reverted powm export. Small fix if Upper-triangular matrix contains I…

…ntegers.

added a few more test for ^;now uses random 10x10 matrix as well.

fixed some ambiguities, corrected the algorithm a bit.

correctly promotes Matrix{Int}^Float64

reverted some wrong changes

reverted logm changes

speedup if A is diagonal, fixed some tests

small changes based on feedback

powm is now scale-invariant

powm scales in-place, fixed testing a bit, corrected some bugs

removed full() calls

added rollback for Matrix^Complex

added schur() support for Symmetric,Hermitian,Triangular,Tridiagonal matrices

base/linalg/dense.jl

@@ -324,53 +324,66 @@ kron(a::AbstractMatrix, b::AbstractVector)=kron(a,reshape(b,length(b),1))
 kron(a::AbstractVector, b::AbstractMatrix)=kron(reshape(a,length(a),1),b)
 
 # Matrix power
-^(A::AbstractMatrix, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A),-p) : Base.power_by_squaring(A,p)
+^{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 issymmetric(A) && T <: Real
- return full(Symmetric(A)^p)
+ if T <: Real && issymmetric(A)
+ return (Symmetric(A)^p)
  end
  if ishermitian(A)
- return full(Hermitian(A)^p)
+ return (Hermitian(A)^p)
  end
 
- # Otherwise, use Schur decomposition
  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 = full(A) ^ floor(p)
+ retmat = A ^ floor(p)
  #Real part
- retmat = retmat * full(powm(UpperTriangular(A), real(p - floor(p))))
- d = diag(A)
+ 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))
- 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
+ #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
- end
- if np_real_eigs
- warn("Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+ retmat = Q * R * Q'
  end
 
- if isreal(A) && !np_real_eigs
+ # if A has nonpositive real eigenvalues, retmat is a nonprincipal matrix power.
+ if isreal(retmat)
  return real(retmat)
  else
  return retmat
  end
 end
+^(A::AbstractMatrix, p::Number) = expm(p*logm(A))
 
 # Matrix exponential
 

base/linalg/linalg.jl

@@ -113,7 +113,6 @@ export
  ordschur,
  peakflops,
  pinv,
- powm,
  qr,
  qrfact!,
  qrfact,
@@ -263,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")
@@ -275,6 +273,8 @@ include("givens.jl")
 include("special.jl")
 include("bitarray.jl")
 include("ldlt.jl")
+include("schur.jl")
+
 
 include("arpack.jl")
 include("arnoldi.jl")

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

base/linalg/symmetric.jl

@@ -412,8 +412,8 @@ function ^{T<:Real}(A::Symmetric{T}, p::Integer)
  end
 end
 function ^{T<:Real}(A::Symmetric{T}, p::Real)
- F = eigfact(full(A))
- if isposdef(F)
+ 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'
@@ -435,7 +435,7 @@ end
 function ^{T}(A::Hermitian{T}, p::Real)
  n = checksquare(A)
  F = eigfact(A)
- if isposdef(F)
+ if all(λ -> λ ≥ 0, F.values)
  retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors'
  if T <: Real
  return Hermitian(retmat)
@@ -451,7 +451,7 @@ function ^{T}(A::Hermitian{T}, p::Real)
  end
 end
 
-function expm{T<:Real}(A::Symmetric{T})
+function expm(A::Symmetric)
  F = eigfact(A)
  return Symmetric((F.vectors * Diagonal(exp.(F.values))) * F.vectors')
 end

base/linalg/triangular.jl

@@ -1668,22 +1668,25 @@ Ac_ldiv_B(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw
 # 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{T}, p::Real) where T<:BlasFloat
+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, m, s = invsquaring(A0, theta)
  A = I - A
 
  # Compute accurate diagonal of I - T
- sqrt_diag!(A0,A,s)
+ sqrt_diag!(A0, A, s)
  for i = 1:n
- A[i,i] = -A[i,i]
+ A[i, i] = -A[i, i]
  end
  # Compute the Padé approximant
  c = 0.5 * (p - m) / (2 * m - 1)
@@ -1694,39 +1697,39 @@ function powm(A0::UpperTriangular{T}, p::Real) where T<:BlasFloat
  j4 = 4 * j
  c = (-p - j) / (j4 + 2)
  for i = 1:n
- @inbounds S[i,i] = S[i,i] + 1
+ @inbounds S[i, i] = S[i, i] + 1
  end
- copy!(Stmp,S)
- scale!(S,A,c)
- A_ldiv_B!(Stmp,S.data)
+ 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
+ @inbounds S[i, i] = S[i, i] + 1
  end
- copy!(Stmp,S)
- scale!(S,A,c)
- A_ldiv_B!(Stmp,S.data)
+ 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
+ S[i, i] = S[i, i] + 1
  end
- copy!(Stmp,S)
- scale!(S,A,-p)
- A_ldiv_B!(Stmp,S.data)
+ 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
+ @inbounds S[i, i] = S[i, i] + 1
  end
 
- blockpower!(A0,S,p/(2^s))
+ 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)))
+ 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
-^(A::LowerTriangular, p::Integer) = ^(A.', p::Integer).'
 powm(A::LowerTriangular, p::Real) = powm(A.', p::Real).'
 
 # Complex matrix logarithm for the upper triangular factor, see:
@@ -1949,9 +1952,11 @@ function sqrt_diag!(A0::UpperTriangular, A::UpperTriangular, s)
  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)
@@ -1964,13 +1969,11 @@ function invsquaring(A0::UpperTriangular, theta)
  dm1 = similar(d, n)
  s = 0
  for i = 1:n
- dm1[i] = d[i] - 1.
+ dm1[i] = d[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
@@ -1986,7 +1989,7 @@ function invsquaring(A0::UpperTriangular, theta)
  alpha2 = max(d2, d3)
  foundm = false
  if alpha2 <= theta[2]
- m = alpha2<=theta[1]?1:2
+ m = alpha2 <= theta[1] ? 1 : 2
  foundm = true
  end
 

test/linalg/dense.jl

@@ -512,37 +512,60 @@ end
  @test_throws ArgumentError diag(zeros(0,1),2)
 end
 
-@testset "New matrix ^ implementation" for elty in (Float64, Complex{Float64})
- A11 = convert(Matrix{elty}, [1 2 3; 4 7 1; 2 1 4])
-
- OLD_STDERR = STDERR
- rd,wr = redirect_stderr()
-
- @test A11^(1/2) ≈ sqrtm(A11)
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- @test A11^(-1/2) ≈ inv(sqrtm(A11))
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- @test A11^(3/4) ≈ sqrtm(A11) * sqrtm(sqrtm(A11))
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- @test A11^(-3/4) ≈ inv(A11) * sqrtm(sqrtm(A11))
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- @test A11^(17/8) ≈ A11^2 * sqrtm(sqrtm(sqrtm(A11)))
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- @test A11^(-17/8) ≈ inv(A11^2 * sqrtm(sqrtm(sqrtm(A11))))
- s = readline(rd)
- @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
-
- redirect_stderr(OLD_STDERR)
+@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