Make the algorithm for real powers of a matrix robust

Author: mfasi

base/linalg/dense.jl

@@ -169,17 +169,50 @@ kron(a::Vector, b::Vector)=vec(kron(reshape(a,length(a),1),reshape(b,length(b),1
 kron(a::Matrix, b::Vector)=kron(a,reshape(b,length(b),1))
 kron(a::Vector, b::Matrix)=kron(reshape(a,length(a),1),b)
 
-^(A::Matrix, p::Integer) = p < 0 ? inv(A^-p) : Base.power_by_squaring(A,p)
-
-function ^(A::Matrix, p::Number)
+# Matrix power
+^(A::AbstractMatrix, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A),-p) : Base.power_by_squaring(A,p)
+function ^(A::AbstractMatrix, p::Real)
+    # For integer powers, use repeated squaring
     if isinteger(p)
         return A^Integer(real(p))
     end
-    chksquare(A)
-    v, X = eig(A)
-    any(v.<0) && (v = complex(v))
-    Xinv = ishermitian(A) ? X' : inv(X)
-    scale(X, v.^p)*Xinv
+
+    # If possible, use diagonalization
+    if issym(A)
+        return full(Symmetric(A)^p)
+    end
+    if ishermitian(A)
+        return full(Hermitian(A)^p)
+    end
+
+    # Otherwise, use Schur decomposition
+    n = chksquare(A)
+    if istriu(A)
+        retmat = full(UpperTriangular(A)^p)
+        d = diag(A)
+    else
+        S,Q,d = schur(complex(A))
+        R = UpperTriangular(S)^p
+        retmat = Q * R * Q'
+    end
+
+    # 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 nonprincipal matrix power will be returned.")
+    end
+
+    if isreal(A) && ~np_real_eigs
+        return real(retmat)
+    else
+        return retmat
+    end
 end
 
 # Matrix exponential
@@ -286,6 +319,9 @@ expm(x::Number) = exp(x)
 
 function logm(A::StridedMatrix)
     # If possible, use diagonalization
+    if issym(A)
+        return full(logm(Symmetric(A)))
+    end
     if ishermitian(A)
         return full(logm(Hermitian(A)))
     end

base/linalg/diagonal.jl

@@ -118,9 +118,14 @@ end
 # identity matrices via eye(Diagonal{type},n)
 eye{T}(::Type{Diagonal{T}}, n::Int) = Diagonal(ones(T,n))
 
-expm(D::Diagonal) = Diagonal(exp(D.diag))
-logm(D::Diagonal) = Diagonal(log(D.diag))
-sqrtm(D::Diagonal) = Diagonal(sqrt(D.diag))
+# Matrix functions
+^(D::Diagonal, p::Integer) = Diagonal((D.diag).^p)
+^(D::Diagonal, p::Real) = Diagonal((D.diag).^p)
+for (funm, func) in ([:expm,:exp], [:sqrtm,:sqrt], [:logm,:log])
+    @eval begin
+        ($funm)(D::Diagonal) = Diagonal(($func)(D.diag))
+    end
+end
 
 #Linear solver
 function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)

base/linalg/symmetric.jl

@@ -131,7 +131,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 ^(A::Symmetric, 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 ^(A::Symmetric, p::Real)
+    F = eigfact(full(A))
+    if isposdef(F)
+        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 = chksquare(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 = chksquare(A)
+    F = eigfact(A)
+    if isposdef(F)
+        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(full(A))
     return Symmetric((F.vectors * Diagonal(exp(F.values))) * F.vectors')
@@ -151,9 +198,7 @@ function expm{T}(A::Hermitian{T})
 end
 
 for (funm, func) in ([:logm,:log], [:sqrtm,:sqrt])
-
     @eval begin
-
         function ($funm)(A::Symmetric)
             F = eigfact(full(A))
             if isposdef(F)
@@ -182,7 +227,5 @@ for (funm, func) in ([:logm,:log], [:sqrtm,:sqrt])
                 return retmat
             end
         end
-
     end
-
 end