RFC: Matrix logarithm function (issue #5840)

base/exports.jl

@@ -659,6 +659,7 @@ export
     linreg,
     logabsdet,
     logdet,
+    logm,
     lu,
     lufact!,
     lufact,

base/linalg.jl

@@ -85,6 +85,7 @@ export
     linreg,
     logabsdet,
     logdet,
+    logm,
     lu,
     lufact,
     lufact!,

base/linalg/dense.jl

@@ -278,6 +278,39 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
     end
 end
 
+function logm(A::StridedMatrix)
+    # If possible, use diagonalization
+    if ishermitian(A)
+        return logm(Hermitian(A))
+    end
+
+    # Use Schur decomposition
+    n = chksquare(A)
+    S,Q,d = schur(complex(A))
+    R = logm(UpperTriangular(S))
+    retmat = Q * R * Q'
+
+    # Check whether the matrix has nonpositive real eigs
+    np_real_eigs = false
+    for i = 1:n
+        if imag(d[i]) < eps() && real(d[i]) <= 0
+            np_real_eigs = true
+            break
+        end
+    end
+    if np_real_eigs
+        warn("Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
+    end
+
+    if isreal(A) && ~np_real_eigs
+        return real(retmat)
+    else
+        return retmat
+    end
+end
+logm(a::Number) = (b = log(complex(a)); imag(b) == 0 ? real(b) : b)
+logm(a::Complex) = log(a)
+
 function sqrtm{T<:Real}(A::StridedMatrix{T})
     if issym(A)
         return sqrtm(Symmetric(A))

base/linalg/symmetric.jl

@@ -113,6 +113,11 @@ end
 
 #Matrix-valued functions
 expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
+function logm{T<:Real}(A::RealHermSymComplexHerm{T})
+    F = eigfact(A)
+    isposdef(F) && return F.vectors*Diagonal(log(F.values))*F.vectors'
+    return F.vectors*Diagonal(log(complex(F.values)))*F.vectors'
+end
 function sqrtm{T<:Real}(A::RealHermSymComplexHerm{T})
     F = eigfact(A)
     isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'

base/linalg/triangular.jl

@@ -889,6 +889,193 @@ for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv
     end
 end
 
+# Complex matrix logarithm for the upper triangular factor, see:
+#   Al-Mohy and Higham, "Improved inverse  scaling and squaring algorithms for
+# the matrix logarithm", SIAM J. Sci. Comput., 34(4), (2012), pp. C153-C169.
+#   Al-Mohy, Higham and Relton, "Computing the Frechet derivative of the matrix
+# logarithm and estimating the condition number", SIAM J. Sci. Comput., 35(4),
+# (2013), C394-C410.
+#   http://eprints.ma.man.ac.uk/1687/02/logm.zip
+function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T})
+    maxsqrt = 100
+    theta = [1.586970738772063e-005,
+         2.313807884242979e-003,
+         1.938179313533253e-002,
+         6.209171588994762e-002,
+         1.276404810806775e-001,
+         2.060962623452836e-001,
+         2.879093714241194e-001]
+    tmax = size(theta, 1)
+    n = size(A0, 1)
+    A = A0
+    p = 0
+    m = 0
+
+    # Compute repeated roots
+    d = diag(A)
+    dm1 = Array(T, n)
+    s = 0
+    for i = 1:n
+        dm1[i] = dm1[i] - 1.
+    end
+    while norm(dm1, Inf) > theta[tmax]
+        for i = 1:n
+            d[i] = sqrt(d[i])
+        end
+        for i = 1:n
+            dm1[i] = d[i] - 1
+        end
+        s = s + 1
+    end
+    s0 = s
+    for k = 1:min(s, maxsqrt)
+        A = sqrtm(A)
+    end
+
+    AmI = A - I
+    d2 = sqrt(norm(AmI^2, 1))
+    d3 = cbrt(norm(AmI^3, 1))
+    alpha2 = max(d2, d3)
+    foundm = false
+    if alpha2 <= theta[2]
+        m = alpha2<=theta[1]?1:2
+        foundm = true
+    end
+
+    while ~foundm
+        more = false
+        if s > s0
+            d3 = cbrt(norm(AmI^3, 1))
+        end
+        d4 = norm(AmI^4, 1)^(1/4)
+        alpha3 = max(d3, d4)
+        if alpha3 <= theta[tmax]
+            for j = 3:tmax
+                if alpha3 <= theta[j]
+                    break
+                end
+            end
+            if j <= 6
+                m = j
+                break
+            elseif alpha3 / 2 <= theta[5] && p < 2
+                more = true
+                p = p + 1
+           end
+        end
+
+        if ~more
+            d5 = norm(AmI^5, 1)^(1/5)
+            alpha4 = max(d4, d5);
+            eta = min(alpha3, alpha4);
+            if eta <= theta[tmax]
+                j = 0
+                for j = 6:tmax
+                    if eta <= theta[j]
+                        m = j
+                        break
+                    end
+                end
+                break
+            end
+        end
+
+        if s == maxsqrt
+            m = tmax
+            break
+        end
+        A = sqrtm(A)
+        AmI = A - I
+        s = s + 1
+    end
+
+    # Compute accurate superdiagonal of T
+    p = 1 / 2^s
+    for k = 1:n-1
+        Ak = A0[k,k]
+        Akp1 = A0[k+1,k+1]
+        Akp = Ak^p
+        Akp1p = Akp1^p
+        A[k,k] = Akp
+        A[k+1,k+1] = Akp1p
+        if Ak == Akp1
+            A[k,k+1] = p * A0[k,k+1] * Ak^(p-1)
+        elseif 2 * abs(Ak) < abs(Akp1) || 2 * abs(Akp1) < abs(Ak)
+            A[k,k+1] = A0[k,k+1] * (Akp1p - Akp) / (Akp1 - Ak)
+        else
+            logAk = log(Ak)
+            logAkp1 = log(Akp1)
+            w = atanh((Akp1 - Ak)/(Akp1 + Ak)) + im*pi*ceil((imag(logAkp1-logAk)-pi)/(2*pi))
+            dd = 2 * exp(p*(logAk+logAkp1)/2) * sinh(p*w) / (Akp1 - Ak);
+            A[k,k+1] = A0[k,k+1] * dd
+        end
+    end
+
+    # Compute accurate diagonal of T
+    for i = 1:n
+        a = A0[i,i]
+        if s == 0
+            r = a - 1
+        end
+        s0 = s
+        if angle(a) >= pi / 2
+            a = sqrt(a)
+            s0 = s - 1
+        end
+        z0 = a - 1
+        a = sqrt(a)
+        r = 1 + a
+        for j = 1:s0-1
+            a = sqrt(a)
+            r = r * (1 + a)
+        end
+        A[i,i] = z0 / r
+    end
+
+    # Compute the Gauss-Legendre quadrature formula
+    R = zeros(Float64, m, m)
+    for i = 1:m - 1
+        R[i,i+1] = i / sqrt((2 * i)^2 - 1)
+        R[i+1,i] = R[i,i+1]
+    end
+    x,V = eig(R)
+    w = Array(Float64, m)
+    for i = 1:m
+        x[i] = (x[i] + 1) / 2
+        w[i] = V[1,i]^2
+    end
+
+    # Compute the Padé approximation
+    Y = zeros(T, n, n)
+    for k = 1:m
+        Y = Y + w[k] * (A / (x[k] * A + I))
+    end
+
+    # Scale back
+    scale!(2^s, Y)
+
+    # Compute accurate diagonal and superdiagonal of log(T)
+    for k = 1:n-1
+        Ak = A0[k,k]
+        Akp1 = A0[k+1,k+1]
+        logAk = log(Ak)
+        logAkp1 = log(Akp1)
+        Y[k,k] = logAk
+        Y[k+1,k+1] = logAkp1
+        if Ak == Akp1
+            Y[k,k+1] = A0[k,k+1] / Ak
+        elseif 2 * abs(Ak) < abs(Akp1) || 2 * abs(Akp1) < abs(Ak)
+            Y[k,k+1] = A0[k,k+1] * (logAkp1 - logAk) / (Akp1 - Ak)
+        else
+            w = atanh((Akp1 - Ak)/(Akp1 + Ak) + im*pi*(ceil((imag(logAkp1-logAk) - pi)/(2*pi))))
+            Y[k,k+1] = 2 * A0[k,k+1] * w / (Akp1 - Ak)
+        end
+    end
+
+    return UpperTriangular(Y)
+
+end
+
 function sqrtm{T}(A::UpperTriangular{T})
     n = size(A, 1)
     TT = typeof(sqrt(zero(T)))

doc/stdlib/linalg.rst

@@ -646,6 +646,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Matrix exponential.
 
+.. function:: logm(A)
+
+   Matrix logarithm.
+
 .. function:: lyap(A, C)
 
    Computes the solution ``X`` to the continuous Lyapunov equation ``AX + XA' + C = 0``, where no eigenvalue of ``A`` has a zero real part and no two eigenvalues are negative complex conjugates of each other.

test/linalg3.jl

@@ -154,29 +154,49 @@ for elty in (Float32, Float64, Complex64, Complex128)
                                  0.135335281175235 0.406005843524598 0.541341126763207]')
     @test_approx_eq expm(A3) eA3
 
-    # issue 5116
-    A4  = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
-    eA4 = [-0.999786072879326  -0.065407069689389   0.0   0.0
-           0.006540706968939  -0.999786072879326   0.0   0.0
-           0.0                 0.0                 1.0   0.0
-           0.013081413937878  -3.999572145758650   0.0   1.0]
-    @test_approx_eq expm(A4) eA4
-
-    # issue 5116
-    A5  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
-    eA5 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
-            0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
-            0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
-            0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
-    @test_approx_eq expm(A5) eA5
-
     # Hessenberg
     @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
                                              [4.000000000000000  -1.414213562373094  -1.414213562373095
                                               -1.414213562373095   4.999999999999996  -0.000000000000000
                                               0  -0.000000000000002   3.000000000000000])
 end
 
+for elty in (Float64, Complex{Float64})
+    A4  = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
+                                 1/3 1/4 1/5 1/6;
+                                 1/4 1/5 1/6 1/7;
+                                 1/5 1/6 1/7 1/8])
+    @test_approx_eq expm(logm(A4)) A4
+
+    A5  = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
+    @test_approx_eq expm(logm(A5)) A5
+
+    A6  = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
+    @test_approx_eq expm(logm(A6)) A6
+
+    A7  = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
+    @test_approx_eq expm(logm(A7)) A7
+end
+
+A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1];
+@test_approx_eq expm(logm(A8)) A8
+
+# issue 5116
+A9  = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
+eA9 = [-0.999786072879326  -0.065407069689389   0.0   0.0
+       0.006540706968939  -0.999786072879326   0.0   0.0
+       0.0                 0.0                 1.0   0.0
+       0.013081413937878  -3.999572145758650   0.0   1.0]
+@test_approx_eq expm(A9) eA9
+
+# issue 5116
+A10  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
+eA10 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
+        0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
+        0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
+        0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
+@test_approx_eq expm(A10) eA10
+
 # matmul for types w/o sizeof (issue #1282)
 A = Array(Complex{Int},10,10)
 A[:] = complex(1,1)