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

A more robust error bound for the determinant. #5902

Merged
merged 1 commit into from
Mar 8, 2014
Merged

A more robust error bound for the determinant. #5902

Show file tree
Hide file tree
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
3 changes: 2 additions & 1 deletion base/linalg/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 #
Expand Down
22 changes: 13 additions & 9 deletions test/linalg1.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand All @@ -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")
Expand All @@ -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
Expand All @@ -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])
Expand Down Expand Up @@ -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
Expand Down
18 changes: 9 additions & 9 deletions test/linalg2.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down