From 3f77e34824e261a88f4bbce0e3c5047ba9696be0 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 7 Jul 2023 16:52:36 +0100 Subject: [PATCH 01/10] Implement reverse cholesky --- src/MatrixFactorizations.jl | 9 +- src/reversecholesky.jl | 252 ++++++++++++++++++++++++ test/runtests.jl | 1 + test/test_reversecholesky.jl | 365 +++++++++++++++++++++++++++++++++++ 4 files changed, 624 insertions(+), 3 deletions(-) create mode 100644 src/reversecholesky.jl create mode 100644 test/test_reversecholesky.jl diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl index b017912..4efb96c 100644 --- a/src/MatrixFactorizations.jl +++ b/src/MatrixFactorizations.jl @@ -2,13 +2,15 @@ module MatrixFactorizations using Base, LinearAlgebra, ArrayLayouts import Base: axes, axes1, getproperty, iterate, tail import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!, - copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym + copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym, + det, logdet, logabsdet, isposdef import LinearAlgebra.LAPACK: chkuplo, chktrans import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen, qr, axpy!, ldiv!, rdiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular, inv, chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, det, logabsdet, AbstractQ, _zeros, _cut_B, _ret_size, require_one_based_indexing, checksquare, - checknonsingular, ipiv2perm, copytri!, issuccess + checknonsingular, ipiv2perm, copytri!, issuccess, RealHermSymComplexHerm, + cholcopy, checkpositivedefinite, char_uplo import Base: getindex, setindex!, *, +, -, ==, <, <=, >, >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds @@ -28,7 +30,7 @@ import ArrayLayouts: reflector!, reflectorApply!, materialize!, @_layoutlmul, @_ -export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, choleskyinv!, choleskyinv +export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, choleskyinv!, choleskyinv, reversecholesky, reversecholesky!, ReverseCholesky const AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint const AbstractQtype = AbstractQ <: AbstractMatrix ? AbstractMatrix : AbstractQ @@ -138,5 +140,6 @@ include("ql.jl") include("rq.jl") include("choleskyinv.jl") include("polar.jl") +include("reversecholesky.jl") end #module diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl new file mode 100644 index 0000000..284825d --- /dev/null +++ b/src/reversecholesky.jl @@ -0,0 +1,252 @@ +# This file was modified from +# a part of Julia. License is MIT: https://julialang.org/license + +########################## +# Reverse Cholesky Factorization # +########################## + +""" + ReverseCholesky <: Factorization + +Matrix factorization type of the reverse Cholesky factorization of a dense symmetric/Hermitian +positive definite matrix `A`. This is the return type of [`reversecholesky`](@ref), +the corresponding matrix factorization function. + +The triangular reverse Cholesky factor can be obtained from the factorization `F::ReverseCholesky` +via `F.L` and `F.U`, where `A ≈ F.U * F.U' ≈ F.L' * F.L`. +``` +""" +struct ReverseCholesky{T,S<:AbstractMatrix} <: Factorization{T} + factors::S + uplo::Char + info::BlasInt + + function ReverseCholesky{T,S}(factors, uplo, info) where {T,S<:AbstractMatrix} + require_one_based_indexing(factors) + new(factors, uplo, info) + end +end +ReverseCholesky(A::AbstractMatrix{T}, uplo::Symbol, info::Integer) where {T} = + ReverseCholesky{T,typeof(A)}(A, char_uplo(uplo), info) +ReverseCholesky(A::AbstractMatrix{T}, uplo::AbstractChar, info::Integer) where {T} = + ReverseCholesky{T,typeof(A)}(A, uplo, info) +ReverseCholesky(U::UpperTriangular{T}) where {T} = ReverseCholesky{T,typeof(U.data)}(U.data, 'U', 0) +ReverseCholesky(L::LowerTriangular{T}) where {T} = ReverseCholesky{T,typeof(L.data)}(L.data, 'L', 0) + +# iteration for destructuring into components +Base.iterate(C::ReverseCholesky) = (C.L, Val(:U)) +Base.iterate(C::ReverseCholesky, ::Val{:U}) = (C.U, Val(:done)) +Base.iterate(C::ReverseCholesky, ::Val{:done}) = nothing + +## Non BLAS/LAPACK element types (generic) +function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) + require_one_based_indexing(A) + n = checksquare(A) + realdiag = eltype(A) <: Complex + @inbounds begin + for k = n:-1:1 + Akk = realdiag ? real(A[k,k]) : A[k,k] + for j = k+1:n + Akk -= realdiag ? abs2(A[k,j]) : A[k,j]'A[k,j] + end + A[k,k] = Akk + Akk, info = _reverse_chol!(Akk, UpperTriangular) + if info != 0 + return UpperTriangular(A), info + end + A[k,k] = Akk + AkkInv = inv(Akk) + for j = k+1:n + @simd for i = 1:k-1 + A[i,k] -= A[i,j]*A[k,j]' + end + end + for i = 1:k-1 + A[i,k] *= AkkInv' + end + end + end + return UpperTriangular(A), convert(BlasInt, 0) +end + +function _reverse_chol!(A::AbstractMatrix, ::Type{LowerTriangular}) + require_one_based_indexing(A) + n = checksquare(A) + realdiag = eltype(A) <: Complex + @inbounds begin + for k = n:-1:1 + Akk = realdiag ? real(A[k,k]) : A[k,k] + for i = k+1:n + Akk -= realdiag ? abs2(A[i,k]) : A[i,k]'A[i,k] + end + A[k,k] = Akk + Akk, info = _reverse_chol!(Akk, LowerTriangular) + if info != 0 + return LowerTriangular(A), info + end + A[k,k] = Akk + AkkInv = inv(copy(Akk')) + for j = 1:k-1 + for i = k+1:n + A[k,j] -= A[i,k]*A[i,j]' + end + A[k,j] = AkkInv*A[k,j] + end + end + return UpperTriangular(A), convert(BlasInt, 0) +end + +## Numbers +_reverse_chol!(x::Number, uplo) = LinearAlgebra._chol!(x, uplo) + +## for StridedMatrices, check that matrix is symmetric/Hermitian + +# cholesky!. Destructive methods for computing Cholesky factorization of real symmetric +# or Hermitian matrix +## No pivoting (default) +function reversecholesky!(A::RealHermSymComplexHerm, ::NoPivot = NoPivot(); check::Bool = true) + C, info = _reverse_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular) + check && checkpositivedefinite(info) + return ReverseCholesky(C.data, A.uplo, info) +end + +### for AbstractMatrix, check that matrix is symmetric/Hermitian +""" + reversecholesky!(A::AbstractMatrix, NoPivot(); check = true) -> ReverseCholesky + +The same as [`reversecholesky`](@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 reversecholesky!(A::AbstractMatrix, ::NoPivot = NoPivot(); check::Bool = true) + checksquare(A) + if !ishermitian(A) # return with info = -1 if not Hermitian + check && checkpositivedefinite(-1) + return ReverseCholesky(A, 'U', convert(BlasInt, -1)) + else + return reversecholesky!(Hermitian(A), NoPivot(); check = check) + end +end +@deprecate reversecholesky!(A::StridedMatrix, ::Val{false}; check::Bool = true) reversecholesky!(A, NoPivot(); check) false +@deprecate reversecholesky!(A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true) reversecholesky!(A, NoPivot(); check) false + + +# reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric +# or Hermitian matrix +## No pivoting (default) +""" + reversecholesky(A, NoPivot(); check = true) -> ReverseCholesky + +Compute the ReverseCholesky factorization of a dense symmetric positive definite matrix `A` +and return a [`ReverseCholesky`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref) or [`Hermitian`](@ref) +[`AbstractMatrix`](@ref) or a *perfectly* symmetric or Hermitian `AbstractMatrix`. + +The triangular ReverseCholesky factor can be obtained from the factorization `F` via `F.L` and `F.U`, +where `A ≈ F.U * F.U' ≈ F.L' * F.L`. +""" +reversecholesky(A::AbstractMatrix, ::NoPivot=NoPivot(); check::Bool = true) = + reversecholesky!(cholcopy(A); check) + +function reversecholesky(A::AbstractMatrix{Float16}, ::NoPivot=NoPivot(); check::Bool = true) + X = reversecholesky!(cholcopy(A); check = check) + return ReverseCholesky{Float16}(X) +end + + +## Number +function reversecholesky(x::Number, uplo::Symbol=:U) + C, info = _reverse_chol!(x, uplo) + xf = fill(C, 1, 1) + ReverseCholesky(xf, uplo, info) +end + + +function ReverseCholesky{T}(C::ReverseCholesky) where T + Cnew = convert(AbstractMatrix{T}, C.factors) + ReverseCholesky{T, typeof(Cnew)}(Cnew, C.uplo, C.info) +end +Factorization{T}(C::ReverseCholesky{T}) where {T} = C +Factorization{T}(C::ReverseCholesky) where {T} = ReverseCholesky{T}(C) + +AbstractMatrix(C::ReverseCholesky) = C.uplo == 'U' ? C.U*C.U' : C.L'*C.L +AbstractArray(C::ReverseCholesky) = AbstractMatrix(C) +Matrix(C::ReverseCholesky) = Array(AbstractArray(C)) +Array(C::ReverseCholesky) = Matrix(C) + +copy(C::ReverseCholesky) = ReverseCholesky(copy(C.factors), C.uplo, C.info) + + +size(C::ReverseCholesky) = size(C.factors) +size(C::ReverseCholesky, d::Integer) = size(C.factors, d) + +function getproperty(C::ReverseCholesky, d::Symbol) + Cfactors = getfield(C, :factors) + Cuplo = getfield(C, :uplo) + if d === :U + return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors')) + elseif d === :L + return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors')) + elseif d === :UL + return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors)) + else + return getfield(C, d) + end +end +Base.propertynames(F::ReverseCholesky, private::Bool=false) = + (:U, :L, :UL, (private ? fieldnames(typeof(F)) : ())...) + + + +issuccess(C::ReverseCholesky) = C.info == 0 + +adjoint(C::ReverseCholesky) = C + +function show(io::IO, mime::MIME{Symbol("text/plain")}, C::ReverseCholesky) + if issuccess(C) + summary(io, C); println(io) + println(io, "$(C.uplo) factor:") + show(io, mime, C.UL) + else + print(io, "Failed factorization of type $(typeof(C))") + end +end + +function ldiv!(C::ReverseCholesky, B::AbstractVecOrMat) + if C.uplo == 'L' + return ldiv!(LowerTriangular(C.factors), ldiv!(adjoint(LowerTriangular(C.factors)), B)) + else + return ldiv!(adjoint(UpperTriangular(C.factors)), ldiv!(UpperTriangular(C.factors), B)) + end +end + +function rdiv!(B::AbstractMatrix, C::ReverseCholesky) + if C.uplo == 'L' + return rdiv!(rdiv!(B, LowerTriangular(C.factors)), adjoint(LowerTriangular(C.factors))) + else + return rdiv!(rdiv!(B, adjoint(UpperTriangular(C.factors))), UpperTriangular(C.factors)) + end +end + +isposdef(C::ReverseCholesky) = C.info == 0 + +function det(C::ReverseCholesky) + dd = one(real(eltype(C))) + @inbounds for i in 1:size(C.factors,1) + dd *= real(C.factors[i,i])^2 + end + return dd +end + +function logdet(C::ReverseCholesky) + dd = zero(real(eltype(C))) + @inbounds for i in 1:size(C.factors,1) + dd += log(real(C.factors[i,i])) + end + dd + dd # instead of 2.0dd which can change the type +end + + + +logabsdet(C::ReverseCholesky) = logdet(C), one(eltype(C)) # since C is p.s.d. diff --git a/test/runtests.jl b/test/runtests.jl index b092d32..176c10b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -469,3 +469,4 @@ end include("test_polar.jl") +include("test_reversecholesky.jl") \ No newline at end of file diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl new file mode 100644 index 0000000..c594519 --- /dev/null +++ b/test/test_reversecholesky.jl @@ -0,0 +1,365 @@ +# This file is modifed from a part of Julia. License is MIT: https://julialang.org/license + +using Test, LinearAlgebra, Random, MatrixFactorizations +using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, + PosDefException, RankDeficientException, chkfullrank + +function unary_ops_tests(a, ca, tol; n=size(a, 1)) + @test inv(ca)*a ≈ Matrix(I, n, n) + @test a*inv(ca) ≈ Matrix(I, n, n) + @test abs((det(ca) - det(a))/det(ca)) <= tol # Ad hoc, but statistically verified, revisit + @test logdet(ca) ≈ logdet(a) + @test logdet(ca) ≈ log(det(ca)) # logdet is less likely to overflow + logabsdet_ca = logabsdet(ca) + logabsdet_a = logabsdet(a) + @test logabsdet_ca[1] ≈ logabsdet_a[1] + @test logabsdet_ca[2] ≈ logabsdet_a[2] + @test isposdef(ca) + @test_throws ErrorException ca.Z + @test size(ca) == size(a) + @test Array(copy(ca)) ≈ a +end + +function factor_recreation_tests(a_U, a_L) + c_U = reversecholesky(a_U) + c_L = reversecholesky(a_L) + cl = c_L.U + ls = c_L.L + @test Array(c_U) ≈ Array(c_L) ≈ a_U + @test ls'*ls ≈ a_U + @test triu(c_U.factors) ≈ c_U.U + @test tril(c_L.factors) ≈ c_L.L + @test istriu(cl) + @test cl*cl' ≈ a_U + @test cl*cl' ≈ a_L +end + +@testset "core functionality" begin + n = 10 + + + + # Split n into 2 parts for tests needing two matrices + n1 = div(n, 2) + n2 = 2*n1 + + Random.seed!(12344) + + areal = randn(n,n)/2 + aimg = randn(n,n)/2 + a2real = randn(n,n)/2 + a2img = randn(n,n)/2 + breal = randn(n,2)/2 + bimg = randn(n,2)/2 + + @testset "basic sanity check" begin + A = [4 2; 2 4] + U = reversecholesky(A).U + @test U*U' ≈ A + @test cholesky(A).L[end:-1:1,end:-1:1] ≈ U + + A = Symmetric(areal) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + end + + for eltya in (Float32, Float64, ComplexF32, ComplexF64, 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) + + ε = εa = eps(abs(float(one(eltya)))) + + # Test of symmetric pos. def. strided matrix + apd = a'*a + @inferred reversecholesky(apd) + capd = reversecholesky(apd) + r = capd.U + κ = cond(apd, 1) #condition number + + unary_ops_tests(apd, capd, ε*κ*n) + if eltya != Int + @test Factorization{eltya}(capd) === capd + if eltya <: Real + @test Array(Factorization{complex(eltya)}(capd)) ≈ Array(factorize(complex(apd))) + @test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya) + end + end + @testset "throw for non-square input" begin + A = rand(eltya, 2, 3) + @test_throws DimensionMismatch reversecholesky(A) + @test_throws DimensionMismatch reversecholesky!(A) + end + + #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 + + #these tests were failing on 64-bit linux when inside the inner loop + #for eltya = ComplexF32 and eltyb = Int. The E[i,j] had NaN32 elements + #but only with Random.seed!(1234321) set before the loops. + E = abs.(apd - r*r') + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + E = abs.(apd - Matrix(capd)) + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + @test LinearAlgebra.issuccess(capd) + @inferred(logdet(capd)) + + apos = apd[1,1] + @test all(x -> x ≈ √apos, cholesky(apos).factors) + + # Test cholesky with Symmetric/Hermitian upper/lower + apds = Symmetric(apd) + apdsL = Symmetric(apd, :L) + apdh = Hermitian(apd) + apdhL = Hermitian(apd, :L) + if eltya <: Real + capds = reversecholesky(apds) + unary_ops_tests(apds, capds, ε*κ*n) + if eltya <: BlasReal + capds = cholesky!(copy(apds)) + unary_ops_tests(apds, capds, ε*κ*n) + end + ulstring = sprint((t, s) -> show(t, "text/plain", s), capds.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capds) == "$(typeof(capds))\nU factor:\n$ulstring" + else + capdh = reversecholesky(apdh) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apdh)) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apd)) + unary_ops_tests(apd, capdh, ε*κ*n) + ulstring = sprint((t, s) -> show(t, "text/plain", s), capdh.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capdh) == "$(typeof(capdh))\nU factor:\n$ulstring" + end + + # test reversecholesky of 2x2 Strang matrix + S = SymTridiagonal{eltya}([2, 2], [-1]) + for uplo in (:U, :L) + @test Matrix(@inferred reversecholesky(Hermitian(S, uplo))) ≈ S + if eltya <: Real + @test Matrix(@inferred reversecholesky(Symmetric(S, uplo))) ≈ S + end + end + @test Matrix(reversecholesky(S).U) ≈ [2 -1; 0 sqrt(eltya(3))] / sqrt(eltya(2)) + @test Matrix(reversecholesky(S)) ≈ S + + # test extraction of factor and re-creating original matrix + if eltya <: Real + factor_recreation_tests(apds, apdsL) + else + factor_recreation_tests(apdh, apdhL) + end + + + + for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + for b in (b, view(b, 1:n, 1)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + x = capd\b + @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + + @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit + + if eltya != BigFloat && eltyb != BigFloat + lapd = reversecholesky(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 + end + end + + for eltyb in (Float64, ComplexF64) + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = (eltya <: Complex || eltyb <: Complex) ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + BB = copy(B) + ldiv!(capd, BB) + @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)*ε)*κ + end + end + + @testset "solve with generic Cholesky" begin + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = eltya <: Complex ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltype(B))))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + cpapd = reversecholesky(eltya <: Complex ? apdh : apds) + 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 = reversecholesky(eltya <: Complex ? apdhL : apdsL) + 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)*ε)*κ + end + end + end +end + +@testset "behavior for non-positive definite matrices" for T in (Float64, ComplexF64) + A = T[1 2; 2 1] + B = T[1 2; 0 1] + # check = (true|false) + for M in (A, Hermitian(A), B) + @test_throws PosDefException reversecholesky(M) + @test_throws PosDefException reversecholesky!(copy(M)) + @test_throws PosDefException reversecholesky(M; check = true) + @test_throws PosDefException reversecholesky!(copy(M); check = true) + @test !LinearAlgebra.issuccess(reversecholesky(M; check = false)) + @test !LinearAlgebra.issuccess(reversecholesky!(copy(M); check = false)) + end + str = sprint((io, x) -> show(io, "text/plain", x), reversecholesky(A; check = false)) +end + +@testset "Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices" begin + X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3] + L = Matrix(MatrixFactorizations._reverse_chol!(X*X', LowerTriangular)[1]) + U = Matrix(MatrixFactorizations._reverse_chol!(X*X', UpperTriangular)[1]) + XX = Matrix(X*X') + + @test sum(sum(norm, L*L' - XX)) < eps() + @test sum(sum(norm, U'*U - XX)) < eps() +end + +@testset "Non-strided Cholesky solves" begin + B = randn(5, 5) + v = rand(5) + @test reversecholesky(Diagonal(v)) \ B ≈ Diagonal(v) \ B + @test B / reversecholesky(Diagonal(v)) ≈ B / Diagonal(v) + @test inv(reversecholesky(Diagonal(v)))::Diagonal ≈ Diagonal(1 ./ v) +end + + +@testset "cholesky Diagonal" begin + # real + d = abs.(randn(3)) .+ 0.1 + D = Diagonal(d) + CD = reversecholesky(D) + CM = reversecholesky(Matrix(D)) + @test CD isa Cholesky{Float64} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + F = reversecholesky(Hermitian(I(3))) + @test F isa Cholesky{Float64,<:Diagonal} + @test Matrix(F) ≈ I(3) + + # real, failing + @test_throws PosDefException reversecholesky(Diagonal([1.0, -2.0])) + Dnpd = reversecholesky(Diagonal([1.0, -2.0]); check = false) + @test Dnpd.info == 2 + + # complex + D = complex(D) + CD = reversecholesky(Hermitian(D)) + CM = reversecholesky(Matrix(Hermitian(D))) + @test CD isa Cholesky{ComplexF64,<:Diagonal} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + # complex, failing + D[2, 2] = 0.0 + 0im + @test_throws PosDefException reversecholesky(D) + Dnpd = reversecholesky(D; check = false) + @test Dnpd.info == 2 + + # InexactError for Int + @test_throws InexactError reversecholesky!(Diagonal([2, 1])) +end + +@testset "Cholesky for AbstractMatrix" begin + S = SymTridiagonal(fill(2.0, 4), ones(3)) + C = reversecholesky(S) + @test C.L * C.U ≈ S +end + +@testset "constructor with non-BlasInt arguments" begin + + x = rand(5,5) + chol = reversecholesky(x'x) + + factors, uplo, info = chol.factors, chol.uplo, chol.info + + @test ReverseCholesky(factors, uplo, Int32(info)) == chol + @test ReverseCholesky(factors, uplo, Int64(info)) == chol +end + + + +@testset "issue #37356, diagonal elements of hermitian generic matrix" begin + B = Hermitian(hcat([one(BigFloat) + im])) + @test Matrix(reversecholesky(B)) ≈ B + C = Hermitian(hcat([one(BigFloat) + im]), :L) + @test Matrix(reversecholesky(C)) ≈ C +end + +@testset "constructing a ReverseCholesky factor from a triangular matrix" begin + A = [1.0 2.0; 3.0 4.0] + let + U = UpperTriangular(A) + C = ReverseCholesky(U) + @test C isa ReverseCholesky{Float64} + @test C.U == U + @test C.L == U' + end + let + L = LowerTriangular(A) + C = ReverseCholesky(L) + @test C isa ReverseCholesky{Float64} + @test C.L == L + @test C.U == L' + end +end + +@testset "adjoint of ReverseCholesky" begin + A = randn(5, 5) + A = A'A + F = reversecholesky(A) + b = ones(size(A, 1)) + @test F\b == F'\b +end + +@testset "Float16" begin + A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.] + B = reversecholesky(A) + B32 = reversecholesky(Float32.(A)) + @test B isa ReverseCholesky{Float16, Matrix{Float16}} + @test B.U isa UpperTriangular{Float16, Matrix{Float16}} + @test B.L isa LowerTriangular{Float16, Matrix{Float16}} + @test B.UL isa UpperTriangular{Float16, Matrix{Float16}} + @test B.U ≈ B32.U + @test B.L ≈ B32.L + @test B.UL ≈ B32.UL + @test Matrix(B) ≈ A +end From b65c2809d54bc158e33a7faa53c0d99b46e9057a Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 8 Jul 2023 20:09:59 +0100 Subject: [PATCH 02/10] tests pass --- src/reversecholesky.jl | 18 ++++++++++++++---- test/test_reversecholesky.jl | 31 +++++++++++++++++-------------- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index 284825d..9ed730e 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -88,10 +88,11 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{LowerTriangular}) AkkInv = inv(copy(Akk')) for j = 1:k-1 for i = k+1:n - A[k,j] -= A[i,k]*A[i,j]' + A[k,j] -= A[i,k]'*A[i,j] end A[k,j] = AkkInv*A[k,j] end + end end return UpperTriangular(A), convert(BlasInt, 0) end @@ -129,9 +130,6 @@ function reversecholesky!(A::AbstractMatrix, ::NoPivot = NoPivot(); check::Bool return reversecholesky!(Hermitian(A), NoPivot(); check = check) end end -@deprecate reversecholesky!(A::StridedMatrix, ::Val{false}; check::Bool = true) reversecholesky!(A, NoPivot(); check) false -@deprecate reversecholesky!(A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true) reversecholesky!(A, NoPivot(); check) false - # reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric # or Hermitian matrix @@ -250,3 +248,15 @@ end logabsdet(C::ReverseCholesky) = logdet(C), one(eltype(C)) # since C is p.s.d. + + +function getproperty(C::ReverseCholesky{<:Any,<:Diagonal}, d::Symbol) + Cfactors = getfield(C, :factors) + if d in (:U, :L, :UL) + return Cfactors + else + return getfield(C, d) + end +end + +inv(C::ReverseCholesky{<:Any,<:Diagonal}) = Diagonal(map(inv∘abs2, C.factors.diag)) \ No newline at end of file diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index c594519..e732798 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -37,8 +37,6 @@ end @testset "core functionality" begin n = 10 - - # Split n into 2 parts for tests needing two matrices n1 = div(n, 2) n2 = 2*n1 @@ -65,6 +63,14 @@ end A = Hermitian(areal + im*aimg) + 10I U = reversecholesky(A).U @test U*U' ≈ A + + A = Symmetric(areal, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A end for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) @@ -146,7 +152,7 @@ end @test Matrix(@inferred reversecholesky(Symmetric(S, uplo))) ≈ S end end - @test Matrix(reversecholesky(S).U) ≈ [2 -1; 0 sqrt(eltya(3))] / sqrt(eltya(2)) + @test Matrix(reversecholesky(S).U) ≈ [sqrt(eltya(3)) -1; 0 2] / sqrt(eltya(2)) @test Matrix(reversecholesky(S)) ≈ S # test extraction of factor and re-creating original matrix @@ -246,8 +252,8 @@ end U = Matrix(MatrixFactorizations._reverse_chol!(X*X', UpperTriangular)[1]) XX = Matrix(X*X') - @test sum(sum(norm, L*L' - XX)) < eps() - @test sum(sum(norm, U'*U - XX)) < eps() + @test_broken sum(sum(norm, L'*L - XX)) < eps() + @test_broken sum(sum(norm, U*U' - XX)) < eps() end @testset "Non-strided Cholesky solves" begin @@ -265,34 +271,32 @@ end D = Diagonal(d) CD = reversecholesky(D) CM = reversecholesky(Matrix(D)) - @test CD isa Cholesky{Float64} + @test CD isa ReverseCholesky{Float64} @test CD.U ≈ Diagonal(.√d) ≈ CM.U @test D ≈ CD.L * CD.U @test CD.info == 0 F = reversecholesky(Hermitian(I(3))) - @test F isa Cholesky{Float64,<:Diagonal} + @test F isa ReverseCholesky{Float64,<:Diagonal} @test Matrix(F) ≈ I(3) # real, failing @test_throws PosDefException reversecholesky(Diagonal([1.0, -2.0])) Dnpd = reversecholesky(Diagonal([1.0, -2.0]); check = false) - @test Dnpd.info == 2 + @test Dnpd.info == 1 # complex D = complex(D) CD = reversecholesky(Hermitian(D)) CM = reversecholesky(Matrix(Hermitian(D))) - @test CD isa Cholesky{ComplexF64,<:Diagonal} + @test CD isa ReverseCholesky{ComplexF64,<:Diagonal} @test CD.U ≈ Diagonal(.√d) ≈ CM.U @test D ≈ CD.L * CD.U @test CD.info == 0 # complex, failing D[2, 2] = 0.0 + 0im - @test_throws PosDefException reversecholesky(D) - Dnpd = reversecholesky(D; check = false) - @test Dnpd.info == 2 + @test_throws ArgumentError reversecholesky(D) # InexactError for Int @test_throws InexactError reversecholesky!(Diagonal([2, 1])) @@ -301,11 +305,10 @@ end @testset "Cholesky for AbstractMatrix" begin S = SymTridiagonal(fill(2.0, 4), ones(3)) C = reversecholesky(S) - @test C.L * C.U ≈ S + @test C.U * C.L ≈ S end @testset "constructor with non-BlasInt arguments" begin - x = rand(5,5) chol = reversecholesky(x'x) From 216b408824b9d40652c85c073c7de0f26e05142c Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 8 Jul 2023 20:16:25 +0100 Subject: [PATCH 03/10] Drop old versions --- .github/workflows/ci.yml | 4 +- Project.toml | 6 +- test/test_reversecholesky.jl | 578 ++++++++++++++++++----------------- 3 files changed, 294 insertions(+), 294 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3239dba..ed5fbd5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,9 +20,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' - - '1' - - '^1.9.0-0' + - '1.9' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index 7ea83cf..d04beb9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MatrixFactorizations" uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" -version = "1.0" +version = "2.0" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" @@ -9,8 +9,8 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -ArrayLayouts = "0.8.11, 1" -julia = "1.6" +ArrayLayouts = "1.0" +julia = "1.9" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index e732798..2ce9b98 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -34,335 +34,337 @@ function factor_recreation_tests(a_U, a_L) @test cl*cl' ≈ a_L end -@testset "core functionality" begin - n = 10 - - # Split n into 2 parts for tests needing two matrices - n1 = div(n, 2) - n2 = 2*n1 - - Random.seed!(12344) - - areal = randn(n,n)/2 - aimg = randn(n,n)/2 - a2real = randn(n,n)/2 - a2img = randn(n,n)/2 - breal = randn(n,2)/2 - bimg = randn(n,2)/2 - - @testset "basic sanity check" begin - A = [4 2; 2 4] - U = reversecholesky(A).U - @test U*U' ≈ A - @test cholesky(A).L[end:-1:1,end:-1:1] ≈ U - - A = Symmetric(areal) + 10I - U = reversecholesky(A).U - @test U*U' ≈ A - - A = Hermitian(areal + im*aimg) + 10I - U = reversecholesky(A).U - @test U*U' ≈ A - - A = Symmetric(areal, :L) + 10I - U = reversecholesky(A).U - @test U*U' ≈ A - - A = Hermitian(areal + im*aimg, :L) + 10I - U = reversecholesky(A).U - @test U*U' ≈ A - end - - for eltya in (Float32, Float64, ComplexF32, ComplexF64, 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) +@testset "ReverseCholesky" begin + @testset "core functionality" begin + n = 10 + + # Split n into 2 parts for tests needing two matrices + n1 = div(n, 2) + n2 = 2*n1 + + Random.seed!(12344) + + areal = randn(n,n)/2 + aimg = randn(n,n)/2 + a2real = randn(n,n)/2 + a2img = randn(n,n)/2 + breal = randn(n,2)/2 + bimg = randn(n,2)/2 + + @testset "basic sanity check" begin + A = [4 2; 2 4] + U = reversecholesky(A).U + @test U*U' ≈ A + @test cholesky(A).L[end:-1:1,end:-1:1] ≈ U + + A = Symmetric(areal) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Symmetric(areal, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + end - ε = εa = eps(abs(float(one(eltya)))) + for eltya in (Float32, Float64, ComplexF32, ComplexF64, 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) + + ε = εa = eps(abs(float(one(eltya)))) + + # Test of symmetric pos. def. strided matrix + apd = a'*a + @inferred reversecholesky(apd) + capd = reversecholesky(apd) + r = capd.U + κ = cond(apd, 1) #condition number + + unary_ops_tests(apd, capd, ε*κ*n) + if eltya != Int + @test Factorization{eltya}(capd) === capd + if eltya <: Real + @test Array(Factorization{complex(eltya)}(capd)) ≈ Array(factorize(complex(apd))) + @test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya) + end + end + @testset "throw for non-square input" begin + A = rand(eltya, 2, 3) + @test_throws DimensionMismatch reversecholesky(A) + @test_throws DimensionMismatch reversecholesky!(A) + end - # Test of symmetric pos. def. strided matrix - apd = a'*a - @inferred reversecholesky(apd) - capd = reversecholesky(apd) - r = capd.U - κ = cond(apd, 1) #condition number + #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 - unary_ops_tests(apd, capd, ε*κ*n) - if eltya != Int - @test Factorization{eltya}(capd) === capd - if eltya <: Real - @test Array(Factorization{complex(eltya)}(capd)) ≈ Array(factorize(complex(apd))) - @test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya) + #these tests were failing on 64-bit linux when inside the inner loop + #for eltya = ComplexF32 and eltyb = Int. The E[i,j] had NaN32 elements + #but only with Random.seed!(1234321) set before the loops. + E = abs.(apd - r*r') + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) end - end - @testset "throw for non-square input" begin - A = rand(eltya, 2, 3) - @test_throws DimensionMismatch reversecholesky(A) - @test_throws DimensionMismatch reversecholesky!(A) - end + E = abs.(apd - Matrix(capd)) + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + @test LinearAlgebra.issuccess(capd) + @inferred(logdet(capd)) - #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 + apos = apd[1,1] + @test all(x -> x ≈ √apos, cholesky(apos).factors) - #these tests were failing on 64-bit linux when inside the inner loop - #for eltya = ComplexF32 and eltyb = Int. The E[i,j] had NaN32 elements - #but only with Random.seed!(1234321) set before the loops. - E = abs.(apd - r*r') - for i=1:n, j=1:n - @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) - end - E = abs.(apd - Matrix(capd)) - for i=1:n, j=1:n - @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) - end - @test LinearAlgebra.issuccess(capd) - @inferred(logdet(capd)) - - apos = apd[1,1] - @test all(x -> x ≈ √apos, cholesky(apos).factors) - - # Test cholesky with Symmetric/Hermitian upper/lower - apds = Symmetric(apd) - apdsL = Symmetric(apd, :L) - apdh = Hermitian(apd) - apdhL = Hermitian(apd, :L) - if eltya <: Real - capds = reversecholesky(apds) - unary_ops_tests(apds, capds, ε*κ*n) - if eltya <: BlasReal - capds = cholesky!(copy(apds)) + # Test cholesky with Symmetric/Hermitian upper/lower + apds = Symmetric(apd) + apdsL = Symmetric(apd, :L) + apdh = Hermitian(apd) + apdhL = Hermitian(apd, :L) + if eltya <: Real + capds = reversecholesky(apds) unary_ops_tests(apds, capds, ε*κ*n) + if eltya <: BlasReal + capds = cholesky!(copy(apds)) + unary_ops_tests(apds, capds, ε*κ*n) + end + ulstring = sprint((t, s) -> show(t, "text/plain", s), capds.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capds) == "$(typeof(capds))\nU factor:\n$ulstring" + else + capdh = reversecholesky(apdh) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apdh)) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apd)) + unary_ops_tests(apd, capdh, ε*κ*n) + ulstring = sprint((t, s) -> show(t, "text/plain", s), capdh.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capdh) == "$(typeof(capdh))\nU factor:\n$ulstring" end - ulstring = sprint((t, s) -> show(t, "text/plain", s), capds.UL) - @test sprint((t, s) -> show(t, "text/plain", s), capds) == "$(typeof(capds))\nU factor:\n$ulstring" - else - capdh = reversecholesky(apdh) - unary_ops_tests(apdh, capdh, ε*κ*n) - capdh = reversecholesky!(copy(apdh)) - unary_ops_tests(apdh, capdh, ε*κ*n) - capdh = reversecholesky!(copy(apd)) - unary_ops_tests(apd, capdh, ε*κ*n) - ulstring = sprint((t, s) -> show(t, "text/plain", s), capdh.UL) - @test sprint((t, s) -> show(t, "text/plain", s), capdh) == "$(typeof(capdh))\nU factor:\n$ulstring" - end - # test reversecholesky of 2x2 Strang matrix - S = SymTridiagonal{eltya}([2, 2], [-1]) - for uplo in (:U, :L) - @test Matrix(@inferred reversecholesky(Hermitian(S, uplo))) ≈ S + # test reversecholesky of 2x2 Strang matrix + S = SymTridiagonal{eltya}([2, 2], [-1]) + for uplo in (:U, :L) + @test Matrix(@inferred reversecholesky(Hermitian(S, uplo))) ≈ S + if eltya <: Real + @test Matrix(@inferred reversecholesky(Symmetric(S, uplo))) ≈ S + end + end + @test Matrix(reversecholesky(S).U) ≈ [sqrt(eltya(3)) -1; 0 2] / sqrt(eltya(2)) + @test Matrix(reversecholesky(S)) ≈ S + + # test extraction of factor and re-creating original matrix if eltya <: Real - @test Matrix(@inferred reversecholesky(Symmetric(S, uplo))) ≈ S + factor_recreation_tests(apds, apdsL) + else + factor_recreation_tests(apdh, apdhL) end - end - @test Matrix(reversecholesky(S).U) ≈ [sqrt(eltya(3)) -1; 0 2] / sqrt(eltya(2)) - @test Matrix(reversecholesky(S)) ≈ S - - # test extraction of factor and re-creating original matrix - if eltya <: Real - factor_recreation_tests(apds, apdsL) - else - factor_recreation_tests(apdh, apdhL) - end - for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) - b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) - ε = max(εa,εb) + for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) - for b in (b, view(b, 1:n, 1)) # Array and SubArray + for b in (b, view(b, 1:n, 1)) # Array and SubArray - # Test error bound on linear solver: LAWNS 14, Theorem 2.1 - # This is a surprisingly loose bound - x = capd\b - @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ - @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + x = capd\b + @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ - @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit + @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit - if eltya != BigFloat && eltyb != BigFloat - lapd = reversecholesky(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 + if eltya != BigFloat && eltyb != BigFloat + lapd = reversecholesky(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 end end - end - for eltyb in (Float64, ComplexF64) - Breal = convert(Matrix{BigFloat}, randn(n,n)/2) - Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) - B = (eltya <: Complex || eltyb <: Complex) ? complex.(Breal, Bimg) : Breal - εb = eps(abs(float(one(eltyb)))) - ε = max(εa,εb) - - for B in (B, view(B, 1:n, 1:n)) # Array and SubArray - - # Test error bound on linear solver: LAWNS 14, Theorem 2.1 - # This is a surprisingly loose bound - BB = copy(B) - ldiv!(capd, BB) - @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)*ε)*κ + for eltyb in (Float64, ComplexF64) + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = (eltya <: Complex || eltyb <: Complex) ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + BB = copy(B) + ldiv!(capd, BB) + @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)*ε)*κ + end end - end - @testset "solve with generic Cholesky" begin - Breal = convert(Matrix{BigFloat}, randn(n,n)/2) - Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) - B = eltya <: Complex ? complex.(Breal, Bimg) : Breal - εb = eps(abs(float(one(eltype(B))))) - ε = max(εa,εb) - - for B in (B, view(B, 1:n, 1:n)) # Array and SubArray - - # Test error bound on linear solver: LAWNS 14, Theorem 2.1 - # This is a surprisingly loose bound - cpapd = reversecholesky(eltya <: Complex ? apdh : apds) - 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 = reversecholesky(eltya <: Complex ? apdhL : apdsL) - 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)*ε)*κ + @testset "solve with generic Cholesky" begin + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = eltya <: Complex ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltype(B))))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + cpapd = reversecholesky(eltya <: Complex ? apdh : apds) + 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 = reversecholesky(eltya <: Complex ? apdhL : apdsL) + 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)*ε)*κ + end end end end -end -@testset "behavior for non-positive definite matrices" for T in (Float64, ComplexF64) - A = T[1 2; 2 1] - B = T[1 2; 0 1] - # check = (true|false) - for M in (A, Hermitian(A), B) - @test_throws PosDefException reversecholesky(M) - @test_throws PosDefException reversecholesky!(copy(M)) - @test_throws PosDefException reversecholesky(M; check = true) - @test_throws PosDefException reversecholesky!(copy(M); check = true) - @test !LinearAlgebra.issuccess(reversecholesky(M; check = false)) - @test !LinearAlgebra.issuccess(reversecholesky!(copy(M); check = false)) + @testset "behavior for non-positive definite matrices" for T in (Float64, ComplexF64) + A = T[1 2; 2 1] + B = T[1 2; 0 1] + # check = (true|false) + for M in (A, Hermitian(A), B) + @test_throws PosDefException reversecholesky(M) + @test_throws PosDefException reversecholesky!(copy(M)) + @test_throws PosDefException reversecholesky(M; check = true) + @test_throws PosDefException reversecholesky!(copy(M); check = true) + @test !LinearAlgebra.issuccess(reversecholesky(M; check = false)) + @test !LinearAlgebra.issuccess(reversecholesky!(copy(M); check = false)) + end + str = sprint((io, x) -> show(io, "text/plain", x), reversecholesky(A; check = false)) end - str = sprint((io, x) -> show(io, "text/plain", x), reversecholesky(A; check = false)) -end -@testset "Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices" begin - X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3] - L = Matrix(MatrixFactorizations._reverse_chol!(X*X', LowerTriangular)[1]) - U = Matrix(MatrixFactorizations._reverse_chol!(X*X', UpperTriangular)[1]) - XX = Matrix(X*X') + @testset "Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices" begin + X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3] + L = Matrix(MatrixFactorizations._reverse_chol!(X*X', LowerTriangular)[1]) + U = Matrix(MatrixFactorizations._reverse_chol!(X*X', UpperTriangular)[1]) + XX = Matrix(X*X') - @test_broken sum(sum(norm, L'*L - XX)) < eps() - @test_broken sum(sum(norm, U*U' - XX)) < eps() -end + @test_broken sum(sum(norm, L'*L - XX)) < eps() + @test_broken sum(sum(norm, U*U' - XX)) < eps() + end -@testset "Non-strided Cholesky solves" begin - B = randn(5, 5) - v = rand(5) - @test reversecholesky(Diagonal(v)) \ B ≈ Diagonal(v) \ B - @test B / reversecholesky(Diagonal(v)) ≈ B / Diagonal(v) - @test inv(reversecholesky(Diagonal(v)))::Diagonal ≈ Diagonal(1 ./ v) -end + @testset "Non-strided Cholesky solves" begin + B = randn(5, 5) + v = rand(5) + @test reversecholesky(Diagonal(v)) \ B ≈ Diagonal(v) \ B + @test B / reversecholesky(Diagonal(v)) ≈ B / Diagonal(v) + @test inv(reversecholesky(Diagonal(v)))::Diagonal ≈ Diagonal(1 ./ v) + end -@testset "cholesky Diagonal" begin - # real - d = abs.(randn(3)) .+ 0.1 - D = Diagonal(d) - CD = reversecholesky(D) - CM = reversecholesky(Matrix(D)) - @test CD isa ReverseCholesky{Float64} - @test CD.U ≈ Diagonal(.√d) ≈ CM.U - @test D ≈ CD.L * CD.U - @test CD.info == 0 - - F = reversecholesky(Hermitian(I(3))) - @test F isa ReverseCholesky{Float64,<:Diagonal} - @test Matrix(F) ≈ I(3) - - # real, failing - @test_throws PosDefException reversecholesky(Diagonal([1.0, -2.0])) - Dnpd = reversecholesky(Diagonal([1.0, -2.0]); check = false) - @test Dnpd.info == 1 - - # complex - D = complex(D) - CD = reversecholesky(Hermitian(D)) - CM = reversecholesky(Matrix(Hermitian(D))) - @test CD isa ReverseCholesky{ComplexF64,<:Diagonal} - @test CD.U ≈ Diagonal(.√d) ≈ CM.U - @test D ≈ CD.L * CD.U - @test CD.info == 0 - - # complex, failing - D[2, 2] = 0.0 + 0im - @test_throws ArgumentError reversecholesky(D) - - # InexactError for Int - @test_throws InexactError reversecholesky!(Diagonal([2, 1])) -end + @testset "cholesky Diagonal" begin + # real + d = abs.(randn(3)) .+ 0.1 + D = Diagonal(d) + CD = reversecholesky(D) + CM = reversecholesky(Matrix(D)) + @test CD isa ReverseCholesky{Float64} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + F = reversecholesky(Hermitian(I(3))) + @test F isa ReverseCholesky{Float64,<:Diagonal} + @test Matrix(F) ≈ I(3) + + # real, failing + @test_throws PosDefException reversecholesky(Diagonal([1.0, -2.0])) + Dnpd = reversecholesky(Diagonal([1.0, -2.0]); check = false) + @test Dnpd.info == 1 + + # complex + D = complex(D) + CD = reversecholesky(Hermitian(D)) + CM = reversecholesky(Matrix(Hermitian(D))) + @test CD isa ReverseCholesky{ComplexF64,<:Diagonal} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + # complex, failing + D[2, 2] = 0.0 + 0im + @test_throws ArgumentError reversecholesky(D) + + # InexactError for Int + @test_throws InexactError reversecholesky!(Diagonal([2, 1])) + end -@testset "Cholesky for AbstractMatrix" begin - S = SymTridiagonal(fill(2.0, 4), ones(3)) - C = reversecholesky(S) - @test C.U * C.L ≈ S -end + @testset "Cholesky for AbstractMatrix" begin + S = SymTridiagonal(fill(2.0, 4), ones(3)) + C = reversecholesky(S) + @test C.U * C.L ≈ S + end -@testset "constructor with non-BlasInt arguments" begin - x = rand(5,5) - chol = reversecholesky(x'x) + @testset "constructor with non-BlasInt arguments" begin + x = rand(5,5) + chol = reversecholesky(x'x) - factors, uplo, info = chol.factors, chol.uplo, chol.info - - @test ReverseCholesky(factors, uplo, Int32(info)) == chol - @test ReverseCholesky(factors, uplo, Int64(info)) == chol -end + factors, uplo, info = chol.factors, chol.uplo, chol.info + @test ReverseCholesky(factors, uplo, Int32(info)) == chol + @test ReverseCholesky(factors, uplo, Int64(info)) == chol + end -@testset "issue #37356, diagonal elements of hermitian generic matrix" begin - B = Hermitian(hcat([one(BigFloat) + im])) - @test Matrix(reversecholesky(B)) ≈ B - C = Hermitian(hcat([one(BigFloat) + im]), :L) - @test Matrix(reversecholesky(C)) ≈ C -end -@testset "constructing a ReverseCholesky factor from a triangular matrix" begin - A = [1.0 2.0; 3.0 4.0] - let - U = UpperTriangular(A) - C = ReverseCholesky(U) - @test C isa ReverseCholesky{Float64} - @test C.U == U - @test C.L == U' + @testset "issue #37356, diagonal elements of hermitian generic matrix" begin + B = Hermitian(hcat([one(BigFloat) + im])) + @test Matrix(reversecholesky(B)) ≈ B + C = Hermitian(hcat([one(BigFloat) + im]), :L) + @test Matrix(reversecholesky(C)) ≈ C end - let - L = LowerTriangular(A) - C = ReverseCholesky(L) - @test C isa ReverseCholesky{Float64} - @test C.L == L - @test C.U == L' + + @testset "constructing a ReverseCholesky factor from a triangular matrix" begin + A = [1.0 2.0; 3.0 4.0] + let + U = UpperTriangular(A) + C = ReverseCholesky(U) + @test C isa ReverseCholesky{Float64} + @test C.U == U + @test C.L == U' + end + let + L = LowerTriangular(A) + C = ReverseCholesky(L) + @test C isa ReverseCholesky{Float64} + @test C.L == L + @test C.U == L' + end end -end -@testset "adjoint of ReverseCholesky" begin - A = randn(5, 5) - A = A'A - F = reversecholesky(A) - b = ones(size(A, 1)) - @test F\b == F'\b -end + @testset "adjoint of ReverseCholesky" begin + A = randn(5, 5) + A = A'A + F = reversecholesky(A) + b = ones(size(A, 1)) + @test F\b == F'\b + end -@testset "Float16" begin - A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.] - B = reversecholesky(A) - B32 = reversecholesky(Float32.(A)) - @test B isa ReverseCholesky{Float16, Matrix{Float16}} - @test B.U isa UpperTriangular{Float16, Matrix{Float16}} - @test B.L isa LowerTriangular{Float16, Matrix{Float16}} - @test B.UL isa UpperTriangular{Float16, Matrix{Float16}} - @test B.U ≈ B32.U - @test B.L ≈ B32.L - @test B.UL ≈ B32.UL - @test Matrix(B) ≈ A -end + @testset "Float16" begin + A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.] + B = reversecholesky(A) + B32 = reversecholesky(Float32.(A)) + @test B isa ReverseCholesky{Float16, Matrix{Float16}} + @test B.U isa UpperTriangular{Float16, Matrix{Float16}} + @test B.L isa LowerTriangular{Float16, Matrix{Float16}} + @test B.UL isa UpperTriangular{Float16, Matrix{Float16}} + @test B.U ≈ B32.U + @test B.L ≈ B32.L + @test B.UL ≈ B32.UL + @test Matrix(B) ≈ A + end +end \ No newline at end of file From a4d13b2d9b91f57630b6b3e280f2184cec897e18 Mon Sep 17 00:00:00 2001 From: Timon Salar Gutleb Date: Sat, 8 Jul 2023 20:18:38 +0100 Subject: [PATCH 04/10] remove choleskyinv (#50) Co-authored-by: Sheehan Olver --- src/MatrixFactorizations.jl | 3 +-- test/runtests.jl | 29 ----------------------------- 2 files changed, 1 insertion(+), 31 deletions(-) diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl index 4efb96c..fa320be 100644 --- a/src/MatrixFactorizations.jl +++ b/src/MatrixFactorizations.jl @@ -29,8 +29,8 @@ import ArrayLayouts: reflector!, reflectorApply!, materialize!, @_layoutlmul, @_ layout_getindex +export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, reversecholesky, reversecholesky!, ReverseCholesky -export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, choleskyinv!, choleskyinv, reversecholesky, reversecholesky!, ReverseCholesky const AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint const AbstractQtype = AbstractQ <: AbstractMatrix ? AbstractMatrix : AbstractQ @@ -138,7 +138,6 @@ include("ul.jl") include("qr.jl") include("ql.jl") include("rq.jl") -include("choleskyinv.jl") include("polar.jl") include("reversecholesky.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 176c10b..77d1c66 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -439,34 +439,5 @@ end end include("test_rq.jl") - -# choleskyinv.jl -@testset "choleskyinv" begin - etol = 1e-9 - n = 20 - # real matrices - A = randn(n, n) - P = A*A' - Pi = inv(P) - C = choleskyinv!(copy(Matrix(P))) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < etol) - C = choleskyinv(P) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < etol) - - # repeat the test for complex matrices - A=randn(ComplexF64, n, n) - P=A*A' - Pi=inv(P) - C=choleskyinv!(copy(Matrix(P))) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < 4etol) - C = choleskyinv(P) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < 4etol) -end - - include("test_polar.jl") include("test_reversecholesky.jl") \ No newline at end of file From 4c7c2c8a8ac6aa31c1dd468a2eb0518621fc50d5 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 8 Jul 2023 20:47:14 +0100 Subject: [PATCH 05/10] use col/rowsupport --- src/reversecholesky.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index 9ed730e..251debe 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -45,8 +45,10 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) realdiag = eltype(A) <: Complex @inbounds begin for k = n:-1:1 + cs = colsupport(A, k) + a,b = first(cs),last(cs) Akk = realdiag ? real(A[k,k]) : A[k,k] - for j = k+1:n + for j = max(a,k+1):min(n,b) Akk -= realdiag ? abs2(A[k,j]) : A[k,j]'A[k,j] end A[k,k] = Akk @@ -56,12 +58,12 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) end A[k,k] = Akk AkkInv = inv(Akk) - for j = k+1:n - @simd for i = 1:k-1 + for j = (k+1:n) ∩ rowsupport(A, k) + @simd for i = (1:k-1) ∩ colsupport(A,j) A[i,k] -= A[i,j]*A[k,j]' end end - for i = 1:k-1 + for i = (1:k-1) ∩ cs A[i,k] *= AkkInv' end end @@ -131,6 +133,8 @@ function reversecholesky!(A::AbstractMatrix, ::NoPivot = NoPivot(); check::Bool end end +reversecholcopy(A) = cholcopy(A) + # reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric # or Hermitian matrix ## No pivoting (default) @@ -145,10 +149,10 @@ The triangular ReverseCholesky factor can be obtained from the factorization `F` where `A ≈ F.U * F.U' ≈ F.L' * F.L`. """ reversecholesky(A::AbstractMatrix, ::NoPivot=NoPivot(); check::Bool = true) = - reversecholesky!(cholcopy(A); check) + reversecholesky!(reversecholcopy(A); check) function reversecholesky(A::AbstractMatrix{Float16}, ::NoPivot=NoPivot(); check::Bool = true) - X = reversecholesky!(cholcopy(A); check = check) + X = reversecholesky!(reversecholcopy(A); check = check) return ReverseCholesky{Float16}(X) end From eed433369203fdf126b11965138b76b00a637508 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 9 Jul 2023 11:35:20 +0100 Subject: [PATCH 06/10] adjust supports --- src/reversecholesky.jl | 38 +++++++++++++++++++++++++++++------- test/test_reversecholesky.jl | 17 +++++++++++++--- 2 files changed, 45 insertions(+), 10 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index 251debe..6195a04 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -46,9 +46,9 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) @inbounds begin for k = n:-1:1 cs = colsupport(A, k) - a,b = first(cs),last(cs) + rs = rowsupport(A, k) Akk = realdiag ? real(A[k,k]) : A[k,k] - for j = max(a,k+1):min(n,b) + for j = (k+1:n) ∩ rs Akk -= realdiag ? abs2(A[k,j]) : A[k,j]'A[k,j] end A[k,k] = Akk @@ -58,8 +58,8 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) end A[k,k] = Akk AkkInv = inv(Akk) - for j = (k+1:n) ∩ rowsupport(A, k) - @simd for i = (1:k-1) ∩ colsupport(A,j) + for j = (k+1:n) ∩ rs + @simd for i = (1:k-1) ∩ cs ∩ colsupport(A,j) A[i,k] -= A[i,j]*A[k,j]' end end @@ -77,8 +77,10 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{LowerTriangular}) realdiag = eltype(A) <: Complex @inbounds begin for k = n:-1:1 + cs = colsupport(A, k) + rs = rowsupport(A, k) Akk = realdiag ? real(A[k,k]) : A[k,k] - for i = k+1:n + for i = (k+1:n) ∩ cs Akk -= realdiag ? abs2(A[i,k]) : A[i,k]'A[i,k] end A[k,k] = Akk @@ -88,8 +90,8 @@ function _reverse_chol!(A::AbstractMatrix, ::Type{LowerTriangular}) end A[k,k] = Akk AkkInv = inv(copy(Akk')) - for j = 1:k-1 - for i = k+1:n + for j = (1:k-1) ∩ rs # colsupport == rowsupport + for i = (k+1:n) ∩ cs ∩ colsupport(A,j) A[k,j] -= A[i,k]'*A[i,j] end A[k,j] = AkkInv*A[k,j] @@ -134,6 +136,10 @@ function reversecholesky!(A::AbstractMatrix, ::NoPivot = NoPivot(); check::Bool end reversecholcopy(A) = cholcopy(A) +function reversecholcopy(A::SymTridiagonal) + T = LinearAlgebra.choltype(A) + Symmetric(Bidiagonal(AbstractVector{T}(A.dv), AbstractVector{T}(A.ev), :U)) +end # reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric # or Hermitian matrix @@ -263,4 +269,22 @@ function getproperty(C::ReverseCholesky{<:Any,<:Diagonal}, d::Symbol) end end +function getproperty(C::ReverseCholesky{<:Any, <:Bidiagonal}, d::Symbol) + Cfactors = getfield(C, :factors) + Cuplo = getfield(C, :uplo) + if d === :U && Cfactors.uplo === Cuplo === 'U' + return Cfactors + elseif d === :L && Cfactors.uplo === Cuplo === 'U' + return Cfactors' + elseif d === :U && Cfactors.uplo === Cuplo === 'L' + return Cfactors' + elseif d === :L && Cfactors.uplo === Cuplo === 'L' + return Cfactors + elseif d === :UL + return Cfactors + else + return getfield(C, d) + end +end + inv(C::ReverseCholesky{<:Any,<:Diagonal}) = Diagonal(map(inv∘abs2, C.factors.diag)) \ No newline at end of file diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index 2ce9b98..888e6f9 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -297,7 +297,8 @@ end # complex, failing D[2, 2] = 0.0 + 0im - @test_throws ArgumentError reversecholesky(D) + @test reversecholesky(D).U ≈ sqrt.(D) + # InexactError for Int @test_throws InexactError reversecholesky!(Diagonal([2, 1])) @@ -319,8 +320,6 @@ end @test ReverseCholesky(factors, uplo, Int64(info)) == chol end - - @testset "issue #37356, diagonal elements of hermitian generic matrix" begin B = Hermitian(hcat([one(BigFloat) + im])) @test Matrix(reversecholesky(B)) ≈ B @@ -367,4 +366,16 @@ end @test B.UL ≈ B32.UL @test Matrix(B) ≈ A end + + @testset "large Sparse" begin + n = 1_000_000 + A = SymTridiagonal(fill(4,n), fill(1,n-1)) + R = reversecholesky(A) + @test R isa ReverseCholesky{Float64, <:Bidiagonal} + @test R.U isa Bidiagonal + @test R.L isa Bidiagonal + @test R.UL isa Bidiagonal + # Bidiagonal multiplication not supported + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + end end \ No newline at end of file From f17fc3dbd63943271fde8c3d1dc40944786891a3 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 10 Jul 2023 17:26:24 +0100 Subject: [PATCH 07/10] add tests, make copy for bidiagonal case --- src/reversecholesky.jl | 22 ++++++++++++++++++++-- test/test_reversecholesky.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index 6195a04..b7dee24 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -34,8 +34,8 @@ ReverseCholesky(U::UpperTriangular{T}) where {T} = ReverseCholesky{T,typeof(U.da ReverseCholesky(L::LowerTriangular{T}) where {T} = ReverseCholesky{T,typeof(L.data)}(L.data, 'L', 0) # iteration for destructuring into components -Base.iterate(C::ReverseCholesky) = (C.L, Val(:U)) -Base.iterate(C::ReverseCholesky, ::Val{:U}) = (C.U, Val(:done)) +Base.iterate(C::ReverseCholesky) = (C.U, Val(:U)) +Base.iterate(C::ReverseCholesky, ::Val{:U}) = (C.L, Val(:done)) Base.iterate(C::ReverseCholesky, ::Val{:done}) = nothing ## Non BLAS/LAPACK element types (generic) @@ -141,6 +141,24 @@ function reversecholcopy(A::SymTridiagonal) Symmetric(Bidiagonal(AbstractVector{T}(A.dv), AbstractVector{T}(A.ev), :U)) end +function reversecholcopy(A::Symmetric{<:Any,<:Tridiagonal}) + T = LinearAlgebra.choltype(A) + if A.uplo == 'U' + Symmetric(Bidiagonal(AbstractVector{T}(parent(A).d), AbstractVector{T}(parent(A).du), :U)) + else + Symmetric(Bidiagonal(AbstractVector{T}(parent(A).d), AbstractVector{T}(parent(A).dl), :L), :L) + end +end + +_copyifsametype(::Type{T}, A::AbstractMatrix{T}) where T = copy(A) +_copyifsametype(_, A) = A + +function reversecholcopy(A::Symmetric{<:Any,<:Bidiagonal}) + T = LinearAlgebra.choltype(A) + B = _copyifsametype(T, AbstractMatrix{T}(parent(A))) + Symmetric{T,typeof(B)}(B, A.uplo) +end + # reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric # or Hermitian matrix ## No pivoting (default) diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index 888e6f9..ddd5a12 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -377,5 +377,31 @@ end @test R.UL isa Bidiagonal # Bidiagonal multiplication not supported @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + + A = Tridiagonal(fill(1.0,n-1), fill(4.0,n), fill(1/2,n-1)) + R = reversecholesky(Symmetric(A)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] + R = reversecholesky(Symmetric(A,:L)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + + A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :U) + R = reversecholesky(Symmetric(A)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] + R = reversecholesky(Symmetric(A,:L)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + + A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :L) + R = reversecholesky(Symmetric(A,:L)) + @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] + end + + @testset "coverage" begin + A = [4 2; 2 4] + R = reversecholesky(A) + @test ReverseCholesky(R.factors, R.uplo, R.info) == R + U,L = R + @test U*U' ≈ A + @test reversecholesky(5).U[1,1] ≈ sqrt(5) + @test propertynames(R) == (:U, :L, :UL) end end \ No newline at end of file From c91afe6d8a2ff6e1d3aa1016aae9a86dff2c40ab Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 10 Jul 2023 17:48:44 +0100 Subject: [PATCH 08/10] fix tests --- src/reversecholesky.jl | 13 +++++++++---- test/test_reversecholesky.jl | 7 +++++-- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index b7dee24..63a9a97 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -290,13 +290,18 @@ end function getproperty(C::ReverseCholesky{<:Any, <:Bidiagonal}, d::Symbol) Cfactors = getfield(C, :factors) Cuplo = getfield(C, :uplo) - if d === :U && Cfactors.uplo === Cuplo === 'U' + @assert Cfactors.uplo === Cuplo + if d === :U && Cuplo === 'U' return Cfactors - elseif d === :L && Cfactors.uplo === Cuplo === 'U' + elseif d === :L && Cuplo === 'U' return Cfactors' - elseif d === :U && Cfactors.uplo === Cuplo === 'L' + elseif d === :U && Cuplo === 'L' return Cfactors' - elseif d === :L && Cfactors.uplo === Cuplo === 'L' + elseif d === :L && Cuplo === 'L' + return Cfactors + elseif d === :U && Cuplo === 'L' + return Cfactors' + elseif d === :L && Cuplo === 'L' return Cfactors elseif d === :UL return Cfactors diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index ddd5a12..e0d1ee5 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -377,22 +377,25 @@ end @test R.UL isa Bidiagonal # Bidiagonal multiplication not supported @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] A = Tridiagonal(fill(1.0,n-1), fill(4.0,n), fill(1/2,n-1)) R = reversecholesky(Symmetric(A)) @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[1,:] R = reversecholesky(Symmetric(A,:L)) @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :U) R = reversecholesky(Symmetric(A)) @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] - R = reversecholesky(Symmetric(A,:L)) - @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[1,:] A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :L) R = reversecholesky(Symmetric(A,:L)) @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] end @testset "coverage" begin From 5ed00ce01256b5190c8f2ddbb4acebcb76c00b13 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 10 Jul 2023 20:10:03 +0100 Subject: [PATCH 09/10] Update test_reversecholesky.jl --- test/test_reversecholesky.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl index e0d1ee5..9be9c7b 100644 --- a/test/test_reversecholesky.jl +++ b/test/test_reversecholesky.jl @@ -396,6 +396,11 @@ end R = reversecholesky(Symmetric(A,:L)) @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] + + A = Bidiagonal(fill(4,n), fill(1,n-1), :L) + R = reversecholesky(Symmetric(A,:L)) + @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] end @testset "coverage" begin From cd6e2ffbbeeb3efd2baeac23e634751cfd02ffc3 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 10 Jul 2023 20:11:29 +0100 Subject: [PATCH 10/10] Update reversecholesky.jl --- src/reversecholesky.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl index 63a9a97..a04802f 100644 --- a/src/reversecholesky.jl +++ b/src/reversecholesky.jl @@ -299,10 +299,6 @@ function getproperty(C::ReverseCholesky{<:Any, <:Bidiagonal}, d::Symbol) return Cfactors' elseif d === :L && Cuplo === 'L' return Cfactors - elseif d === :U && Cuplo === 'L' - return Cfactors' - elseif d === :L && Cuplo === 'L' - return Cfactors elseif d === :UL return Cfactors else