diff --git a/NEWS.md b/NEWS.md index 5c415e6c703ee0..5941161625b072 100644 --- a/NEWS.md +++ b/NEWS.md @@ -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 ========================== diff --git a/base/deprecated.jl b/base/deprecated.jl index 66d4af48e4cf82..2a8a9de828fc56 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -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 diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl index 2b5c89cad03efc..1e9fa98ea28a8c 100644 --- a/base/linalg/cholesky.jl +++ b/base/linalg/cholesky.jl @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 0364b4bfb9fa27..531f19ee15230f 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -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) @@ -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) @@ -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 @@ -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 - 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] @@ -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] @@ -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) @@ -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 @@ -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] @@ -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)))) @@ -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