DESP (version 0.1-5)

scsSOCP: solve a second-order cone program using SCS

Description

This function solves a second-order cone program using the C package SCS that solves convex cone problems via operator splitting. We use the direct method of the C package SCS that also proposes an indirect method. This method solves an optimization problem of the form : minimize c^T x subject to b - A x = s, s in K, where K is a product of zero, linear and second-order cones, containing the following informations (in this order) : f : number of linear equality constraints, l : length of LP cone, q : array of second-order cone constraints, qsize : length of SOC array (number of second-order cones).

Usage

scsSOCP(model)

Arguments

model

The model representing the convex problem to solve is a list containing the components that follow. The sparse matrix A is stored using the compressed sparse column (CSC) format. [object Object],[object Object],[object Object],[object Object],[object

Value

• scsSOCP returns a list containing the optimal solution, with components:
• x
  The value of the primal variable for the best solution.
• status
  The status of the optimization, returned as an integer : 1 in case of a solved program. A negative status corresponds to an infeasible, unbounded, inaccurate, interrupted, failed or indeterminate problem. Refer to the SCS development pages for more documentation on status codes.

Details

More documentation on SCS can be found at https://github.com/cvxgrp/scs.

References

O'Donoghue, B. and Chu, E. and Parikh, N. and Boyd, S. (2015): Conic Optimization via Operator Splitting and Homogeneous Self-Dual Embedding.

Examples

## This example is the same that the one given for sqR_Lasso() function, 
## with the difference that the optimization problem in R is set up in R,
## before calling scsSOCP() function.
##
## set the design matrix
X <- matrix(c(1,0,2,2,1,0,-1,1,1,2,0,1),4,3,byrow=TRUE)
## set the vector of observations
Y <- c(1,0,2,1)
## set the penalty level
lambda <- 1
## compute the square-root Lasso estimate using SCS
## computation of beta that minimize |Y-X*beta|_2 + lambda |beta|_1
##
## read the sample size and the number of variables
D <- dim(X)
n <- D[1]               # n is the sample size
p <- D[2]               # p is the dimension
## normalize the columns of X by dividing them by their norm
vn <- 1/sqrt(apply(X^2,2,mean))
X_ <- tcrossprod(X,diag(vn))
## minimize |Y-X*beta|_2 + lambda |beta|_1 (square-root Lasso)
## using the variables : x=c(z=abs(beta),beta,v=Y-X*beta,t=norm(v))
## the optimization problem could be rewritten :
##  minimize
##        c^T x
##  subject to
##        v + X beta  = Y    (1)
##        z + beta  >= 0     (2)
##        z - beta  >= 0     (3)
##        |v|_2^ 2 <= t ^ 2  (4) (a second-order cone constraint)
##
## to solve this problem using SCS, we write
##  minimize
##        c^T x
##  subject to
##        b - A x  = s
##        s in K
## where K is a product of zero, linear and second-order cones, 
## containing the following informations (in this order) :
## f = n : number of linear equality constraints (cf. (1))
## l = 2*p : length of LP cone (cf. (2) and (3))
## *q = [1+n] : array of second-order cone constraints (cf. (4))
## qsize = 1 : length of SOC array (number of second-order cones)
##
require(Matrix)
## A matrix
## these lines introduce the constraints (1) : v+X*beta=Y
A <- cBind(Matrix(0,n,p), X_, Diagonal(n), Matrix(0,n,1))
## these lines introduce the constraints (2) and (3) : |beta_{j}| = z_{j}
A <- rBind(A,cBind( -Diagonal(p), -Diagonal(p), Matrix(0, p, n+1)), 
cBind( -Diagonal(p) , Diagonal(p), Matrix(0, p, n+1)))
## these lines introduce the quadratic constraint (4) : |v|_2 <= t
A <- rBind(A, cBind(Matrix(0, 1, 2*p+n),-1), cBind(Matrix(0, n, 2*p), Diagonal(n), Matrix(0, n, 1)))
##
## b vector
b <- c(Y,(1:(2*p+1+n))*0)
##
## c vector : the cost function
c <- c((1:p)*0+lambda, (1:(p+n))*0, 1)
##
## K vector of second-order cone constraints
K <- c(1+n)
##
model <- list()
model$Ax <- A@x # access to the slot containing the non-zero elements of the matrix
model$Ai <- A@i # access to the slot containing the row index of class dgCMatrix,
		# subclass of CsparseMatrix
model$Ap <- A@p # access to the slot containing the column pointer of class dgCMatrix,
		# subclass of CsparseMatrix
model$Am <- 2*n+2*p+1 # A number of rows
model$An <- 2*p+n+1 # A number of columns
model$b <- b
model$c <- c
model$Kf <- n
model$Kl <- 2*p
model$Kq <- K
model$Kqsize <- 1
model$n <- n
model$p <- p
model$settings <- list(normalize=0,verbose=0)
##
r <- try(scsSOCP(model), silent = TRUE)
if (r$status != 1) {
  stop (paste("SCS status :",r$status))
}
if ( inherits (r , "try-error")) {
  stop ("SCS failed somehow !")
}
## get beta, the vector of the coefficients of regression
as.matrix(r$x[p+1:p] * vn)