Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use testsets in cholesky test #22138

Merged
merged 1 commit into from
May 31, 2017
Merged

Use testsets in cholesky test #22138

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
324 changes: 159 additions & 165 deletions test/linalg/cholesky.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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])
Expand All @@ -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'
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down