Fixed the algorithm for powers of a matrix. #21184

iamnapo · 2017-03-27T19:15:52Z

I did a few small bug-fixes and also changed the ^(A::AbstractMatrix, p::Integer) algorithm a bit to make use of powm(A0::UpperTriangular{T}, p::Real).

tkelman · 2017-03-27T19:38:17Z

base/linalg/linalg.jl

@@ -113,6 +113,7 @@ export
    ordschur,
    peakflops,
    pinv,
+    powm,


this is an internal helper function, it shouldn't be exported

iamnapo · 2017-03-27T19:52:44Z

fixed it.

iamnapo · 2017-03-28T00:41:20Z

I changed some and added a few more tests based on:
Deadman & Higham [March 2014] , Higham & Lin [September 2013]

kshyatt · 2017-03-28T19:49:52Z

Is this good to go?

andreasnoack · 2017-03-28T19:51:33Z

I don't know yet. I haven't looked carefully. What is the easiest way to resume? Should I open a PR from my branch or should this PR be against master?

tkelman · 2017-03-28T19:59:17Z

let's change this ~~branch~~ PR to be against master, so we can get a more complete view of the diff including the few commits from Andreas' branch

tkelman · 2017-03-28T20:00:14Z

and close/reopen to trigger CI, apparently changing the target branch doesn't do that

tkelman · 2017-03-28T20:04:23Z

base/linalg/dense.jl

-function ^(A::AbstractMatrix, p::Number)
+# Matrix power
+^(A::AbstractMatrix, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A),-p) : Base.power_by_squaring(A,p)
+function ^{T}(A::AbstractMatrix{T}, p::Real)


what's the right approach for complex p ? is there a different method for that?

Not sure. Maybe we can ask @higham what the right approach is.

Only approach I found is by definition: expm(p * logm(A)). But @stevengj pointed out correctly that it's very computationally expensive, so I didn't add it.

I think that's all there is for complex p. If you do an initial Schur decomp of A then expm and logm are applied to triangular matrices, which reduces the cost.

It would be reasonable to have ^(A::AbstractMatrix, p::Number) = expm(p * logm(A)) as a fallback. We can optimize it later if needed.

iamnapo · 2017-03-29T15:55:05Z

base/intfuncs.jl

@@ -193,7 +193,6 @@ end

 ^{T<:Integer}(x::T, p::T) = power_by_squaring(x,p)
 ^(x::Number, p::Integer)  = power_by_squaring(x,p)
-^(x, p::Integer)          = power_by_squaring(x,p)


not sure if this should be here.

This fallback method shouldn't be removed by this PR.

stevengj · 2017-03-29T15:55:39Z

base/linalg/dense.jl

    if isinteger(p)
-        return A^Integer(real(p))
+        return p < 0 ? Base.power_by_squaring(inv(A),-Integer(p)) : Base.power_by_squaring(A,Integer(p))


The promotion here is still wrong if p is not an integer, as explained in #21143. A^2.0 should produce a floating-point matrix even if A is integer.

The return type should be something like Matrix{Base.promote_op(^, eltype(A), typeof(p))}.

Also, note that power_by_squaring is pretty inefficient for matrices because it generates a lot of temporary matrices. We should really have a specialized _power_by_squaring!(Aᵖ::Matrix, A::Matrix, p::Integer) function that minimizes the number of temporary arrays and uses a pre-allocated output array. But that can be a separate PR, since it is just a performance optimization.

If A is an integer matrix, A^2 should be an integer matrix, but A^2.0 should be a floating-point matrix.

We definitely should not use diagonalization for integer p, because it will be way slower even than what we are doing now, despite the extra temporary matrices (which mostly waste memory, not time, for large A). e.g.

julia> A = randn(1000,1000); A = A'*A; # random SPD matrix julia> @time A^4; 0.428707 seconds (114.28 k allocations: 20.806 MiB, 20.50% gc time) julia> @time A^4.5; 3.669913 seconds (2.70 M allocations: 172.916 MiB, 2.36% gc time)

stevengj · 2017-03-29T17:10:41Z

base/linalg/dense.jl

@@ -310,15 +310,23 @@ 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
+function ^{T}(A::AbstractMatrix{T}, p::Integer)
+    # non-inversible matrix to negative power is undefined.
+    if (p < 0) && (det(A) == 0)


This check does not seem like a good idea. First, inv(A) will already throw an error if A is singular, and just as quickly as det will. Second, computing det(A) means that you compute the LU factorization of A twice. Third, you could get false positives because of underflow (e.g. det([1e-300 0; 0 1e-300]) == 0 even though the matrix is not singular).

(Basically, you should forget everything you learned in high school about the determinant as a computational tool.)

stevengj · 2017-03-29T17:16:13Z

base/linalg/dense.jl

-    if isinteger(p)
-        return p < 0 ? Base.power_by_squaring(inv(A),-Integer(p)) : Base.power_by_squaring(A,Integer(p))
+    if isinteger(p) && !(typeof(p) <: Int64)
+        return float.(A^Integer(p))


This is not right. First, Matrix{Integer}^Integer should be an integer matrix (even if p is not Int64). Second, for Matrix{Int}^Float64 we need to convert to floating point and then exponentiate in order to avoid overflow problems. Third, calling float.(X) will make an extra copy even if the matrix is already floating point.

We want something like:

convert(Matrix{promote_op(^, eltype(A), typeof(p)}, A)^Integer(p)

although that isn't quite right either: if A is not a Matrix, but is some other AbstractMatrix type, then we probably want the same type. Maybe

T = Base.promote_op(^, eltype(A), typeof(p)) return (T == eltype(A) ? A : copy!(similar(A, T), A))^Integer(p)

iamnapo · 2017-03-29T20:04:40Z

One thing I'm not sure is what should happen is A is singular and p is negative real

stevengj · 2017-03-29T20:34:30Z

base/linalg/diagonal.jl

@@ -285,7 +285,6 @@ logm(D::Diagonal) = Diagonal(log.(D.diag))
 logm(D::Diagonal{<:AbstractMatrix}) = Diagonal(logm.(D.diag))
 sqrtm(D::Diagonal) = Diagonal(sqrt.(D.diag))
 sqrtm(D::Diagonal{<:AbstractMatrix}) = Diagonal(sqrtm.(D.diag))
-^(D::Diagonal, p::Real) = Diagonal((D.diag).^p)


Why was this method removed?

stevengj · 2017-03-29T20:39:25Z

base/linalg/dense.jl

+        end
+    end
+    if np_real_eigs
+        warn("Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")


Maybe pass once=true here so that people don't see repeated warnings?

you mean in something like (A^0.5) * (A^0.5)? If so, I'm not sure how to come up with a right key for the warning.

The default key is the warning message, which seems fine here.

stevengj · 2017-03-29T20:41:00Z

base/linalg/dense.jl

+        warn("Matrix with nonpositive real eigenvalues, a nonprincipal matrix power will be returned.")
+    end
+
+    if isreal(A) && ~np_real_eigs


!np_real_eigs

stevengj · 2017-03-29T20:45:25Z

base/linalg/triangular.jl

+    sqrt_diag!(A0, A, s)
+
+    # Compute the Gauss-Legendre quadrature formula
+    x, w = Base.QuadGK.gauss(Float64, m)


This won't work now that QuadGK was removed from Base, no?

Since it looks like m ≤ length(theta), we can just precompute these and store them in a const array.

iamnapo · 2017-03-30T18:55:45Z

@stevengj
If I put once=true then the warning will only show once, even if later I have a different matrix.
I removed ^(D::Diagonal, p::Real) because it created ambigiuties with ^(x, p::Integer), for some reason.
Regarding logm: #21179

jebej · 2017-04-04T20:23:08Z

test/linalg/dense.jl

@@ -512,6 +512,49 @@ end
    @test_throws ArgumentError diag(zeros(0,1),2)
 end

+@testset "New matrix ^ implementation" for elty in (Float64, Complex{Float64})
+    A11 = convert(Matrix{elty}, rand(10, 10))


You might want to just have rand(elty,10,10) here so that you actually make complex matrices.

tkelman · 2017-04-05T11:05:15Z

base/linalg/dense.jl

+    # 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


should this be the eps of the element type? not sure how likely it would be to matter

There's also typemin() for example, but I thought default eps() is generally a safe condition to produce a warning. You think I should change it nevertheless?

This is not a good way to do the test because it depends on overall scaling of the matrix (if you multiply the whole matrix by 1e±100 the result will change).

Also, shouldn't you be checking abs(imag), not imag? And shouldn't the warning be for real < 0, not real ≤ 0?

I'm not sure what a robust test would be here. Maybe abs(imag(d[i])) ≤ ɛ && real(d[i]) < 0, where ɛ = eps(maxabs(d))?

tkelman · 2017-04-05T11:15:53Z

It evidently isn't tested, but tightening ^(A::AbstractMatrix, p::Number) to ^(A::AbstractMatrix, p::Real) and not having a method for other subtypes of Number could be breaking

tkelman · 2017-04-05T11:19:42Z

base/linalg/symmetric.jl

+        retmat = Base.power_by_squaring(A, p)
+    end
+    for i = 1:n
+        retmat[i,i] = real(retmat[i,i])


I know there are separate issues where this has been discussed, but just checking - any imaginary components here should only ever be round-off, right? Does blas hemm touch the diagonal imaginary components or wouldn't it leave them alone? Is the issue round-off for julia-generic or non-blas element types?

Yes, the imaginary part is guaranteed to be zero, so any nonzero imaginary part is purely due to roundoff errors.

The errors arise for any values, including blas floats, because multiplying by the matrix of eigenvectors uses the general matmul routine. I don't think there is a specialized BLAS routine for matrix * real diagonal * matrix', which is the operation here.

Sorry, for Hermitian^p the operation is hermitian * hermitian, but again there is no specialized BLAS method for this because in general the result would not be Hermitian. (The result is Hermitian here because the matrices being multiplied are A^p and A^q, i.e. different powers of the same Hermitian matrix.) So the roundoff errors still arise for any values, including BLAS floats.

For example, note the nonzero imaginary parts on the diagonal of:

julia> A = rand(Complex128,10,10); A = Hermitian(A + A') 10×10 Hermitian{Complex{Float64},Array{Complex{Float64},2}}: 0.0861346-0.0im 0.837547-0.743436im … 1.09493-0.163298im 0.837547+0.743436im 1.96064-0.0im 0.491744-0.792954im 0.869364-0.131544im 1.06965-0.252339im 0.370136+0.0155668im 0.83727+0.303993im 1.18724+0.108695im 1.38671+0.191601im 1.50475+0.357721im 1.00663-0.421604im 1.02963-0.13435im 1.1999+0.0748608im 0.370293+0.503285im … 1.13011-0.0552782im 1.56458+0.546159im 0.512528-0.0406777im 0.721986+0.445042im 0.394002+0.645965im 0.935305-0.473984im 0.971132+0.493203im 0.741598-0.175156im 1.91631+0.111909im 0.868631+0.834707im 1.09493+0.163298im 0.491744+0.792954im 0.00985939-0.0im julia> imag(diag(A^2 * A^3)) 10-element Array{Float64,1}: 0.0 5.68434e-14 0.0 0.0 0.0 4.54747e-13 0.0 5.68434e-14 0.0 0.0

tkelman · 2017-04-05T11:23:34Z

base/linalg/symmetric.jl

+    end
+end
+
+function expm{T<:Real}(A::Symmetric{T})


would this, logm, and sqrtm for Symmetric of non-Real element types become no-method errors here?

No, these are covered by Hermitian methods

what about complex symmetric, non hermitian?

oooh I'm sorry, for some reason I thought symmetric matrices are always hermitian

it's a bit of a niche type, a fair number of routines even in blas or lapack aren't always implemented for c/z symmetric

stevengj · 2017-04-21T01:56:48Z

LGTM. The Travis failure seems unrelated.

tkelman · 2017-04-21T02:08:48Z

base/linalg/symmetric.jl

+    end
+end
+function ^{T<:Real}(A::Symmetric{T}, p::Real)
+    F = eigfact(full(A))


why is it we need the full call here?

(In general, a lot of the full calls seem odd to me; we seem better off adding additional methods elsewhere, but I thought that could be left for another PR. This particular one is especially egregious though because eigfact should be especially efficient for real-symmetric matrices. What happens if you just delete all of the full calls?)

jebej · 2017-04-21T13:05:08Z

Thanks for pushing with this. It would be really great to have it for 0.6.0

jebej · 2017-04-21T13:08:04Z

base/linalg/dense.jl

-
-function ^(A::AbstractMatrix, p::Number)
+# Matrix power
+^{T}(A::AbstractMatrix{T}, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A),-Integer(p)) : Base.power_by_squaring(A,Integer(p))


Is Integer(p) needed here?

tkelman · 2017-04-21T14:30:16Z

The tightening of the allowed power type here from Number to Real would currently break [1 2; 3 4]^im

iamnapo · 2017-04-21T15:07:51Z

@tkelman I'm not sure where are you referring to

tkelman · 2017-04-21T15:09:47Z

it's collapsed right now, but the line at #21184 (comment)

iamnapo · 2017-04-21T15:16:29Z

Currently [1 2; 3 4]^im produces a MethodError. What do you think should happen?

(I plan on working at function ^{T}(A::AbstractMatrix{T}, p::Complex) at another PR)

tkelman · 2017-04-21T15:19:40Z

On master it gives a consistent answer with http://www.wolframalpha.com/input/?i=N%5B%7B%7B1,+2%7D,+%7B3,+4%7D%7D%5EI%5D. I don't know of a use case for this off the top of my head, but that doesn't seem like a good reason to change it from consistent with Wolfram to an error.

stevengj · 2017-04-21T16:34:04Z

@iamnapo, we should just add the ^(A::AbstractMatrix, p::Number) = expm(p*logm(A)) fallback routine and a test.

tkelman · 2017-04-21T23:46:46Z

base/linalg/dense.jl

+            retmat = retmat * powm!(UpperTriangular(float.(A)), real(p - floor(p)))
+        end
+    else
+        S,Q,d = schur(complex(full(A)))


another full that may have been missed, or was this one needed?

This one is needed

good to know - what was happening without it?

if A is symmetric with complex values, but not hermitian

@iamnapo, do you mean that there was a methoderror for schur(::Symmetric{<:Complex}}? It would be better to solve that by adding a method to shur or shurfact that calls full.

Looks like we just need:

schur(A::Symmetric{<:Complex}) = schur(full(A))

Yes. I'll add it.

tkelman · 2017-04-21T23:57:27Z

base/linalg/symmetric.jl

    @eval begin
-
-        function ($funm)(A::Symmetric)
+        function ($funm){T<:Real}(A::Symmetric{T})


is the problem for Symmetric{<:Complex} that isposdef doesn't work?

No, the issue is that Symmetric{<:Complex} may have ill-conditioned eigenvectors F.vectors, so it is not reliable to use diagonalization to compute matrix functions. Whereas for real-symmetric matrices the eigenvectors are unitary, which we are exploiting here: F.vectors' == inv(F.vectors). The complex analogue is a Hermitian matrix, not complex-symmetric.

Cool. Useful to write that down somewhere, I wouldn't have remembered that detail from a linear algebra textbook.

stevengj · 2017-04-22T14:29:07Z

base/linalg/schur.jl

@@ -107,6 +107,7 @@ function schur(A::StridedMatrix)
    SchurF = schurfact(A)
    SchurF[:T], SchurF[:Z], SchurF[:values]
 end
+schur(A::AbstractMatrix) = schur(full(A))


The problem with this definition is that it will produce an infinite loop for any AbstractMatrix type where full(A) === A and there is no more-specific schur method.

I would just add schur(A::Triangular{<:Complex}) = schur(full(A)) (and anything else as needed).

Isn't full(A) always of type StridedMatrix?

stevengj · 2017-04-22T22:06:27Z

base/linalg/schur.jl

@@ -111,8 +111,7 @@ 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))


These methods can obviously be made vastly more efficient. For example, an upper-triangular matrix is already in Schur form. But that can be left for another PR.

Of course. I just added these definitions to make it work for now. I plan to work on schur() after this PR is closed.

tkelman · 2017-04-22T22:40:24Z

Do the original 2 commits by a different author pass on their own? Any failing commits should be squashed here, but git is bad at attribution if you squash together commits by different authors. Otherwise this lgtm now, thanks for the epic persistence @iamnapo !

iamnapo · 2017-04-22T23:18:04Z

I don't think there will be a problem, but in the case there is, you can use command line to --amend authorship, if that's what you mean.

fixed Matrix^Real_Number functionality reverted logm changes remove ambiguous Diagonal^Real method

…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

and be uniform about space after # in comments

Sacha0 · 2017-04-23T16:54:22Z

thanks for the epic persistence @iamnapo !

Epic persistence indeed! Thanks @iamnapo! :)

iamnapo · 2017-04-24T16:46:48Z

is this good to go?

tkelman reviewed Mar 27, 2017

View reviewed changes

ararslan added domain:linear algebra Linear algebra needs tests Unit tests are required for this change labels Mar 27, 2017

ararslan removed the needs tests Unit tests are required for this change label Mar 28, 2017

iamnapo changed the base branch from anj/powm to master March 28, 2017 19:59

tkelman closed this Mar 28, 2017

tkelman reopened this Mar 28, 2017

tkelman reviewed Mar 28, 2017

View reviewed changes

iamnapo commented Mar 29, 2017

View reviewed changes

stevengj reviewed Mar 29, 2017

View reviewed changes

tkelman mentioned this pull request Apr 4, 2017

No method for isless(::Complex128,::Int64) #21273

Closed

jebej reviewed Apr 4, 2017

View reviewed changes

tkelman reviewed Apr 5, 2017

View reviewed changes

stevengj approved these changes Apr 21, 2017

View reviewed changes

tkelman reviewed Apr 21, 2017

View reviewed changes

jebej reviewed Apr 21, 2017

View reviewed changes

tkelman reviewed Apr 21, 2017

View reviewed changes

iamnapo changed the title ~~Fixed the algorithm for real powers of a matrix.~~ Fixed the algorithm for powers of a matrix. Apr 22, 2017

stevengj reviewed Apr 22, 2017

View reviewed changes

mfasi and others added 3 commits April 23, 2017 08:02

Make the algorithm for real powers of a matrix robust

0935284

fixed Matrix^Real_Number functionality reverted logm changes remove ambiguous Diagonal^Real method

Remove some blank lines at start or end of blocks

2fbeba3

and be uniform about space after # in comments

tkelman merged commit f56147d into JuliaLang:master Apr 24, 2017

iamnapo deleted the anj/powm branch April 24, 2017 18:43

jiahao mentioned this pull request May 1, 2017

Implement more matrix functions #5840

Closed

17 tasks

This was referenced May 1, 2017

Error in matrix^power if eigenvalues are complex #16930

Closed

WIP: Improve accuracy of ^ for ill-conditioned matrices #12584

Closed

fredrikekre mentioned this pull request Sep 30, 2017

change order of exponent for neg. matrix powers #20071

Closed

oxinabox mentioned this pull request Apr 30, 2020

Matrix Powers are type-unstable: sometime Real sometimes Complex #35667

Open

Fixed the algorithm for powers of a matrix. #21184

Fixed the algorithm for powers of a matrix. #21184

Conversation

iamnapo commented Mar 27, 2017 • edited Loading

Choose a reason for hiding this comment

iamnapo commented Mar 27, 2017

iamnapo commented Mar 28, 2017

kshyatt commented Mar 28, 2017

andreasnoack commented Mar 28, 2017

tkelman commented Mar 28, 2017 • edited Loading

tkelman commented Mar 28, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iamnapo Mar 28, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Mar 29, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Mar 29, 2017 • edited Loading

Choose a reason for hiding this comment

stevengj Mar 29, 2017 • edited Loading

Choose a reason for hiding this comment

stevengj Mar 29, 2017 • edited Loading

Choose a reason for hiding this comment

iamnapo commented Mar 29, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iamnapo commented Mar 30, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Apr 5, 2017 • edited Loading

Choose a reason for hiding this comment

tkelman commented Apr 5, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Apr 13, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj commented Apr 21, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jebej commented Apr 21, 2017

Choose a reason for hiding this comment

tkelman commented Apr 21, 2017

iamnapo commented Apr 21, 2017

tkelman commented Apr 21, 2017

iamnapo commented Apr 21, 2017

tkelman commented Apr 21, 2017

stevengj commented Apr 21, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iamnapo Apr 22, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Apr 22, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

iamnapo Apr 22, 2017 • edited Loading

Choose a reason for hiding this comment

iamnapo commented Mar 27, 2017 •

edited

Loading

tkelman commented Mar 28, 2017 •

edited

Loading

iamnapo Mar 28, 2017 •

edited

Loading

stevengj Mar 29, 2017 •

edited

Loading

stevengj Mar 29, 2017 •

edited

Loading

stevengj Mar 29, 2017 •

edited

Loading

stevengj Mar 29, 2017 •

edited

Loading

stevengj Apr 5, 2017 •

edited

Loading

stevengj Apr 13, 2017 •

edited

Loading

iamnapo Apr 22, 2017 •

edited

Loading

stevengj Apr 22, 2017 •

edited

Loading

iamnapo Apr 22, 2017 •

edited

Loading