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**: key directions, matrix powers, stability 5. **SVD**: works for any matrix; rank, conditioning, low rank 6. **Applications**: PCA, image compression, pseudo-inverse 7. **Iterative methods**: power method, QR algorithm 8. **Networks and text**: PageRank, latent semantic analysis 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. If several right-hand sides share the same `\(\mathbf A\)`, factor once, then reuse the factorization for cheap triangular solves. --- # Check the residual Given a computed solution `\(\widehat{\mathbf x}\)`, check the residual vector `$$\mathbf r=\mathbf b-\mathbf A\widehat{\mathbf x}.$$` To compare across scales, use the relative residual $$ \frac{\lVert\mathbf A\widehat{\mathbf x}-\mathbf b\rVert_2} {\lVert\mathbf A\rVert_2\lVert\widehat{\mathbf x}\rVert_2+\lVert\mathbf b\rVert_2}. $$ ```r relative_residual(A, x_direct, b) #> [1] 8.0646e-17 ``` --- # Residual vs conditioning A small relative residual says: "the reported solution nearly satisfies the equation you handed the computer". The **algorithm** has done its job. It does **not** guarantee that the solution is stable under small perturbations of `\(\mathbf A\)` or `\(\mathbf b\)` --- that is a property of the **problem**, not the algorithm. -- - **Residual**: did the algorithm solve the given problem accurately? - **Condition number**: is the problem itself sensitive to perturbation? For now, treat `\(\kappa\)` as a sensitivity measure of the *problem*, and use `kappa(A, exact = TRUE)` in R. Its precise meaning will emerge once we have the SVD: `\(\kappa_2(\mathbf A)=\sigma_1/\sigma_p\)`. --- 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()`. --- # Example: covariance matrix of stock returns Using the built-in `EuStockMarkets` data (for matrix practice, not investment advice): .scroll-output[ ```r stock_prices <- as.matrix(EuStockMarkets) stock_returns <- diff(log(stock_prices)) Sigma <- cov(stock_returns) b_stock <- colMeans(stock_returns) R_stock <- chol(Sigma) x_chol <- backsolve(R_stock, forwardsolve(t(R_stock), b_stock)) x_solve <- solve(Sigma, b_stock) c( factorization_error = norm(Sigma - crossprod(R_stock), "F") / norm(Sigma, "F"), solution_difference = max(abs(x_chol - x_solve)), relative_residual = relative_residual(Sigma, x_chol, b_stock) ) #> factorization_error solution_difference relative_residual #> 0.0000e+00 1.7764e-15 1.7875e-17 ``` ] --- # Cholesky gives a stable log-determinant Multivariate normal likelihoods and Bayesian computations need `\(\log\det(\mathbf A)\)`. From the factorization, `$$\log\det(\mathbf A)=2\sum_{j=1}^p\log R_{jj}.$$` ```r c( log_det_chol = 2 * sum(log(diag(R_stock))), log_det_direct = as.numeric(determinant(Sigma, logarithm = TRUE)$modulus) ) #> log_det_chol log_det_direct #> -39.388 -39.388 ``` No eigenvalues, no determinant of a huge matrix --- just the diagonal of `\(\mathbf R\)`. --- # 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. Do **not** "fix" this by silently adding a large diagonal constant --- that changes the problem you are solving. --- # Lab session **Cholesky and residuals.** Build a `\(4\times4\)` SPD matrix; verify `\(\mathbf A=\mathbf R^\top\mathbf R\)`; solve via two triangular solves; report the relative residual and compare with `solve(A, b)`. Then make the matrix indefinite and interpret the `chol()` error. --- 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 data usually do not look like that: with `\(n\)` observations and `\(p\)` parameters, the design matrix `\(\mathbf X\)` is `\(n\times p\)` and `$$\mathbf y=\mathbf X\boldsymbol\beta$$` generally has **no exact solution**. We minimize the squared residual instead: `$$\min_{\boldsymbol\beta}\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2.$$` --- # The tempting route is a trap A seemingly natural route: form the **normal equations** `$$\mathbf X^\top\mathbf X\widehat{\boldsymbol\beta}=\mathbf X^\top\mathbf y,$$` note that `\(\mathbf X^\top\mathbf X\)` is SPD (full column rank), and hand them to the Cholesky routine you just learned. -- Mathematically correct. Numerically dangerous --- and the reason is one identity: `$$\kappa_2(\mathbf X^\top\mathbf X)=\kappa_2(\mathbf X)^2.$$` Forming `\(\mathbf X^\top\mathbf X\)` **squares** the conditioning. The problem is not Cholesky --- it is the matrix you feed it. --- # 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\)` orthogonal and `\(\mathbf R\)` upper triangular. 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 `\((a_1,\ldots,a_n)\)` can be turned into an orthonormal basis `\((q_1,\ldots,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> (In production code, QR is computed by Householder reflections or Givens rotations, which are more stable than classical Gram--Schmidt; the idea, however, is the same.) --- # Least squares via QR With the thin QR of an `\(n\times p\)` full-column-rank design matrix, `\(\mathbf X=\mathbf Q_1\mathbf R_1\)`, the least squares problem becomes a triangular system: `$$\mathbf R_1\widehat{\boldsymbol\beta}=\mathbf Q_1^\top\mathbf y,$$` solved by one back substitution. (For a square system `\(\mathbf A\mathbf x=\mathbf b\)`, the same manipulation gives `\(\mathbf R_1\mathbf x=\mathbf Q_1^\top\mathbf b\)`.) - No `\(\mathbf X^\top\mathbf X\)` is ever formed. - This is what `lm()`, `lm.fit()` and `qr.solve()` do internally (with pivoting). --- # 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, normal_eq = 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 .scroll-output[ ```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. --- # When it actually breaks .scroll-output[ ```r set.seed(2026) n <- 1000 Xw <- matrix(rnorm(n * 20), n, 20) yw <- rnorm(n) W <- cbind(Xw, Xw[, 1] + rnorm(n, sd = 1e-12)) # near-duplicate column solve(crossprod(W), crossprod(W, yw)) # normal equations #> [,1] #> [1,] 16.0636068 #> [2,] -0.0240405 #> [3,] -0.0113174 #> [4,] -0.0081561 #> [5,] 0.0303891 #> [6,] 0.0050347 #> [7,] -0.0399401 #> [8,] 0.0338903 #> [9,] 0.0046889 #> [10,] -0.0189772 #> [11,] -0.0192437 #> [12,] -0.0331134 #> [13,] 0.0069001 #> [14,] 0.0109173 #> [15,] -0.0036603 #> [16,] 0.0080537 #> [17,] -0.0329856 #> [18,] 0.0200623 #> [19,] 0.0335841 #> [20,] -0.0287181 #> [21,] -16.0704783 ``` ] .scroll-output[ ```r qr.coef(qr(W), yw) |> round(3) #> [1] -0.007 -0.024 -0.011 -0.008 0.030 0.005 -0.040 0.034 0.005 -0.019 -0.019 -0.033 0.007 #> [14] 0.011 -0.004 0.008 -0.033 0.020 0.034 -0.029 NA ``` ] A tiny perturbation makes `\(\mathbf W^\top\mathbf W\)` computationally singular, so `solve()` refuses. Pivoted QR **detects** the aliasing and flags the redundant coefficient as `NA` --- exactly the failure predicted by the `\(\kappa^2\)` identity. --- # Residual orthogonality: a free diagnostic At the least squares solution, the residual is orthogonal to the columns of `\(\mathbf X\)`: `$$\mathbf X^\top(\mathbf y-\mathbf X\widehat{\boldsymbol\beta})\approx\mathbf 0.$$` ```r r_hat <- as.vector(y - X %*% beta_qr) c( relative_fit_residual = sqrt(sum(r_hat^2)) / sqrt(sum(y^2)), orthogonality = sqrt(sum(crossprod(X, r_hat)^2)) / (norm(X, type = "2") * sqrt(sum(r_hat^2))) ) #> relative_fit_residual orthogonality #> 3.5018e-03 2.2608e-13 ``` Both should be at rounding level. Together with the condition number, they are the core diagnostics for least squares. --- # Practical defaults in R - Fitted regression: `lm()` or `lm.fit()` --- pivoted QR inside. - Solving with a general square matrix: `solve(A, b)`. - Least squares by hand: `qr.solve(X, y)`. - SPD system or log-determinant: `chol()` + `forwardsolve()`/`backsolve()`. - Rank deficient or nearly so: pivoted QR or SVD, plus singular values and numerical rank (coming next). Explicit inverses and hand-built normal equations are fine on paper, and rarely the right default in code. --- # Lab session **Ill-conditioned least squares.** With `longley`, verify `\(\kappa_2(\mathbf X^\top\mathbf X)=\kappa_2(\mathbf X)^2\)`; compare normal equations, QR and `lm.fit()`. Add a near-duplicate column and watch numerical rank and coefficient stability. --- class: inverse, center, middle # Eigenanalysis --- # What is `\(\mathbf A\mathbf x=\lambda\mathbf x\)` asking? - In general, a matrix acts on a vector by changing both its magnitude and its direction. - Some special directions are **stable**: the matrix only stretches/compresses (possibly flips) the vector along the same line. <center> <img src="./figs/eigen.png" height="180px"/> </center> - `\(|\lambda|>1\)`: stretched; `\(0<|\lambda|<1\)`: compressed; `\(\lambda<0\)`: flipped; `\(\lambda=0\)`: collapsed. - Real matrices may have complex eigenvalues (rotation enters); real **symmetric** matrices --- covariances! --- have real eigenvalues and orthogonal eigenvectors. --- # The eigendecomposition If `\(\mathbf A\)` has `\(n\)` linearly independent eigenvectors, it is diagonalizable: `$$\mathbf A=\mathbf Q\boldsymbol\Lambda\mathbf Q^{-1}.$$` Read it as three steps: change coordinates to the eigenbasis, scale each direction by `\(\lambda_i\)`, change back. If `\(\mathbf A\)` is real **symmetric** (covariance, correlation, many quadratic forms): - all `\(n\)` eigenvalues are real, and the eigenvectors are orthogonal; - `\(\mathbf A=\mathbf Q\boldsymbol\Lambda\mathbf Q^\top\)` with `\(\mathbf Q^{-1}=\mathbf Q^\top\)` (spectral decomposition). If no diagonalization exists, a triangularization always does (Schur decomposition, `\(\mathbf A=\mathbf Q\mathbf S\mathbf Q^\top\)`) --- what numerical software actually computes. --- # Eigenanalysis in R ```r A_eig <- matrix(c(2, 1, 1, 2), 2, 2, byrow = TRUE) eig <- eigen(A_eig) eig$values #> [1] 3 1 round(eig$vectors, 3) #> [,1] [,2] #> [1,] 0.707 -0.707 #> [2,] 0.707 0.707 v1 <- eig$vectors[, 1] round(cbind(A_v1 = A_eig %*% v1, lambda_v1 = eig$values[1] * v1), 3) #> lambda_v1 #> [1,] 2.121 2.121 #> [2,] 2.121 2.121 ``` The last two columns agree: `\(\mathbf A\mathbf v_1=\lambda_1\mathbf v_1\)`. Eigenvector **signs** are arbitrary ($\mathbf v$ and `\(-\mathbf v\)` span the same direction) --- never over-interpret a sign. --- # Why care? Matrix powers Consider `\(\mathbf A=\left(\begin{array}{cc}5 & -1\\-2 & 4 \end{array}\right)\)` and `\(\mathbf x=(1,-1)^\top\)`. What is `\(\mathbf A^{20}\mathbf x\)`? Diagonalization turns matrix powers into scalar powers: `$$\mathbf A^k=\mathbf Q\boldsymbol\Lambda^k\mathbf Q^{-1}.$$` In Markov chains, population and industry transition, and web ranking, a state is repeatedly updated by the same matrix; the long-run behavior is governed by the **largest** eigenvalues and their eigenvectors. --- # Why care? Inverses, stability, collinearity If `\(\mathbf A\)` is diagonalizable with no zero eigenvalue, `$$\mathbf A^{-1}=\mathbf Q\boldsymbol\Lambda^{-1}\mathbf Q^{-1}.$$` `\(\boldsymbol\Lambda^{-1}\)` reciprocates the eigenvalues, so: - eigenvalues **near zero** `\(\Rightarrow\)` huge reciprocals `\(\Rightarrow\)` results highly sensitive to rounding; - multicollinearity in regression *is* small eigenvalues of `\(\mathbf X^\top\mathbf X\)`. Eigenanalysis is not only a way to *compute*; it tells you whether a problem is *easy to compute*. More generally, `\(f(\mathbf A)=\mathbf Q f(\boldsymbol\Lambda)\mathbf Q^{-1}\)` defines matrix functions: eigenvalues describe the strength of the matrix along its key directions. --- 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, PCA, 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\in\mathbb R^{m\times m}\)` and `\(\mathbf V\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> Geometrically: rotate/reflect, stretch along orthogonal directions by `\(\sigma_j\)`, rotate/reflect again. --- # 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`. --- # Properties you should remember - The non-zero singular values of `\(\mathbf A\)` are the square roots of the non-zero eigenvalues of `\(\mathbf A^\top\mathbf A\)` and `\(\mathbf A\mathbf A^\top\)`. - `\(\operatorname{rank}(\mathbf A)\)` = number of non-zero singular values. - For full column rank matrices, the 2-norm **condition number** --- promised earlier --- is `$$\kappa_2(\mathbf A)=\frac{\sigma_1}{\sigma_p}.$$` Large `\(\kappa\)`: computations are sensitive to rounding and data perturbation. If `\(\sigma_p=0\)` the standard condition number is infinite; `\(\sigma_1/\sigma_r\)` computed on the nonzero singular values is an **effective** condition number --- a different quantity, report the threshold used. --- # Economy SVD and the outer-product form If `\(\operatorname{rank}(\mathbf A)=r\)`, only `\(r\)` singular values are positive and `$$\mathbf A=\mathbf U_r\mathbf D_r\mathbf V_r^\top =\sum_{j=1}^r\sigma_j\mathbf u_j\mathbf v_j^\top =\sigma_1\mathbf u_1\mathbf v_1^\top+\cdots+\sigma_r\mathbf u_r\mathbf v_r^\top.$$` <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\}.$$` .scroll-output[ ```r A_rd <- cbind(c(1, 2, 3, 4), c(2, 4, 6, 8)) # exactly collinear columns d <- svd(A_rd, nu = 0, nv = 0)$d c(singular_values = signif(d, 3), ratio = signif(d[2] / d[1], 3)) #> singular_values1 singular_values2 ratio #> 1.22e+01 8.13e-16 6.64e-17 ``` ] --- class: inverse, center, middle # Applications of SVD --- # PCA is SVD in disguise Center (and usually standardize) the data matrix `\(\mathbf X_c\in\mathbb R^{n\times p}\)`, take its thin SVD `\(\mathbf X_c=\mathbf U\mathbf D\mathbf V^\top\)`. Then: - `\(\mathbf V\)`'s columns are the **loading** directions; - `\(\mathbf U\mathbf D=\mathbf X_c\mathbf V\)` are the **scores**; - the `\(j\)`-th pc has sample variance `\(d_j^2/(n-1)\)`, and variance share `\(d_j^2\big/\sum_\ell d_\ell^2\)`. PCA is not a separate algorithm: it is the **statistical interpretation** of the SVD of a centered data matrix. --- # PCA example: US state indicators `state.x77` (1970s data --- algorithm practice only): income, illiteracy, life expectancy, murder rate, HS graduation. Units differ, so standardize. .scroll-output[ ```r socio <- state.x77[, c("Income", "Illiteracy", "Life Exp", "Murder", "HS Grad")] pca <- prcomp(socio, center = TRUE, scale. = TRUE) variance_share <- pca$sdev^2 / sum(pca$sdev^2) round(cbind(share = variance_share, cumulative = cumsum(variance_share)), 3) #> share cumulative #> [1,] 0.640 0.640 #> [2,] 0.188 0.828 #> [3,] 0.080 0.908 #> [4,] 0.062 0.969 #> [5,] 0.031 1.000 round(pca$rotation[, 1:2], 3) #> PC1 PC2 #> Income 0.347 0.732 #> Illiteracy -0.480 0.069 #> Life Exp 0.469 -0.324 #> Murder -0.459 0.492 #> HS Grad 0.467 0.336 ``` ] --- # Reading the output - PC1 loads with the same sign on income, life expectancy, HS graduation, and the opposite sign on illiteracy and murder: a **socioeconomic gradient**. - The overall **sign of a PC is arbitrary**; what matters is the direction and the pattern of oppositions. - The largest-variance direction is **not** automatically "development", "wellbeing", or a causal effect. - If PCA enters a prediction pipeline, centering/scaling/loadings must be estimated on training data only --- otherwise data leakage. --- # Verify: `prcomp()` is the SVD .scroll-output[ ```r Z <- scale(socio) svd_socio <- svd(Z) c( sdev_difference = max(abs(svd_socio$d / sqrt(nrow(Z) - 1) - pca$sdev)), loading_difference_up_to_sign = max( abs(abs(svd_socio$v) - abs(pca$rotation)) ), rank2_relative_error = norm( Z - svd_socio$u[, 1:2] %*% diag(svd_socio$d[1:2]) %*% t(svd_socio$v[, 1:2]), "F") / norm(Z, "F") ) #> sdev_difference loading_difference_up_to_sign rank2_relative_error #> 0.00000 0.00000 0.41486 ``` ] --- # PCA output: variance shares and state scores .scroll-output[ ```r old_par <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1)) plot(variance_share, type = "b", pch = 19, ylim = c(0, 1), xlab = "Component", ylab = "Variance share") lines(cumsum(variance_share), type = "b", pch = 1, lty = 2, col = "steelblue") plot(pca$x[, 1], pca$x[, 2], type = "n", xlab = "PC1 score", ylab = "PC2 score") abline(h = 0, v = 0, col = "gray80") text(pca$x[, 1], pca$x[, 2], labels = state.abb, cex = 0.6) ``` <img src="BSC-L2-matrices_files/figure-html/pca-plot-1.png" width="90%" style="display: block; margin: auto;" /> ```r par(old_par) ``` ] --- # Lab session **PCA.** Run `prcomp()` on `state.x77` on both scales; reproduce the standardized scores via SVD; compute the rank-2 reconstruction error. --- # Image compression = low-rank approximation A grayscale image is a matrix (an RGB image: three matrices). Storing `\(\mathbf A\)` takes `\(mn\)` numbers; storing the rank `\(k\)` factors takes `$$k(m+n+1)\quad\text{numbers per channel}.$$` Since `\(\lVert\mathbf A-\mathbf A_k\rVert_2=\sigma_{k+1}\)`, a fast decay of singular values means visually faithful reconstructions from few components. ```r img <- jpeg::readJPEG("./figs/ykang.jpg") s_ch <- lapply(1:3, function(j) svd(img[, , j])) k <- 116 rec <- sapply(s_ch, function(s) { s$u[, 1:k] %*% diag(s$d[1:k]) %*% t(s$v[, 1:k]) }, simplify = "array") ``` --- # Rank 200, 116, 32 <center> <img src="./figs/ykang_svd_rank_200.jpg" height="230px"/> <img src="./figs/ykang_svd_rank_116.jpg" height="230px"/> <img src="./figs/ykang_svd_rank_32.jpg" height="230px"/> </center> Higher rank: more detail, larger factors. Lower rank: smaller factors, visible loss. Report *which* metric you mean by "compression": visual quality, matrix error, factor count, and file size are four different things. --- # The pseudo-inverse For `\(\mathbf A=\mathbf U_r\mathbf D_r\mathbf V_r^\top\)`, the Moore--Penrose **pseudo-inverse** is `$$\mathbf A^+=\mathbf V_r\mathbf D_r^{-1}\mathbf U_r^\top.$$` It exists for **any** matrix, and solves `$$\min_{\mathbf x}\lVert\mathbf A\mathbf x-\mathbf b\rVert_2^2$$` giving the usual least squares solution when it is unique, and the **minimum-norm** solution otherwise. For least squares this says: a small `\(\sigma_j\)` contributes `\(1/\sigma_j\)` --- tiny singular values mean trouble, and the threshold you keep determines the answer. --- # Pseudo-inverse in action Back to the near-duplicate-column matrix `\(\mathbf W\)` where `solve()` refused and `qr.coef()` flagged `NA`: .scroll-output[ ```r svd_W <- svd(W) keep <- svd_W$d > 1e-10 W_ginv <- svd_W$v[, keep, drop = FALSE] %*% (1 / svd_W$d[keep] * t(svd_W$u[, keep, drop = FALSE])) c(numerical_rank = sum(keep), columns = ncol(W)) #> numerical_rank columns #> 20 21 round(W_ginv %*% yw, 3) #> [,1] #> [1,] -0.003 #> [2,] -0.024 #> [3,] -0.011 #> [4,] -0.008 #> [5,] 0.030 #> [6,] 0.005 #> [7,] -0.040 #> [8,] 0.034 #> [9,] 0.005 #> [10,] -0.019 #> [11,] -0.019 #> [12,] -0.033 #> [13,] 0.007 #> [14,] 0.011 #> [15,] -0.004 #> [16,] 0.008 #> [17,] -0.033 #> [18,] 0.020 #> [19,] 0.034 #> [20,] -0.029 #> [21,] -0.003 ``` ] --- # Many "models" are matrix problems Curve fitting with a quadratic polynomial, `\(y_i\approx\beta_0+\beta_1x_i+\beta_2x_i^2\)`, is a least squares problem with design rows `\((1,x_i,x_i^2)\)`: .scroll-output[ ```r set.seed(123) x <- seq(-1, 1, length.out = 30) yc <- 1 + 2 * x - 1.5 * x^2 + rnorm(30, sd = 0.15) Xq <- cbind(1, x, x^2) beta_hat <- qr.solve(Xq, yc) round(beta_hat, 3) #> x #> 0.984 1.946 -1.475 ``` ] --- # Many "models" are matrix problems ```r plot(x, yc, pch = 19, col = "gray40", xlab = "x", ylab = "y") curve(beta_hat[1] + beta_hat[2] * x + beta_hat[3] * x^2, add = TRUE, col = "steelblue", lwd = 2) ``` <img src="BSC-L2-matrices_files/figure-html/curve-fit-plot-1.png" width="50%" style="display: block; margin: auto;" /> With high-degree polynomials the design matrix becomes ill-conditioned: centering, standardization, and stable factorizations matter again. --- 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. Rule of this section: **iteration count is not evidence of convergence**; algorithms must report residuals and convergence status. --- # The power method Start with `\(\mathbf x_0\)` and repeatedly apply `\(\mathbf A\)`: `$$\mathbf x_k=\mathbf A^k\mathbf x_0.$$` Expand `\(\mathbf x_0=c_1\mathbf q_1+\cdots+c_n\mathbf q_n\)` in the eigenbasis with `\(|\lambda_1|>|\lambda_2|\geq\cdots\)`: `$$\mathbf A^k\mathbf x_0 =\lambda_1^k\left\{c_1\mathbf q_1+ \sum_{i=2}^n c_i\left(\frac{\lambda_i}{\lambda_1}\right)^k\mathbf q_i\right\} \;\rightarrow\; c_1\lambda_1^k\mathbf q_1.$$` Every term with `\(i\geq 2\)` decays geometrically at rate `\(|\lambda_2/\lambda_1|\)` --- the direction of `\(\mathbf x_k\)` converges to the dominant eigenvector `\(\mathbf q_1\)` (provided `\(c_1\neq 0\)`). --- # The normalized power method Since only the direction matters, normalize at every step; estimate the eigenvalue with the **Rayleigh quotient**: 1. Set `\(\mathbf q_0=\mathbf x_0/\lVert\mathbf x_0\rVert_2\)`. 2. For `\(k=1,2,\ldots\)`: - `\(\mathbf z_k=\mathbf A\mathbf q_{k-1}\)` - `\(\mathbf q_k=\mathbf z_k/\lVert\mathbf z_k\rVert_2\)` - `\(\lambda_k=\mathbf q_k^\top\mathbf A\mathbf q_k\)` - stop when the scaled eigen-residual `\(\lVert\mathbf A\mathbf q_k-\lambda_k\mathbf q_k\rVert\)` is small enough. Each step costs one matrix--vector product --- that is why sparse, huge problems are tractable. --- # Power method in R .scroll-output[ ```r A_pm <- matrix(c(2, -1, -1, 2), 2, 2, byrow = TRUE) set.seed(2026) q <- rnorm(2) q <- q / sqrt(sum(q^2)) for (k in 1:50) { z <- as.vector(A_pm %*% q) q <- z / sqrt(sum(z^2)) } lambda <- as.numeric(t(q) %*% A_pm %*% q) c( power_value = lambda, exact = eigen(A_pm, symmetric = TRUE)$values[1], eigen_residual = sqrt(sum((A_pm %*% q - lambda * q)^2)) ) #> power_value exact eigen_residual #> 3 3 0 round(q, 3) # direction (1, -1)/sqrt(2), up to sign #> [1] 0.707 -0.707 ``` ] --- # Comments - Convergence is fast when `\(\lambda_1\)` clearly dominates; slow when `\(|\lambda_2|\approx|\lambda_1|\)`. - Fails if the start vector is orthogonal to `\(\mathbf q_1\)` (here: an all-ones start is orthogonal to `\((1,-1)\)`!). Random starts avoid this with probability 1, but always report the residual. - Not applicable to non-diagonalizable matrices in general; it gives only the dominant pair. - Used at enormous scale: Google's PageRank and Twitter's follow recommendations. - For non-dominant eigenpairs, or all eigenvalues: the QR algorithm. --- # The QR algorithm Idea: iterate QR factorizations to drive the matrix toward (quasi-)triangular, eigenvalue-revealing form. 1. Set `\(\mathbf A_0:=\mathbf A\)`. 2. For `\(k=1,2,\ldots\)`: - factor `\(\mathbf A_{k-1}=\mathbf Q_{k-1}\mathbf R_{k-1}\)`; - set `\(\mathbf A_k:=\mathbf R_{k-1}\mathbf Q_{k-1}\)`. Since `\(\mathbf A_k=\mathbf Q_{k-1}^\top\mathbf A_{k-1}\mathbf Q_{k-1}\)`, every iterate is **similar** to `\(\mathbf A\)`: eigenvalues never change. For symmetric `\(\mathbf A\)`, the iterates converge to a diagonal matrix of eigenvalues; for general real matrices, to real Schur form (quasi-triangular, `\(2\times2\)` blocks for complex conjugate pairs). --- # QR algorithm in R .scroll-output[ ```r Ak <- matrix(c(4, 1, 1, 3), 2, 2, byrow = TRUE) for (k in 1:30) { # no column pivoting: keeps the iteration a similarity transform fa <- qr(Ak, LAPACK = FALSE, tol = 0) Ak <- qr.R(fa) %*% qr.Q(fa) } c( QR_iteration = sort(diag(Ak), decreasing = TRUE), eigen = sort(eigen(matrix(c(4, 1, 1, 3), 2, 2, byrow = TRUE))$values, decreasing = TRUE) ) #> QR_iteration1 QR_iteration2 eigen1 eigen2 #> 4.618 2.382 4.618 2.382 ``` ] "QR the factorization" and "QR the eigenvalue algorithm" share one tool but solve different tasks. Production code adds shifts and Hessenberg reduction --- never hand-roll your eigenvalue solver. --- # Lab session - Go to R and code up the QR algorithm. - Use the QR algorithm on $$ \mathbf A= \begin{pmatrix} 1 & 2\\ 3 & 4 \end{pmatrix}. $$ - **Power method diagnostics.** Choose an initial vector orthogonal to the dominant eigenvector and watch the failure; then compare the power method and QR algorithm on a matrix with nearly equal leading eigenvalues. --- class: inverse, center, middle # PageRank: # ranking the web with eigenvectors --- # Application of eigenanalysis: Google <center> <img src="./figs/google.png" height="380px"/> </center> --- # The idea: links as weighted votes A page is important if **important** pages link to it. If page `\(i\)` has `\(N_i\)` outlinks, each link casts a vote of weight `\(1/N_i\)`. Write the link structure as a transition matrix `\(\mathbf B\)`: `$$B_{ij}:=\left\{\begin{array}{ll} 1/N_{i} & \text{if page } i \text{ links to page } j,\\ 0 & \text{otherwise.} \end{array}\right.$$` At a stationary importance distribution, one step of link-following must not change it: `$$\mathbf r=\mathbf B^\top\mathbf r.$$` The PageRank vector `\(\mathbf r\)` is the eigenvector of `\(\mathbf B^\top\)` for eigenvalue `\(1\)` --- a ranking problem becomes an eigenvalue problem. (The real web: billions of pages, extremely sparse matrix.) --- # A small example .pull-left[ Five pages with these links: <br> <img src="./figs/googleeg.png" height="330px"/> ] -- .pull-right[ `$$\mathbf B=\left(\begin{array}{ccccc} 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)$$` The ranking of page 4, say, collects votes: `$$r_4=\tfrac13 r_1+\tfrac13 r_2+\tfrac12 r_3.$$` `\(n=5\)` linear equations; rescale so `\(\sum_i r_i=1\)`. ] --- # PageRank in R .scroll-output[ ```r A_pr <- matrix(c( 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0 ), 5, 5, byrow = TRUE) B <- sweep(A_pr, 1, rowSums(A_pr), "/") er <- eigen(t(B)) r <- Re(er$vectors[, which.min(Mod(er$values - 1))]) if (sum(r) < 0) r <- -r r <- r / sum(r) round(r, 3) #> [1] 0.150 0.050 0.300 0.217 0.283 ``` ] --- # PageRank from the power-method view A surfer starts uniformly at random, `\(\mathbf r_0=(1/5,\ldots,1/5)^\top\)`, and follows a random outlink each step: `$$\mathbf r_1=\mathbf B^\top\mathbf r_0,\quad \mathbf r_2=\mathbf B^\top\mathbf r_1,\quad\ldots$$` This is **exactly the power method** applied to `\(\mathbf B^\top\)`. ```r r_iter <- rep(1 / 5, 5) for (k in 1:50) r_iter <- as.vector(t(B) %*% r_iter) c(power_iter = round(r_iter, 3), eigen = round(r, 3)) #> power_iter1 power_iter2 power_iter3 power_iter4 power_iter5 eigen1 eigen2 eigen3 #> 0.150 0.050 0.300 0.217 0.283 0.150 0.050 0.300 #> eigen4 eigen5 #> 0.217 0.283 ``` --- # Lab session **PageRank.** Add a dangling page and a two-cycle to the 5-page example; repair dangling rows; compare `\(\alpha=1,0.85,0.5\)` (iterations, residual, ranking). --- class: inverse, center, middle # Wrap-up --- # The full picture | Task | Matrix structure | Method | Diagnostics | |:---|:---|:---|:---| | Solve `\(\mathbf A\mathbf x=\mathbf b\)` | general invertible | `solve(A, b)` | relative residual, `\(\kappa\)` | | Solve `\(\mathbf A\mathbf x=\mathbf b\)` | symmetric positive definite | Cholesky + triangular solves | symmetry, PD, residual | | Minimize `\(\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2\)` | full column rank | QR | numerical rank, residual | | Least squares | rank deficient / nearly so | pivoted QR or SVD | singular values, rank | | Dimension reduction | centered data matrix | SVD / PCA | variance shares, error | | Dominant eigenpair at scale | huge, sparse | power method | eigen-residual, convergence | --- # 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).