diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index d70f7c2549f08..70222f62844c0 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -29,7 +29,7 @@ function scale!(X::Array{T}, s::Real) where T<:BlasComplex X end -#Test whether a matrix is positive-definite +# Test whether a matrix is positive-definite isposdef!(A::StridedMatrix{<:BlasFloat}, UL::Symbol) = LAPACK.potrf!(char_uplo(UL), A)[2] == 0 """ @@ -323,18 +323,67 @@ kron(a::AbstractVector, b::AbstractVector)=vec(kron(reshape(a,length(a),1),resha kron(a::AbstractMatrix, b::AbstractVector)=kron(a,reshape(b,length(b),1)) kron(a::AbstractVector, b::AbstractMatrix)=kron(reshape(a,length(a),1),b) -^(A::AbstractMatrix, p::Integer) = p < 0 ? inv(A^-p) : Base.power_by_squaring(A,p) - -function ^(A::AbstractMatrix, p::Number) +# Matrix power +^{T}(A::AbstractMatrix{T}, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A), -p) : Base.power_by_squaring(A, p) +function ^{T}(A::AbstractMatrix{T}, p::Real) + # For integer powers, use repeated squaring if isinteger(p) - return A^Integer(real(p)) + TT = Base.promote_op(^, eltype(A), typeof(p)) + return (TT == eltype(A) ? A : copy!(similar(A, TT), A))^Integer(p) + end + + # If possible, use diagonalization + if T <: Real && issymmetric(A) + return (Symmetric(A)^p) + end + if ishermitian(A) + return (Hermitian(A)^p) + end + + n = checksquare(A) + + # Quicker return if A is diagonal + if isdiag(A) + retmat = copy(A) + for i in 1:n + retmat[i, i] = retmat[i, i] ^ p + end + return retmat + end + + # Otherwise, use Schur decomposition + if istriu(A) + # Integer part + retmat = A ^ floor(p) + # Real part + if p - floor(p) == 0.5 + # special case: A^0.5 === sqrtm(A) + retmat = retmat * sqrtm(A) + else + retmat = retmat * powm!(UpperTriangular(float.(A)), real(p - floor(p))) + end + else + S,Q,d = schur(complex(A)) + # Integer part + R = S ^ floor(p) + # Real part + if p - floor(p) == 0.5 + # special case: A^0.5 === sqrtm(A) + R = R * sqrtm(S) + else + R = R * powm!(UpperTriangular(float.(S)), real(p - floor(p))) + end + retmat = Q * R * Q' + end + + # if A has nonpositive real eigenvalues, retmat is a nonprincipal matrix power. + if isreal(retmat) + return real(retmat) + else + return retmat end - checksquare(A) - v, X = eig(A) - any(v.<0) && (v = complex(v)) - Xinv = ishermitian(A) ? X' : inv(X) - (X * Diagonal(v.^p)) * Xinv end +^(A::AbstractMatrix, p::Number) = expm(p*logm(A)) # Matrix exponential @@ -466,7 +515,7 @@ function rcswap!(i::Integer, j::Integer, X::StridedMatrix{<:Number}) end """ - logm(A::StridedMatrix) + logm(A{T}::StridedMatrix{T}) If `A` has no negative real eigenvalue, compute the principal matrix logarithm of `A`, i.e. the unique matrix ``X`` such that ``e^X = A`` and ``-\\pi < Im(\\lambda) < \\pi`` for all @@ -497,8 +546,11 @@ julia> logm(A) 0.0 1.0 ``` """ -function logm(A::StridedMatrix) +function logm{T}(A::StridedMatrix{T}) # If possible, use diagonalization + if issymmetric(A) && T <: Real + return full(logm(Symmetric(A))) + end if ishermitian(A) return full(logm(Hermitian(A))) end diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index 34e1d632d6df6..9c8647fe2c22f 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -288,6 +288,7 @@ end # identity matrices via eye(Diagonal{type},n) eye(::Type{Diagonal{T}}, n::Int) where {T} = Diagonal(ones(T,n)) +# Matrix functions expm(D::Diagonal) = Diagonal(exp.(D.diag)) expm(D::Diagonal{<:AbstractMatrix}) = Diagonal(expm.(D.diag)) logm(D::Diagonal) = Diagonal(log.(D.diag)) diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index b4004f008536e..d92d2c39d8baa 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -262,7 +262,6 @@ include("hessenberg.jl") include("lq.jl") include("eigen.jl") include("svd.jl") -include("schur.jl") include("symmetric.jl") include("cholesky.jl") include("lu.jl") @@ -274,6 +273,8 @@ include("givens.jl") include("special.jl") include("bitarray.jl") include("ldlt.jl") +include("schur.jl") + include("arpack.jl") include("arnoldi.jl") diff --git a/base/linalg/schur.jl b/base/linalg/schur.jl index ad6cd33b800cb..da63220b5128a 100644 --- a/base/linalg/schur.jl +++ b/base/linalg/schur.jl @@ -107,6 +107,12 @@ function schur(A::StridedMatrix) SchurF = schurfact(A) SchurF[:T], SchurF[:Z], SchurF[:values] end +schur(A::Symmetric) = schur(full(A)) +schur(A::Hermitian) = schur(full(A)) +schur(A::UpperTriangular) = schur(full(A)) +schur(A::LowerTriangular) = schur(full(A)) +schur(A::Tridiagonal) = schur(full(A)) + """ ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index f0163356afc15..f417ce2ead09f 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -1,6 +1,6 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -#Symmetric and Hermitian matrices +# Symmetric and Hermitian matrices struct Symmetric{T,S<:AbstractMatrix} <: AbstractMatrix{T} data::S uplo::Char @@ -181,7 +181,7 @@ trace(A::Hermitian) = real(trace(A.data)) Base.conj(A::HermOrSym) = typeof(A)(conj(A.data), A.uplo) Base.conj!(A::HermOrSym) = typeof(A)(conj!(A.data), A.uplo) -#tril/triu +# tril/triu function tril(A::Hermitian, k::Integer=0) if A.uplo == 'U' && k <= 0 return tril!(A.data',k) @@ -235,7 +235,7 @@ end ## Matvec A_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::Symmetric{T,<:StridedMatrix}, x::StridedVector{T}) = BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y) A_mul_B!{T<:BlasComplex}(y::StridedVector{T}, A::Hermitian{T,<:StridedMatrix}, x::StridedVector{T}) = BLAS.hemv!(A.uplo, one(T), A.data, x, zero(T), y) -##Matmat +## Matmat A_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::Symmetric{T,<:StridedMatrix}, B::StridedMatrix{T}) = BLAS.symm!('L', A.uplo, one(T), A.data, B, zero(T), C) A_mul_B!{T<:BlasFloat}(C::StridedMatrix{T}, A::StridedMatrix{T}, B::Symmetric{T,<:StridedMatrix}) = BLAS.symm!('R', B.uplo, one(T), B.data, A, zero(T), C) A_mul_B!{T<:BlasComplex}(C::StridedMatrix{T}, A::Hermitian{T,<:StridedMatrix}, B::StridedMatrix{T}) = BLAS.hemm!('L', A.uplo, one(T), A.data, B, zero(T), C) @@ -403,7 +403,54 @@ function svdvals!{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{ return sort!(vals, rev = true) end -#Matrix-valued functions +# Matrix functions +function ^{T<:Real}(A::Symmetric{T}, p::Integer) + if p < 0 + return Symmetric(Base.power_by_squaring(inv(A), -p)) + else + return Symmetric(Base.power_by_squaring(A, p)) + end +end +function ^{T<:Real}(A::Symmetric{T}, p::Real) + F = eigfact(A) + if all(λ -> λ ≥ 0, F.values) + retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors' + else + retmat = (F.vectors * Diagonal((complex(F.values)).^p)) * F.vectors' + end + return Symmetric(retmat) +end +function ^(A::Hermitian, p::Integer) + n = checksquare(A) + if p < 0 + retmat = Base.power_by_squaring(inv(A), -p) + else + retmat = Base.power_by_squaring(A, p) + end + for i = 1:n + retmat[i,i] = real(retmat[i,i]) + end + return Hermitian(retmat) +end +function ^{T}(A::Hermitian{T}, p::Real) + n = checksquare(A) + F = eigfact(A) + if all(λ -> λ ≥ 0, F.values) + retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors' + if T <: Real + return Hermitian(retmat) + else + for i = 1:n + retmat[i,i] = real(retmat[i,i]) + end + return Hermitian(retmat) + end + else + retmat = (F.vectors * Diagonal((complex(F.values).^p))) * F.vectors' + return retmat + end +end + function expm(A::Symmetric) F = eigfact(A) return Symmetric((F.vectors * Diagonal(exp.(F.values))) * F.vectors') @@ -423,10 +470,8 @@ function expm{T}(A::Hermitian{T}) end for (funm, func) in ([:logm,:log], [:sqrtm,:sqrt]) - @eval begin - - function ($funm)(A::Symmetric) + function ($funm){T<:Real}(A::Symmetric{T}) F = eigfact(A) if isposdef(F) retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors' @@ -454,7 +499,5 @@ for (funm, func) in ([:logm,:log], [:sqrtm,:sqrt]) return retmat end end - end - end diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index cd517d1985c22..e940ff2c138e8 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1662,13 +1662,81 @@ At_ldiv_B(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw Ac_ldiv_B(::Union{UpperTriangular,LowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) Ac_ldiv_B(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +# Complex matrix power for upper triangular factor, see: +# Higham and Lin, "A Schur-Padé algorithm for fractional powers of a Matrix", +# SIAM J. Matrix Anal. & Appl., 32 (3), (2011) 1056–1078. +# Higham and Lin, "An improved Schur-Padé algorithm for fractional powers of +# a matrix and their Fréchet derivatives", SIAM. J. Matrix Anal. & Appl., +# 34(3), (2013) 1341–1360. +function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real) + if abs(p) >= 1 + ArgumentError("p must be a real number in (-1,1), got $p") + end + + normA0 = norm(A0, 1) + scale!(A0, 1/normA0) + + theta = [1.53e-5, 2.25e-3, 1.92e-2, 6.08e-2, 1.25e-1, 2.03e-1, 2.84e-1] + n = checksquare(A0) + + A, m, s = invsquaring(A0, theta) + A = I - A + + # Compute accurate diagonal of I - T + sqrt_diag!(A0, A, s) + for i = 1:n + A[i, i] = -A[i, i] + end + # Compute the Padé approximant + c = 0.5 * (p - m) / (2 * m - 1) + triu!(A) + S = c * A + Stmp = similar(S) + for j = m-1:-1:1 + j4 = 4 * j + c = (-p - j) / (j4 + 2) + for i = 1:n + @inbounds S[i, i] = S[i, i] + 1 + end + copy!(Stmp, S) + scale!(S, A, c) + A_ldiv_B!(Stmp, S.data) + + c = (p - j) / (j4 - 2) + for i = 1:n + @inbounds S[i, i] = S[i, i] + 1 + end + copy!(Stmp, S) + scale!(S, A, c) + A_ldiv_B!(Stmp, S.data) + end + for i = 1:n + S[i, i] = S[i, i] + 1 + end + copy!(Stmp, S) + scale!(S, A, -p) + A_ldiv_B!(Stmp, S.data) + for i = 1:n + @inbounds S[i, i] = S[i, i] + 1 + end + + blockpower!(A0, S, p/(2^s)) + for m = 1:s + A_mul_B!(Stmp.data, S, S) + copy!(S, Stmp) + blockpower!(A0, S, p/(2^(s-m))) + end + scale!(S, normA0^p) + return S +end +powm(A::LowerTriangular, p::Real) = powm(A.', p::Real).' # 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. +# 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. +# logarithm and estimating the condition number", SIAM J. Sci. Comput., +# 35(4), (2013), C394–C410. # # Based on the code available at http://eprints.ma.man.ac.uk/1851/02/logm.zip, # Copyright (c) 2011, Awad H. Al-Mohy and Nicholas J. Higham @@ -1809,7 +1877,7 @@ function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T}) A[i,i] = z0 / r end - # Compute the Gauss-Legendre quadrature formula + # Get the Gauss-Legendre quadrature points and weights R = zeros(Float64, m, m) for i = 1:m - 1 R[i,i+1] = i / sqrt((2 * i)^2 - 1) @@ -1850,10 +1918,169 @@ function logm{T<:Union{Float64,Complex{Float64}}}(A0::UpperTriangular{T}) end return UpperTriangular(Y) - end logm(A::LowerTriangular) = logm(A.').' +# Auxiliary functions for logm and matrix power + +# Compute accurate diagonal of A = A0^s - I +# Al-Mohy, "A more accurate Briggs method for the logarithm", +# Numer. Algorithms, 59, (2012), 393–402. +function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s) + n = checksquare(A0) + @inbounds for i = 1:n + a = complex(A0[i,i]) + if s == 0 + A[i,i] = a - 1 + else + s0 = s + if imag(a) >= 0 && real(a) <= 0 && a != 0 + 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 + end +end + +# Used only by powm at the moment +# Repeatedly compute the square roots of A so that in the end its +# eigenvalues are close enough to the positive real line +function invsquaring(A0::UpperTriangular, theta) + # assumes theta is in ascending order + maxsqrt = 100 + tmax = size(theta, 1) + n = checksquare(A0) + A = complex(copy(A0)) + p = 0 + m = 0 + + # Compute repeated roots + d = complex(diag(A)) + dm1 = similar(d, n) + s = 0 + for i = 1:n + dm1[i] = d[i] - 1 + end + while norm(dm1, Inf) > theta[tmax] + for i = 1:n + d[i] = sqrt(d[i]) + 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 + elseif alpha3 / 2 <= theta[5] && p < 2 + more = true + p = p + 1 + end + end + if j <= 6 + m = j + foundm = true + 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 + break + end + end + if s == maxsqrt + m = tmax + break + end + A = sqrtm(A) + AmI = A - I + s = s + 1 + end + end + + # Compute accurate superdiagonal of T + p = 1 / 2^s + A = complex(A) + blockpower!(A, A0, p) + return A,m,s +end + +# Compute accurate diagonal and superdiagonal of A = A0^p +function blockpower!(A::UpperTriangular, A0::UpperTriangular, p) + n = checksquare(A0) + @inbounds for k = 1:n-1 + Ak = complex(A0[k,k]) + Akp1 = complex(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 * unw(logAkp1-logAk) + dd = 2 * exp(p*(logAk+logAkp1)/2) * sinh(p*w) / (Akp1 - Ak); + A[k,k+1] = A0[k,k+1] * dd + end + end +end + +# Unwinding number +unw(x::Real) = 0 +unw(x::Number) = ceil((imag(x) - pi) / (2 * pi)) + +# End of auxiliary functions for logm and matrix power + function sqrtm(A::UpperTriangular) realmatrix = false if isreal(A) @@ -1907,7 +2134,7 @@ end sqrtm(A::LowerTriangular) = sqrtm(A.').' sqrtm(A::UnitLowerTriangular) = sqrtm(A.').' -#Generic eigensystems +# Generic eigensystems eigvals(A::AbstractTriangular) = diag(A) function eigvecs{T}(A::AbstractTriangular{T}) TT = promote_type(T, Float32) @@ -1938,7 +2165,7 @@ end eigfact(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A)) -#Generic singular systems +# Generic singular systems for func in (:svd, :svdfact, :svdfact!, :svdvals) @eval begin ($func)(A::AbstractTriangular) = ($func)(full(A)) diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index fdf107d743469..29a7659caad0f 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -512,6 +512,61 @@ end @test_throws ArgumentError diag(zeros(0,1),2) end +@testset "Matrix to real power" for elty in (Float64, Complex{Float64}) +# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014 + #Aa : only positive real eigenvalues + Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2]) + + @test Aa^(1/2) ≈ sqrtm(Aa) + @test Aa^(-1/2) ≈ inv(sqrtm(Aa)) + @test Aa^(3/4) ≈ sqrtm(Aa) * sqrtm(sqrtm(Aa)) + @test Aa^(-3/4) ≈ inv(Aa) * sqrtm(sqrtm(Aa)) + @test Aa^(17/8) ≈ Aa^2 * sqrtm(sqrtm(sqrtm(Aa))) + @test Aa^(-17/8) ≈ inv(Aa^2 * sqrtm(sqrtm(sqrtm(Aa)))) + @test (Aa^0.2)^5 ≈ Aa + @test (Aa^(2/3))*(Aa^(1/3)) ≈ Aa + @test (Aa^im)^(-im) ≈ Aa + + #Ab : both positive and negative real eigenvalues + Ab = convert(Matrix{elty}, [1 2 3; 4 7 1; 2 1 4]) + + @test Ab^(1/2) ≈ sqrtm(Ab) + @test Ab^(-1/2) ≈ inv(sqrtm(Ab)) + @test Ab^(3/4) ≈ sqrtm(Ab) * sqrtm(sqrtm(Ab)) + @test Ab^(-3/4) ≈ inv(Ab) * sqrtm(sqrtm(Ab)) + @test Ab^(17/8) ≈ Ab^2 * sqrtm(sqrtm(sqrtm(Ab))) + @test Ab^(-17/8) ≈ inv(Ab^2 * sqrtm(sqrtm(sqrtm(Ab)))) + @test (Ab^0.2)^5 ≈ Ab + @test (Ab^(2/3))*(Ab^(1/3)) ≈ Ab + @test (Ab^im)^(-im) ≈ Ab + + #Ac : complex eigenvalues + Ac = convert(Matrix{elty}, [5 4 2 1;0 1 -1 -1;-1 -1 3 6;1 1 -1 5]) + + @test Ac^(1/2) ≈ sqrtm(Ac) + @test Ac^(-1/2) ≈ inv(sqrtm(Ac)) + @test Ac^(3/4) ≈ sqrtm(Ac) * sqrtm(sqrtm(Ac)) + @test Ac^(-3/4) ≈ inv(Ac) * sqrtm(sqrtm(Ac)) + @test Ac^(17/8) ≈ Ac^2 * sqrtm(sqrtm(sqrtm(Ac))) + @test Ac^(-17/8) ≈ inv(Ac^2 * sqrtm(sqrtm(sqrtm(Ac)))) + @test (Ac^0.2)^5 ≈ Ac + @test (Ac^(2/3))*(Ac^(1/3)) ≈ Ac + @test (Ac^im)^(-im) ≈ Ac + + #Ad : defective Matrix + Ad = convert(Matrix{elty}, [3 1; 0 3]) + + @test Ad^(1/2) ≈ sqrtm(Ad) + @test Ad^(-1/2) ≈ inv(sqrtm(Ad)) + @test Ad^(3/4) ≈ sqrtm(Ad) * sqrtm(sqrtm(Ad)) + @test Ad^(-3/4) ≈ inv(Ad) * sqrtm(sqrtm(Ad)) + @test Ad^(17/8) ≈ Ad^2 * sqrtm(sqrtm(sqrtm(Ad))) + @test Ad^(-17/8) ≈ inv(Ad^2 * sqrtm(sqrtm(sqrtm(Ad)))) + @test (Ad^0.2)^5 ≈ Ad + @test (Ad^(2/3))*(Ad^(1/3)) ≈ Ad + @test (Ad^im)^(-im) ≈ Ad +end + @testset "Least squares solutions" begin a = [ones(20) 1:20 1:20] b = reshape(eye(8, 5), 20, 2)