Ax=b

From Wikimization

(Difference between revisions)
Jump to: navigation, search
(Seven ways Matlab can solve thin-matrix linear equality)
Current revision (20:46, 20 December 2017) (edit) (undo)
 
(4 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>
+
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
Line 107: Line 107:
</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) can increase 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