Skip to content

Commit 4d9bffa

Browse files
mfasimfasi
committed
Matrix logarithm function logm added
1 parent c75f4d4 commit 4d9bffa

File tree

7 files changed

+272
-20
lines changed

7 files changed

+272
-20
lines changed

‎base/exports.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -659,6 +659,7 @@ export
659659
linreg,
660660
logabsdet,
661661
logdet,
662+
logm,
662663
lu,
663664
lufact!,
664665
lufact,

‎base/linalg.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,7 @@ export
8585
linreg,
8686
logabsdet,
8787
logdet,
88+
logm,
8889
lu,
8990
lufact,
9091
lufact!,

‎base/linalg/dense.jl

Lines changed: 35 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -278,9 +278,42 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
278278
end
279279
end
280280

281+
function logm{T}(A::StridedMatrix{T})
282+
# If possible, use diagonalization
283+
if ishermitian(A)
284+
return funm(log, Hermitian(A))
285+
end
286+
287+
# Use Schur decomposition
288+
n = chksquare(A)
289+
S,Q,d = schur(complex(A))
290+
R = full(logm(UpperTriangular(S)))
291+
retmat = Q*R*Q'
292+
293+
# Check whether the matrix has nonpositive real eigs
294+
np_real_eigs = false
295+
for i = 1:n
296+
if imag(d[i]) < eps() && real(d[i]) <= 0
297+
np_real_eigs = true
298+
break
299+
end
300+
end
301+
if np_real_eigs
302+
warn("Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.")
303+
end
304+
305+
if isreal(A) && ~np_real_eigs
306+
return real(retmat)
307+
else
308+
return retmat
309+
end
310+
end
311+
logm(a::Number) = (b = log(complex(a)); imag(b) == 0 ? real(b) : b)
312+
logm(a::Complex) = log(a)
313+
281314
function sqrtm{T<:Real}(A::StridedMatrix{T})
282315
if issym(A)
283-
return sqrtm(Symmetric(A))
316+
return funm(sqrt, Symmetric(A))
284317
end
285318
n = chksquare(A)
286319
SchurF = schurfact(complex(A))
@@ -290,7 +323,7 @@ function sqrtm{T<:Real}(A::StridedMatrix{T})
290323
end
291324
function sqrtm{T<:Complex}(A::StridedMatrix{T})
292325
if ishermitian(A)
293-
return sqrtm(Hermitian(A))
326+
return funm(sqrt, Hermitian(A))
294327
end
295328
n = chksquare(A)
296329
SchurF = schurfact(A)

‎base/linalg/symmetric.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -113,8 +113,8 @@ end
113113

114114
#Matrix-valued functions
115115
expm{T<:Real}(A::RealHermSymComplexHerm{T}) = (F = eigfact(A); F.vectors*Diagonal(exp(F.values))*F.vectors')
116-
function sqrtm{T<:Real}(A::RealHermSymComplexHerm{T})
116+
function funm{T<:Real}(f::Function, A::RealHermSymComplexHerm{T})
117117
F = eigfact(A)
118-
isposdef(F) && return F.vectors*Diagonal(sqrt(F.values))*F.vectors'
119-
return F.vectors*Diagonal(sqrt(complex(F.values)))*F.vectors'
118+
isposdef(F) && return F.vectors*Diagonal(f(F.values))*F.vectors'
119+
return F.vectors*Diagonal(f(complex(F.values)))*F.vectors'
120120
end

‎base/linalg/triangular.jl

Lines changed: 195 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -889,6 +889,201 @@ for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv
889889
end
890890
end
891891

892+
# Complex matrix logarithm for the upper triangular factor, see:
893+
# Al-Mohy and Higham, "Improved inverse scaling and squaring algorithms for
894+
# the matrix logarithm", SIAM J. Sci. Comput., 34(4), (2012), pp. C153-C169.
895+
# Al-Mohy, Higham and Relton, "Computing the Frechet derivative of the matrix
896+
# logarithm and estimating the condition number", SIAM J. Sci. Comput., 35(4),
897+
# (2013), C394-C410.
898+
function logm{T}(A0::UpperTriangular{T})
899+
maxsqrt = 100
900+
theta = [1.586970738772063e-005,
901+
2.313807884242979e-003,
902+
1.938179313533253e-002,
903+
6.209171588994762e-002,
904+
1.276404810806775e-001,
905+
2.060962623452836e-001,
906+
2.879093714241194e-001]
907+
tmax = size(theta, 1)
908+
n = size(A0, 1)
909+
A = A0
910+
p = 0
911+
m = 0
912+
913+
# Compute repeated roots
914+
d = diag(A)
915+
dm1 = Array(T, n)
916+
s = 0
917+
for i = 1:n
918+
dm1[i] = dm1[i] - 1.
919+
end
920+
while norm(dm1, Inf) > theta[tmax]
921+
for i = 1:n
922+
d[i] = sqrt(d[i])
923+
end
924+
for i = 1:n
925+
dm1[i] = d[i] - 1
926+
end
927+
s = s + 1
928+
end
929+
s0 = s
930+
for k = 1:min(s, maxsqrt)
931+
A = sqrtm(A)
932+
end
933+
934+
AmI = A - UpperTriangular(eye(n))
935+
d2 = norm(AmI^2, 1)^(1/2)
936+
d3 = norm(AmI^3, 1)^(1/3)
937+
alpha2 = max(d2, d3)
938+
foundm = false
939+
if alpha2 <= theta[2]
940+
m = alpha2<=theta[1]?1:2
941+
foundm = true
942+
end
943+
944+
while ~foundm
945+
more = false
946+
if s > s0
947+
d3 = norm(AmI^3, 1)^(1/3)
948+
end
949+
d4 = norm(AmI^4, 1)^(1/4)
950+
alpha3 = max(d3, d4)
951+
if alpha3 <= theta[tmax]
952+
for j = 3:tmax
953+
if alpha3 <= theta[j]
954+
break
955+
end
956+
end
957+
if j <= 6
958+
m = j
959+
break
960+
elseif alpha3/2 <= theta[5] && p < 2
961+
more = true
962+
p = p + 1
963+
end
964+
end
965+
966+
if ~more
967+
d5 = norm(AmI^5, 1)^(1/5)
968+
alpha4 = max(d4, d5);
969+
eta = min(alpha3, alpha4);
970+
if eta <= theta[tmax]
971+
j = 0
972+
for j = 6:tmax
973+
if eta <= theta[j]
974+
m = j
975+
break
976+
end
977+
end
978+
break
979+
end
980+
end
981+
982+
if s == maxsqrt
983+
m = tmax
984+
break
985+
end
986+
A = sqrtm(A)
987+
AmI = A - UpperTriangular(eye(n))
988+
s = s + 1
989+
end
990+
991+
# Compute accurate superdiagonal of T
992+
p = (1/2^s)
993+
for k = 1:n-1
994+
Ak = A0[k,k]
995+
Akp1 = A0[k+1,k+1]
996+
Akp = Ak^p
997+
Akp1p = Akp1^p
998+
A[k,k] = Akp
999+
A[k+1,k+1] = Akp1p
1000+
if Ak == Akp1
1001+
A[k,k+1] = p * A0[k,k+1] * Ak^(p-1)
1002+
elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2
1003+
A[k,k+1] = A0[k,k+1] * (Akp1p - Akp) / (Akp1 - Ak)
1004+
else
1005+
logAk = log(Ak)
1006+
logAkp1 = log(Akp1)
1007+
w = atanh((Akp1 - Ak)/(Akp1 + Ak)) + im*pi*ceil((imag(logAkp1-logAk)-pi)/(2*pi))
1008+
dd = 2 * exp(p*(logAk+logAkp1)/2) * sinh(p*w) / (Akp1 - Ak);
1009+
A[k,k+1] = A0[k,k+1] * dd
1010+
end
1011+
end
1012+
1013+
# Compute accurate diagonal of T
1014+
for i = 1:n
1015+
a = A0[i,i]
1016+
if s == 0
1017+
r = a - 1
1018+
end
1019+
s0 = s
1020+
if angle(a) >= pi/2
1021+
a = sqrt(a)
1022+
s0 = s - 1
1023+
end
1024+
z0 = a - 1
1025+
a = sqrt(a)
1026+
r = 1 + a
1027+
for j = 1:s0-1
1028+
a = sqrt(a)
1029+
r = r * (1 + a)
1030+
end
1031+
A[i,i] = z0 / r
1032+
end
1033+
1034+
# Compute the Gauss-Legendre quadrature formula
1035+
R = zeros(Float64, m, m)
1036+
for i = 1:m - 1
1037+
R[i,i+1] = i / sqrt((2 * i).^2-1)
1038+
R[i+1,i] = R[i,i+1]
1039+
end
1040+
x,V = eig(R)
1041+
w = Array(Float64, m)
1042+
for i = 1:m
1043+
x[i] = (x[i] + 1) / 2
1044+
w[i] = (V[1,i])^2
1045+
end
1046+
1047+
# Compute the Padé approximation
1048+
B = zeros(T, n, n)
1049+
Y = zeros(T, n, n)
1050+
for k = 1:m
1051+
for j = 1:n
1052+
B[j,j] = 1. + A[j,j] * x[k]
1053+
for i = 1:j-1
1054+
B[i,j] = A[i,j] * x[k]
1055+
end
1056+
end
1057+
B = A * inv(B)
1058+
scale!(w[k], B)
1059+
Y = Y + B
1060+
end
1061+
1062+
# Scale back
1063+
scale!(2^s, Y)
1064+
1065+
# Compute accurate diagonal and superdiagonal of log(T)
1066+
for k = 1:n-1
1067+
Ak = A0[k,k]
1068+
Akp1 = A0[k+1,k+1]
1069+
logAk = log(Ak)
1070+
logAkp1 = log(Akp1)
1071+
Y[k,k] = logAk
1072+
Y[k+1,k+1] = logAkp1
1073+
if Ak == Akp1
1074+
Y[k,k+1] = A0[k,k+1] / Ak
1075+
elseif abs(Ak) < abs(Akp1)/2 || abs(Akp1) < abs(Ak)/2
1076+
Y[k,k+1] = A0[k,k+1] * (logAkp1 - logAk) / (Akp1 - Ak)
1077+
else
1078+
w = atanh((Akp1 - Ak)/(Akp1 + Ak) + im*pi*(ceil((imag(logAkp1-logAk) - pi)/(2*pi))))
1079+
Y[k,k+1] = 2 * A0[k,k+1] * w / (Akp1 - Ak)
1080+
end
1081+
end
1082+
1083+
return UpperTriangular(Y)
1084+
1085+
end
1086+
8921087
function sqrtm{T}(A::UpperTriangular{T})
8931088
n = size(A, 1)
8941089
TT = typeof(sqrt(zero(T)))

‎doc/stdlib/linalg.rst

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -646,6 +646,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
646646

647647
Matrix exponential.
648648

649+
.. function:: logm(A)
650+
651+
Matrix logarithm.
652+
649653
.. function:: lyap(A, C)
650654

651655
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

Lines changed: 33 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -154,21 +154,20 @@ for elty in (Float32, Float64, Complex64, Complex128)
154154
0.135335281175235 0.406005843524598 0.541341126763207]')
155155
@test_approx_eq expm(A3) eA3
156156

157-
# issue 5116
158-
A4 = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
159-
eA4 = [-0.999786072879326 -0.065407069689389 0.0 0.0
160-
0.006540706968939 -0.999786072879326 0.0 0.0
161-
0.0 0.0 1.0 0.0
162-
0.013081413937878 -3.999572145758650 0.0 1.0]
163-
@test_approx_eq expm(A4) eA4
164-
165-
# issue 5116
166-
A5 = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
167-
eA5 = [ 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
168-
0.0+0.0im 1.543080634815244+0.0im 0.0-1.175201193643801im 0.0+0.0im
169-
0.0+0.0im 0.0+1.175201193643801im 1.543080634815243+0.0im 0.0+0.0im
170-
0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im]
171-
@test_approx_eq expm(A5) eA5
157+
A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5;
158+
1/3 1/4 1/5 1/6;
159+
1/4 1/5 1/6 1/7;
160+
1/5 1/6 1/7 1/8])
161+
@test_approx_eq expm(logm(A4)) A4
162+
163+
A5 = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
164+
@test_approx_eq expm(logm(A5)) A5
165+
166+
A6 = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
167+
@test_approx_eq expm(logm(A6)) A6
168+
169+
A7 = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
170+
@test_approx_eq expm(logm(A7)) A7
172171

173172
# Hessenberg
174173
@test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
@@ -177,6 +176,25 @@ for elty in (Float32, Float64, Complex64, Complex128)
177176
0 -0.000000000000002 3.000000000000000])
178177
end
179178

179+
A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1];
180+
@test_approx_eq expm(logm(A8)) A8
181+
182+
# issue 5116
183+
A9 = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
184+
eA9 = [-0.999786072879326 -0.065407069689389 0.0 0.0
185+
0.006540706968939 -0.999786072879326 0.0 0.0
186+
0.0 0.0 1.0 0.0
187+
0.013081413937878 -3.999572145758650 0.0 1.0]
188+
@test_approx_eq expm(A9) eA9
189+
190+
# issue 5116
191+
A10 = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
192+
eA10 = [ 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
193+
0.0+0.0im 1.543080634815244+0.0im 0.0-1.175201193643801im 0.0+0.0im
194+
0.0+0.0im 0.0+1.175201193643801im 1.543080634815243+0.0im 0.0+0.0im
195+
0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im]
196+
@test_approx_eq expm(A10) eA10
197+
180198
# matmul for types w/o sizeof (issue #1282)
181199
A = Array(Complex{Int},10,10)
182200
A[:] = complex(1,1)

0 commit comments

Comments
 (0)