Perform mldivide between 3x3 matrix M and every RGB pixel in a image in GPU

I have a rgb image (2048x2048x3), and I wan to perform mldivide on every rgb (3,1) pixels with a M matrix (3x3).

So each pixel will do a M\pixel operation. The easy way to do this is using two for loops to go over each pixels. However, I am wondering if I can use arrayfun or pagefun for this so I can use GPU.

Right now in order to use arrayfun, I have to hard code mldivide so M becomes 9x4194304(2048*2048) array which takes lots of vram. I am wondering if there is a way I can run this without having a big M gpu array as input.

[KSSV]

But 3x3\1x3 will throw an error right?

M = rand(3) ;
pixel = rand(1,3) ;
M\pixel

Error using \
Matrix dimensions must agree.

[Chi Huang]

Sorry I made a mistake, it it will be a 3x1 matrix, rand(3,1). I forgot matlab has y,x convention.

[Stephen23]

"Sorry I made a mistake, it it will be a 3x1 matrix, rand(3,1). I forgot matlab has y,x convention."

... I guess you mean all of mathematics has that convention: https://en.wikipedia.org/wiki/Matrix_(mathematics)

Answer 1

[Stephen23]

format compact
S = 5;
M = rand(3,3);
I = rand(S,S,3);
% reference:
A = nan(S,S,3);
for ii = 1:S
    for jj = 1:S
        pixel = reshape(I(ii,jj,:),3,1);
        A(ii,jj,:) = M\pixel;
    end
end
A

A = 
A(:,:,1) =
    4.2226   15.3976    9.7294    0.6686   -8.5764
   15.5251   -0.8818    3.5689    4.1730    2.4676
    7.1266   -6.9632    9.8103   10.5705   15.3470
    7.5849    0.9940   16.1186   -0.2064   16.8039
   -7.8875   -7.1051    0.2250    2.7415   -7.6742
A(:,:,2) =
    0.7674    3.0061    1.4504    1.4626    0.3589
    2.4699    2.0878    2.4253    1.5477    2.3345
    1.0364    0.5316    1.4698    2.3068    2.5930
    1.5519    2.5403    3.6515    1.6914    2.3027
    0.6423    0.3699    1.7695    2.2491    0.0037
A(:,:,3) =
   -1.0156   -6.3149   -3.2293   -0.7714    2.6126
   -5.7948   -0.6883   -2.1061   -2.1869   -1.4797
   -2.2634    1.9633   -3.2872   -4.0813   -5.4468
   -2.4617   -1.5182   -6.7929   -0.7229   -5.7964
    2.5722    2.1096   -1.0220   -1.5365    2.8622

Remember that RESHAPE is a very fast operation (it does not change your data in memory, only the shape meta-data in its header). In contrast you want to avoid operations that rearrange your data in memory (e.g. TRANSPOSE, PEMUTE). Note that using PAGEFUN would require permutimh your data.

% RESHAPE() and TRANSPOSE() and MLDIVIDE():
B = reshape((M\reshape(I,S*S,3).').',S,S,3)

B = 
B(:,:,1) =
    4.2226   15.3976    9.7294    0.6686   -8.5764
   15.5251   -0.8818    3.5689    4.1730    2.4676
    7.1266   -6.9632    9.8103   10.5705   15.3470
    7.5849    0.9940   16.1186   -0.2064   16.8039
   -7.8875   -7.1051    0.2250    2.7415   -7.6742
B(:,:,2) =
    0.7674    3.0061    1.4504    1.4626    0.3589
    2.4699    2.0878    2.4253    1.5477    2.3345
    1.0364    0.5316    1.4698    2.3068    2.5930
    1.5519    2.5403    3.6515    1.6914    2.3027
    0.6423    0.3699    1.7695    2.2491    0.0037
B(:,:,3) =
   -1.0156   -6.3149   -3.2293   -0.7714    2.6126
   -5.7948   -0.6883   -2.1061   -2.1869   -1.4797
   -2.2634    1.9633   -3.2872   -4.0813   -5.4468
   -2.4617   -1.5182   -6.7929   -0.7229   -5.7964
    2.5722    2.1096   -1.0220   -1.5365    2.8622

% we can simplify the above into RESHAPE() and MRDIVIDE():
C = reshape((reshape(I,S*S,3)/M.'),S,S,3) % recommended

C = 
C(:,:,1) =
    4.2226   15.3976    9.7294    0.6686   -8.5764
   15.5251   -0.8818    3.5689    4.1730    2.4676
    7.1266   -6.9632    9.8103   10.5705   15.3470
    7.5849    0.9940   16.1186   -0.2064   16.8039
   -7.8875   -7.1051    0.2250    2.7415   -7.6742
C(:,:,2) =
    0.7674    3.0061    1.4504    1.4626    0.3589
    2.4699    2.0878    2.4253    1.5477    2.3345
    1.0364    0.5316    1.4698    2.3068    2.5930
    1.5519    2.5403    3.6515    1.6914    2.3027
    0.6423    0.3699    1.7695    2.2491    0.0037
C(:,:,3) =
   -1.0156   -6.3149   -3.2293   -0.7714    2.6126
   -5.7948   -0.6883   -2.1061   -2.1869   -1.4797
   -2.2634    1.9633   -3.2872   -4.0813   -5.4468
   -2.4617   -1.5182   -6.7929   -0.7229   -5.7964
    2.5722    2.1096   -1.0220   -1.5365    2.8622

[Chi Huang]

wow thx,

Shuld have think about side of the out of the diea of looping the pixels.

Answer 2

[Joss Knight]

I feel like I'm missing something - this is just a single backslash with multiple right-hand sides, or to avoid permutation a single mrdivide (edit: you still have to transpose M, but that's very quick):

[h,w] = size(im,[1 2]);
imout = reshape(reshape(im,[],3)/(M.'), h, w, 3);

[Stephen23]

"I feel like I'm missing something"

Yes, the relationship between MLDIVIDE and MRDIVIDE.

"this is just a single backslash with multiple right-hand sides, or to avoid permutation a single mrdivide:"

Almost, in fact the difference is already explained in the MRDIVIDE and MLDIVIDE documentation (which I strongly recommend reading). But in any case, lets compare:

format compact
S = 5;
M = rand(3,3);
I = rand(S,S,3);
% reference:
A = nan(S,S,3);
for ii = 1:S
    for jj = 1:S
        pixel = reshape(I(ii,jj,:),3,1);
        A(ii,jj,:) = M\pixel;
    end
end
A

A = 
A(:,:,1) =
    0.7193    0.8243    0.5550   -0.0297   -0.0509
    0.5016    0.7826    0.0961    0.6376    0.6199
   -0.1862   -0.0229   -0.0054   -0.5392    0.1112
    0.4610    0.2288    1.0311    1.0132   -0.0740
    0.0678    0.7926   -0.1280    0.2075    0.8620
A(:,:,2) =
   -0.2069   -0.3187    1.3863    2.3973    0.7306
   -2.8763   -1.0749    2.4043    1.1967   -0.6490
    3.0769    2.7335    0.8889    2.8601   -0.0210
    0.0231   -1.3710   -0.7220   -0.2142    0.6677
   -0.1136   -0.8085    0.6678    0.0182   -1.8847
A(:,:,3) =
    0.6982    0.2779   -0.2764   -1.4822    0.2922
    2.7752    1.5345   -1.4558   -0.1879    1.2744
   -1.7259   -1.6548   -0.0137   -1.6193    0.7842
    0.2625    1.8839    1.1231    0.5188    0.2445
    0.4455    0.6565   -0.2994    0.9861    2.0038

% Joss Knight's untested code which gives a different result:
[h,w] = size(I,[1 2]);
B = reshape(reshape(I,[],3)/M, h, w, 3)

B = 
B(:,:,1) =
   -0.6455   -2.4000    0.1715    0.5034    2.4232
    0.8047   -0.1307    0.1996   -0.2022    0.4604
    1.6167    0.6411    1.7098    2.6048    1.8655
   -0.6421    2.0287   -1.4028   -2.0945    2.2354
    0.8471   -2.1673    0.8408    2.2212   -0.6079
B(:,:,2) =
    0.7500    0.6188    0.6066   -0.0697    0.2641
    0.8322    0.9388    0.0313    0.6563    0.8174
   -0.1187   -0.0569    0.2048   -0.3734    0.4047
    0.4400    0.6268    1.0242    0.8784    0.2126
    0.2087    0.6386   -0.0554    0.5671    0.9978
B(:,:,3) =
    1.1052    2.4277    0.7798    0.2389   -1.5846
   -0.8369    0.5733    0.5795    1.0767    0.1053
   -0.5286    0.2353   -0.9909   -1.6324   -1.2260
    0.9117   -1.5806    1.8133    2.4211   -1.4896
   -0.5670    2.1032   -0.5488   -1.3754    0.7788

Similar? Not at all. In fact, very very different results.

I already showed the solution using MRDIVIDE two days ago (the line is clearly marked with "Recommended"), it is very similar to what you wrote, with some small differences:

• it actually gives the same result (within FP error) as the OP's reference approach using two loops.
• it follows the relationship given in the MATLAB documentation.

[Joss Knight]

Right, but that's just because I forgot to transpose M for the mrdivide version (fixed in the original answer).

[h,w] = size(I,[1 2]);
B = reshape(reshape(I,[],3)/(M.'), h, w, 3);

max(abs(A(:)-B(:)))

ans =

8.8818e-16

The two versions are dramatically different in performance so there's no way you should use the looping version if you can avoid it:

im = im2double(imresize(imread('peppers.png'),[2048 2048]));
M = randn(3);

imoutStephen = Stephen23sVersion(im,M);
imoutJoss = JosssVersion(im,M);
differenceBetweenAnswers = max(abs(imoutStephen(:)-imoutJoss(:)))

differenceBetweenAnswers =

3.3307e-16

timeForStephen = timeit(@()Stephen23sVersion(im,M))

timeForStephen =

31.1067

timeForJoss = timeit(@()JosssVersion(im,M))

timeForJoss =

0.1709

So your solution is 182x slower for Chi's Huang's problem size.

[Bruno Luong]

Despite of what people keep saying not using INV (the argument never totally convince me), I think this is where INV can be used without any issue

I = im2double(imresize(imread('peppers.png'),[2048 2048]));
M = randn(3);

tic
[h,w] = size(I,[1 2]);
B = reshape(reshape(I,[],3)/(M.'), h, w, 3);
toc % Joss

Elapsed time is 0.184491 seconds.

tic
[h,w] = size(I,[1 2]);
B = reshape(reshape(I,[],3)*inv(M.'), h, w, 3);
toc % Bruno

Elapsed time is 0.035146 seconds.

[Stephen23]

"Right, but that's just because I forgot to transpose M for the mrdivide version (fixed in the original answer)."

You mean that after you posted the wrong code (which you did not bother to test and someone else had to test on your behalf) you have now managed to copy the approach that I gave two days before you even posted your incorrect answer? (namely MRDIVIDE with a transposed M matrix)

So... nearly three days to copy someone's existing approach, without attribution... very impressive work.

"So your solution is 182x slower for Chi's Huang's problem size."

Ummm... whose solution would that be? Why are you hiding these supposed "solutions", why not let us see them?

Making claims about speed using test functions that you did not upload or show anyone... very convincing.

"The two versions are dramatically different in performance so there's no way you should use the looping version if you can avoid it:"

Luckily I didn't recommend a loop.

[Joss Knight]

It just looks like an honest misunderstanding Stephen. You posted looping code which I now see is a reference version. Then you reposted it when you were explaining to me the difference between MLDIVIDE and MRDIVIDE, both of which are functions I personally author in MATLAB.

So I sincerely apologise for not reading your post properly and noticing your code matching mine and I concede that you got to the answer first and have no intention of stealing that accolade.

I also apologise for the honest mistake in not transposing M. It is common practice to edit mistakes in Answers given in this forum to ensure that nobody gets confused, but I did clearly flag the edit.

It didn't seem necessary to repost the two implementations I used, since they were clearly available in the post and your prior one. However, to ensure full disclosure, here they are:

function imout = Stephen23sVersion(im,M)
[h,w] = size(im,[1 2]);
imout = nan(size(im));
for ii = 1:h
    for jj = 1:w
        pixel = reshape(im(ii,jj,:),3,1);
        imout(ii,jj,:) = M\pixel;
    end
end

function imout = JosssVersion(im,M)
[h,w] = size(im,[1 2]);
imout = reshape(reshape(im,[],3)/(M.'), h, w, 3);
end

[Joss Knight]

Bruno, we generally advise against inv for reasons of numerical stability and accuracy. inv(A) is (literally) the solution to A\eye(size(A,1)), which means there is no pivoting during backsubstitution to ensure numerical stability. Essentially, what you are getting out of mldivide is a separate solution to each system of equations, numerically balanced for that right-hand side. When you know your system is well-conditioned, however, it's probably fine to use inv, especially when it's giving you performance benefits.

[Bruno Luong]

Correct me if I'm wrong but

As I understand the backslash x = A\b does that internally

% [Q,R,P] = qr(A);

% c = Q'*b
% y = R\b; % back subsititution where the effective rank is cuts from i such that abs(R(i,i)) >= tol*R(1,1)
% x = P'*y;

So actually it mainly depends on A in my point of view.

I make this test to see if I can see any difference between "\" and inv() even for ill-conditioned matrxi A and honestly I don't see anything you pretent to be a potential issue, or to me "\" is not more robust than "inv()" (sorry for 300 warning messages, you need to scroll down to see the graph):

n = 5;
numtest = 100;
relerrx12 = zeros(1,numtest);
relerrx13 = zeros(1,numtest);
condtarget = 1e10;
for k = 1:numtest
    % Generate random A with condition number ~ 1e10
    [~,U]=qr(rand(n));
    [~,V]=qr(rand(n));
    S=diag(logspace(0,-log(condtarget),n));
    A=U*S*V';
    % fprint('cond(A) = %f\n', cond(A));
    
    % Generate random rhs
    b=randn(n,1);
    
    x1=A\b;
    
    x2=inv(A)*b;
   
    [Q,R,P] = qr(A);
    iA = P*inv(R)*Q'; % EDIT: bug fix
    x3=iA*b;
    
    relerrx12(k) = norm(x1-x2,2)/norm(x1,2);
    relerrx13(k) = norm(x1-x3,2)/norm(x1,2);
end

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.351927e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.351927e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.351929e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.015422e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.015422e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.015424e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.273102e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.273102e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.273104e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.718002e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.718002e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.718012e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.717390e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.717390e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.717391e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.282214e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.282214e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.282222e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.913614e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.913614e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.913625e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.385603e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.385603e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.385614e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.961429e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.961429e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.961433e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.632683e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.632683e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.632685e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.722869e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.722869e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.722870e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.664052e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.664052e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.664058e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.047569e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.047569e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.047570e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.915664e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.915664e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.915672e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.025172e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.025172e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.025173e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.311005e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.311005e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.311007e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.129441e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.129441e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.129443e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.349298e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.349298e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.349302e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.358772e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.358772e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.358773e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.508930e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.508930e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.508931e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.304268e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.304268e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.304272e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.742637e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.742637e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.742646e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.500471e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.500471e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.500477e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.208903e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.208903e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.208905e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.220874e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.220874e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.220875e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.639888e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.639888e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.639890e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.373887e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.373887e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.373892e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.286178e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.286178e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.286180e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.016567e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.016567e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.016569e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.814374e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.814374e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.814380e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.767163e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.767163e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.767167e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.897976e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.897976e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.897977e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.710342e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.710342e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.710344e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.285219e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.285219e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.285220e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.050556e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.050556e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.050565e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.042485e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.042485e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.042489e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.201109e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.201109e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.201110e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.231589e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.231589e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.231591e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.339934e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.339934e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.339935e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.265787e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.265787e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.265788e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.492448e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.492448e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.492451e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.536611e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.536611e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.536613e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.670475e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.670475e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.670476e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.488415e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.488415e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.488416e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.006400e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.006400e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.006407e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.087629e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.087629e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.087629e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.293629e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.293629e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.293634e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.088514e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.088514e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.088516e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.749474e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.749474e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.749476e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.576523e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.576523e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.576532e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.108253e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.108253e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.108256e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.584905e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.584905e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.584906e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.783064e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.783064e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.783066e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.368129e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.368129e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.368135e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.377319e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.377319e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.377320e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.962485e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.962485e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.962489e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.546735e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.546735e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.546736e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.885805e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.885805e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.885813e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.895576e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.895576e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.895578e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.515101e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.515101e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.515112e-27.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.859788e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.859788e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.859795e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.127998e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.127998e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.127999e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.099736e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.099736e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.099737e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.387907e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.387907e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.387908e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.549194e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.549194e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.549197e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.597175e-28.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.597175e-28.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.597177e-28.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.587253e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.587253e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.587257e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.230982e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.230982e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.230983e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.177299e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.177299e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.177301e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.266159e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.266159e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.266161e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.110531e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.110531e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.110535e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.898803e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.898803e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.898814e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.877006e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.877006e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.877009e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.983905e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.983905e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.983911e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.877873e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.877873e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.877876e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.782964e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.782964e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.782976e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.478958e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.478958e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.478960e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.528100e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.528100e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.528107e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.233357e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.233357e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.233358e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.128693e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.128693e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.128694e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.595941e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.595941e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.595947e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.607702e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.607702e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.607703e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.650565e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.650565e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.650565e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.193406e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.193406e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.193407e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.081795e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.081795e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.081796e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.137365e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.137365e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.137372e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.819046e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.819046e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.819048e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.587779e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.587779e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.587782e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.310237e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.310237e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.310237e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.194711e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.194711e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.194712e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.497856e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.497856e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.497862e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.491136e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.491136e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.491139e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.393771e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.393771e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.393774e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.306608e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.306608e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.306611e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.997209e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.997209e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.997212e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.713936e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.713936e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.713939e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.763256e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.763256e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.763259e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.258626e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.258626e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.258626e-26.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.078047e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.078047e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.078048e-24.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.551038e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.551038e-25.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.551041e-25.

figure
h = semilogy(1:numtest, relerrx12, 'b', 1:numtest, relerrx13, 'r');
legend(h, 'Relative Error x1 and x2', 'Relative Error x1 and x3');

[Joss Knight]

Thanks for the detailed analysis! Apologies if I don't spend time picking apart your code - it looks sound. But to summarize, it's not about the conditioning of A alone but the whole system including the RHS. It's the same reason that

isequal(2^53+1+1, 1+1+2^53)

ans =

logical

0

Which I'm sure you understand is a basic principle of floating point numbers. If your RHS contains a wide range of values (in terms of scale) then you benefit from pivoting such that as a general rule you will get a better result using backslash than inv:

for i = 1:1000
    A = randn(1000);
    b = randn(1000,1);
    res1(i) = norm(A*(A\b)-b,2);
    res2(i) = norm(A*(inv(A)*b)-b,2);
end
mean(res1)
mean(res2)

ans =

1.3011e-10

ans =

2.5798e-10

For any system matrix it would be possible to contrive a pathological RHS or set of RHSs that give worse precision for inv, but it's harder the other way around.

For what it's worth, the MATLAB documentation's example is probably way more convincing than mine.

The performance is easier explain without a toy example. Wherever the number of right-hand sides is lower than the number of variables you should expect backslash to perform better, since it is performing triangular solves fewer times.

A = randn(10000);
b = randn(10000,1);
timeit(@()A\b)
timeit(@()inv(A)*b)

ans =

2.5859

ans =

5.9034

Nevertheless, your point is not lost on me. Of course, in this situation (Chi Huang's situation) we have a great many more right-hand sides than variables.

[Bruno Luong]

"Which I'm sure you understand is a basic principle of floating point numbers. If your RHS contains a wide range of values (in terms of scale) then you benefit from pivoting such that as a general rule you will get a better result using backslash than inv:"

No I don' totally agree with it. When one deals with matrix with large condition number (are we agree that there is no difference between inv and "\" in terme of accuracy when matrix is well conditioned?, I hope we do) we are dealing with ill-conditioning problem and it is well-known that reduce the residue to 0 at all cost is a very bad thing to do in the presence of noise on the rhs (the TMW example completely ignores this possibility). Any good book on inverse-problem and ill-conditioned system will tell you that, and the first recommend method is using cutting rank (as with pinv) or use regularization, both methods with give a bigger residu, but far robust in term of accuracy of the solution. Look at the residu is thus not a relevant thing to do.

The example of TMW shows kind of extrem condition number of 1e10 and inv gives indeed the smaller residu and barely better accuracy of the solution (err_bs = 4.2512e-06 vs err_inv 5.3997e-06). For smaller condition number and smaller dimension, the difference is more like numerical noise as I show in my script above.

Of course I know that inv(A) is O(n^3) computation and A\b is O(n^2). The computation cost is only advantage for "\" when n >= 100 for a single rhs, for multiple rhs I have shown inv can be 10 time faster with OP case.

The fact is INV can have advantage in some case, and should not be systematically considered as evil as I often seen people do.

(On other hand I also seen often people use INV for no good reasons).

But I really want to rehabilitate inv that can be a good function to be used in some circumstance.

[Joss Knight]

(are we agree that there is no difference between inv and "\" in terme of accuracy when matrix is well conditioned?, I hope we do)

No, I can't agree with that, indeed my example code demonstrates it is false. inv(A) solves the system for the set of right-hand sides eye(n), which will result in a good solution for those numbers, but those numbers may have a completely different magnitude to the numbers in b and every time two floating point numbers of different magnitudes meet, you risk losing precision, hence my making the point that isequal(2^53+1+1, 1+1+2^53) is false. The subsequent matrix multiplication is doing reductions and this will necessarily lose precision. I may not be as versed in the academics of numerical analytics as you but I'm confident with my floating point numerics!

I concede the point that sometimes inv is useful. I started on this journey merely to attempt to explain why the advice is as it is. I'm pretty sure that nowhere does MathWorks advise to avoid inv at all costs. In my example above I make no attempt to game the result by selecting a poorly conditioned system and yet I still get better residuals for backslash than for inv. Feels pretty convincing to me. Although I can see you have brought into question the validity of the residual as a measurement of accuracy.

[Bruno Luong]

Then please provide a concrete example of well conditioned matrix A (says cond A < 1000) and b (pick whatever you like, arbitray large or small) where

inv(A)*b
A\b

are significan't different in the sense that;

norm(inv(A)*b - A\b) / norm(A\b) >= 1e-6 % I just fix 1e-6 as significant difference

I really don't understand quite your point on 2^53 + 1 + 1... We are talking about matrix and linear algebra, not problem with finite precision floating point. It might linked in the sense that the numercial error will propagete but I don't see clearly what is your point in the discussion, sorry.

[Joss Knight]

2^53+1+1 is [1 1 1]*[2^53;1;1] which is what a matrix multiplication is doing. inv(A) introduces these reductions and misses out on the pivoting inherent in the triangular solve, since there is no pivoting for A\eye(n).

norm(inv(A)*b - A\b) / norm(A\b)

Is this really what you mean here? This isn't the residual error, it's the relative difference between the two answers.

As for your request, I thought that I had done that, but I can make it more explicit:

n = 1000;
for i = 1:100
    Q = orth(randn(n,n));
    d = logspace(0,-3,n);
    A = Q*diag(d)*Q';
    cond_num(i) = cond(A);
    b = randn(1000,1);
    res1(i) = norm(A*(A\b)-b);
    res2(i) = norm(A*(inv(A)*b)-b);
end
max(cond_num)
mean(res1)
mean(res2)

ans =

1.0000e+03

ans =

7.0523e-11

ans =

1.2468e-10

[Bruno Luong]

Not the residual, I want the error norm on the solution

The residual consideration is a wrong thing to focus on IMO. You should read on overfitting that arises on ill-donditioning problem as I suggested ealier.

At the end what we want is the accuracy of the solution.

[Bruno Luong]

"2^53+1+1 is [1 1 1]*[2^53;1;1] which is what a matrix multiplication is doing. inv(A) introduces these reductions and misses out on the pivoting inherent in the triangular solve, since there is no pivoting for A\eye(n)."

This is not about pivoting, it's about order of summing. You can pivot and still run into such difference of results, which I call "finite precision numerical error noise".

Such error is amplified at most by condition number of the matrix, this is well known :

"The condition number is defined more precisely to be the maximum ratio of the relative error in x to the relative error in b. "

https://en.wikipedia.org/wiki/Condition_number

So if the matrix has condition number = kappa, the relative error of the solution is

norm(inv(A)*b) <= kappa * norm(xtrue) * eps(1)

I challenge to find the example because this niquelity tell that for A with cond(1000) the relative error is NEVER below

eps*1000

ans = 2.2204e-13

The solution obtained by inv(A) (or any algorithm) will provide at least 13 digits of precision. That is largely enough for most application. Not a big deal to malke a fuss about avoiding it.

[Joss Knight]

The residual is the error norm on the solution. There's no way to understand the relative importance of each element of b for a general problem, we do not have that information. We could look at the maximum error? But that just gives the same result.

When you requested this:

norm(inv(A)*b - A\b) / norm(A\b) >= 1e-6

You are essentially saying that you'll only agree that inv gives a worse solution if it is potentially 5 or 6 orders of magnitude worse, because that's what 1e-6 is compared to 1e-11. I don't really agree with that evaluation. I don't agree that the two answers are approximately the same unless their difference is down at machine precision for FP64, which means 1e-11 or better for this size of problem.

[Bruno Luong]

"The residual is the error norm on the solution"

No the residual is the error of the rhs, or the fit, it leaves in b-space, not x-space.

As example : If you fit a polynomial to a data point, the residual norm is the error of the data; what I'm interested is the error of the polynomial coefficients.

If you don't like my 1e-6 I'll make it tighter 1e-12.

norm(inv(A)*b - A\b) / norm(A\b) >= 1e-12

with cond(A) <= 1000;

Can you find A and b with such property?

At the end you might be able to find some exemple with 1e-14 relative difference but not larger. That means inv(x) is able to give 13 digits precision on the solution. You might be call it bad, but it good enough for most case to me.

[Joss Knight]

To say that you only care about x-space error is being a bit naive about real-world use cases. In most real engineering problems, in my experience, I don't care as much about the precision of the input variables as I do about the precision of the output. For instance, in a control system, I care that the input parameters give me the output that I want. If a slightly different control input is equidistant from the true solution but means my car is further from the centre of the lane, it is a worse solution.

Still, I reran my example code using error norm as the measure and still got the same results:

n = 1000;
for i = 1:100
    Q = orth(randn(n,n));
    d = logspace(0,-3,n);
    A = Q*diag(d)*Q';
    x = randn(n,1);
    cond_num(i) = cond(A);
    b = A*x;
    res1(i) = norm((A\b)-x);
    res2(i) = norm((inv(A)*b)-x);
end
max(cond_num)
mean(res1)
mean(res2)

ans =

1.0000e+03

ans =

7.0544e-11

ans =

1.2458e-10

inv is still nearly 2x worse than backslash. I don't really know how else to say it. I guess what you're saying is, it's quite good enough for your purposes.

[Bruno Luong]

In general one never have ill-posed system in control system, because otherwise your system is not really controllable. (If that is the case it is quite easy to change earth climate, unfortunately it doesn't work likke that.)

If you have expreince in control you should know the the measurement (b-space) would rarely have a precision up to few digits (I think a record of physics measurement is about 11 or 12 digits).

So yeah in REAL world 13 digits is enough IMO.

The ill-posed problem occurs in other circumtance other than control, and I invite you to read about overfitting and convince yourself that reduce the residu in this case is a bad thing to do, because you blow up the solutions to fit to the noise.

You can stick with your mesure of residu, but you can never convince me that is the right thing yo do.

[Joss Knight]

Excellent. Well, this has been an education, thank you!

[Paul]

Setting aside accuracy considerations for the moment, is one of A\b or inv(A)*b more computationally efficient than the other?

[Bruno Luong]

Also a side note, the pivoting order perfomed by A\B is fully dictated by A, and not B as showed in this example where all the 0s of X are at the same indices:

A=rand(3,5)

A = 3×5

0.7398 0.4359 0.9511 0.9418 0.8110 0.6184 0.0447 0.5276 0.4571 0.3606 0.4980 0.9869 0.5944 0.0123 0.1086

B=rand(3,10)

B = 3×10

0.5389 0.6791 0.1235 0.5291 0.2913 0.2468 0.4506 0.0579 0.3858 0.1901 0.0083 0.7803 0.6447 0.8726 0.2500 0.8429 0.2431 0.1708 0.0072 0.6199 0.9360 0.3331 0.1131 0.5518 0.5435 0.5719 0.0392 0.1861 0.7293 0.3692

X = A\B

X = 5×10

0 0 0 0 0 0 0 0 0 0 1.2893 -1.4892 -2.0578 -1.9928 -0.1708 -2.3693 -0.0741 -0.4377 0.9427 -1.7485 -0.5776 3.0676 3.6609 4.2960 1.2152 4.9716 0.1821 1.0564 -0.3450 3.5780 0.5587 -1.6877 -2.6135 -2.8544 -0.8389 -3.6621 0.3288 -0.8028 0.3217 -2.6022 0 0 0 0 0 0 0 0 0 0

[Bruno Luong]

For b that is vector,

• A\b complexity O(n^2),
• inv(A)*b has O(n^3),

and you find various tic/toc along the discussion above.

For b that is n x m

• A\b complexity O(n^2)+O(n*m),
• inv(A)*b has O(n^3)+O(n*m),

with the O(n*m) much favarable to the inv(), possibly 6-10 time faster for small n and large m as for the OP's question here.

[Paul]

@Bruno Luong

No worries. I just deleted my comment that flagged your edit and it looks like Stephen23 deleted his related comment. If you delete your comment and can delete this comment as well, we'll go back to having a nice, clean thread as if it never happened.

Answer 3

[Bruno Luong]

B = reshape(reshape(I,[],3)*inv(M.'),size(I));