# Numerically stable way to implement Cramer's rule analog

## 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$$.
• Re "I get crazy results that change every time I run it": If this happens with identical input data, there is a high likelihood that an out-of-bounds access or uninitialized variable occurs somewhere in the code. While results may consists entirely of noise in extremely ill-conditioned cases, floating-point arithmetic using the same order of operations on the same hardware should yield identical results for every run (that is, there should be complete reproducibility). Jun 8 at 2:24

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.

• Thanks! No I am not sure that my SVD method works. You said you have an example where it fails? Would you mind providing it?
– Joe
May 2 at 4:19
• @Joe I just generated matrices at random, and noted that $\det S_1\neq \det A_1$, and also $\det S_1 / \det S \neq \det A_1 / \det A$, in case you were just comparing ratios. May 2 at 6:21
• Is $S$ what I called $\Sigma$? In that case, indeed the determinants that you mentioned should not match. Instead, one needs to rotate the vector $(\det S_1, … , \det S_n)$ by the matrix $V$. And to be clear, $A_i$ is gotten from $A$ by replacing the $i$th row with $c$, while $S_i$ is gotten from $S$ by replacing the $i$th row by $U^T c$. Does this still give you incorrect results?
– Joe
May 2 at 8:38
• Sorry, things seems to match now, up to a sign; I confirm that your method works in practice! Yes, I mis-described (and possibly also mis-computed) things. May 2 at 14:33
• I have also noticed that sometimes overall signs don’t match with the SVD method. I’m a bit confused why that would be.
– Joe
May 2 at 15:20