Nonlinear least squares and regularization

Consider the nonlinear least-squares minimization of a vector of $$n$$ residuals $$\mathbf{f}$$ in $$p$$ parameters $$\mathbf{x}$$: $$\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2$$ This can be done with a standard Gauss-Newton or Levenberg-Marquardt approach, using the $$n \times p$$ Jacobian matrix $$J_{ij} = \frac{\partial f_i}{\partial x_j}$$ Now suppose I want to regularize the solution, so I want to solve $$\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2 + \lambda^2 x^T C x$$ for some symmetric $$p \times p$$ matrix $$C$$. It seems to me there are two ways to proceed.

Option 1

If $$C$$ happens to be SPD, then we can do a Cholesky decomposition $$C = U^T U$$ where $$U$$ is upper triangular. Then we get $$\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2 + \lambda^2 || U \mathbf{x}||^2$$ Here we can form the augmented residual vector and Jacobian: \begin{align} \tilde{\mathbf{f}} &= \pmatrix{ \mathbf{f} \\ \lambda U \mathbf{x}} \\ \tilde{J} &= \pmatrix{ J \\ \lambda U} \end{align} and the problem becomes $$\min_{\mathbf{x}} || \tilde{\mathbf{f}}(\mathbf{x}) ||^2$$ for which we can apply a standard Gauss-Newton or LM code to solve. The augmented system will have $$n + p$$ residuals and a $$(n+p) \times p$$ Jacobian matrix. For a large number of parameters $$p$$ (say several thousand), this could add significant cost to the nonlinear iteration.

Option 2

Alternatively, suppose $$C$$ is symmetric but not positive definite, then we can't perform a Cholesky factorization. In this case we can still proceed by forming the augmented system: \begin{align} \tilde{\mathbf{f}} &= \pmatrix{ \mathbf{f} \\ \lambda \sqrt{x^T C x}} \\ \tilde{J} &= \pmatrix{ J \\ \lambda \frac{x^T C}{x^T C x}} \end{align} In this case, we now have $$n + 1$$ residuals and the Jacobian is $$(n+1) \times p$$.

My question is whether Option 2 is superior to Option 1? Option 2 requires factoring a smaller Jacobian matrix which could result in significant speed improvements when $$p$$ is very large. However I wonder if Option 1 might converge faster to the solution, since more "information" is stored in the residual vector and Jacobian. Does anyone have any insights into this?

For a long time I used Option 1, but recently I encountered a problem where $$C$$ is symmetric but not positive definite, so I am forced to consider Option 2. Just wondering if Option 2 might have worse convergence properties.

EDIT: it seems there is a third option which may or may not be useful:

Option 3

The minimization problem can be put in the form:

$$\min_{\mathbf{x}} \tilde{\mathbf{f}}^T W \tilde{\mathbf{f}}$$ where $$W = \pmatrix{ I & 0 \\ 0 & C}$$ and $$\tilde{\mathbf{f}} = \pmatrix{ \mathbf{f} \\ \lambda x}$$ In this form, it may not be possible to use an off-the-shelf LM software due to the $$W$$ matrix, but certainly a specialized algorithm could be developed to solve this problem. Again I am not clear if this offers any advantage over Option 1 or 2.

• When you say that $C$ may not be positive definite, are you saying that it might be indefinite, or just positive semidefinite? (The former is much, much harder than the latter.) – Richard Zhang Dec 12 '18 at 16:43
• It is positive semi-definite (eigenvalues are either 0 or positive) – vibe Dec 12 '18 at 17:03
• So for my problem, the choice is between Option 2 and 3 - I cannot use Option 1 since the Cholesky factorization doesn't exist for my $C$. – vibe Dec 12 '18 at 17:11
• Another option is to use the symmetric square root of your positive semidefinite matrix. – Brian Borchers Dec 12 '18 at 17:49
• The Cholesky factorization does exist, but square-root based routines can run into stability issues. Instead, you can also just compute the LDL decomposition and round all negative elements in D to zero. – Richard Zhang Dec 12 '18 at 18:24