layout: true --- class: inverse, center, middle background-image: url(../figs/titlepage16-9.png) background-size: cover <br> <br> # Bayesian Statistics and Computing ## Lecture 2: Linear Algebra in Statistical Computing <img src="../figs/slides.png" width="150px"/> #### *Yanfei Kang | BSC | Beihang University* --- class: inverse, center, middle # Why linear algebra # in statistical computing? --- # Matrices are everywhere in statistics Once data are organized as a table, you are already doing matrix computations. - Regression: the design matrix `\(\mathbf X\)` and the response `\(\mathbf y\)` - Covariance and correlation structure of multivariate data - Transition matrices in Markov chains; adjacency matrices in networks - Term--document matrices in text analysis; pixel matrices in images The question of this lecture is **not** how to *prove* theorems about matrices, but a practical one: > Once a statistical problem is written in matrix form, how do we compute it **fast**, **stably**, and **verifiably**? --- # A correct formula can be a bad algorithm The least squares estimator is usually written as `$$\widehat{\boldsymbol\beta} =(\mathbf X^\top\mathbf X)^{-1}\mathbf X^\top\mathbf y.$$` This formula is essential for *understanding* the estimator. -- But as *code*, it suggests two things we should almost never do: - explicitly form the inverse `\((\mathbf X^\top\mathbf X)^{-1}\)`; - explicitly form the matrix `\(\mathbf X^\top\mathbf X\)`. Statistical computing is not "translating formulas into code": it is choosing a numerical method that respects the structure of the problem. --- # What this lecture covers 1. **Linear systems**: `solve()`, residuals, conditioning 2. **Symmetric positive definite systems**: Cholesky 3. **Least squares**: the normal equations trap, QR 4. **Eigenanalysis**: matrix powers, PageRank 5. **SVD**: works for any matrix; rank, conditioning, low rank 6. **Applications**: image compression, pseudo-inverse 7. **Iterative methods**: power method, QR algorithm As you listen, keep three questions in mind: 1. What is the **goal** of the computation? 2. What **structure** does the matrix have? 3. Which **quantities** verify the result? --- class: inverse, center, middle # Linear systems: # solve, don't invert --- # Where do linear systems come from? Many statistical computations eventually require solving `$$\mathbf A\mathbf x=\mathbf b.$$` Sources in statistics: - conditional distributions of multivariate normals - Newton directions in optimization (positive definite Hessians) - ridge-type penalized estimates - covariance structures `\(\Rightarrow\)` symmetric positive definite systems --- # Solve the equation; do not form the inverse Mathematically `\(\mathbf x=\mathbf A^{-1}\mathbf b\)`. Computationally, treat the task as "solve a system". - Forming `\(\mathbf A^{-1}\)` costs more flops and storage and adds rounding error. - In R write `solve(A, b)`, **not** `solve(A) %*% b`. ```r set.seed(2026) A <- matrix(rnorm(25), 5, 5) + 5 * diag(5) b <- rnorm(5) x_direct <- solve(A, b) x_via_inverse <- solve(A) %*% b max(abs(x_direct - x_via_inverse)) #> [1] 1.1102e-16 ``` Same answer here --- but the first route is cheaper and more stable. --- class: inverse, center, middle # Symmetric positive definite systems: # Cholesky --- # SPD matrices are native to statistics Covariance matrices, positive definite Hessians, and ridge-type matrices are symmetric positive definite (SPD). If `\(\mathbf A\in\mathbb R^{p\times p}\)` is SPD, there is a unique upper-triangular `\(\mathbf R\)` with positive diagonal such that `$$\mathbf A=\mathbf R^\top\mathbf R.$$` This is the **Cholesky factorization**. Solving `\(\mathbf A\mathbf x=\mathbf b\)` becomes two triangular solves: `$$\mathbf R^\top\mathbf z=\mathbf b \qquad\text{(forward solve)}, \qquad \mathbf R\mathbf x=\mathbf z \qquad\text{(back solve)}.$$` In R: `chol()`, `forwardsolve()`, `backsolve()`. --- # When `chol()` fails ```r A_ind <- matrix(c(1, 2, 2, 1), 2, 2) # eigenvalues 3 and -1 eigen(A_ind, only.values = TRUE)$values #> [1] 3 -1 chol(A_ind) #> Error in chol.default(A_ind): the leading minor of order 2 is not positive ``` `chol()` fails when the matrix is **not** symmetric positive definite: exact collinearity, or rounding pushing a tiny eigenvalue below zero. --- class: inverse, center, middle # Least squares by QR --- # Data give overdetermined systems The systems so far (general square, SPD) have a unique exact solution. Regression is the same equation in its statistical spelling: `\(\mathbf A\to\)` design matrix `\(\mathbf X\)`, `\(\mathbf x\to\)` coefficients `\(\boldsymbol\beta\)`, `\(\mathbf b\to\)` response `\(\mathbf y\)`. What changes is the shape: with `\(n\)` observations and `\(p\)` parameters, `\(\mathbf X\)` is `\(n\times p\)`, and `$$\mathbf y=\mathbf X\boldsymbol\beta$$` generally has **no exact solution** --- `\(\mathbf y\)` rarely lies in the column space of `\(\mathbf X\)`. "Solving" is then redefined as projection: minimize the squared residual `$$\min_{\boldsymbol\beta}\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2.$$` --- # Least squares via the normal equations Route 1: square the system up, then reuse the SPD machinery of the previous section. For a full-column-rank design matrix, the minimizer satisfies the normal equations `$$\mathbf X^\top\mathbf X\widehat{\boldsymbol\beta}=\mathbf X^\top\mathbf y,$$` and `\(\mathbf X^\top\mathbf X\)` is symmetric positive definite (in exact arithmetic): 1. Cholesky factorization `\(\mathbf X^\top\mathbf X=\mathbf R^\top\mathbf R\)`. 2. Forward solve `\(\mathbf R^\top\mathbf w=\mathbf X^\top\mathbf y\)`. 3. Back solve `\(\mathbf R\widehat{\boldsymbol\beta}=\mathbf w\)`. No explicit inverse --- but forming `\(\mathbf X^\top\mathbf X\)` squares the condition number, `\(\kappa_2(\mathbf X^\top\mathbf X)=\kappa_2(\mathbf X)^2\)`, and rounding can make `chol(crossprod(X))` fail even when `\(\mathbf X\)` has full column rank mathematically. --- # The QR factorization Any `\(m\times n\)` matrix `\(\mathbf A\)` with `\(m\geq n\)` can be written as `$$\mathbf A=\mathbf Q\mathbf R,$$` with `\(\mathbf Q\)` an `\(m\times m\)` orthogonal matrix and `\(\mathbf R\)` an `\(m\times n\)` upper trapezoidal matrix. For tall matrices the last `\(m-n\)` rows of `\(\mathbf R\)` are zero, so we only need the **thin** (economy) form: `$$\mathbf A=\mathbf Q\mathbf R=(\mathbf Q_1,\mathbf Q_2) \left(\begin{array}{cc}\mathbf R_1 \\ \mathbf 0 \end{array}\right)=\mathbf Q_1\mathbf R_1,$$` where `\(\mathbf Q_1\)` has orthonormal columns and `\(\mathbf R_1\)` is `\(n\times n\)` upper triangular. --- # How is QR computed? Gram--Schmidt Any basis `\((\mathbf a_1,\ldots,\mathbf a_n)\)` can be turned into an orthonormal basis `\((\mathbf q_1,\ldots,\mathbf q_n)\)` by projecting out, from each vector, the components along the previous directions. <center> <img src="./figs/gso.png" height="300px"/> </center> --- # From Gram--Schmidt to QR Collecting the projection coefficients of `\(\mathbf A=\mathbf Q\mathbf R\)` gives exactly the Gram--Schmidt quantities: `\(\mathbf R\)` holds the coefficients, `\(\mathbf Q\)` the orthonormalized directions. <center> <img src="./figs/QR.png" height="230px"/> </center> --- # Least squares via QR Route 2: do not square up --- factor the design matrix itself (the QR theorem with `\(\mathbf A=\mathbf X\)`): `$$\mathbf X=\mathbf Q_1\mathbf R_1, \qquad \mathbf Q_1\in\mathbb R^{n\times p},\ \mathbf R_1\in\mathbb R^{p\times p}\ \text{upper triangular}.$$` 1. Thin QR factorization `\(\mathbf X=\mathbf Q_1\mathbf R_1\)`. 2. Compute `\(\mathbf c=\mathbf Q_1^\top\mathbf y\)`. 3. Solve `\(\mathbf R_1\widehat{\boldsymbol\beta}=\mathbf c\)` by back substitution. --- # Advantages of QR - No need to form `\(\mathbf X^\top\mathbf X\)`: avoid squaring the condition number. - QR is generally more numerically reliable than Cholesky applied to the normal equations. - `lm()`, `lm.fit()` and the default `qr()` route use QR with pivoting to handle numerical rank. QR improves the computation; it does not remove the sensitivity caused by nearly dependent columns. --- # Longley data: three routes compared `longley`: 16 annual observations, highly collinear macro variables. .scroll-output[ ```r X <- model.matrix( Employed ~ GNP + Unemployed + Armed.Forces + Population + Year, data = longley ) y <- longley$Employed beta_qr <- qr.coef(qr(X), y) beta_lm <- lm.fit(X, y)$coefficients beta_ne <- solve(crossprod(X), crossprod(X, y)) round(cbind(QR = beta_qr, lm = beta_lm, normaleq = beta_ne), 6) #> QR lm #> (Intercept) -3.4499e+03 -3.4499e+03 -3.4499e+03 #> GNP -3.1961e-02 -3.1961e-02 -3.1961e-02 #> Unemployed -1.9721e-02 -1.9721e-02 -1.9721e-02 #> Armed.Forces -1.0200e-02 -1.0200e-02 -1.0200e-02 #> Population -7.7537e-02 -7.7537e-02 -7.7537e-02 #> Year 1.8141e+00 1.8141e+00 1.8141e+00 ``` ] --- # The condition numbers tell the story ```r X_scaled <- cbind(`(Intercept)` = 1, scale(X[, -1])) data.frame( matrix = c("X", "crossprod(X)", "scaled X"), kappa_2 = c( kappa(X, exact = TRUE), kappa(crossprod(X), exact = TRUE), kappa(X_scaled, exact = TRUE) ) ) #> matrix kappa_2 #> 1 X 2.3311e+07 #> 2 crossprod(X) 5.4341e+14 #> 3 scaled X 8.0668e+01 ``` - `\(\kappa_2(\mathbf X^\top\mathbf X)\)` is orders of magnitude larger than `\(\kappa_2(\mathbf X)\)`: the squared-conditioning trap in action. - Standardization fixes *scaling* problems, **not** genuine near-linear dependence among variables. --- # Example: a nearly collinear design Generate `\(n=5000\)` observations with `\(p=100\)` predictors, then add an almost duplicate of the first column. ```r set.seed(2023-09-08) n_qr <- 5000; p_qr <- 100 X_qr <- matrix(rnorm(n_qr * p_qr), n_qr, p_qr) y_qr <- rnorm(n_qr) W_qr <- cbind(X_qr, X_qr[, 1] + rnorm(n_qr, sd = 1e-14)) ``` --- # Collinearity: normal equations ```r solve(crossprod(W_qr), crossprod(W_qr, y_qr)) #> Error in solve.default(crossprod(W_qr), crossprod(W_qr, y_qr)): system is computationally singular: reciprocal condition number = 1.89233e-16 ``` The almost duplicate column makes `\(\mathbf W^\top\mathbf W\)` computationally singular: the normal-equation solve fails. --- # Collinearity: QR decomposition .scroll-output[ ```r qr_W <- qr(W_qr) beta_qr_W <- qr.coef(qr_W, y_qr) round(beta_qr_W, 3) #> [1] 0.021 -0.011 -0.005 -0.014 -0.017 -0.025 -0.006 -0.005 -0.014 0.001 0.011 0.024 -0.001 #> [14] 0.040 -0.005 0.016 0.037 -0.006 0.000 0.000 0.011 0.003 -0.010 0.009 -0.006 -0.010 #> [27] 0.000 0.008 0.010 -0.021 -0.005 0.006 -0.007 -0.004 0.007 -0.001 -0.026 -0.005 -0.020 #> [40] -0.011 0.010 0.002 0.010 -0.026 0.000 -0.020 0.020 0.025 0.014 -0.008 -0.033 0.007 #> [53] -0.025 -0.015 0.012 0.010 -0.004 0.012 0.006 0.002 0.001 -0.007 -0.010 -0.009 0.005 #> [66] -0.022 0.018 0.014 0.005 0.017 -0.025 -0.001 -0.011 -0.002 0.008 0.011 -0.014 -0.008 #> [79] 0.011 -0.001 -0.002 -0.009 0.023 -0.009 0.010 -0.019 -0.015 0.013 -0.009 0.005 -0.021 #> [92] -0.013 0.020 0.003 -0.009 0.006 -0.002 -0.015 -0.002 -0.016 NA c(columns = ncol(W_qr), numerical_rank = qr_W$rank, aliased_coefficients = sum(is.na(beta_qr_W))) #> columns numerical_rank aliased_coefficients #> 101 100 1 ``` ] QR detects a numerical rank of 100 among 101 columns and marks the redundant coefficient as `NA`. It does not claim that all 101 coefficients are identifiable. --- class: inverse, center, middle # Eigenanalysis --- # Eigenvectors and eigenvalues Everything so far asked: given `\(\mathbf A\)` and `\(\mathbf b\)`, find a usable solution `\(\mathbf x\)` --- exact or least squares. Eigenanalysis asks a different question about `\(\mathbf A\)` itself. - In general, a matrix acts on a vector by changing both its magnitude and its direction. - However, a matrix may act on certain vectors by changing only their magnitude, and leaving their direction unchanged (or reversing it). - these vectors are the eigenvectors of the matrix - the scaling factor is the eigenvalue. <center> <img src="./figs/eigen.png" height="200px" /> </center> --- # The symmetric case - If `\(\mathbf{A}\)` is a real symmetric matrix - all eigenvalues are real - `\(n\)` orthonormal eigenvectors - `\(\mathbf{A}=\mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^\top\)` with `\(\mathbf Q\)` orthogonal - Not every square matrix diagonalises. A triangularisation always exists, though: the Schur decomposition `\(\mathbf{A}=\mathbf{Q}\mathbf S\mathbf{Q}^\top\)`. - What if `\(\mathbf{A}\)` is not even square? Hold that thought for the SVD. --- # What does eigenanalysis help to do? Work in a coordinate system built from the eigenvectors, and `\(\mathbf{Ax} = \lambda \mathbf{x}\)` lets us replace the matrix by a scalar. - Consider `\(\mathbf{A} = \left(\begin{array}{cc}5 & -1\\-2 & 4 \end{array}\right)\)` and `\(\mathbf{x} = \left(\begin{array}{c}1 \\-1 \end{array}\right)\)`. Check that `\(\mathbf{x}\)` is an eigenvector --- what is the eigenvalue? - What if we now want to calculate `\(\mathbf{A^{20}x}\)`? - What if we now want to calculate `\(\mathbf{A^{-1}x}\)`? - Any other vector is a combination of eigenvectors, so the scalar rule applies term by term. --- # Application of Eigenanalysis: Google <center> <img src="./figs/google.png" height="500px" /> </center> --- # How does Google rank pages? - How does the search engine know which pages are the most important? - Google assigns each page a number expressing its importance: the **PageRank**, computed via the eigenvalue problem `\(\mathbf A\mathbf w=\lambda\mathbf w\)`, where `\(\mathbf A\)` is based on the link structure of the Internet. --- # Google's pagerank - Suppose we have a set of webpages `\(W\)`, with `\(|W|=n\)` as the number of webpages. - We define a connectivity (adjacency) matrix `\(\mathbf A\)` as follows: `$$\mathbf A_{i j}:=\left\{\begin{array}{ll} 1 / N_{i} & \text{if Page } i \text{ links to Page } j, \\ 0 & \text { otherwise.} \end{array}\right.$$` where `\(N_i\)` is the number of outlinks on Page `\(i\)`. - Typically `\(\mathbf A\)` is huge (27.5 billion `\(\times\)` 27.5 billion in 2011), but extremely sparse (lots of zero values) --- # A small example .pull-left[ Consider the small set of five webpages: <br>  ] -- .pull-right[ The connectivity matrix is `$$\mathbf A =\left(\begin{array}{lllll} 0 & 1/3 & 0 & 1/3 & 1/3 \\ 0 & 0 & 1/3 & 1/3 & 1/3 \\ 1/2 & 0 & 0 & 1/2 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \end{array}\right)$$` ] --- # A small example - The underlying idea for PageRank: a page is important, if other important pages link to it. - The ranking of page `\(i\)` is proportional to the weighted sum of the rankings of all pages that link to `\(i\)`: `\(r_4 = \frac{1}{3}r_1 + \frac{1}{3}r_2 + \frac{1}{2}r_3\)`. - In vector form, `\(\mathbf r=\mathbf A^\top\mathbf r\)`: a system of `\(n = 5\)` linear equations, and an eigenvector equation for `\(\mathbf A^\top\)` with eigenvalue 1 (please write them out). - The ranking vector is treated like a probability of relevance, so we rescale it to satisfy `\(\sum_{i=1}^n r_i = 1\)`. - Go to R and compute page ranks in this example. --- # PageRank at real scale - Finding eigenvalues via the characteristic polynomial is hopeless at this scale (don't even try). - We need methods that use only matrix--vector products with a sparse `\(\mathbf A\)`, and that only look for the dominant eigenvector. - The **power method**, at the end of this lecture, is exactly such a method. --- class: inverse, center, middle # Singular Value Decomposition --- # Non-square matrices Eigendecomposition has an obvious limitation: it targets **square** matrices, and not even all of them diagonalize nicely. Statistical data matrices are rarely square: `\(n\)` observations `\(\times\)` `\(p\)` variables, image pixels, term--document matrices, user--item ratings. The SVD is the most useful "diagonal" decomposition that works for **any** real matrix. It underlies rank, conditioning, low-rank approximation, least squares with rank deficiency, and recommender systems. --- # Singular values and singular vectors For `\(\mathbf A\in\mathbb R^{m\times n}\)`, the singular values are `$$\sigma_j=\sqrt{\lambda_j(\mathbf A^\top\mathbf A)},\qquad \sigma_1\geq\sigma_2\geq\cdots\geq\sigma_p\geq 0,\quad p=\min\{m,n\}.$$` For each `\(\sigma_j>0\)` there is a pair `\((\mathbf u_j,\mathbf v_j)\)` with `$$\mathbf A\mathbf v_j=\sigma_j\mathbf u_j, \qquad \mathbf A^\top\mathbf u_j=\sigma_j\mathbf v_j.$$` `\(\mathbf v_j\)` lives in the input space, `\(\mathbf u_j\)` in the output space: the matrix maps direction `\(\mathbf v_j\)` onto direction `\(\mathbf u_j\)`, scaling lengths by `\(\sigma_j\)`. Theoretically `\(\sigma_j\)` comes from eigenvalues of `\(\mathbf A^\top\mathbf A\)`; numerically, **never form** that product because it squares the condition number. --- # The SVD For any `\(\mathbf A\in\mathbb R^{m\times n}\)`, there exist orthogonal matrices `\(\mathbf U = [\mathbf u_1, \cdots, \mathbf u_m]\in\mathbb R^{m\times m}\)` and `\(\mathbf V = [\mathbf v_1, \cdots, \mathbf v_n]\in\mathbb R^{n\times n}\)` such that `$$\mathbf U^\top\mathbf A\mathbf V=\boldsymbol\Sigma,\qquad \mathbf A=\mathbf U\boldsymbol\Sigma\mathbf V^\top,$$` with `\(\boldsymbol\Sigma\in\mathbb R^{m\times n}\)` diagonal, entries `\(\sigma_1\geq\cdots\geq\sigma_p\geq 0\)`. <center> <img src="https://sthalles.github.io/assets/svd-for-regression/full-svd-matrices.png" height="180px"/> </center> --- # SVD in R ```r A_svd <- matrix(c(3, 2, 2, 3, 1, 1), 3, 2, byrow = TRUE) s <- svd(A_svd) s$d #> [1] 5.1962 1.0000 round(s$u %*% diag(s$d) %*% t(s$v), 6) # reconstruction #> [,1] [,2] #> [1,] 3 2 #> [2,] 2 3 #> [3,] 1 1 ``` The reconstruction equals `\(\mathbf A\)` up to rounding. `svd()` returns `d` (singular values), `u`, `v`. --- # Some properties of SVD - `\(\sigma_j\)` are the singular values of `\(\mathbf A\)`. - The non-zero singular values of `\(\mathbf A\)` are the square roots of the non-zero eigenvalues of both `\(\mathbf A^\top\mathbf A\)` and `\(\mathbf A\mathbf A^\top\)`. - The rank of `\(\mathbf A\)` equals the number of non-zero singular values. - The condition number measures how close `\(\mathbf A\)` is to rank deficiency. For a full-column-rank `\(\mathbf A\)`: `$$\kappa_2(\mathbf A)=\frac{\sigma_1}{\sigma_p}.$$` --- # Economy SVD and the outer-product form `$$\mathbf A=\mathbf U\boldsymbol\Sigma\mathbf V^\top =\left(\begin{array}{cc}\mathbf U_1 & \mathbf U_2\end{array}\right) \left(\begin{array}{cc}\mathbf D & \mathbf 0\\ \mathbf 0 & \mathbf 0\end{array}\right) \left(\begin{array}{c}\mathbf V_1^\top \\ \mathbf V_2^\top\end{array}\right) =\mathbf U_1\mathbf D\mathbf V_1^\top,$$` where `\(\mathbf D=\operatorname{diag}\{\sigma_1,\cdots,\sigma_r\}\)` with `\(\sigma_r>0\)` and `\(r=\operatorname{rank}(\mathbf A)\)`. <center> <img src="https://sthalles.github.io/assets/svd-for-regression/economy-svd-matrices.png" height="170px"/> </center> The matrix is a **sum of rank-one pieces**, ordered by contribution. Truncating the sum is the basis of compression and noise filtering. --- # Best low-rank approximation Keep the first `\(k\)` terms: `$$\mathbf A_k=\sum_{j=1}^k\sigma_j\mathbf u_j\mathbf v_j^\top.$$` In the spectral norm, the **optimal rank `\(k\)`** approximation error is `$$\lVert\mathbf A-\mathbf A_k\rVert_2=\sigma_{k+1}.$$` If singular values decay fast, a few directions carry most of the information. In floating point, singular values are rarely exactly zero, so we define the **numerical rank** with a relative threshold `\(\tau\)`: `$$r_\tau=\#\{j:\ \sigma_j>\tau\sigma_1\}.$$` --- class: inverse, center, middle # Applications of SVD --- # Images are matrices - A grayscale image is one matrix of pixel intensities; a colour image has three matrices of the same size, one for each RGB channel. - Each channel uses intensities from 0 to 255 in an 8-bit image; `jpeg::readJPEG()` returns values scaled to [0, 1]. - Compressing the image is then exactly the low-rank approximation from the previous slides: keep the leading `\(k\)` singular values of each channel. --- # Reconstructing the approximated image - Choose `\(k\)`, rebuild each channel from its first `\(k\)` singular values and singular vectors, and clip the intensities back to [0, 1]. - For colour images, reconstruct the R, G, B matrices separately and then reassemble. Storing an `\(m \times n\)` channel takes `\(mn\)` numbers; its rank-$k$ factors take `\(k(m+n+1)\)` numbers. --- # Data compression .scroll-output[ ```r library(jpeg) img <- readJPEG('./figs/svdimg.jpg') r <- img[,,1] g <- img[,,2] b <- img[,,3] img.r.svd <- svd(r) img.g.svd <- svd(g) img.b.svd <- svd(b) # Function to reconstruct the image with a limited number of singular values reconstruct_image <- function(svd_r, svd_g, svd_b, k) { # Reconstruct the individual color channels using the first k singular values r_reconstructed <- svd_r$u[, 1:k] %*% diag(svd_r$d[1:k]) %*% t(svd_r$v[, 1:k]) g_reconstructed <- svd_g$u[, 1:k] %*% diag(svd_g$d[1:k]) %*% t(svd_g$v[, 1:k]) b_reconstructed <- svd_b$u[, 1:k] %*% diag(svd_b$d[1:k]) %*% t(svd_b$v[, 1:k]) # Clip the values to be within the valid range [0, 1] r_reconstructed <- pmin(pmax(r_reconstructed, 0), 1) g_reconstructed <- pmin(pmax(g_reconstructed, 0), 1) b_reconstructed <- pmin(pmax(b_reconstructed, 0), 1) # Combine the channels back into a single image reconstructed_image <- array(0, dim = c(nrow(svd_r$u), nrow(svd_r$v), 3)) reconstructed_image[,,1] <- r_reconstructed reconstructed_image[,,2] <- g_reconstructed reconstructed_image[,,3] <- b_reconstructed return(reconstructed_image) } # Set the number of singular values to test singular_values <- seq(5,200,10) # Loop over different k values and plot the reconstructed images for (k in singular_values) { reconstructed <- reconstruct_image(img.r.svd, img.g.svd, img.b.svd, k) plot(1, type = "n", xlab = "", ylab = "", xlim = c(0, 1), ylim = c(0, 1), axes = FALSE, ann = FALSE, frame.plot = FALSE) rasterImage(reconstructed, 0, 0, 1, 1) title(paste("Reconstructed Image with", k, "Singular Values")) } ``` ] --- # Data compression <center> <img src="./figs/image-compression.gif" height="420px" alt="SVD image reconstruction with increasing numbers of singular values"/> </center> --- # Pseudo-inverse Write the SVD of `\(\mathbf A\in\mathbb R^{m\times n}\)` in partitioned form, `$$\mathbf A=\mathbf U\boldsymbol\Sigma\mathbf V^\top =\left(\begin{array}{cc}\mathbf U_1 & \mathbf U_2\end{array}\right) \left(\begin{array}{cc}\mathbf D & \mathbf 0\\ \mathbf 0 & \mathbf 0\end{array}\right) \left(\begin{array}{c}\mathbf V_1^\top \\ \mathbf V_2^\top\end{array}\right) =\mathbf U_1\mathbf D\mathbf V_1^\top,$$` where `\(\mathbf D=\operatorname{diag}\{\sigma_1,\cdots,\sigma_r\}\)` with `\(\sigma_r>0\)` and `\(r=\operatorname{rank}(\mathbf A)\)`. - The pseudo-inverse, or generalised inverse, of `\(\mathbf A\)` can always be calculated, even if `\(\mathbf A\)` is not invertible: `$$\mathbf A^{+}=\mathbf V_1\mathbf D^{-1}\mathbf U_1^\top$$` - `\(\mathbf A^{+}\)` is `\(n\times m\)`. --- # Linear regression with SVD For `\(\mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon\)`, minimize `\(\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2\)`. With the economy SVD `\(\mathbf X=\mathbf U_1\mathbf D\mathbf V_1^\top\)` and `\(r=\operatorname{rank}(\mathbf X)\)`, `$$\widehat{\boldsymbol\beta}=\mathbf X^+\mathbf y =\mathbf V_1\mathbf D^{-1}\mathbf U_1^\top\mathbf y =\sum_{i=1}^{r}\frac{\mathbf u_i^\top\mathbf y}{\sigma_i}\mathbf v_i.$$` **Least squares via SVD (Route 3):** 1. Compute the economy SVD, retaining the positive singular values. 2. Project and rescale: `\(\mathbf c=\mathbf U_1^\top\mathbf y\)`, then `\(w_i=c_i/\sigma_i\)`. 3. Transform back: `\(\widehat{\boldsymbol\beta}=\mathbf V_1\mathbf w\)`. --- # Collinear design, revisited by SVD .scroll-output[ ```r # Perform SVD on W svd_W <- svd(W_qr) # Extract components (svd() returns the economy U and V by default) U <- svd_W$u # 5000 x 101 D <- svd_W$d # singular values, length 101 V <- svd_W$v # 101 x 101 # Define a threshold to filter out small singular values threshold <- 1e-10 Positive <- D > threshold # Calculate ginv W_ginv <- V[, Positive, drop = FALSE] %*% ((1/D[Positive]) * t(U[, Positive, drop = FALSE])) # You can also use # W_ginv <- MASS::ginv(W_qr) # Estimation W_ginv %*% y_qr #> [,1] #> [1,] 0.01035248 #> [2,] -0.01077988 #> [3,] -0.00533464 #> [4,] -0.01390428 #> [5,] -0.01656258 #> [6,] -0.02473690 #> [7,] -0.00559365 #> [8,] -0.00506709 #> [9,] -0.01418955 #> [10,] 0.00146504 #> [11,] 0.01058719 #> [12,] 0.02368550 #> [13,] -0.00091239 #> [14,] 0.03956990 #> [15,] -0.00454597 #> [16,] 0.01564260 #> [17,] 0.03689864 #> [18,] -0.00622416 #> [19,] 0.00037195 #> [20,] -0.00010454 #> [21,] 0.01054372 #> [22,] 0.00320691 #> [23,] -0.00991130 #> [24,] 0.00920917 #> [25,] -0.00562838 #> [26,] -0.01018367 #> [27,] 0.00044099 #> [28,] 0.00790167 #> [29,] 0.01016269 #> [30,] -0.02084090 #> [31,] -0.00549245 #> [32,] 0.00642719 #> [33,] -0.00655908 #> [34,] -0.00378638 #> [35,] 0.00721263 #> [36,] -0.00136867 #> [37,] -0.02615884 #> [38,] -0.00452603 #> [39,] -0.02046003 #> [40,] -0.01087911 #> [41,] 0.01005100 #> [42,] 0.00218466 #> [43,] 0.00960766 #> [44,] -0.02567444 #> [45,] -0.00020722 #> [46,] -0.02030924 #> [47,] 0.01994909 #> [48,] 0.02454570 #> [49,] 0.01413681 #> [50,] -0.00758212 #> [51,] -0.03319998 #> [52,] 0.00707289 #> [53,] -0.02516811 #> [54,] -0.01455012 #> [55,] 0.01249593 #> [56,] 0.01017242 #> [57,] -0.00416696 #> [58,] 0.01183680 #> [59,] 0.00587026 #> [60,] 0.00173228 #> [61,] 0.00085731 #> [62,] -0.00685357 #> [63,] -0.01012746 #> [64,] -0.00936435 #> [65,] 0.00451944 #> [66,] -0.02244896 #> [67,] 0.01754686 #> [68,] 0.01355087 #> [69,] 0.00487726 #> [70,] 0.01684452 #> [71,] -0.02545877 #> [72,] -0.00127009 #> [73,] -0.01063024 #> [74,] -0.00177621 #> [75,] 0.00774909 #> [76,] 0.01109049 #> [77,] -0.01447771 #> [78,] -0.00797451 #> [79,] 0.01120052 #> [80,] -0.00126976 #> [81,] -0.00240462 #> [82,] -0.00942650 #> [83,] 0.02340285 #> [84,] -0.00878272 #> [85,] 0.00972327 #> [86,] -0.01883327 #> [87,] -0.01511776 #> [88,] 0.01294354 #> [89,] -0.00860643 #> [90,] 0.00480738 #> [91,] -0.02096349 #> [92,] -0.01338732 #> [93,] 0.02034945 #> [94,] 0.00318856 #> [95,] -0.00882699 #> [96,] 0.00641799 #> [97,] -0.00169024 #> [98,] -0.01498294 #> [99,] -0.00181189 #> [100,] -0.01615946 #> [101,] 0.01035248 ``` ] --- # Discussion If SVD handles rank deficiency gracefully, why did R still make QR the default inside `lm()`? --- class: inverse, center, middle # Iterative computation of eigenvalues --- # Why iterate at all? `eigen()` and `svd()` are fine for small dense matrices. At scale, matrices can have millions of rows/columns --- full factorization is infeasible, and often **unnecessary**: - PageRank needs only the dominant eigenvector; - low-rank approximation needs only the top singular values. Iterative methods start from an initial guess and repeatedly apply cheap operations --- typically matrix--vector products --- until a convergence criterion is met. --- # The Power Method - The Power Method finds the dominant eigenvalue `\(\lambda_1\)` and corresponding dominant eigenvector `\(\mathbf q_1\)` of a matrix `\(\mathbf A\)`. - The dominant eigenvalue is the one with the largest modulus (absolute value for real eigenvalues). - The Power Method is an iterative approach that generates a sequence of scalars that converge to `\(\lambda_1\)` and a sequence of vectors that converge to `\(\mathbf q_1\)`. - The Power Method works well (converges quickly) when the dominant eigenvalue is clearly dominant. --- # The Power Method - It works by starting with an initial vector `\(\mathbf x_0\)`, transforming to `\(\mathbf x_1=\mathbf A\mathbf x_0\)`, transforming `\(\mathbf x_1\)` to `\(\mathbf x_2=\mathbf A\mathbf x_1\)`, etc. `$$\mathbf x_1=\mathbf A\mathbf x_0,\quad \mathbf x_2=\mathbf A\mathbf x_1=\mathbf A^2\mathbf x_0,\ \cdots,\ \mathbf x_k=\mathbf A^k\mathbf x_0.$$` - As `\(k\rightarrow \infty\)`, the direction of `\(\mathbf x_k\)` converges to `\(\mathbf q_1\)`. --- # Estimation of eigenvalues So once we have an eigenvector estimate we can quickly estimate the corresponding eigenvalue `$$\mathbf A\mathbf x=\lambda\mathbf x\ \Rightarrow\ \mathbf x^\top\mathbf A\mathbf x=\lambda\,\mathbf x^\top\mathbf x \ \Rightarrow\ \lambda=\frac{\mathbf x^\top\mathbf A\mathbf x}{\mathbf x^\top\mathbf x} =\frac{\mathbf x^\top\mathbf A\mathbf x}{\lVert\mathbf x\rVert_2^2} =\mathbf q^\top\mathbf A\mathbf q,$$` where `\(\mathbf q=\mathbf x/\lVert\mathbf x\rVert_2\)` is the normalized vector; `\(\mathbf q^\top\mathbf A\mathbf q\)` is the **Rayleigh quotient**. --- # Practical Power method (normalised) Since the components of `\(\mathbf x_k\)` just get larger and larger as the Power Method iterates, and we really just want to know the direction of `\(\mathbf q_1\)` (not its magnitude), we normalise each `\(\mathbf x_k\)`. So the Power Method can be summarised as: - Set initial vector `\(\mathbf q_0 = \mathbf x_0/\lVert\mathbf x_0\rVert_2\)`. - Repeat for `\(k\)`: - Compute `\(\mathbf x_k = \mathbf A\mathbf q_{k-1}\)` - Normalise `\(\mathbf q_k = \mathbf x_k/\lVert\mathbf x_k\rVert_2\)` - Estimate `\(\lambda_k = \mathbf q_k^\top\mathbf A\mathbf q_k\)` --- # Comments - The Power method is not expected to converge if the matrix `\(\mathbf A\)` is not diagonalisable. - Convergence rate depends on how dominant `\(\lambda_1\)` is. - Google uses it to calculate the PageRank and Twitter uses it to show users recommendations of who to follow. - And for non-dominant eigenvalues/vectors? --- # Lab session Now use the power method to redo your google pagerank problem. --- # QR algorithm This algorithm computes an upper triangular matrix `\(\mathbf S\)` and an orthogonal matrix `\(\mathbf Q\)` such that `\(\mathbf A=\mathbf Q\mathbf S\mathbf Q^\top\)` is the Schur decomposition of `\(\mathbf A\)`. 1. Set `\(\mathbf A_0 := \mathbf A\)`. 2. for `\(k = 1, 2, \cdots\)`, - `\(\mathbf A_{k-1} = \mathbf Q_{k-1}\mathbf R_{k-1}.\)` - `\(\mathbf A_k := \mathbf R_{k-1}\mathbf Q_{k-1}.\)` 3. Set `\(\mathbf S := \mathbf A_{\infty}.\)` --- # Lab session - Go to R code up QR algorithm. - Use QR algorithm on `\(\mathbf A = \left( \begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array}\right)\)`. --- class: inverse, center, middle # Wrap-up --- # The full picture One question ran through the lecture: you are handed a matrix --- what does it look like? Structure picks the method; diagnostics verify the result. .small[ | Task | Matrix structure | Method | Diagnostics | |:--|:--|:--|:--| | Solve `\(\mathbf A\mathbf x=\mathbf b\)` | general invertible | `solve(A, b)` | residual, `\(\kappa\)` | | Solve `\(\mathbf A\mathbf x=\mathbf b\)` | SPD | Cholesky + triangular solves | symmetry, PD, residual | | Least squares | full column rank | QR | numerical rank, residual | | Least squares | rank deficient / nearly so | pivoted QR or SVD | singular values, rank | | Low-rank approximation | any matrix | truncated SVD | singular values, error | | Dominant eigenpair | huge, sparse | power method | eigen-residual | ] --- # Further readings - Golub, G. H. and Van Loan, C. F. (2013). *Matrix Computations* (4th edition). Johns Hopkins University Press. - Higham, N. J. (2002). *Accuracy and Stability of Numerical Algorithms* (2nd edition). SIAM. - Gentle, J. E. (2009). *Computational Statistics*. Springer. - [Google's PageRank](https://yanfei.site/docs/bsc/pagerank.pdf); [Iterative Methods for Computing Eigenvalues and Eigenvectors](http://yanfei.site/docs/bsc/svdnum.pdf). When reading software documentation, distinguish between the mathematical definition, the teaching implementation, and the production algorithm (pivoting, shifts, scaling, sparse storage, error handling).