Pivoted Cholesky vs Modified Cholesky

Question

I am solving nonlinear least squares problems with the normal equations approach, so on each iteration, I need to solve: $$ J^T J \delta = -J^T f $$ for the step $\delta$, where $J$ is a large (millions of rows and about $15,000$ columns), dense, Jacobian matrix, and $f$ is the residual vector. In my problem, $J$ is often ill-conditioned, and so the normal equations matrix $J^T J$ is even worse, and sometimes the standard Cholesky decomposition fails due to this during the iteration.

I am trying to find a robust method of solving this linear system and continuing the iteration, even if $J$ is badly conditioned. I have found two possibilities:

1. Pivoted Cholesky (description here), which is available in LAPACK using Level 3 BLAS calls, so it is presumably quite fast. This algorithm factors the matrix $A = J^T J$ as: $$ P^T A P = L L^T $$ and stops the decomposition when one of the diagonal entries of $A$ becomes smaller than some specified tolerance, so the algorithm is able to provide an estimate of rank($A$).

2. Modified Cholesky: this algorithm factors the matrix as

$$ P^T (A + E) P = LL^T $$ where $E$ is a small diagonal perturbation matrix designed to make $A+E$ positive definite and well-conditioned. This method appears to be specifically designed for nonlinear optimization problems, and is what I am currently using. The main issue is I cannot find an implementation of this method which is based on Level 3 BLAS. I have found a Level 2 BLAS based implementation, but it is significantly slower than the Level 3 Cholesky algorithm.

A reference for this is:

Schnabel, R. B., & Eskow, E. (1990). A new modified Cholesky factorization. SIAM Journal on Scientific and Statistical Computing, 11(6), 1136-1158.

My question is: what is a fast, robust approach to this problem? I want a reasonably accurate solution $\delta$ for each iteration, and I don't want the factorization to fail on any iteration, since it takes several days for the full nonlinear iteration to complete (I'm currently using the Level 2 Modified Cholesky approach and am looking to speed it up). I am tempted to switch to the Pivoted Cholesky approach, but I don't know if it is as robust as Modified Cholesky. Anyone with experience with these problems, I would be grateful for advice.

Answer 1

Neither of these approaches is recommended. Although pivoted Cholesky factorization can help with badly conditioned matrices, it ultimately won't help with a singular matrix. The modified Cholesky factorization could be used, but it's quite expensive computationally in comparison with an efficient implementation of the Cholesky factorization and isn't necessary to stabilize the algorithm.

There are several issues here that you should consider carefully.

1. Is the nonlinear least squares problem that you're solving badly conditioned in the sense that rather than having a unique optimal solution there is large set of solutions with the same (or nearly the same) optimal objective value? If there's nonuniqueness of the optimal solution to the nonlinear least squares problem then you will want to regularize the problem in some way.

This is a higher level issue than any difficulty in using the Cholesky factorization to factor $J^{T}J$.

2. The performance of the naive Gauss-Newton method that you're using can typically be improved substantially by using the Levenberg-Marquardt method, which modifies the equations to

$(J^{T}J+\lambda I) \delta = -J^{T}f$

where $\lambda$ is a parameter that is adjusted to keep the matrix effectively non-singular. You should be using the LM method.

3. How to solve the normal equations in each iteration. It is generally preferable to set this up as a linear least squares problem

$\min \left\| \left[ \begin{array}{c} J \\ \sqrt{\lambda} I \end{array} \right] \delta -\left[ \begin{array}{c} J^{T}f \\ 0 \end{array} \right] \right\|_{2} $

and then use a numerically stable technique involving the QR factorization or SVD to solve the least squares problem. However, since $J$ has a million or more rows and only 15,000 columns in your problem, you will want to use the Cholesky factorization of the normal equations for this problem.