Skip to content

Support extracting factors from cholfact(::SparseMatrixCSC) #10402

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

Merged
timholy merged 2 commits into from
Jun 6, 2015
Merged

Support extracting factors from cholfact(::SparseMatrixCSC) #10402

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6b0e706
Support extracting factors from cholfact(::SparseMatrixCSC)
timholy Jun 5, 2015
208fd76
Add documentation on extracting factors from sparse cholfact/ldltfact
timholy Jun 5, 2015
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
181 changes: 171 additions & 10 deletions base/sparse/cholmod.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ import Base: (*), convert, copy, eltype, getindex, show, size,
linearindexing, LinearFast, LinearSlow

import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
cholfact, cholfact!, det, diag, ishermitian, isposdef,
cholfact, det, diag, ishermitian, isposdef,
issym, ldltfact, logdet

import Base.SparseMatrix: sparse, nnz
Expand Down Expand Up @@ -132,10 +132,10 @@ end
# Type definitions #
####################

# The three core data types for CHOLMOD: Dense, Sparse and Factor. CHOLMOD
# manages the memory, so the Julia versions only wrap a pointer to a struct.
# Therefore finalizers should be registret each time a pointer is returned from
# CHOLMOD.
# The three core data types for CHOLMOD: Dense, Sparse and Factor.
# CHOLMOD manages the memory, so the Julia versions only wrap a
# pointer to a struct. Therefore finalizers should be registered each
# time a pointer is returned from CHOLMOD.

# Dense
immutable C_Dense{T<:VTypes}
Expand Down Expand Up @@ -273,6 +273,28 @@ type Factor{Tv,Ti} <: Factorization{Tv}
p::Ptr{C_Factor{Tv,Ti}}
end

# FactorComponent, for encoding particular factors from a factorization
type FactorComponent{Tv,Ti,S} <: AbstractMatrix{Tv}
F::Factor{Tv,Ti}

function FactorComponent(F::Factor{Tv,Ti})
s = unsafe_load(F.p)
if s.is_ll != 0
S == :L || S == :U || S == :PtL || S == :UP || throw(CHOLMODException(string(S, " not supported for sparse LLt matrices; try :L, :U, :PtL, or :UP")))
else
S == :L || S == :U || S == :PtL || S == :UP ||
S == :D || S == :LD || S == :DU || S == :PtLD || S == :DUP ||
throw(CHOLMODException(string(S, " not supported for sparse LDLt matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")))
end
new(F)
end
end
function FactorComponent{Tv,Ti}(F::Factor{Tv,Ti}, sym::Symbol)
FactorComponent{Tv,Ti,sym}(F)
end

Factor(FC::FactorComponent) = Factor(FC.F)

#################
# Thin wrappers #
#################
Expand Down Expand Up @@ -478,7 +500,7 @@ for Ti in IndexTypes
s
end

function transpose{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer)
function transpose_{Tv<:VTypes}(A::Sparse{Tv,$Ti}, values::Integer)
s = Sparse(ccall((@cholmod_name("transpose", $Ti),:libcholmod), Ptr{C_Sparse{Tv,$Ti}},
(Ptr{C_Sparse{Tv,$Ti}}, Cint, Ptr{UInt8}),
A.p, values, common($Ti)))
Expand Down Expand Up @@ -637,6 +659,12 @@ for Ti in IndexTypes
finalizer(f, free!)
f
end
function factorize!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8})
@isok ccall((@cholmod_name("factorize", $Ti),:libcholmod), Cint,
(Ptr{C_Sparse{Tv,$Ti}}, Ptr{C_Factor{Tv,$Ti}}, Ptr{UInt8}),
A.p, F.p, cmmn)
F
end
function factorize_p!{Tv<:VTypes}(A::Sparse{Tv,$Ti}, β::Real, F::Factor{Tv,$Ti}, cmmn::Vector{UInt8})
# note that β is passed as a complex number (double beta[2]),
# but the CHOLMOD manual says that only beta[0] (real part) is used
Expand Down Expand Up @@ -691,6 +719,13 @@ for Ti in IndexTypes
end
end

function get_perm(F::Factor)
s = unsafe_load(F.p)
p = pointer_to_array(s.Perm, s.n, false)
p+1
end
get_perm(FC::FactorComponent) = get_perm(Factor(FC))

#########################
# High level interfaces #
#########################
Expand Down Expand Up @@ -879,9 +914,37 @@ function sparse{Ti}(A::Sparse{Complex{Float64},Ti}) # Notice! Cannot be type sta
end
return convert(Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, A)
end
sparse(L::Factor) = sparse(Sparse(L))
function sparse(F::Factor)
s = unsafe_load(F.p)
if s.is_ll != 0
L = Sparse(F)
A = sparse(L*L')
else
LD = sparse(F[:LD])
L, d = getLd!(LD)
A = scale(L, d)*L'
end
p = get_perm(F)
if p != [1:s.n;]
pinv = Array(Int, length(p))
for k = 1:length(p)
pinv[p[k]] = k
end
A = A[pinv,pinv]
Copy link
Member

@andreasnoack andreasnoack Mar 8, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not a good idea as you also anticipated in your comment. The problem LPx=b has to be solved Yy=b first and then Px=y, i.e. apply the inverse permutation to the solution. The operation A[pinv,pinv] is probably quite expensive for a sparse matrix, but destroying the triangular structure will require a new factorization before you can solve your problem.

I can see two solutions. Either F[:LP] returns a PivotedCholeskyFactor or we only support extracting F[:L] and F[:p] and it would be the responsibility of the user to apply the permutations after the solve. I prefer the last solution for simplicity for the users as well as the maintainers.

Copy link
Member Author

@timholy timholy Mar 8, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But note this function isn't about extracting individual factors: this is just so that sparse(cholfact(A)) == A.

Copy link
Member Author

@timholy timholy Mar 8, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(The ones for extracting components are just a few lines down from here, and indeed they follow the prescription you suggest.)

Copy link
Member

@andreasnoack andreasnoack Mar 8, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for not reading carefully enough. I can see that FactorComponent acts like the suggested PivotedCholeskyFactor. I'm still in doubt if it is a good idea to combine L and P in one type. It is already on the todo to make a permutation matrix type and we could support extracting P from a non-pivoted dense Cholesky where it would then be almost a noop in Px=y. Then code could still be generic for dense and sparse input.

Copy link
Member Author

@timholy timholy Mar 8, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm simply glad you're reading it at all, no apologies needed.

To be clear, a FactorComponent is any type of factor: it can represent just :L or :D (meaning, it also represents single factors) or :PtL or :LD or ...

Rather than mushing them together into a single type, I think it would be fine to have separated factors using an expanded type hierarchy: L = cholfact[:L]; P = cholfact[:P]. But I do wonder whether ready-made combinations are still desirable, in much the same way that a Cholesky is a nice convenience that guards against dumb mistakes like x = L\(L'\b) (instead of L'\(L\b) as it should be). Or stated another way, what advantages do you see that our cholfact types offer that's different from L, p = chol(A)?

Copy link
Contributor

@tkelman tkelman Jun 5, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The cholfact types don't actually require constructing the mathematically-valid representation of the factors, they can keep the internal data represented however the library chooses, right? So you don't have to construct L unless requested, which is potentially an expensive step that isn't needed for solving a system of equations.

Copy link
Member Author

@timholy timholy Jun 5, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, F[:PtL] returns an opaque object.

end
A
end

sparse(D::Dense) = sparse(Sparse(D))

function sparse{Tv,Ti}(FC::FactorComponent{Tv,Ti,:L})
F = Factor(FC)
s = unsafe_load(F.p)
s.is_ll != 0 || throw(CHOLMODException("sparse: supported only for :LD on LDLt factorizations"))
sparse(Sparse(F))
end
sparse{Tv,Ti}(FC::FactorComponent{Tv,Ti,:LD}) = sparse(Sparse(Factor(FC)))

# Calculate the offset into the stype field of the cholmod_sparse_struct and
# change the value
let offidx=findfirst(fieldnames(C_Sparse) .== :stype)
Expand All @@ -905,14 +968,25 @@ eltype{T<:VTypes}(A::Sparse{T}) = T
nnz(F::Factor) = nnz(Sparse(F))

function show(io::IO, F::Factor)
s = unsafe_load(F.p)
println(io, typeof(F))
showfactor(io, F)
end

function show(io::IO, FC::FactorComponent)
println(io, typeof(FC))
showfactor(io, Factor(FC))
end

function showfactor(io::IO, F::Factor)
s = unsafe_load(F.p)
@printf(io, "type: %12s\n", s.is_ll!=0 ? "LLt" : "LDLt")
@printf(io, "method: %10s\n", s.is_super!=0 ? "supernodal" : "simplicial")
@printf(io, "maxnnz: %10d\n", Int(s.nzmax))
@printf(io, "nnz: %13d\n", nnz(F))
end



isvalid(A::Dense) = check_dense(A)
isvalid(A::Sparse) = check_sparse(A)
isvalid(A::Factor) = check_factor(A)
Expand All @@ -936,7 +1010,22 @@ function size(F::Factor, i::Integer)
return 1
end


linearindexing(::Dense) = LinearFast()

size(FC::FactorComponent, i::Integer) = size(FC.F, i)
size(FC::FactorComponent) = size(FC.F)

ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:L}) = FactorComponent{Tv,Ti,:U}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:U}) = FactorComponent{Tv,Ti,:L}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtL}) = FactorComponent{Tv,Ti,:UP}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:UP}) = FactorComponent{Tv,Ti,:PtL}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:D}) = FC
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:LD}) = FactorComponent{Tv,Ti,:DU}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DU}) = FactorComponent{Tv,Ti,:LD}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:PtLD}) = FactorComponent{Tv,Ti,:DUP}(FC.F)
ctranspose{Tv,Ti}(FC::FactorComponent{Tv,Ti,:DUP}) = FactorComponent{Tv,Ti,:PtLD}(FC.F)

function getindex(A::Dense, i::Integer)
s = unsafe_load(A.p)
0 < i <= s.nrow*s.ncol || throw(BoundsError())
Expand All @@ -957,6 +1046,27 @@ function getindex{T}(A::Sparse{T}, i0::Integer, i1::Integer)
((r1 > r2) || (unsafe_load(s.i, r1) + 1 != i0)) ? zero(T) : unsafe_load(s.x, r1)
end

function getindex(F::Factor, sym::Symbol)
sym == :p && return get_perm(F)
FactorComponent(F, sym)
end

function getLd!(S::SparseMatrixCSC)
d = Array(eltype(S), size(S, 1))
fill!(d, 0)
col = 1
for k = 1:length(S.nzval)
while k >= S.colptr[col+1]
col += 1
end
if S.rowval[k] == col
d[col] = S.nzval[k]
S.nzval[k] = 1
end
end
S, d
end

## Multiplication
(*)(A::Sparse, B::Sparse) = ssmult(A, B, 0, true, true)
(*)(A::Sparse, B::Dense) = sdmult!(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2)))
Expand All @@ -966,7 +1076,7 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti})
cm = common(Ti)

if !is(A,B)
aa1 = transpose(B, 2)
aa1 = transpose_(B, 2)
## result of ssmult will have stype==0, contain numerical values and be sorted
return ssmult(A, aa1, 0, true, true)
end
Expand All @@ -984,7 +1094,7 @@ function A_mul_Bc{Tv<:VRealTypes,Ti<:ITypes}(A::Sparse{Tv,Ti}, B::Sparse{Tv,Ti})
end

function Ac_mul_B(A::Sparse, B::Sparse)
aa1 = transpose(A, 2)
aa1 = transpose_(A, 2)
if is(A,B)
return A_mul_Bc(aa1, aa1)
end
Expand Down Expand Up @@ -1071,6 +1181,57 @@ update!{T<:VTypes}(F::Factor{T}, A::SparseMatrixCSC{T}; kws...) = update!(F, Spa

## Solvers

for (T, f) in ((:Dense, :solve), (:Sparse, :spsolve))
@eval begin
# Solve Lx = b and L'x=b where A = L*L'
function (\){T,Ti}(L::FactorComponent{T,Ti,:L}, B::$T)
($f)(CHOLMOD_L, Factor(L), B)
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:U}, B::$T)
($f)(CHOLMOD_Lt, Factor(L), B)
end
# Solve PLx = b and L'P'x=b where A = P*L*L'*P'
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtL}, B::$T)
F = Factor(L)
($f)(CHOLMOD_L, F, ($f)(CHOLMOD_P, F, B)) # Confusingly, CHOLMOD_P solves P'x = b
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:UP}, B::$T)
F = Factor(L)
($f)(CHOLMOD_Pt, F, ($f)(CHOLMOD_Lt, F, B))
end
# Solve various equations for A = L*D*L' and A = P*L*D*L'*P'
function (\){T,Ti}(L::FactorComponent{T,Ti,:D}, B::$T)
($f)(CHOLMOD_D, Factor(L), B)
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:LD}, B::$T)
($f)(CHOLMOD_LD, Factor(L), B)
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:DU}, B::$T)
($f)(CHOLMOD_DLt, Factor(L), B)
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:PtLD}, B::$T)
F = Factor(L)
($f)(CHOLMOD_LD, F, ($f)(CHOLMOD_P, F, B))
end
function (\){T,Ti}(L::FactorComponent{T,Ti,:DUP}, B::$T)
F = Factor(L)
($f)(CHOLMOD_Pt, F, ($f)(CHOLMOD_DLt, F, B))
end
end
end

function (\)(L::FactorComponent, b::Vector)
reshape(convert(Matrix, L\Dense(b)), length(b))
end
function (\)(L::FactorComponent, B::Matrix)
convert(Matrix, L\Dense(B))
end
function (\)(L::FactorComponent, B::SparseMatrixCSC)
sparse(L\Sparse(B,0))
end

Ac_ldiv_B(L::FactorComponent, B) = ctranspose(L)\B

(\)(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B)
(\)(L::Factor, b::Vector) = reshape(convert(Matrix, solve(CHOLMOD_A, L, Dense(b))), length(b))
(\)(L::Factor, B::Matrix) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B)))
Expand Down
34 changes: 32 additions & 2 deletions doc/stdlib/linalg.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,13 +104,27 @@ Linear algebra functions in Julia are largely implemented by calling functions f

.. function:: cholfact(A; shift=0, perm=Int[]) -> CHOLMOD.Factor

Compute the Cholesky factorization of a sparse positive definite matrix ``A``. A fill-reducing permutation is used. The main application of this type is to solve systems of equations with ``\``, but also the methods ``diag``, ``det``, ``logdet`` are defined. The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
Compute the Cholesky factorization of a sparse positive definite
matrix ``A``. A fill-reducing permutation is used. ``F =
cholfact(A)`` is most frequently used to solve systems of equations
with ``F\b``, but also the methods ``diag``, ``det``, ``logdet``
are defined for ``F``. You can also extract individual factors
from ``F``, using ``F[:L]``. However, since pivoting is on by
default, the factorization is internally represented as ``A ==
P'*L*L'*P`` with a permutation matrix ``P``; using just ``L``
without accounting for ``P`` will give incorrect answers. To
include the effects of permutation, it's typically preferable to
extact "combined" factors like ``PtL = F[:PtL]`` (the equivalent of
``P'*L``) and ``LtP = F[:UP]`` (the equivalent of ``L'*P``).

Setting optional ``shift`` keyword argument computes the factorization
of ``A+shift*I`` instead of ``A``. If the ``perm`` argument is nonempty,
it should be a permutation of `1:size(A,1)` giving the ordering to use
(instead of CHOLMOD's default AMD ordering).

The function calls the C library CHOLMOD and many other functions
from the library are wrapped but not exported.

.. function:: cholfact!(A [,LU=:U [,pivot=Val{false}]][;tol=-1.0]) -> Cholesky

``cholfact!`` is the same as :func:`cholfact`, but saves space by overwriting the input ``A``, instead of creating a copy. ``cholfact!`` can also reuse the symbolic factorization from a different matrix ``F`` with the same structure when used as: ``cholfact!(F::CholmodFactor, A)``.
Expand All @@ -121,13 +135,29 @@ Linear algebra functions in Julia are largely implemented by calling functions f

.. function:: ldltfact(A; shift=0, perm=Int[]) -> CHOLMOD.Factor

Compute the LDLt factorization of a sparse symmetric or Hermitian matrix ``A``. A fill-reducing permutation is used. The main application of this type is to solve systems of equations with ``\``, but also the methods ``diag``, ``det``, ``logdet`` are defined. The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
Compute the LDLt factorization of a sparse symmetric or Hermitian
matrix ``A``. A fill-reducing permutation is used. ``F =
ldltfact(A)`` is most frequently used to solve systems of equations
with ``F\b``, but also the methods ``diag``, ``det``, ``logdet``
are defined for ``F``. You can also extract individual factors from
``F``, using ``F[:L]``. However, since pivoting is on by default,
the factorization is internally represented as ``A == P'*L*D*L'*P``
with a permutation matrix ``P``; using just ``L`` without
accounting for ``P`` will give incorrect answers. To include the
effects of permutation, it's typically preferable to extact
"combined" factors like ``PtL = F[:PtL]`` (the equivalent of
``P'*L``) and ``LtP = F[:UP]`` (the equivalent of ``L'*P``). The
complete list of supported factors is ``:L, :PtL, :D, :UP, :U, :LD,
:DU, :PtLD, :DUP``.

Setting optional ``shift`` keyword argument computes the factorization
of ``A+shift*I`` instead of ``A``. If the ``perm`` argument is nonempty,
it should be a permutation of `1:size(A,1)` giving the ordering to use
(instead of CHOLMOD's default AMD ordering).

The function calls the C library CHOLMOD and many other functions
from the library are wrapped but not exported.

.. function:: qr(A [,pivot=Val{false}][;thin=true]) -> Q, R, [p]

Compute the (pivoted) QR factorization of ``A`` such that either ``A = Q*R`` or ``A[:,p] = Q*R``. Also see ``qrfact``. The default is to compute a thin factorization. Note that ``R`` is not extended with zeros when the full ``Q`` is requested.
Expand Down
Loading