3
$\begingroup$

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.

$\endgroup$
  • $\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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.