Under what circumstances can two (nearly) identical sparse matrices give different solutions to Mx = b?

Question

Suppose I have two sparse matrices, $A$ and $B$, of size $N \times N$. They each have the same sparcity pattern ("footprint"). They each also have values which in theory should be identical, but aren't due to machine precision error.

A simple example of this in MATLAB syntax:

% Dimensions of square, sparse matrices
N = 10;
A = sparse(N,N);
B = sparse(N,N);

% Add a dozen nonzero entries at random locations
numEntries = 12;
rows = randi( N, numEntries,1 );
cols = randi( N, numEntries,1 );

% Populate the entries, one at a time
for entry = 1:numEntries
    i = rows(entry);
    j = cols(entry);

    % Generate random numbers
    x = rand(1,1);
    y = rand(1,1);

    % Populate matrices
    % Note: (x^2-y^2) / (x-y) == (x+y)*(x-y)/(x-y) == x+y
    A(i,j) = (x^2-y^2) / (x-y);
    B(i,j) = x + y;

end

diffMat = A - B;
disp(diffMat);

Obviously the values of $A$ and $B$ should be identical, but they aren't. The output of this script is as follows:

   1.0e-15 *

   (1,2)     -0.222044604925031
   (2,2)      0.222044604925031
  (10,3)      0.111022302462516
   (8,4)     -0.222044604925031
   (3,5)      0.222044604925031
  (10,5)      0.555111512312578
   (3,6)      0.222044604925031

Suppose a pair of matrices like this that give wildly different results $x_A$ and $x_B$ for the matrix equations $A x_A = b$ and $B x_B = b$ (same $b$ in each case). Suppose I have confirmed that the values of $A$ and $B$ are identical to less than one part in 1e15 (that is, abs( (A(A~=0)-B(B~=0)) ./ A(A~=0) ) < 1e-15). Further suppose that the estimated condition number (condest()) of each matrix is identical.

My question: Under what circumstances could this occur?

Disclaimer

I originally asked this question because I had this problem in my code. It turned out that my $b$ vectors were not identical due to a small typo in the code that generated them (thanks @JesseChan!).

I am still interested if and how this can occur. My guess is that very high condition numbers can produce these kinds of results, but I'm curious if anyone can produce and explain a minimum working example of this behavior.

Answer 1

This cannot happen, unless the matrices are very ill-behaved. If you have two nonsingular matrices $A, B$ with $$ A-B = O(\epsilon), $$ then the relative error in solving $Ax=b$ is $$ \begin{aligned} \frac{\|A^{-1}b-B^{-1}b\|}{\|A^{-1}b\|} &\leq \sup_x \frac{\|A^{-1}x-B^{-1}x\|}{\|A^{-1}x\|} \\&\leq \sup_y \frac{\|y - B^{-1}Ay\|}{\|y\|} \\&\leq \|I - B^{-1}A\| \\&\leq \|B^{-1}(B-A)\| \\&\leq \|B^{-1}\| \|A-B\| \\&\leq \frac{\kappa(B)}{\|B\|}\|A-B\| \end{aligned} $$

So if $A$ and $B$ are identical to within $\epsilon$, and also have, as you say, identical condition numbers $\kappa(A)=\kappa(B)$, then the relative error can be at most $ \kappa(A)\epsilon/\|A\|. $

In particular, for the error to be on the order of $O(1)$, the condition number of $A$ must be on the order of $\|A\|\epsilon^{-1}$.

Answer 2

Glad the answer was in the $b$ vector. The answer to the original equation depends on exactly what method you use to solve the system (backslash has a ton of different options built into it). You usually have that the relative difference in your solutions is $O(\kappa \epsilon)$, where $\epsilon$ is the difference between the matrices and $\kappa$ is the condition number of $A,B$. An easier way of thinking about it is that the order of magnitude of your conditioning illustrates roughly how many digits difference you can expect (for a stable method of solving the system).