Question

Levenberg-Marquardt (Damped Least Squares)

Recall the least squares problem:

Given a set of data points \((x_i, y_i)\) where \(x_i\)‘s are independent variables (in \(\mathbb{R}^n\) and the \(y_i\)‘s are response variables (in \(\mathbb{R}\)), find the parameter values of \(\beta\) for the model \(f(x;\beta)\) so that

\[S(\beta) = \sum\limits_{i=1}^m \left(y_i - f(x_i;\beta)\right)^2\]

is minimized.

If we were to use Newton’s method, our update step would look like:

\[\beta_{k+1} = \beta_k - H^{-1}\nabla S(\beta_k)\]

Gradient descent, on the other hand, would yield:

\[\beta_{k+1} = \beta_k - \gamma\nabla S(\beta_k)\]

Levenberg-Marquardt adaptively switches between Newton’s method and gradient descent.

\[\beta_{k+1} = \beta_k - (H + \lambda I)^{-1}\nabla S(\beta_k)\]

When \(\lambda\) is small, the update is essentially Newton-Gauss, while for \(\lambda\) large, the update is gradient descent.

Newton-Krylov

The notion of a Krylov space comes from the Cayley-Hamilton theorem (CH). CH states that a matrix \(A\) satisfies its characteristic polynomial. A direct corollary is that \(A^{-1}\) may be written as a linear combination of powers of the matrix (where the highest power is \(n-1\)).

The Krylov space of order \(r\) generated by an \(n\times n\) matrix \(A\) and an \(n\)-dimensional vector \(b\) is given by:

\[\mathcal{K}_r(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{r-1}b \}\]

Thes are actually the subspaces spanned by the conjugate vectors we mentioned in Newton-CG, so, technically speaking, Newton-CG is a Krylov method.

Now, the scipy.optimize newton-krylov solver is what is known as a ‘Jacobian Free Newton Krylov’. It is a very efficient algorithm for solving large \(n\times n\) non-linear systems. We won’t go into detail of the algorithm’s steps, as this is really more applicable to problems in physics and non-linear dynamics.