# Pivoted Cholesky vs Modified Cholesky

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.

• It is well-known that the normal equations are ill-conditioned and so are not typically used to solve least squares problems. The standard algorithms are QR and SVD decompositions. Why did you not consider those algorithms? – Bill Greene May 29 '19 at 16:41
• Because $J$ has a million rows QR or SVD will be impractical. – Brian Borchers May 29 '19 at 17:02
• Speed and memory. I do about 50 iterations to solve the nonlinear problem, so I need to solve 50 least squares problem with a matrix of about 10 million rows by 15,000 columns. To use QR/SVD, I need to store the matrix in memory, which would require about 1200 GB (which I don't have). I know there are fancier routines like TSQR but I don't have a fast implementation readily available. With normal eqs I don't need to store the full matrix, I can update $J^T J$ with blocks of rows at a time. – vibe May 29 '19 at 17:03
• Also, despite the ill-conditioning, my normal equations approach (with Modified Cholesky) does converge to a reasonable solution. I am just trying to speed it up by switching to a robust Level 3 BLAS algorithm. – vibe May 29 '19 at 17:06

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

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

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

• That is a good point about Levenberg-Marquardt. I have been experimenting with LM vs dogleg (dogleg requires calculating the Gauss-Newton solution). Perhaps I should restrict myself to LM only and then use the Level 3 Cholesky – vibe May 29 '19 at 17:15
• You'll find that LAPACK's Cholesky factorization with LM will do the job well. I'd encourage you to use an optimized BLAS/LAPACK library (like Intel's MKL or OpenBLAS) to get good performance. – Brian Borchers May 29 '19 at 17:32