From e16666627f6ead25f58b7504c9f5465b1dc9007d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 14 Dec 2021 15:54:59 +0100 Subject: [PATCH] Replace Val-types in cholesky[!] (#41640) --- NEWS.md | 4 ++ stdlib/LinearAlgebra/src/cholesky.jl | 49 ++++++++++++--------- stdlib/LinearAlgebra/src/diagonal.jl | 8 ++-- stdlib/LinearAlgebra/test/cholesky.jl | 50 +++++++++++----------- stdlib/LinearAlgebra/test/factorization.jl | 4 +- 5 files changed, 64 insertions(+), 51 deletions(-) diff --git a/NEWS.md b/NEWS.md index 408dce72dd322..3ca8f0aa540c4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -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 diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl index 1941220914e6f..22e9fcfdbca83 100644 --- a/stdlib/LinearAlgebra/src/cholesky.jl +++ b/stdlib/LinearAlgebra/src/cholesky.jl @@ -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. @@ -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` @@ -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}}: @@ -261,7 +261,7 @@ 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) @@ -269,7 +269,7 @@ 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 @@ -289,41 +289,45 @@ 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), @@ -331,15 +335,16 @@ function cholesky!(A::StridedMatrix, ::Val{true}; tol = 0.0, check::Bool = true) 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) @@ -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) @@ -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}}: @@ -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) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 09562615621a1..ef01bb6cd3ba9 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -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 @@ -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) diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl index 1d258aab8ad3a..8a83ba768f2f6 100644 --- a/stdlib/LinearAlgebra/test/cholesky.jl +++ b/stdlib/LinearAlgebra/test/cholesky.jl @@ -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 @@ -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 @@ -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)*ε)*κ @@ -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)*ε)*κ @@ -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) @@ -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)))) @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/stdlib/LinearAlgebra/test/factorization.jl b/stdlib/LinearAlgebra/test/factorization.jl index 740a134aa5303..d200eff2f17bf 100644 --- a/stdlib/LinearAlgebra/test/factorization.jl +++ b/stdlib/LinearAlgebra/test/factorization.jl @@ -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, @@ -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,