Fixed-structure CHOLMOD with in-place `set_nonzeros!`

General Usage Performance
, ,
1

I have a large sparse matrix. Its nonzeros change over the course of the problem, but its structure doesn’t. I want to do repeated in-place Cholesky factorizations, re-using the symbolic factorization and modifying the numeric factorization in-place. I can do the first factorization just fine:

using SparseArrays
using LinearAlgebra
A = sprand(10, 10, 0.3)
B = A + A' + 10 * LinearAlgebra.I(10)
B_cholmod = SparseArrays.CHOLMOD.Sparse(B)
sym = SparseArrays.CHOLMOD.symbolic(B_cholmod)
SparseArrays.CHOLMOD.cholesky!(sym, B_cholmod)

However, Julia’s interface does not allow me to write to the array of nonzero values of B_cholmod. So I’m implementing my own wrapper. Here’s what I have:

using SparseArrays
using LinearAlgebra
const VTypes = SparseArrays.CHOLMOD.VRealTypes
const ITypes = SparseArrays.CHOLMOD.ITypes
struct FixedStructureCHOLMOD{T, I} <: AbstractSparseMatrix{T, I}
    _mat::SparseArrays.CHOLMOD.Sparse{T, I}
    _values::Vector{T}
end

function FixedStructureCHOLMOD(A::SparseMatrixCSC{T, I}) where {T<:VTypes, I<:ITypes}
    cholmod_mat = SparseArrays.CHOLMOD.Sparse(A, -1) # -1 for hermitian
    chol = unsafe_load(cholmod_mat.ptr)
    x_ptr = chol.x
    # max capacity may not be the same as the number of nonzeros:
    # instead of chol.nzmax, we want chol.p[chol.ncol], in C indexing.
    # (what about packed matrices, though?)
    nonzero_count = unsafe_load(Ptr{I}(chol.p), chol.ncol + 1)
    values = unsafe_wrap(Vector{T}, Ptr{T}(x_ptr), nonzero_count)
    return FixedStructureCHOLMOD(cholmod_mat, values)
end


"""
    set_values!(mat::FixedStructureCHOLMOD, new_vals::AbstractVector{Float64})
In-place update of the numeric values in the CHOLMOD matrix.
"""
function set_values!(mat::FixedStructureCHOLMOD{T, I}, new_vals::AbstractVector{T}) where {T<:VTypes, I<:ITypes}
    @assert length(new_vals) == length(mat._values) "New values must be same length as allocated storage!"
    copyto!(mat._values, new_vals)
    return mat
end

function symbolic_factor(mat::FixedStructureCHOLMOD{T, I}) where {T<:VTypes, I<:ITypes}
    return SparseArrays.CHOLMOD.symbolic(mat._mat)
end

function numeric_factor!(F::SparseArrays.CHOLMOD.Factor{T, I}, mat::FixedStructureCHOLMOD{T, I}) where {T<:VTypes, I<:ITypes}
    SparseArrays.CHOLMOD.factorize!(mat._mat, F)
    return
end

My question: should this struct be mutable, with a proper finalizer, or is it memory-safe as is? edit: I’m asking this in the context described above. So constant sparse structure, one call to symbolic_factor, then many calls to numeric_refactor!.

Useful links: SuiteSparse’s cholmod.h, CHOLMOD’s User Guide, and Julia’s cholmod.jl.

Sample test code 
const K = 10
A = sprand(10, 10, 0.3)
B = A + A' + K * LinearAlgebra.I(10)
test_struct = FixedStructureCHOLMOD(B)
sym = symbolic_factor(test_struct)
numeric_factor!(sym, test_struct)
b = rand(10)
x = sym \ b
@assert isapprox(B * x, b, atol=1e-6)

# Test the in-place modification of values.
A2 = deepcopy(A)
nonzeros(A2) = rand(Float64, size(nonzeros(A2)))
B2 = A2 + A2' +  K * LinearAlgebra.I(10) # same structure
set_values!(test_struct, SparseArrays.nonzeros(B2))
numeric_factor!(sym, test_struct)
x2 = sym \ b
@assert isapprox(B2 * x2, b, atol=1e-6)
2

Why not just do F = cholesky(B) the first time and then do cholesky!(F, B) thereafter?

This is the documented way to re-use the symbolic factorization, and then you can just modify B and F in-place instead of messing with low-level CHOLMOD data structures.

3

I expect that the conversion to Sparse here is allocating and has a negative impact on performance. I wrote up a quick-and-dirty script designed to test that assumption:

using LinearAlgebra
using SparseArrays
using SuiteSparse
function testSolve(F, A)
    b = rand(size(A, 1))
    x = F \ b
    @assert isapprox(A * x, b, atol = 1e-6)
end

function choleskyAllocations(A1, A2)
    F = cholesky(A1)
    testSolve(F, A1)
    result = @allocated cholesky!(F, A2)
    testSolve(F, A2)
    return result
end

for N in (100, 1000, 10000)
    density = 0.001
    B = sprand(N, N, density)
    A1 = B + B' + 2*N/100 * LinearAlgebra.I(N)
    copyto!(nonzeros(B), rand(size(nonzeros(B), 1)))
    A2 = B + B' + N/100 * LinearAlgebra.I(N)
    println("N = $N:")
    println("   cholesky! memory usage: $(choleskyAllocations(A1, A2))")
    println("   Sparse memory usage: $(@allocated SuiteSparse.CHOLMOD.Sparse(A2))")
end

Output:

N = 100:
   cholesky! memory usage: 5936
   Sparse memory usage: 4640
N = 1000:
   cholesky! memory usage: 113152
   Sparse memory usage: 57024
N = 10000:
   cholesky! memory usage: 250266272
   Sparse memory usage: 3506880

Huh. So for large matrices, converting to Sparse has a pretty insignificant memory footprint compared to the library call itself, from half the memory footprint at N = 1,000 to around 1% at N = 10,000.

4

Yes, I think you are right — looks like it makes a copy of the data in the constructor here.

It’s certainly true that the sparse Cholesky factor typically requires a lot more memory than the matrix, due to fill-in. However, be aware that this depends strongly on the sparsity pattern, so using sprand may be misleading if your actual matrix is more structured.

(I’m also not sure if Julia is tracking all of the allocations that happen in the external library here?)