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;
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.

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;


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;
Dxx = spdiags([e -2*e e], -1:1, N+1, N+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:


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:

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