Writing Fast Matlab Code #8: Solving a linear system (Ax = b)
Many problems involve solving a linear system, that is, solving the Ax = b problem
The condition number cond(A) describes how sensitive the problem is to small perturbations. Matlab can efficiently approximate cond(A) with condest. Smaller condition number is better, and in this case A is said to be well-conditioned. An infinite condition number implies that A is singular—the problem cannot be solved, though it may be solved approximately with the pseudoinverse :
x = pinv(A)*b
Small- to moderately-sized problems are solved efficiently by the backslash operator:
x = A\b; % Solves A*x = bWe shall focus on large problems.
Iterative Methods
For large problems, Matlab is well-equipped with iterative solvers: bicg, bicgstab, cgs, gmres, lsqr,
minres, pcg, symmlq, qmr. The general calling syntax (except for gmres) is
[x,flag,relres,iter,resvec] = method (A,b,tol,maxit,M1,M2,x0)
where the required arguments are the matrix A, the right-hand side b, and the solution x. The minres, pcg, and symmlq methods require A to be symmetric (A = A’) and pcg additionally requires positive definiteness (all eigenvalues are positive). The optional input variables are
- tol accuracy tolerance (default 10−6 )
- maxit maximum number of iterations (default 20)
- M1, M2 the method uses preconditioner M = M1*M2; it solves M −1 Ax = M −1 b
- x0 initial guess
The optional output variables are
- flag convergence flag (0 = success, nonzero values are error codes)
- relres relative residual kb − Axk/kbk
- iter number of iterations
- resvec a vector listing the residual norm kb − Axk at each iteration
The gmres method includes one additional input,
[x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,...)
The method restarts every restart iterations, or if restart=[], the method does not restart.
Which method works best depends on the particular problem. This diagram (adapted from Demmel [3])
provides a reasonable starting point.
Functional Representation
Rather than storing the matrix explicitly, it is often convenient to represent A as a linear function acting on x. For example,
can be expressed as A*x = cumsum(x). In addition to being conceptually convenient, functional repre- sentation avoids the memory overhead of storing the matrix.
All the iterative solvers support functional representation, where a function handle afun is given in place of matrix A. The function afun(x) should accept a vector x and return A*x. The help for lsqr (“help lsqr”) has a functional example where A is the central differencing operator. The bicg, lsqr, and qmr methods also require the transpose: afun(x,’notransp’) should return A*x and afun(x,’transp’) should return A’*x. Similarly, the preconditioner may be given as a function handle mfun.
Special Matrices
For certain special matrix structures, the Ax = b problem can be solved very quickly.
- Circulant matrices. Matrices corresponding to cyclic convolution Ax = h ∗ x are diagonalized in the Fourier domain, and can be solved by x = ifft(fft(b)./fft(h));
- Triangular and banded. Triangular matrices and diagonally-dominant banded matrices are solved efficiently by sparse LU factorization: [L,U] = lu(sparse(A)); x = U\(L\b);
- Poisson problem. If A is a finite difference approximation of the Laplacian, the problem is effi- ciently solved by multigrid methods (e.g., web search for “matlab multigrid”).
Inlined PCG
As discussed in Section 6, it is sometimes advantageous to inline a simple piece of code to tailor it to a specific application. While most iterative methods are too lengthy for inlining, conjugate gradients has a short implementation. The value of inlining pcg is that calls to afun and mfun may be inlined. Beware that in the following code, A must be symmetric positive definite—otherwise, it will not work!
%% Use PCG to solve Ax = b %% tol = 1e−6; % Convergence tolerance maxit = 20; % Maximum number of iterations x = x0; % Set x to initial guess (if no guess is available, use x0 = 0) Compute w = A x r = b − w; % r = residual vector Set p = r, or if using a preconditioner, p = M −1 r d = r'*p; for iter = 1:maxit if norm(p) < tol break; % PCG converged successfully (small change in x) end Compute w = A p alpha = e/(p'*w); x = x + alpha*p; r = r − alpha*w; Set w = r, or if using a preconditioner, w = M −1 r if norm(w) < tol && norm(r) < tol break; % PCG converged successfully (small residual) end end dlast = d; d = r'*w; p = w + (d/dlast)*p;
After the loop, the final solution is in x. The method converged successfully with accuracy tol if it encountered a break statement, otherwise, it failed or needed more iterations. The relative residual may be computed by relres = norm(r)/norm(b) after the loop and resvec by storing the value of norm(r) each iteration.
Popularity: 1% [?]
Linear Systems Solver (for small Systems)
SolverNxN is designed to solve many small linear system in a vectorized way in order to improve time perfomaces. It is in practice a faster alternative to the matlab commands:
for k = 1:number_of_systems</span></span>
x=A\b;
end
The Matlab looping speed can visibly slow down this process so most of the time is spent in the for loop and not to solve the system. This is particulary true when A is a small matrix. In this case SolverNxN can be a useful tool.
The speed improvement is more visibile the greater number are the systems and the smaller is their dimension.
The algorithm works assembling all the small elementary systems into a sparse matrix so the \ command is called just one time.
VERY IMPORTANT:
Because of algorithm structure is very important that no system is singular. One singular system can ruin the whole solution. So use this routine only if you are sure all the systems are solvable.
Let’s say we want to solve Nsyst m x m linear systems.
The function call is:
x=SolverNxN(Amxn,b);
INPUTS:
- Amxn: is a matrix of dimensions [m ,(m x Nsyst)], so it is just an horizontal concatenation of the small elementary systems.
- b: is a [(m x Nsyst),1] vector formed by the vertical concatenation of the elematry system vectors.
OUTPUT:
- x: is a [(m x Nsyst),1] vector formed by the vertical concatenation of the elematry system solutions.
Here it is an example:
It is inlcuded a Test for speed comparison, here it is some results:
Program started We are going to solve 1000000 3x3 systems Looping took: 19.6869 s Solver NxN took: 3.1457 s
Residue is 0.0000 Solution is correct!
For question,suggestion,bugs report: giaccariluigi@msn.com
It is brand new so please report problems to my e-mail
Author: Giaccari Luigi
Created: 12/02/2009
Last Update: 12/02/2009
Download Now
| for i=1:number_of_systems | |
| x=A\b; | |
| end | |
Popularity: 6% [?]
Translator
Follow Us
Video Posts
