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

Jump to bottom

deprecate cholfact[!](A, uplo) #22188

Merged
merged 1 commit into from
Jun 5, 2017
+57 −70
Merged

deprecate cholfact[!](A, uplo) #22188

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
3 changes: 3 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1354,6 +1354,14 @@ end
@deprecate versioninfo(verbose::Bool) versioninfo(verbose=verbose)
@deprecate versioninfo(io::IO, verbose::Bool) versioninfo(io, verbose=verbose)

# 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
View file
Open in desktop
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
55 changes: 27 additions & 28 deletions test/linalg/cholesky.jl
View file
Open in desktop
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,13 @@ 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
Copy link
Member

Choose a reason for hiding this comment

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

Comment worth preserving/revising?

Copy link
Member Author

Choose a reason for hiding this comment

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

I thought it only applied to the tests which are now deprecated? Both Hermitian(A, :U) and Hermitian(A, :L) are tested in the tests just below.

Copy link
Member

Choose a reason for hiding this comment

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

My suggestion to preserve/revise the comment associated with that block was independent of whether that comment applies primarily or strictly to the deprecated tests: There is no other comment describing the surviving contents of that block, and having an up to date comment describing the remaining tests in that block would be worthwhile :).

Copy link
Member Author

Choose a reason for hiding this comment

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

Makes sense! I will add something descriptive :)

Copy link
Member

Choose a reason for hiding this comment

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

Awesome! Thanks! :)

Copy link
Member Author

Choose a reason for hiding this comment

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

Done! Moved around some things, I think it is clearer now.

Copy link
Member

Choose a reason for hiding this comment

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

If you mean that the comment added on line 61 should also apply to this block, perhaps move this block adjacent to the block to which that comment is attached (rather than leave those blocks separated by the block dedicated to Strang matrices)? Forty lines, paragraph breaks, and other comments interceding obscures the relationship :).

Copy link
Member Author

Choose a reason for hiding this comment

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

Hmm, I confused myself. Actually, there is no reason not to merge the two blocks, I will do that.

Copy link
Member Author

Choose a reason for hiding this comment

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

Ended up adding a comment instead.

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]
# test extraction of factor and re-creating original matrix
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 +111,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 +122,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 +157,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 +211,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 +241,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 +252,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