2
$\begingroup$

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$.
$\endgroup$
1
  • $\begingroup$ 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). $\endgroup$
    – njuffa
    Jun 8, 2023 at 2:24

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – Joe
    May 2, 2023 at 4:19
  • $\begingroup$ @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. $\endgroup$ May 2, 2023 at 6:21
  • $\begingroup$ 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? $\endgroup$
    – Joe
    May 2, 2023 at 8:38
  • $\begingroup$ 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. $\endgroup$ May 2, 2023 at 14:33
  • $\begingroup$ I have also noticed that sometimes overall signs don’t match with the SVD method. I’m a bit confused why that would be. $\endgroup$
    – Joe
    May 2, 2023 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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