Efficiency improvement of matrix exponential

[jarlebring]

The matrix exponential exp!(A::StridedMatrix{T}) can be improved by appropriate use of Inn=UniformScaling(one(T)) instead of

https://github.com/JuliaLang/julia/blob/697e782ab86bfcdd7fd15550241fe162c51d9f98/stdlib/LinearAlgebra/src/dense.jl#L561

This should ....

• Lessen the cost of some matrix additions (e.g. in U = A * (A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ CC[8].*A6 .+ CC[6].*A4 .+ CC[4].*A2 .+ CC[2].*Inn)
• Lessen the cost of some scalar multiplication e.g. P = copy(Inn); U = C[2] * P
• (Most importantly) Reduce the number of matrix-matrix multiplications since in the first round in this loop, P=I:
  https://github.com/JuliaLang/julia/blob/697e782ab86bfcdd7fd15550241fe162c51d9f98/stdlib/LinearAlgebra/src/dense.jl#L581-L583

[jarlebring]

The nof mat-mat mult is important. The implementation really has one too many mat-mats:

julia> using LinearAlgebra; import Base.*;
julia> function *(A::Matrix{Float64},B::Matrix{Float64}) ; 
    if (norm(B-I)==0 || norm(A-I)==0) 
         println("Identity mult"); 
   end; 
   C=deepcopy(B);
   return mul!(C,A,B)  # Call original function
end
julia> A=randn(100,100); A=A/opnorm(A,1);
julia> LinearAlgebra.exp!(A);
Identity mult
julia> 

[musm]

@jarlebring looks like a reasonable change. Do you intend to open a PR with these changes?

[jarlebring]

From my perspective it's up for grabs. I would probably do it eventually if nobody else provides a PR.

[jebej]

The problem with making P = I is it would make the function type unstable as soon as you multiply it with A2. It might be possible to avoid the unnecessary multiplication but it won't be as simple as that.

For the matrix addition, the issue is that for the moment I doesn't work with fusion/broadcasting.

Since you are looping over many matrices for the addition anyway, looping through the zero entries of Inn ends up being faster (I think) than either 1. doing a separate loop for addition, or 2. implementing the branching logic for I to only add to diagonal elements in the fusion.

[jarlebring]

Since you are looping over many matrices for the addition anyway, looping through the zero entries of Inn ends up being faster (I think) than either 1. doing a separate loop for addition, or 2. implementing the branching logic for I to only add to diagonal elements in the fusion.

I find that hard to believe. Looping through all zero elements is O(n^2), and both 1 and 2 are O(n). It sounds difficult to beat that complexity difference even with proper for-loops. Compare for instance with all the variants here: https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/5

[jebej]

I find that hard to believe. Looping through all zero elements is O(n^2), and both 1 and 2 are O(n). It sounds difficult to beat that complexity difference even with proper for-loops.

This complexity argument would apply only if each addition operation was truly independent and had a fixed cost. But as you might have noticed in your PR, the cost of adding additional matrices becomes smaller the more matrices you add (due to fusion and SIMD). Thus, adding 4 matrices is only 1.5x slower (on my computer) than adding 2 matrices, despite there being 3x the number of additions:

julia> @btime $A .+ $B;
  12.000 μs (2 allocations: 78.20 KiB)

julia> @btime $A .+ $B .+ $C;
  15.700 μs (2 allocations: 78.20 KiB)

julia> @btime $A .+ $B .+ $C .+ $D;
  17.799 μs (2 allocations: 78.20 KiB)

[jebej]

Though with your background you would know this already, sorry about that! It might be the case that for larger matrices separating the O(n) additions does end up being faster.

[jarlebring]

It might be the case that for larger matrices separating the O(n) additions does end up being faster.

Indeed. It was faster to do the O(n) operations for larger matrices, although not very much faster in a relative sense. This is due to the fact that the function contains several matrix-matrix products which are dominating.