4
$\begingroup$

I have an equation $Ax=b$ for which I need to solve for numerous $x$ matrices given $b$. Both $x$ and $b$ are nx1 matrices. Unfortunately, $A$ is a 32x32 matrix and inversion gives highly unstable values. I am wondering what the appropriate algorithm/technique is for obtaining a stable $A$ matrix.

My searches for an alternative calculation approach or way to restore stability to obtain a solution have so far led me nowhere so I would appreciate any input.

Thanks,

Jason

$\endgroup$
  • 4
    $\begingroup$ The $A$ matrix comes from your problem, so I can't help you on how to make that. As far as solving the system, the fully pivoted LU decomposition is numerically stable, but you will still lose precision if $A$ has a large condition number. What have you tried? $\endgroup$ – Tyler Olsen May 17 '15 at 5:51
  • 2
    $\begingroup$ Since you have already tagged this as an inverse problem, presumably you also know the standard regularization techniques (truncated singular value decomposition and Tikhonov regularization) to obtain a stable inversion. Are these not working for you? $\endgroup$ – Christian Clason May 17 '15 at 9:27
  • $\begingroup$ Matrix A is roughly known in this situation so not a problem. To get x, I have only tried LU decomposition to obtain the inverse and Cholesky decomposition to avoid calculating the inverse. I am using the Apache math commons libraries for this. I am not at all familiar with the techniques and am just doing a bunch of reading to try to find an approach. I am not familiar with Tikhonov regularization or truncated singular value decomposition so I will have a read on these to see if I make use of them. Thanks for the ideas so far. $\endgroup$ – Jason K. May 17 '15 at 14:56
  • $\begingroup$ Right, then I assumed incorrectly. Tikhonov and TSVD are exactly what you need here; I'll give a brief description below. $\endgroup$ – Christian Clason May 17 '15 at 15:21
13
$\begingroup$

If the solution of $Ax=b$ is unstable, the matrix is very ill-conditioned (i.e., has a very large condition number), and (paraphrasing Lanczos) no amount of mathematical trickery can make it stable. The best you can hope for is to solve a different problem that is a) stable and b) gives you a solution that is sufficiently close; this is called regularization.

For linear ill-conditioned problems, there are two classical approaches:

  1. Tikhonov regularization replaces $Ax=b$ by the stabilized least-squares problem $$ \min_x \|Ax-b\|^2 + \alpha \|x\|^2$$ for some regularization parameter $\alpha>0$. The minimizer can of course be computed by solving the normal equations $$ (A^TA+\alpha I)x_\alpha = A^Tb.$$

  2. Truncated singular value decomposition computes the singular value decomposition $$ U \Sigma V^T = A,$$ where the columns of $U$ and $V$ contain the left and right singular vectors, respectively, and $\Sigma$ is a diagonal matrix containing the singular values (usually in order of descending magnitude). Since $U$ and $V$ are unitary, the inverse of $A$ (if it exists) is given by $A^{-1}=V\Sigma^{-1} U^T$. Ill-conditionedness of $A$ manifests in the existence of very small singular values, such that the corresponding entry in $\Sigma^{-1}$ would be very large, amplifying small perturbations (e.g., due to finite numerical precision). As the name indicates, stability is restored by ignoring these small entries, i.e., replacing $\Sigma^{-1}$ by $\Sigma_{\alpha}^{-1}$ where $$[\Sigma_{\alpha}^{-1}]_{ii} = \begin{cases} \Sigma_{ii}^{-1} & \text{if } |\Sigma_{ii}| > \alpha \\ 0 &\text{else} \end{cases}$$ and setting $$x_\alpha = V\Sigma_{\alpha}^{-1} U^T b.$$

In both cases, you need to choose $\alpha$ specifically for your problem to get good results. You could start by computing the SVD of a representative $A$ and looking at the singular values to see if there's a clear threshold.

$\endgroup$
  • $\begingroup$ I am off to try these approaches. This is what I needed. Thanks. $\endgroup$ – Jason K. May 17 '15 at 15:59

Your Answer

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

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