1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
#' @title Coordinate descent algorithm for SparseStep
#'
#'
#' @export
#'
sparsestep.cd <- function(x, y, lambdas=NULL, epsilon=1e-5, intercept=TRUE,
...)
{
nm <- dim(X)
n <- nm[1]
m <- nm[2]
one <- rep(1, n)
if (intercept) {
meanx <- drop(one %*% x)/n
x <- scale(x, meanx, FALSE)
mu <- mean(y)
y <- drop(y - mu)
} else {
meanx <- rep(0, m)
mu <- 0
y <- drop(y)
}
XX <- t(x) %*% x
Xy <- t(x) %*% y
#gammas <- 2^(seq(log(1e6)/log(2), log(1e-8)/log(2)))
num.lambdas <- length(lambdas)
#num.gammas <- length(gammas)
#betas <- array(0, dim=c(num.gammas, num.lambdas, m))
betas <- array(0, dim=c(num.lambdas, m))
for (l in num.lambdas:1) {
lambda <- lambdas[l]
# initialize beta
if (l == num.lambdas) {
beta <- as.vector(matrix(0, 1, m))
} else {
#beta <- betas[num.gammas, l+1, ]
beta <- betas[l+1, ]
}
j <- 1
last.beta <- as.vector(matrix(0, 1, m))
while (TRUE) {
# code
b <- -2 * x[, j] %*% (y - x[, -j] %*% last.beta[-j])
a <- x[, j] %*% x[, j]
if (abs(last.beta[j]) > epsilon) {
beta[j] <- b/a
} else {
beta[j] <- (2*b - lambda*sign(last.beta[j]))/
(2*a)
}
# check convergence
if (sum(abs(beta - last.beta)) < 1e-10) {
break
} else {
last.beta <- beta
}
# continue
j <- j %% m + 1
}
betas[l, ] <- beta
}
return(betas)
}
|
