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