Matrix inversion differences between versions

Question

Adam Ward 2020-6-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/539991-matrix-inversion-differences-between-versions

评论： Bjorn Gustavsson 2020-6-2

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bjorn Gustavsson 2020-6-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/539991-matrix-inversion-differences-between-versions#answer_444431

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Bjorn Gustavsson 2020-6-1

在 MATLAB Online 中打开

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)

Adam Ward 2020-6-1

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

Bjorn Gustavsson 2020-6-2

在 MATLAB Online 中打开

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:

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.
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
99885
99496
99317
0027
0222
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.

请先登录，再进行评论。

Matrix inversion differences between versions

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Matrix inversion differences between versions

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论