# Numerically stable approach for calculating x in Ax=b

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

• 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? – Tyler Olsen May 17 '15 at 5:51
• 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? – Christian Clason May 17 '15 at 9:27
• 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. – Jason K. May 17 '15 at 14:56
• Right, then I assumed incorrectly. Tikhonov and TSVD are exactly what you need here; I'll give a brief description below. – Christian Clason May 17 '15 at 15:21

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.
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.