acroy / expokit.jl Goto Github PK

In PR #30 , we find the variable anorm = norm(A, Inf) in expmv function does not affect the correctness and performance of the program.

As norm function is not designed for a general Linear Operator, it is natural to ask what it is used for and whether it should be required?

I was using Expokit to compute some unitary transformations. However, I found out that the output is wrong "sometime".

using Expokit

badk = zeros(Int64, 0)
badx = Array{Array{Float64,1}}(0)
for k in 1:1000
    x = expmv(-0.1, sparse(9:-1:1, 1:9, -4:4), ones(9))
    if norm(x) > 4
        push!(badk, k)
        push!(badx, x)
    end
end

length(badk)

sparse(9:-1:1, 1:9, -4:4) is just a 9x9 skew-symmetric matrix, then its exponential is a unitary/orthonormal matrix. Therefore norm(x) should be approximately norm(ones(9)) = 3.

Depending on you luck, you will get non-empty badk above.
It happens on my Macbook as well as my office linux computing server.

The function signature in https://github.com/acroy/Expokit.jl/blob/master/src/expmv.jl#L44
has w::Vector. This could be relaxted to w::AbstractVector to apply expmv! to views, I believe.

In Julia v0.3.5 I get:

In[41]: expmv(π/4,[0 -1im; 1im 0]|>sparse,[1.;0.])
InexactError()
while loading In[41], in expression starting on line 1

 in setindex! at array.jl:307
 in A_mul_B! at linalg/sparse.jl:24
 in expmv! at /home/msilva/.julia/v0.3/Expokit/src/Expokit.jl:71
 in expmv! at /home/msilva/.julia/v0.3/Expokit/src/Expokit.jl:20

In Julia 1.0 the matrix exponential is done via exp and not expm as before. To have a consistent naming I would suggest to rename expmv and phimv.

It looks like may be using norm(matrix). In Julia 0.7, this will compute the Frobenius norm (vecnorm in Julia 0.6), due to JuliaLang/julia#27401. If you want the induced/operator norm as in Julia 0.6, use opnorm(matrix) instead, or Compat.opnorm(matrix) to work in 0.6 and 0.7 (JuliaLang/Compat.jl#577).

Note that, for testing purposes, rather than @test norm(A - B) ≤ tol, it is usually preferred to do @test A ≈ B or @test A ≈ B rtol=... (which uses isapprox).

A bug that was fixed probably in #25 still exists in the current release version of Expokit.jl that people download with Pkg.add("Expokit"). I suggest to make a new release soon.

(the error was found again here)

As briefly discussed in #9 and #18 having benchmarks of a couple of "real world problems" would be nice. Some suitable examples are provided in this paper, among them the 3D quantum harmonic oscillator. Additionally, @mforets suggested to "use Hurwitz-stable matrices arising from control systems".

New benchmarks should go into the folder benchmarks/ and this issue can be used for discussions.

It would be nice if there could be an optional argument for the Krylov vectors so that way if this is called in a loop the caches don't need to be created each time.

as a step forward what do you think about:

• add more functionality:
  • phimv #11
  • padm (a preliminary version is here, and some results here) #14
  • chbv #15
• make docs
• register the package
• benchmarks #21

CC: @acroy, @garrison

The phimv function, phiv in Expokit, computes the solution of a nonhomogeneous linear ODE without direct evaluation of a matrix exponential.

As started in another thread, the expmv() function can also be written in the Al-Mohy's algorithm (A. H. Al-Mohy and N. J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. MIMS EPrint 2010.30, The University of Manchester, 2010; SIAM J. Sci. Comput.), beyond the Krylov, Chebyshev and other usual methods. The MATLAB implementation of the Al-Mohy algorithm (with intensive for-loops) and a C rewritten can be found here and here. It would be good to compare those different methods in Julia and provide parametric interface to competitive algorithms. See also the roadmap of IterativeSolvers.jl.

This happens even for 2x2 matrices, and affects complex and real matrices equally

julia> expm(pi/4*[0 1; 1 0])*[1.;0.]
2-element Array{Float64,1}:
 1.32461 
 0.868671

julia> expmv(pi/4,[0 1; 1 0]|>sparse, [1.;0.])
2-element Array{Float64,1}:
 1.19636
 0.65672

julia> expm(-pi/4*1im*[0 1; 1 0])*[1.;0.]
2-element Array{Complex{Float64},1}:
 0.707107+0.0im     
      0.0-0.707107im

julia> expmv(pi/4,-1im*[0 1; 1 0]|>sparse, [1.;0.])
2-element Array{Complex{Float64},1}:
 0.815705+0.0im     
      0.0-0.578469im

The function padm from Expokit computes the matrix exponential exp(A) of a large sparse matrix using Pade approximants.

A preliminary version is in this branch. Some results are here.

When trying to avoid allocations, note that \ actually allocates. You will instead want to use A_ldiv_B! for the non-allocating form.

Maybe we should also add chbv for completeness. Seems to be straight-forward.

Like for expmv only mul! should be used.