-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fixed the algorithm for powers of a matrix. #21184
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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}) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. is the problem for There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. No, the issue is that There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Cool. Useful to write that down somewhere, I wouldn't have remembered that detail from a linear algebra textbook. |
||
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 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For example, note the nonzero imaginary parts on the diagonal of:
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
#17367