Ax=b

(Difference between revisions)
 Revision as of 15:29, 17 December 2017 (edit)← Previous diff 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 $A$== - + For $b\notin\mathcal{R}(A)\,,\,$ find best fit $x$ by least squares such that $A_{}x\approx b\,$: - When $b\notin\mathcal{R}(A)$ + - +

%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

- [[http://www.convexoptimization.com/wikimization/index.php/Accumulator_Error_Feedback|csum() routine]] + [[Accumulator_Error_Feedback | csum() routine]] - (with presorting) increases precision by orders of magnitude. + (with presorting) can increase precision by orders of magnitude.

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.