## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----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)

## ----non-lazy-approach--------------------------------------------------------
scaled_a <- base::scale(A)
non_lazy_beta <- stats::lm.fit(scaled_a, y)$coefficients

## ----result-------------------------------------------------------------------
isTRUE(all.equal(as.vector(lazy_beta), as.vector(non_lazy_beta)))

## ----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)
}

## ----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

## -----------------------------------------------------------------------------
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])

## ----svd, warning=FALSE-------------------------------------------------------
#3. SVD
svd_lazy <- svd(X)

# 4. Print singular values
print("The singular values are: ")
svd_lazy$d

print("The first 6 rows of the left singular vectors are: ")
head(svd_lazy$u)

print("The frist 6 rows of the right singular vectors are: ")
head(svd_lazy$v)

## ----prcomp, warning=FALSE----------------------------------------------------
pca_lazy <- prcomp(X)
print("The first 6 rows of the principal components are: ")
head(pca_lazy$rotation)