From 4d9bffae8ca7e9d42a977dbaa4c3502edc3ffb4c Mon Sep 17 00:00:00 2001 From: Massimiliano Fasi Date: Sun, 19 Jul 2015 18:10:23 +0100 Subject: [PATCH 1/2] Matrix logarithm function logm added --- base/exports.jl | 1 + base/linalg.jl | 1 + base/linalg/dense.jl | 37 +++++++- base/linalg/symmetric.jl | 6 +- base/linalg/triangular.jl | 195 ++++++++++++++++++++++++++++++++++++++ doc/stdlib/linalg.rst | 4 + test/linalg3.jl | 48 +++++++--- 7 files changed, 272 insertions(+), 20 deletions(-) diff --git a/base/exports.jl b/base/exports.jl index 7960ff9bc8f7e..94bcd18396cf4 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -659,6 +659,7 @@ export linreg, logabsdet, logdet, + logm, lu, lufact!, lufact, diff --git a/base/linalg.jl b/base/linalg.jl index 66fbd72440bad..ef57347fc3105 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -85,6 +85,7 @@ export linreg, logabsdet, logdet, + logm, lu, lufact, lufact!, diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index b05b533114ba7..1646dc09486f1 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -278,9 +278,42 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T}) end end +function logm{T}(A::StridedMatrix{T}) + # If possible, use diagonalization + if ishermitian(A) + return funm(log, Hermitian(A)) + end + + # Use Schur decomposition + n = chksquare(A) + S,Q,d = schur(complex(A)) + R = full(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)) + return funm(sqrt, Symmetric(A)) end n = chksquare(A) SchurF = schurfact(complex(A)) @@ -290,7 +323,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T}) end function sqrtm{T<:Complex}(A::StridedMatrix{T}) if ishermitian(A) - return sqrtm(Hermitian(A)) + return funm(sqrt, Hermitian(A)) end n = chksquare(A) SchurF = schurfact(A) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 9a209007f3510..9c2addbf7b4a1 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -113,8 +113,8 @@ end #Matrix-valued functions expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors') -function sqrtm{T<:Real}(A::RealHermSymComplexHerm{T}) +function funm{T<:Real}(f::Function, A::RealHermSymComplexHerm{T}) F = eigfact(A) - isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors' - return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors' + isposdef(F) && return F.vectors*Diagonal(f(F.values))*F.vectors' + return F.vectors*Diagonal(f(complex(F.values)))*F.vectors' end diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 4c5645e1bf340..607316e4b50f0 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -889,6 +889,201 @@ 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. +function logm{T}(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 - UpperTriangular(eye(n)) + d2 = norm(AmI^2, 1)^(1/2) + d3 = norm(AmI^3, 1)^(1/3) + 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 = norm(AmI^3, 1)^(1/3) + 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 - UpperTriangular(eye(n)) + 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 abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2 + 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 + B = zeros(T, n, n) + Y = zeros(T, n, n) + for k = 1:m + for j = 1:n + B[j,j] = 1. + A[j,j] * x[k] + for i = 1:j-1 + B[i,j] = A[i,j] * x[k] + end + end + B = A * inv(B) + scale!(w[k], B) + Y = Y + B + 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 abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2 + 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))) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 4b03ba664d62d..ccb5cc80f56b0 100644 --- a/doc/stdlib/linalg.rst +++ b/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. diff --git a/test/linalg3.jl b/test/linalg3.jl index 1df5e5a25dcdc..becc3b5daca1d 100644 --- a/test/linalg3.jl +++ b/test/linalg3.jl @@ -154,21 +154,20 @@ 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 + A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5; + 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 # Hessenberg @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty}, @@ -177,6 +176,25 @@ for elty in (Float32, Float64, Complex64, Complex128) 0 -0.000000000000002 3.000000000000000]) 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) From 059d8704b7ff4317ee18b5f6c39587d52d30c3ba Mon Sep 17 00:00:00 2001 From: Massimiliano Fasi Date: Sun, 19 Jul 2015 22:54:14 +0100 Subject: [PATCH 2/2] A few minor improvements, high order functions removed, added type checking to ensure the precision to be double --- base/linalg/dense.jl | 12 ++++++------ base/linalg/symmetric.jl | 11 ++++++++--- base/linalg/triangular.jl | 38 +++++++++++++++----------------------- test/linalg3.jl | 16 +++++++++------- 4 files changed, 38 insertions(+), 39 deletions(-) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 1646dc09486f1..b5940c6dfd5b2 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -278,17 +278,17 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T}) end end -function logm{T}(A::StridedMatrix{T}) +function logm(A::StridedMatrix) # If possible, use diagonalization if ishermitian(A) - return funm(log, Hermitian(A)) + return logm(Hermitian(A)) end # Use Schur decomposition n = chksquare(A) S,Q,d = schur(complex(A)) - R = full(logm(UpperTriangular(S))) - retmat = Q*R*Q' + R = logm(UpperTriangular(S)) + retmat = Q * R * Q' # Check whether the matrix has nonpositive real eigs np_real_eigs = false @@ -313,7 +313,7 @@ logm(a::Complex) = log(a) function sqrtm{T<:Real}(A::StridedMatrix{T}) if issym(A) - return funm(sqrt, Symmetric(A)) + return sqrtm(Symmetric(A)) end n = chksquare(A) SchurF = schurfact(complex(A)) @@ -323,7 +323,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T}) end function sqrtm{T<:Complex}(A::StridedMatrix{T}) if ishermitian(A) - return funm(sqrt, Hermitian(A)) + return sqrtm(Hermitian(A)) end n = chksquare(A) SchurF = schurfact(A) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 9c2addbf7b4a1..9396ba27f8e7d 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -113,8 +113,13 @@ end #Matrix-valued functions expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors') -function funm{T<:Real}(f::Function, A::RealHermSymComplexHerm{T}) +function logm{T<:Real}(A::RealHermSymComplexHerm{T}) F = eigfact(A) - isposdef(F) && return F.vectors*Diagonal(f(F.values))*F.vectors' - return F.vectors*Diagonal(f(complex(F.values)))*F.vectors' + 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' + return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors' end diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 607316e4b50f0..199b59eb71605 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -895,7 +895,8 @@ end # 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. -function logm{T}(A0::UpperTriangular{T}) +# 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, @@ -931,9 +932,9 @@ function logm{T}(A0::UpperTriangular{T}) A = sqrtm(A) end - AmI = A - UpperTriangular(eye(n)) - d2 = norm(AmI^2, 1)^(1/2) - d3 = norm(AmI^3, 1)^(1/3) + 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] @@ -944,7 +945,7 @@ function logm{T}(A0::UpperTriangular{T}) while ~foundm more = false if s > s0 - d3 = norm(AmI^3, 1)^(1/3) + d3 = cbrt(norm(AmI^3, 1)) end d4 = norm(AmI^4, 1)^(1/4) alpha3 = max(d3, d4) @@ -957,7 +958,7 @@ function logm{T}(A0::UpperTriangular{T}) if j <= 6 m = j break - elseif alpha3/2 <= theta[5] && p < 2 + elseif alpha3 / 2 <= theta[5] && p < 2 more = true p = p + 1 end @@ -984,12 +985,12 @@ function logm{T}(A0::UpperTriangular{T}) break end A = sqrtm(A) - AmI = A - UpperTriangular(eye(n)) + AmI = A - I s = s + 1 end # Compute accurate superdiagonal of T - p = (1/2^s) + p = 1 / 2^s for k = 1:n-1 Ak = A0[k,k] Akp1 = A0[k+1,k+1] @@ -999,7 +1000,7 @@ function logm{T}(A0::UpperTriangular{T}) A[k+1,k+1] = Akp1p if Ak == Akp1 A[k,k+1] = p * A0[k,k+1] * Ak^(p-1) - elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2 + 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) @@ -1017,7 +1018,7 @@ function logm{T}(A0::UpperTriangular{T}) r = a - 1 end s0 = s - if angle(a) >= pi/2 + if angle(a) >= pi / 2 a = sqrt(a) s0 = s - 1 end @@ -1034,29 +1035,20 @@ function logm{T}(A0::UpperTriangular{T}) # 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,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 + w[i] = V[1,i]^2 end # Compute the Padé approximation - B = zeros(T, n, n) Y = zeros(T, n, n) for k = 1:m - for j = 1:n - B[j,j] = 1. + A[j,j] * x[k] - for i = 1:j-1 - B[i,j] = A[i,j] * x[k] - end - end - B = A * inv(B) - scale!(w[k], B) - Y = Y + B + Y = Y + w[k] * (A / (x[k] * A + I)) end # Scale back @@ -1072,7 +1064,7 @@ function logm{T}(A0::UpperTriangular{T}) Y[k+1,k+1] = logAkp1 if Ak == Akp1 Y[k,k+1] = A0[k,k+1] / Ak - elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2 + 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)))) diff --git a/test/linalg3.jl b/test/linalg3.jl index becc3b5daca1d..d17b10799e701 100644 --- a/test/linalg3.jl +++ b/test/linalg3.jl @@ -154,7 +154,15 @@ for elty in (Float32, Float64, Complex64, Complex128) 0.135335281175235 0.406005843524598 0.541341126763207]') @test_approx_eq expm(A3) eA3 - A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5; + # 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]) @@ -168,12 +176,6 @@ for elty in (Float32, Float64, Complex64, Complex128) 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 - - # 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 A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1];