RE: gsl_multifit

[Mikael Adlers <mikael at mathcore dot com>]

Hi,
the algorithm in gsl_multifit_linear uses the singular value
decomposition (SVD). This routine gives the best stability and
accuracy when solving a least squares problem. However, in most
cases it is like shooting a fly with a cannon. The complexity 
of the SVD is about 4mn^2+11/3n^3 compared to 1/2mn^2 +n^3/6 
using the normal equations + Cholesky factorization (or twice the 
amount using QR factorization instead (better stability than
the normal equations)). 

Further, gsl_multifit_linear computes a covariance matrix and
some other nifty stuff eating flops away.

The algorithm in Matlab (and Octave?) uses the QR factorization 
(together with one step of iterative refinement I believe) to 
compute the solution.

You could try computing the solution using the following 
functions instead and see if you get any speed improvement.
gsl_linalg_QR_decomp
gsl_linalg_QR_lssolve

/Mikael Adlers




------------------------------------------------------------------ 
 Mikael Adlers, Ph.D.          email: mikael@mathcore.com 
 MathCore AB                   phone: +4613 32 85 07 
 Wallenbergsgata 4             fax:         21 27 01
 SE-583 35 Linköping, Sweden   http://www.mathcore.com




> -----Original Message-----
> From: Kai Trukenmueller [mailto:trukenm@ag2.mechanik.tu-darmstadt.de] 
> Sent: den 16 oktober 2001 22:51
> To: gsl-discuss@sources.redhat.com
> Subject: gsl_multifit
> 
> 
> Hi,
>
> The problem A*x=b with A not square leeds to
> x = (A^T *A)^-1 A^T *b
> at least analytically.
> That should be equivalent to a multidimensional fitting, as in
> `gsl_multifit_linear'.
> I implemented an example and it works good.
>
> Than I compared it to a simple octave (a free matlab clone) routine,
> using
> x=A\b;
>
> The results are the same (as it should be), but it surprised me a lot,
> that the octave-script works _much_ faster for high dimensions than the
> compiled C programm. It seems, that these algorithems are better.
>
> Is that a known problem, or did I use the wrong routines.
> Is there any better way for solving A x = b?
> -- 
> :wq `Kai Trukenmueller'


> -----Original Message-----
> From: Kai Trukenmueller [mailto:trukenm@ag2.mechanik.tu-darmstadt.de] 
> Sent: den 17 oktober 2001 01:36
> To: gsl-discuss@sources.redhat.com
> Subject: Re: gsl_multifit
> 
> 
> Hi,
> 
> I benchmarked the programs on system with Athlon 1GHz, 640MB 
> running Linux-2.4.10 with the following results:
> 
> Dimensions	(120,80)  (300,200)  (600,400)  (900,600)  (1200,800)
> Octave-Script	0.060s    0.410s      3.390s     12.260s    39.790s
> Compiled gsl    0.110s    2.240s      31.310s    2m7.590s   seg.fault
> 
> The gsl routine not just turns out to be much slower, it also seems to
> be `instable' for higher orders (the results do not convege -> inft..)
> Maybe sth. is wrong in my code. For low-orders (~<500) both 
> results are
> equivalent.
> 
> I'm not using the atlas-blas but gslcblas.
> 
> Code attached;
> The octave-script `sinft.m' should be executable, and can be started
> directely from the shell (iff /usr/bin/octave exists).
> The matrix dimensions must be edited in the source code (in the c
> programm, they are arguments).
> 
> 
> -- 
> :wq `Kai Trukenmueller'
>