Any example of `cholmod_updown_solve()` (Updating in CHOLMOD)?

Question

I'm working on a project using CHOLMOD in C++ to do Cholesky factorization updating. The only reference I can find is the User Guide. And cholmod_updown_solve() seems to be the right function for me. But there is no example on this function and I can't get a correct result.

Below is the essential part of my code. According to the user guide, the only thing to notice seems to be sorting the update matrix C in advance. But that doesn't work. Could anyone tell me if there's any important step I missed?

Also, I'm confused because the user guide says that the solution phi will be given "in the permuted ordering, not your original ordering". How could I restore the order without knowing the permutation matrix? (L->Perm doesn't work.)

// The original system is At*A*phi = At*b
// Updates: C*Ct is added to At*A (Here C is a sparse column vector)

size_t n = m_pMesh->numVertices();
double w = 1e3;

cholmod_sparse *C;
cholmod_triplet *C_coefficients;
cholmod_dense *Delta_Atb;

cholmod_common common;
cholmod_common *cm = &common;
cholmod_start(cm);

C_coefficients = cholmod_allocate_triplet(n, 1, 2, 0, CHOLMOD_REAL, cm);
Delta_Atb = cholmod_zeros(n, 1, CHOLMOD_REAL, cm);

// updates: two more constraints
CViewerVertex *pNew = stroke_ends.start;
CViewerVertex *qNew = stroke_ends.end;
cholmodEntry(C_coefficients, pNew->sid(), pNew->sid(), w, cm);
cholmodEntry(C_coefficients, qNew->sid(), qNew->sid(), w, cm);

// change to At*b is Delta_Atb
((double*)Delta_Atb->x)[pNew->sid()] = w * w;

C = cholmod_triplet_to_sparse(C_coefficients, C_coefficients->nnz, cm);
cholmod_sort(C, cm);

// phi is the given solution to the original system At*A*phi = At*b
// L is the Cholesky factor to modify
// Both phi and L should be overwritten here
cholmod_updown_solve(1, C, L, phi, Delta_Atb, cm);

cholmod_free_sparse(&C, cm);
cholmod_free_triplet(&C_coefficients, cm);
cholmod_free_dense(&Delta_Atb, cm);

cholmod_finish(cm);

Answer 1

I also had several issues in getting a correct implementation using cholmod_updown_solve(). First of all, it might help understanding what the advantage of using this function is. A typical implementation using cholmod_updown() will look like this:

cholmod_updown(update, C, L, cm);
cholmod_dense *x = cholmod_solve(CHOLMOD_A, L, b, cm);

Running cholmod_solve() for the system CHOLMOD_A will run both the forward and backward subsitution algorithms, each requiring 2nnz(L) float point operations.

Alternatively, you can compute an initial solution to the Ly=Pb system as follows and reuse it later:

cholmod_dense *Pb = cholmod_solve(CHOLMOD_P, L, b, cm);
cholmod_dense *y = cholmod_solve(CHOLMOD_L, L, Pb, cm);
cholmod_dense *zeros = cholmod_zeros(L->n, 1, CHOLMOD_REAL, cm);
...
cholmod_updown_solve(update, C, L, y, zeros, cm);
cholmod_dense *Px = cholmod_solve(CHOLMOD_DLt, L, y, cm);
cholmod_dense *x = cholmod_solve(CHOLMOD_Pt, L, Px, cm);

This will be faster as cholmod_solve() for the system CHOLMOD_DLt will only run the backward substitution algorithm so if your update changes a small number of the columns of L you could get a speed improvement of almost 2x.

You have to be careful though, as cholmod_updown_solve() will convert any supernodal or LL' factorization into a simplicial LDL' factorization so before you initialize a solution for Ly=Pb you need to make sure you are working with an LDL' factorization. This can be achieved by running:

if (L->is_super || L->is_ll) cholmod_change_factor(CHOLMOD_REAL, 0, 0, 0, 0, L, cm);

before calling cholmod_solve() to initialize a solution to the Ly=Pb system.

As for constructing the update sparse matrix C to your matrix A=P'LDL'P, as explained in the source code, you have to actually provide cholmod_updown() or cholmod_updown_solve() with the sparse matrix PC which can be generated with the command:

cholmod sparse *PC = cholmod_submatrix(C, L->Perm, L->n, NULL, -1, 1, 1, cm);

Alternatively you can build an inverse permutation array:

int *IPerm = (int *)malloc(sizeof(int) * L->n);
for (int i = 0; i < L->n; i++) IPerm[((int *)L->Perm)[i]] = i;

And then update the row indexes of C using IPerm:

for (int i = 0; i < C->nzmax; i++) ((int *)C->i)[i] = IPerm[((int *)C->i)[i]];