From b695c7a1131fe898f4297e02160de2d2c6754b9f Mon Sep 17 00:00:00 2001 From: Andreas Noack Jensen Date: Fri, 21 Feb 2014 19:28:01 +0100 Subject: [PATCH] A more robust error bound for the determinant. Remove Float16 from linalg tests. Add fallback condition number method. --- base/linalg/factorization.jl | 3 ++- test/linalg1.jl | 22 +++++++++++++--------- test/linalg2.jl | 18 +++++++++--------- 3 files changed, 24 insertions(+), 19 deletions(-) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 479b1470042af..c710c1fe81842 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -250,7 +250,8 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular inv{T<:BlasFloat}(A::LU{T})=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info -cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))) +cond{T<:BlasFloat}(A::LU{T}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))) +cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p) #################### # QR Factorization # diff --git a/test/linalg1.jl b/test/linalg1.jl index ba59b8842555c..a2956c47b3ff1 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -20,14 +20,16 @@ areal = randn(n,n)/2 aimg = randn(n,n)/2 breal = randn(n,2)/2 bimg = randn(n,2)/2 -for eltya in (Float16, Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int) - for eltyb in (Float16, Float32, Float64, Complex32, Complex64, Complex128, Int) - a = eltya == Int ? rand(1:5, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) +for eltya in (Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int) + for eltyb in (Float32, Float64, Complex32, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) asym = a'+a # symmetric indefinite apd = a'*a # symmetric positive-definite - ε = max(eps(abs(float(one(eltya)))),eps(abs(float(one(eltyb))))) + εa = eps(abs(float(one(eltya)))) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") @@ -50,8 +52,8 @@ debug && println("(Automatic) upper Cholesky factor") @test norm(apd*x-b)/norm(b) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ @test_approx_eq apd * inv(capd) eye(n) - @test_approx_eq_eps a*(capd\(a'*b)) b 8000ε - @test_approx_eq det(capd) det(apd) + @test norm(a*(capd\(a'*b)) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit + @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow debug && println("lower Cholesky factor") @@ -64,7 +66,7 @@ debug && println("pivoted Choleksy decomposition") cpapd = cholfact(apd, pivot=true) @test rank(cpapd) == n @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing - @test_approx_eq_eps b apd * (cpapd\b) 15000ε + @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit if isreal(apd) @test_approx_eq apd * inv(cpapd) eye(n) end @@ -83,12 +85,13 @@ debug && println("Bunch-Kaufman factors of a pos-def matrix") end debug && println("(Automatic) Square LU decomposition") + κ = cond(a,1) lua = factorize(a) l,u,p = lua[:L], lua[:U], lua[:p] @test_approx_eq l*u a[p,:] @test_approx_eq l[invperm(p),:]*u a @test_approx_eq a * inv(lua) eye(n) - @test_approx_eq_eps a*(lua\b) b 80ε + @test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns debug && println("Thin LU") lua = lufact(a[:,1:5]) @@ -188,10 +191,11 @@ debug && println("Generalized svd") end debug && println("Solve square general system of equations") + κ = cond(a,1) x = a \ b - @test_approx_eq_eps a*x b 80ε @test_throws b'\b @test_throws b\b' + @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit! debug && println("Solve upper triangular system") x = triu(a) \ b diff --git a/test/linalg2.jl b/test/linalg2.jl index 9ddb296788171..3e1410196c7c6 100644 --- a/test/linalg2.jl +++ b/test/linalg2.jl @@ -43,6 +43,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) V = convert(Matrix{elty}, V) C = convert(Matrix{elty}, C) end + ε = eps(abs2(float(one(elty)))) T = Tridiagonal(dl, d, du) @test size(T, 1) == n @test size(T) == (n, n) @@ -95,8 +96,8 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) F = full(W) @test_approx_eq W*v F*v iFv = F\v - @test_approx_eq W\v iFv - @test_approx_eq det(W) det(F) + @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong + @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong iWv = similar(iFv) if elty != Int Base.LinAlg.solve!(iWv, W, v) @@ -208,7 +209,7 @@ end #Triangular matrices n=12 -for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) A = convert(Matrix{elty}, randn(n, n)) b = convert(Matrix{elty}, randn(n, 2)) if elty <: Complex @@ -245,7 +246,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt end #Tridiagonal matrices -for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty}) +for relty in (Float32, Float64), elty in (relty, Complex{relty}) a = convert(Vector{elty}, randn(n-1)) b = convert(Vector{elty}, randn(n)) c = convert(Vector{elty}, randn(n-1)) @@ -263,11 +264,10 @@ for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty}) end #SymTridiagonal (symmetric tridiagonal) matrices -for relty in (Float16, Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work +for relty in (Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work a = convert(Vector{elty}, randn(n)) b = convert(Vector{elty}, randn(n-1)) if elty <: Complex - relty==Float16 && continue a += im*convert(Vector{elty}, randn(n)) b += im*convert(Vector{elty}, randn(n-1)) end @@ -306,7 +306,7 @@ for elty in (Float32, Float64) end #Bidiagonal matrices -for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) dv = convert(Vector{elty}, randn(n)) ev = convert(Vector{elty}, randn(n-1)) b = convert(Matrix{elty}, randn(n, 2)) @@ -361,7 +361,7 @@ end #Diagonal matrices n=12 -for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) +for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty}) d=convert(Vector{elty}, randn(n)) v=convert(Vector{elty}, randn(n)) U=convert(Matrix{elty}, randn(n,n)) @@ -529,7 +529,7 @@ end nnorm = 1000 mmat = 100 nmat = 80 -for elty in (Float16, Float32, Float64, BigFloat, Complex{Float16}, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt) +for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt) debug && println(elty) ## Vector