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.
