diff --git a/base/linalg.jl b/base/linalg.jl index 96eb19ceea0de..175632561b2ee 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -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") diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 5218d38359d4e..611df8b838004 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -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 +# 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 @@ -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 @@ -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 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 @@ -78,19 +105,8 @@ 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")) @@ -98,35 +114,80 @@ function chol!(x::Number, uplo) 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 """ 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 @@ -134,27 +195,62 @@ 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. @@ -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. @@ -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) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index ddd832be36c12..71c2594b573ea 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -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 @@ -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 @@ -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 diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 028dca616279d..3b1120a810a82 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -117,8 +117,8 @@ 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") - if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted cholesky decomposition in julia +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 @test norm(apd * (cpapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n @@ -136,8 +136,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() @@ -152,8 +152,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 @@ -196,3 +196,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)))