Skip to content
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.

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

Add sparse matrix support #4

Open
ivan-pi opened this issue Feb 11, 2022 · 1 comment
Open

Add sparse matrix support #4

ivan-pi opened this issue Feb 11, 2022 · 1 comment

Comments

@ivan-pi
Copy link
Owner

ivan-pi commented Feb 11, 2022

The biggest advantage of libdogleg is it supports problems with sparse Jacobians. This can result in huge savings for large least squares problems.

Generally speaking, we just need to recover the sparse storage buffers of the compressed-column sparse (CCS) matrix of type cholmod_sparse and fill them with our Jacobian values. In C, the unpacking is done as follows:

static void optimizerCallback(const double*   p,
                              double*         x,
                              cholmod_sparse* Jt,
                              void*           cookie __attribute__ ((unused)) )
{

  // These are convenient so that I only apply the casts once
  int*    Jrowptr = (int*)Jt->p;
  int*    Jcolidx = (int*)Jt->i;
  double* Jval    = (double*)Jt->x;

Notably, the Jacobian matrix is referenced as J^T. The array Jrowptr is of size Nmeas + 1 (recall the output vector is of size Nmeas). The arrays Jcolidx and Jval are of size Jnnz, which is the number of non-zero values in the sparse matrix.

The number of non-zeros must be given explicitly to the libdogleg optimizer:

    optimum = dogleg_optimize2(p, Nstate, Nmeasurements, Jnnz,
                               &optimizerCallback, NULL,
                               &dogleg_parameters, NULL);

A full description of the CHOLMOD sparse storage format is given in the user guide: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/master/CHOLMOD/Doc/CHOLMOD_UserGuide.pdf

@ivan-pi
Copy link
Owner Author

ivan-pi commented May 29, 2022

One more observation, in the Fortran interface we should add support for 1-based indexing.

Since CHOLMOD is a C library, it expects 0-based indexing. This should be handled as an (optional) argument in the call to the optimizer.

By default we'd assume Fortran indexing is used.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant