nls.lm

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

nls.lm

Mike meyer
>From my reading of the above cited document I get the impression that the algorithm
(algorithm 3.16, p27) can easily be adapted to handle the case m<n.
In this case the Jacobian Jf(x) is mxn and the matrix A(x)=Jf(x)'Jf(x) is nxn and rank deficient.

This seems to be a problem since the search direction h at iterate x is computed from
       (A(x)+mu*I)h=-J(x)'f(x)
where mu -> 0 and so the system becomes ill conditioned.

Why can we not get around this as follows: as soon as mu is below some threshold
we solve instead the seemingly worse A(x)h=-J(x)'f(x) which however is of the form

      J(x)'J(x)h=-J(x)'f(x)

and can be solved by solving J(x)h=-f(x) and the condition for this to be solvable
is that J(x) have full rank. In fact this is the Gauss-Newton step
and to switch to that in the case of small mu is exactly in keeping with the idea of the
LM-algorithm.

If h is _any_ solution we have h'F'(x)=-||J(x)h||²<=0 (document, p21, eq(3.10))
and is a descent direction if this is in fact < 0.
If it is equal to zero, then f(x)=J(x)h=0 and we have a global minimum at x.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: nls.lm

Berend Hasselman

> On 19 Oct 2016, at 19:12, Mike meyer <[hidden email]> wrote:
>
>> From my reading of the above cited document I get the impression that the algorithm
> (algorithm 3.16, p27) can easily be adapted to handle the case m<n.
> In this case the Jacobian Jf(x) is mxn and the matrix A(x)=Jf(x)'Jf(x) is nxn and rank deficient.
>
> This seems to be a problem since the search direction h at iterate x is computed from
>       (A(x)+mu*I)h=-J(x)'f(x)
> where mu -> 0 and so the system becomes ill conditioned.
>
> Why can we not get around this as follows: as soon as mu is below some threshold
> we solve instead the seemingly worse A(x)h=-J(x)'f(x) which however is of the form
>
>      J(x)'J(x)h=-J(x)'f(x)
>
> and can be solved by solving J(x)h=-f(x) and the condition for this to be solvable
> is that J(x) have full rank. In fact this is the Gauss-Newton step



Oh?

You are solving a system

Ax=b where A is a matrix mxn with m<n with full rank (m).

Let Aplus be the generalized inverse of A.
Let z <- runif(nrow(Aplus))

The general solution of that linear system is   Aplus %*% b + (diag(nrow(Aplus)) - Aplus %*% A) %*% z
where the minimum norm solution is  Aplus %*% b.

Since z is any random vector there are many solutions.
So which one is the Gauss-Newton step?
 
Berend

> and to switch to that in the case of small mu is exactly in keeping with the idea of the
> LM-algorithm.
>
> If h is _any_ solution we have h'F'(x)=-||J(x)h||²<=0 (document, p21, eq(3.10))
> and is a descent direction if this is in fact < 0.
> If it is equal to zero, then f(x)=J(x)h=0 and we have a global minimum at x.
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.