Implement matrix logarithm

[eschnett]

Closes #993.

[eschnett]

The tests currently fail because they test 0x0 matrices as well, and this is only implemented in https://github.com/JuliaArrays/StaticArrays.jl/pull/991. Waiting for this to be applied, then the tests here should succeed.

[eschnett]

This PR is ready to be merged.

[hyrodium]

The log function seems not type-stable, and sometimes throws DomainError.

julia> m0 = one(SMatrix{3,3})
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> log(m0)
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0

julia> m1 = rand(SMatrix{3,3})
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.273815   0.156067  0.538108
 0.0792263  0.149036  0.6332
 0.940688   0.6957    0.550016

julia> log(m1)
ERROR: DomainError with -0.5292055332168919:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:304
  [3] log
    @ ./special/log.jl:269 [inlined]
  [4] macro expansion
    @ ~/.julia/dev/StaticArrays/src/broadcast.jl:126 [inlined]
  [5] _broadcast
    @ ~/.julia/dev/StaticArrays/src/broadcast.jl:100 [inlined]
  [6] copy
    @ ~/.julia/dev/StaticArrays/src/broadcast.jl:27 [inlined]
  [7] materialize
    @ ./broadcast.jl:860 [inlined]
  [8] _log(#unused#::Size{(3, 3)}, A::SMatrix{3, 3, Float64, 9})
    @ StaticArrays ~/.julia/dev/StaticArrays/src/logm.jl:19
  [9] log(A::SMatrix{3, 3, Float64, 9})
    @ StaticArrays ~/.julia/dev/StaticArrays/src/logm.jl:1
 [10] top-level scope
    @ REPL[46]:1

julia> m2 = rand(SMatrix{3,3})
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.582469   0.912389   0.435272
 0.0275857  0.765842   0.944147
 0.73664    0.0436018  0.405992

julia> log(m2)
3×3 SMatrix{3, 3, ComplexF64, 9} with indices SOneTo(3)×SOneTo(3):
 -0.280447+0.0im    1.10906+0.0im   -0.20155+0.0im
 -0.696253+0.0im   0.239951+0.0im    1.33066+0.0im
   1.18863+0.0im  -0.643787+0.0im  -0.582311+0.0im

julia> m3 = rand(SMatrix{3,3})
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.239263  0.566921  0.172025
 0.219623  0.874731  0.0596628
 0.436071  0.930055  0.503457

julia> log(m3)
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -2.85822    1.32526    0.790975
  0.691332  -0.388829  -0.0421343
  1.41868    1.03391   -1.04496

julia> log(@SMatrix randn(3,3))
3×3 SMatrix{3, 3, ComplexF64, 9} with indices SOneTo(3)×SOneTo(3):
  -4.89227+3.05795im     -1.89937+0.792861im    -0.568245+0.491583im
 -0.469905+0.313737im   -0.424589+0.0813454im    -1.86939+0.0504351im
 -0.470829+0.0142845im    4.29088+0.00370366im    1.46013+0.00229632im

[eschnett]

@hyrodium I confirm. Apparently the eigendecomposition is not type stable. I expected it to be.

The domain error is intended: If there is no real-valued matrix logarithm, then there should be a domain error, because returning a complex result would lead to a type instability.

The implementation special cases small matrices. This is where the domain error is generated. Larger matrices use the eigendecomposition which is causing the type instability.

[eschnett]

How should we handle type instabilities? For example, currently:

• sqrt returns an error for indefinite matrices
• eigen is type unstable and might convert to complex numbers

I think there is much value in type-stable algorithms for small matrices. For example, I want to store small matrices in arrays, and type instabilities would lead to very slow code.

I suggest the following:

• all matrix operations are type-stable
• "standard" calls return errors for real arguments if the result would be complex
• we introduce an additional, optional first argument that defines the return type, allowing conversion to complex numbers, similar to round

This would look like:

• log(A::SMatrix{...,<:Complex}) always succeeds
• log(A::SMatrix{...,<:Real}) might fail
• log(Real, A::SMatrix{...<:Real}) has the same behaviour
• log(Complex, A::SMatrix{...<:Real}) always succeeds, always returning a complex result

The return type of log(A::SMatrix) is typeof(log(one(eltype(A)))). The return type of log(::Type{T}, A::SMatrix) is typeof(log(T(one(eltype(A))))).

[hyrodium]

• all matrix operations are type-stable
• "standard" calls return errors for real arguments if the result would be complex

I agree with that :+1:

• we introduce an additional, optional first argument that defines the return type, allowing conversion to complex numbers, similar to round

I'm not sure we need this method here. These are what I thought:

• We don't have log(::Type, ::Number) method in Base.
• The round(T,x) is almost same as T(round(x)), but log(T,A) will not be like that; it converts eltype.
• We can just use log(complex(A::SMatrix{3,3,<:Real})) instead of log(Complex, A:SMatrix{3,3,<:Real}). Or, is there any performance advantage with log(Complex, A::SMatrix{...<:Real})?

[eschnett]

I thought there might be performance advantages. However, in a preliminary implementation, I'm just calling Complex.(A) anyway, so let's not introduce that new signature.

Currently, the eigendecomposition is not type stable. Should we make it type stable in the same way? This could be a breaking change.

[thchr]

I think having a lot of StaticArrays-specific matrix functions is undesirable (e.g. with return-type as a part of the signature): a good chunk of the appeal of StaticArrays is that they can be very nearly "plug-and-play" swapped for standard arrays.

Would it be crazy to simply promote all results to the corresponding complex type? If someone knows the output element type should really be <:Real, they could always wrap the input in Hermitian (e.g. for eigen).

[eschnett]

There are several related matrix functions, e.g. exp, log, sqrt. I think they should all behave in a similar manner regarding the return type. I would prefer that the return type to be the same as the input type, as is currently the case for exp and sqrt. Having log return complex numbers by default would be strange.

My current experimental implementation does just as you suggest: convert everything to complex, then calculate the logarithm, then convert back to real. The problem is that there seem to be many cases where the conversion to real "should" work, but doesn't, purely due to round-off. These are matrices with all real entries (even integer entries), and Mathematica easily finds an all-real logarithm. However, the eigenvalues and eigenvectors are complex. The imaginary part that we find (due to round-off) is about 1e-10, which is too large for my taste to discard silently. log([4 2; -2 1]) is an example of this.

I think that my approach to calculate the matrix log via eigenvalues is too naive to work well in a numerical setting. A Schur decomposition might be necessary.

[MasonProtter]

The problem here is basically that base makes these functions type unstable because it's known that the overhead of working with regular heavy arrays makes that instability okay to work with.

I think it's important for static arrays that they be type stable, because their whole reason to exist is for getting top-notch performance on small array operations. To that end, I'd advocate in favour of making these methods either error, or return NaNs in the type unstable case (I'd actually prefer NaN, but I know that's likely an upopular opinion so I'd rather separate it out from the rest of this point which I think is less controversial).

I think that it just needs to be documented that these things can error (or return nan) if the user does not ask for a suitably wide (i.e. complex) type. I do not think that the matrix logarithm should be complex by default.

[eschnett]

I agree regarding type stability.

The current problem (that made me stop working on this) is that the Base algorithms require allocating regular dense Matrix objects. I don't think that's reasonable, I want to avoid allocation. This would force me to re-implement quite a few algorithms to handle the cases that interest me (symmetric or not, real or complex).

As far as I can see, calculating the matrix logarithm essentially requires an eigenvalue decomposition. These would be good to have anyway.