Required R packages¶

Computer methods for solving linear systems of equations¶

We consider computer algorithms for solving linear equations $\mathbf{A} \mathbf{x} = \mathbf{b}$, a ubiquitous task in statistics.

Our mantra¶

The form of a mathematical expression and the way the expression should be evaluated in actual practice may be quite different.

-- James Gentle, Matrix Algebra, Springer, New York (2007).

Algorithms¶

• Algorithm is loosely defined as a set of instructions for doing something, which terminates in finite time. An algorithm requires input and output.

Measure of algorithm efficiency¶

• A basic unit for measuring algorithmic efficiency is flop.

A flop (floating point operation) consists of a floating point addition, subtraction, multiplication, division, or comparison, and the usually accompanying fetch and store.

• How to measure efficiency of an algorithm? Big O notation. If $n$ is the size of a problem, an algorithm has order $O(f(n))$, where the leading term in the number of flops is $c \cdot f(n)$. For example,

  • matrix-vector multiplication A * b, where A is $m \times n$ and b is $n \times 1$, takes $2mn$ or $O(mn)$ flops
  • matrix-matrix multiplication A * B, where A is $m \times n$ and B is $n \times p$, takes $2mnp$ or $O(mnp)$ flops

• Typical orders of computational complexity:

  • Exponential order: $O(b^n)$ ("horrible")
  • Polynomial order: $O(n^q)$ (doable)
  • $O(n \log n )$ (fast)
  • Linear order $O(n)$ (fast)
  • Logarithmic order $O(\log n)$ (super fast)

• Difference of $O(n^2)$ and $O(n\log n)$ on massive data. Suppose we have a teraflops supercomputer capable of doing $10^{12}$ flops per second. For a problem of size $n=10^{12}$, $O(n \log n)$ algorithm takes about

$O(n^2)$ algorithm takes about $10^{12}$ seconds, which is approximately 31710 years!

• The Quicksort and Fast Fourier Transform (invented by John Tukey--recall "bit"--) are celebrated algorithms that turn $O(n^2)$ operations into $O(n \log n)$.

Triangular systems¶

Idea: turning original problem into an easy one, e.g., triangular system.

Lower triangular system¶

To solve $\mathbf{A} \mathbf{x} = \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{n \times n}$ is lower triangular

• Forward substitution:

• $1 + 3 + 5 + \cdots + (2n-1) = n^2$ flops.

Upper triangular system¶

To solve $\mathbf{A} \mathbf{x} = \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{n \times n}$ is upper triangular
$$ \begin{bmatrix} a_{11} & \cdots & a_{1,n-1} & a_{1n} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & a_{n-1,n-1} & a_{n-1,n} \\ 0 & 0 & 0 & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix} = \begin{bmatrix} b_1 \\ \vdots \\ b_{n-1} \\ b_n \end{bmatrix}. $$

• Back substitution (backsolve):

• $n^2$ flops.

Implementation¶

• Most numerical computing software packages (including R) are built on top of the BLAS (basic linear algebra subprograms).

  • Netlib provides a reference implementation.
  • Vanila R (the binary that you can download from https://www.r-project.org) uses its own implementation of BLAS (see sessionInfo() above)

• The BLAS triangular system solve is done in place, i.e., $\mathbf{b}$ is overwritten by $\mathbf{x}$.

# forward substitution
    for (i in seq_len(n)) {
        for (j in seq_len(i-1)) {
            b[i] <- b[i] - A[i, j] * b[j]
        }
        b[i] <- b[i] / A[i, i]
    }
    # backsolve
    for (i in rev(seq_len(n))) {
        for (j in i + seq_len(max(n - i, 0))) {
            b[i] <- b[i] - A[i, j] * b[j]
        }
        b[i] <- b[i] / A[i, i]
    }

• R
  • ?backsolve, ?forwardsolve
  • This is a wrapper for the level-3 BLAS routine dtrsm
  • The d stands for double precision.

Some algebraic facts of triangular system¶

• Eigenvalues of a triangular matrix $\mathbf{A}$ are diagonal entries $\lambda_i = a_{ii}$.

• Determinant $\det(\mathbf{A}) = \prod_i a_{ii}$.

• The product of two upper (lower) triangular matrices is upper (lower) triangular.

• The inverse of an upper (lower) triangular matrix is upper (lower) triangular.

• The product of two unit upper (lower) triangular matrices is unit upper (lower) triangular.

• The inverse of a unit upper (lower) triangular matrix is unit upper (lower) triangular.

Gaussian Elimination and LU Decomposition¶

• A toy example.

• History: the method is named after Carl Friedrich Gauss (1777–1855), although it stems from the notes of Isaac Newton. If fact, GE was known to Chinese mathematicians as early as 179 AD.

What happens when we call solve(A, b)?

Elementary operator matrix¶

• Elementary operator matrix is the identity matrix with the 0 in position $(j,k)$ replaced by $c$

• $\mathbf{E}_{jk}(c)$ is unit triangular, full rank, and its inverse is

• $\mathbf{E}_{jk}(c)$ left-multiplies an $n \times m$ matrix $\mathbf{X}$ effectively replace its $j$-th row $\mathbf{x}_{j\cdot}$ by $c \mathbf{x}_{k \cdot} + \mathbf{x}_{j \cdot}$

$2m$ flops (why?).

• Gaussian elimination can be understood as applying a sequence of elementary operator matrices to transform the linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$ to an upper triangular system

Let's get back to our toy example.

Column 1¶

Column 2¶

Thus we have arrived at an upper triangular matrix $\mathbf{U}$.

Gauss transformations¶

• We can combine elementary operations by column.

• For column 1,

• For column $k$,

• Matrices $\mathbf{M}_1, \ldots, \mathbf{M}_{n-1}$ are called the Gauss transformations.

• Gauss transformations $\mathbf{M}_k$ are unit triangular, full rank, with inverse

LU decomposition¶

Gaussian elimination does $$ \mathbf{M}_{n-1} \cdots \mathbf{M}_1 \mathbf{A} = \mathbf{U}. $$
Let

Thus we have the LU decomposition $$ \mathbf{A} = \mathbf{L} \mathbf{U}, $$ where $\mathbf{L}$ is unit lower triangular and $\mathbf{U}$ is upper triangular.

• LU decomposition exists if the principal sub-matrix $\mathbf{A}[1:k, 1:k]$ is non-singular for $k=1,\ldots,n-1$.

• If the LU decomposition exists and $\mathbf{A}$ is non-singular, then the LU decmpositon is unique and

- This is basically how R (and MATLAB and Julia) compute determinants.

• The whole LU algorithm is done in place, i.e., $\mathbf{A}$ is overwritten by $\mathbf{L}$ and $\mathbf{U}$.

for (k in seq_len(n-1)) {
    for (i in k +seq_len(max(n - k, 0))) {
        A[i, k] <- A[i, k] / A[k, k]   # computes L
        for (j in k +seq_len(max(n - k, 0))) {
            A[i, j] <- A[i, j] - A[i, k] * A[k, j]  # computes U
        }
    }
}

Source: http://www.netlib.org

• The LU decomposition costs (ignoring the cost of division)

• Given LU decomposition $\mathbf{A} = \mathbf{L} \mathbf{U}$, solving $(\mathbf{L} \mathbf{U}) \mathbf{x} = \mathbf{b}$ for one right hand side costs $2n^2$ flops:

  • One forward substitution ($n^2$ flops) to solve

  $$ \mathbf{L} \mathbf{y} = \mathbf{b} $$

  • One back substitution ($n^2$ flops) to solve

  $$ \mathbf{U} \mathbf{x} = \mathbf{y} $$

• Therefore to solve $\mathbf{A} \mathbf{x} = \mathbf{b}$ via LU, we need a total of

• If there are multiple right hand sides, LU only needs to be done once.

Matrix inversion¶

• For matrix inversion, there are $n$ right hand sides $\mathbf{e}_1, \ldots, \mathbf{e}_n$. Taking advantage of 0s reduces $2n^3$ flops to $\frac 43 n^3$ flops. So matrix inversion via LU costs

• No inversion mentality:

  Whenever we see matrix inverse, we should think in terms of solving linear equations. In short, no solve(A) %*% b, but solve(A, b).

  We do not compute matrix inverse unless

  1. it is necessary to compute standard errors
  2. number of right hand sides is much larger than $n$
  3. $n$ is small

Pivoting¶

• Let

Is there a solution to $\mathbf{A} \mathbf{x} = \mathbf{b}$ for an arbitrary $\mathbf{b}$? Does GE/LU work for $\mathbf{A}$?

• What if, during LU procedure, the pivot $a_{kk}$ is 0 or nearly 0 due to underflow?
  Example:

with solution

• With two significant digits, GE via LU yields

or
$$ \begin{split} 0.001x_1 + x_2 &= 1 \\ -1000 x_2 &= -1000 \end{split} , $$ yielding $$ \begin{bmatrix} x_1 \\ x_ 2 \end{bmatrix} = \begin{bmatrix} 0.0 \\ 1.0 \end{bmatrix} ! $$

Condition number and stability¶

• The problem of solving $\mathbf{A}\mathbf{x}=\mathbf{b}$ for fixed $b$ can be viewed as a map $f: \mathbb{R}^{n\times n} \ni \mathbf{A} \mapsto \mathbf{A}^{-1}\mathbf{b} \in \mathbb{R}^n$.

• A well-conditioned problem is one such that all small perturbations of a given input $\mathbf{A}$ lead to only small changes in $f(\mathbf{A})$.

  • An ill-conditioned problem is one such that some small perturbation of $\mathbf{A}$ leads to a large change in $f(\mathbf{A})$.

• The conditioning of the problem of computing $\mathbf{A}^{-1}\mathbf{b}$ given $\mathbf{A}$ is usually measured by the condition number

+ $\sigma_{\max}(\mathbf{A}) = \max_{\mathbf{x}\neq 0} \|\mathbf{A}\mathbf{x}\|/\|\mathbf{x}\|$ is the maximum singular value of $\mathbf{A}$;
+ $\sigma_{\min}(\mathbf{A}) = \min_{\mathbf{x}\neq 0} \|\mathbf{A}\mathbf{x}\|/\|\mathbf{x}\|$ is the minimum singular value of $\mathbf{A}$.

• If $\mathbf{A}$ is singular, then $\sigma_{\min}(\mathbf{A}) = 0$, hence $\kappa(\mathbf{A}) = \infty$. So a large condition number indicates that $\mathbf{A}$ is close to singular and the problem is inherently hard.

Stability¶

• An algorithm for computing $\mathbf{A}^{-1}\mathbf{b}$ using floating-point arithmetric can be viewed as another map $\tilde{f}: \mathbb{R}^{n\times n} \to \mathbb{R}^{n}$.

• For given input $\mathbf{A}$, the solution computed by algorithm $\tilde{f}$ is $\hat{\mathbf{x}} = \tilde{f}(\mathbf{A})$.

  • Example 1: solve $\mathbf{A}\mathbf{x}=\mathbf{b}$ by GE/LU followed by forward and backward substitutions.
  • Example 2: solve $\mathbf{A}\mathbf{x}=\mathbf{b}$ by Gauss-Seidel (an iterative method to come).

• Algorithms will be affected by at least rounding errors.

• Algorithm $\tilde{f}$ is called stable if for any perturbed input $\tilde{\mathbf{A}}$ from $\mathbf{A}$ with relative error bounded by $\epsilon$, the relative error between $\hat{\mathbf{x}}$ and $\tilde{\mathbf{A}}^{-1}\mathbf{b}$ is bounded by $O(\epsilon)$.

Analysis of the example¶

• The condition number of matrix $\mathbf{A} = \begin{bmatrix} 0.001 & 1 \\ 1 & 1 \end{bmatrix}$ is 2.262, so the problem itself is well-conditioned: a small change in $\mathbf{A}$ won't yield a large change in the solution $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$.

• It is the instability of the GE/LU algorithm that yields the large error between the true solution $\mathbf{x}\approx(1, 1)$ and the computed solution (or the output of the algorithm) $\hat{\mathbf{x}}=(0, 1)$. This example suggests that with GE/LU, the worst-case relative error is $\|\mathbf{x}-\hat{\mathbf{x}}\| / \|\mathbf{x}\| \ge 1/\sqrt{2}$ (zero relative error in input $\tilde{\mathbf{A}}$), not small.

• The root cause is the vastly different scale of $a_{11}$ with other coefficients and the order of ellimination. In particular, any $a_{22}$ such that $[a_{22} - 1000] = -1000$ and $b_2$ such that $[b_2 - 1000] = 1000$ will yield the same solution $\hat{\mathbf{x}}=(0, 1)$.

That is, the $1000 = a_{21}/a_{11}$ causes the loss of information in $a_{22}$ and $b_2$. (Recall Scenario 1 of Catastrophic Cancellation in Lecture 1.)

• Indeed, the LU decomposition of $\mathbf{A}$ within the precision is

resulting in $\mathbf{L}\mathbf{U} = \begin{bmatrix} 0.001 & 1 \\ 1 & 0 \end{bmatrix}$, so $\mathbf{L}\mathbf{U}$ is not even closed to $\mathbf{A} = \begin{bmatrix} 0.001 & 1 \\ 1 & 1 \end{bmatrix}$!

Pivoting (cont'd)¶

• Solution: fix the algorithm by pivoting (i.e., change the order of elimination)

• With two significant digits, GE yields

or $$ \begin{split} x_1 + x_2 &= 2 \\ 1.0 x_2 &= 1.0 \end{split} , $$ yielding $$ \begin{bmatrix} x_1 \\ x_ 2 \end{bmatrix} = \begin{bmatrix} 1.0 \\ 1.0 \end{bmatrix} . $$

• Note that $0.001=a_{11}/a_{22}$ does not cause a loss of information in $a_{12}$ and $b_1$. The LU decomposition within the precision is

• Partial pivoting: before zeroing the $k$-th column, move the row with $\max_{i=k}^n |a_{ik}|$ is into the $k$-th row (called GEPP).

  • Cost: $(n-1) + (n-2) + \dotsb + 1 = O(n^2)$.

• LU with partial pivoting yields

where $\mathbf{P}$ is a permutation matrix, $\mathbf{L}$ is unit lower triangular with $|\ell_{ij}| \le 1$, and $\mathbf{U}$ is upper triangular.

• Complete pivoting (GECP): do both row and column interchanges so that the largest entry in the sub matrix A[k:n, k:n] is permuted to the $(k,k)$-th entry. This yields the decomposition

where $|\ell_{ij}| \le 1$.

• GEPP is the most commonly used methods for solving general linear systems. GECP is more stable but costs more computation. Partial pivoting is stable most of times (which is partially because GEPP guarantees $|\ell_{ij}| \le 1$).

Implementation¶

• Except for the iterative methods, most of these numerical linear algebra tasks are implemented in the BLAS and LAPACK libraries. LAPACK is built on top of BLAS. They form the building blocks of most statistical computing tasks.

LAPACK (Linear Algebra PACKage) is a standard software library for numerical linear algebra. It provides routines for solving systems of linear equations and linear least squares, eigenvalue problems, and singular value decomposition. -Wikipedia

• LAPACK routine GETRF does $\mathbf{P} \mathbf{A} = \mathbf{L} \mathbf{U}$ (LU decomposition with partial pivoting) in place.

• R: solve() implicitly performs LU decomposition (wrapper of LAPACK routine DGESV). solve() allows specifying a single or multiple right hand sides. If none, it computes the matrix inverse. The Matrix package contains lu() function that outputs L, U, and pivot.