COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

### Weighted least squares with minpack?

• Follow

```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

```
 0

```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)

```
 0
Reply highegg (245) 9/19/2007 6:39:23 PM

```> 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

```
 0

```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

```
 0
Reply highegg (245) 9/19/2007 7:16:52 PM

```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.
```
 0
Reply NOSPAM3879 (120) 9/21/2007 1:48:56 AM

```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.

definition of residual was opposite in sign to the conventional one.

```
 0
Reply NOSPAM3879 (120) 9/21/2007 4:16:26 AM

5 Replies
30 Views

5/22/2013 8:34:52 PM