Today's Topics

• Singular value decomposition

• Sparse matrices

• Case study: regression and PCA for sparse data

Singular Value Decomposition

• Singular value decomposition (SVD) can be viewed as an extension to the eigenvalue decomposition of symmetric matrices

• SVD applies to general rectangular matrices

• Wide applications in dimension reduction, image compression, recommender system, etc.

SVD Theorem ^[2]

Every matrix $A\in R^{m\times n}$ of rank $r$ can be decomposed as

$$A=U\Sigma V^{\prime }=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{\Sigma }_1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}{V}_1^{\prime }\\ V_2^{\prime }\end{array}\right),$$

where $U_{m\times m}=(u_1,\dots ,u_m)$ and $V_{n\times n}=(v_1,\dots ,v_n)$ are orthogonal matrices, $U_1\in R^{m\times r}$, $V_1\in R^{n\times r}$, and ${\Sigma }_1=diag({\sigma }_1,\dots ,{\sigma }_r)$ is a nonnegative diagonal matrix.

${\sigma }_1\ge \cdots \ge {\sigma }_r>0$ are called singular values of $A$. $u_i,i=1,\dots ,m$ and $v_j,j=1,\dots ,n$ are called left and right singular vectors of $A$, respectively.

Compact SVD

According to the factorization theorem, we can also write $$A=U_1{\Sigma }_1{V}_1^{\prime },$$ where $A$ is $m\times n$, $U_1$ is $m\times r$, ${\Sigma }_1$ is $r\times r$, and $V_1$ is $n\times r$.

$U_1$ and $V_1$ are column-orthonormal matrices.

This is typically called the compact SVD.

Relation to Eigenvalues

Suppose we have the (compact) SVD $A=U\Sigma V^{\prime }$, then

• $A^{\prime }A=V\Sigma U^{\prime }U\Sigma V^{\prime }=V{\Sigma }^2{V}^{\prime }$

• $A{A}^{\prime }=U\Sigma V^{\prime }V\Sigma U^{\prime }=U{\Sigma }^2{U}^{\prime }$

• Singular values of $A$ are square roots of positive eigenvalues of $A^{\prime }A$ and $A{A}^{\prime }$

For generality, we usually also allow singular values to be zero, by extending the $U$, $\Sigma$, and $V$ matrices to proper dimensions.

Partial SVD

If we only extract the largest $k$ singular values and associated singular vectors of $A$, we call the result $U_{\left(k\right)}{\Sigma }_{\left(k\right)}{V}_{\left(k\right)}^{\prime }$ the partial SVD of $A$.

• ${\Sigma }_{\left(k\right)}=diag({\sigma }_1,\dots ,{\sigma }_k)$, ${\sigma }_1\ge \cdots \ge {\sigma }_k$

• $U_{\left(k\right)}=(u_1,\dots ,u_k)$

• $V_{\left(k\right)}=(v_1,\dots ,v_k)$

Relation to Low-rank Approximation

Partial SVD is useful due to an interesting theorem.

Suppose that a matrix $A_{m\times n}$, $m\ge n$ has SVD $A=U\Sigma V^{\prime }$, where $\Sigma =diag({\sigma }_1,\dots ,{\sigma }_n)$ and ${\sigma }_1\ge \cdots \ge {\sigma }_n\ge 0$. Then the best rank-k approximation to $A$ in the Frobenius norm is given by

$$A_k=\sum _{i=1}^k{\sigma }_i{u}_i{v}_i^{\prime }.$$

The claim also holds with the operator norm.

SVD Computation

Similar to eigenvalue computation, there are different use cases of SVD:

• Computing all singular values (full SVD)

• Computing the largest $k$ singular values (partial SVD)

It is not common to request the smallest singular values.

Computing Full SVD

Similar to full eigenvalue decomposition, computing full SVD also typically has two steps:

1. Looking for orthogonal matrices $U_1$ and $V_1$ such that $U_1^{\prime }A{V}_1$ is bidiagonal $$U_1^{\prime }A{V}_1=\left[\begin{array}{c}B\\ 0\end{array}\right],B=\left[\begin{array}{ccccc}{\rho }_1& {\theta }_2& & & \\ & {\rho }_2& {\theta }_3& & \\ & & \ddots & \ddots & \\ & & & {\rho }_{n-1}& {\theta }_n\\ & & & & {\rho }_n\end{array}\right]$$

2. Applying the QR algorithm to the matrix $B$ to compute all singular values

The Golub-Kahan Householder (Bidiagonalization) Algorithm ^[2]

Any $A_{m\times n}$ with $m\ge n$ can be reduced to $$U_1^{\prime }A{V}_1=\left[\begin{array}{c}B\\ 0\end{array}\right]$$ by orthogonal matrices $U_1\in R^{m\times m}$ and $V_1\in R^{n\times n}$. $U_1$ and $V_1$ are product of simple orthogonal matrices:

• $U_1=Q_1{Q}_2\cdots Q_n$
• $V_1=P_0{P}_1\cdots P_{n-2}$

The singular values of $B$ are the same as those of $A$.

The QR Algorithm

Since singular values of $B$ are square roots of eigenvalues of $B^{\prime }B$, the QR algorithm can proceed as follows:

• Compute the QR decomposition: $B^{\left(k\right)\prime }{B}^{\left(k\right)}=Q_k{R}_k$

• Update $B^{\left(k\right)}$ to $B^{(k+1)}$ such that $B^{(k+1)\prime }{B}^{(k+1)}=R_k{Q}_k=Q_k^{\prime }{B}^{\left(k\right)\prime }{B}^{\left(k\right)}{Q}_k$

$B^{\left(k\right)}$ is bidiagonal for all $k$.

The QR decomposition and the update of $B^{\left(k\right)}$ would use the special structure of $B^{\left(k\right)}$, so $B^{\left(k\right)\prime }{B}^{\left(k\right)}$ will not be explicitly formed.

Advanced Algorithms

• In practice, there are many improvements made to this basic method

• A good summary of those advanced algorithms can be found at https://www.cs.utexas.edu/users/inderjit /public_papers/HLA_SVD.pdf

Computing Partial SVD

• Partial SVD seeks the largest $k$ singular values of a matrix $A_{m\times n}$

• Can be reduced to computing the largest $k$ eigenvalues of $A^{\prime }A$ (if $m\ge n$) or $A{A}^{\prime }$ (if $m<n$)

• Using the power method or other iterative methods, and assuming $m\ge n$, we need to realize the operation $v\to A^{\prime }Av$

• Be cautious of the computing order! You should compute $x\leftarrow Av$ first and then do $y\leftarrow A^{\prime }x$

Advanced Algorithms

• Another algorithm to compute partial SVD is the augmented implicitly restarted Lanczos bidiagonalization method ^[4]

• It works on the original matrix $A$ instead of $A^{\prime }A$

• Potentially more numerical stable than computing the eigenvalues of $A^{\prime }A$

Implementation

• Simply call svd() for full (but compact) SVD
• You can specify how many left/right singular vectors to return

set.seed(123)
n = 1000
p = 500
x = matrix(rnorm(n * p), n, p)
str(svd(x))

## List of 3
##  $ d: num [1:500] 53.6 53.3 53.3 53 52.6 ...
##  $ u: num [1:1000, 1:500] 0.00265 -0.03029 0.00451 -0.00763 -0.0034 ...
##  $ v: num [1:500, 1:500] 0.01699 0.00386 0.00971 -0.01931 0.01455 ...

str(svd(x, nu = 10, nv = 10))

## List of 3
##  $ d: num [1:500] 53.6 53.3 53.3 53 52.6 ...
##  $ u: num [1:1000, 1:10] 0.00265 -0.03029 0.00451 -0.00763 -0.0034 ...
##  $ v: num [1:500, 1:10] 0.01699 0.00386 0.00971 -0.01931 0.01455 ...

Benchmark

• Use nu = 0 and nv = 0 if you do not need singular vectors

library(dplyr)
bench::mark(
    svd(x),
    svd(x, nu = 100, nv = 100),
    svd(x, nu = 1, nv = 1),
    svd(x, nu = 0, nv = 0),
    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 svd(x)                        785ms    795ms
## 2 svd(x, nu = 100, nv = 100)    784ms    786ms
## 3 svd(x, nu = 1, nv = 1)        786ms    788ms
## 4 svd(x, nu = 0, nv = 0)        272ms    275ms

Implementation

• There are various options for partial SVD
  • irlba::irlba
  • RSpectra::svds
  • svd::propack.svd, svd::trlan.svd, and svd::ztrlan.svd

res1 = irlba::irlba(x, 10)
res1$d

##  [1] 53.56025 53.30589 53.25349 53.03145 52.56782 52.30170 52.13693 51.85140
##  [9] 51.64052 51.42326

res2 = RSpectra::svds(x, k = 10)
res2$d

##  [1] 53.56025 53.30589 53.25349 53.03145 52.56782 52.30170 52.13693 51.85140
##  [9] 51.64052 51.42326

Implementation

res3 = svd::propack.svd(x, neig = 10)
res3$d

##  [1] 53.56025 53.30589 53.25349 53.03145 52.56782 52.30170 52.13693 51.85140
##  [9] 51.64052 51.42326

res4 = svd::trlan.svd(x, neig = 10)
res4$d

##  [1] 53.56025 53.30589 53.25349 53.03145 52.56782 52.30170 52.13693 51.85140
##  [9] 51.64052 51.42326

res5 = svd::ztrlan.svd(x, neig = 10)
res5$d

##  [1] 53.56025 53.30589 53.25349 53.03145 52.56782 52.30170 52.13693 51.85140
##  [9] 51.64052 51.42326

Benchmark

• irlba::irlba, RSpectra::svds, and svd::propack.svd are roughly at the same level

bench::mark(
    irlba::irlba(x, 10, tol = 1e-8),
    RSpectra::svds(x, k = 10, opts = list(tol = 1e-8)),
    svd::propack.svd(x, neig = 10, opts = list(tol = 1e-8)),
    svd::trlan.svd(x, neig = 10, opts = list(tol = 1e-8)),
    svd::ztrlan.svd(x, neig = 10, opts = list(tol = 1e-8)),
    min_iterations = 10, max_iterations = 10, check = FALSE
) %>% select(expression, min, median)

## # A tibble: 5 × 3
##   expression                                                    min   median
##   <bch:expr>                                               <bch:tm> <bch:tm>
## 1 irlba::irlba(x, 10, tol = 1e-08)                           97.7ms  109.8ms
## 2 RSpectra::svds(x, k = 10, opts = list(tol = 1e-08))        88.4ms   89.8ms
## 3 svd::propack.svd(x, neig = 10, opts = list(tol = 1e-08))   99.2ms  102.3ms
## 4 svd::trlan.svd(x, neig = 10, opts = list(tol = 1e-08))    241.6ms  248.7ms
## 5 svd::ztrlan.svd(x, neig = 10, opts = list(tol = 1e-08))   400.8ms  406.1ms

Implementation

• We will see how partial SVD can be used to efficiently compute PCA of sparse data

Sparse Matrix

A matrix that has very few nonzero elements.

• Mathematically, not very different from dense matrices
• But sparsity can make a huge impact on computing

Problems

We are mostly interested in the following problems:

• How to store sparse data/matrices

• How to implement basic matrix operations on sparse matrices

• How to do statistical computing on sparse matrices

Storage

Two major storage formats:

• Coordinate (COO) format

• Compressed sparse column (CSC) format

COO Format

• Simplest storage scheme

• Row indices, column indices, nonzero elements

$$A=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 3& 0& 0& 5& 0\\ 0& 0& 7& 8& 9\\ 0& 0& 2& 1& 0\\ 0& 0& 0& 0& 0\end{array}\right]$$

$$\begin{array}{cc}i& =\left[\begin{array}{cccccccc}1& 2& 3& 4& 2& 3& 4& 3\end{array}\right]\\ j& =\left[\begin{array}{cccccccc}1& 1& 3& 3& 4& 4& 4& 5\end{array}\right]\\ x& =\left[\begin{array}{cccccccc}1& 3& 7& 2& 5& 8& 1& 9\end{array}\right]\end{array}$$

COO Format

• For an $m\times n$ matrix with $nnz$ nonzero elements

• COO format needs to store $3\cdot nnz$ numbers

• However, the $j$ vector contains redundant information

• The CSC format is a more economic scheme to represent sparse matrices

CSC Format

• The $x$ vector stores elements in $A$ column by column

• Use a "pointer vector" $p$, where $p_{j+1}-p_j$ is the number of nonzero elements in column $j$

$$A=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 3& 0& 0& 5& 0\\ 0& 0& 7& 8& 9\\ 0& 0& 2& 1& 0\\ 0& 0& 0& 0& 0\end{array}\right]$$

$$\begin{array}{cc}i& =\left[\begin{array}{cccccccc}1& 2& 3& 4& 2& 3& 4& 3\end{array}\right]\\ j& =\left[\begin{array}{cccccccc}1& 1& 3& 3& 4& 4& 4& 5\end{array}\right]\\ p& =\left[\begin{array}{cccccc}0& 2& 2& 4& 7& 8\end{array}\right]\\ x& =\left[\begin{array}{cccccccc}1& 3& 7& 2& 5& 8& 1& 9\end{array}\right]\end{array}$$

CSC Format

• The CSC format only needs to store the $i$, $p$, and $x$ vectors

• $2\cdot nnz+n+1$ numbers in total

• Question: If we also know that the sparse matrix is binary, i.e., every element is either 0 or 1, then what is a good scheme for storage?

Basic Operations

• Elementwise operations are relatively simple

• One of the most important operations for linear systems and eigenvalue computation is the matrix-vector multiplication

• $v\to Av$ and/or $v\to A^{\prime }v$

Computing $A^{\prime }v$

• The key point is to extract each column of $A$

• Each column can be viewed as a sparse vector

$$A=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 3& 0& 0& 5& 0\\ 0& 0& 7& 8& 9\\ 0& 0& 2& 1& 0\\ 0& 0& 0& 0& 0\end{array}\right]=\left[\begin{array}{ccccc}{a}_1& a_2& a_3& a_4& a_5\end{array}\right]$$

$$A^{\prime }v={\left[\begin{array}{ccccc}{a}_1^{\prime }v& a_2^{\prime }v& a_3^{\prime }v& a_4^{\prime }v& a_5^{\prime }v\end{array}\right]}^{\prime }$$

Extract $a_j$

• $a_j$ has $(p_{j+1}-p_j)$ nonzero elements

• The $(p_j+1)$-th element of $x$ is the first nonzero element of $a_j$
• The $(p_j+1)$-th element of $i$ is the location of this element in the sparse vector

• The $(p_j+2)$-th element of $x$ is the second nonzero element of $a_j$
• The $(p_j+2)$-th element of $i$ is the location of this element in the sparse vector

• ...

Extract $a_j$

Computing $A^{\prime }v$

• Similarly, we can get $$\begin{array}{cc}{a}_1& =\{1\to 1,2\to 3\}\\ a_2& =\left\{\right\}\\ a_3& =\{3\to 7,4\to 2\}\\ a_4& =\{2\to 5,3\to 8,4\to 1\}\\ a_5& =\{3\to 9\}\end{array}$$

• Then the inner products only involve nonzero elements $$\begin{array}{cc}{a}_1^{\prime }v& =1\cdot v_1+3\cdot v_2\\ a_2^{\prime }v& =0\\ a_3^{\prime }v& =7\cdot v_3+2\cdot v_4\\ & \cdots \end{array}$$

Computing $Av$

• We can again use the sparse vectors $a_j$ to compute $Av$

$$A=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 3& 0& 0& 5& 0\\ 0& 0& 7& 8& 9\\ 0& 0& 2& 1& 0\\ 0& 0& 0& 0& 0\end{array}\right]=\left[\begin{array}{ccccc}{a}_1& a_2& a_3& a_4& a_5\end{array}\right]$$

$$Av=v_1{a}_1+\cdots +v_5{a}_5$$

• Each $v_j{a}_j$ is still a sparse vector formed by multiplying each nonzero element of $a_j$ by $v_j$

Computing $Av$

• Therefore, computing $Av$ reduces to:

  1. Initialize $r\leftarrow 0_m$

  2. For $j=1,\dots ,n$ do $r\leftarrow r+v_j{a}_j$

• For example, $a_3=\{3\to 7,4\to 2\}$

• Then $r\leftarrow r+v_3{a}_3$ expands to

$$\begin{array}{cc}{r}_3& \leftarrow r_3+7\cdot v_3\\ r_4& \leftarrow r_4+2\cdot v_4\end{array}$$

Computational Complexity

• In general, for matrix $A_{m\times n}$ with $nnz$ nonzero elements

• The complexity of $v\to Av$ and $v\to A^{\prime }v$ are both $O\left(nnz\right)$

• This makes many iterative methods very efficient on sparse matrices

Implementation

• The Matrix R package has nice support for different kinds of sparse matrices
• A simple way to construct sparse matrices is to create from an (i, j, x) triplet (COO format)
• The result is typically in CSC format (dgCMatrix)

library(Matrix)
i = c(1, 2, 3, 4, 2, 3, 4, 3)
j = c(1, 1, 3, 3, 4, 4, 4, 5)
x = c(1, 3, 7, 2, 5, 8, 1, 9)
xsp = sparseMatrix(i = i, j = j, x = x, dims = c(5, 5))
xsp

## 5 x 5 sparse Matrix of class "dgCMatrix"
##               
## [1,] 1 . . . .
## [2,] 3 . . 5 .
## [3,] . . 7 8 9
## [4,] . . 2 1 .
## [5,] . . . . .

Implementation

• Note that internally, the $i$ vector uses 0-based indices

xsp

## 5 x 5 sparse Matrix of class "dgCMatrix"
##               
## [1,] 1 . . . .
## [2,] 3 . . 5 .
## [3,] . . 7 8 9
## [4,] . . 2 1 .
## [5,] . . . . .

str(xsp)

## Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   ..@ i       : int [1:8] 0 1 2 3 1 2 3 2
##   ..@ p       : int [1:6] 0 2 2 4 7 8
##   ..@ Dim     : int [1:2] 5 5
##   ..@ Dimnames:List of 2
##   .. ..$ : NULL
##   .. ..$ : NULL
##   ..@ x       : num [1:8] 1 3 7 2 5 8 1 9
##   ..@ factors : list()

Implementation

• Another way to construct sparse matrices is to convert from dense matrices

set.seed(123)
n = 1000
nnz = floor(0.01 * n^2)
xdata = numeric(n^2)
xdata[sample(n^2, nnz)] = rnorm(nnz)
x = matrix(xdata, n, n)
class(x)

## [1] "matrix" "array"

xsp = as(x, "sparseMatrix")
class(xsp)

## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"

Benchmark

• Matrix multiplication operations are already defined for sparse matrices

bench::mark(
    x %*% x, crossprod(x),
    xsp %*% xsp, crossprod(xsp),
    min_iterations = 10, max_iterations = 10, check = FALSE
) %>% select(expression, min, median)

## # A tibble: 4 × 3
##   expression          min   median
##   <bch:expr>     <bch:tm> <bch:tm>
## 1 x %*% x        488.19ms 493.53ms
## 2 crossprod(x)    482.8ms 493.61ms
## 3 xsp %*% xsp       1.3ms   1.53ms
## 4 crossprod(xsp)   1.48ms   1.58ms

Direct Methods

• Direct methods for sparse matrices are typically difficult problems

• Preservation of sparsity is a big challenge

• If $A$ is sparse, $A=LU$, $L$ and $U$ are not necessarily sparse

• Therefore, proper ordering of rows/columns of $A$ is critical

• With some permutation matrices $P$ and $Q$, $PAQ=\tilde{L}\tilde{U}$, $\tilde{L}$ and $\tilde{U}$ may be much more sparse than $L$ and $U$

Direct Methods

• There are indeed sparse decomposition methods, e.g.,

  • Sparse LU decomposition

  • Sparse Cholesky decomposition

• They are typically very sophisticated

• Software support such as SuiteSparse

Iterative Methods

• We are mostly interested in iterative methods for sparse matrices

• Small memory cost

• Computationally efficient

• Easy to implement

Core Idea

• Iterative methods that rely on $v\to Av$, $v\to A^{\prime }v$, or some other minor operations are readily available for sparse matrices

• Conjugate gradient method

• Iterative eigenvalue algorithms such as power method

• Partial SVD

Statistical Computing on Sparse Matrix

• Regression

• Ridge regression

• PCA

• Lasso (preview)

Regression on Sparse Data

• Suppose the data matrix $X$ is sparse

• Target is $\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$

• However, $X^{\prime }X$ may no longer be sparse

• Also, $X$ typically does not have a "specialized" QR decomposition

• Therefore, there is no obvious benefit of the sparsity (except for computing $X^{\prime }X$)

Regression on Sparse Data

• In contrast, CG can effectively utilize the sparsity

• CG solves $Ax=b$, where $A=X^{\prime }X$ and $b=X^{\prime }Y$

• $A$ is never explictly formed

• CG requires the operation $v\to Av=X^{\prime }Xv$

• Compute two sparse matrix-vector multiplications

  1. $u\leftarrow Xv$

  2. $w\leftarrow X^{\prime }u$

Implementation

• First simulate sparse data and use the dense solver to get the result

library(Matrix)
set.seed(123)
n = 10000
p = 500
xsp = rsparsematrix(n, p, density = 0.001)
x = as.matrix(xsp)
y = rnorm(n)
bhat = solve(crossprod(x), crossprod(x, y))
bhat[1:5]

## [1] -0.003835525 -0.211645950  0.102156056  0.294562285  0.031367565

Implementation

• Slight modification of the CG code

reg_cg = function(Xsp, y, x0 = rep(0, ncol(Xsp)), eps = 1e-6)
{
    b = as.numeric(crossprod(Xsp, y))
    x = x0
    p = r = b - as.numeric(crossprod(Xsp, Xsp %*% x))
    r2 = sum(r^2)
    errs = c()
    for(i in seq_along(b))
    {
        Ap = as.numeric(crossprod(Xsp, Xsp %*% 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(beta_hat = x, errs = errs, niter = i)
}

Implementation

• Verify the result

cg = reg_cg(xsp, y, eps = 1e-6)
bhat_cg = cg$beta_hat
cg$niter

## [1] 56

bhat_cg[1:5]

## [1] -0.003835524 -0.211645950  0.102156055  0.294562283  0.031367563

max(abs(bhat_cg - bhat))

## [1] 1.548608e-07

Benchmark

• Significant speedup over dense operations

bench::mark(
    solve(crossprod(x), crossprod(x, y)),
    reg_cg(xsp, y),
    min_iterations = 3, max_iterations = 10, check = FALSE
) %>% select(expression, min, median)

## # A tibble: 2 × 3
##   expression                                min   median
##   <bch:expr>                           <bch:tm> <bch:tm>
## 1 solve(crossprod(x), crossprod(x, y))    1.37s    1.37s
## 2 reg_cg(xsp, y)                         3.83ms   5.06ms

Ridge Regression on Sparse Data

• Ridge regression is similar, but now $A=X^{\prime }X+\lambda I$ and $b=X^{\prime }Y$

• Again, $A$ should not be explictly formed

• Now we need the operation $v\to Av=(X^{\prime }X+\lambda I)v$

• Compute two sparse matrix-vector multiplications and one vector addition

  1. $u\leftarrow Xv$

  2. $w\leftarrow X^{\prime }u$

  3. $z\leftarrow w+\lambda v$

Excercises

1. Implement the ridge regression for a sparse $X$ with $n=10000$ and $p=500$

2. What if $p\gg n$? For example, $X$ is a sparse matrix with $n=500$ and $p=10000$

3. We can find that the majority part of the CG code is unchanged in different problems. Design an implementation of CG that is as general as possible. Hint: pass a function Ax as the first argument of CG, where Ax has the signature

   Ax = function(x, args) {...}

   Ax should implement the operation $v\to Av$, and args contains information about $A.$

PCA on Sparse Data

• Another common task for large data is computing PCA

• Equivalent to computing the largest $k$ eigenvalues of the sample covariance matrix

• An obvious way is to first compute $$S=\frac{1}{n-1}\sum _{i=1}^n(x_{i\cdot }-\overline{x})(x_{i\cdot }-\overline{x}{)}^{\prime },$$ where $x_{i\cdot }\in R^p$ is the $i$-th observation, and $\overline{x}\in R^p$ is the sample mean vector

• Then compute eigenvalues of $S$

PCA on Sparse Data

• However, this does not fully utilize the sparsity

• In fact, $S$ is generally a dense matrix

• A better scheme is to observe that $$S=\frac{1}{n-1}(X-1_n{\overline{x}}^{\prime }{)}^{\prime }(X-1_n{\overline{x}}^{\prime }),$$ where $X$ is the data matrix, and $1_n$ is a vector of ones

• Ignoring the $(n-1{)}^{-1}$ factor, we only need to compute the eigenvalues of $A^{\prime }A$, where $A=X-1_n{\overline{x}}^{\prime }$

• Also equivalent to the partial SVD of $A$

PCA on Sparse Data

• Then we need to implement the operations $v\to Av$ and $v\to A^{\prime }v$

• As a first step, we need to compute $\overline{x}$. This is relatively easy since we have shown how to extract columns of $X$

• Then to do $v\to Av=(X-1_n{\overline{x}}^{\prime })v$:

  1. $u\leftarrow Xv$

  2. $s\leftarrow {\overline{x}}^{\prime }v$

  3. $w\leftarrow u-s{1}_n$

PCA on Sparse Data

• Similarly, to do $v\to A^{\prime }v=(X-1_n{\overline{x}}^{\prime }{)}^{\prime }v$:

  1. $u\leftarrow X^{\prime }v$

  2. $s\leftarrow 1^{\prime }v$

  3. $w\leftarrow u-s\overline{x}$

Implementation and Benchmark

• See https://statr.me/2019/11/rspectra-center-scale/ for more details

• Especially the center and scale parameters in the RSpectra::svds function

References