Skip to content

Commit

Permalink
deprecate cholfact[!](A, uplo)
Browse files Browse the repository at this point in the history
restructure the cholesky test
  • Loading branch information
fredrikekre committed Jun 2, 2017
1 parent ef9ab60 commit 1ae54ff
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 70 deletions.
3 changes: 3 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,9 @@ Deprecated or removed

* The method `srand(rng, filename, n=4)` has been deprecated ([#21359]).

* The `cholfact`/`cholfact!` methods that accepted an `uplo` symbol have been deprecated
in favor of using `Hermitian` (or `Symmetric`) views ([#22187], [#22188]).


Julia v0.6.0 Release Notes
==========================
Expand Down
8 changes: 8 additions & 0 deletions base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1350,6 +1350,14 @@ end
@deprecate srand(filename::AbstractString, n::Integer=4) srand(read!(filename, Array{UInt32}(Int(n))))
@deprecate MersenneTwister(filename::AbstractString) srand(MersenneTwister(0), read!(filename, Array{UInt32}(Int(4))))

# PR #22188
@deprecate cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) cholfact!(Hermitian(A, uplo), Val{false})
@deprecate cholfact!(A::StridedMatrix, uplo::Symbol) cholfact!(Hermitian(A, uplo))
@deprecate cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) cholfact(Hermitian(A, uplo), Val{false})
@deprecate cholfact(A::StridedMatrix, uplo::Symbol) cholfact(Hermitian(A, uplo))
@deprecate cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) cholfact!(Hermitian(A, uplo), Val{true}, tol = tol)
@deprecate cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0) cholfact(Hermitian(A, uplo), Val{true}, tol = tol)

# END 0.7 deprecations

# BEGIN 1.0 deprecations
Expand Down
61 changes: 19 additions & 42 deletions base/linalg/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
# The internal structure is as follows
# - _chol! returns the factor and info without checking positive definiteness
# - chol/chol! returns the factor and checks for positive definiteness
# - cholfact/cholfact! returns Cholesky with checking positive definiteness
# - cholfact/cholfact! returns Cholesky without checking positive definiteness

# 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
Expand Down Expand Up @@ -207,8 +207,8 @@ chol(x::Number, args...) = ((C, info) = _chol!(x, nothing); @assertposdef C info

# cholfact!. Destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting
function cholfact!(A::RealHermSymComplexHerm, ::Type{Val{false}})
## No pivoting (default)
function cholfact!(A::RealHermSymComplexHerm, ::Type{Val{false}}=Val{false})
if A.uplo == 'U'
CU, info = _chol!(A.data, UpperTriangular)
Cholesky(CU.data, 'U', info)
Expand All @@ -220,7 +220,7 @@ end

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact!(A, [uplo::Symbol,] Val{false}) -> Cholesky
cholfact!(A, Val{false}) -> Cholesky
The same as [`cholfact`](@ref), but saves space by overwriting the input `A`,
instead of creating a copy. An [`InexactError`](@ref) exception is thrown if
Expand All @@ -239,18 +239,9 @@ julia> cholfact!(A)
ERROR: InexactError()
```
"""
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}})
function cholfact!(A::StridedMatrix, ::Type{Val{false}}=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::RealHermSymComplexHerm) = cholfact!(A, 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))
return cholfact!(Hermitian(A), Val{false})
end


Expand All @@ -270,37 +261,34 @@ cholfact!(A::RealHermSymComplexHerm{<:Real}, ::Type{Val{true}};

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact!(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted
cholfact!(A, Val{true}; tol = 0.0) -> CholeskyPivoted
The same as [`cholfact`](@ref), but saves space by overwriting the input `A`,
instead of creating a copy. An [`InexactError`](@ref) exception is thrown if the
factorization produces a number not representable by the element type of `A`,
e.g. for integer types.
"""
function cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0)
function cholfact!(A::StridedMatrix, ::Type{Val{true}}; tol = 0.0)
ishermitian(A) || non_hermitian_error("cholfact!")
return cholfact!(Hermitian(A, uplo), Val{true}; tol = tol)
return cholfact!(Hermitian(A), Val{true}; tol = tol)
end

# cholfact. Non-destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting
cholfact(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}, ::Type{Val{false}}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)), Val{false})
## No pivoting (default)
cholfact(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}, ::Type{Val{false}}=Val{false}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)))

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact(A, [uplo::Symbol,] Val{false}) -> Cholesky
cholfact(A, Val{false}) -> Cholesky
Compute the Cholesky factorization of a dense symmetric positive definite matrix `A`
and return a `Cholesky` factorization. The matrix `A` can either be a [`Symmetric`](@ref) or [`Hermitian`](@ref)
`StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`. In the latter case,
the optional argument `uplo` may be `:L` for using the lower part or `:U` for the upper part of `A`.
The default is to use `:U`.
`StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`.
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`](@ref), [`\\`](@ref),
[`inv`](@ref), and [`det`](@ref).
A `PosDefException` exception is thrown in case the matrix is not positive definite.
# Example
Expand Down Expand Up @@ -331,18 +319,9 @@ julia> C[:L] * C[:U] == A
true
```
"""
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::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = cholfact(A, Val{false})

#### for StridedMatrices, check that matrix is symmetric/Hermitian
function cholfact(A::StridedMatrix, uplo::Symbol = :U)
function cholfact(A::StridedMatrix, ::Type{Val{false}}=Val{false})
ishermitian(A) || non_hermitian_error("cholfact")
return cholfact(Hermitian(A, uplo))
return cholfact(Hermitian(A))
end


Expand All @@ -353,22 +332,20 @@ cholfact(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}, ::Type{Val{true}}; t

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholfact(A, [uplo::Symbol,] Val{true}; tol = 0.0) -> CholeskyPivoted
cholfact(A, 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. The matrix `A` can either be a [`Symmetric`](@ref)
or [`Hermitian`](@ref) `StridedMatrix` or a *perfectly* symmetric or Hermitian `StridedMatrix`.
In the latter case, the optional argument `uplo` 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 `PivotedCholesky` objects:
[`size`](@ref), [`\\`](@ref), [`inv`](@ref), [`det`](@ref), and [`rank`](@ref).
The argument `tol` determines the tolerance for determining the rank.
For negative values, the tolerance is the machine precision.
"""
function cholfact(A::StridedMatrix, uplo::Symbol, ::Type{Val{true}}; tol = 0.0)
function cholfact(A::StridedMatrix, ::Type{Val{true}}; tol = 0.0)
ishermitian(A) || non_hermitian_error("cholfact")
return cholfact(Hermitian(A, uplo), Val{true}; tol = tol)
return cholfact(Hermitian(A), Val{true}; tol = tol)
end

## Number
Expand Down
54 changes: 26 additions & 28 deletions test/linalg/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,11 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
apd = a'*a # symmetric positive-definite

apds = Symmetric(apd)
apdsL = Symmetric(apd, :L)
apdh = Hermitian(apd)
apdhL = Hermitian(apd, :L)
ε = εa = eps(abs(float(one(eltya))))

# Test of symmetric pos. def. strided matrix
apd = a'*a
@inferred cholfact(apd)
@inferred chol(apd)
capd = factorize(apd)
Expand Down Expand Up @@ -61,6 +58,11 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
@test all(x -> x apos, cholfact(apos).factors)
@test_throws PosDefException chol(-one(eltya))

# Test cholfact with Symmetric/Hermitian upper/lower
apds = Symmetric(apd)
apdsL = Symmetric(apd, :L)
apdh = Hermitian(apd)
apdhL = Hermitian(apd, :L)
if eltya <: Real
capds = cholfact(apds)
@test inv(capds)*apds eye(n)
Expand All @@ -82,9 +84,6 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
capdh = cholfact!(copy(apd))
@test inv(capdh)*apdh eye(n)
@test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
capdh = cholfact!(copy(apd), :L)
@test inv(capdh)*apdh eye(n)
@test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
ulstring = sprint(show,capdh[:UL])
@test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring"
end
Expand All @@ -94,18 +93,12 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
U = Bidiagonal([2,sqrt(eltya(3))],[-1],true) / sqrt(eltya(2))
@test full(chol(S)) full(U)

#lower Cholesky factor
lapd = cholfact(apd, :L)
@test full(lapd) apd
l = lapd[:L]
@test l*l' apd
@test triu(capd.factors) lapd[:U]
@test tril(lapd.factors) capd[:L]
if eltya <: Real
capds = cholfact(apds)
lapds = cholfact(apdsL)
cl = chol(apdsL)
ls = lapds[:L]
@test full(capds) full(lapds) apd
@test ls*ls' apd
@test triu(capds.factors) lapds[:U]
@test tril(lapds.factors) capds[:L]
Expand All @@ -117,6 +110,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
lapdh = cholfact(apdhL)
cl = chol(apdhL)
ls = lapdh[:L]
@test full(capdh) full(lapdh) apd
@test ls*ls' apd
@test triu(capdh.factors) lapdh[:U]
@test tril(lapdh.factors) capdh[:L]
Expand All @@ -127,9 +121,9 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException

#pivoted upper Cholesky
if eltya != BigFloat
cz = cholfact(zeros(eltya,n,n), :U, Val{true})
cz = cholfact(Hermitian(zeros(eltya,n,n)), Val{true})
@test_throws Base.LinAlg.RankDeficientException Base.LinAlg.chkfullrank(cz)
cpapd = cholfact(apd, :U, Val{true})
cpapd = cholfact(apdh, Val{true})
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing
if isreal(apd)
Expand Down Expand Up @@ -162,17 +156,19 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
@test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit

if eltya != BigFloat && eltyb != BigFloat
lapd = cholfact(apdhL)
@test norm(apd * (lapd\b) - b)/norm(b) <= ε*κ*n
@test norm(apd * (lapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n
end
@test_throws DimensionMismatch lapd\RowVector(ones(n))
@test_throws DimensionMismatch cholfact(apdhL)\RowVector(ones(n))

if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted Cholesky decomposition in julia

cpapd = cholfact(apdh, Val{true})
@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

lpapd = cholfact(apd, :L, Val{true})
lpapd = cholfact(apdhL, Val{true})
@test norm(apd * (lpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
@test norm(apd * (lpapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n

Expand Down Expand Up @@ -214,8 +210,8 @@ end
AcA = A'A
BcB = AcA + v*v'
BcB = (BcB + BcB')/2
F = cholfact(AcA, uplo)
G = cholfact(BcB, uplo)
F = cholfact(Hermitian(AcA, uplo))
G = cholfact(Hermitian(BcB, uplo))
@test LinAlg.lowrankupdate(F, v)[uplo] G[uplo]
@test_throws DimensionMismatch LinAlg.lowrankupdate(F, ones(eltype(v), length(v)+1))
@test LinAlg.lowrankdowndate(G, v)[uplo] F[uplo]
Expand Down Expand Up @@ -244,7 +240,7 @@ end
0.25336108035924787 + 0.975317836492159im 0.0628393808469436 - 0.1253397353973715im
0.11192755545114 - 0.1603741874112385im 0.8439562576196216 + 1.0850814110398734im
-1.0568488936791578 - 0.06025820467086475im 0.12696236014017806 - 0.09853584666755086im]
cholfact(apd, :L, Val{true}) \ b
cholfact(Hermitian(apd, :L), Val{true}) \ b
r = factorize(apd)[:U]
E = abs.(apd - r'*r)
ε = eps(abs(float(one(Complex64))))
Expand All @@ -255,12 +251,14 @@ end
end

@testset "throw if non-Hermitian" begin
@test_throws ArgumentError cholfact(randn(5,5))
@test_throws ArgumentError cholfact(complex.(randn(5,5), randn(5,5)))
@test_throws ArgumentError Base.LinAlg.chol!(randn(5,5))
@test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{false})
@test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{true})
@test_throws ArgumentError cholfact(randn(5,5),:U,Val{false})
R = randn(5, 5)
C = complex.(R, R)
for A in (R, C)
@test_throws ArgumentError cholfact(A)
@test_throws ArgumentError cholfact!(copy(A))
@test_throws ArgumentError chol(A)
@test_throws ArgumentError Base.LinAlg.chol!(copy(A))
end
end

@testset "fail for non-BLAS element types" begin
Expand Down

0 comments on commit 1ae54ff

Please sign in to comment.