>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. |
> 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. |
Free forum by Nabble - Free Resume Builder | Edit this page |