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.

  • $\begingroup$ 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.) $\endgroup$ – Richard Zhang Dec 12 '18 at 16:43
  • $\begingroup$ It is positive semi-definite (eigenvalues are either 0 or positive) $\endgroup$ – vibe Dec 12 '18 at 17:03
  • $\begingroup$ 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$. $\endgroup$ – vibe Dec 12 '18 at 17:11
  • $\begingroup$ Another option is to use the symmetric square root of your positive semidefinite matrix. $\endgroup$ – Brian Borchers Dec 12 '18 at 17:49
  • 1
    $\begingroup$ 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. $\endgroup$ – Richard Zhang Dec 12 '18 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.