2
$\begingroup$

I have read that matrix regularisation can improve the stability of LU or Cholesky decomposition of ill conditioned problems. The idea is to add a value to the diagonals of a matrix:
$B=A+cI$
In the above mentioned case $A$ is the initial matrix, $c$ is the constant, $I$ is the identity matrix and $B$ is the regulized matrix.
Two questions arise here from my side:

  1. How to select the value of $c$?
  2. If I want to solve $Ax=b$ but have regulized the matrix leading to $Bx=b$, how should I now solve for $x$ as my initial problem is not the same anymore?
Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

matrix regularisation can improve the stability of LU or Cholesky decomposition of ill conditioned problems.

Not really, at least in the way the word "stability" is typically used in the numerical literature. Cholesky is always backward stable, and LU (with pivoting) is backward stable when the entries of $U$ do not grow too much, a property that usually holds even for ill-conditioned matrices.

Matrix regularization has a different goal. In short, it helps you produce a prettier "pseudo-solution" with a small residual but also small norm. This is useful, for instance, if you have reason to believe that the original problem has a small-norm solution, and noise (combined with the ill-conditioning) has perturbed the problem so much that this solution is impossible to recover. A full treatment would fill a book on inverse problems, and I am not a domain expert so I wouldn't be able to give it.

If I want to solve ๐ด๐‘ฅ=๐‘ but have regulized the matrix leading to ๐ต๐‘ฅ=๐‘, how should I now solve for ๐‘ฅ as my initial problem is not the same anymore?

Typically, you can't: noise (combined with the ill-conditioning) has perturbed your problem. You can compute the exact solution to $Ax=b$ with the noisy values of $A$, $b$ that you have on your machine, but typically it will be an unusable mess with huge, imbalanced entries.

How to select the value of ๐‘?

It's impossible to tell based on math alone, since the goal here is impossible to define without referencing the original application. There are several common heuristics, such as looking for a 'bend' in the plot of the solution norm $\|(A+cI)^{-1}b\|$ vs. $c$, or estimating the amount of noise in your data. Machine learning people love grid searches and cross-validation, which are a way to do so.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for the detailed answer. Basically I was looking for some pre-processing step of a sparse matrix coming from FEM to avoid pivoting in direct solution methods. $\endgroup$
    – vydesaster
    Jul 11, 2022 at 16:46
  • $\begingroup$ So you have a well-conditioned problem, and you are adding a perturbation to avoid pivoting? Interesting set-up, I had never heard about this use case. In this case, you are indeed perturbing to improve stability, then. And you might be able to recover the original solution, at least with low accuracy, using extrapolation. $\endgroup$
    – Federico Poloni
    Jul 11, 2022 at 16:50
  • $\begingroup$ Do you know any other pre-processing steps which enables me to avoid pivoting during factorization? $\endgroup$
    – vydesaster
    Jul 11, 2022 at 17:03
  • 2
    $\begingroup$ If it's a SPD system that is eligible for sparse Cholesky factorization, you really shouldn't have to do this. If it's indefinite or nonsymmetric, then it's almost inevitable that some form of pivoting will be required. You might consider a new question that's more focused on your problem (ie that clarifies you're interested in sparse matrices, and talks about the origins/spectral properties of your matrix). $\endgroup$
    – rchilton1980
    Jul 11, 2022 at 17:08
  • 1
    $\begingroup$ @vydesaster you may want to look at threshold pivoting (which allows smaller pivots to be selected if they're good for another metric, e.g., fill-in) and static-pivoting (which doesn't pivot during the factorization and instead just adds a small factor to small diagonal elements). Both approaches usually have larger errors that conventional pivoting, so iterative refinement or another preconditioned-method (e.g. gmres) is needed if you want a full precision solution. $\endgroup$
    – Neil Lindquist
    Jul 11, 2022 at 22:25

Your Answer

By clicking โ€œPost Your Answerโ€, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.