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 ea755f5
Show file tree
Hide file tree
Showing 4 changed files with 182 additions and 61 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
4 changes: 2 additions & 2 deletions doc/stdlib/linalg.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}) -> CholeskyPivoted

.. Docstring generated from Julia source
Expand Down
14 changes: 9 additions & 5 deletions test/linalg/cholesky.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
@test norm(apd * (lapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n
end

debug && println("pivoted Choleksy decomposition")
debug && println("pivoted Cholesky decomposition")
if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted cholesky decomposition in julia

@test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
Expand All @@ -127,8 +127,8 @@ begin
# Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices

X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3]
L = full(Base.LinAlg.chol!(X*X', LowerTriangular))
U = full(Base.LinAlg.chol!(X*X', UpperTriangular))
L = full(Base.LinAlg._chol!(X*X', LowerTriangular))
U = full(Base.LinAlg._chol!(X*X', UpperTriangular))
XX = full(X*X')

@test sum(sum(norm, L*L' - XX)) < eps()
Expand All @@ -143,8 +143,8 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
A = randn(5,5)
end
A = convert(Matrix{elty}, A'A)
@test_approx_eq full(cholfact(A)[:L]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular))
@test_approx_eq full(cholfact(A)[:U]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular))
@test_approx_eq full(cholfact(A)[:L]) full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular))
@test_approx_eq full(cholfact(A)[:U]) full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular))
end

# Test up- and downdates
Expand Down Expand Up @@ -187,3 +187,7 @@ let apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.
@test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j]))
end
end

# Fail if non-Hermitian
@test_throws ArgumentError cholfact(randn(5,5))
@test_throws ArgumentError cholfact(complex(randn(5,5), randn(5,5)))

0 comments on commit ea755f5

Please sign in to comment.