Ax=b
From Wikimization
(Difference between revisions)
Line 1: | Line 1: | ||
=Seven ways Matlab can solve thin-matrix linear equality= | =Seven ways Matlab can solve thin-matrix linear equality= | ||
- | + | For <math>b\notin\mathcal{R}(A)</math> | |
- | + | ||
- | + | ||
<pre> | <pre> | ||
%test backslash timing | %test backslash timing | ||
clc %clear all; close all; fclose all; slow execution by order of magnitude | clc %clear all; close all; fclose all; slow execution by order of magnitude | ||
- | + | m=1e6; n=52; | |
- | + | A = randn(m,n); | |
- | xact = randn( | + | xact = randn(n,1); |
- | b = A*xact; | + | b = A*xact; |
opt.SYM = true; opt.POSDEF = true; | opt.SYM = true; opt.POSDEF = true; | ||
opt1.LT = true; opt2.UT = true; | opt1.LT = true; opt2.UT = true; | ||
- | %form A'*A | ||
- | [m n] = size(A); | ||
AA = zeros(n,n); | AA = zeros(n,n); | ||
bb = zeros(n,1); | bb = zeros(n,1); | ||
+ | %form A'A for accuracy, precision, and symmetry | ||
for i=1:n | for i=1:n | ||
for j=1:i | for j=1:i | ||
Line 23: | Line 20: | ||
AA(j,i) = AA(i,j); | AA(j,i) = AA(i,j); | ||
end | end | ||
- | + | bb(i) = csum(A(:,i).*b); | |
end | end | ||
+ | % AA = A'*A; | ||
+ | % AA = (AA + AA')/2; | ||
spAA = sparse(AA); | spAA = sparse(AA); | ||
- | for i=1: | + | % spA = sparse(A); |
+ | for i=1:3 %Manually disable unnecessary OS processes. Runs faster on subsequent loops. | ||
if i > 1, disp(' ');disp(' '); end | if i > 1, disp(' ');disp(' '); end | ||
% disp(' backslash on A') | % disp(' backslash on A') | ||
% tic | % tic | ||
- | % | + | % x1a = A\b; |
% toc | % toc | ||
- | % disp(['error = ' num2str(norm( | + | % disp(['error = ' num2str(norm(x1a - xact)/norm(xact))]) |
disp(' backslash on A''A') | disp(' backslash on A''A') | ||
tic | tic | ||
- | + | x1 = AA\bb; | |
toc | toc | ||
- | disp(['error = ' num2str(norm( | + | disp(['error = ' num2str(norm(x1 - xact)/norm(xact))]) |
% disp(' pinv on A') | % disp(' pinv on A') | ||
% tic | % tic | ||
- | % | + | % x2a = pinv(A)*b; |
% toc | % toc | ||
- | % disp(['error = ' num2str(norm( | + | % disp(['error = ' num2str(norm(x2a - xact)/norm(xact))]) |
disp(' pinv on A''A') | disp(' pinv on A''A') | ||
tic | tic | ||
- | + | x2 = pinv(AA)*bb; | |
toc | toc | ||
- | disp(['error = ' num2str(norm( | + | disp(['error = ' num2str(norm(x2 - xact)/norm(xact))]) |
% disp(' QR on A') | % disp(' QR on A') | ||
% tic | % tic | ||
% [c,R] = qr(spA,b,0); | % [c,R] = qr(spA,b,0); | ||
- | % | + | % x3a = R\c; |
% toc | % toc | ||
- | % disp(['error = ' num2str(norm( | + | % disp(['error = ' num2str(norm(x3a - xact)/norm(xact))]) |
disp(' QR on A''A') | disp(' QR on A''A') | ||
tic | tic | ||
[c,R] = qr(spAA,bb,0); | [c,R] = qr(spAA,bb,0); | ||
- | + | x3 = R\c; | |
toc | toc | ||
- | disp(['error = ' num2str(norm( | + | disp(['error = ' num2str(norm(x3 - xact)/norm(xact))]) |
disp(' U\(U''\A''*b) on U=chol(A''A)') | disp(' U\(U''\A''*b) on U=chol(A''A)') | ||
Line 83: | Line 83: | ||
% disp(' linsolve on A') | % disp(' linsolve on A') | ||
% tic | % tic | ||
- | % | + | % x6a = linsolve(A,b); |
% toc | % toc | ||
- | % disp(['error = ' num2str(norm( | + | % disp(['error = ' num2str(norm(x6a - xact)/norm(xact))]) |
disp(' linsolve on A''A') | disp(' linsolve on A''A') | ||
tic | tic | ||
- | + | x6 = linsolve(AA,bb,opt); | |
toc | toc | ||
- | disp(['error = ' num2str(norm( | + | disp(['error = ' num2str(norm(x6 - xact)/norm(xact))]) |
% disp(' lscov on A') | % disp(' lscov on A') | ||
% tic | % tic | ||
- | % | + | % x7a = lscov(A,b); |
% toc | % toc | ||
- | % disp(['error = ' num2str(norm( | + | % disp(['error = ' num2str(norm(x7a - xact)/norm(xact))]) |
disp(' lscov on A''A') | disp(' lscov on A''A') | ||
tic | tic | ||
- | + | x7 = lscov(AA,bb); | |
toc | toc | ||
- | disp(['error = ' num2str(norm( | + | disp(['error = ' num2str(norm(x7 - xact)/norm(xact))]) |
end | end | ||
</pre> | </pre> | ||
[[http://www.convexoptimization.com/wikimization/index.php/Accumulator_Error_Feedback|<tt>csum()</tt> routine]] | [[http://www.convexoptimization.com/wikimization/index.php/Accumulator_Error_Feedback|<tt>csum()</tt> routine]] | ||
- | (with presorting) | + | (with presorting) can increase precision by orders of magnitude. |
Revision as of 16:16, 17 December 2017
Seven ways Matlab can solve thin-matrix linear equality
For
%test backslash timing clc %clear all; close all; fclose all; slow execution by order of magnitude m=1e6; n=52; A = randn(m,n); xact = randn(n,1); b = A*xact; opt.SYM = true; opt.POSDEF = true; opt1.LT = true; opt2.UT = true; AA = zeros(n,n); bb = zeros(n,1); %form A'A for accuracy, precision, and symmetry for i=1:n for j=1:i AA(i,j) = csum(A(:,i).*A(:,j)); AA(j,i) = AA(i,j); end bb(i) = csum(A(:,i).*b); end % AA = A'*A; % AA = (AA + AA')/2; spAA = sparse(AA); % spA = sparse(A); for i=1:3 %Manually disable unnecessary OS processes. Runs faster on subsequent loops. if i > 1, disp(' ');disp(' '); end % disp(' backslash on A') % tic % x1a = A\b; % toc % disp(['error = ' num2str(norm(x1a - xact)/norm(xact))]) disp(' backslash on A''A') tic x1 = AA\bb; toc disp(['error = ' num2str(norm(x1 - xact)/norm(xact))]) % disp(' pinv on A') % tic % x2a = pinv(A)*b; % toc % disp(['error = ' num2str(norm(x2a - xact)/norm(xact))]) disp(' pinv on A''A') tic x2 = pinv(AA)*bb; toc disp(['error = ' num2str(norm(x2 - xact)/norm(xact))]) % disp(' QR on A') % tic % [c,R] = qr(spA,b,0); % x3a = R\c; % toc % disp(['error = ' num2str(norm(x3a - xact)/norm(xact))]) disp(' QR on A''A') tic [c,R] = qr(spAA,bb,0); x3 = R\c; toc disp(['error = ' num2str(norm(x3 - xact)/norm(xact))]) disp(' U\(U''\A''*b) on U=chol(A''A)') tic; U = chol(AA); x4 = U\(U'\bb); toc disp(['error = ' num2str(norm(x4 - xact)/norm(xact))]) disp(' linsolve on chol(A''A)') %winner on Dec.14 2017 tic U = chol(AA); x5 = linsolve(U, linsolve(U',bb,opt1), opt2); toc disp(['error = ' num2str(norm(x5 - xact)/norm(xact))]) % disp(' linsolve on A') % tic % x6a = linsolve(A,b); % toc % disp(['error = ' num2str(norm(x6a - xact)/norm(xact))]) disp(' linsolve on A''A') tic x6 = linsolve(AA,bb,opt); toc disp(['error = ' num2str(norm(x6 - xact)/norm(xact))]) % disp(' lscov on A') % tic % x7a = lscov(A,b); % toc % disp(['error = ' num2str(norm(x7a - xact)/norm(xact))]) disp(' lscov on A''A') tic x7 = lscov(AA,bb); toc disp(['error = ' num2str(norm(x7 - xact)/norm(xact))]) end
[csum() routine] (with presorting) can increase precision by orders of magnitude.