gsl_multifit_fdfsolver_lsmder memory usage

Peter Verveer
Sun Mar 16 19:54:00 GMT 2003

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:
<   gsl_matrix *q = state->q;
<   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:
<     gsl_matrix *q;
<   q = gsl_matrix_calloc (n, n);
<   if (q == 0)
<     {
<       GSL_ERROR ("failed to allocate space for q", GSL_ENOMEM);
<     }
<   state->q = q;
<   gsl_matrix_free (state->q);

In lmiterate replace compute_qtf() and gsl_linalg_QRPT_decomp2():
<   gsl_matrix *q = state->q;
<   compute_qtf (q, f, qtf);
>   gsl_vector_memcpy (qtf, f);
>   gsl_linalg_QR_QTvec (r, tau, qtf);
<         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():
< 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);
<     }
< }

More information about the Gsl-discuss mailing list