File: factorize_demo.html

package info (click to toggle)
suitesparse 1%3A5.12.0%2Bdfsg-2
links: PTS, VCS
area: main
in suites: bookworm
size: 176,720 kB
sloc: ansic: 1,193,914; cpp: 31,704; makefile: 6,638; fortran: 1,927; java: 1,826; csh: 765; ruby: 725; sh: 529; python: 333; perl: 225; sed: 164; awk: 35
file content (633 lines) | stat: -rw-r--r-- 30,228 bytes
parent folder | download | duplicates (7)

<!DOCTYPE html
  PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   <!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
      --><title>THE FACTORIZE OBJECT for solving linear systems</title><meta name="generator" content="MATLAB 7.9"><meta name="date" content="2012-07-09"><meta name="m-file" content="factorize_demo"><style type="text/css">

body {
  background-color: white;
  margin:10px;
}

h1 {
  color: #990000; 
  font-size: x-large;
}

h2 {
  color: #990000;
  font-size: medium;
}

/* Make the text shrink to fit narrow windows, but not stretch too far in 
wide windows. */ 
p,h1,h2,div.content div {
  max-width: 600px;
  /* Hack for IE6 */
  width: auto !important; width: 600px;
}

pre.codeinput {
  background: #EEEEEE;
  padding: 10px;
}
@media print {
  pre.codeinput {word-wrap:break-word; width:100%;}
} 

span.keyword {color: #0000FF}
span.comment {color: #228B22}
span.string {color: #A020F0}
span.untermstring {color: #B20000}
span.syscmd {color: #B28C00}

pre.codeoutput {
  color: #666666;
  padding: 10px;
}

pre.error {
  color: red;
}

p.footer {
  text-align: right;
  font-size: xx-small;
  font-weight: lighter;
  font-style: italic;
  color: gray;
}

  </style></head><body><div class="content"><h1>THE FACTORIZE OBJECT for solving linear systems</h1><!--introduction--><p>Copyright 2011-2012, Timothy A. Davis, <a href="http://www.suitesparse.com">&lt;mailto:DrTimothyAldenDavis@gmail.com DrTimothyAldenDavis@gmail.com</a>&gt;</p><p>This is a demonstration of the FACTORIZE object for solving linear systems and least-squares problems, and for computations with the matrix inverse and pseudo-inverse.</p><!--/introduction--><h2>Contents</h2><div><ul><li><a href="#1">Rule Number One: never multiply by the inverse, inv(A)</a></li><li><a href="#2">Rule Number Two:  never break Rule Number One</a></li><li><a href="#3">How to use BACKSLASH to solve A*x=b</a></li><li><a href="#4">BACKSLASH versus INV ... let the battle begin</a></li><li><a href="#8">LU and LINSOLVE are fast and accurate but complicated to use</a></li><li><a href="#9">INV is easy to use, but slow and inaccurate</a></li><li><a href="#11">So the winner is ... nobody</a></li><li><a href="#12">The FACTORIZE object to the rescue</a></li><li><a href="#13">Least-squares problems</a></li><li><a href="#15">Underdetermined systems</a></li><li><a href="#16">Computing selected entries in the inverse or pseudo-inverse</a></li><li><a href="#17">Computing the entire inverse or pseudo-inverse</a></li><li><a href="#18">Update/downdate of a dense Cholesky factorization</a></li><li><a href="#19">Caveat Executor</a></li><li><a href="#20">Summary</a></li></ul></div><h2>Rule Number One: never multiply by the inverse, inv(A)<a name="1"></a></h2><p>Use backslash or a matrix factorization instead (LU, CHOL, or QR).</p><h2>Rule Number Two:  never break Rule Number One<a name="2"></a></h2><p>However, the problem with Rule Number One is that it can be hard to figure out which matrix factorization to use and how to use it.  Using LU, CHOL, or QR is complicated, particularly if you want the best performance.  BACKSLASH (MLDIVIDE) is great, but it can't be reused when solving multiple systems (x=A\b and y=A\c).  Its syntax doesn't match the use of the inverse in mathematical expressions, either.</p><p>The goal of the FACTORIZE object is to solve this problem ...</p><p>"Don't let that INV go past your eyes; to solve that system, FACTORIZE!"</p><h2>How to use BACKSLASH to solve A*x=b<a name="3"></a></h2><p>First, let's create a square matrix A and a right-hand-side b for a linear system A*x=b.  There are many ways to solve this system.  The best way is to use x=A\b.  The residual r is a vector of what's left over in each equation, and its norm tells you how accurately the system was solved.</p><pre class="codeinput">format <span class="string">compact</span> ;
A = rand (3)
b = rand (3,1)
x = A\b
r = b-A*x ;
norm (r)
</pre><pre class="codeoutput">A =
    0.7221    0.6055    0.5715
    0.1602    0.9760    0.1407
    0.0346    0.4527    0.2428
b =
    0.6253
    0.4261
    0.3978
x =
   -0.2529
    0.3238
    1.0707
ans =
     0
</pre><h2>BACKSLASH versus INV ... let the battle begin<a name="4"></a></h2><p>The backslash operation x=A\b is mathematically the same as x=inv(A)*b. However, backslash is faster and more accurate since it uses a matrix factorization instead of multiplying by the inverse.  Even though your linear algebra textbook might write x=A^(-1)*b as the solution to the system A*x=b, your textbook author never means for you to compute the inverse.</p><p>These next statements give the same answer, so what's the big deal?</p><pre class="codeinput">S = inv(A) ;
x = S*b
x = A\b
</pre><pre class="codeoutput">x =
   -0.2529
    0.3238
    1.0707
x =
   -0.2529
    0.3238
    1.0707
</pre><p>The big deal is that you should care about speed and you should care even more about accuracy.  BACKSLASH relies on matrix factorization (LU, CHOL, QR, or other specialized methods).  It's faster and more reliable than multiplying by the inverse, particularly for large matrices and sparse matrices.  Here's an illustration of how pathetic inv(A)*b can be.</p><pre class="codeinput">A = gallery (<span class="string">'frank'</span>,16) ; xtrue = ones (16,1) ; b = A*xtrue ;

x = inv(A)*b ; norm (b-A*x)
x = A\b      ; norm (b-A*x)
</pre><pre class="codeoutput">ans =
    0.0145
ans =
   1.7764e-15
</pre><p>The performance difference between BACKSLASH and INV for even small sparse matrices is striking.</p><pre class="codeinput">load <span class="string">west0479</span> ;
A = west0479 ;
n = size (A,1)
b = rand (n,1) ;
tic ; x = A\b ; toc
norm (b-A*x)
tic ; x = inv(A)*b ; toc
norm (b-A*x)
</pre><pre class="codeoutput">n =
   479
Elapsed time is 0.002701 seconds.
ans =
   1.3584e-10
Elapsed time is 0.118024 seconds.
ans =
   2.4936e-09
</pre><p>What if you want to solve multiple systems?  Use a matrix factorization. But which one?  And how do you use it?  Here are some alternatives using LU for the sparse west0479 matrix, but some are faster than others.</p><pre class="codeinput">tic ; [L,U]     = lu(A) ; x1 = U \ (L \ b)         ; t1=toc ; nz1=nnz(L+U);
tic ; [L,U,P]   = lu(A) ; x2 = U \ (L \ P*b)       ; t2=toc ; nz2=nnz(L+U);
tic ; [L,U,P,Q] = lu(A) ; x3 = Q * (U \ (L \ P*b)) ; t3=toc ; nz3=nnz(L+U);

fprintf (<span class="string">'1: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz1,t1, norm(b-A*x1));
fprintf (<span class="string">'2: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz2,t2, norm(b-A*x2));
fprintf (<span class="string">'3: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz3,t3, norm(b-A*x3));
</pre><pre class="codeoutput">1: nnz(L+U): 16151 time:   0.0056 resid: 1.059567e-10
2: nnz(L+U): 15826 time:   0.0147 resid: 7.104000e-11
3: nnz(L+U):  3704 time:   0.0051 resid: 6.371471e-11
</pre><h2>LU and LINSOLVE are fast and accurate but complicated to use<a name="8"></a></h2><p>A quick look at ``help lu'' will scroll off your screen.  For full matrices, [L,U,p] = lu (A,'vector') is fastest.  Then for the forward/backsolves, use LINSOLVE instead of BACKSLASH for even faster performance.  But for sparse matrices, use the optional 'Q' output of LU so you get a good fill-reducing ordering.  But you can't use 'Q' if the matrix is full.  But LINSOLVE doesn't work on sparse matrices.</p><p>But ... Ack!  That's getting complicated ...</p><p>Here's the best way to solve A*x=b and A*y=c when A is full and unsymmetric:</p><pre class="codeinput">n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;
tic ; [L,U,p] = lu (A, <span class="string">'vector'</span>) ; LUtime = toc

tic ; x = U \ (L \ b (p,:)) ;
      y = U \ (L \ c (p,:)) ; toc

tic ; opL = struct (<span class="string">'LT'</span>, true) ;
      opU = struct (<span class="string">'UT'</span>, true) ;
      x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
      y = linsolve (U, linsolve (L, c(p,:), opL), opU) ; toc
</pre><pre class="codeoutput">LUtime =
    0.2897
Elapsed time is 0.013282 seconds.
Elapsed time is 0.006304 seconds.
</pre><h2>INV is easy to use, but slow and inaccurate<a name="9"></a></h2><p>Oh bother!  Using LU and LINSOLVE is too complicated.  You just want to solve your system.  Let's just compute inv(A) and use it twice.  Easy to write, but slower and less accurate ...</p><pre class="codeinput">S = inv (A) ;
x = S*b ; norm (b-A*x)
y = S*c ; norm (c-A*y)
</pre><pre class="codeoutput">ans =
   2.0238e-11
ans =
   1.8024e-11
</pre><p>Sometimes using the inverse seems inevitable.  For example, your textbook might show the Schur complement formula as S = A-B*inv(D)*C.  This can be done without inv(D) in one of two ways: SLASH or BACKSLASH (MRDIVIDE or MLDIVIDE to be precise).</p><p>inv(A)*B and A\B are mathematically equivalent, as are B*inv(A) and B/A, so these three methods give the same results (ignoring computational errors, which are worse for inv(D)).  Only the first equation looks like the equation in your textbook, however.</p><pre class="codeinput">A = rand (200) ; B = rand (200) ; C = rand (200) ; D = rand (200) ;

tic ; S1 = A - B*inv(D)*C ; toc ;
tic ; S2 = A - B*(D\C) ;    toc ;
tic ; S3 = A - (B/D)*C ;    toc ;
</pre><pre class="codeoutput">Elapsed time is 0.019368 seconds.
Elapsed time is 0.015578 seconds.
Elapsed time is 0.017005 seconds.
</pre><h2>So the winner is ... nobody<a name="11"></a></h2><p>BACKSLASH: mostly simple to use (except remember that Schur complement       formula?).  Fast and accurate ... but slow if you want to solve       two linear systems with the same matrix A.</p><p>LU, QR, CHOL: fast and accurate.  Awful syntax to use.  Drag out your       linear algebra textbook if you want to use these in MATLAB.       Whenever I use them I have to derive them from scratch, even       though I <b>wrote</b> most of the sparse factorizations used in MATLAB!</p><p>INV: slow and inaccurate.  Wins big on ease-of-use, though, since it's a       direct plug-in for all your nice mathematical formulas.</p><p>No method is best on all three criterion: speed, accuracy, and ease of use.</p><p>Is there a solution?  Yes ... keeping reading ...</p><h2>The FACTORIZE object to the rescue<a name="12"></a></h2><p>The FACTORIZE method is just as easy to use as INV, but just as fast and accurate as BACKSLASH, LU, QR, CHOL, and LINSOLVE.</p><p>F = factorize(A) computes the factorization of A and returns it as an object that you can reuse to solve a linear system with x=F\b.  It picks LU, QR, or Cholesky for you, just like BACKSLASH.</p><p>S = inverse(A) is simpler yet.  It does NOT compute inv(A), but factorizes A.  When multiplying S*b, it doesn't mulitply by the inverse, but uses the correct forward/backsolve equations to solve the linear system.</p><pre class="codeinput">n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;

tic ;                       x = A\b ; y = A\c ; toc
tic ; S = inv(A) ;          x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ;    x = F\b ; y = F\c ; toc
tic ; S = inverse(A) ;      x = S*b ; y = S*c ; toc
</pre><pre class="codeoutput">Elapsed time is 0.631824 seconds.
Elapsed time is 0.790060 seconds.
Elapsed time is 0.326646 seconds.
Elapsed time is 0.321894 seconds.
</pre><h2>Least-squares problems<a name="13"></a></h2><p>Here are some different methods for solving a least-squares problem when your system is over-determined.  The last two methods are the same.</p><pre class="codeinput">A = rand (1000,200) ;
b = rand (1000,1) ;

tic ; x = A\b            ; toc, norm (A'*A*x-A'*b)
tic ; x = pinv(A)*b      ; toc, norm (A'*A*x-A'*b)
tic ; x = inverse(A)*b   ; toc, norm (A'*A*x-A'*b)
tic ; x = factorize(A)\b ; toc, norm (A'*A*x-A'*b)
</pre><pre class="codeoutput">Elapsed time is 0.056806 seconds.
ans =
   1.9098e-12
Elapsed time is 0.340721 seconds.
ans =
   2.0025e-12
Elapsed time is 0.079824 seconds.
ans =
   2.6834e-12
Elapsed time is 0.078512 seconds.
ans =
   2.6834e-12
</pre><p>FACTORIZE is better than BACKSLASH because you can reuse the factorization for different right-hand-sides.  For full-rank matrices, it's better than PINV because it's faster (and PINV fails for sparse matrices).</p><pre class="codeinput">A = rand (1000,200) ;
b = rand (1000,1) ;
c = rand (1000,1) ;

tic ;                  ; x = A\b ; y = A\c ; toc
tic ; S = pinv(A)      ; x = S*b ; y = S*c ; toc
tic ; S = inverse(A)   ; x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ; x = F\b ; y = F\c ; toc
</pre><pre class="codeoutput">Elapsed time is 0.112498 seconds.
Elapsed time is 0.341902 seconds.
Elapsed time is 0.080206 seconds.
Elapsed time is 0.082501 seconds.
</pre><h2>Underdetermined systems<a name="15"></a></h2><p>The under-determined system A*x=b where A has more columns than rows has many solutions.  x=A\b finds a basic solution (some of the entries in x are zero).  pinv(A)*b finds a minimum 2-norm solution, but it's slow.  QR factorization will do the same if A has full rank.  That's what the factorize(A) and inverse(A) methods do.</p><pre class="codeinput">A = rand (200,1000) ;
b = rand (200,1) ;

tic ; x = A\b            ; toc, norm (x)
tic ; x = pinv(A)*b      ; toc, norm (x)
tic ; x = inverse(A)*b   ; toc, norm (x)
tic ; x = factorize(A)\b ; toc, norm (x)
</pre><pre class="codeoutput">Elapsed time is 0.075132 seconds.
ans =
    3.9615
Elapsed time is 0.344866 seconds.
ans =
    0.5288
Elapsed time is 0.082219 seconds.
ans =
    0.5288
Elapsed time is 0.081563 seconds.
ans =
    0.5288
</pre><h2>Computing selected entries in the inverse or pseudo-inverse<a name="16"></a></h2><p>If you want just a few entries from the inverse, it's still better to formulate the problem as a system of linear equations and use a matrix factorization instead of computing inv(A).  The FACTORIZE object does this for you, by overloading the subsref operator.</p><pre class="codeinput">A = rand (1000) ;

tic ; S = inv (A)     ; S (2:3,4), toc
tic ; S = inverse (A) ; S (2:3,4), toc
</pre><pre class="codeoutput">ans =
    0.0622
   -0.0401
Elapsed time is 0.788354 seconds.
ans =
    0.0622
   -0.0401
Elapsed time is 0.325303 seconds.
</pre><h2>Computing the entire inverse or pseudo-inverse<a name="17"></a></h2><p>Rarely, and I mean RARELY, you really do need the inverse.  More frequently what you want is the pseudo-inverse.  You can force a factorization to become a plain matrix by converting it to double.  Note that inverse(A) only handles full-rank matrices (either dense or sparse), whereas pinv(A) works for all dense matrices (not sparse).</p><p>The explicit need for inv(A) (or S=A\eye(n), which is the same thing) is RARE.  If you ever find yourself multiplying by the inverse, then you know one thing for sure.  You know with certainty that you don't know what you're doing.</p><pre class="codeinput">A = rand (500) ;
tic ; S1 = inv (A) ;            ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)

A = rand (500,400) ;
tic ; S1 = pinv (A)             ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)
</pre><pre class="codeoutput">Elapsed time is 0.105186 seconds.
Elapsed time is 0.146987 seconds.
ans =
   3.1250e-13
Elapsed time is 1.685084 seconds.
Elapsed time is 0.212152 seconds.
ans =
   3.5864e-14
</pre><h2>Update/downdate of a dense Cholesky factorization<a name="18"></a></h2><p>Wilkinson considered the update/downdate of a matrix factorization to be a key problem in computational linear algebra.  The idea is that you first factorize a matrix.  Next, make a low-rank change to A, and patch up (or down...) the factorization so that it becomes the factorization of the new matrix.  In MATLAB, this only works for dense symmetric positive definite matrices, via cholupdate.  This is much faster than computing the new factorization from scratch.</p><pre class="codeinput">n = 1000 ;
A = rand (n) ;
A = A*A' + n*eye (n) ;
w = rand (n,1) ; t = rand (n,1) ; b = rand (n,1) ;
F = factorize (A) ;

tic ; F = cholupdate (F,w,<span class="string">'+'</span>) ; x = F\b ; toc
tic ; y = (A+w*w')\b ;      toc
norm (x-y)

tic ; F = cholupdate (F,t,<span class="string">'-'</span>) ; x = F\b ; toc
tic ; y = (A+w*w'-t*t')\b ; toc
norm (x-y)
</pre><pre class="codeoutput">Elapsed time is 0.035651 seconds.
Elapsed time is 0.197350 seconds.
ans =
   3.4517e-17
Elapsed time is 0.037860 seconds.
Elapsed time is 0.203125 seconds.
ans =
   3.5007e-17
</pre><h2>Caveat Executor<a name="19"></a></h2><p>One caveat:  If you have a large number of very small systems to solve, the object-oriented overhead of creating and using an object can dominate the run time, at least in MATLAB R2011a.  For this case, if you want the best performance, stick with BACKSLASH, or LU and LINSOLVE (just extract the appropriate formulas from the M-files in the FACTORIZE package).</p><p>Hopefully the object-oriented overhead will drop in future versions of MATLAB, and you can ignore this caveat.</p><pre class="codeinput">A = rand (10) ; b = rand (10,1) ; F = factorize (A) ;

tic ; <span class="keyword">for</span> k = 1:10000, x = F\b ; <span class="keyword">end</span> ; toc

tic ; <span class="keyword">for</span> k = 1:10000, x = A\b ; <span class="keyword">end</span> ; toc

[L,U,p] = lu (A, <span class="string">'vector'</span>) ;
opL = struct (<span class="string">'LT'</span>, true) ;
opU = struct (<span class="string">'UT'</span>, true) ;
tic ;
<span class="keyword">for</span> k = 1:10000
    x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
<span class="keyword">end</span>
toc
</pre><pre class="codeoutput">Elapsed time is 2.953968 seconds.
Elapsed time is 0.350816 seconds.
Elapsed time is 0.125989 seconds.
</pre><h2>Summary<a name="20"></a></h2><p>So ... don't use INV, and don't worry about how to use LU, CHOL, or QR factorization.  Just install the FACTORIZE package, and you're on your way.  Assuming you are now in the Factorize/ directory, cut-and-paste these commands into your command window:</p><pre> addpath (pwd)
 savepath</pre><p>And remember ...</p><p>"Don't let that INV go past your eyes; to solve that system, FACTORIZE!"</p><p class="footer"><br>
      Published with MATLAB&reg; 7.9<br></p></div><!--
##### SOURCE BEGIN #####
%% THE FACTORIZE OBJECT for solving linear systems
%
% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com
% DrTimothyAldenDavis@gmail.com
%
% This is a demonstration of the FACTORIZE object for solving linear
% systems and least-squares problems, and for computations with the
% matrix inverse and pseudo-inverse.

%% Rule Number One: never multiply by the inverse, inv(A)
%
% Use backslash or a matrix factorization instead (LU, CHOL, or QR).

%% Rule Number Two:  never break Rule Number One
%
% However, the problem with Rule Number One is that it can be hard to
% figure out which matrix factorization to use and how to use it.  Using
% LU, CHOL, or QR is complicated, particularly if you want the best
% performance.  BACKSLASH (MLDIVIDE) is great, but it can't be reused when
% solving multiple systems (x=A\b and y=A\c).  Its syntax doesn't match
% the use of the inverse in mathematical expressions, either.
%
% The goal of the FACTORIZE object is to solve this problem ...
%
% "Don't let that INV go past your eyes; to solve that system, FACTORIZE!"

%% How to use BACKSLASH to solve A*x=b
%
% First, let's create a square matrix A and a right-hand-side b for a
% linear system A*x=b.  There are many ways to solve this system.  The
% best way is to use x=A\b.  The residual r is a vector of what's left
% over in each equation, and its norm tells you how accurately the system
% was solved.

format compact ;
A = rand (3)
b = rand (3,1)
x = A\b
r = b-A*x ;
norm (r)

%% BACKSLASH versus INV ... let the battle begin
%
% The backslash operation x=A\b is mathematically the same as x=inv(A)*b.
% However, backslash is faster and more accurate since it uses a matrix
% factorization instead of multiplying by the inverse.  Even though your
% linear algebra textbook might write x=A^(-1)*b as the solution to the
% system A*x=b, your textbook author never means for you to compute the
% inverse.
%
% These next statements give the same answer, so what's the big deal?

S = inv(A) ;
x = S*b
x = A\b

%%
% The big deal is that you should care about speed and you should care even
% more about accuracy.  BACKSLASH relies on matrix factorization (LU, CHOL,
% QR, or other specialized methods).  It's faster and more reliable than
% multiplying by the inverse, particularly for large matrices and sparse
% matrices.  Here's an illustration of how pathetic inv(A)*b can be.

A = gallery ('frank',16) ; xtrue = ones (16,1) ; b = A*xtrue ;

x = inv(A)*b ; norm (b-A*x)
x = A\b      ; norm (b-A*x)

%%
% The performance difference between BACKSLASH and INV for even small
% sparse matrices is striking.

load west0479 ;
A = west0479 ;
n = size (A,1)
b = rand (n,1) ;
tic ; x = A\b ; toc
norm (b-A*x)
tic ; x = inv(A)*b ; toc
norm (b-A*x)

%%
% What if you want to solve multiple systems?  Use a matrix factorization.
% But which one?  And how do you use it?  Here are some alternatives using
% LU for the sparse west0479 matrix, but some are faster than others.

tic ; [L,U]     = lu(A) ; x1 = U \ (L \ b)         ; t1=toc ; nz1=nnz(L+U);
tic ; [L,U,P]   = lu(A) ; x2 = U \ (L \ P*b)       ; t2=toc ; nz2=nnz(L+U);
tic ; [L,U,P,Q] = lu(A) ; x3 = Q * (U \ (L \ P*b)) ; t3=toc ; nz3=nnz(L+U);

fprintf ('1: nnz(L+U): %5d time: %8.4f resid: %e\n', nz1,t1, norm(b-A*x1));
fprintf ('2: nnz(L+U): %5d time: %8.4f resid: %e\n', nz2,t2, norm(b-A*x2));
fprintf ('3: nnz(L+U): %5d time: %8.4f resid: %e\n', nz3,t3, norm(b-A*x3));

%% LU and LINSOLVE are fast and accurate but complicated to use
%
% A quick look at ``help lu'' will scroll off your screen.  For full
% matrices, [L,U,p] = lu (A,'vector') is fastest.  Then for the
% forward/backsolves, use LINSOLVE instead of BACKSLASH for even faster
% performance.  But for sparse matrices, use the optional 'Q' output of LU
% so you get a good fill-reducing ordering.  But you can't use 'Q' if the
% matrix is full.  But LINSOLVE doesn't work on sparse matrices.
%
% But ... Ack!  That's getting complicated ...
%
% Here's the best way to solve A*x=b and A*y=c when A is full and
% unsymmetric:

n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;
tic ; [L,U,p] = lu (A, 'vector') ; LUtime = toc

tic ; x = U \ (L \ b (p,:)) ;
      y = U \ (L \ c (p,:)) ; toc

tic ; opL = struct ('LT', true) ;
      opU = struct ('UT', true) ;
      x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
      y = linsolve (U, linsolve (L, c(p,:), opL), opU) ; toc

%% INV is easy to use, but slow and inaccurate
%
% Oh bother!  Using LU and LINSOLVE is too complicated.  You just want to
% solve your system.  Let's just compute inv(A) and use it twice.  Easy to
% write, but slower and less accurate ...

S = inv (A) ;
x = S*b ; norm (b-A*x)
y = S*c ; norm (c-A*y)

%%
% Sometimes using the inverse seems inevitable.  For example, your textbook
% might show the Schur complement formula as S = A-B*inv(D)*C.  This can be
% done without inv(D) in one of two ways: SLASH or BACKSLASH (MRDIVIDE or
% MLDIVIDE to be precise).
%
% inv(A)*B and A\B are mathematically equivalent, as are B*inv(A) and B/A,
% so these three methods give the same results (ignoring computational
% errors, which are worse for inv(D)).  Only the first equation looks like
% the equation in your textbook, however.

A = rand (200) ; B = rand (200) ; C = rand (200) ; D = rand (200) ;

tic ; S1 = A - B*inv(D)*C ; toc ;
tic ; S2 = A - B*(D\C) ;    toc ;
tic ; S3 = A - (B/D)*C ;    toc ;

%% So the winner is ... nobody
%
% BACKSLASH: mostly simple to use (except remember that Schur complement
%       formula?).  Fast and accurate ... but slow if you want to solve
%       two linear systems with the same matrix A.
%
% LU, QR, CHOL: fast and accurate.  Awful syntax to use.  Drag out your
%       linear algebra textbook if you want to use these in MATLAB.
%       Whenever I use them I have to derive them from scratch, even
%       though I *wrote* most of the sparse factorizations used in MATLAB!
%
% INV: slow and inaccurate.  Wins big on ease-of-use, though, since it's a
%       direct plug-in for all your nice mathematical formulas.
%
% No method is best on all three criterion: speed, accuracy, and ease of
% use.
%
% Is there a solution?  Yes ... keeping reading ...

%% The FACTORIZE object to the rescue
%
% The FACTORIZE method is just as easy to use as INV, but just as fast and
% accurate as BACKSLASH, LU, QR, CHOL, and LINSOLVE.
%
% F = factorize(A) computes the factorization of A and returns it as an
% object that you can reuse to solve a linear system with x=F\b.  It picks
% LU, QR, or Cholesky for you, just like BACKSLASH.
%
% S = inverse(A) is simpler yet.  It does NOT compute inv(A), but
% factorizes A.  When multiplying S*b, it doesn't mulitply by the inverse,
% but uses the correct forward/backsolve equations to solve the linear
% system.

n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;

tic ;                       x = A\b ; y = A\c ; toc
tic ; S = inv(A) ;          x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ;    x = F\b ; y = F\c ; toc
tic ; S = inverse(A) ;      x = S*b ; y = S*c ; toc

%% Least-squares problems
%
% Here are some different methods for solving a least-squares problem when
% your system is over-determined.  The last two methods are the same.

A = rand (1000,200) ;
b = rand (1000,1) ;

tic ; x = A\b            ; toc, norm (A'*A*x-A'*b)
tic ; x = pinv(A)*b      ; toc, norm (A'*A*x-A'*b)
tic ; x = inverse(A)*b   ; toc, norm (A'*A*x-A'*b)
tic ; x = factorize(A)\b ; toc, norm (A'*A*x-A'*b)

%%
% FACTORIZE is better than BACKSLASH because you can reuse the
% factorization for different right-hand-sides.  For full-rank matrices,
% it's better than PINV because it's faster (and PINV fails for sparse
% matrices).

A = rand (1000,200) ;
b = rand (1000,1) ;
c = rand (1000,1) ;

tic ;                  ; x = A\b ; y = A\c ; toc
tic ; S = pinv(A)      ; x = S*b ; y = S*c ; toc
tic ; S = inverse(A)   ; x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ; x = F\b ; y = F\c ; toc

%% Underdetermined systems
%
% The under-determined system A*x=b where A has more columns than rows has
% many solutions.  x=A\b finds a basic solution (some of the entries in x
% are zero).  pinv(A)*b finds a minimum 2-norm solution, but it's slow.  QR
% factorization will do the same if A has full rank.  That's what the
% factorize(A) and inverse(A) methods do.

A = rand (200,1000) ;
b = rand (200,1) ;

tic ; x = A\b            ; toc, norm (x)
tic ; x = pinv(A)*b      ; toc, norm (x)
tic ; x = inverse(A)*b   ; toc, norm (x)
tic ; x = factorize(A)\b ; toc, norm (x)

%% Computing selected entries in the inverse or pseudo-inverse
%
% If you want just a few entries from the inverse, it's still better to
% formulate the problem as a system of linear equations and use a matrix
% factorization instead of computing inv(A).  The FACTORIZE object does
% this for you, by overloading the subsref operator.

A = rand (1000) ;

tic ; S = inv (A)     ; S (2:3,4), toc
tic ; S = inverse (A) ; S (2:3,4), toc

%% Computing the entire inverse or pseudo-inverse
%
% Rarely, and I mean RARELY, you really do need the inverse.  More
% frequently what you want is the pseudo-inverse.  You can force a
% factorization to become a plain matrix by converting it to double.  Note
% that inverse(A) only handles full-rank matrices (either dense or
% sparse), whereas pinv(A) works for all dense matrices (not sparse).
%
% The explicit need for inv(A) (or S=A\eye(n), which is the same thing) is
% RARE.  If you ever find yourself multiplying by the inverse, then you
% know one thing for sure.  You know with certainty that you don't know
% what you're doing.

A = rand (500) ;
tic ; S1 = inv (A) ;            ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)

A = rand (500,400) ;
tic ; S1 = pinv (A)             ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)

%% Update/downdate of a dense Cholesky factorization
%
% Wilkinson considered the update/downdate of a matrix factorization to be
% a key problem in computational linear algebra.  The idea is that you
% first factorize a matrix.  Next, make a low-rank change to A, and patch
% up (or down...) the factorization so that it becomes the factorization of
% the new matrix.  In MATLAB, this only works for dense symmetric positive
% definite matrices, via cholupdate.  This is much faster than computing
% the new factorization from scratch.

n = 1000 ;
A = rand (n) ;
A = A*A' + n*eye (n) ;
w = rand (n,1) ; t = rand (n,1) ; b = rand (n,1) ;
F = factorize (A) ;

tic ; F = cholupdate (F,w,'+') ; x = F\b ; toc
tic ; y = (A+w*w')\b ;      toc
norm (x-y)

tic ; F = cholupdate (F,t,'-') ; x = F\b ; toc
tic ; y = (A+w*w'-t*t')\b ; toc
norm (x-y)

%% Caveat Executor
%
% One caveat:  If you have a large number of very small systems to solve,
% the object-oriented overhead of creating and using an object can dominate
% the run time, at least in MATLAB R2011a.  For this case, if you want the
% best performance, stick with BACKSLASH, or LU and LINSOLVE (just extract
% the appropriate formulas from the M-files in the FACTORIZE package).
%
% Hopefully the object-oriented overhead will drop in future versions of
% MATLAB, and you can ignore this caveat.

A = rand (10) ; b = rand (10,1) ; F = factorize (A) ;

tic ; for k = 1:10000, x = F\b ; end ; toc

tic ; for k = 1:10000, x = A\b ; end ; toc

[L,U,p] = lu (A, 'vector') ;
opL = struct ('LT', true) ;
opU = struct ('UT', true) ;
tic ; 
for k = 1:10000
    x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
end
toc

%% Summary
%
% So ... don't use INV, and don't worry about how to use LU, CHOL, or QR
% factorization.  Just install the FACTORIZE package, and you're on your
% way.  Assuming you are now in the Factorize/ directory, cut-and-paste
% these commands into your command window:
%
%   addpath (pwd)
%   savepath
%
% And remember ...
%
% "Don't let that INV go past your eyes; to solve that system, FACTORIZE!"


##### SOURCE END #####
--></body></html>