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