ivan-pi / libdogleg-f

Fortran bindings to libdogleg - a large-scale nonlinear least-squares optimization library

libdogleg-f

Fortran bindings to libdogleg - a large-scale nonlinear least-squares optimization library.

Currently only the dense optimizer calls are supported.

Installation

To use libdogleg-f you will need:

• a Fortran compiler supporting Fortran 2018 (GCC 5 or newer or Intel Fortran)
• the Fortran package manager (fpm)
• libdogleg
• CHOLMOD and SuiteSparse
• LAPACK

Include libdogleg-f in your fpm package manifest using

[dependencies]
libdogleg-f.git = "https://github.com/ivan-pi/libdogleg-f"

License and copyright

Copyright (C) 2022 Ivan Pribec

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with this program (see COPYING.LESSER and COPYING). If not, see https://www.gnu.org/licenses/.

libdogleg-f's Issues

Add support for different versions

The latest version of libdogleg slightly expanded the API with the parameter type. The Fortran wrapper here will only work with the latest version (0.15), and not with the older version (0.14) available in older Debian/Ubuntu distros.

We can support both with some fpp preprocessor statements.

Build systems

Since the library is relatively small, support for additional build systems including Make and CMake should be added.

Documentation following Diataxis

The documentation should follow the Diataxis framework (https://diataxis.fr/).

Here are a few possible sub-topics:

reference

• description
  • which method is used
  • which factorizations are used
• additional links
  • original library
  • method used
  • data structures
• other nonlinear least squares solvers

tutorial

• learn to use libdogleg-f via a curve-fitting example
• determining the equivalent point of a titration
  • (not the perfect example, because we need to find the derivative, and minimum)
• camera calibration
• computer vision, fitting a rectangle to an image
• variable scaling for ill-conditioned problems

explanation

• stating the problem dimensions (explicit vs implicit)
  • generic functions (retrieve size from dummy arrays)
  • retrieve sizes, via host association
  • global variables
• sparse vs dense jacobians
• libdogleg-f internals and the adaptor pattern

how to's

• how to re-use an existing callback in libdogleg
• installing libdogleg dependencies
• converting from a minpack callback
• how to retrieve a solver context
• controlling the solver parameters

Add support for retrieval and handling of solver context

The return Context is handled as an opaque C pointer. To actually access it's data members we need to write some C functions.

This is only for advanced users, who may want to analyze the Jacobian factorization.

add topics

I suggest adding the topics nonlinear-least-square, numerical-optimization. I added the repo to my categorized list.

Refactor examples

The subroutines to print the state vector before and after, and the subroutine to display the dogleg parameters should go to a dedicated module, reused by all the examples.

Add sparse matrix support

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