linalg1.jl
debug = false
import Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted
n = 10
srand(1234321)
a = rand(n,n)
for elty in (Float32, Float64, Complex64, Complex128)
a = convert(Matrix{elty}, a)
# cond
@test_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01
@test_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01
@test_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01
@test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01
end
areal = randn(n,n)/2
aimg = randn(n,n)/2
breal = randn(n,2)/2
bimg = randn(n,2)/2
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
ε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")
debug && println("(Automatic) upper Cholesky factor")
if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement cholesky decomposition in julia
capd = factorize(apd)
r = capd[:U]
κ = cond(apd) #condition number
#Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1
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
#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)/norm(x) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test norm(apd*x-b)/norm(b) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test_approx_eq apd * inv(capd) eye(n)
@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")
l = cholfact(apd, :L)[:L]
@test_approx_eq l*l' apd
end
debug && println("pivoted Choleksy decomposition")
if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted cholesky decomposition in julia
cpapd = cholfact(apd, pivot=true)
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
@test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
if isreal(apd)
@test_approx_eq apd * inv(cpapd) eye(n)
end
end
debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix")
if eltya != BigFloat && eltyb != BigFloat # Not implemented for BigFloat and I don't think it will.
bc1 = factorize(asym)
@test_approx_eq inv(bc1) * asym eye(n)
@test_approx_eq_eps asym * (bc1\b) b 600ε
debug && println("Bunch-Kaufman factors of a pos-def matrix")
bc2 = bkfact(apd)
@test_approx_eq inv(bc2) * apd eye(n)
@test_approx_eq_eps apd * (bc2\b) b 60000ε
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 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])
@test_approx_eq lua[:L]*lua[:U] lua[:P]*a[:,1:5]
debug && println("Fat LU")
lua = lufact(a[1:5,:])
@test_approx_eq lua[:L]*lua[:U] lua[:P]*a[1:5,:]
debug && println("QR decomposition (without pivoting)")
qra = qrfact(a, pivot=false)
q,r = qra[:Q], qra[:R]
@test_approx_eq q'*full(q, thin=false) eye(n)
@test_approx_eq q*full(q, thin=false)' eye(n)
@test_approx_eq q*r a
@test_approx_eq_eps a*(qra\b) b 3000ε
debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
qrpa = factorize(a[1:5,:])
q,r = qrpa[:Q], qrpa[:R]
if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
@test_approx_eq q'*full(q, thin=false) eye(5)
@test_approx_eq q*full(q, thin=false)' eye(5)
@test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:5,p] : a[1:5,:]
@test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:5,:]
@test_approx_eq_eps a[1:5,:]*(qrpa\b[1:5]) b[1:5] 5000ε
debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
qrpa = factorize(a[:,1:5])
q,r = qrpa[:Q], qrpa[:R]
if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
@test_approx_eq q'*full(q, thin=false) eye(n)
@test_approx_eq q*full(q, thin=false)' eye(n)
@test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:5]
@test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:5]
debug && println("symmetric eigen-decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
d,v = eig(asym)
@test_approx_eq asym*v[:,1] d[1]*v[:,1]
@test_approx_eq v*scale(d,v') asym
@test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
end
debug && println("non-symmetric eigen decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
d,v = eig(a)
for i in 1:size(a,2) @test_approx_eq a*v[:,i] d[i]*v[:,i] end
end
debug && println("symmetric generalized eigenproblem")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
a610 = a[:,6:10]
f = eigfact(asym[1:5,1:5], a610'a610)
@test_approx_eq asym[1:5,1:5]*f[:vectors] scale(a610'a610*f[:vectors], f[:values])
@test_approx_eq f[:values] eigvals(asym[1:5,1:5], a610'a610)
@test_approx_eq_eps prod(f[:values]) prod(eigvals(asym[1:5,1:5]/(a610'a610))) 200ε
end
debug && println("Non-symmetric generalized eigenproblem")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
f = eigfact(a[1:5,1:5], a[6:10,6:10])
@test_approx_eq a[1:5,1:5]*f[:vectors] scale(a[6:10,6:10]*f[:vectors], f[:values])
@test_approx_eq f[:values] eigvals(a[1:5,1:5], a[6:10,6:10])
@test_approx_eq_eps prod(f[:values]) prod(eigvals(a[1:5,1:5]/a[6:10,6:10])) 50000ε
end
debug && println("Schur")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
f = schurfact(a)
@test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a
@test_approx_eq sort(real(f[:values])) sort(real(d))
@test_approx_eq sort(imag(f[:values])) sort(imag(d))
@test istriu(f[:Schur]) || iseltype(a,Real)
end
debug && println("Generalized Schur")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
f = schurfact(a[1:5,1:5], a[6:10,6:10])
@test_approx_eq f[:Q]*f[:S]*f[:Z]' a[1:5,1:5]
@test_approx_eq f[:Q]*f[:T]*f[:Z]' a[6:10,6:10]
@test istriu(f[:S]) || iseltype(a,Real)
@test istriu(f[:T]) || iseltype(a,Real)
end
debug && println("singular value decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
usv = svdfact(a)
@test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a
end
debug && println("Generalized svd")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
gsvd = svdfact(a,a[1:5,:])
@test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a
@test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a[1:5,:]
end
debug && println("Solve square general system of equations")
κ = cond(a,1)
x = a \ b
@test_throws DimensionMismatch b'\b
@test_throws DimensionMismatch 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
#Test forward error [JIN 5705] if this is not a BigFloat
γ = n*ε/(1-n*ε)
if eltya != BigFloat
bigA = big(triu(a))
̂x = bigA \ b
for i=1:size(b, 2)
@test norm(̂x[:,i]-x[:,i], Inf)/norm(x[:,i], Inf) <= abs(condskeel(bigA, x[:,i])*γ/(1-condskeel(bigA)*γ))
end
end
#Test backward error [JIN 5705]
for i=1:size(b, 2)
@test norm(abs(b[:,i] - triu(a)*x[:,i]), Inf) <= γ * norm(triu(a), Inf) * norm(x[:,i], Inf)
end
debug && println("Solve lower triangular system")
x = tril(a)\b
#Test forward error [JIN 5705] if this is not a BigFloat
γ = n*ε/(1-n*ε)
if eltya != BigFloat
bigA = big(tril(a))
̂x = bigA \ b
for i=1:size(b, 2)
@test norm(̂x[:,i]-x[:,i], Inf)/norm(x[:,i], Inf) <= abs(condskeel(bigA, x[:,i])*γ/(1-condskeel(bigA)*γ))
end
end
#Test backward error [JIN 5705]
for i=1:size(b, 2)
@test norm(abs(b[:,i] - tril(a)*x[:,i]), Inf) <= γ * norm(tril(a), Inf) * norm(x[:,i], Inf)
end
debug && println("Test null")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
a15null = null(a[:,1:5]')
@test rank([a[:,1:5] a15null]) == 10
@test_approx_eq_eps norm(a[:,1:5]'a15null, Inf) zero(eltya) 300ε
@test_approx_eq_eps norm(a15null'a[:,1:5], Inf) zero(eltya) 400ε
@test size(null(b), 2) == 0
end
debug && println("Test pinv")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
pinva15 = pinv(a[:,1:5])
@test_approx_eq a[:,1:5]*pinva15*a[:,1:5] a[:,1:5]
@test_approx_eq pinva15*a[:,1:5]*pinva15 pinva15
end
# if isreal(a)
debug && println("Matrix square root")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
asq = sqrtm(a)
@test_approx_eq asq*asq a
asymsq = sqrtm(asym)
@test_approx_eq asymsq*asymsq asym
end
end
end