Skip to content
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

Check that matrix is symmetric/Hermitian in dense Cholesky #16799

Merged
merged 2 commits into from
Jun 11, 2016
Merged

Check that matrix is symmetric/Hermitian in dense Cholesky #16799

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
597d945
Check that matrix is symmetric/Hermitian in dense Cholesky
andreasnoack Jun 3, 2016
9900fb6
Factor copy-pasted non-hermitian error into a helper function
tkelman Jun 10, 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
3 changes: 1 addition & 2 deletions base/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,11 +234,10 @@ include("linalg/lq.jl")
include("linalg/eigen.jl")
include("linalg/svd.jl")
include("linalg/schur.jl")
include("linalg/symmetric.jl")
include("linalg/cholesky.jl")
include("linalg/lu.jl")

include("linalg/bunchkaufman.jl")
include("linalg/symmetric.jl")
include("linalg/diagonal.jl")
include("linalg/bidiag.jl")
include("linalg/uniformscaling.jl")
Expand Down
226 changes: 174 additions & 52 deletions base/linalg/cholesky.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,30 @@
##########################
# Cholesky Factorization #
##########################

# The dispatch structure in the chol!, chol, cholfact, and cholfact! methods is a bit
# complicated and some explanation is therefore provided in the following
#
# In the methods below, LAPACK is called when possible, i.e. StridedMatrices with Float32,
# Float64, Complex{Float32}, and Complex{Float64} element types. For other element or
# matrix types, the unblocked Julia implementation in _chol! is used. For cholfact
# and cholfact! pivoting is supported through a Val{Bool} argument. A type argument is
# necessary for type stability since the output of cholfact and cholfact! is either
# Cholesky or PivotedCholesky. The latter is only
# supported for the four LAPACK element types. For other types, e.g. BigFloats Val{true} will
# give an error. It is required that the input is Hermitian (including real symmetric) either
# through the Hermitian and Symmetric views or exact symmetric or Hermitian elements which
# is checked for and an error is thrown if the check fails. The dispatch
# is further complicated by a limitation in the formulation of Unions. The relevant union
# would be Union{Symmetric{T<:Real,S}, Hermitian} but, right now, it doesn't work in Julia
# so we'll have to define methods for the two elements of the union separately.

# FixMe? The dispatch below seems overly complicated. One simplification could be to
# merge the two Cholesky types into one. It would remove the need for Val completely but
# the cost would be extra unnecessary/unused fields for the unpivoted Cholesky and runtime
Copy link
Contributor

Choose a reason for hiding this comment

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

So this renames the two-argument chol! method with the upper/lower type second argument, making it no longer available under the chol! name at all, right? I did a quick grep of packages and luckily it looks like no one was calling that in anything that's registered.

# checks of those fields before calls to LAPACK to check which version of the Cholesky
# factorization the type represents.

immutable Cholesky{T,S<:AbstractMatrix} <: Factorization{T}
factors::S
uplo::Char
Expand All @@ -23,24 +47,27 @@ function CholeskyPivoted{T}(A::AbstractMatrix{T}, uplo::Char, piv::Vector{BlasIn
CholeskyPivoted{T,typeof(A)}(A, uplo, piv, rank, tol, info)
end

function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{UpperTriangular})

# _chol!. Internal methods for calling unpivoted Cholesky
## BLAS/LAPACK element types
function _chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{UpperTriangular})
C, info = LAPACK.potrf!('U', A)
return @assertposdef UpperTriangular(C) info
end
function chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{LowerTriangular})
function _chol!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{LowerTriangular})
C, info = LAPACK.potrf!('L', A)
return @assertposdef LowerTriangular(C) info
end
chol!(A::StridedMatrix) = chol!(A, UpperTriangular)

function chol!{T}(A::AbstractMatrix{T}, ::Type{UpperTriangular})
## Non BLAS/LAPACK element types (generic)
function _chol!(A::AbstractMatrix, ::Type{UpperTriangular})
n = checksquare(A)
@inbounds begin
for k = 1:n
for i = 1:k - 1
A[k,k] -= A[i,k]'A[i,k]
end
Akk = chol!(A[k,k], UpperTriangular)
Akk = _chol!(A[k,k], UpperTriangular)
A[k,k] = Akk
AkkInv = inv(Akk')
for j = k + 1:n
Expand All @@ -53,14 +80,14 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{UpperTriangular})
end
return UpperTriangular(A)
end
function chol!{T}(A::AbstractMatrix{T}, ::Type{LowerTriangular})
function _chol!(A::AbstractMatrix, ::Type{LowerTriangular})
n = checksquare(A)
@inbounds begin
for k = 1:n
for i = 1:k - 1
Copy link
Contributor

Choose a reason for hiding this comment

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

if I'm reading this right, is there no longer a definition for AbstractMatrix at all?

Copy link
Member Author

Choose a reason for hiding this comment

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

I'd restricted it to StridedMatrix to avoid overlap with Hermitian but it's not necessary in the version of this I ended up with so I've changed it back to AbstractMatrix.

A[k,k] -= A[k,i]*A[k,i]'
end
Akk = chol!(A[k,k], LowerTriangular)
Akk = _chol!(A[k,k], LowerTriangular)
A[k,k] = Akk
AkkInv = inv(Akk)
for j = 1:k
Expand All @@ -78,83 +105,152 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{LowerTriangular})
return LowerTriangular(A)
end

"""
chol(A::AbstractMatrix) -> U

Compute the Cholesky factorization of a positive definite matrix `A`
and return the UpperTriangular matrix `U` such that `A = U'U`.
"""
function chol{T}(A::AbstractMatrix{T})
S = promote_type(typeof(chol(one(T))), Float32)
AA = similar(A, S, size(A))
copy!(AA, A)
chol!(AA)
end
function chol!(x::Number, uplo)
## Numbers
function _chol!(x::Number, uplo)
rx = real(x)
if rx != abs(x)
throw(ArgumentError("x must be positive semidefinite"))
end
rxr = sqrt(rx)
convert(promote_type(typeof(x), typeof(rxr)), rxr)
end

non_hermitian_error(f) = throw(ArgumentError("matrix is not symmetric/" *
"Hermitian. This error can be avoided by calling $f(Hermitian(A)) " *
"which will ignore either the upper or lower triangle of the matrix."))

# chol!. Destructive methods for computing Cholesky factor of real symmetric or Hermitian
# matrix
chol!(A::Hermitian) = _chol!(A.data, UpperTriangular)
chol!{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}) = _chol!(A.data, UpperTriangular)
function chol!(A::StridedMatrix)
ishermitian(A) || non_hermitian_error("chol!")
return _chol!(A, UpperTriangular)
end



# chol. Non-destructive methods for computing Cholesky factor of a real symmetric or
# Hermitian matrix. Promotes elements to a type that is stable under square roots.
function chol(A::Hermitian)
T = promote_type(typeof(chol(one(eltype(A)))), Float32)
AA = similar(A, T, size(A))
copy!(AA, A.data)
chol!(Hermitian(AA))
end
function chol{T<:Real,S<:AbstractMatrix}(A::Symmetric{T,S})
TT = promote_type(typeof(chol(one(T))), Float32)
AA = similar(A, TT, size(A))
copy!(AA, A.data)
chol!(Hermitian(AA))
end

## for StridedMatrices, check that matrix is symmetric/Hermitian
"""
chol(A) -> U

Compute the Cholesky factorization of a positive definite matrix `A`
and return the UpperTriangular matrix `U` such that `A = U'U`.
"""
function chol(A::AbstractMatrix)
ishermitian(A) || non_hermitian_error("chol")
return chol(Hermitian(A))
end

## Numbers
"""
Copy link
Contributor

Choose a reason for hiding this comment

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

Isn't the function signature here wrong?

Copy link
Member Author

Choose a reason for hiding this comment

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

Are you talking about the args... part?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yeah, does that need to be in the sig?

Copy link
Member Author

Choose a reason for hiding this comment

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

Often the documentation signatures cheat a bit for simplification. I think that having the args... part in the docstring would confuse more than help. Maybe you could write that "extra arguments are ignored". However, this is not part of the change in this PR so if we want that, it could be done in a separate PR.

Copy link
Contributor

Choose a reason for hiding this comment

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

Sure; I just want to make sure the docs still work! On the other hand I don't think the masses are really clamoring for chol(Number) that we need to worry about it :).

chol(x::Number) -> y

Compute the square root of a non-negative number `x`.
"""
chol(x::Number, args...) = chol!(x, nothing)
chol(x::Number, args...) = _chol!(x, nothing)


cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) =
_cholfact!(A, Val{false}, uplo)
cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
_cholfact!(A, Val{true}, uplo; tol = tol)
cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol = :U) =
_cholfact!(A, Val{false}, uplo)

function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U)
# cholfact!. Destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting
function cholfact!(A::Hermitian, uplo::Symbol, ::Type{Val{false}})
if uplo == :U
return Cholesky(chol!(A, UpperTriangular).data, uplo)
Cholesky(_chol!(A.data, UpperTriangular).data, 'U')
else
return Cholesky(chol!(A, LowerTriangular).data, uplo)
Cholesky(_chol!(A.data, LowerTriangular).data, 'L')
end
end
function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{true}}, uplo::Symbol=:U; tol=0.0)
uplochar = char_uplo(uplo)
A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
return CholeskyPivoted{T,typeof(A)}(A, uplochar, piv, rank, tol, info)
function cholfact!{T<:Real,S}(A::Symmetric{T,S}, uplo::Symbol, ::Type{Val{false}})
if uplo == :U
Cholesky(_chol!(A.data, UpperTriangular).data, 'U')
else
Cholesky(_chol!(A.data, LowerTriangular).data, 'L')
end
end

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact!(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
cholfact!(A, uplo::Symbol, Val{false}) -> Cholesky

The same as `cholfact`, but saves space by overwriting the input `A`, instead
of creating a copy. An `InexactError` exception is thrown if the factorisation
produces a number not representable by the element type of `A`, e.g. for
integer types.
"""
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}})
if uplo == :U
Cholesky(chol!(A, UpperTriangular).data, 'U')
else
Cholesky(chol!(A, LowerTriangular).data, 'L')
end
ishermitian(A) || non_hermitian_error("cholfact!")
return cholfact!(Hermitian(A), uplo, Val{false})
end

### Default to no pivoting (and storing of upper factor) when not explicit
cholfact!(A::Hermitian, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false})
cholfact!{T<:Real,S}(A::Symmetric{T,S}, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false})
#### for StridedMatrices, check that matrix is symmetric/Hermitian
function cholfact!(A::StridedMatrix, uplo::Symbol = :U)
ishermitian(A) || non_hermitian_error("cholfact!")
return cholfact!(Hermitian(A), uplo)
end


## With pivoting
### BLAS/LAPACK element types
function cholfact!{T<:BlasReal,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S},
uplo::Symbol, ::Type{Val{true}}; tol = 0.0)
uplochar = char_uplo(uplo)
AA, piv, rank, info = LAPACK.pstrf!(uplochar, A.data, tol)
return CholeskyPivoted{eltype(AA),typeof(AA)}(AA, uplochar, piv, rank, tol, info)
end

### Non BLAS/LAPACK element types (generic). Since generic fallback for pivoted Cholesky
### is not implemented yet we throw an error
cholfact!{T<:Real,S}(A::RealHermSymComplexHerm{T,S}, uplo::Symbol, ::Type{Val{true}};
tol = 0.0) =
throw(ArgumentError("generic pivoted Cholesky factorization is not implemented yet"))

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact!(A::StridedMatrix, uplo::Symbol, Val{true}) -> PivotedCholesky
cholfact!(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted

The same as `cholfact`, but saves space by overwriting the input `A`, instead
of creating a copy. An `InexactError` exception is thrown if the
factorisation produces a number not representable by the element type of `A`,
e.g. for integer types.
"""
cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
throw(ArgumentError("generic pivoted Cholesky fectorization is not implemented yet"))
cholfact!(A::StridedMatrix, uplo::Symbol = :U) = cholfact!(A, uplo, Val{false})
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0)
ishermitian(A) || non_hermitian_error("cholfact!")
return cholfact!(Hermitian(A), uplo, Val{true}; tol = tol)
end



# cholfact. Non-destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting
cholfact(A::Hermitian, uplo::Symbol, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), uplo, Val{false})
cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}, uplo::Symbol, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), uplo, Val{false})

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
cholfact(A, uplo::Symbol, Val{false}) -> Cholesky

Compute the Cholesky factorization of a dense symmetric positive definite matrix `A`
and return a `Cholesky` factorization.
Expand All @@ -164,11 +260,34 @@ The triangular Cholesky factor can be obtained from the factorization `F` with:
The following functions are available for `Cholesky` objects: `size`, `\\`, `inv`, `det`.
A `PosDefException` exception is thrown in case the matrix is not positive definite.
"""
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{false})
function cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}})
ishermitian(A) || non_hermitian_error("cholfact")
return cholfact(Hermitian(A), uplo, Val{false})
end

### Default to no pivoting (and storing of upper factor) when not explicit
cholfact(A::Hermitian, uplo::Symbol = :U) = cholfact(A, uplo, Val{false})
cholfact{T<:Real,S<:StridedMatrix}(A::Symmetric{T,S}, uplo::Symbol = :U) =
cholfact(A, uplo, Val{false})
#### for StridedMatrices, check that matrix is symmetric/Hermitian
function cholfact(A::StridedMatrix, uplo::Symbol = :U)
ishermitian(A) || non_hermitian_error("cholfact")
return cholfact(Hermitian(A), uplo)
end


## With pivoting
cholfact(A::Hermitian, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)),
uplo, Val{true}; tol = tol)
cholfact{T<:Real,S<:StridedMatrix}(A::RealHermSymComplexHerm{T,S}, uplo::Symbol,
::Type{Val{true}}; tol = 0.0) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)),
uplo, Val{true}; tol = tol)

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact(A::StridedMatrix, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted
cholfact(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted

Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix `A`
and return a `CholeskyPivoted` factorization.
Expand All @@ -179,16 +298,19 @@ The following functions are available for `PivotedCholesky` objects: `size`, `\\
The argument `tol` determines the tolerance for determining the rank.
For negative values, the tolerance is the machine precision.
"""
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{true}; tol = tol)

cholfact{T}(A::StridedMatrix{T}, uplo::Symbol = :U) = cholfact(A, uplo, Val{false})
function cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0)
ishermitian(A) || non_hermitian_error("cholfact")
return cholfact(Hermitian(A), uplo, Val{true}; tol = tol)
end

## Number
function cholfact(x::Number, uplo::Symbol=:U)
xf = fill(chol(x), 1, 1)
Cholesky(xf, uplo)
end



function convert{Tnew,Told,S}(::Type{Cholesky{Tnew}}, C::Cholesky{Told,S})
Cnew = convert(AbstractMatrix{Tnew}, C.factors)
Cholesky{Tnew, typeof(Cnew)}(Cnew, C.uplo)
Expand Down
10 changes: 5 additions & 5 deletions doc/stdlib/linalg.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f

``lufact!`` is the same as :func:`lufact`\ , but saves space by overwriting the input ``A``\ , instead of creating a copy. An ``InexactError`` exception is thrown if the factorisation produces a number not representable by the element type of ``A``\ , e.g. for integer types.

.. function:: chol(A::AbstractMatrix) -> U
.. function:: chol(A) -> U

.. Docstring generated from Julia source

Expand All @@ -169,13 +169,13 @@ Linear algebra functions in Julia are largely implemented by calling functions f

Compute the square root of a non-negative number ``x``\ .

.. function:: cholfact(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
.. function:: cholfact(A, uplo::Symbol, Val{false}) -> Cholesky

.. Docstring generated from Julia source

Compute the Cholesky factorization of a dense symmetric positive definite matrix ``A`` and return a ``Cholesky`` factorization. The ``uplo`` argument may be ``:L`` for using the lower part or ``:U`` for the upper part of ``A``\ . The default is to use ``:U``\ . The triangular Cholesky factor can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``\ . The following functions are available for ``Cholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ . A ``PosDefException`` exception is thrown in case the matrix is not positive definite.

.. function:: cholfact(A::StridedMatrix, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted
.. function:: cholfact(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted

.. Docstring generated from Julia source

Expand Down Expand Up @@ -205,13 +205,13 @@ Linear algebra functions in Julia are largely implemented by calling functions f

This method uses the CHOLMOD library from SuiteSparse, which only supports doubles or complex doubles. Input matrices not of those element types will be converted to ``SparseMatrixCSC{Float64}`` or ``SparseMatrixCSC{Complex128}`` as appropriate.

.. function:: cholfact!(A::StridedMatrix, uplo::Symbol, Val{false}) -> Cholesky
.. function:: cholfact!(A, uplo::Symbol, Val{false}) -> Cholesky

.. Docstring generated from Julia source

The same as ``cholfact``\ , but saves space by overwriting the input ``A``\ , instead of creating a copy. An ``InexactError`` exception is thrown if the factorisation produces a number not representable by the element type of ``A``\ , e.g. for integer types.

.. function:: cholfact!(A::StridedMatrix, uplo::Symbol, Val{true}) -> PivotedCholesky
.. function:: cholfact!(A, uplo::Symbol, Val{true}; tol = 0.0) -> CholeskyPivoted

.. Docstring generated from Julia source

Expand Down
Loading