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]