Matrix inversion differences between versions

I have been working with some relatively large arrays (approx. 32000x32000, with a sparsity of >0.99) as part of a modelling project. I first wrote a version of my code in 2018a, and taking the inverse of my matrix, A, and multiplying this with a vector, b, worked fine using the "A\b" commmand.

I got a new laptop today and finally got around to updating to version 2020a. The same matrix inversion now doesn't work - I get the warning that the "matrix is singular to working precision". After this warning has been displayed, the operation "A\b" returns an array comprised entirely of NaN.

I tried running the exact same code on my desktop to check (using 2018a), and it worked as normal. I then installed version 2018a on my new laptop, and again the code ran fine. So I am quite sure that it is to do with some difference between 2018a and 2020a.

Does anybody know anything about this difference between the two versions? I would like to be using the most up to date release of MATLAB, but it seems that my code won't run. Any insight would be much appreciated.

**Note: I have not provided my code here as it is part of an academic study and it would not be appropriate for me to release it.

Answer 1

[Bjorn Gustavsson]

Have you tried to look at the smallest singular values? You can calculate those with svds. If they are small or very small your calculation A\b is not necessarily all that meaningful.

In my fields we are now required to give open access to our data.

[Adam Ward]

Hi Bjorn,

Thank you for your speedy response! I hope you are keeping safe and well wherever in the world you are.

The svds are indeed very small, they're on the order of 10^-13. I'm guessing this is due to both the fact that the matrix is very large, as well as being highly sparse. Although, I was of the impression that as long as they are non-zero then the matrix should be invertible. Indeed, 2018a seems able to calculate this inverse and the results do give a reasonable physical interpretation in our problem.

My question was more focussed on the difference between 2018a and 2020a. Specifically, why 2018a would be able to calculate this inverse, but 2020a would not.

[Bjorn Gustavsson]

OK, that looks traditionally troublesome. Unless the largest eigenvalues are similarly small you have a problem. The size of the smallest eigenvalues has "nothing" to do with the size or sparsity of your matrix A. It has to do with the matrix being close to singular. In a reality-based context it is not enough that the eigenvalues are non-zero - since you have data with noise or some finite precision (which is the same as discretization-noise.)

Linear discrete inverse-problem theory in 21 lines:

b = A*x,

We can factorize A in two orthonormal base-vector-matrices, U and V as: A = U*S*V.' where S is a diagonal matrix with the eigenvalues along the diagonal. This gives us:

b = U*S*V.'*x,

This allows us to solve for x even in the case where A is singular (in a regularized sense) by stepwise multiply from the left on both sides:

U.'*b = U.'*U*S*V.'*x == S*V.'*x => 
inv(S)*U.'*b = inv(S)*S*V.'*x == V.'*x => 
V*inv(S)*U.'*b = V*V.'*x == x.

The problem with real-world data in b is that the first equation is wrong - it should look like this:

b = A*x + n,  

where n is your noise (measurement, discretization or other)

this changes everything to give us:

V*inv(S)*U.'*b - V*inv(S)*U.'*n = V*V.'*x == x.

Here the first term on the left-hand-side is what we want, but the second term is the problematic, there the noise components that happen to be multiplied by the inverse of thhe smaller eigenvalues will be amplified - in your case by 1/1e-13.

Assuming that the larger eigenvalues are on the order of unity I think I would not be rude to suggest that you don't have your data in b with that small uncertainties. In your case you might have found that matlab 2018 used some sensible iterative method with some implicit damping (or stopping-criteria) such that you did not get noise-only solution in x, while matlab2020 might have started taking a more conservative approach and bailed.

My suggestion is that you preferably look for some of the explicit solvers to systems of sparse linear equations:

sparse matrices (including eq-solvers)

[John D'Errico]

Note that what is important is to consider the RATIO of the smallest singular values with the largest ones. That is, if the largest singular value is 1e-5, and the smallest is 1e-13, then the matrix is not close to singular. But if the largest singular value is 1e+5, then that ratio spans 18 orders of magnitude. That matrix is singular. The most extreme case is of course

M = 1e-13*eye(32000);

which is perfectly well behaved.

[Adam Ward]

Thank you for the detailed response, it is very helpful.

[Bjorn Gustavsson]

To follow up on Johns comment: I should have made that a bit more clear. However, it is still just the beginning of the next chapter in practical application of linear inverse-problem theory. There are two fundamental facts that are often neglected in introductory linear algebra courses:

1. There will always be some kind of noise in the data, b in this case, even if there is no measurement noise the AD-conversion will have some limited number of bits worth of accuracy.
2. For the case where we have a discretized version of a continuous unknown, x in my answer, and discrete data (for example tomographic reconstruction-problems) the problem is underdetermined as soon as we have our discrete data, and the "measurement" or "theory" matrix has infinite condition-number. If we hide that by under-discretize the model-space, i.e. the dimension of x, we have implicitly constrained/regularized the solution but this does not change the fundamental character of the problem.

That we need to take a few steps beyond the condition-number is obvious if we look at two modified version of John's example. Let us have two data-vectors:

b1 = [0.99876
      0.99885
      0.99496
      0.99317
      1.0027
      1.0222
      6.6246e-14];

and

b1p = [7.0761e-14
       8.738e-14
       9.8103e-14
       1.2058e-13
       8.1381e-14
       9.5989e-14
       7.4862e-14];

Both produced from the same unknowns array:

x = ones(7,1);

but with two different theory-matrices:

M7 = diag([ones(1,6),1e-13]);
M7e = diag([ones(1,6)*1e-13,1e-13]);

The first has a condition-number of 1e13 the second has a condition-number of 1. If we have normal-distributed measurement noise with widths that is max of 1% and 4/2^16 (just to model something with finite AD-precision):

b1 = M7*ones(7,1);
b1 = b1 + [0.01*randn(6,1);0.25*randn(1)*1e-13];
b1p = M7e*ones(7,1);
b1p = b1p + [0.25(randn(7,1))].*b1p;  

Our x-estimates will be far better from the first problem (x1 = M7\b1) than from the second (x1p = M7e\b1p), even though the condition-number points in the other direction. However, the last component of x1 will have much larger variance than the first 6.

For the case where we do not under-discretize the theory-matrix it will now be singular, i.e. with some number of non-zero eigenvalues, lets say p. Then the condition-number is infinity, but this is both expected and irellevant since the problem was underdetermined the moment we got our finite set of measurements. To get an estimate of our unknown we can naively try the general (Moore-Penrose) svd-inverse:

x = V_p*inv(S_p)*U_p.'*b_obs;

Where U_p is the p rows of U corresponding to the non-zero eigenvalues etc. But since our observations, b_obs, are the sum of the ideal measurements and some unavoidable noise:

b_obs = b_ideal + n

we see that with

x = V*inv(S)*U.'*b_ideal + V*inv(S)*U.'*n;

a problem with the components of U.'*n that corresponds to the really small eigenvalues will be hugely amplified, and here the amplified noise component will, in general, spread to most components in x (in contrast to the previous nice diagonal problem). Instead of using the "raw" inverse of S we can resort to either a further truncation of S (and U and V) or use a damped version (Tikhonov-inverse) where we use

alpha = 1; % Or some other suitable value
iZ_alpha = diag(diag(S)./(diag(S).^2+alpha));
x = V_p*iZ*U_p.'*b_obs;

Here we can continue and calculate the covariance matrix for x and the resolution-matrix, but not at this moment.

If the theory-matrix is too big for SVD-decomposition we have to resort to iterative methods and hope that they behave well - at least some do.