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

Question

Chi Huang 2022-7-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754945-perform-mldivide-between-3x3-matrix-m-and-every-rgb-pixel-in-a-image-in-gpu

评论： Paul 2022-7-12

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 个评论
显示 1更早的评论隐藏 1更早的评论

Chi Huang 2022-7-8

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

Stephen23 2022-7-8

"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 2022-7-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754945-perform-mldivide-between-3x3-matrix-m-and-every-rgb-pixel-in-a-image-in-gpu#answer_1002735

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Chi Huang 2022-7-8

wow thx,

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

请先登录，再进行评论。

Answer 2

Joss Knight 2022-7-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754945-perform-mldivide-between-3x3-matrix-m-and-every-rgb-pixel-in-a-image-in-gpu#answer_1003770

编辑：Joss Knight 2022-7-10

在 MATLAB Online 中打开

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 个评论
显示 21更早的评论隐藏 21更早的评论

Stephen23 2022-7-10

编辑：Stephen23 2022-7-10

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

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

Bruno Luong 2022-7-10

编辑：Bruno Luong 2022-7-11

在 MATLAB Online 中打开

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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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 = 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.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.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.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.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 = 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.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.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.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.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.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.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.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 = 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.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 2022-7-10

编辑：Joss Knight 2022-7-10

在 MATLAB Online 中打开

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

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

编辑：Joss Knight 2022-7-11

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

Joss Knight 2022-7-11

在 MATLAB Online 中打开

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

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

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

请先登录，再进行评论。