Solving Linear Systems
This script is a tutorial for how to solve linear systems using both direct and iterative methods in Matlab.
A linear system of equations, of the form Ax=b, can be solved both directly and iteratively. Matlab has a number of different methods built into the function "\". To solve the system Ax=b you would call "A\b". This is by far the quickest way to directly solve a linear system since Matlab chooses the best method depending on the classification of the matrix.
N = 10; A1=rand(N); b1=rand(N,1); tic x=A1\b1; toc
Elapsed time is 0.014921 seconds.
The other way of solving a linear system is by using iterative methods. This minimizes the residual or error to within a specified tolerance, solving the system. We can use the built in function "bicgstab" to use Biconjugate gradients stabilized method.
tic bicgstab(A1,b1,1e-5,1000); toc
bicgstab converged at iteration 17.5 to a solution with relative residual 7e-06. Elapsed time is 0.056772 seconds.
Consider now if we change N to be much larger.
N = 1000; A2=rand(N); b2=rand(N,1); tic x=A2\b2; toc tic bicgstab(A2,b2,1e-6,1000); toc
Elapsed time is 0.034424 seconds. bicgstab stopped at iteration 195 without converging to the desired tolerance 1e-06 because a scalar quantity became too small or too large to continue computing. The iterate returned (number 1) has relative residual 0.5. Elapsed time is 0.211840 seconds.
As you can see, for large dense matrices, the iterative method does not perform anywhere close to the matrix left divide function (in fact, in many cases it does not even converge). However, let's now consider a sparse matrix.
xmin = 0; xmax = 50; N = 2^10; x = linspace(xmin,xmax,N+1); dx=x(2)-x(1); dx2=dx*dx; e=ones(N+1,1); Dxx = spdiags([e -2*e e], -1:1, N+1, N+1); Dxx=(1/dx2)*Dxx; b3=ones(size(Dxx,1),1);
This is the second order differentiation matrix from an earlier script. It is sparse with only entries along the diagonal. Now let's try again:
tic x=Dxx\b3; toc tic bicgstab(Dxx,b3,1e-6,1000); toc
Elapsed time is 0.000257 seconds. bicgstab converged at iteration 600.5 to a solution with relative residual 5.3e-07. Elapsed time is 0.061603 seconds.
In general iterative methods work better on sparse matrices and ones with diagonal dominance. In specific cases and for particular problems iterative methods will be better, but in most cases you should use the built in "backslash" solve since it has been optimized over decades to use the best method for a given matrix.
If you want to see exactly what "\" is doing under the hood, you can change the verbosity so that it prints its choices:
spparms('spumoni',2) A2\rand(size(A2,1),1); Dxx\b3; sparse(rand(N).*round(rand(N)-0.2))\rand(N,1);
sp\: bandwidth = 1+1+1. sp\: is A diagonal? no. sp\: is it ok to solve real tridiagonal A without pivoting? yes. sp\: bandwidth = 1023+1+1022. sp\: is A diagonal? no. sp\: is band density (0.299638) > bandden (0.5) to try banded solver? no. sp\: is A triangular? no. sp\: is A morally triangular? no. sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? no. sp\: use Unsymmetric MultiFrontal PACKagewith automatic reordering. UMFPACK V5.4.0 (May 20, 2009), Control: Matrix entry defined as: double Int (generic integer) defined as: UF_long 0: print level: 2 1: dense row parameter: 0.2 "dense" rows have > max (16, (0.2)*16*sqrt(n_col) entries) 2: dense column parameter: 0.2 "dense" columns have > max (16, (0.2)*16*sqrt(n_row) entries) 3: pivot tolerance: 0.1 4: block size for dense matrix kernels: 32 5: strategy: 0 (auto) 6: initial allocation ratio: 0.7 7: max iterative refinement steps: 2 12: 2-by-2 pivot tolerance: 0.01 13: Q fixed during numerical factorization: 0 (auto) 14: AMD dense row/col parameter: 10 "dense" rows/columns have > max (16, (10)*sqrt(n)) entries Only used if the AMD ordering is used. 15: diagonal pivot tolerance: 0.001 Only used if diagonal pivoting is attempted. 16: scaling: 1 (divide each row by sum of abs. values in each row) 17: frontal matrix allocation ratio: 0.5 18: drop tolerance: 0 19: AMD and COLAMD aggressive absorption: 1 (yes) The following options can only be changed at compile-time: 8: BLAS library used: Fortran BLAS. size of BLAS integer: 8 9: compiled for MATLAB 10: CPU timer is POSIX times ( ) routine. 11: compiled for normal operation (debugging disabled) computer/operating system: unknown size of int: 4 UF_long: 8 Int: 8 pointer: 8 double: 8 Entry: 8 (in bytes) sp\: UMFPACK's factorization was successful. sp\: UMFPACK's solve was successful. UMFPACK V5.4.0 (May 20, 2009), Info: matrix entry defined as: double Int (generic integer) defined as: UF_long BLAS library used: Fortran BLAS. size of BLAS integer: 8 MATLAB: yes. CPU timer: POSIX times ( ) routine. number of rows in matrix A: 1024 number of columns in matrix A: 1024 entries in matrix A: 314193 memory usage reported in: 16-byte Units size of int: 4 bytes size of UF_long: 8 bytes size of pointer: 8 bytes size of numerical entry: 8 bytes strategy used: unsymmetric ordering used: colamd on A modify Q during factorization: yes prefer diagonal pivoting: no pivots with zero Markowitz cost: 0 submatrix S after removing zero-cost pivots: number of "dense" rows: 1024 number of "dense" columns: 1024 number of empty rows: 0 number of empty columns 0 submatrix S square and diagonal preserved pattern of square submatrix S: number rows and columns 1024 symmetry of nonzero pattern: 0.297084 nz in S+S' (excl. diagonal): 534552 nz on diagonal of matrix S: 289 fraction of nz on diagonal: 0.282227 symbolic factorization defragmentations: 0 symbolic memory usage (Units): 703584 symbolic memory usage (MBytes): 10.7 Symbolic size (Units): 2629 Symbolic size (MBytes): 0 symbolic factorization CPU time (sec): 0.02 symbolic factorization wallclock time(sec): 0.02 matrix scaled: yes (divided each row by sum of abs values in each row) minimum sum (abs (rows of A)): 1.28869e+02 maximum sum (abs (rows of A)): 1.82775e+02 symbolic/numeric factorization: upper bound actual % variable-sized part of Numeric object: initial size (Units) 642729 641704 100% peak size (Units) 1731846 1504811 87% final size (Units) 540112 520516 96% Numeric final size (Units) 545780 525672 96% Numeric final size (MBytes) 8.3 8.0 96% peak memory usage (Units) 1749383 1522348 87% peak memory usage (MBytes) 26.7 23.2 87% numeric factorization flops 7.10325e+08 6.96778e+08 98% nz in L (incl diagonal) 522362 519916 100% nz in U (incl diagonal) 524800 518172 99% nz in L+U (incl diagonal) 1046138 1037064 99% largest front (# entries) 1029120 984966 96% largest # rows in front 1005 1004 100% largest # columns in front 1024 983 96% initial allocation ratio used: 0.7 # of forced updates due to frontal growth: 3 nz in L (incl diagonal), if none dropped 519916 nz in U (incl diagonal), if none dropped 518172 number of small entries dropped 0 nonzeros on diagonal of U: 1024 min abs. value on diagonal of U: 1.16e-03 max abs. value on diagonal of U: 2.45e-01 estimate of reciprocal of condition number: 4.74e-03 indices in compressed pattern: 2607 numerical values stored in Numeric object: 1038406 numeric factorization defragmentations: 2 numeric factorization reallocations: 2 costly numeric factorization reallocations: 2 numeric factorization CPU time (sec): 0.18 numeric factorization wallclock time (sec): 0.14 numeric factorization mflops (CPU time): 3870.99 numeric factorization mflops (wallclock): 4976.98 symbolic + numeric CPU time (sec): 0.20 symbolic + numeric mflops (CPU time): 3483.89 symbolic + numeric wall clock time (sec): 0.16 symbolic + numeric mflops (wall clock): 4354.86 solve flops: 6.99340e+06 iterative refinement steps taken: 1 iterative refinement steps attempted: 1 sparse backward error omega1: 1.60e-16 sparse backward error omega2: 0.00e+00 solve CPU time (sec): 0.01 solve wall clock time (sec): 0.00 solve mflops (CPU time): 699.34 total symbolic + numeric + solve flops: 7.03771e+08 total symbolic + numeric + solve CPU time: 0.21 total symbolic + numeric + solve mflops (CPU): 3351.29