diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index b6d92d89ba9f7..0364b4bfb9fa2 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -1,165 +1,157 @@ # This file is a part of Julia. License is MIT: https://julialang.org/license -debug = false - using Base.Test using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException -n = 10 +@testset "core functionality" begin + n = 10 -# Split n into 2 parts for tests needing two matrices -n1 = div(n, 2) -n2 = 2*n1 + # Split n into 2 parts for tests needing two matrices + n1 = div(n, 2) + n2 = 2*n1 -srand(1234321) + srand(1234321) -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 + 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 -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 + 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)))) + apds = Symmetric(apd) + apdsL = Symmetric(apd, :L) + apdh = Hermitian(apd) + apdhL = Hermitian(apd, :L) + ε = εa = eps(abs(float(one(eltya)))) - @inferred cholfact(apd) - @inferred chol(apd) - capd = factorize(apd) - r = capd[:U] - κ = cond(apd, 1) #condition number + @inferred cholfact(apd) + @inferred chol(apd) + capd = factorize(apd) + r = capd[:U] + κ = cond(apd, 1) #condition number - #getindex - @test_throws KeyError capd[:Z] + #getindex + @test_throws KeyError capd[:Z] - #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 + #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 = Complex64 and eltyb = Int. The E[i,j] had NaN32 elements - #but only with srand(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 - full(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 apd*inv(capd) ≈ eye(n) - @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit - @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow - - apos = apd[1,1] # test chol(x::Number), needs x>0 - @test all(x -> x ≈ √apos, cholfact(apos).factors) - @test_throws PosDefException chol(-one(eltya)) - - if eltya <: Real - capds = cholfact(apds) - @test inv(capds)*apds ≈ eye(n) - @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n - if eltya <: BlasReal - capds = cholfact!(copy(apds)) + #these tests were failing on 64-bit linux when inside the inner loop + #for eltya = Complex64 and eltyb = Int. The E[i,j] had NaN32 elements + #but only with srand(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 - full(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 apd*inv(capd) ≈ eye(n) + @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit + @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow + + apos = apd[1,1] # test chol(x::Number), needs x>0 + @test all(x -> x ≈ √apos, cholfact(apos).factors) + @test_throws PosDefException chol(-one(eltya)) + + if eltya <: Real + capds = cholfact(apds) @test inv(capds)*apds ≈ eye(n) @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n + if eltya <: BlasReal + capds = cholfact!(copy(apds)) + @test inv(capds)*apds ≈ eye(n) + @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n + end + ulstring = sprint(show,capds[:UL]) + @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring" + else + capdh = cholfact(apdh) + @test inv(capdh)*apdh ≈ eye(n) + @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n + capdh = cholfact!(copy(apdh)) + @test inv(capdh)*apdh ≈ eye(n) + @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n + 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 - ulstring = sprint(show,capds[:UL]) - @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring" - else - capdh = cholfact(apdh) - @test inv(capdh)*apdh ≈ eye(n) - @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n - capdh = cholfact!(copy(apdh)) - @test inv(capdh)*apdh ≈ eye(n) - @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n - 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 - # test chol of 2x2 Strang matrix - S = convert(AbstractMatrix{eltya},full(SymTridiagonal([2,2],[-1]))) - 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] - if eltya <: Real - capds = cholfact(apds) - lapds = cholfact(apdsL) - cl = chol(apdsL) - ls = lapds[:L] - @test ls*ls' ≈ apd - @test triu(capds.factors) ≈ lapds[:U] - @test tril(lapds.factors) ≈ capds[:L] - @test istriu(cl) - @test cl'cl ≈ apds - @test cl'cl ≈ apdsL - else - capdh = cholfact(apdh) - lapdh = cholfact(apdhL) - cl = chol(apdhL) - ls = lapdh[:L] - @test ls*ls' ≈ apd - @test triu(capdh.factors) ≈ lapdh[:U] - @test tril(lapdh.factors) ≈ capdh[:L] - @test istriu(cl) - @test cl'cl ≈ apdh - @test cl'cl ≈ apdhL - end + # test chol of 2x2 Strang matrix + S = convert(AbstractMatrix{eltya},full(SymTridiagonal([2,2],[-1]))) + 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] + if eltya <: Real + capds = cholfact(apds) + lapds = cholfact(apdsL) + cl = chol(apdsL) + ls = lapds[:L] + @test ls*ls' ≈ apd + @test triu(capds.factors) ≈ lapds[:U] + @test tril(lapds.factors) ≈ capds[:L] + @test istriu(cl) + @test cl'cl ≈ apds + @test cl'cl ≈ apdsL + else + capdh = cholfact(apdh) + lapdh = cholfact(apdhL) + cl = chol(apdhL) + ls = lapdh[:L] + @test ls*ls' ≈ apd + @test triu(capdh.factors) ≈ lapdh[:U] + @test tril(lapdh.factors) ≈ capdh[:L] + @test istriu(cl) + @test cl'cl ≈ apdh + @test cl'cl ≈ apdhL + end - #pivoted upper Cholesky - if eltya != BigFloat - cz = cholfact(zeros(eltya,n,n), :U, Val{true}) - @test_throws Base.LinAlg.RankDeficientException Base.LinAlg.chkfullrank(cz) - cpapd = cholfact(apd, :U, Val{true}) - @test rank(cpapd) == n - @test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing - if isreal(apd) - @test apd*inv(cpapd) ≈ eye(n) + #pivoted upper Cholesky + if eltya != BigFloat + cz = cholfact(zeros(eltya,n,n), :U, Val{true}) + @test_throws Base.LinAlg.RankDeficientException Base.LinAlg.chkfullrank(cz) + cpapd = cholfact(apd, :U, Val{true}) + @test rank(cpapd) == n + @test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing + if isreal(apd) + @test apd*inv(cpapd) ≈ eye(n) + end + @test full(cpapd) ≈ apd + #getindex + @test_throws KeyError cpapd[:Z] + + @test size(cpapd) == size(apd) + @test full(copy(cpapd)) ≈ apd + @test det(cpapd) ≈ det(apd) + @test logdet(cpapd) ≈ logdet(apd) + @test cpapd[:P]*cpapd[:L]*cpapd[:U]*cpapd[:P]' ≈ apd end - @test full(cpapd) ≈ apd - #getindex - @test_throws KeyError cpapd[:Z] - @test size(cpapd) == size(apd) - @test full(copy(cpapd)) ≈ apd - @test det(cpapd) ≈ det(apd) - @test logdet(cpapd) ≈ logdet(apd) - @test cpapd[:P]*cpapd[:L]*cpapd[:U]*cpapd[:P]' ≈ apd - end + for eltyb in (Float32, Float64, Complex64, Complex128, 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, Complex64, Complex128, 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) - -debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") - let Bs = b - for atype in ("Array", "SubArray") - if atype == "Array" - b = Bs - else - b = view(Bs, 1:n, 1) - end + 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 @@ -175,7 +167,6 @@ debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") end @test_throws DimensionMismatch lapd\RowVector(ones(n)) -debug && println("pivoted Cholesky decomposition") if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted Cholesky decomposition in julia @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit @@ -192,8 +183,7 @@ debug && println("pivoted Cholesky decomposition") end end -begin - # Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices +@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 = full(Base.LinAlg._chol!(X*X', LowerTriangular)[1]) U = full(Base.LinAlg._chol!(X*X', UpperTriangular)[1]) @@ -204,19 +194,22 @@ begin end # Test generic cholfact! -for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) - if elty <: Complex - A = complex.(randn(5,5), randn(5,5)) - else - A = randn(5,5) +@testset "generic cholfact!" begin + for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) + if elty <: Complex + A = complex.(randn(5,5), randn(5,5)) + else + A = randn(5,5) + end + A = convert(Matrix{elty}, A'A) + @test full(cholfact(A)[:L]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular)[1]) + @test full(cholfact(A)[:U]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular)[1]) end - A = convert(Matrix{elty}, A'A) - @test full(cholfact(A)[:L]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular)[1]) - @test full(cholfact(A)[:U]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular)[1]) end -# Test up- and downdates -let A = complex.(randn(10,5), randn(10, 5)), v = complex.(randn(5), randn(5)) +@testset "cholesky up- and downdates" begin + A = complex.(randn(10,5), randn(10, 5)) + v = complex.(randn(5), randn(5)) for uplo in (:U, :L) AcA = A'A BcB = AcA + v*v' @@ -230,8 +223,8 @@ let A = complex.(randn(10,5), randn(10, 5)), v = complex.(randn(5), randn(5)) end end -# issue #13243, unexpected nans in complex cholfact -let apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.3824869f0im 2.619998f0 + 1.8532984f0im -1.8306153f0 - 1.2336911f0im 0.32275113f0 + 0.015575029f0im 2.1968813f0 + 1.0640624f0im 0.27894387f0 + 0.97911835f0im 3.0476584f0 + 0.18548489f0im 0.3842994f0 + 0.7050991f0im +@testset "issue #13243, unexpected nans in complex cholfact" begin + apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.3824869f0im 2.619998f0 + 1.8532984f0im -1.8306153f0 - 1.2336911f0im 0.32275113f0 + 0.015575029f0im 2.1968813f0 + 1.0640624f0im 0.27894387f0 + 0.97911835f0im 3.0476584f0 + 0.18548489f0im 0.3842994f0 + 0.7050991f0im -0.79540455f0 - 0.7066077f0im 8.313246f0 + 0.0f0im -1.8076122f0 - 0.8882447f0im 0.47806996f0 + 0.48494184f0im 0.5096429f0 - 0.5395974f0im -0.7285097f0 - 0.10360408f0im -1.1760061f0 - 2.7146957f0im -0.4271084f0 + 0.042899966f0im -1.7228563f0 + 2.8335886f0im 1.8942566f0 + 0.6389735f0im 0.98274714f0 - 1.3824869f0im -1.8076122f0 + 0.8882447f0im 9.367975f0 + 0.0f0im -0.1838578f0 + 0.6468568f0im -1.8338387f0 + 0.7064959f0im 0.041852742f0 - 0.6556877f0im 2.5673025f0 + 1.9732997f0im -1.1148382f0 - 0.15693812f0im 2.4704504f0 - 1.0389464f0im 1.0858271f0 - 1.298006f0im 2.619998f0 - 1.8532984f0im 0.47806996f0 - 0.48494184f0im -0.1838578f0 - 0.6468568f0im 3.1117508f0 + 0.0f0im -1.956626f0 + 0.22825956f0im 0.07081801f0 - 0.31801307f0im 0.3698375f0 - 0.5400855f0im 0.80686307f0 + 1.5315914f0im 1.5649154f0 - 1.6229297f0im -0.112077385f0 + 1.2014246f0im @@ -261,22 +254,23 @@ let apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1. end end -# Fail if non-Hermitian -@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}) +@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}) +end -# Fail for non-BLAS element types -@test_throws ArgumentError cholfact!(Hermitian(rand(Float16, 5,5)), Val{true}) +@testset "fail for non-BLAS element types" begin + @test_throws ArgumentError cholfact!(Hermitian(rand(Float16, 5,5)), Val{true}) +end -@testset "throw for non positive matrix" begin +@testset "throw for non positive definite matrix" begin for T in (Float32, Float64, Complex64, Complex128) A = T[1 2; 2 1]; B = T[1, 1] C = cholfact(A) - @show typeof(A), typeof(B), typeof(C.factors) @test_throws PosDefException C\B @test_throws PosDefException det(C) @test_throws PosDefException logdet(C)