Skip to content

RFC : isposdef, bkfact, eigvals! for Symmetric Matrix #15458

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.

Sign up for GitHub

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

Closed
wants to merge 2 commits into from
Closed

RFC : isposdef, bkfact, eigvals! for Symmetric Matrix #15458

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d590582
RFC:isposdef,bkfact,eigvals! for Symmetric Matrix
amitjamadagni Mar 11, 2016
487d63b
Revamped isposdef
amitjamadagni Mar 12, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions base/linalg/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,17 @@ issymmetric{T<:Complex,S}(A::Hermitian{T,S}) = all(imag(A.data) .== 0)
issymmetric(A::Symmetric) = true
transpose(A::Symmetric) = A
ctranspose{T<:Real}(A::Symmetric{T}) = A

function isposdef{T<:Real,S}(A::HermOrSym{T,S})
try
f = cholfact(A)
catch e
isa(e, LinAlg.PosDefException) || rethrow(e)
return false
end
true
end

function ctranspose(A::Symmetric)
AC = ctranspose(A.data)
return Symmetric(AC, ifelse(A.uplo == 'U', :L, :U))
Expand Down Expand Up @@ -138,7 +149,7 @@ A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{
*(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
*(A::StridedMatrix, B::HermOrSym) = A*full(B)

bkfact(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issymmetric(A))
bkfact(A::HermOrSym) = bkfact(full(A.data), symbol(A.uplo), issymmetric(A))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't want to call full implicitly on sparse matrices since it might not fit in memory. It will be a big task to write a sparse Bunch-Kaufman factorization so we'll probably not be able to support this method for a while. Meanwhile, we could consider an error message suggesting ldltfact with a warning that pivoting is only done as part of the symbolic factorization to reduce fill in.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't want to call full implicitly on sparse matrices since it might not fit in memory. It will be a big task to write a sparse Bunch-Kaufman factorization so we'll probably not be able to support this method for a while.

Is there some reference which I could look into to start working on this ?!

Meanwhile, we could consider an error message suggesting ldltfact with a warning that
pivoting is only done as part of the symbolic factorization to reduce fill in.

I understand the synopsis that it should be an error message but I guess my knowledge is pretty limited and it would be really helpful if you point me to some sources on the understanding of this ! Thanks :-)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't point to any specific literature but searching for "sparse", "indefinite", and "pivoting" would probably give you something. @tkelman is more into this so he might be able to help.

factorize(A::HermOrSym) = bkfact(A)

# Is just RealHermSymComplexHerm, but type alias seems to be broken
Expand All @@ -163,7 +174,8 @@ eigfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, vl::Real,
# Because of #6721 it is necessary to specify the parameters explicitly here.
eigfact{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}, vl::Real, vh::Real) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigfact!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh))

eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
# eigvals!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}) = LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)[1]
eigvals!{T<:BlasReal,S<:AbstractMatrix}(A::RealHermSymComplexHerm{T,S}) = LAPACK.syevr!('N', 'A', A.uplo, full(A.data), 0.0, 0.0, 0, 0, -1.0)[1]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as for bkfact. Maybe we could use eigs to return some of the eigenvalues, but so far we have followed MATLAB and kept the iterative methods in separate functions, i.e. svds and eigs.

# Because of #6721 it is necessary to specify the parameters explicitly here.
eigvals{T1<:Real,T2}(A::RealHermSymComplexHerm{T1,T2}) = (T = eltype(A); S = promote_type(Float32, typeof(zero(T)/norm(one(T)))); eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A)))

Expand Down