Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate chol(A, Val{:U/:L}) and move documentation to method definitions #13680

Merged
merged 1 commit into from
Oct 22, 2015
Merged

Deprecate chol(A, Val{:U/:L}) and move documentation to method definitions #13680

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions base/deprecated.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -858,3 +858,8 @@ end

#13496
@deprecate A_ldiv_B!(A::SparseMatrixCSC, B::StridedVecOrMat) A_ldiv_B!(factorize(A), B)

@deprecate chol(A::Number, ::Type{Val{:U}}) chol(A)
@deprecate chol(A::AbstractMatrix, ::Type{Val{:U}}) chol(A)
@deprecate chol(A::Number, ::Type{Val{:L}}) ctranspose(chol(A))
@deprecate chol(A::AbstractMatrix, ::Type{Val{:L}}) ctranspose(chol(A))
7 changes: 0 additions & 7 deletions base/docs/helpdb.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5734,13 +5734,6 @@ as ``v0``. In general, this cannot be used with empty collections
"""
foldr(op, itr)

doc"""
chol(A, [LU]) -> F

Compute the Cholesky factorization of a symmetric positive definite matrix `A` and return the matrix `F`. If `LU` is `Val{:U}` (Upper), `F` is of type `UpperTriangular` and `A = F'*F`. If `LU` is `Val{:L}` (Lower), `F` is of type `LowerTriangular` and `A = F*F'`. `LU` defaults to `Val{:U}`.
"""
chol

doc"""
ParseError(msg)

Expand Down
96 changes: 50 additions & 46 deletions base/linalg/cholesky.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,24 +23,24 @@ 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{Val{:U}})
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{Val{:L}})
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, Val{:U})
chol!(A::StridedMatrix) = chol!(A, UpperTriangular)

function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:U}})
function chol!{T}(A::AbstractMatrix{T}, ::Type{UpperTriangular})
n = chksquare(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], Val{:U})
Akk = chol!(A[k,k], UpperTriangular)
A[k,k] = Akk
AkkInv = inv(Akk')
for j = k + 1:n
Expand All @@ -53,14 +53,14 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:U}})
end
return UpperTriangular(A)
end
function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:L}})
function chol!{T}(A::AbstractMatrix{T}, ::Type{LowerTriangular})
n = chksquare(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], Val{:L})
Akk = chol!(A[k,k], LowerTriangular)
A[k,k] = Akk
AkkInv = inv(Akk)
for j = 1:k
Expand All @@ -78,14 +78,15 @@ function chol!{T}(A::AbstractMatrix{T}, ::Type{Val{:L}})
return LowerTriangular(A)
end

doc"""
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)
chol!(copy_oftype(A, S))
end
function chol{T}(A::AbstractMatrix{T}, uplo::Union{Type{Val{:L}}, Type{Val{:U}}})
S = promote_type(typeof(chol(one(T))), Float32)
chol!(copy_oftype(A, S), uplo)
end
function chol!(x::Number, uplo)
rx = real(x)
if rx != abs(x)
Expand All @@ -94,50 +95,50 @@ function chol!(x::Number, uplo)
rxr = sqrt(rx)
convert(promote_type(typeof(x), typeof(rxr)), rxr)
end
chol(x::Number, uplo::Symbol=:U) = chol!(x, uplo)

function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U,
pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{false}; tol=0.0)
_cholfact!(A, pivot, uplo, tol=tol)
end
function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U; tol=0.0)
return Cholesky(chol!(A, Val{uplo}).data, uplo)
doc"""
chol(x::Number) -> y

Compute the square root of a non-negative number `x`.
"""
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)
if uplo == :U
return Cholesky(chol!(A, UpperTriangular).data, uplo)
else
return Cholesky(chol!(A, LowerTriangular).data, uplo)
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,StridedMatrix{T}}(A, uplochar, piv, rank, tol, info)
end
function cholfact!(A::StridedMatrix, uplo::Symbol=:U,
pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{false}; tol=0.0)
Cholesky(chol!(A, Val{uplo}).data, uplo)
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}})
if uplo == :U
Cholesky(chol!(A, UpperTriangular).data, 'U')
else
Cholesky(chol!(A, UpperTriangular).data, 'U')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Whoops, was this supposed to be Cholesky(chol!(A, LowerTriangular).data, 'L') ? Definitely missing test coverage if so.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. I'm really surprised that it is not covered. I'll fix that now.

end
end
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!{T<:BlasFloat,S}(C::Cholesky{T,S})
# _, info = LAPACK.potrf!(C.uplo, C.factors)
# info[1]>0 && throw(PosDefException(info[1]))
# C
# end

# cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}; tol=0.0) =
# cholfact!(copy(A), uplo, pivot, tol=tol)


cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}; tol=0.0) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot; tol = tol)
# _cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{true}}, uplo::Symbol=:U; tol=0.0) =
# cholfact!(A, uplo, pivot, tol = tol)
# _cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{false}}, uplo::Symbol=:U; tol=0.0) =
# cholfact!(A, uplo, pivot, tol = tol)

# _cholfact{T}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U; tol=0.0) =
# cholfact!(A, uplo)
# _cholfact{T}(A::StridedMatrix{T}, ::Type{Val{true}}, uplo::Symbol=:U; tol=0.0) =
# throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}"))

cholfact{T}(A::StridedMatrix{T}, uplo::Symbol, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, Val{false})
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(x::Number, uplo::Symbol=:U)
xf = fill(chol!(x, uplo), 1, 1)
xf = fill(chol(x), 1, 1)
Cholesky(xf, uplo)
end

Expand All @@ -155,7 +156,10 @@ convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) =
convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C)

full{T,S}(C::Cholesky{T,S}) = C.uplo == 'U' ? C[:U]'C[:U] : C[:L]*C[:L]'
full(F::CholeskyPivoted) = (ip=invperm(F[:p]); (F[:L] * F[:U])[ip,ip])
function full(F::CholeskyPivoted)
ip = invperm(F[:p])
return (F[:L] * F[:U])[ip,ip]
end

copy(C::Cholesky) = Cholesky(copy(C.factors), C.uplo)
copy(C::CholeskyPivoted) = CholeskyPivoted(copy(C.factors), C.uplo, C.piv, C.rank, C.tol, C.info)
Expand Down
10 changes: 8 additions & 2 deletions doc/stdlib/linalg.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,11 +111,17 @@ 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. For sparse ``A`` the ``nzval`` field is not overwritten but the index fields, ``colptr`` and ``rowval`` are decremented in place, converting from 1-based indices to 0-based indices.

.. function:: chol(A, [LU]) -> F
.. function:: chol(A::AbstractMatrix) -> U

.. Docstring generated from Julia source

Compute the Cholesky factorization of a symmetric positive definite matrix ``A`` and return the matrix ``F``\ . If ``LU`` is ``Val{:U}`` (Upper), ``F`` is of type ``UpperTriangular`` and ``A = F'*F``\ . If ``LU`` is ``Val{:L}`` (Lower), ``F`` is of type ``LowerTriangular`` and ``A = F*F'``\ . ``LU`` defaults to ``Val{:U}``\ .
Compute the Cholesky factorization of a positive definite matrix ``A`` and return the UpperTriangular matrix ``U`` such that ``A = U'U``\ .

.. function:: chol(x::Number) -> y

.. Docstring generated from Julia source

Compute the square root of a non-negative number ``x``\ .

.. function:: cholfact(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) -> Cholesky

Expand Down
9 changes: 4 additions & 5 deletions test/linalg/cholesky.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,8 +125,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', Val{:L}))
U = full(Base.LinAlg.chol!(X*X', Val{:U}))
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 @@ -141,7 +141,6 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
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!, Tuple{AbstractMatrix, Type{Val{ul}}},copy(A), Val{ul}))
end
@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
4 changes: 2 additions & 2 deletions test/linalg/triangular.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
(UnitLowerTriangular, :L))

# Construct test matrix
A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't, Val{uplo1})))
A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't) |> t -> uplo1 == :U ? t : ctranspose(t)))

debug && println("elty1: $elty1, A1: $t1")

Expand Down Expand Up @@ -243,7 +243,7 @@ for elty1 in (Float32, Float64, Complex64, Complex128, BigFloat, Int)

debug && println("elty1: $elty1, A1: $t1, elty2: $elty2")

A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t-> chol(t't, Val{uplo2})))
A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't) |> t -> uplo2 == :U ? t : ctranspose(t)))

# Convert
if elty1 <: Real && !(elty2 <: Integer)
Expand Down