Ax=b

From Wikimization

(Difference between revisions)
Jump to: navigation, search
Current revision (20:46, 20 December 2017) (edit) (undo)
 
(6 intermediate revisions not shown.)
Line 1: Line 1:
-
=Seven ways Matlab can solve thin-matrix linear equality=
+
==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> by least squares such that <math>A_{}x\approx b\,</math>:
-
When <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
-
A = randn(1e6,52);
+
m=1e6; n=52;
-
% spA = sparse(A);
+
A = randn(m,n);
-
xact = randn(52,1);
+
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);
+
bb(i) = csum(A(:,i).*b);
end
end
 +
% AA = A'*A;
 +
% AA = (AA + AA')/2;
spAA = sparse(AA);
spAA = sparse(AA);
-
for i=1:4
+
% 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
-
% x1 = A\b;
+
% x1a = A\b;
% toc
% toc
-
% disp(['error = ' num2str(norm(x1 - xact)/norm(xact))])
+
% disp(['error = ' num2str(norm(x1a - xact)/norm(xact))])
disp(' backslash on A''A')
disp(' backslash on A''A')
tic
tic
-
x1a = AA\bb;
+
x1 = AA\bb;
toc
toc
-
disp(['error = ' num2str(norm(x1a - xact)/norm(xact))])
+
disp(['error = ' num2str(norm(x1 - xact)/norm(xact))])
% disp(' pinv on A')
% disp(' pinv on A')
% tic
% tic
-
% x2 = pinv(A)*b;
+
% x2a = pinv(A)*b;
% toc
% toc
-
% disp(['error = ' num2str(norm(x2 - xact)/norm(xact))])
+
% disp(['error = ' num2str(norm(x2a - xact)/norm(xact))])
disp(' pinv on A''A')
disp(' pinv on A''A')
tic
tic
-
x2a = pinv(AA)*bb;
+
x2 = pinv(AA)*bb;
toc
toc
-
disp(['error = ' num2str(norm(x2a - xact)/norm(xact))])
+
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);
-
% x3 = R\c;
+
% x3a = R\c;
% toc
% toc
-
% disp(['error = ' num2str(norm(x3 - xact)/norm(xact))])
+
% 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);
-
x3a = R\c;
+
x3 = R\c;
toc
toc
-
disp(['error = ' num2str(norm(x3a - xact)/norm(xact))])
+
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
-
% x6 = linsolve(A,b);
+
% x6a = linsolve(A,b);
% toc
% toc
-
% disp(['error = ' num2str(norm(x6 - xact)/norm(xact))])
+
% disp(['error = ' num2str(norm(x6a - xact)/norm(xact))])
disp(' linsolve on A''A')
disp(' linsolve on A''A')
tic
tic
-
x6a = linsolve(AA,bb,opt);
+
x6 = linsolve(AA,bb,opt);
toc
toc
-
disp(['error = ' num2str(norm(x6a - xact)/norm(xact))])
+
disp(['error = ' num2str(norm(x6 - xact)/norm(xact))])
% disp(' lscov on A')
% disp(' lscov on A')
% tic
% tic
-
% x7 = lscov(A,b);
+
% x7a = lscov(A,b);
% toc
% toc
-
% disp(['error = ' num2str(norm(x7 - xact)/norm(xact))])
+
% disp(['error = ' num2str(norm(x7a - xact)/norm(xact))])
disp(' lscov on A''A')
disp(' lscov on A''A')
tic
tic
-
x7a = lscov(AA,bb);
+
x7 = lscov(AA,bb);
toc
toc
-
disp(['error = ' num2str(norm(x7a - xact)/norm(xact))])
+
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]]
+
[[Accumulator_Error_Feedback | <tt>csum()</tt> routine]]
-
(with presorting) increases precision by orders of magnitude.
+
(with presorting) can increase precision by orders of magnitude.

Current revision

Seven ways Matlab can optimally solve normal equations given thin matrix LaTeX: A

For LaTeX: b\notin\mathcal{R}(A)\,,\, find best fit LaTeX: x by least squares such that LaTeX: A_{}x\approx b\,:

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

Personal tools