1 files changed, 52 insertions, 0 deletions
diff --git a/R/run.sparsestep.R b/R/run.sparsestep.R
new file mode 100644
index 0000000..43017c9
--- /dev/null
+++ b/R/run.sparsestep.R
@@ -0,0 +1,52 @@
+# Core SparseStep routine. This could be implemented in a low-level language 
+# in the future, if necessary.
+#
+run.sparsestep <- function(XX, Xy, nvars, lambda, gamma0, gammastep, gammastop,
+			   IMsteps, force.zero, threshold) {
+	# Start solving SparseStep
+	gamma <- gamma0
+	beta <- matrix(0.0, nvars, 1)
+
+	it <- 0
+	while (gamma > gammastop)
+	{
+		for (i in 1:IMsteps) {
+			alpha <- beta
+			omega <- gamma^2/(alpha^2 + gamma^2)^2
+			Omega <- diag(as.vector(omega), nvars, nvars)
+			beta <- solve(XX + lambda * Omega, Xy)
+		}
+		it <- it + 1
+		gamma <- gamma / gammastep
+	}
+
+	# perform IM until convergence
+	epsilon <- 1e-14
+	loss <- get.loss(x, y, gamma, beta, lambda)
+	lbar <- (1.0 + 2.0 * epsilon) * loss
+	while ((lbar - loss)/loss > epsilon)
+	{
+		alpha <- beta
+		omega <- gamma^2/(alpha^2 + gamma^2)^2
+		Omega <- diag(as.vector(omega), nvars, nvars)
+		beta <- solve(XX + lambda * Omega, Xy)
+		lbar <- loss
+		loss <- get.loss(x, y, gamma, beta, lambda)
+	}
+
+	# postprocessing
+	if (force.zero) {
+		beta[which(abs(beta) < threshold)] <- 0
+	}
+	return(beta)
+}
+
+get.loss <- function(x, y, gamma, beta, lambda)
+{
+	Xb <- x %*% beta
+	diff <- y - Xb
+	b2 <- beta^2
+	binv <- 1/(b2 + gamma^2)
+	loss <- t(diff) %*% diff + lambda * t(b2) %*% binv
+	return(loss)
+}