## Sparse matrix should be operated on with suitable data structure

From wikipedia:

When storing and manipulating sparse matrices on a computer, it is beneficial and often necessary to use specialized algorithms and data structures that take advantage of the sparse structure of the matrix. Operations using standard dense-matrix structures and algorithms are relatively slow and consume large amounts of memory when applied to large sparse matrices. Sparse data is by nature easily compressed, and this compression almost always results in significantly less computer data storage usage. Indeed, some very large sparse matrices are infeasible to manipulate using standard dense algorithms.

R has the Matrix package with special data structures for dealing with this. Here is an demonstration in R:

Performance test of sparse matrix in R
n = 9000
k = 100
u = matrix(rnorm(n*k), k)
v = matrix(rnorm(n*k), n)
w = rnorm(n)
W1 = diag(w)
W2 = Diagonal(x=w)
system.time(u %*% W1 %*% v)
system.time(u %*% W2 %*% v)

#////////////////////////////////////////////////////
system.time(u %*% W1 %*% v)
user  system elapsed
1.956   0.000   1.956
#////////////////////////////////////////////////////
system.time(u %*% W2 %*% v)
user  system elapsed
0.081   0.004   0.086



This brings about an staggering improvement of over 20 times! Let's revisit one previous example: logistic regression implemented with batch gradient descent (in that post we didn't calculate the standard error of coefficients, if we do, glm will beat us by a large margin!), but this time take advantage of sparse matrix:

Logistic regression implementation improved
require(Matrix)
logreg = function(y, x) {
x = as.matrix(x)
x = apply(x, 2, scale)
x = cbind(1, x)
m = nrow(x)
n = ncol(x)
alpha = 2/m

# b = matrix(rnorm(n))
# b = matrix(summary(lm(y~x))$coef[, 1]) b = matrix(rep(0, n)) v = exp(-x %*% b) h = 1 / (1 + v) J = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h)) derivJ = t(x) %*% (h-y) niter = 0 while(1) { niter = niter + 1 newb = b - alpha * derivJ v = exp(-x %*% newb) h = 1 / (1 + v) newJ = -(t(y) %*% log(h) + t(0-y) %*% log(1 -h)) while((newJ - J) >= 0) { print("inner while...") # step adjust alpha = alpha / 1.15 newb = b - alpha * derivJ v = exp(-x %*% newb) h = 1 / (1 + v) newJ = -(t(y) %*% log(h) + t(1-y) %*% log(1 -h)) } if(max(abs(b - newb)) < 0.001) { break } b = newb J = newJ derivJ = t(x) %*% (h-y) } b w = h^2 * v # # hessian matrix of cost function hess = t(x) %*% Diagonal(x = as.vector(w)) %*% x seMat = sqrt(diag(solve(hess))) zscore = b / seMat cbind(b, zscore) } nr = 5000 nc = 5 # set.seed(17) x = matrix(rnorm(nr*nc, 0, 999), nr) x = apply(x, 2, scale) # y = matrix(sample(0:1, nr, repl=T), nr) h = 1/(1 + exp(-x %*% rnorm(nc))) y = round(h) y[1:round(nr/2)] = sample(0:1, round(nr/2), repl=T) ntests = 300 testglm = function() { for(i in 1:ntests) { res = summary(glm(y~x, family=binomial))$coef
}
print(res)
}

testlogreg = function() {
for(i in 1:ntests) {
res = logreg(y, x)
}
print(res)
}

print(system.time(testlogreg()))
print(system.time(testglm()))

# benchmark results
#////////////////////////////////////////////////////
print(system.time(testlogreg()))
[,1]      [,2]
[1,] -0.00639276 -0.205752
[2,] -0.23057674 -7.327096
[3,] -0.05069555 -1.634037
[4,] -0.00793098 -0.254736
[5,]  0.47984767 14.795614
[6,]  0.84068549 23.592259
user  system elapsed
8.671   0.000   8.673
#////////////////////////////////////////////////////
print(system.time(testglm()))
Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -0.00641490  0.0310827 -0.206382  8.36493e-01
x1          -0.23164275  0.0314839 -7.357510  1.87372e-13
x2          -0.05082648  0.0310369 -1.637614  1.01502e-01
x3          -0.00789098  0.0311468 -0.253348  7.99999e-01
x4           0.48154867  0.0324508 14.839341  8.15527e-50
x5           0.84393705  0.0356694 23.659992 9.31454e-124
user  system elapsed
7.571   0.000   7.572


Not so bad, sometimes our algorithm can even beat glm.