3

Based on the wonderful MATLAB documentation, it is clear that when solving systems of linear equations, it is preferable to use X\b to inv(X)*b for reasons of numerical stability. It also notes that X^(-1) is equivalent to inv(X). But it makes no mention of the relationship between X\eye(size(X)) and inv(X)! I would like to use the more succinct notation, and it makes sense that inv(X) should be the preferred form, but cannot find confirmation of this anywhere in the documentation.

The docs mention that \ uses Gaussian elimination whereas inv uses LU decomposition.

When I tested them using an example similar to that provided in the documenation, I find that the two are not equivalent:

n = 500;
Q = orth(randn(n,n));
d = logspace(0,-10,n);
A = Q*diag(d)*Q';

y = inv(A);
z = A\eye(n);
disp(['Difference: ', num2str(norm(inv(A) - A\eye(n)))]);

The results:

Difference: 2.2957e-05

Can anyone explain which of the two methods is more numerically stable?

TimD1
  • 982
  • 15
  • 26
  • 5
    Use `timeit` to make meaningful comparisons between runtimes. Single `tic` and `toc` calls are essentially useless for benchmarking purposes. – sco1 Jul 25 '18 at 18:17
  • Thanks, I did not know that! Just removed that, since I'm more interested in numerical stability anyways. – TimD1 Jul 25 '18 at 18:22
  • I guess that you have already ran this a lot of times to find the numerical stability and timing differences? Please update the question with those data. – Adriaan Jul 25 '18 at 18:25
  • My problem here is that I can get the differences, but that information is useless since I don't know which one is "more correct". I would need to use a matrix where the inverse is known beforehand to benchmark, but then I fear my results would only be applicable to how those algorithms perform on that specific type of matrix... – TimD1 Jul 25 '18 at 18:27
  • 2
    `\ ` uses a huge amount of algorithms, not Gaussian elimination. However, the real question you need to ask is: do you really need the matrix inverse? Or are you using this inverse somewhere else? Because in the second case, you want to use `\ `, in that "somwhere else" place, and that would be by far the most accurate and stable method – Ander Biguri Jul 25 '18 at 18:35
  • 2
    Does it really matter in practice, or is this a question of academic nature? My advice to you is - whenever you have a choice between code that "does the job better" versus code that "looks better", you should keep the more performant code active and possibly keep the other option as a comment nearby. Maybe both of them are "good enough"? How did you reach the conclusion that `\\` uses Gaussian elimination? [It depends on the input](https://www.mathworks.com/help/matlab/ref/mldivide.html#bt4jslc-6)! – Dev-iL Jul 25 '18 at 18:36
  • @Dev-iL I found it on another MathWorks documentation page: "A better way, from the standpoint of both execution time and numerical accuracy, is to use the matrix backslash operator x = A\b. This produces the solution using Gaussian elimination, without explicitly forming the inverse." – TimD1 Jul 25 '18 at 18:42
  • 1
    You can alternative go to the main documentation page of `\ ` and look at the algorithms: https://uk.mathworks.com/help/matlab/ref/mldivide.html#bt4jslc-6 – Ander Biguri Jul 25 '18 at 18:45
  • There are various ways to go about this, it's a question of what you want to achieve at the end of the day. See for example: [`decomposition`](https://www.mathworks.com/help/matlab/ref/decomposition.html). – Dev-iL Jul 25 '18 at 18:46

0 Answers0