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

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.

• The plots look like they might be from an error in the imposition of a boundary condition. Do you happen to have Matlab code generating the actual $A,B$ matrices and $b$ vector? Is the $b$ vector the same in both cases as well? – Jesse Chan Aug 11 '15 at 3:49
• What's the condition number of A? B? – Brian Borchers Aug 11 '15 at 3:53
• Are you plotting the entries of xA and xB (as you seem to state on top of the plots) or some function of them (as you seem to state below the plot)? – Federico Poloni Aug 11 '15 at 9:49
• @JesseChan: Turns out my $b$ vector was not the same for each, and that was what caused the problems. I fixed that and everything's working again. That said, I'm curious about this question for general knowledge! – jvriesem Aug 11 '15 at 19:00
• @FedericoPoloni: I plotted the value of the normal derivative, both the fixed value (fixed via the Neumann boundary conditions) and the value I calculated from the solution. Moot point now though. – jvriesem Aug 11 '15 at 19:11

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}$.
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).