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

restructure the cholesky test

Author: fredrikekre

NEWS.md

@@ -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
 ==========================

base/deprecated.jl

@@ -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

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

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