Hi all,

I'm fumbling around trying to figure out how to use minpack (lmstr)
for a weighted non-linear least squares fit and I know it shouldn't be
as tough as I'm making it.  For weights, I'm using a full covariance
matrix, S, so in Bx = y form, the linear problem is:

  (K^T S^(-1) K) x = K^T S^(-1) y

where K is the Jacobian matrix.  My confusion is that lmstr (and most
similar non-linear least squares routines) all expect you to provide:

  m, n = rows, cols of K
  x = initial guess (overwritten with solution)
  y = data to fit
  fcn = function to evaluate F(x) and Jacobian matrix (dF/dx)

It isn't obvious to me how to weight the data or write the model
function 'fcn' when weights are used.

Does anyone have a basic fortran example showing how to set this up or
can point me to a text with practical examples?

Many thanks,
Dave

Confused Scientist napsal(a):
> Hi all,
>
> I'm fumbling around trying to figure out how to use minpack (lmstr)
> for a weighted non-linear least squares fit and I know it shouldn't be
> as tough as I'm making it.  For weights, I'm using a full covariance
> matrix, S, so in Bx = y form, the linear problem is:
>
>   (K^T S^(-1) K) x = K^T S^(-1) y
>
> where K is the Jacobian matrix.  My confusion is that lmstr (and most
> similar non-linear least squares routines) all expect you to provide:
>
>   m, n = rows, cols of K
>   x = initial guess (overwritten with solution)
>   y = data to fit
>   fcn = function to evaluate F(x) and Jacobian matrix (dF/dx)
>
> It isn't obvious to me how to weight the data or write the model
> function 'fcn' when weights are used.
>
> Does anyone have a basic fortran example showing how to set this up or
> can point me to a text with practical examples?
>
> Many thanks,
> Dave

A good solution is to factorize (cholesky) the covariance matrix as S
= L*L' and then,
instead of fitting f(x) to y, fit inv(L) * f(x) to  inv(L) * y. (' is
transpose, inv is inverse)
(i.e. additional TRSM after function or jacobian evaluation)

> A good solution is to factorize (cholesky) the covariance matrix as S
> = L*L' and then,
> instead of fitting f(x) to y, fit inv(L) * f(x) to  inv(L) * y. (' is
> transpose, inv is inverse)
> (i.e. additional TRSM after function or jacobian evaluation)

What does 'TRSM' mean?

If I understand what you suggest, I should pre-multiply y and fcn so
would pass the following to lmstr:

  x = first guess
  fcn = inv(L) * fcn_orig(x)
  y = inv(L) * y_orig

where fcn_orig is my original function and y_orig is my original
data.  NR gives a trivial Cholesky factorization routine so that
should be pretty easy.

Thanks, I'll see how that works... heck, I might even post my
attempt ;)
Dave

On Sep 19, 9:08 pm, Confused  Scientist <Confused.Scient...@gmail.com>
wrote:
> > A good solution is to factorize (cholesky) the covariance matrix as S
> > = L*L' and then,
> > instead of fitting f(x) to y, fit inv(L) * f(x) to  inv(L) * y. (' is
> > transpose, inv is inverse)
> > (i.e. additional TRSM after function or jacobian evaluation)
>
> What does 'TRSM' mean?


xTRSM are BLAS routines to multiply a matrix by the inversion of a
triangular matrix.

>
> If I understand what you suggest, I should pre-multiply y and fcn so
> would pass the following to lmstr:
>
>   x = first guess
>   fcn = inv(L) * fcn_orig(x)
>   y = inv(L) * y_orig
>
> where fcn_orig is my original function and y_orig is my original
> data.  NR gives a trivial Cholesky factorization routine so that
> should be pretty easy.

Yeah, that's it.
Nitpick: I think you do not give y and fcn separately to MINPACK
routines (lmder or lmstr), instead you provide fcn - y (the algorithms
simply minimize norm(fcn(x))).


> Thanks, I'll see how that works... heck, I might even post my
> attempt ;)
> Dave

Confused Scientist wrote:
> Hi all,
> 
> I'm fumbling around trying to figure out how to use minpack (lmstr)
> for a weighted non-linear least squares fit and I know it shouldn't be
> as tough as I'm making it.  For weights, I'm using a full covariance
> matrix, S, so in Bx = y form, the linear problem is:
> 
>   (K^T S^(-1) K) x = K^T S^(-1) y
> 
> where K is the Jacobian matrix.  My confusion is that lmstr (and most
> similar non-linear least squares routines) all expect you to provide:
> 
>   m, n = rows, cols of K
>   x = initial guess (overwritten with solution)
>   y = data to fit
>   fcn = function to evaluate F(x) and Jacobian matrix (dF/dx)
> 
> It isn't obvious to me how to weight the data or write the model
> function 'fcn' when weights are used.
> 
> Does anyone have a basic fortran example showing how to set this up or
> can point me to a text with practical examples?
> 
> Many thanks,
> Dave
> 
If I remember correctly, Minpack uses as input the residuals deltay = 
y_calculated - y_observed.  Simply replace each deltay(i) in the input 
array by deltay(i)/sigma(i), where sigma(i) is the standard deviation of 
measurement i, and you have weighted regression.

Charles Russell wrote:
> Confused Scientist wrote:
>> Hi all,
>>
>> I'm fumbling around trying to figure out how to use minpack (lmstr)
>> for a weighted non-linear least squares fit and I know it shouldn't be
>> as tough as I'm making it.  For weights, I'm using a full covariance
>> matrix, S, so in Bx = y form, the linear problem is:
>>
>>   (K^T S^(-1) K) x = K^T S^(-1) y
>>
>> where K is the Jacobian matrix.  My confusion is that lmstr (and most
>> similar non-linear least squares routines) all expect you to provide:
>>
>>   m, n = rows, cols of K
>>   x = initial guess (overwritten with solution)
>>   y = data to fit
>>   fcn = function to evaluate F(x) and Jacobian matrix (dF/dx)
>>
>> It isn't obvious to me how to weight the data or write the model
>> function 'fcn' when weights are used.
>>
>> Does anyone have a basic fortran example showing how to set this up or
>> can point me to a text with practical examples?
>>
>> Many thanks,
>> Dave
>>
> If I remember correctly, Minpack uses as input the residuals deltay = 
> y_calculated - y_observed.  Simply replace each deltay(i) in the input 
> array by deltay(i)/sigma(i), where sigma(i) is the standard deviation of 
> measurement i, and you have weighted regression.

On second reading, your weighting is more complicated.  Sorry.  Also, my 
definition of residual was opposite in sign to the conventional one.