Skip to content

Commit

Permalink
Check that matrix is symmetric/Hermitian in dense Cholesky
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasnoack
andreasnoack committed Jun 7, 2016
1 parent 44d778a commit 2d5d3f9
Show file tree
Hide file tree
Showing 4 changed files with 188 additions and 66 deletions.
3 changes: 1 addition & 2 deletions base/linalg.jl
View file
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
222 changes: 170 additions & 52 deletions base/linalg/cholesky.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,29 @@
##########################
# 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 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 then checked 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 simplifacation could be to merge
# the two Cholesky types into a single. It would remove the need for Val completely but
# the cost would be extra unnecessary/ununsed fields for the unpivoted Cholesky and runtime # 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 +46,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 in LAPACK
## 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 +79,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
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 +104,149 @@ 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



# chol!. Destructive methods for computing Cholesky factor of real symmetric of 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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling chol!(Hermitian(A)) instead."))
return _chol!(A, UpperTriangular)
end



# chol. Non-destructive methods for computing Cholesky factor of real symmetric of
# Hermitian matrix. Promotes elements to type that is stable under sqrt 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::StridedMatrix)
ishermitian(A) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling chol(Hermitian(A)) instead."))
return chol(Hermitian(A))
end

## Numbers
"""
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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact!(Hermitian(A)) instead."))
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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact!(Hermitian(A)) instead."))
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 fectorization 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}) -> 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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact!(Hermitian(A)) instead."))
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 +256,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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact(Hermitian(A)) instead."))
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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact(Hermitian(A)) instead."))
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 +294,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) || throw(ArgumentError("matrix is not symmetric/Hermitian. Consider calling cholfact(Hermitian(A)) instead."))
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
15 changes: 8 additions & 7 deletions doc/stdlib/linalg.rst
View file
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,17 +169,17 @@ 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.
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}; tol = 0.0) -> CholeskyPivoted
.. function:: cholfact!(A, uplo::Symbol, Val{true}) -> CholeskyPivoted

.. Docstring generated from Julia source
Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix ``A`` and return a ``CholeskyPivoted`` 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 ``PivotedCholesky`` objects: ``size``\ , ``\``\ , ``inv``\ , ``det``\ , and ``rank``\ . The argument ``tol`` determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision.
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(::Union{SparseMatrixCSC{<:Real},SparseMatrixCSC{Complex{<:Real}},Symmetric{<:Real,SparseMatrixCSC{<:Real,SuiteSparse_long}},Hermitian{Complex{<:Real},SparseMatrixCSC{Complex{<:Real},SuiteSparse_long}}}; shift = 0.0, perm = Int[]) -> CHOLMOD.Factor

Expand All @@ -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}) -> CholeskyPivoted

.. Docstring generated from Julia source
Expand Down Expand Up @@ -2212,3 +2212,4 @@ set of functions in future releases.
Solves the Sylvester matrix equation ``A * X +/- X * B = scale*C`` where ``A`` and ``B`` are both quasi-upper triangular. If ``transa = N``\ , ``A`` is not modified. If ``transa = T``\ , ``A`` is transposed. If ``transa = C``\ , ``A`` is conjugate transposed. Similarly for ``transb`` and ``B``\ . If ``isgn = 1``\ , the equation ``A * X + X * B = scale * C`` is solved. If ``isgn = -1``\ , the equation ``A * X - X * B = scale * C`` is solved.

Returns ``X`` (overwriting ``C``\ ) and ``scale``\ .

Loading

0 comments on commit 2d5d3f9

Please sign in to comment.