We present the package sdpt3r, an R implementation of the Matlab
package SDPT3 (Toh et al. 1999). The purpose of the software is to solve
semidefinite quadratic linear programming (SQLP) problems, which
encompasses problems such as D-optimal experimental design, the
nearest correlation matrix problem, and distance weighted
discrimination, as well as problems in graph theory such as finding
the maximum cut or Lovasz number of a graph.
Current optimization packages in R include Rdsdp, Rcsdp, scs,
cccp, and Rmosek. Of these, scs and Rmosek solve a similar
suite of problems. In addition to these solvers, the R packages CXVR
and ROI provide sophisticated modelling interfaces to these solvers.
As a point of difference from the current solvers in R, sdpt3r
allows for log-barrier terms in the objective function, which allows
for problems such as the D-optimal design of experiments to be solved
with minimal modifications . The sdpt3r package also provides helper
functions, which formulate the required input for several well-known
problems, an additional perk not present in the other R packages.
Convex optimization is a well traversed field with far reaching applications. While perhaps unfamiliar to those in the statistical sciences, many problems important to statisticians can be formulated as a convex optimization, perhaps the most well known of which would be the least squares problem. More specifically, many problems in statistics can be formulated as a subset of these convex optimization problems, known as conic linear optimization problems.
One such example would be the nearest correlation matrix problem (Higham 2002), which was first considered when attempting to find correlations between stocks, where incomplete data on daily stock returns are not unusual. Pairwise correlations are only computed when data is available for both pairs of stocks under consideration, resulting in a correlation matrix that contains pairwise correlations, but is not necessarily positive semidefinite - an approximate correlation matrix. The goal is to then find the correlation matrix that is nearest to the approximate correlation matrix in some way.
Other examples of problems that can be formulated in terms of a conic linear optimization problem include D-optimal experimental design (Smith 1918), classification using distance weighted discrimination (Marron et al. 2007), minimum volume ellipsoids (John 2014), and problems in educational testing (Chu and Wright 1995).
Problems in related fields can also be solved, including finding the maximum cut (or maximum k-cut) of a graph, finding the upper bound of the Shannon entropy of a graph, also known as the Lovasz number (Vandenberghe et al. 1998), as well as problems in control theory, Toeplitz matrix approximation, and Chebyshev approximation.
For the purpose of solving these conic linear optimization problems, we introduce the R package sdpt3r, an implementation of the Matlab package SDPT3 by Toh et al. (1999). Of the R packages available to perform conic optimization, sdpt3r is among the most general. Rdsdp (Zhu and Ye 2016) & Rcsdp (Bravo 2016) are capable of solving semidefinite conic optimization problems, while cccp (Pfaff 2015) solves linear and quadratic conic optimization problems. The sdpt3r package allows for all of linear, quadratic, and semidefinite conic optimization to be solved simultaneously (i.e., a problem with any combination of semidefinite, quadratic, or linear cones can be solved). Two comparable packages, scs (O’Donoghue et al. 2017) and Rmosek (Friberg 2012), solve a similar suite of problems. Additionally, the R packages CXVR (Fu et al. 2017) and ROI (Theußl et al. 2017) provide sophisticated modelling interfaces to these solvers.
As a point of difference, scs and Rmosek allow for the exponential and power cones to be included in the constraints, while sdpt3r handles log-barrier terms in the objective function directly. The inclusion of log-barrier terms allows for the D-optimal design of experiments and minimum volume ellipsoid problems to be solved with minimal modifications. In addition, sdpt3r provides helper functions which directly solve a number of well known problems (such as the “max cut” or nearest correlation matrix problems) with minimal input burden on the user. This additional functionality is not found in either scs or Rmosek (although scs can be used with the CVXR package).
This paper is structured as follows. In Section 2 we
discuss in greater detail the mathematical formulation of the linear
conic optimization problem, and introduce three examples to explore the
increasing generality of the problem to be solved. Section
3 discusses the R implementation of sdpt3r, and the main
function by which conic linear optimization problems are solved, sqlp
,
including the required input, and the output generated. The same
examples used in Section 2 will be used to demonstrate how
a standard conic linear optimization problem can be converted to a form
solvable by sqlp
. Section 4 presents the classic form of
several other well known problems that can be solved using sdpt3r, as
well as the helper functions available to convert them to the
appropriate form. Finally Section 5 provides some closing
remarks.
At its simplest, a conic linear optimization problem has the following standard form (Tütüncü et al. 2003):
where
Here,
One of the simplest problems that can be formulated in terms of a conic
linear optimization problem is finding the maximum cut of a graph. Let
Finding the maximum cut is referred to as the Max-Cut Problem, and was one of the first problems found to be NP-complete, and is also one of the 21 algorithms on Karp’s 21 NP-complete problems (Karp 1972). The Max-Cut problem is also known to be APX hard (Papadimitriou and Yannakakis 1991), meaning in addition to there being no polynomial time solution, there is also no polynomial time approximation.
Using the semidefinite programming approximation formulation of
Goemans and Williamson (1995), the Max-Cut problem can be approximated to within
an approximation constant. For a weighted adjacency matrix
where
To see that the Max-Cut problem is a conic linear optimization problem
it needs to be written in the same form as Equation (1). The
objective function is already in a form identical to that of Equation
(1), with minimization occurring over
Now
Allowing for optimization to occur over only one variable at a time is quite restrictive, as only a small number of problems can be formulated in this form. Allowing optimization to occur over multiple variables simultaneously would allow for a broader range of problems to be solved.
The conic linear optimization problem actually covers a much wider class
of problems than those expressible as in Equation (1).
Variables can be separated into those which are constrained to a
semidefinite cone,
Here, svec takes the upper triangular elements of a matrix (including
the diagonal) in a column-wise fashion and vectorizes them. In general
for an
Some important problems in statistics can be formulated to fit this form of the optimization problem.
First addressed by Higham (2002) in dealing with correlations between stock prices, difficulty arises when data is not available for all stocks on each day, which is unfortunately a common occurrence. To help address this situation, correlations are calculated for pairs of stocks only when data is available for both stocks on any given day. The resulting correlation matrix is only approximate in that it is not necessarily positive semidefinite.
This problem was cast by Higham (2002) as
where
Since the matrix
Now, introduce a variable
Here,
To see this as a conic linear optimization problem, notice that
Re-writing the first constraint as
we can easily define the constraint matrices and right hand side of the first constraint as
The second constraint is identical to the constraint from the Max-Cut
problem, where each diagonal element of
yielding
Since
The final step is to combine the individual constraint matrices from
each constraint to form one constraint matrix for each variable, which
is done by defining
While Equation (2) allows for additional variables to be present, it can be made more general still to allow even more problems to be solved. We will refer to this general form as a semidefinite quadratic linear programming (SQLP) problem.
The first generality afforded by an SQLP is the addition of an
unconstrained variable
Recall that for any linear optimization problem, there exists two formulations - the primal formulation and the dual formulation. For the purposes of a semidefinite quadratic linear programming problem, the primal problem will always be defined as a minimization, and the associated dual problem will therefore be a maximization
The primal formulation of the SQLP problem is
For each
As before, the matrices
where
The dual problem associated with the semidefinite quadratic linear programming formulation is
where
Consider the problem of estimating a vector
The variance-covariance matrix of such an estimator is proportional to
There are many different ways in which
Perhaps the most commonly used of these optimality criteria is
D-Optimality, which is equivalent to maximizing the determinant of
Given that the matrix
where
Due to the inequality constraint, this primal formulation cannot be
interpreted as an SQLP of the form of Equation (3). By
defining
Keeping in mind that this is a dual configuration, and thus follows
Equation (4), we proceed with writing the D-Optimal design
problem as an SQLP by first considering the objective function. The
objective function depends only on the determinant of the matrix
variable
The constraint matrices
Then, the constraint matrix is
Finally, define the right hand side of each constraint
which fully specifies the D-Optimal design problem as an SQLP.
In the next section, we will demonstrate using R how these definitions can be translated for use in the main function of sdpt3r so an SQLP problem can be solved.
Each of the problems presented in Section 2 can be solved using the sdpt3r package, an R implementation of the Matlab program SDPT3. The algorithm is an infeasible primal-dual predictor-corrector path-following method, utilizing either an HKM (Helmberg et al. 1996) or NT (Nesterov and Todd 1997) search direction. The interested reader is directed to Tütüncü et al. (2003) for further details surrounding the implementation.
The main function available in sdpt3r is sqlp
, which takes a number
of inputs (or an sqlp_input
object) specifying the problem to be
solved, and executes the optimization, returning both the primal and
dual solution to the problem. This function will be thoroughly discussed
in Section 3.1, and examples will be provided. In addition
to sqlp
, a prospective user will also have access to a number of
helper functions for well known problems that can be solved using
sdpt3r. For example, the function maxcut
takes as input an adjacency
matrix sqlp
. These functions
will be discussed in Sections 3.3, 3.4,
3.4.2, and 4.
For sdpt3r, each optimization variable will be referred to as a
block in the space in which it is restricted. For instance, if we have
an optimization variable
The main function call in sdpt3r is sqlp
, which takes the following
input variables
blk |
A named-list object describing the block structure of the optimization variables. |
At |
A list object containing constraint matrices |
for the primal-dual problem. | |
b |
A vector containing the right hand side of the equality constraints, |
in the primal problem, or equivalently the constant vector in the dual. | |
C |
A list object containing the constant |
function or equivalently the corresponding right hand side of the equality | |
constraints in the dual problem. | |
X0 , y0 , Z0 |
Matrix objects containing an initial iterate for the |
the SQLP problem. If not provided, an initial iterate is computed internally. | |
control |
A list object providing additional parameters for use in sqlp . |
If not provided, default values are used. |
The input variable blk
describes the block structure of the problem.
Letting L
be the total number of semidefinite, quadratic, linear, and
unrestricted blocks in the SQLP problem, define blk
to be a
named-vector object of length
Block type | Name | Value |
---|---|---|
Semidefinite | s | |
Quadratic | q | |
Linear | l | |
Unrestricted | u |
The input variable At
corresponds to the constraint matrices in
Equation (3), and C
the constant matrices in the
objective function. The size of these input variables depends on the
block they are representing, summarized in Table 2 for
each block type.
Block type | ||||
Semidefinite | Quadratic | Linear | Unrestricted | |
At |
||||
C |
Note that in Table 2, At
in the semidefinite block reflects the upper-triangular
input format that has been discussed previously. In a semidefinite
block, the optimization variable
It is important to note that both input variables At
and C
are lists
containing constraint and constant matrices for each optimization
variable. In general, the user need not supply initial iterates X0
,
y0
, and Z0
for a solution to be found using sqlp
. The infeasible
starting point generated internally by sqlp
tends to be sufficient to
find a solution. If the user wishes to provide a starting point however,
the size parameters in Table 3 must be met for each block.
Block type | ||||
Semidefinite | Quadratic | Linear | Unrestricted | |
X0 |
||||
y0 |
||||
Z0 |
The user may choose to depart from the default values of several
parameters which could affect the optimization by specifying alternative
values in the control
list. A complete list of all parameters that can
be altered can be found in Appendix 6.
An important example is the specification of the parbarrier
parameter
in control
, which specifies the presence of a log-barrier in the
objective function. The default case in control
assumes that the
parameters controlparbarrier
to store the values of these parameters
(including zeros). If the
When executed, sqlp
simultaneously solves both the primal and dual
problems, meaning solutions for both problems are returned. The
relevance of each output therefore depends on the problem being solved.
The following object of class sqlp_output
is returned upon completion
pobj |
the value of the primary objective function |
dobj |
the value of the dual objective function |
X |
A list object containing the optimal matrix |
y |
A vector object containing the optimal vector |
Z |
A list object containing the optimal matrix |
The examples in subsequent subsections will demonstrate the output
provided by sqlp
.
Before moving on to more complex problems, consider first some very simple example to illustrate the functionality of the sdpt3r package. First, consider the following simple linear programming problem:
This problem can be solved using sdpt3r in very straightforward
fashion. First, this is a linear programming problem with two variables,
blk = c("l" = 2)
. Next the
objective function can be written as C = matrix(c(1,1),nrow=1)
. The constraints can be summarized in matrix
form as:
so A = matrix(c(1,3,4,-1), nrow=2))
and At = t(A)
. Finally the right
hand side can be written in vector form as b = c(12,10)
.
Pulling these all together, the problem is solved using sqlp
:
blk = c("l" = 2)
C = matrix(c(1,1),nrow=1)
A = matrix(c(1,3,4,-1), nrow=2)
At = t(A)
b = c(12,10)
out = sqlp(blk,list(At),list(C),b)
out
$X
$X[[1]]
2 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 4
[2,] 2
$y
[,1]
[1,] 0.3076923
[2,] 0.2307692
$Z
$Z[[1]]
2 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 6.494441e-10
[2,] 1.234448e-09
$pobj
[1] 6
$dobj
[1] 6
which returns the solution
This problem can be solved using sdpt3r by formulating the input variables in a similar fashion as the linear programming problem:
blk = c("q" = 2)
C = matrix(c(0.5,-1),nrow=1)
A = matrix(c(2,1,-1,1), nrow=2)
At = t(A)
b = c(5,4)
out = sqlp(blk,list(At),list(C),b)
out
$X
$X[[1]]
2 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 3
[2,] 1
$y
[,1]
[1,] 0.5
[2,] -0.5
$Z
$Z[[1]]
2 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 2.186180e-09
[2,] -3.522956e-10
$pobj
[1] 0.5
$dobj
[1] 0.5
which returns the solution
This problem is written almost exactly in the language used by sdpt3, and so can be easily solved by taking:
blk = c("s" = 3)
C = list(matrix(c(1,2,3,2,9,0,3,0,7), nrow=3))
A1 = matrix(c(1,0,1,0,3,7,1,7,5), nrow=3)
A2 = matrix(c(0,2,8,2,6,0,8,0,4), nrow=3)
At = svec(blk,list(A1,A2))
b = c(11,9)
out = sqlp(blk,At,C,b)
out
$X
$X[[1]]
3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 0.08928297 0.1606827 0.2453417
[2,] 0.16068265 0.2891815 0.4415426
[3,] 0.24534167 0.4415426 0.6741785
$y
[,1]
[1,] 0.5172462
[2,] 0.4262486
$Z
$Z[[1]]
3 x 3 Matrix of class "dsyMatrix"
[,1] [,2] [,3]
[1,] 0.4827538 1.147503 -0.9272352
[2,] 1.1475028 4.890770 -3.6207235
[3,] -0.9272352 -3.620723 2.7087744
$pobj
[1] 9.525946
$dobj
[1] 9.525946
which provides the optimal matrix solution svec
is used since
the problem is a semidefinite programming problem, and thus each A
matrix is necessarily symmetric.
Recall that the maximum cut of a graph
where
where we defined
To convert this to a form usable by sqlp
, we begin by noting that we
have one optimization variable, B
of dimension n
for which we would like to
determine the Max-Cut, B
(as in Figure 3.1), blk
is specified as
B <- rbind(c(0, 0, 0, 1, 0, 0, 1, 1, 0, 0),
c(0, 0, 0, 1, 0, 0, 1, 0, 1, 1),
c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0),
c(1, 1, 0, 0, 0, 0, 0, 1, 0, 1),
c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0),
c(1, 1, 0, 0, 1, 0, 0, 1, 1, 1),
c(1, 0, 1, 1, 1, 0, 1, 0, 0, 0),
c(0, 1, 0, 0, 1, 1, 1, 0, 0, 1),
c(0, 1, 0, 1, 1, 0, 1, 0, 1, 0))
n <- max(dim(B))
blk <- c("s" = n)
With the objective function in the form
C
as
one <- matrix(1, nrow = n, ncol = 1)
C <- -(diag(c(B %*% one)) - B) / 4
where, again,
The matrix At
is constructed using the upper triangular portion of the
svec
is made
available in sdpt3r.
A <- list()
for(k in 1:n){
A[[k]] <- Matrix(0,n,n)
A[[k]][k,k] <- 1
}
At <- svec(blk[1],A,1)
Having each of the diagonal elements of b
is a
b <- matrix(1, nrow = n, ncol = 1)
With all the input variables now defined, we can now call sqlp
to
solve the Max-Cut problem
sqlp(blk, At, list(C), b)
The built-in function maxcut
takes as input a (weighted) adjacency
matrix B
and returns the optimal solution directly. If we wish to find
to the maximum cut of the graph in Figure 3.1, given the
adjacency matrix maxcut
as
out <- maxcut(B)
out
$pobj
[1] -14.67622
$X
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
V1 1.000 0.987 -0.136 -0.858 0.480 0.857 -0.879 0.136 -0.857 0.597
V2 0.987 1.000 0.026 -0.763 0.616 0.929 -0.791 -0.026 -0.929 0.459
V3 -0.136 0.026 1.000 0.626 0.804 0.394 0.592 -1.000 -0.394 -0.876
V4 -0.858 -0.763 0.626 1.000 0.039 -0.469 0.999 -0.626 0.470 -0.925
V5 0.480 0.616 0.804 0.039 1.000 0.864 -0.004 -0.804 -0.864 -0.417
V6 0.857 0.929 0.394 -0.469 0.864 1.000 -0.508 -0.394 -1.000 0.098
V7 -0.879 -0.791 0.592 0.999 -0.004 -0.508 1.000 -0.592 0.508 -0.907
V8 0.136 -0.026 -1.000 -0.626 -0.804 -0.394 -0.592 1.000 0.394 0.876
V9 -0.857 -0.929 -0.394 0.470 -0.864 -1.000 0.508 0.394 1.000 -0.098
V10 0.597 0.459 -0.876 -0.925 -0.417 0.098 -0.907 0.876 -0.098 1.000
Note that the value of the primary objective function is negative as we
have defined
As an interesting aside, we can show that the matrix
eigen(out$X[[1]])
$values
[1] 5.59e+00 4.41e+00 2.07e-07 1.08e-07 4.92e-08 3.62e-08 3.22e-08
[8] 1.90e-08 1.66e-08 9.38e-09
The fact that
Recall that the nearest correlation matrix is found as the solution to
In Section 2.1 we wrote this as the following SQLP
for
where
and
To be solved using sqlp
, we first define blk
. There are two
optimization variables in the formulation of the nearest correlation
matrix -
data(Hnearcorr)
X = Hnearcorr
n = max(dim(X))
n2 = n * (n + 1) / 2
blk <- c("s" = n, "q" = n2+1)
Note that C
entry corresponding to the block variable C
as
C1 <- matrix(0, nrow = n, ncol = n)
C2 <- rbind(1, matrix(0, nrow = n2, ncol = 1))
C <- list(C1,C2)
Next comes the constraint matrix for
Aks <- matrix(list(), nrow = 1, ncol = n)
for(k in 1:n){
Aks[[k]] <- matrix(0, nrow = n, ncol = n)
diag(Aks[[k]])[k] <- 1
}
A1s <- svec(blk[1], Aks)[[1]]
A2s <- diag(1, nrow = n2, ncol = n2)
At1 <- cbind(A1s,A2s)
then the constraint matrix for
A1q <- matrix(0, nrow = n, ncol = n2 + 1)
A2q1 <- matrix(0, nrow = n2, ncol = 1)
A2q2 <- diag(1, nrow = n2, ncol = n2)
A2q <- cbind(A2q1, A2q2)
At2 <- rbind(A1q, A2q)
and the right hand side vector b
b <- rbind(matrix(1, n, 1),svec(blk[1], X))
The nearest correlation matrix problem is now solved by
sqlp(blk, list(At1,At2), C, b)
To demonstrate the nearest correlation matrix problem, we will modify an existing correlation matrix by exploring the effect of changing the sign of just one of the pairwise correlations. In the context of stock correlations, we make use of tools available in the R package quantmod (Ryan and Ulrich 2017) to access stock data from five tech firms (Microsoft, Apple, Amazon, Alphabet/Google, and IBM) beginning in 2007.
library("quantmod")
getSymbols(c("MSFT", "AAPL", "AMZN", "GOOGL", "IBM"))
stock.close <- as.xts(merge(MSFT, AAPL, AMZN,
GOOGL, IBM))[, c(4, 10, 16, 22, 28)]
The correlation matrix for these five stocks is
stock.corr <- cor(stock.close)
stock.corr
MSFT.Close AAPL.Close AMZN.Close GOOGL.Close IBM.Close
MSFT.Close 1.0000000 -0.2990463 0.9301085 0.5480033 0.2825698
AAPL.Close -0.2990463 1.0000000 -0.1514348 0.3908624 0.6887127
AMZN.Close 0.9301085 -0.1514348 1.0000000 0.6228299 0.3870390
GOOGL.Close 0.5480033 0.3908624 0.6228299 1.0000000 0.5885146
IBM.Close 0.2825698 0.6887127 0.3870390 0.5885146 1.0000000
Next, consider the effect of having a positive correlation between Microsoft and Apple
stock.corr[1, 2] <- -1 * stock.corr[1, 2]
stock.corr[2, 1] <- stock.corr[1, 2]
stock.corr
MSFT.Close AAPL.Close AMZN.Close GOOGL.Close IBM.Close
MSFT.Close 1.0000000 0.2990463 0.9301085 0.5480033 0.2825698
AAPL.Close 0.2990463 1.0000000 -0.1514348 0.3908624 0.6887127
AMZN.Close 0.9301085 -0.1514348 1.0000000 0.6228299 0.3870390
GOOGL.Close 0.5480033 0.3908624 0.6228299 1.0000000 0.5885146
IBM.Close 0.2825698 0.6887127 0.3870390 0.5885146 1.0000000
Unfortunately, this correlation matrix is not positive semidefinite
eigen(stock.corr)$values
[1] 2.8850790 1.4306393 0.4902211 0.3294150 -0.1353544
Given the approximate correlation matrix stock.corr
, the built-in
function nearcorr
solves the nearest correlation matrix problem using
sqlp
out <- nearcorr(stock.corr)
Since this is a minimization problem, and thus a primal formulation of
the SQLP, the output X
from sqlp
will provide the optimal solution
to the problem - that is, X
will be the nearest correlation matrix to
stock.corr
.
out$X
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 0.25388359 0.86150833 0.5600734 0.3126420
[2,] 0.2538836 1.00000000 -0.09611382 0.3808981 0.6643566
[3,] 0.8615083 -0.09611382 1.00000000 0.6115212 0.3480430
[4,] 0.5600734 0.38089811 0.61152116 1.0000000 0.5935021
[5,] 0.3126420 0.66435657 0.34804303 0.5935021 1.0000000
The matrix above is symmetric with unit diagonal and all entries in
eigen(out$X)
$values
[1] 2.846016e+00 1.384062e+00 4.570408e-01 3.128807e-01 9.680507e-11
we can see that X
is indeed a correlation matrix.
Recall from Section 2.2 that the D-Optimal experimental design problem was stated as the following dual SQLP
which we wrote as
where
Here,
To convert this to a form usable by sdpt3r, we first declare the three
blocks in blk
. For a matrix The first block is semidefinite containing
the matrix
data(DoptDesign)
V = DoptDesign
n = nrow(V)
p = ncol(V)
blk = c("s" = n, "l" = p, "u" = 1)
Next, by noting the variable b
as a vector of zeros
b <- matrix(0, nrow = p, ncol = 1)
Next, looking at the right-hand side of the constraints, we define the
matrices C
C1 <- matrix(0, nrow = n, ncol = n)
C2 <- matrix(0, nrow = p, ncol = 1)
C3 <- 1
C = list(C1,C2,C3)
Finally, we construct At
for each variable
A <- matrix(list(), nrow = p, ncol = 1)
for(k in 1:p){
A[[k]] <- -V[,k] %*% t(V[,k])
}
At1 <- svec(blk[1], A)[[1]]
At2 <- diag(-1, nrow = p, ncol = p)
At3 <- matrix(1, nrow = 1, ncol = p)
At = list(At1,At2,At3)
The final hurdle necessary to address in this problem is the existence
of the log-barrier. Recall that it is assumed that control
. In this
case, we can see that is not true, as we have a log term containing
sqlp
, we define the controlparbarrier
variable as
control <- list(parbarrier = matrix(list(),3,1))
control$parbarrier[[1]] <- 1
control$parbarrier[[2]] <- 0
control$parbarrier[[3]] <- 0
The D-Optimal experimental design problem can now be solved using sqlp
sqlp(blk, At, C, b, control)
To demonstrate the output generated from a D-optimal experimental design
problem, we consider a simple sqlp
package). To solve the problem usingsqlp
, we use the
function doptimal
, which takes as input an y
, corresponding to
data("DoptDesign")
out <- doptimal(DoptDesign)
out$y
[,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.000
[6,] 0.000
[7,] 0.154
[8,] 0.000
[9,] 0.000
[10,] 0.000
[11,] 0.000
[12,] 0.000
[13,] 0.319
[14,] 0.000
[15,] 0.000
[16,] 0.240
[17,] 0.000
[18,] 0.000
[19,] 0.000
[20,] 0.000
[21,] 0.000
[22,] 0.000
[23,] 0.287
[24,] 0.000
[25,] 0.000
The information matrix y
above.
The sdpt3r package considerably broadens the set of optimization problems that can be solved in R. In addition to those problems presented in detail in Section 3, there are a large number of well known problems that can also be formulated as an SQLP.
Each problem presented will be described briefly, with appropriate
references for the interested reader, and presented mathematically in
its classical form, not as an SQLP as in Equation (3) or
(4). Accompanying each problem will be an R helper function,
which will solve the corresponding problem using sqlp
. Each helper
function in sdpt3r (including those for the max-cut, D-optimal
experimental design, and nearest correlation matrix) is an R
implementation of the helper functions that are available to the user in
the Matlab SDPT3 package (Toh et al. 1999).
The problem of finding the ellipsoid of minimum volume containing a set
of points
The function minelips
takes as input an sqlp
.
data(Vminelips)
out <- minelips(Vminelips)
Given two sets of points in a matrix
Of course, it may be impossible to have a perfect separation of points
using a linear hyperplane, so an error term
Distance weighted discrimination (Marron et al. 2007) solves the following optimization problem to find the optimal hyperplane.
where
The function dwd
takes as input two X1
and
X2
containing the points to be separated, as well as a penalty term
C
sqlp
.
data(Andwd)
data(Apdwd)
C <- 0.5
out <- dwd(Apdwd,Andwd,penalty)
Similar to the Max-Cut problem, the Max-kCut problem asks, given a graph
Here,
The function maxkcut
takes as input an adjacency matrix B
and an
integer k
, and returns the optimal solution using sqlp
.
data(Bmaxkcut)
k = 2
out <- maxkcut(Bmaxkcut,k)
The graph partitioning problem can be formulated as the following primal optimization problem
Here,
The function gpp
, takes as input a weighted adjacency matrix B
and a
real number alpha
and returns the optimal solution using sqlp
.
data(Bgpp)
alpha <- nrow(Bgpp)
out <- gpp(Bgpp, alpha)
The Lovasz Number of a graph
The function lovasz
takes as input an adjacency matrix sqlp
.
data(Glovasz)
out <- lovasz(Glovasz)
Given a symmetric matrix
is a general Toeplitz matrix. The problem is formulated as the following optimization problem
where
The function toep
takes as input a symmetric matrix F
for which we
would like to find the nearest Toeplitz matrix, and returns the optimal
solution using sqlp
.
data(Ftoep)
out <- toep(Ftoep)
The educational testing problem arises in measuring the reliability of a
student’s total score in an examination consisting of a number of
sub-tests (Fletcher 1981). In terms of formulation as an
optimization problem, the problem is to determine how much can be
subtracted from the diagonal of a given symmetric positive definite
matrix
The Educational Testing Problem (ETP) is formulated as the following dual problem
where
The corresponding primal problem is
The function etp
takes as input an A
, and returns the optimal solution using sqlp
.
data(Betp)
out <- etp(Betp)
For a
where
The function logcheby
takes as input a matrix B
and vector f
, and
returns the optimal solution to the Logarithmic Chebyshev Approximation
problem using sqlp
.
data(Blogcheby)
data(flogcheby)
out <- logcheby(Blogcheby, flogcheby)
We consider three distinct linear matrix inequality problems, all
written in the form of a dual optimization problem. The first linear
matrix inequality problem we will consider is defined by the following
optimization equation for some
The function lmi1
takes as input a matrix sqlp
.
B <- matrix(c(-1,5,1,0,-2,1,0,0,-1), nrow=3)
out <- lmi1(B)
The second linear matrix inequality problem is
Here, the matrices
The function lmi2
takes the matrices A1
, A2
, and B
as input, and
returns the optimal solution using sqlp
.
A1 <- matrix(c(-1,0,1,0,-2,1,0,0,-1),3,3)
A2 <- A1 + 0.1*t(A1)
B <- matrix(c(1,3,5,2,4,6),3,2)
out <- lmi2(A1,A2,B)
The final linear matrix inequality problem originates from a problem in
control theory (Boyd et al. 1994) and requires three matrices be known in
advance,
The function lmi3
takes as input the matrices A
, B
, and G
, and
returns the optimal solution using sqlp
.
A <- matrix(c(-1,0,1,0,-2,1,0,0,-1),3,3)
B <- matrix(c(1,2,3,4,5,6), 2, 3)
G <- matrix(1,3,3)
out <- lmi3(A,B,G)
In Section 2, we introduced the problem of conic linear optimization. Using the Max-Cut, Nearest Correlation Matrix, and D-Optimal Experimental Design problems as examples, we demonstrated the increasing generality of the problem, culminating in a general form of the conic linear optimization problem, known as the semidefinite quadratic linear program, in Section 2.2.
In Section 3, we introduced the R package sdpt3r, and
the main function call available in the package, sqlp
. The specifics
of the necessary input variables, the optional input variables, and the
output variables provided by sqlp
were presented. Using the examples
from Section 2, we showed how a problem written as a
semidefinite quadratic linear program could be solved in R using
sdpt3r.
Finally, in Section 4, we presented a number of additional
problems that can be solved using the sdpt3r package, and presented
the helper functions available so these problems could be easily solved
using sqlp
.
The sdpt3r package broadens the range of problems that can be solved
using R. Here, we discussed a number of problems that can be solved
using sdpt3r, including problems in the statistical sciences, graph
theory, classification, control theory, and general matrix theory. The
sqlp
function in sdpt3r is in fact even more general, and users may
apply it to any other conic linear optimization problem that can be
written in the form of Equation (3) or (4) by
specifying the input variables blk
, At
, C
, and b
for their
particular problem.
vers |
specifies the search direction |
0, HKM if semidefinite blocks present, NT otherwise (default) | |
1, HKM direction | |
2, NT direction | |
predcorr |
TRUE, use Mehrotra prediction-correction (default) |
FALSE, otherwise | |
gam |
step-length (default 0) |
expon |
exponent used to decrease sigma (default 1) |
gaptol |
tolerance for duality gap as a fraction of the objective function (default |
inftol |
tolerance for stopping due to infeasibility (default 1e-8) |
steptol |
tolerance for stopping due to small steps (default 1e-6) |
maxit |
maximum number of iterations (default 100) |
stoplevel |
0, continue until successful completion, maximum iteration, or numerical failure |
1, automatically detect termination, restart if small steps is cause (default) | |
2, automatically detect termination | |
scale_data |
TRUE, scale data prior to solving |
FALSE, otherwise (default) | |
rmdepconstr |
TRUE, remove nearly dependent constraints |
FALSE, otherwise (default) | |
parbarrier |
declare the existence of a log barrier term |
default value is 0 (i.e., no log barrier) |
sdpt3r, Rdsdp, Rcsdp, cccp, scs, Rmosek, quantmod
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Rahman, "sdpt3r: Semidefinite Quadratic Linear Programming in R", The R Journal, 2018
BibTeX citation
@article{RJ-2018-063, author = {Rahman, Adam}, title = {sdpt3r: Semidefinite Quadratic Linear Programming in R}, journal = {The R Journal}, year = {2018}, note = {https://rjournal.github.io/}, volume = {10}, issue = {2}, issn = {2073-4859}, pages = {371-394} }