--- title: "Use Cases for LazyMatrix: Statistical Algorithms" output: rmarkdown::html_vignette: number_sections: true bibliography: references.bib vignette: > %\VignetteIndexEntry{Use Cases for LazyMatrix: Statistical Algorithms} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Linear Regression with LSQR ## Introduction Several algorithms for iterative least squares algorithms have been proposed when working with sparse matrices. We have chosen to work with the method proposed by @paige1982lsqr that is called LSQR. This iterative algorithm allows us to use custom methods for matrix operations which allows for our lazy computations. We will see that the operations `%*%`, `crossprod` and `norm` are enough for getting stable estimates of the regression coefficients. ## Notation The authors proposes an algorithm as a set of initalisation parameters and the design matrix $A$. As mentioned by the authors, the matrix $A$ is used only to compute products of the form $Av$ and $A^Tu$. This object will be our `LazyMatrix` objects and hence, all computations will be done lazily at every step of the iteration. For clarity, matrices will be denoted by capital letters $A, B,\cdots$, vectors by $v,w,\cdots$ and scalars by greek letters $\alpha, \beta,\cdots$. For vectors, the norm is the Euclidean norm and for matrices we assume the Frobenius norm $||A||_F=\sqrt{\sum |a_{ij}|^2}$. ## Mathematical Background We will make use of the [Lanczos process](https://en.wikipedia.org/wiki/Lanczos_algorithm), for which we will not go into further detail. The main idea is that with the help of a set of initialization parameters, we can with the help of the design matrix $A$ update the regression coefficients for every step of the algorithm. The algorithm is explained as 1. Initialization: \begin{align*} \beta_1 &= \sqrt{\sum b^2}, \\ \mathbf{u}_1 &= \frac{\mathbf{b}}{\beta_1}, \\ \mathbf{A}^\top \mathbf{u}_1 &= \mathbf{A}^\top \mathbf{u}_1, \\ \alpha_1 &= \sqrt{\sum (\mathbf{A}^\top \mathbf{u}_1)^2}, \\ \mathbf{v}_1 &= \frac{\mathbf{A}^\top \mathbf{u}_1}{\alpha_1}, \\ \mathbf{w}_1 &= \mathbf{v}_1, \\ \mathbf{x}_0 &= \mathbf{0}, \\ \bar{\phi_1} &= \beta_1, \\ \bar{\rho_1} &= \alpha_1. \end{align*} 2. For $i = 1, 2, \dots$, repeat steps 3--6: 3. Continue the bidiagonalization: \begin{align*} \mathbf{\beta}_u &= \mathbf{A} \mathbf{v}_1 - \alpha_1 \mathbf{u}_1, \\ \beta_2 &= \sqrt{\sum \mathbf{\beta}_u^2}, \\ \mathbf{u}_2 &= \frac{\mathbf{\beta}_u}{\beta_2}, \\ \mathbf{\alpha}_v &= \mathbf{A}^\top \mathbf{u}_2 - \beta_2 \mathbf{v}_1, \\ \alpha_2 &= \sqrt{\sum \mathbf{\alpha}_v^2}, \\ \mathbf{v}_2 &= \frac{\mathbf{\alpha}_v}{\alpha_2}. \end{align*} 4. Construct and apply the next orthogonal transformation: \begin{align*} \rho_1 &= \sqrt{(\bar{\rho_1})^2 + \beta_2^2}, \\ c_1 &= \frac{\bar{\rho_1}}{\rho_1}, \\ s_1 &= \frac{\beta_2}{\rho_1}, \\ \theta_2 &= s_1 \alpha_2, \\ \bar{\rho_2} &= -c_1 \alpha_2, \\ \phi_1 &= c_1\bar{\phi_1}, \\ \bar{\phi_2} &= s_1 \bar{\phi_1}. \end{align*} 5. Update $\mathbf{x}$ and $\mathbf{w}$: \begin{align*} \mathbf{x}_1 &= \mathbf{x}_0 + \frac{\phi_1}{\rho_1} \mathbf{w}_1, \\ \mathbf{w}_2 &= \mathbf{v}_2 - \frac{\theta_2}{\rho_1} \mathbf{w}_1. \end{align*} 6. Reset the loop variables: \begin{align*} \beta_1 &= \beta_2, \\ \mathbf{u}_1 &= \mathbf{u}_2, \\ \alpha_1 &= \alpha_2, \\ \mathbf{v}_1 &= \mathbf{v}_2, \\ \bar{\rho_1} &= \bar{\rho_2}, \\ \bar{\phi_1} &= \bar{\phi_2}, \\ \mathbf{x}_0 &= \mathbf{x}_1, \\ \mathbf{w}_1 &= \mathbf{w}_2. \end{align*} 7. Check for convergence. At this point, the authors do not propose a criterion for when we should stop iterating. However, we have chosen to use the criterion that [LSQR in scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsqr.html) [@virtanen2020scipy; @scipy_lsqr] uses. First, compute the residual between the response vector $y$ and the estimate as $$ r = y-Ax_0 $$ and then we stop iterating when $$ \|r\| \leq \text{tol}\cdot \|A\| \cdot \|x_0\| + \text{tol} \cdot \|b\|. $$ `tol` is a parameter for tolerance level of the difference that we accept. In our algorithm, it is by default $1e-6$. ## Implementation Assuming we want to fit regression coefficient vector $\beta$ on the scaled matrix $X$ with response vector $y$. We can compute these estimates using lazy computation. Assuming the original sparse matrix `A`, create an instance of `LazyMatrix` and run the method `lazymatrix::lsqr`. ```{r setup} library(lazymatrix) # 1. Define sparseMatrix base::set.seed(123) n_row <- 50 n_col <- 10 i <- c( 1:50, base::sample(1:50, 20, replace = TRUE), base::sample(1:50, 15, replace = TRUE) ) j <- c( base::rep_len(1:10, 50), base::sample(1:10, 20, replace = TRUE), base::sample(1:10, 15, replace = TRUE) ) pairs <- base::unique(data.frame(i = i, j = j)) i <- pairs$i j <- pairs$j x <- stats::rnorm(length(i)) A <- Matrix::sparseMatrix(i = i, j = j, x = x, dims = c(n_row, n_col)) # 2. Define response vector y <- stats::rnorm(nrow(A)) # 3. Define LazyMatrix X <- LazyMatrix(A, "sd", "mean") # 4. Perform lsqr lazy_beta <- lazymatrix::lsqr(X, y) ``` Now,the non-lazy approach would be to scale `A` using `base::scale()`and then use something like `stats::lm.fit()` to get the estimates. This would imply creating a copy of `A` with non-zero elements taking up unneccessary memory and being much slower generally. ```{r non-lazy-approach} scaled_a <- base::scale(A) non_lazy_beta <- stats::lm.fit(scaled_a, y)$coefficients ``` The estimates are equivalent as is shown below. ```{r result} isTRUE(all.equal(as.vector(lazy_beta), as.vector(non_lazy_beta))) ``` # Gradient Descent ## Introduction In statistical modeling, the classical gradient descent algorithm remains a fundamental part in any kind of modern machine learning approach. As an example of the power of `lazymatrix`, we propose The approach we propose the simplest form of gradient descent on a scaled matrix. As many of our implementations show, it is sufficient to use a small set of matrix operations in order to perform this algorithm. ## Mathematical Background The algorithm follows the approach as explained by @ng2025cs229. Assume the cost function $J(\theta)$ defined as $$ J(\theta)=\frac{1}{2}\sum_{i=1}^n(h_\theta(x^{(i)})-y^{(i)})^2 $$ where $h_\theta(x)=\theta^Tx$, $x$ being the elements of the design matrix and $\theta$ the parameters. In gradient descent, we usually talk about two types of parameters, weights and biases. In the notation above, they are both included within the parameter vector $\theta$. This is done by letting $x_0=1$ so that the first term in $\sum_{i=0}^n\theta_ix_i$ includes the additive bias term. $y$ is the response vector, so we recognize $J$ to be the cost function representing the residual sum of squares. Hence, the goal is to choose $\theta$ to minimize $J$. We start with an initial guess of the parameter and repeatedly update it to make the cost function smaller. This is done with the gradient descent algorithm defined as $$ \theta_j:=\theta_j - \alpha\frac{\partial}{\partial\theta_j}J(\theta) $$ where the partial derivative is the derivative of the cost function with regards to the parameter. $\alpha$ is the learning rate that usually takes a small value to control how fast the algorithm converges. Computing this derivative is straightforward and we get the expression $$ \frac{\partial}{\partial\theta_j}J(\theta)=(h_\theta(x)-y)x_j, $$ providing us with the update rule $$ \theta_j:=\theta_j+\alpha(y^{(i)}-h_\theta(x^{(i)}))x_j^{(i)}. $$ ## Implementation The following gradient descent algorithm works as helper function for which we will compare the lazy-approach to the non-lazy. ```{r gradient descent algorithm} gradient_descent_helper <- function( x, y, w_init, b_init, learning_rate, n_epochs ) { w <- w_init b <- Matrix::Matrix(b_init, nrow = nrow(x)) n <- nrow(x) for (epoch in 1:n_epochs) { y_pred <- x %*% w + b error <- y_pred - y w <- w - (learning_rate * crossprod(x, error)) / n b <- b - (learning_rate * sum(error)) / n } list(w = w, b = b) } ``` Compared to the definition above, in this algorithm we include the bias term as `b`. This results in a list of the two parameter vectors, weights and biases. `n_epochs` is the number of time we choose to iterate. As the result of `LazyMatrix %*% vector` results in a `Matrix`-object, the addition of the bias term can be handled completely within the `Matrix` framework. The key of `lazymatrix` is to handle the heavy matrix operations lazily which is done in this case using `%*%` and `crossprod`. Now, we use the sparse matrix `A` from section #1 and it's `LazyMatrix` `X` and set initial values randomly for the weights and biases and the response vector `y_true`. ```{r initialization} set.seed(4567) p <- ncol(X) w_init <- stats::rnorm(p) b_init <- stats::rnorm(1) y_true <- stats::rnorm(nrow(A)) learning_rate <- 0.01 n_epochs <- 10 ``` Note that the function returns the parameters from the iterative algorithm. ```{r} pars_lazy <- gradient_descent_helper( x = X, y = y_true, w_init = w_init, b_init = b_init, learning_rate = learning_rate, n_epochs = n_epochs ) print("First 5 parameter weights and the bias term.") print(pars_lazy$w[1:5]) print(pars_lazy$b[1]) ``` These parameters can now be used for predictions and are comparible with the non-lazy way of using `scale(A)`and then running gradient descent. ## Discussion Nonetheless there exists much more complex versions of the gradient descent algorithm, this example rather proves that with a few matrix operations being performed lazily, we can avoid storing heavy copies of large matrices and perform operations solely when they are needed. The idea is to implement methods to `lazymatrix` which allows for more complex algorithms such as Stochastic gradient descent. # Principal Component Analysis ## Introduction A fundamental statistical algorithm for dimension reduction is principal component analysis (PCA). In base R, the function used to perform PCA is `stats::prcomp()`, which does so by breaking up the matrix `A` into it's singular value decomposition (SVD) using `base::svd()`. This implies that we require a lazy implementation of svd in order for principal component analysis to work. Moreover, many methods have been proposed to perform svd on sparse matrices and in `lazymatrix` we have chosen to follow the iterative algorithm proposed by @baglama2005augmented, which the authors have also implemented in the R package `irlba`. As this method uses iterative methods, we can easily include `lazymatrix`'s methods for `%*%` and `crossprod` in order to get a working algorithm for svd. However, as the authors themselves say in the paper, the `irlba` function will not compute all singular values and proposes the user to use `base::svd` if one wishes to compute all singular values. This would however involve materializing `A` breaking the lazy computation in the process. Therefore, we work solely with the algorithm proposed by `irlba`. ## Mathematical Background ### Partial Singular Value Decomposition with `irlba` Assume the matrix $\tilde{A}$ being the scaled version of `A` being $n\times p$. We wish to perform a partial singular value decomposition, where the first step is using a Lanczos process which can be described as 1. Choose a random unit vector $p_1\in R^p$ 2. Compute $\tilde{A}p_1=z\in R^n$, transforming $p_1$ from $R^p$ to $R^n$ 3. Normalize the vector $q_1=\frac{z}{\lVert z \rVert}$ and let $\alpha_1=\lVert z \rVert$. This builds up the matrix $B$ as it's first diagonal entry. 4. Compute the residual vector by first returning to $R^p$ with $\tilde{A}^Tq_1=z\in R^p$ and let $$ r=z-\alpha_1p_1 \in R^p. $$ 5. Compute the length of the residual vector $\beta_1=\lVert r \rVert$ and let $$ p_2 = \frac{r}{\beta_1}. $$ 6. Continue to build up $B$ with $\beta_1$ as $$ B=\begin{pmatrix} \alpha_1 & \beta_1 & \end{pmatrix} $$ At the next iteration, $p_2$ is used as input vector and we start over until we reach convergence after $m$ steps, resulting in the bidiagonal form of $B\in R^{m\times m}$. At every point when matrix multiplication is being computed with $\tilde{A}$ or $\tilde{A}^T$, lazy computation through `lazymatrix` is being implemented. Once we have performed every iteration in the Lanczos process, we can stack the right basis vectors $p_1, p_2, \ldots, p_m$ into the matrix $P_m\in R^{p\times m}$ and the left basis vectors $q_1, q_2,...,q_m$ into $Q_m\in R^{n\times m}$ resulting in the equations $$ \begin{align*} \tilde{A}P_m &= Q_mB_m\\ \tilde{A}^TQ_m &= P_mB_m^T+r_me_m^T \end{align*} $$ where we know that $B_m$ is bidiagonal and since $m \ll \min(n,p)$, this is a smaller matrix and the decomposition can be written as $$ \tilde{A}P_m=Q_mB_m $$ where we know that $B_m$ is bidiagonal and since $m<<\min(n,p)$, this is a smaller matrix and the decomposition can be written as $$ B_m=U_B\Sigma_B V_B^T $$ where $$ \begin{align*} \Sigma_B&=\text{diag}(\sigma_1^B, \sigma_2^B,\ldots,\sigma_m^B)\newline U_B&\in R^{m\times m} \\ V_B&\in R^{m\times m}. \end{align*} $$ This decompoisition can be used to get approximate solutions to the svd of $\tilde{A}$ with the relationships $$ \begin{align*} \tilde{\sigma}_j&=\sigma_j^B\\ \tilde{u}_j &= Q_mu_j^B \in R^n\\ \tilde{v}_j &= P_mv_j^B \in R^p.\\ \end{align*} $$ These equations can be validated easily by looking first at the left singular vectors $$ \tilde{A}\tilde{v}_j=\tilde{A}P_mv_j^B=Q_mB_mv_j^B=\sigma_j^BQ_mu_j^B=\tilde{\sigma}_j\tilde{u}_j $$ which implies that $(\tilde{\sigma}_j,\tilde{u}_j, \tilde{v}_j)$ satisfies the left singular vector equation. Instead, the second equation $\tilde{u}_j = Q_mu_j^B$, has a residual error if we look at $$ \tilde{A}^TQ_m=P_mB_m^T+r_me_m^T $$ where it can be shown that $$ \tilde{A}^T\tilde{u_j}=\tilde{\sigma}_j\tilde{v_j}+(e_m^Tu_j^B)r_m. $$ This error term is proportional to the last computation of $u_j^B$ and residual $r_m$. Now, we are not neccessarily interested in all $m$ iterations, but it is sufficient to choose $k