Sparse linear regression

sparse-linreg.R

library(Matrix)

rm(list=ls())

set.seed(123)


## parameters

n <- 1e6
p1 <- 10
p2 <- 1e4
s <- 1e3


## generate X

nn <- n*p2/s

X1 <- Matrix(rnorm(n*p1), n, p1, sparse = TRUE)

system.time( ix <- sample(n, nn, replace = TRUE) )
jx <- sample(p2, nn, replace = TRUE)
system.time( X2 <- sparseMatrix(i = ix, j = jx, dims = c(n,p2), x = rnorm(nn)) )

system.time( X <- cBind(X1, X2) )


## generate y

beta0 <- rep(1, p1+p2)
beta0[1] <- 2
beta0[p1+1] <- 2.5

eps <- rnorm(n, sd = 0.1)
system.time( y <- X %*% beta0 + eps )


## solve lin.regr. y ~ X

system.time( XX <- crossprod(X) )
system.time( Xy <- crossprod(X,y) )
system.time( beta <- solve(XX, Xy) )


## verify beta (should be approx beta0)

head(beta0,20)
head(beta,20)
plot(c(beta[2:10],beta[12:20]))