Ax=b
From Wikimization
(Difference between revisions)
| Line 1: | Line 1: | ||
==Seven ways Matlab can optimally solve normal equations given thin matrix <math>A</math>== | ==Seven ways Matlab can optimally solve normal equations given thin matrix <math>A</math>== | ||
| - | For <math>b\notin\mathcal{R}(A)\,,\,</math> find best fit <math>x</math> such that <math>A_{}x\approx b\,</math>: | + | For <math>b\notin\mathcal{R}(A)\,,\,</math> find best fit <math>x</math> by least squares such that <math>A_{}x\approx b\,</math>: |
<pre> | <pre> | ||
%test backslash timing | %test backslash timing | ||
Current revision
Seven ways Matlab can optimally solve normal equations given thin matrix 
For find best fit
by least squares such that
:
%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.