In article <1992May11.070411.15451@dutrun2.tudelft.nl> rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:
>In article <1992May8.181301.9243@news.cs.indiana.edu> pete shirley:
>> [ Linear eqs. ] What is the fast way to solve this?  For 2 by 2
>> I use Cramer's rule, but several people have told me for N > 2 that
>> is stupid.  What is fastest?
>
>For small matrices up to and including 4 x 4, Cramer's rule is fastest.
>-- 
>* Han de Bruijn; Applications&Graphics | "A little bit of Physics * No
>* TUD Computing Centre; P.O. Box 354   | would be NO idleness in  * Oil
>* 2600 AJ  Delft; The Netherlands.     | Mathematics" (HdB).      * for
>* E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood

Cramer's rule is not a good idea, even for 2x2 systems.  I am
attaching a posting I made some time ago which shows why.  For garden
variety problems, Gaussian elimination with partial pivoting is the
algorithm of choice.  If you expect ill-conditioned problems, you can
use a condition estimator or any one of a number of rank revealing
decompositions.  A bibliography is appended. (Incidentally, my
bibliographic data base along with a program to search on keywords is
available by anonymous ftp from thales.cs.umd.edu in the directory
pub/references.)

Pete Stewart

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An earlier posting, in which I noted that Cramer's rule is unstable,
drew some queries.  Since the question is a hardy perennial, I am
appending a short matlab code for people to play with.  It generates
an ill-conditioned system whose exact solution is x = [1,2]', rounds
it a bit, and solves it twice--once by Cramer's rule and again by
Gaussian elimination.  The printout consists of the condition number,
and the solution and its scaled residual for each algorithm.
Here is sample output (slightly edited).

conda = 5.7268e+14

       xc           rc           xg           rg

   1.0164e+00  -1.5707e-03   1.0302e+00  -3.7064e-17
   2.0000e+00  -4.9346e-03   1.9697e+00            0

Note that the solution by Cramer's rule (xc) is slightly more accurate
than the solution by Gaussian elimination (xg).  Nonetheless, the
residuals (rc and rg) show that the solution by Cramer's rule
satisfies the equations to only three figures, whereas the solution by
Gaussian elimination satisfies it right up to the rounding unit.  The
latter is a consequence of the backward stability of Gaussian
elimination, a property that Cramer's rule does not share.

One caveat.  The above output is not entirely typical.  Because a two
by two system gives little opportunity for rounding error to dig in,
the residuals for Gaussian elimination are often zero.  Moreover,
Gaussian elimination sometimes produces the exact solution, as does
Cramer's rule--though less frequently.

Pete Stewart

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a Matlab M-file.  To execute it, start up Matlab in this
% directory and type 
%
%	>> Cramer
%
% at the Matlab prompt.
%
% Cramer.m generates an ill-conditioned two by two system whose exact
% solution is x = [1,2]'.  Then it rounds the matrix a bit, and solves
% the system twice--once by Cramer's rule and again by Gaussian
% elimination.  It prints out the condition number, and the solution and
% its scaled residual for each algorithm.

format short e;
rand('normal');

fac = 1.e13;        % Condition control

a = rand(2);        % Get a random a with the second row equal
a(2,:) = pi*a(1,:); % to pi times the first, and fix it to
a = fix(fac*a);     % log10(fac) digits.

conda = cond(a)     % Print the condition number of a.

x = [1,2]';         % Form a right hand side with solution
b = a*x;            % whose components are 1 and 2.

a = a/fac;          % Introduce some initial rounding error
b = b/fac;          % by scaling back.

% Solution and residual by Cramer's rule

d = a(1,1)*a(2,2) - a(2,1)*a(1,2);
xc = [b(1)*a(2,2)-b(2)*a(1,2), a(1,1)*b(2)-a(2,1)*b(1)]'/d;
rc = (b - a*xc)/(norm(a)*norm(xc));

% Solution and residual by Gaussian elimination

xg = a\b;
rg = (b - a*xg)/(norm(a)*norm(xg));

[xc, rc, xg, rg]    % Print out the results

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{bisc:90}
C.~H. Bischof.
\newblock Incremental condition estimation.
\newblock {\it SIAM Journal on Matrix Analysis and Applications}, 11:312--322,
  1990.

\bibitem{biha:89}
C.~H. Bischof and P.~C. Hansen.
\newblock {\it Structure-Preserving and Rank-Revealing {QR}-Factorizations}.
\newblock Preprint~MCS-P100-0989, Mathematics and Computer Science Division,
  Argonne National Laboratory, 1989.

\bibitem{bilp:89}
C.~H. Bischof, J.~G. Lewis, and D.~J. Pierce.
\newblock {\it Incremental Condition Estimation for General Matrices}.
\newblock Preprint~MCS-P106-0989, Mathematics and Computer Science Division,
  Argonne National Laboratory, 1989.

\bibitem{bipl:90}
C.~H. Bischof, D.~J. Pierce, and J.~G. Lewis.
\newblock Incremental condition estimation for sparse matrices.
\newblock {\it SIAM Journal on Matrix Analysis and Applications}, 11:644--659,
  1990.
\newblock Cited in {\AA ke Bj\"orck's} bibliography on least squares. Available
  by anonymous ftp from math.liu.se in pub/references.

\bibitem{bishr:89}
C.~H. Bischof and G.~M. Shroff.
\newblock {\it On Updating Signal Subspaces}.
\newblock Preprint~MCS-P101-0989, Mathematics and Computer Science Division,
  Argonne National Laboratory, 1989.

\bibitem{chfp:77}
S.~P. Chan, R. Feldman, and B.~N. Parlett.
\newblock Algorithm 517: a program for computing the condition numbers of
  matrix eigenvalues without computing eigenvectors.
\newblock {\it ACM Transactions on Mathematical Software}, 3:186--203, 1977.

\bibitem{chan:87}
T.~F. Chan.
\newblock Rank revealing {QR}~factorizations.
\newblock {\it Linear Algebra and Its Applications}, 88/89:67--82, 1987.

\bibitem{chha:90}
T.~F. Chan and P.~C. Hansen.
\newblock Computing truncated singular value decomposition least squares
  solutions by rank revealing {QR}-factorizations.
\newblock {\it SIAM Journal on Scientific and Statistical Computing},
  11:519--530, 1990.

\bibitem{chha:91}
T.~F. Chan and P.~C. Hansen.
\newblock Low-rank revealing {QR}~factorizations.
\newblock 1990.
\newblock Manuscript, Department of Mathematics, University of California, Los
  Angeles.

\bibitem{chha:92}
T.~F. Chan and P.~C. Hansen.
\newblock Some applications of the rank revealing {QR}~factorization.
\newblock {\it SIAM Journal on Scientific and Statistical Computing},
  13:727--741, 1992.

\bibitem{clcv:82}
A.~K. Cline, A.~R. Conn, and C.~F. {Van Loan}.
\newblock Generalizing the {LINPACK} condition estimator.
\newblock In J.P. Hennart, editor, {\it Numerical Analysis}, Springer Verlag,
  1982.
\newblock Cited in {\AA ke Bj\"orck's} bibliography on least squares. Available
  by anonymous ftp from math.liu.se in pub/references.

\bibitem{cmsw:79}
A.~K. Cline, C.~B. Moler, G.~W. Stewart, and J.~H. Wilkinson.
\newblock An estimate for the condition number of a matrix.
\newblock {\it SIAM Journal on Numerical Analysis}, 16:368--375, 1979.

\bibitem{goks:76}
G.~H. Golub, V. Klema, and G.~W. Stewart.
\newblock {\it Rank Degeneracy and Least Squares Problems}.
\newblock Technical Report~TR-456, Department of Computer Science, University
  of Maryland, 1976.

\bibitem{gore:70}
G.~H. Golub and C. Reinsch.
\newblock Singular value decomposition and least squares solution.
\newblock {\it Numerische Mathematik}, 14:403--420, 1970.
\newblock Also in \cite[pp.134--151]{wire:71}.

\bibitem{grst:76}
W.~B. Gragg and G.~W. Stewart.
\newblock A stable variant of the secant method for solving nonlinear
  equations.
\newblock {\it SIAM Journal on Numerical Analysis}, 13:880--903, 1976.

\bibitem{hage:84}
W.~W. Hager.
\newblock Condition estimators.
\newblock {\it SIAM Journal on Scientific and Statistical Computing},
  5:311--316, 1984.
\newblock Cited in \cite{high:87}.

\bibitem{high:87}
N.~J. Higham.
\newblock A survey of condition number estimation for triangular matrices.
\newblock {\it SIAM Review}, 29:575--596, 1987.

\bibitem{hopa:92}
Y.~P. Hong and C.-T. Pan.
\newblock Rank-revealing {QR}~factorizations and the singular value
  decomposition.
\newblock {\it Mathematics of Computation}, 58:213--232, 1992.

\bibitem{kaha:72}
W. Kahan.
\newblock {\it Conserving Confluence Curbs Ill-Conditioning}.
\newblock Technical Report~6, Computer Science Department, University of
  California, Berkeley, 1972.

\bibitem{nats:83}
M. Natori and A. Tsukamoto.
\newblock A fast method for estimating the condition number of a matrix.
\newblock {\it Journal of Information Processing}, 6:138--140, 1983.

\bibitem{pipl:92}
D.~J. Pierce and R.~J. Plemmons.
\newblock Fast adaptive condition estimation.
\newblock {\it SIAM Journal on Matrix Analysis and Applications}, 13:274--291,
  1992.

\bibitem{shbi:90}
G.~M. Shroff and C.~H. Bischof.
\newblock {\it Adaptive Condition Estimation for Rank-One Updates of
  QR~Factorizations}.
\newblock Preprint~MCS-P166-0790, Mathematics and Computer Science Division,
  Argonne National Laboratory, 1990.

\bibitem{stew:80}
G.~W. Stewart.
\newblock The efficient generation of random orthogonal matrices with an
  application to condition estimators.
\newblock {\it SIAM Journal on Numerical Analysis}, 17:403--404, 1980.

\bibitem{stew:84}
G.~W. Stewart.
\newblock Rank degeneracy.
\newblock {\it SIAM Journal on Scientific and Statistical Computing},
  5:403--413, 1984.

\bibitem{stew:91a}
G.~W. Stewart.
\newblock {\it Updating a Rank-Revealing {ULV}~Decomposition}.
\newblock Technical Report~CS-TR 2627, Department of Computer Science,
  University of Maryland, 1991.

\bibitem{stew:90b}
G.~W. Stewart.
\newblock {\it An Updating Algorithm for Subspace Tracking}.
\newblock Technical Report~CS-TR 2494, Department of Computer Science,
  University of Maryland, 1990.
\newblock To appear in IEEE Transactions on Signal Processing.

\bibitem{vdwe:79}
A. van~der~Sluis and G.~W. Veltkamp.
\newblock Restoring rank and consistency by orthogonal projection.
\newblock {\it Linear Algebra and Its Applications}, 28:257--278, 1979.

\bibitem{wire:71}
J.~H. Wilkinson and C. Reinsch.
\newblock {\it Handbook for Automatic Computation. Vol.~II Linear Algebra}.
\newblock Springer, New York, 1971.

\bibitem{xuka:90}
G. Xu and T. Kailath.
\newblock Fast signal-subspace decompsotion\,---\,part i: {Ideal} covariance
  matrices.
\newblock 1990.
\newblock Manuscript submitted to ASSP. Information Systems Laboratory,
  Stanford University.

\bibitem{xuka:90a}
G. Xu and T. Kailath.
\newblock Fast signal-subspace decomposition\,---\,part {II}: {Sample}
  covariance matrices.
\newblock 1990.
\newblock Manuscript submitted to ASSP. Information Systems Laboratory,
  Stanford University.