From 5e880747d8c0e444863a6676a4ea43c758858999 Mon Sep 17 00:00:00 2001 From: Andreas Noack Jensen Date: Thu, 12 Jun 2014 22:54:38 +0200 Subject: [PATCH] Add generic Cholesky decomposition and make Cholesky parametric on matrix type plus some fixes. Define tril! and triu! for AbstractMatrix and use elements instead of type to construct zeros. Also use elements for constructing zeros in getindex(Triangular). --- base/linalg.jl | 1 + base/linalg/cholesky.jl | 181 +++++++++++++++++++++++++++++++++++ base/linalg/dense.jl | 8 +- base/linalg/factorization.jl | 120 ----------------------- base/linalg/triangular.jl | 10 +- test/linalg1.jl | 55 ++++++----- test/linalg4.jl | 13 +++ 7 files changed, 231 insertions(+), 157 deletions(-) create mode 100644 base/linalg/cholesky.jl diff --git a/base/linalg.jl b/base/linalg.jl index 1b0e422a84b5c..560e4325acca6 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -198,6 +198,7 @@ include("linalg/dense.jl") include("linalg/tridiag.jl") include("linalg/triangular.jl") include("linalg/factorization.jl") +include("linalg/cholesky.jl") include("linalg/lu.jl") include("linalg/bunchkaufman.jl") diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl new file mode 100644 index 0000000000000..d9a83d1a1e9df --- /dev/null +++ b/base/linalg/cholesky.jl @@ -0,0 +1,181 @@ +########################## +# Cholesky Factorization # +########################## +immutable Cholesky{T,S<:AbstractMatrix{T},UpLo} <: Factorization{T} + UL::S +end +immutable CholeskyPivoted{T} <: Factorization{T} + UL::Matrix{T} + uplo::Char + piv::Vector{BlasInt} + rank::BlasInt + tol::Real + info::BlasInt +end + +function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U) + C, info = LAPACK.potrf!(string(uplo)[1], A) + return @assertposdef Triangular(C, uplo, false) info +end + +function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U) + n = chksquare(A) + @inbounds begin + if uplo == :L + for k = 1:n + for i = 1:k - 1 + A[k,k] -= A[k,i]*A[k,i]' + end + A[k,k] = chol!(A[k,k], uplo) + AkkInv = inv(A[k,k]') + for j = 1:k + for i = k + 1:n + j == 1 && (A[i,k] = A[i,k]*AkkInv) + j < k && (A[i,k] -= A[i,j]*A[k,j]'*AkkInv) + end + end + end + elseif uplo == :U + for k = 1:n + for i = 1:k - 1 + A[k,k] -= A[i,k]'A[i,k] + end + A[k,k] = chol!(A[k,k], uplo) + AkkInv = inv(A[k,k]) + for j = k + 1:n + for i = 1:k - 1 + A[k,j] -= A[i,k]'A[i,j] + end + A[k,j] = A[k,k]'\A[k,j] + end + end + else + throw(ArgumentError("uplo must be either :U or :L but was $(uplo)")) + end + end + return Triangular(A, uplo, false) +end + +function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) + uplochar = string(uplo)[1] + if pivot + A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol) + return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info) + end + return Cholesky{T,typeof(A),uplo}(chol!(A, uplo).data) +end +cholfact!(A::AbstractMatrix, uplo::Symbol=:U) = Cholesky{eltype(A),typeof(A),uplo}(chol!(A, uplo).data) + +cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol) +function cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) + S = promote_type(typeof(chol(one(T))),Float32) + S <: BlasFloat && return cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol) + pivot && throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}")) + S != T && return cholfact!(convert(AbstractMatrix{S}, A), uplo) + return cholfact!(copy(A), uplo) +end +function cholfact(x::Number, uplo::Symbol=:U) + xf = fill(chol!(x), 1, 1) + Cholesky{:U, eltype(xf), typeof(xf)}(xf) +end + +chol(A::AbstractMatrix, uplo::Symbol=:U) = Triangular(chol!(copy(A), uplo), uplo, false) +function chol!(x::Number, uplo::Symbol=:U) + rx = real(x) + rx == abs(x) || throw(DomainError()) + rxr = sqrt(rx) + convert(promote_type(typeof(x), typeof(rxr)), rxr) +end +chol(x::Number, uplo::Symbol=:U) = chol!(x, uplo) + +function convert{Tnew,Told,S,UpLo}(::Type{Cholesky{Tnew}},C::Cholesky{Told,S,UpLo}) + Cnew = convert(AbstractMatrix{Tnew}, C.UL) + Cholesky{Tnew, typeof(Cnew), UpLo}(Cnew) +end +function convert{T,S,UpLo}(::Type{Cholesky{T,S,UpLo}},C::Cholesky) + Cnew = convert(AbstractMatrix{T}, C.UL) + Cholesky{T, typeof(Cnew), UpLo}(Cnew) +end +convert{T}(::Type{Factorization{T}}, C::Cholesky) = convert(Cholesky{T}, C) +convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = CholeskyPivoted(convert(AbstractMatrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info) +convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C) + +full{T,S}(C::Cholesky{T,S,:U}) = C[:U]'C[:U] +full{T,S}(C::Cholesky{T,S,:L}) = C[:L]*C[:L]' + +size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL) +size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d) + +function getindex{T,S,UpLo}(C::Cholesky{T,S,UpLo}, d::Symbol) + d == :U && return Triangular(UpLo == d ? C.UL : C.UL',:U) + d == :L && return Triangular(UpLo == d ? C.UL : C.UL',:L) + d == :UL && return Triangular(C.UL, UpLo) + throw(KeyError(d)) +end +function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol) + d == :U && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :U) + d == :L && return Triangular(symbol(C.uplo) == d ? C.UL : C.UL', :L) + d == :p && return C.piv + if d == :P + n = size(C, 1) + P = zeros(T, n, n) + for i=1:n + P[C.piv[i],i] = one(T) + end + return P + end + throw(KeyError(d)) +end + +show{T,S<:AbstractMatrix,UpLo}(io::IO, C::Cholesky{T,S,UpLo}) = (println("$(typeof(C)) with factor:");show(io,C[UpLo])) + +A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat{T}) = LAPACK.potrs!('U', C.UL, B) +A_ldiv_B!{T<:BlasFloat,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat{T}) = LAPACK.potrs!('L', C.UL, B) +A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:L}, B::StridedVecOrMat) = Ac_ldiv_B!(Triangular(C.UL, :L, false), A_ldiv_B!(Triangular(C.UL, :L, false), B)) +A_ldiv_B!{T,S<:AbstractMatrix}(C::Cholesky{T,S,:U}, B::StridedVecOrMat) = A_ldiv_B!(Triangular(C.UL, :U, false), Ac_ldiv_B!(Triangular(C.UL, :U, false), B)) + +function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T}) + chkfullrank(C) + ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv) +end +function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T}) + chkfullrank(C) + n = size(C, 1) + for i=1:size(B, 2) + permute!(sub(B, 1:n, i), C.piv) + end + LAPACK.potrs!(C.uplo, C.UL, B) + for i=1:size(B, 2) + ipermute!(sub(B, 1:n, i), C.piv) + end + B +end +A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv]))[invperm(C.piv)] +A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), Ac_ldiv_B!(Triangular(C.UL, symbol(C.uplo), false), B[C.piv,:]))[invperm(C.piv),:] + +function det{T,S,UpLo}(C::Cholesky{T,S,UpLo}) + dd = one(T) + for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end + dd +end + +det{T}(C::CholeskyPivoted{T}) = C.rank 0) - -chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo] -chol(A::Union(Number, AbstractMatrix)) = triu!(cholfact(A, :U).UL) - -convert{T}(::Type{Cholesky{T}},C::Cholesky) = Cholesky(convert(AbstractMatrix{T},C.UL),C.uplo) -convert{T}(::Type{Factorization{T}}, C::Cholesky) = convert(Cholesky{T}, C) -convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = CholeskyPivoted(convert(AbstractMatrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info) -convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C) - -function full{T<:BlasFloat}(C::Cholesky{T}) - if C.uplo == 'U' - BLAS.trmm!('R', C.uplo, 'N', 'N', one(T), C.UL, tril!(C.UL')) - else - BLAS.trmm!('L', C.uplo, 'N', 'N', one(T), C.UL, triu!(C.UL')) - end -end - -size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL) -size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d) - -function getindex(C::Cholesky, d::Symbol) - d == :U && return Triangular(triu!(symbol(C.uplo) == d ? C.UL : C.UL'),:U) - d == :L && return Triangular(tril!(symbol(C.uplo) == d ? C.UL : C.UL'),:L) - d == :UL && return Triangular(C.UL, symbol(C.uplo)) - throw(KeyError(d)) -end -function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol) - d == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL') - d == :L && return tril!(symbol(C.uplo) == d ? C.UL : C.UL') - d == :p && return C.piv - if d == :P - n = size(C, 1) - P = zeros(T, n, n) - for i=1:n - P[C.piv[i],i] = one(T) - end - return P - end - throw(KeyError(d)) -end - -show(io::IO, C::Cholesky) = (println(io,"$(typeof(C)) with factor:");show(io,C[symbol(C.uplo)])) - -A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B) -A_ldiv_B!(C::Cholesky, B::StridedVecOrMat) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B)) : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B)) - -function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T}) - chkfullrank(C) - ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv) -end -function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T}) - chkfullrank(C) - n = size(C, 1) - for i=1:size(B, 2) - permute!(sub(B, 1:n, i), C.piv) - end - LAPACK.potrs!(C.uplo, C.UL, B) - for i=1:size(B, 2) - ipermute!(sub(B, 1:n, i), C.piv) - end - B -end -A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)] -A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:] - -function det{T}(C::Cholesky{T}) - dd = one(T) - for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end - dd -end - -det{T}(C::CholeskyPivoted{T}) = C.rank j ? A.data[i,j] : zero(T)) -getindex{T,S}(A::Triangular{T,S,:L,false}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(T) -getindex{T,S}(A::Triangular{T,S,:U,true}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(T)) -getindex{T,S}(A::Triangular{T,S,:U,false}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(T) -getindex{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}, i::Integer) = ((m, n) = divrem(i - 1, size(A,1)); A[m + 1, n + 1]) +getindex{T,S}(A::Triangular{T,S,:L,true}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j])) +getindex{T,S}(A::Triangular{T,S,:L,false}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j]) +getindex{T,S}(A::Triangular{T,S,:U,true}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j])) +getindex{T,S}(A::Triangular{T,S,:U,false}, i::Integer, j::Integer) = i <= j ? A.data[i,j] : zero(A.data[i,j]) +getindex(A::Triangular, i::Integer) = ((m, n) = divrem(i - 1, size(A,1)); A[m + 1, n + 1]) istril{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :L istriu{T,S,UpLo,IsUnit}(A::Triangular{T,S,UpLo,IsUnit}) = UpLo == :U diff --git a/test/linalg1.jl b/test/linalg1.jl index 2b4b76113128c..d75e903c63116 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -37,38 +37,37 @@ for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") debug && println("(Automatic) upper Cholesky factor") - if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement cholesky decomposition in julia - capd = factorize(apd) - r = capd[:U] - κ = cond(apd) #condition number - - #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 - E = abs(apd - r'*r) - for i=1:n, j=1:n - @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) - end - E = abs(apd - full(capd)) - for i=1:n, j=1:n - @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) - end + + capd = factorize(apd) + r = capd[:U] + κ = cond(apd, 1) #condition number + + #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 + E = abs(apd - r'*r) + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + E = abs(apd - full(capd)) + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end - #Test error bound on linear solver: LAWNS 14, Theorem 2.1 - #This is a surprisingly loose bound... - x = capd\b - @test norm(x-apd\b)/norm(x) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ - @test norm(apd*x-b)/norm(b) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + #Test error bound on linear solver: LAWNS 14, Theorem 2.1 + #This is a surprisingly loose bound... + x = capd\b + @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ - @test_approx_eq apd * inv(capd) eye(n) - @test norm(a*(capd\(a'*b)) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit - @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit - @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow + @test_approx_eq apd * inv(capd) eye(n) + @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit + @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit + @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow debug && println("lower Cholesky factor") - lapd = cholfact(apd, :L) - @test_approx_eq full(lapd) apd - l = lapd[:L] - @test_approx_eq l*l' apd - end + lapd = cholfact(apd, :L) + @test_approx_eq full(lapd) apd + l = lapd[:L] + @test_approx_eq l*l' apd debug && println("pivoted Choleksy decomposition") if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted cholesky decomposition in julia diff --git a/test/linalg4.jl b/test/linalg4.jl index c2d8b78c60c6a..1a9f6a9898c57 100644 --- a/test/linalg4.jl +++ b/test/linalg4.jl @@ -331,6 +331,19 @@ for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Triangular, Matrix] @test full(convert(newtype, A)) == full(A) end +# Test generic cholfact! +for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) + if elty <: Complex + A = complex(randn(5,5), randn(5,5)) + else + A = randn(5,5) + end + A = convert(Matrix{elty}, A'A) + for ul in (:U, :L) + @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, (AbstractMatrix,Symbol),copy(A), ul)) + end +end + # Issue #7886 x, r = LAPACK.gelsy!([0 1; 0 2; 0 3.], [2, 4, 6.]) @test_approx_eq x [0,2]