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% [?]
Related Posts
Comments
Translator
Follow Us
Video Posts
