Matrix `exp` and `log` are not inverse of each other

[olivierverdier]

1. Output of versioninfo()

Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin22.4.0)
CPU: 8 × Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_EDITOR = emacsclient
JULIA_INPUT_COLOR = bold
JULIA_ANSWER_COLOR = bold

2. How Julia was installed

With juliaup

3. Minimal working example:

M = [0.9949357359852791 -0.015567763143324862 -0.09091193493947397 -0.03994428739762443 0.07338356301650806; 0.011813655598647289 0.9968988574699793 -0.06204555000202496 0.04694097614450692 0.09028834462782365; 0.092737943594701 0.059546719185135925 0.9935850721633324 0.025348893985651405 -0.018530261590167685; 0.0369187299165628 -0.04903571106913449 -0.025962938675946543 0.9977767446862031 0.12901494726320517; 0.0 0.0 0.0 0.0 1.0]
sdiff = exp(log(M)) - M
-log10(maximum(abs.(sdiff))) # < 3

The same test in Python gives

import numpy as np
import scipy as sp


M = np.array([[0.9949357359852791,-0.015567763143324862,-0.09091193493947397,-0.03994428739762443,0.07338356301650806],[0.011813655598647289,0.9968988574699793,-0.06204555000202496,0.04694097614450692,0.09028834462782365],[0.092737943594701,0.059546719185135925,0.9935850721633324,0.025348893985651405,-0.018530261590167685],[0.0369187299165628,-0.04903571106913449,-0.025962938675946543,0.9977767446862031,0.12901494726320517],[0.0,0.0,0.0,0.0,1.0]])

sdiff = sp.linalg.expm(sp.linalg.logm(M)) - M

print(-np.log10(np.max(np.abs(sdiff)))) # > 15

[dkarrasch]

Could you please quickly check which one seems to be wrong: log or exp?

[oscardssmith]

reproduced on x86 linux. The issue appears to be in log (although it's somewhat hard to prove because matrix logs are not unique).

[olivierverdier]

A simpler example:

N = [0.0 -0.8 1.1 0.2 -0.2; 0.8 0.0 0.8 -0.4 -0.1; -1.1 -0.8 0.0 0.4 -0.1; -0.2 0.4 -0.4 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]

k = 10 # also bad for k = 100 or 1000

sdiff = log(exp(N/k)) - N/k
@show -log10(maximum(abs.(sdiff))) # 2.9

I agree with @oscardssmith , I suspect the matrix logarithm to be at fault.

---

Edit: for completeness, the Python code

N = np.array([[0.0, -0.8, 1.1, 0.2, -0.2], [0.8, 0.0, 0.8, -0.4, -0.1], [-1.1, -0.8, 0.0, 0.4, -0.1], [-0.2, 0.4, -0.4, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0],])

k = 10

sdiff = sp.linalg.logm(sp.linalg.expm(N/k)) - N/k
print(-np.log10(np.max(np.abs(sdiff)))) # 15.0

[oscardssmith]

This second example isn't as good because iiuc, log(exp(M)) == M isn't expected for matrices do to branch cuts.

[dkarrasch]

On v1.6.7 I get

julia> -log10(maximum(abs.(sdiff))) # < 3
15.175451669208671

on v1.9.4 and v1.8.5

julia> -log10(maximum(abs.(sdiff))) # < 3
2.233375639545474

so this must be an "old" bug.

[olivierverdier]

This second example isn't as good because iiuc, log(exp(M)) == M isn't expected for matrices do to branch cuts.

That's why I included the parameter k. If M is small enough, I would expect log(exp(M)) == M to be true no matter what, right?

[olivierverdier]

so this must be an "old" bug.

I suspect the problem to be in log_quasitriu, which was introduced in 6913f9ccb156230f54830c8b7b8ef82b41cac27e, in 2021.

[giordano]

I suspect the problem to be in log_quasitriu, which was introduced in JuliaLang/julia@6913f9c, in 2021.

That's more recent than Julia v1.6 though

[olivierverdier]

More evidence that the problem is the logarithm:

N = [0.0 -0.8 1.1 0.2 -0.2; 0.8 0.0 0.8 -0.4 -0.1; -1.1 -0.8 0.0 0.4 -0.1; -0.2 0.4 -0.4 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]
k = 10
M = exp(N/k)

sdiff = log(M*M) - 2*log(M)

@show -log10(maximum(abs.(sdiff))) # 2.6

The equivalent code in Python gives 14.7.

[olivierverdier]

I suspect the problem to be in log_quasitriu, which was introduced in 6913f9c, in 2021.

That's more recent than Julia v1.6 though

If it were the case that PR JuliaLang/julia#39973 is to blame (wild guess), then the bug would have appeared first in v1.7.0.

[olivierverdier]

For completeness, I checked the code above (i.e., with log(M*M) - 2*log(M)) with various versions, and indeed, the bug appears between v1.6.7 and v1.7.3.
My money is still on PR JuliaLang/julia#39973.

[olivierverdier]

An even simpler test:

function test_log(N, k=10)
    M = exp(N / k)
    sdiff = log(M*M) - 2*log(M)
    return -log10(maximum(abs.(sdiff)))
end

Now simply with

N = [0.0 -1.0 1.0; 1.0 0.0 0.0; 0.0 0.0 0.0]

run

test_log(N) # 2.1

and the same in Python returns 15.4.

[dkarrasch]

It's not clear to me what is the cause. I tried a specific test from scipy's test suite, and this is the result:

julia> VERSION
v"1.10.4"

julia> A = [3.2346e-01 3.0000e+04 3.0000e+04 3.0000e+04;
              0.0000e+00 3.0089e-01 3.0000e+04 3.0000e+04;
              0.0000e+00 0.0000e+00 3.2210e-01 3.0000e+04;
              0.0000e+00 0.0000e+00 0.0000e+00 3.0744e-01];

julia> logA = [ -1.12867982029050462e+00  9.61418377142025565e+04 -4.52485573953179264e+09  2.92496941103871812e+14;
         0.00000000000000000e+00 -1.20101052953082288e+00 9.63469687211303099e+04 -4.68104828911105442e+09;
         0.00000000000000000e+00  0.00000000000000000e+00 -1.13289322264498393e+00  9.53249183094775653e+04;
         0.00000000000000000e+00  0.00000000000000000e+00 0.00000000000000000e+00 -1.17947533272554850e+00];

julia> (logA - log(A)) / norm(logA)
4×4 Matrix{Float64}:
 0.0  5.97008e-25  9.78138e-21  -1.06839e-15
 0.0  0.0          4.97507e-25   1.30418e-20
 0.0  0.0          0.0           2.98504e-25
 0.0  0.0          0.0           0.0

julia> (exp(logA) - A) / norm(A)
4×4 Matrix{Float64}:
 -1.81534e-8  0.000297088  50.5783       -2.79971e6
  0.0         2.33977e-8    0.00116877    8.26782
  0.0         0.0           3.53834e-10  -0.00160443
  0.0         0.0           0.0          -3.34387e-8

julia> (exp(log(A)) - A) / norm(A)
4×4 Matrix{Float64}:
 -1.81534e-8  0.000297088  50.5783       -2.79971e6
  0.0         2.33977e-8    0.00116877    8.26782
  0.0         0.0           3.53834e-10  -0.00160443
  0.0         0.0           0.0          -3.34387e-8

julia> (exp(log(A)) - A)
4×4 Matrix{Float64}:
 -0.001334  21.8314       3.71673e6     -2.05736e11
  0.0        0.00171937  85.8868         6.07558e5
  0.0        0.0          2.60014e-5  -117.901
  0.0        0.0          0.0           -0.00245723

OTOH, we have:

>>> (sp.linalg.expm(A_logm) - A) / sp.linalg.norm(A)
array([[ 0.00000000e+00,  1.33667877e-15,  3.49612787e-10,
        -5.04994630e-07],
       [ 0.00000000e+00,  0.00000000e+00,  1.23766552e-15,
         4.78558722e-11],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         1.38618539e-15],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00]])

So, the computed log is really close to the "known" solution, but the exp of it is completely wrong in the upper right corner.

ADDENDUM: It seems, however, that this specific example fails on all versions starting from 1.6.x.

[mikmoore]

I tried a specific test from scipy's test suite

It's worth noting that the condition number of that matrix is awful:

julia> cond(A)
1.8884133342622014e20

julia> cond(logA)
4.086839610043594e28

I have an independent myexp implementation (our actual exp implementation is better but has limited type support JuliaLang/julia#51008) and it similarly struggles with this:

function myexp(A)
	# compute exp(A) == exp(A/2^n)^(2^n)
	iszero(A) && return float(one(A))
	# TODO: decide the best norm to use here: op-1, op-Inf, or frobenius are probably good candidates - maybe take the min of all three
	nA = float(norm(A,2))
	fr,ex = frexp(nA)
	logscale = max(4 + ex,0) # ensure nA^(4+1)/factorial(4+1) < eps for 4 or whatever power our approximation is
	p = ntuple(i->inv(prod(1.0:i-1)),Val(9))
	scale = exp2(-logscale)
	sA = A * scale
	if sA isa AbstractMatrix
		expA = evalpoly(sA,UniformScaling.(p))
	else
		expA = evalpoly(sA,p)
	end
	for _ in 1:logscale
		expA = expA*expA
	end
	return expA
end

And with this and BigFloat:

julia> norm(exp(logA) - A) / norm(A)
2.799714043033207e6

julia> norm(myexp(logA) - A) / norm(A) # myexp is generally worse
5.134964610902018e6

julia> norm(myexp(big.(logA)) - A) / norm(A) |> Float64 # despite being worse, BigFloat makes this version work
8.141481015780612e-7

Note that the squaring part of our exp routine is probably very-much to blame, given that BigFloat resolves that problem even though the exp kernel is not expanded beyond what is roughly needed for Float64. I don't think we should use that example for this particular issue.

However, myexp does nothing to resolve the original example:

julia> cond(M) # extremely well conditioned
1.1906431494194285

julia> norm(exp(log(M)) - M) / norm(M) |> Float64
0.003052024023511132

julia> norm(myexp(big.(log(M))) - M) / norm(M) |> Float64
0.0030520240235111305

Given that the matrix exponential is reasonably "simple" to compute (modulo stability issues, which even the literature has not fully settled), I am much more inclined towards expecting that log is the culprit.