Skip to content

Commit

Permalink
Replace Val-types in cholesky[!] (#41640)
Browse files Browse the repository at this point in the history
  • Loading branch information
dkarrasch authored Dec 14, 2021
1 parent afdd6bd commit 490d78c
Show file tree
Hide file tree
Showing 5 changed files with 64 additions and 51 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,10 @@ Standard library changes
#### LinearAlgebra
* The BLAS submodule now supports the level-2 BLAS subroutine `spr!` ([#42830]).

* `cholesky[!]` now supports `LinearAlgebra.PivotingStrategy` (singleton type) values
as its optional `pivot` argument: the default is `cholesky(A, NoPivot())` (vs.
`cholesky(A, RowMaximum())`); the former `Val{true/false}`-based calls are deprecated. ([#41640])

#### Markdown

#### Printf
Expand Down
49 changes: 28 additions & 21 deletions stdlib/LinearAlgebra/src/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,10 @@
# In the methods below, LAPACK is called when possible, i.e. StridedMatrices with Float32,
# Float64, ComplexF32, and ComplexF64 element types. For other element or
# matrix types, the unblocked Julia implementation in _chol! is used. For cholesky
# and cholesky! pivoting is supported through a Val(Bool) argument. A type argument is
# and cholesky! pivoting is supported through a RowMaximum() argument. A type argument is
# necessary for type stability since the output of cholesky and cholesky! is either
# Cholesky or CholeskyPivoted. The latter is only
# supported for the four LAPACK element types. For other types, e.g. BigFloats Val(true) will
# supported for the four LAPACK element types. For other types, e.g. BigFloats RowMaximum() will
# give an error. It is required that the input is Hermitian (including real symmetric) either
# through the Hermitian and Symmetric views or exact symmetric or Hermitian elements which
# is checked for and an error is thrown if the check fails.
Expand Down Expand Up @@ -106,7 +106,7 @@ Base.iterate(C::Cholesky, ::Val{:done}) = nothing
CholeskyPivoted
Matrix factorization type of the pivoted Cholesky factorization of a dense symmetric/Hermitian
positive semi-definite matrix `A`. This is the return type of [`cholesky(_, Val(true))`](@ref),
positive semi-definite matrix `A`. This is the return type of [`cholesky(_, ::RowMaximum)`](@ref),
the corresponding matrix factorization function.
The triangular Cholesky factor can be obtained from the factorization `F::CholeskyPivoted`
Expand All @@ -126,7 +126,7 @@ julia> X = [1.0, 2.0, 3.0, 4.0];
julia> A = X * X';
julia> C = cholesky(A, Val(true), check = false)
julia> C = cholesky(A, RowMaximum(), check = false)
CholeskyPivoted{Float64, Matrix{Float64}}
U factor with rank 1:
4×4 UpperTriangular{Float64, Matrix{Float64}}:
Expand Down Expand Up @@ -261,15 +261,15 @@ end
# cholesky!. Destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting (default)
function cholesky!(A::RealHermSymComplexHerm, ::Val{false}=Val(false); check::Bool = true)
function cholesky!(A::RealHermSymComplexHerm, ::NoPivot = NoPivot(); check::Bool = true)
C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
check && checkpositivedefinite(info)
return Cholesky(C.data, A.uplo, info)
end

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholesky!(A::StridedMatrix, Val(false); check = true) -> Cholesky
cholesky!(A::StridedMatrix, NoPivot(); check = true) -> Cholesky
The same as [`cholesky`](@ref), but saves space by overwriting the input `A`,
instead of creating a copy. An [`InexactError`](@ref) exception is thrown if
Expand All @@ -289,57 +289,62 @@ Stacktrace:
[...]
```
"""
function cholesky!(A::StridedMatrix, ::Val{false}=Val(false); check::Bool = true)
function cholesky!(A::StridedMatrix, ::NoPivot = NoPivot(); check::Bool = true)
checksquare(A)
if !ishermitian(A) # return with info = -1 if not Hermitian
check && checkpositivedefinite(-1)
return Cholesky(A, 'U', convert(BlasInt, -1))
else
return cholesky!(Hermitian(A), Val(false); check = check)
return cholesky!(Hermitian(A), NoPivot(); check = check)
end
end
@deprecate cholesky!(A::StridedMatrix, ::Val{false}; check::Bool = true) cholesky!(A, NoPivot(); check) false
@deprecate cholesky!(A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true) cholesky!(A, NoPivot(); check) false

## With pivoting
### BLAS/LAPACK element types
function cholesky!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix},
::Val{true}; tol = 0.0, check::Bool = true)
::RowMaximum; tol = 0.0, check::Bool = true)
AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
check && chkfullrank(C)
return C
end
@deprecate cholesky!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, ::Val{true}; kwargs...) cholesky!(A, RowMaximum(); kwargs...) false

### Non BLAS/LAPACK element types (generic). Since generic fallback for pivoted Cholesky
### is not implemented yet we throw an error
cholesky!(A::RealHermSymComplexHerm{<:Real}, ::Val{true}; tol = 0.0, check::Bool = true) =
cholesky!(A::RealHermSymComplexHerm{<:Real}, ::RowMaximum; tol = 0.0, check::Bool = true) =
throw(ArgumentError("generic pivoted Cholesky factorization is not implemented yet"))
@deprecate cholesky!(A::RealHermSymComplexHerm{<:Real}, ::Val{true}; kwargs...) cholesky!(A, RowMaximum(); kwargs...) false

### for StridedMatrices, check that matrix is symmetric/Hermitian
"""
cholesky!(A::StridedMatrix, Val(true); tol = 0.0, check = true) -> CholeskyPivoted
cholesky!(A::StridedMatrix, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivoted
The same as [`cholesky`](@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 cholesky!(A::StridedMatrix, ::Val{true}; tol = 0.0, check::Bool = true)
function cholesky!(A::StridedMatrix, ::RowMaximum; tol = 0.0, check::Bool = true)
checksquare(A)
if !ishermitian(A)
C = CholeskyPivoted(A, 'U', Vector{BlasInt}(),convert(BlasInt, 1),
tol, convert(BlasInt, -1))
check && chkfullrank(C)
return C
else
return cholesky!(Hermitian(A), Val(true); tol = tol, check = check)
return cholesky!(Hermitian(A), RowMaximum(); tol = tol, check = check)
end
end
@deprecate cholesky!(A::StridedMatrix, ::Val{true}; kwargs...) cholesky!(A, RowMaximum(); kwargs...) false

# cholesky. Non-destructive methods for computing Cholesky factorization of real symmetric
# or Hermitian matrix
## No pivoting (default)
"""
cholesky(A, Val(false); check = true) -> Cholesky
cholesky(A, NoPivot(); check = true) -> Cholesky
Compute the Cholesky factorization of a dense symmetric positive definite matrix `A`
and return a [`Cholesky`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref) or [`Hermitian`](@ref)
Expand Down Expand Up @@ -391,17 +396,18 @@ true
```
"""
cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}},
::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
::NoPivot=NoPivot(); check::Bool = true) = cholesky!(cholcopy(A); check = check)
@deprecate cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}}, ::Val{false}; check::Bool = true) cholesky(A, NoPivot(); check) false

function cholesky(A::Union{StridedMatrix{Float16},RealHermSymComplexHerm{Float16,<:StridedMatrix}}, ::Val{false}=Val(false); check::Bool = true)
function cholesky(A::Union{StridedMatrix{Float16},RealHermSymComplexHerm{Float16,<:StridedMatrix}}, ::NoPivot=NoPivot(); check::Bool = true)
X = cholesky!(cholcopy(A); check = check)
return Cholesky{Float16}(X)
end

@deprecate cholesky(A::Union{StridedMatrix{Float16},RealHermSymComplexHerm{Float16,<:StridedMatrix}}, ::Val{false}; check::Bool = true) cholesky(A, NoPivot(); check) false

## With pivoting
"""
cholesky(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted
cholesky(A, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivoted
Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix `A`
and return a [`CholeskyPivoted`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref)
Expand Down Expand Up @@ -431,7 +437,7 @@ julia> X = [1.0, 2.0, 3.0, 4.0];
julia> A = X * X';
julia> C = cholesky(A, Val(true), check = false)
julia> C = cholesky(A, RowMaximum(), check = false)
CholeskyPivoted{Float64, Matrix{Float64}}
U factor with rank 1:
4×4 UpperTriangular{Float64, Matrix{Float64}}:
Expand All @@ -456,8 +462,9 @@ true
```
"""
cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}},
::Val{true}; tol = 0.0, check::Bool = true) =
cholesky!(cholcopy(A), Val(true); tol = tol, check = check)
::RowMaximum; tol = 0.0, check::Bool = true) =
cholesky!(cholcopy(A), RowMaximum(); tol = tol, check = check)
@deprecate cholesky(A::Union{StridedMatrix,RealHermSymComplexHerm{<:Real,<:StridedMatrix}}, ::Val{true}; tol = 0.0, check::Bool = true) cholesky(A, RowMaximum(); tol, check) false

## Number
function cholesky(x::Number, uplo::Symbol=:U)
Expand Down
8 changes: 5 additions & 3 deletions stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -715,7 +715,7 @@ function _mapreduce_prod(f, x, D::Diagonal, y)
end
end

function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
function cholesky!(A::Diagonal, ::NoPivot = NoPivot(); check::Bool = true)
info = 0
for (i, di) in enumerate(A.diag)
if isreal(di) && real(di) > 0
Expand All @@ -729,9 +729,11 @@ function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
end
Cholesky(A, 'U', convert(BlasInt, info))
end
@deprecate cholesky!(A::Diagonal, ::Val{false}; check::Bool = true) cholesky!(A::Diagonal, NoPivot(); check) false

cholesky(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) =
cholesky!(cholcopy(A), Val(false); check = check)
cholesky(A::Diagonal, ::NoPivot = NoPivot(); check::Bool = true) =
cholesky!(cholcopy(A), NoPivot(); check = check)
@deprecate cholesky(A::Diagonal, ::Val{false}; check::Bool = true) cholesky(A::Diagonal, NoPivot(); check) false

function getproperty(C::Cholesky{<:Any,<:Diagonal}, d::Symbol)
Cfactors = getfield(C, :factors)
Expand Down
50 changes: 25 additions & 25 deletions stdlib/LinearAlgebra/test/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ end

#pivoted upper Cholesky
if eltya != BigFloat
cpapd = cholesky(apdh, Val(true))
cpapd = cholesky(apdh, RowMaximum())
unary_ops_tests(apdh, cpapd, ε*κ*n)
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing
Expand Down Expand Up @@ -167,11 +167,11 @@ end

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

cpapd = cholesky(apdh, Val(true))
cpapd = cholesky(apdh, RowMaximum())
@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 = cholesky(apdhL, Val(true))
lpapd = cholesky(apdhL, RowMaximum())
@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
end
Expand All @@ -194,7 +194,7 @@ end
@test norm(apd \ B - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test norm(apd * BB - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
if eltya != BigFloat
cpapd = cholesky(apdh, Val(true))
cpapd = cholesky(apdh, RowMaximum())
BB = copy(B)
ldiv!(cpapd, BB)
@test norm(apd \ B - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
Expand Down Expand Up @@ -225,12 +225,12 @@ end
@test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
if eltya != BigFloat
cpapd = cholesky(eltya <: Complex ? apdh : apds, Val(true))
cpapd = cholesky(eltya <: Complex ? apdh : apds, RowMaximum())
BB = copy(B)
rdiv!(BB, cpapd)
@test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
cpapd = cholesky(eltya <: Complex ? apdhL : apdsL, Val(true))
cpapd = cholesky(eltya <: Complex ? apdhL : apdsL, RowMaximum())
BB = copy(B)
rdiv!(BB, cpapd)
@test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
Expand Down Expand Up @@ -266,15 +266,15 @@ end
@test !LinearAlgebra.issuccess(cholesky!(copy(M); check = false))
end
for M in (A, Hermitian(A), B)
@test_throws RankDeficientException cholesky(M, Val(true))
@test_throws RankDeficientException cholesky!(copy(M), Val(true))
@test_throws RankDeficientException cholesky(M, Val(true); check = true)
@test_throws RankDeficientException cholesky!(copy(M), Val(true); check = true)
@test !LinearAlgebra.issuccess(cholesky(M, Val(true); check = false))
@test !LinearAlgebra.issuccess(cholesky!(copy(M), Val(true); check = false))
C = cholesky(M, Val(true); check = false)
@test_throws RankDeficientException cholesky(M, RowMaximum())
@test_throws RankDeficientException cholesky!(copy(M), RowMaximum())
@test_throws RankDeficientException cholesky(M, RowMaximum(); check = true)
@test_throws RankDeficientException cholesky!(copy(M), RowMaximum(); check = true)
@test !LinearAlgebra.issuccess(cholesky(M, RowMaximum(); check = false))
@test !LinearAlgebra.issuccess(cholesky!(copy(M), RowMaximum(); check = false))
C = cholesky(M, RowMaximum(); check = false)
@test_throws RankDeficientException chkfullrank(C)
C = cholesky!(copy(M), Val(true); check = false)
C = cholesky!(copy(M), RowMaximum(); check = false)
@test_throws RankDeficientException chkfullrank(C)
end
@test !isposdef(A)
Expand Down Expand Up @@ -339,7 +339,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]
cholesky(Hermitian(apd, :L), Val(true)) \ b
cholesky(Hermitian(apd, :L), RowMaximum()) \ b
r = factorize(apd).U
E = abs.(apd - r'*r)
ε = eps(abs(float(one(ComplexF32))))
Expand All @@ -350,7 +350,7 @@ end
end

@testset "fail for non-BLAS element types" begin
@test_throws ArgumentError cholesky!(Hermitian(rand(Float16, 5,5)), Val(true))
@test_throws ArgumentError cholesky!(Hermitian(rand(Float16, 5,5)), RowMaximum())
end

@testset "cholesky Diagonal" begin
Expand Down Expand Up @@ -398,7 +398,7 @@ end
@test Cholesky(factors, uplo, Int32(info)) == chol
@test Cholesky(factors, uplo, Int64(info)) == chol

cholp = cholesky(x'x, Val(true))
cholp = cholesky(x'x, RowMaximum())

factors, uplo, piv, rank, tol, info =
cholp.factors, cholp.uplo, cholp.piv, cholp.rank, cholp.tol, cholp.info
Expand All @@ -417,25 +417,25 @@ end
@testset "issue #33704, casting low-rank CholeskyPivoted to Matrix" begin
A = randn(1,8)
B = A'A
C = cholesky(B, Val(true), check=false)
C = cholesky(B, RowMaximum(), check=false)
@test B Matrix(C)
end

@testset "CholeskyPivoted and Factorization" begin
A = randn(8,8)
B = A'A
C = cholesky(B, Val(true), check=false)
C = cholesky(B, RowMaximum(), check=false)
@test CholeskyPivoted{eltype(C)}(C) === C
@test Factorization{eltype(C)}(C) === C
@test Array(CholeskyPivoted{complex(eltype(C))}(C)) Array(cholesky(complex(B), Val(true), check=false))
@test Array(Factorization{complex(eltype(C))}(C)) Array(cholesky(complex(B), Val(true), check=false))
@test Array(CholeskyPivoted{complex(eltype(C))}(C)) Array(cholesky(complex(B), RowMaximum(), check=false))
@test Array(Factorization{complex(eltype(C))}(C)) Array(cholesky(complex(B), RowMaximum(), check=false))
@test eltype(Factorization{complex(eltype(C))}(C)) == complex(eltype(C))
end

@testset "REPL printing of CholeskyPivoted" begin
A = randn(8,8)
B = A'A
C = cholesky(B, Val(true), check=false)
C = cholesky(B, RowMaximum(), check=false)
cholstring = sprint((t, s) -> show(t, "text/plain", s), C)
rankstring = "$(C.uplo) factor with rank $(rank(C)):"
factorstring = sprint((t, s) -> show(t, "text/plain", s), C.uplo == 'U' ? C.U : C.L)
Expand All @@ -444,10 +444,10 @@ end
end

@testset "destructuring for Cholesky[Pivoted]" begin
for val in (true, false)
for val in (NoPivot(), RowMaximum())
A = rand(8, 8)
B = A'A
C = cholesky(B, Val(val), check=false)
C = cholesky(B, val, check=false)
l, u = C
@test l == C.L
@test u == C.U
Expand Down Expand Up @@ -507,7 +507,7 @@ end
2048 1920 2940 1008 2240 2740;
4470 4200 6410 2240 4875 6015;
5490 5140 7903 2740 6015 7370]
B = cholesky(A, Val(true), check=false)
B = cholesky(A, RowMaximum(), check=false)
@test det(B) == 0.0
@test det(B) det(A) atol=eps()
@test logdet(B) == -Inf
Expand Down
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/test/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using Test, LinearAlgebra
@testset "equality for factorizations - $f" for f in Any[
bunchkaufman,
cholesky,
x -> cholesky(x, Val(true)),
x -> cholesky(x, RowMaximum()),
eigen,
hessenberg,
lq,
Expand Down Expand Up @@ -45,7 +45,7 @@ end
@testset "size for factorizations - $f" for f in Any[
bunchkaufman,
cholesky,
x -> cholesky(x, Val(true)),
x -> cholesky(x, RowMaximum()),
hessenberg,
lq,
lu,
Expand Down

0 comments on commit 490d78c

Please sign in to comment.