**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

7 views (last 30 days)

Show older comments

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.

##### 3 Comments

Chi Huang
on 8 Jul 2022

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

Stephen23
on 8 Jul 2022

"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)

### Accepted Answer

Stephen23
on 8 Jul 2022

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

##### 1 Comment

Chi Huang
on 8 Jul 2022

wow thx,

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

### More Answers (2)

Joss Knight
on 9 Jul 2022

Edited: Joss Knight
on 10 Jul 2022

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);

##### 23 Comments

Stephen23
on 10 Jul 2022

Edited: Stephen23
on 10 Jul 2022

"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
on 10 Jul 2022

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
on 10 Jul 2022

Edited: Bruno Luong
on 10 Jul 2022

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
on 10 Jul 2022

Edited: Stephen23
on 10 Jul 2022

"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
on 10 Jul 2022

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
on 10 Jul 2022

Bruno Luong
on 10 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

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
on 10 Jul 2022

Edited: Joss Knight
on 10 Jul 2022

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.

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
on 10 Jul 2022

"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
on 11 Jul 2022

Edited: Joss Knight
on 11 Jul 2022

(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
on 11 Jul 2022

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
on 11 Jul 2022

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
on 11 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

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
on 11 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

"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. "

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
on 11 Jul 2022

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
on 11 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

"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
on 11 Jul 2022

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
on 11 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

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.

Bruno Luong
on 12 Jul 2022

Edited: Bruno Luong
on 12 Jul 2022

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
on 12 Jul 2022

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
on 12 Jul 2022

Bruno Luong
on 11 Jul 2022

Edited: Bruno Luong
on 11 Jul 2022

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

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)