Numerically stable way to implement Cramer's rule analog

Question

Problem statement

Let $A$ be an $n\times n$ matrix and $b$ an $n$-dimensional vector. For $j\in \{1, \dots, n \}$, let $A_j$ be the matrix where we take $A$ and replace the $j^{\rm th}$ column with $b$. I want to find the vector $z = \frac{\tilde z}{\max_j |\tilde z_j|}$, where $\tilde z = (\det A_1, \dots, \det A_n)$.

If $A$ is invertible, then by Cramer's rule it follows that $\tilde z = (\det A) A^{-1}b$. So to compute $z$, I can use my favorite linear solver to solve the linear system $Ax=b$ and set $z = \frac{x}{\max_j |x_j|}$.

However, I am interested in the case when $A$ is not invertible. In this case, I cannot just solve the linear system. Instead, as far as I can tell, I have two options:

1. Compute all $n$ determinants explicitly to determine $\tilde z$, or
2. Let $U \Sigma V^T$ be the singular value decomposition of $A$. Then $\tilde z = V (\det \Sigma_1, \dots, \det \Sigma_n)$, where $\Sigma_j$ is the matrix where we take $\Sigma$ and replace the $j^{\rm th}$ column by the vector $U^T b$.

The second option is probably more computationally efficient since computing the determinant of $\Sigma_i$ is very easy since $\Sigma$ is diagonal.

Question 1

Is there another (better) way to find $z$?

Question 2 (my main question)

When I implement the above procedures in Python using Numpy, I get crazy results that change every time I run it. The reason for this (as far as I can tell) is that the whole procedure is very numerically unstable (please let me know if you know of a better reason why this is happening!). $A$ is not invertible, so up to numerical precision $\det A=0$. But furthermore, generally $\det A_j$ seems to also be very close to zero, albeit not actually zero in general. Thus, $\tilde z$ is filled with very very small numbers that are essentially zero up machine precision. However, to find $z$ I have to divide by $\max_j |\tilde z_j|$! So in the end, the vector $z$ should have perfectly reasonable not-essentially-zero entries, but numerically I have to get these reasonable numbers by dividing really small numbers by a really small number. Is there any way that I can get around this?

Misc notes

• A possible third way of approximately computing $\tilde z$ is to perturb $A$ by a very small but nonzero matrix $B$ such that $A+B$ is invertible even though $A$ is not invertible. But the problem with this is that:
• Even in the case that $A$ is invertible but close to being singular, I run into these numerical issues when calculating $z$.

Answer 1

A partial answer to your point 2, from a comment I wrote to a now-deleted answer: there is nothing wrong about $\det A_j$ being very small: determinants are notoriously poorly scaled. For instance the determinant of $0.8I$, where $I$ is a $500\times 500$ identity matrix, is $3.5\times 10^{-49}$; yet $0.8I$ is a perfectly reasonable matrix to work with.

So this is not the cause of your problem. And, to me, your SVD-based method looks good from the point of view of stability.

Can you tell us more on why your example fails? And I suggest to look at the condition numbers of the $A_j$ rather than their determinants, to understand if they really are singular up to numerical precision.

EDIT: on a second look, are you sure that your SVD method works? Have you tested it on an easy non-singular example? It seems strange that those column replacements can be done after orthogonal transformations.