From Wikimization

Seven ways Matlab can solve thin-matrix linear equality

When

%test backslash timing
clc  %clear all; close all; fclose all; slow execution by order of magnitude
A = randn(600*96000,52);
% spA = sparse(A);
xact = randn(52,1); 
b = A*xact;

opt.SYM = true;  opt.POSDEF = true;
opt1.LT = true;  opt2.UT = true;

[m n] = size(A);
AA = zeros(n,n);
bb = zeros(n,1);
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
spAA = sparse(AA);
for i=1:4
   if i > 1, disp(' ');disp(' '); end

%    disp('   backslash on A')
%    tic
%    x1 = A\b;
%    toc
%    disp(['error = ' num2str(norm(x1 - xact)/norm(xact))])

   disp('   backslash on A''A')
   tic
   x1a = AA\bb;
   toc
   disp(['error = ' num2str(norm(x1a - xact)/norm(xact))])

%    disp('   pinv on A')
%    tic
%    x2 = pinv(A)*b;
%    toc
%    disp(['error = ' num2str(norm(x2 - xact)/norm(xact))])

   disp('   pinv on A''A')
   tic
   x2a = pinv(AA)*bb;
   toc
   disp(['error = ' num2str(norm(x2a - xact)/norm(xact))])

%    disp('   QR on A')
%    tic
%    [c,R] = qr(spA,b,0);  
%    x3 = R\c;
%    toc
%    disp(['error = ' num2str(norm(x3 - xact)/norm(xact))])

   disp('   QR on A''A')
   tic
   [c,R] = qr(spAA,bb,0);  
   x3a = R\c;
   toc
   disp(['error = ' num2str(norm(x3a - 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
%    x6 = linsolve(A,b);
%    toc
%    disp(['error = ' num2str(norm(x6 - xact)/norm(xact))])

   disp('   linsolve on A''A')
   tic
   x6a = linsolve(AA,bb,opt);
   toc
   disp(['error = ' num2str(norm(x6a - xact)/norm(xact))])

%    disp('   lscov on A')
%    tic
%    x7 = lscov(A,b);
%    toc
%    disp(['error = ' num2str(norm(x7 - xact)/norm(xact))])

   disp('   lscov on A''A')
   tic
   x7a = lscov(AA,bb);
   toc
   disp(['error = ' num2str(norm(x7a - xact)/norm(xact))])
end