Re: gsl_multifit_fdfsolver_lsmder memory usage

[Peter Verveer <verveer at embl-heidelberg dot de>]

Hi Brian,

>  > I recently started converting some numerical code to GSL and
>  > replaced a Levenberg-Marquardt routine which was a C translation of
>  > lmdif from MINPACK.  However, I found that the GSL solver
>  > gsl_multifit_fdfsolver_lsmder uses a lot of memory, since it
>  > appears to allocate a matrix of m times m, where m is the number of
>  > functions. This causes me many problems since I work with large
>  > data sets (10,000-100.000 data points is not unusual) and obviously
>  > memory requirements quickly get out of hand then. I had a look at
>  > the original fortran MINPACK which does not seem to require
>  > this. Is there a particular reason gsl_multifit_fdfsolver_lsmder
>  > needs so much memory, and could this eventually be changed?
>
> When I implemented the algorithm I believe I stored the full Q matrix
> from the QR decomposition because it was easier to implement that way.
>
> Maybe you could look at what is needed to avoid storing the full Q
> matrix, as in MINPACK. I don't think it's too complicated now that the
> basic algorithm is there.

You are right. It turns out that the changes are very minor. I replaced the QR 
decomposition by its in-place version, and used the resulting encoded QR 
matrix to calculate qtf. Then the q matrix in the solver can be eliminated, 
saving huge amounts of memory in my case.

Below are the patches that I applied. They are quite small, so maybe they can 
be included in future versions of GSL? I hope they are correct and don't 
introduce any problems I did not discover.

I should note that I do not actually unpack the QR matrix to get R. The packed 
matrix QR can be used directly instead of R since its upper triangular matrix 
is equal to that of R (at least I think...). This seems to work since the 
rest of the code does not access the lower triangular matrix of R. If that 
seems a problem, R should be extracted, or the lower triangular part of QR 
set to zero to get R.

Thanks for the hint, the code works now fine for me.

Cheers, Peter
-------------------------------------------------------------------------------------------------------------------------------
Here are the patches:

In lmset.c just replace the QR decomposition with its in-place version:
6d5
<   gsl_matrix *q = state->q;
41c40,41
<   gsl_linalg_QRPT_decomp2 (J, q, r, tau, perm, &signum, work1);
---
>   gsl_matrix_memcpy (r, J);
>   gsl_linalg_QRPT_decomp (r, tau, perm, &signum, work1);

In lmder.c, remove the allocation and deallocation of q:
43d42
<     gsl_matrix *q;
84,92d82
<   q = gsl_matrix_calloc (n, n);
< 
<   if (q == 0)
<     {
<       GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM);
<     }
< 
<   state->q = q;
< 
389d378
<   gsl_matrix_free (state->q);

In lmiterate replace compute_qtf() and gsl_linalg_QRPT_decomp2():
6d5
<   gsl_matrix *q = state->q;
31,32c30,31
< 
<   compute_qtf (q, f, qtf);
---
>   gsl_vector_memcpy (qtf, f);
>   gsl_linalg_QR_QTvec (r, tau, qtf);
187c186,187
<         gsl_linalg_QRPT_decomp2 (J, q, r, tau, perm, &signum, work1);
---
> 	gsl_matrix_memcpy (r, J);
> 	gsl_linalg_QRPT_decomp (r, tau, perm, &signum, work1);

In lmutils, remove the definition of compute_qtf():
100,113d99
< static void
< compute_qtf (const gsl_matrix * q, const gsl_vector * f, gsl_vector * qtf)
< {
<   size_t i, j, N = f->size;
< 
<   for (j = 0; j < N; j++)
<     {
<       double sum = 0;
<       for (i = 0; i < N; i++)
< 	sum += gsl_matrix_get (q, i, j) * gsl_vector_get (f, i);
< 
<       gsl_vector_set (qtf, j, sum);
<     }
< }