Today's Topics

• Conjugate gradient method

• Eigenvalue computation

Direct Methods

• All the methods introduced so far can be categorized as direct methods to solve linear systems

• Meaning that the exact solution to $Ax=b$ can be computed within a finite number of operations (in exact arithmetic)

• Cons: too expensive to handle very large linear systems

Iterative Methods

• Computes a sequence of approximate solutions $x^{\left(k\right)}$ that converges to the exact solution

• Stops when the precision is sufficient

• Do not request more precision than is needed

• In each iteration, the computation is typically cheap, e.g., matrix-vector multiplication $v\to Av$

• Especially efficient for sparse matrices

Conjugate Gradient Method

• Conjugate gradient method (CG) is a very special linear system solver

• It is a direct method used as an iterative one

• Aims to solve the linear system $Ax=b$ when $A$ is positive definite (p.d.)

Conjugate Gradient Method

• Only uses the matrix-vector multiplication $v\to Av$

• Obtains the exact solution after $n$ steps

• Converges fast under some conditions

Algorithm ^[2]

Convergence ^[2]

In exact arithmetic, the algorithm converges to the solution of the linear system $Ax=b$ in at most $n$ iterations.

Convergence ^[2]

Let $A$ be a p.d. matrix and denote by $x^{\ast }$ the solution to $Ax=b$. Let $x^{\left(k\right)}$ be the sequence of approximate solutions produced by the conjugate gradient method. Then

$$\parallel x^{\left(k\right)}-x^{\ast }{\parallel }_2\le 2\sqrt{\kappa }{\left(\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\right)}^k\parallel x_0-x^{\ast }{\parallel }_2,$$ where $\kappa$ is the condition number of $A$ (defined later).

Condition Number

• For any matrix $A$, the condition number $\kappa \left(A\right)$ is defined as $\kappa \left(A\right)={\mu }_{max}\left(A\right)/{\mu }_{min}\left(A\right)$, where ${\mu }_{max}$ and ${\mu }_{min}$ are the largest and smallest singular values, respectively

• For a p.d. matrix $A$, $\kappa \left(A\right)={\lambda }_{max}\left(A\right)/{\lambda }_{min}\left(A\right)$, where ${\lambda }_{max}$ and ${\lambda }_{min}$ are the largest and smallest eigenvalues, respectively

• $\kappa \left(A\right)\ge 1$

• For any orthogonal matrix $Q$, $\kappa \left(Q\right)=1$ and $\kappa \left(QA\right)=\kappa \left(AQ\right)=\kappa \left(A\right)$

Insights

• We say a matrix $A$ is well-conditioned if $\kappa \left(A\right)$ is small (i.e., $\kappa \left(A\right)\approx 1$)
• Otherwise we say $A$ is ill-conditioned (i.e., $\kappa \left(A\right)\gg 1$)

• If a p.d. matrix $A$ has almost equal eigenvalues, then $\kappa \left(A\right)\approx 1$ and $$\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\approx 0.$$ CG will have extremely fast convergence

• If $A$ almost loses positive definiteness, ${\lambda }_{min}\left(A\right)\approx 0$, then $\kappa \left(A\right)\gg 1$, and CG almost has no progress

Implementation

• A direct "translation" of the algorithm gives

cg = function(A, b, x0 = rep(0, length(b)), eps = 1e-6)
{
    m = length(b)
    x = x0
    p = r = b - A %*% x
    r2 = sum(r^2)
    errs = c()
    for(i in 1:m)
    {
        Ap = A %*% p
        alpha = r2 / sum(p * Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        r2_new = sum(r^2)
        err = sqrt(r2_new)
        errs = c(errs, err)
        if(err < eps)
            break
        beta = r2_new / r2
        p = r + beta * p
        r2 = r2_new
    }
    list(x = x, errs = errs, niter = i)
}

Implementation

• We test on a simulated matrix

• However, the algorithm does not seem to converge, and the error is quite large

set.seed(123)
n = 100
M = matrix(0.1 * rnorm(n^2), n)
A = crossprod(M)
b = rnorm(n)
sol = cg(A, b, eps = 1e-12)
sol$niter

## [1] 100

max(abs(A %*% sol$x - b))

## [1] 0.5763024

Implementation

• We can find that $A$ is ill-conditioned

kappa(A)

## [1] 1164091

• Adding a positive number to the diagonal results in much smaller condition number

kappa(A1 <- A + diag(rep(0.1, n)))

## [1] 134.8415

kappa(A2 <- A + diag(rep(1, n)))

## [1] 15.01074

Implementation

sol1 = cg(A1, b, eps = 1e-12)
sol1$niter

## [1] 78

max(abs(A1 %*% sol1$x - b))

## [1] 1.588243e-13

sol2 = cg(A2, b, eps = 1e-12)
sol2$niter

## [1] 29

max(abs(A2 %*% sol2$x - b))

## [1] 1.858513e-13

Implementation

• Plot the residual norms (in log-scale) against iterations

Implementation

• We will see a more useful example when we consider linear regression on sparse matrices

When to Use

• $A$ is p.d. and is well-conditioned

• The matrix-vector operation $v\to Av$ is cheap

• One example is sparse matrix, or product of sparse matrices

• We will elaborate this point when we study sparse matrices

Eigenvalues

• Eigenvalues are important characteristics of matrices

• Wide applications in statistics, e.g., principal component analysis (PCA), spectral clustering, matrix norm, random matrix theory, etc.

• Eigenvalue computation is one of the core topics of numerical linear algebra

Definition

• The definition of eigenvalue is actually quite simple

• For a square matrix $A$, we call the solution $\lambda$ to the equation $Ax=\lambda x$ an eigenvalue of $A$

• The corresponding $x$ vector is an eigenvector

• Here we focus on real symmetric matrices

Decomposition Theorem ^[2]

A matrix $A_{n\times n}$ is real symmetric if and only if there exists a real orthogonal matrix $\Gamma$ and real eigenvalues ${\lambda }_1,\dots ,{\lambda }_n$ such that $$A=\Gamma D{\Gamma }^{\prime },$$ where $D=diag({\lambda }_1,\dots ,{\lambda }_n)$.

This is typically called the eigenvalue decomposition, or spectral decomposition, of $A$.

Basic Properties

Suppose that $A$ is a real symmetric matrix

• $A$ is p.d. if all its eigenvalues are strictly positive

• $A^k=\Gamma D^k{\Gamma }^{\prime }$, $k=\dots ,-2,-1,0,1,2,\dots$

• $det\left(A\right)={\lambda }_1\cdots {\lambda }_n$

• $tr\left(A\right)={\lambda }_1+\cdots +{\lambda }_n$

Basic Properties

• Min-Max principle $${\lambda }_{min}\left(A\right)=\underset{x\ne 0}{min}\frac{x^{\prime }Ax}{x^{\prime }x}=\underset{\parallel x\parallel =1}{min}{x}^{\prime }Ax$$ $${\lambda }_{max}\left(A\right)=\underset{x\ne 0}{max}\frac{x^{\prime }Ax}{x^{\prime }x}=\underset{\parallel x\parallel =1}{max}{x}^{\prime }Ax$$

• This gives PCA a nice interpretation: the first principal component maximizes the explained variance

Computing Eigenvalues

There are different cases of computing eigenvalues

• Computing all eigenvalues (eigenvalue decomposition)

• Computing the largest eigenvalue

• Computing the eigenvalue that is closest to a given number

• ...

Computing Eigenvalues

• Eigenvalue computation is a large field in numerical linear algebra

• The actual algorithms used in practice are usually quite involved

• For the purpose of statistical computing, we mainly introduce the basic ideas behind these algorithms and some recommended implementations

• In what follows, we focus on finding eigenvalues of real symmetric matrices

Computing All Eigenvalues

The most commonly-used method to compute all eigenvalues consists of two steps:

1. Looking for an orthogonal matrix $Q$ such that $Q^{\prime }AQ$ is tridiagonal

$$Q^{\prime }AQ=T=\left[\begin{array}{ccccc}{\alpha }_1& {\beta }_2& & & \\ {\beta }_2& {\alpha }_2& {\beta }_3& & \\ & \ddots & \ddots & \ddots & \\ & & {\beta }_{n-1}& {\alpha }_{n-1}& {\beta }_n\\ & & & {\beta }_n& {\alpha }_n\end{array}\right]$$

1. Applying the QR algorithm (introduced later) to the matrix $T$ to compute all eigenvalues

The Householder (Tridiagonalization) Algorithm

• The Householder algorithm finds a sequence of (simple) orthogonal matrices $Q_k$ such that $$T=(Q_1\cdots Q_{n-2}{)}^{\prime }A(Q_1\cdots Q_{n-2}),$$ where $T$ is tridiagonal

• $Q=Q_1\cdots Q_{n-2}$ is also orthogonal (why?)

• $A$ and $T$ have the same eigenvalues (why?)

• The motivation is that $T$ has cheaper matrix operations

The Householder (Tridiagonalization) Algorithm

• A very nice intro orthogonal transformatioduction to the Householder algorithm can be found at Section 4.3 and Section 4.6.1 of https://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

The QR Algorithm

The QR algorithm iteratively reduces the tridiagonal $T$ into a diagonal matrix by orthogonal transformations:

1. Initialize $T^{\left(0\right)}=T$, and begin the loop $k=0,1,\dots$

2. Compute the QR decomposition of the matrix $T^{\left(k\right)}$: $$T^{\left(k\right)}=Q_k{R}_k$$

3. Update $T^{(k+1)}=R_k{Q}_k=Q_k^{\prime }{T}^{\left(k\right)}{Q}_k$

Assume that $|{\lambda }_1|>\cdots >|{\lambda }_n|>0$, and then $T^{\left(k\right)}$ converges to a diagonal matrix, whose diagonal elements are eigenvalues of $T$ and $A$.

Practical Algorithm

• The lower triangular elements $T_{ij}^{\left(k\right)}$, $i>j$ converge to zero at the rate of $O\left(|{\lambda }_i/{\lambda }_j|\right)$ ^[2]
• To accelerate the convergence, a shift can be applied to the original matrix:

1. Compute the QR decomposition of the shifted matrix, $$T^{\left(k\right)}-{\sigma }_k{I}_n=Q_k{R}_k$$
2. Update $T^{(k+1)}=R_k{Q}_k+{\sigma }_k{I}_n=Q_k^{\prime }{T}^{\left(k\right)}{Q}_k$

• We can take ${\sigma }_k=(T^{\left(k\right)}{)}_{nn}$ until $T_{n-1,n}^{\left(k\right)}\approx 0$, and then switch to ${\sigma }_k=(T^{\left(k\right)}{)}_{n-1,n-1}$

Complexity

• The Householder algorithm costs $O\left(n^3\right)$

• The QR decomposition of tridiagonal matrices has a specialized algorithm, which costs $O\left(n^2\right)$ instead of $O\left(n^3\right)$ (the latter for general dense matrices)

• It can be proved that all $T^{\left(k\right)}$ matrices are tridiagonal

• The number of iterations required in the QR algorithm depends on the distribution of eigenvalues

Computing the Largest Eigenvalue

• In many applications, we only need the largest eigenvalue (in absolute value)

• It can be extended to computing the largest $k$ eigenvalues

• If ${\lambda }_1,\dots ,{\lambda }_n$ are eigenvalues of $A$ with $|{\lambda }_1|>\cdots >|{\lambda }_n|$, and ${\gamma }_1,\dots ,{\gamma }_n$ are the associated eigenvectors, then:

  • ${\lambda }_2$ is the largest eigenvalue of $A-{\lambda }_1{\gamma }_1{\gamma }_1^{\prime }$
  • ${\lambda }_3$ is the largest eigenvalue of $A-{\lambda }_1{\gamma }_1{\gamma }_1^{\prime }-{\lambda }_2{\gamma }_2{\gamma }_2^{\prime }$
  • ...

Power Method ^[2]

One widely-used and possibly the simplest method to compute the largest eigenvalue (in absolute value) is the power method

Convergence ^[2]

Assume $A$ has eigenvalues ${\lambda }_1>|{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$ and associated eigenvectors ${\gamma }_1,\dots ,{\gamma }_n$. Also assume the initial vector satisfies $x_0={\sum }_{i=1}^n{\beta }_i{\gamma }_i$ for some ${\beta }_1\ne 0$. Then

$$|\parallel y_k\parallel -{\lambda }_1|\le C{\left(\frac{|{\lambda }_2|}{{\lambda }_1}\right)}^k,\parallel x_k-{\tilde{\gamma }}_1\parallel \le C{\left(\frac{|{\lambda }_2|}{{\lambda }_1}\right)}^k,$$ where $\widetilde{{\gamma }_1}=\pm {\gamma }_1$.

Pros and Cons

• The power method only requires the matrix-vector multiplication $v\to Av$ to compute eigenvalues of $A$ (recall the CG method)

• This is important and useful for sparse matrices

• However, it does not effectively use the information of past iterations

• Also, it can only compute one eigenvalue at one time

Advanced Algorithms

• Mainstream software packages typically use other advanced algorithms to compute the largest eigenvalues, e.g.,

  • Implicitly restarted Lanczos method ^[4]
  • Jacobi–Davidson algorithm ^[5]

• They also only need the $v\to Av$ operation

• Can compute multiple eigenvalues together

• These algorithms themselves are quite sophisticated, however

Computing Eigenvalue Closest to $\mu$

• Another type of problem is to find the eigenvalue that is closest to a given number $\mu$

• For example, if $\mu =0$, then this is equivalent to finding the smallest eigenvalue (in absolute value)

• Fortunately, we can make use of existing methods that compute largest eigenvalues for this task, via a technique called spectral transformation

Spectral Transformation

Spectral transformation is based on the following finding.

If ${\lambda }_i$, $i=1,\dots ,n$ are eigenvalues of $A$, then $A{x}_i={\lambda }_i{x}_i$ for the eigenvectors $x_i$. Choose a shift $\mu \ne {\lambda }_i$, then the eigenvalues of $B=(A-\mu I_n{)}^{-1}$ are ${\theta }_i=({\lambda }_i-\mu {)}^{-1}$, and $B{\theta }_i={\theta }_i{x}_i$.

Therefore, the largest eigenvalue of $B$ corresponds to the eigenvalue of $A$ that is closest to $\mu$.

Inverse Iteration

• If we want to compute the eigenvalue of $A$ closest to $\mu$

• We can apply the power method (or other methods based on matrix-vector multiplication) to $B=(A-\mu I_n{)}^{-1}$

• The iteration becomes $$\begin{array}{cc}{y}_k& =(A-\mu I_n{)}^{-1}{x}_{k-1}\\ x_k& =y_k/\parallel y_k\parallel \end{array}$$

• By the convergence theorem, we have $\parallel y_k\parallel \to {\theta }_1$

• Then recover the eigenvalue of $A$ by ${\lambda }_{\mu }=\mu +1/{\theta }_1$

Solving Linear Systems

• In actual implementation, we do not directly compute the matrix inverse $(A-\mu I_n{)}^{-1}$

• Instead, we factorize $A-\mu I_n$ once using LU decomposition or $LD{L}^T$ decomposition

• And then solve the linear systems $(A-\mu I_n)y_k=x_{k-1}$

Implementation

• The eigen() function in base R is a very stable and efficient function to compute all eigenvalues

• Specify symmetric = TRUE if you know your matrix is symmetric (by default it will detect symmetry)

• Specify only.values = TRUE if you do not need eigenvectors

Benchmark

library(dplyr)
set.seed(123)
n = 1000
M = matrix(rnorm(n^2), n)
A = M + t(M)
bench::mark(
    eigen(A),
    eigen(A, symmetric = FALSE),
    eigen(A, symmetric = TRUE),
    eigen(A, symmetric = TRUE, only.values = TRUE),
    min_iterations = 3, max_iterations = 10, check = FALSE
) %>% select(expression, min, median)

## # A tibble: 4 × 3
##   expression                                          min   median
##   <bch:expr>                                     <bch:tm> <bch:tm>
## 1 eigen(A)                                          1.15s    1.17s
## 2 eigen(A, symmetric = FALSE)                       4.13s    4.14s
## 3 eigen(A, symmetric = TRUE)                        1.11s    1.13s
## 4 eigen(A, symmetric = TRUE, only.values = TRUE) 314.12ms 314.19ms

Implementation

• To compute a subset of eigenvalues up to some selection rule, we have several options

• The RSpectra package implements the implicitly restarted Lanczos method
• The PRIMME package uses the Jacobi-Davidson algorithm
• The irlba package is designed for singular value decomposition (next class), but can also be used to compute the largest eigenvalues

• My personal favorite is the RSpectra package (this is a highly biased choice, of course)

Implementation

• The code below computes three eigenvalues with largest absolute values

library(RSpectra)
e = eigs_sym(A, k = 3, which = "LM")
e$values

## [1]  88.83776  88.31500 -88.98945

head(e$vectors)

##             [,1]         [,2]         [,3]
## [1,]  0.01506291  0.029602478 -0.015855700
## [2,] -0.05124211  0.016404205 -0.014749483
## [3,] -0.01542738 -0.002350839  0.005520141
## [4,]  0.02921002 -0.006653817 -0.003387408
## [5,] -0.02406133  0.039802584  0.028103112
## [6,] -0.01309102 -0.026105929  0.073813795

Implementation

• If you want the largest three eigenvalues without taking absolute values, use which = "LA"

e = eigs_sym(A, k = 3, which = "LA")
e$values

## [1] 88.83776 88.31500 87.18646

• If no eigenvectors are needed, add the following option

e = eigs_sym(A, k = 3, which = "LA", opts = list(retvec = FALSE))
e$vectors

## NULL

Benchmark

• If only part of the eigenvalues are needed, eigs_sym() is much faster than eigen()

bench::mark(
    eigen(A, symmetric = TRUE),
    eigen(A, symmetric = TRUE, only.values = TRUE),
    eigs_sym(A, k = 3, which = "LM"),
    eigs_sym(A, k = 3, which = "LM", opts = list(retvec = FALSE)),
    min_iterations = 3, max_iterations = 10, check = FALSE
) %>% select(min, median)

## # A tibble: 4 × 2
##        min   median
##   <bch:tm> <bch:tm>
## 1    1.12s    1.12s
## 2 310.66ms 310.66ms
## 3 149.93ms 150.31ms
## 4 149.43ms 150.85ms

• The difference will be more evident on sparse matrices (next class)

Implementation

• If one needs the smallest three eigenvalues (in absolute value)

e = eigs_sym(A, k = 3, which = "LM", sigma = 0)
e$values

## [1]  0.01589015 -0.14248784 -0.26307129

• In general, if sigma is not NULL, then the selection rule is applied to $1/({\lambda }_i-\sigma )$

• This only affects the selection rule

• The returned eigenvalues are still for $A$

References