Skip to content

Instantly share code, notes, and snippets.

@mschauer
Last active August 29, 2015 14:20
Show Gist options
  • Save mschauer/c325ff89cd378fe7ecd6 to your computer and use it in GitHub Desktop.
Save mschauer/c325ff89cd378fe7ecd6 to your computer and use it in GitHub Desktop.
Shape preservation for blocked cholesky
import Base: ctranspose, *, -, *, \, zero, showlimited, print, show, writemime, showcompact, showcompact_lim, inv, one
import Base.display
import Base.LinAlg.chol!
import Base.LinAlg.chol
const Dadd = Dict(
('L','L') => 'L',
('U','U') => 'U',
('L','U') => 'A',
('A','L') => 'A',
('A','U') => 'A',
('A','A') => 'A',
('0','L') => 'L',
('0','U') => 'U',
('0','A') => 'A',
('0','F') => 'F',
('0','0') => '0',
('A', 'F') => 'F',
('F', 'F') => 'F'
)
const Dmul = Dict(
('L','L') => 'L',
('U','U') => 'U',
('L','U') => 'A',
('L','A') => 'A',
('A','U') => 'A',
('A', '0') => '0',
('F', '0') => '0',
('L', '0') => '0',
('U', '0') => '0',
('0', '0') => '0',
('0', 'A') => '0',
('0', 'L') => '0',
('0', 'F') => '0',
('0', 'U') => '0',
('F', 'F') => 'F',
('A', 'A') => 'F',
)
const Dldiv = Dict(
('L','A') => 'A',
('L','0') => '0',
('L','L') => 'L',
('L','F') => 'F',
('F','F') => 'F',
('F','0') => '0',
('F','1') => 'F',
('L','1') => 'L',
('A','1') => 'F',
('U','1') => 'U',
('A','A') => 'A',
)
type Ring
x::Float64
c::Char
end
zero(::Type{Ring}) = Ring(0.0, '0')
one(::Type{Ring}) = Ring(1, 'A')
print(io::IO, r::Ring) = print(io,r.c)
showlimited(io::IO, r::Ring) = print(io, r)
writemime(io::IO, ::MIME"text/plain", r::Ring) = print(io, r)
showcompact_lim(io::IO, r::Ring) = print(io, r)
showcompact(io::IO, r::Ring) = print(io, r)
ctranspose(r::Ring) = Ring(r.x, Dict('L'=>'U', 'U'=>'L', 'A'=>'A', '0'=>'0', 'F'=>'F')[r.c])
+(r1::Ring, r2::Ring) = Ring(r1.x + r2.x, Dadd[minmax(r1.c, r2.c)])
-(r1::Ring, r2::Ring) = Ring(r1.x - r2.x, Dadd[minmax(r1.c, r2.c)])
*(r1::Ring, r2::Ring) = Ring(r1.x * r2.x, Dmul[r1.c, r2.c])
*(f1::Float64, r2::Ring) = Ring(f1 * r2.x, r2.c)
\(r1::Ring, r2::Ring) = Ring(r1.x \ r2.x, Dldiv[r1.c, r2.c])
inv(r::Ring) = Ring(inv(r.x), Dldiv[r.c, '1'])
function root(r::Ring)
# assert(r.c == 'A' || r.c == '0')
Ring(sqrt(r.x), Dict('A'=>'U', '0'=>'0', 'F'=>'F')[r.c])
end
chol(r::Ring, v) = chol!(r, v)
chol(r::Ring) = chol!(r, Val{:L})
function chol!(r::Ring, ::Type{Val{:U}})
Ring(sqrt(r.x), Dict('A'=>'U', '0'=>'0', 'F'=>'F')[r.c])
end
function chol!(r::Ring, ::Type{Val{:L}})
Ring(sqrt(r.x), Dict('A'=>'L', '0'=>'0', 'F'=>'F')[r.c])
end
L = Ring[Ring(3.,'L') Ring(1.,'L')
Ring(1.,'L') Ring(3.,'L')]
A = L*(L')
chol!(copy(A))
F = Float64[
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
1 2 1 1 1 1 1 1 1 2 1 1 2 1 9 0 0
2 1 1 1 1 1 1 1 2 1 1 2 1 2 7 9 0
1 1 1 1 1 1 2 1 1 2 1 9 1 1 7 4 9
]
n = size(F,1)
R = randn(n,n)
F = R .* F
U = F * F' + n*eye(n)
A2 = map(x-> Ring(x, x != 0 ? 'A':'0'), U)
L2 = chol!(copy(A2))
display(A2)
display(L2)
UU = rand(n*n, n*n)
1(abs(chol!(UU*UU' .* kron(U, U) + kron(U,U) + 100*n*n*eye(n*n))) .> eps() )
@andreasnoack
Copy link

Thanks for the link. I like the example and surprisingly enough, in this week, I've actually looked at the shape preserving property of the arrow matrix under Cholesky/LDLt factorizations. We might need the backward Cholesky at some point such that we can get the same property when the arrowhead is to the upper left.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment