%% THE FACTORIZE OBJECT for solving linear systems
%
% Copyright 2011, Timothy A. Davis, University of Florida.
% davis@cise.ufl.edu
% http://www.cise.ufl.edu/~davis
%
% 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 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!"
