%\VignetteIndexEntry{sdpt3r}
% JSS Article Template
\documentclass[nojss]{jss}

%%% Additional Packages %%%

% Algorithms
% http://tex.stackexchange.com/questions/229355/algorithm-algorithmic-algorithmicx-algorithm2e-algpseudocode-confused
\usepackage{algorithm}% http://ctan.org/pkg/algorithms
\usepackage{algpseudocode}% http://ctan.org/pkg/algorithmicx


% Font stuff
\usepackage[T1]{fontenc}% for correct hyphenation and T1 encoding
%\usepackage{lmodern}% latin modern font
\usepackage[american]{babel}% for American English

% Math fonts/symbols
\usepackage{amsmath}% sophisticated mathematical formulas with amstex (includes \text{})
\usepackage{mathtools}% fix amsmath deficiencies
\usepackage{amssymb}% sophisticated mathematical symbols with amstex
\usepackage{amstext}% ams mathematical fonts
\usepackage{amsfonts}% ams mathematical fonts

% Figures
\usepackage{graphicx}% for including figures
%\usepackage{epstopdf}
\ifpdf
  \DeclareGraphicsExtensions{.eps,.pdf,.png,.jpg}
\else
  \DeclareGraphicsExtensions{.eps}
\fi

%%% Commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% vetors
\newcommand{\ve}[1]{\mathbf{#1}}           % for vetors
\newcommand{\sv}[1]{\boldsymbol{#1}}   % for greek letters
\newcommand{\m}[1]{\mathbf{#1}}               % for matrices
\newcommand{\sm}[1]{\boldsymbol{#1}}   % for greek letters
\newcommand{\tr}[1]{{#1}^{\mkern-1.5mu\mathsf{T}}}              % for transpose
\newcommand{\norm}[1]{||{#1}||}              % for transpose
\newcommand*{\mve}{\operatorname{ve}}
\newcommand*{\trace}{\operatorname{trace}}
\newcommand*{\rank}{\operatorname{rank}}
\newcommand*{\diag}{\operatorname{diag}}
\newcommand*{\vspan}{\operatorname{span}}
\newcommand*{\rowsp}{\operatorname{rowsp}}
\newcommand*{\colsp}{\operatorname{colsp}}
\newcommand*{\svd}{\operatorname{svd}}
\newcommand*{\edm}{\operatorname{edm}}  % euclidean distance matrix (D * D)


% statistical
\newcommand{\widebar}[1]{\overline{#1}}  

% 
% operators
\newcommand{\Had}{\circ}
\DeclareMathOperator*{\lmin}{Minimize}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\arginf}{arg\,inf}
\DeclareMathOperator*{\argsup}{arg\,sup}
%\newcommand*{\arginf}{\operatorname*{arginf}}
%\newcommand*{\argsup}{\operatorname*{argsup}}

% Sets
\newcommand*{\intersect}{\cap}
\newcommand*{\union}{\cup}
\let\oldemptyset\emptyset
\let\emptyset\varnothing

% Fields, Reals, etc. etc
\newcommand{\field}[1]{\mathbb{#1}}
\newcommand{\Reals}{\field{R}}
\newcommand{\Integers}{\field{Z}}
\newcommand{\Naturals}{\field{N}}
\newcommand{\Complex}{\field{C}}
\newcommand{\Rationals}{\field{Q}}

% Hyphenation
\hyphenation{Ar-chi-me-dean}

% Editorial
\newcommand*{\TODO}[1]{\textcolor{red}{TODO: #1}}
\newcommand*{\NOTE}[1]{\textcolor{blue}{Note: #1}}
\newtheorem{theorem}{Theorem}
\newtheorem{Proof}{Proof}

% Misc
%\makeatletter
%\newcommand\myisodate{\number\year-\ifcase\month\or 01\or 02\or 03\or 04\or 05\or 06\or 07\or 08\or 09\or 10\or 11\or 12\fi-\ifcase\day\or 01\or 02\or 03\or 04\or 05\or 06\or 07\or 08\or 09\or 10\or 11\or 12\or 13\or 14\or 15\or 16\or 17\or 18\or 19\or 20\or 21\or 22\or 23\or 24\or 25\or 26\or 27\or 28\or 29\or 30\or 31\fi}% create iso date
\makeatother
%\newcommand*{\abstractnoindent}{}% define abstract such that it has no indent
%\let\abstractnoindent\abstract
%\renewcommand*{\abstract}{\let\quotation\quote\let\endquotation\endquote
%  \abstractnoindent}
%\deffootnote[1em]{1em}{1em}{\textsuperscript{\thefootnotemark}}% setting for footnote
\pdfstringdefDisableCommands{\let\bm\relax}% allow \bm{} in ((sub)sub)section headings

%%% Document %%%

\author{Adam Rahman\\University of Waterloo}
\Plainauthor{Adam Rahman}
\title{\pkg{sdpt3r} - Semidefinite Quadratic Linear Programming in \proglang{R}}
\Plaintitle{sdpt3r - Semidefinite Quadratic Linear Programming in R}
\Shorttitle{sdpt3r - SQLP in R}
\Abstract{
We present the package \pkg{sdpt3r}, an \proglang{R} implementation of the \proglang{MATLAB} package \pkg{SDPT3} (\cite{toh1999sdpt3}). The purpose of the software is to solve semidefinite quadratic linear programming (SQLP) problems, a framework in which many important statistical problems can be cast. These include 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, among others.

While some tools exist in \proglang{R} to solve individual problems, such as the \code{nearPD} function in the \pkg{Matrix} package (\cite{bates2017package}) which solves the nearest correlation matrix problem, there currently does not exist a general solver for the SQLP class of problems, making \pkg{sdpt3r} a novel addition to the \proglang{R} library.
}

\Keywords{R, semidefinite programming, linear programming, quadratic programming, nearest correlation matrix, D-optimal experimental design, maximum cut, Lovasz number}

\Address{
  Adam Rahman\\
  Department of Statistics and Actuarial Science\\
  Faculty of Mathematics\\
  University of Waterloo\\
  200 University Ave W, Waterloo, Ontario, Canada\\
  E-mail: \email{a45rahma@uwaterloo.ca} \\
}

%% need no \usepackage{Sweave.sty}

\begin{document}
\SweaveOpts{concordance=FALSE}

\section{Introduction}\label{sec:intro}

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 \emph{conic linear optimization problems}.

One such example would be the nearest correlation matrix problem (\cite{higham2002computing}), 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 \emph{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 (\cite{smith1918standard}), classification using distance weighted discrimination (\cite{marron2007distance}), minimum volume ellipsoids (\cite{john2014extremum}), and problems in educational testing (\cite{chu1995educational}).

Important 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 (\cite{vandenberghe1998determinant}), 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 \proglang{R} package \pkg{sdpt3r}, an implementation of the \proglang{MATLAB} package \pkg{SDPT3} by \cite{toh1999sdpt3}. While there are currently functions in \proglang{R} available to solve some of the specific problems mentioned above - for instance the function \code{nearPD} in the \pkg{Matrix} package by \cite{bates2017package} solves the nearest correlation matrix problem - there currently does not exist a general solver for conic linear optimization in \proglang{R}, making \pkg{sdpt3r} a novel addition to the \proglang{R} library.

This paper is structured as follows. In Section \ref{sec:math} 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 \ref{sec:usage} discusses the \proglang{R} implementation of \pkg{sdpt3r}, and the main function by which conic linear optimization problems are solved, \code{sqlp}, including the required input, and the output generated. The same examples used in Section \ref{sec:math} will be used to demonstrate how a standard conic linear optimization problem can be converted to a form solvable by \code{sqlp}. Section \ref{sec:morex} presents the classic form of several other well known problems that can be solved using \pkg{sdpt3r}, as well as the helper functions available to convert them to the appropriate form. Finally Section \ref{sec:conc} provides some closing remarks.

\section{Conic linear optimization}
\label{sec:math}

At its simplest, a conic linear optimization problem has the following standard form (\cite{tutuncu2003solving}):

\begin{equation}
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & \langle \m{C},~\m{X} \rangle \\
\text{subject to} & \\
& 
\begin{array}{rll}
\langle \m{A}_{k},~\m{X} \rangle &=~~\ve{b}_{k}, & k=1,...,m \\
\m{X} &\in~~\mathcal{K} &
\end{array}
\end{array}
\label{eq:con1}
\end{equation}

\noindent where $\mathcal{K}$ is a cone. Generally, $\mathcal{K}$ is either a 

\begin{itemize}
\item Semidefinite Cone - $\mathcal{S}^{n} = \{\m{X} \in \mathcal{R}^{n\times n} : \m{X} \succeq 0, \m{X}_{ij} = \m{X}_{ji}~\forall~i\neq j\}$
\item Quadratic Cone - $\mathcal{Q}^{n} = \{\ve{x}=[x_{0};\tilde{\ve{x}}] \in \mathcal{R}^{n}: x_{0} \geq \sqrt{\tr{\tilde{\ve{x}}}\tilde{\ve{x}}}\}$
\item Linear Cone - $\mathcal{L}^{n}$ - non-negative orthant of $\mathcal{R}^{n}$
\end{itemize}

Here, $\tilde{\ve{x}} = [x_{1},\dots,x_{n-1}]$, and $\langle\cdot{,}\cdot\rangle$ represents the standard inner product in the appropriate space. In the semidefinite cone the inner product is $\langle \m{X},~\m{Y} \rangle = \tr{vec(\m{X})}vec(\m{Y})$, where the operator $vec$ is the by-column vector version of the matrix $\m{X}$, that is, for the $n\times n$ matrix $\m{X} = [x_{ij}]$, $vec(\m{X})$ is the $n^{2}\times 1$ vector $\tr{[x_{11},x_{12},x_{13},\dots,x_{(n-1)n},x_{nn}]}$. Note that $vec$ does not require a square matrix in general.

While not inherently statistical, 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 $\m{G} = [\m{V},\m{E}]$ be a graph with vertices $\m{V}$ and edges $\m{E}$. A \emph{cut} of the graph $\m{G}$ is a partition of the vertices of $\m{G}$ into two disjoint subsets $\m{G}_{1} = [\m{V}_{1},\m{E}_{1}]$, $\m{G}_{2} = [\m{V}_{2},\m{E}_{2}]$, with $\m{V}_{1} \cap \m{V}_{2} = \emptyset$. The size of the cut is defined to be the number of edges connecting the two subsets. The \emph{maximum cut} is defined to be the cut of a graph $\m{G}$ whose size is at least as large as any other cut. For a weighted graph object, we can also define the maximum cut to be the cut with weight at least as large as any other cut. 

Finding the maximum cut is referred to as the \textbf{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 (\cite{karp1972reducibility}). The Max-Cut problem is also known to be \emph{APX hard} (\cite{papadimitriou1991optimization}), meaning in addition to there being no polynomial time solution, there is also no polynomial time approximation. 

Using the semidefinite programming approximation formulation of \cite{goemans1995improved}, the Max-Cut problem can be approximated to within an \emph{approximation constant}. For a weighted adjacency matrix $\m{B}$, the objective function can be stated as

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & \langle \m{C},\m{X} \rangle \\
\text{subject to} & \\
& 
\begin{array}{rl}
diag(\m{X}) &=~~ \ve{1} \\
\m{X} & \in~~\mathcal{S}^{n} \\
\end{array}
\end{array}
\]

\noindent where $\mathcal{S}^{n}$ is the cone of symmetric positive semidefinite matrices of size $n$, and $\m{C} = -(diag(\m{B}\ve{1}) - \m{B})/4$. Here, we define $diag(\ve{a})$ for an $n \times 1$ vector $\ve{a}$ to be the diagonal matrix $\m{A} = [A_{ij}]$ of size $n \times n$ with $A_{ii} = a_{i},~~ i=1,\dots,n$. For a matrix $\m{X}$, $diag(\m{X})$ extracts the diagonal elements from $\m{X}$ and places them in a column-vector.

To see that the Max-Cut problem is a conic linear optimization problem it needs to be written in the same form as Equation \ref{eq:con1}. The objective function is already in a form identical to that of Equation \ref{eq:con1}, with minimization occurring over $\m{X}$ of its inner product with a constant matrix $\m{C} = -(diag(\m{B}\ve{1}) - \m{B})/4$. There are $n$ equality constraints of the form $x_{kk} = 1,~~k=1,...,n$, where $x_{kk}$ is the $k^{th}$ diagonal element of $\m{X}$, and $b_{k} = 1$ in Equation \ref{eq:con1}. To represent this in the form $\langle \m{A}_{k},~\m{X} \rangle = x_{kk}$,  take $\m{A}_{k}$ to be

\[
\m{A}_{k} = [a_{ij}] = 
\begin{cases}
1, & i = j = k \\
0, & \text{otherwise}
\end{cases}
\]

Now $\langle \m{A}_{k},~\m{X} \rangle = \tr{vec(\m{A}_{k})}vec(\m{X}) = x_{kk}$ as required, and the Max-Cut problem is specified as a conic linear optimization problem.

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.  

\subsection{A separable set of variables}\label{sec:gen2}

The conic linear optimization problem actually covers a much wider class of problems than those expressible as in Equation \ref{eq:con1}. Variables can be separated into those which are constrained to a semidefinite cone, $\mathcal{S}$, a quadratic cone, $\mathcal{Q}$, or a linear cone, $\mathcal{L}$. The objective function is a sum of the corresponding inner products of each set of variables. The linear constraint is simply a sum of variables of linear functions of each set. This more general  version of the conic linear optimization problem is

\begin{equation}
\begin{array}{ll}
\underset{\m{X}^{s},\m{X}^{q},\m{X}^{l}}{\text{minimize}} & \sum_{j=1}^{n_{s}} \langle \m{C}_{j}^{s},~\m{X}_{j}^{s} \rangle + \sum_{i=1}^{n_{q}} \langle \m{C}_{i}^{q},~\m{X}_{i}^{q} \rangle +  \langle \m{C}^{l},~\m{X}^{l} \rangle\\
\text{subject to} & \\
& 
\begin{array}{rl}
\sum_{j=1}^{n_{s}} \tr{(\m{A}_{j}^{s})}svec(\m{X}_{j}^{s}) + \sum_{i=1}^{n_{q}} \tr{(\m{A}_{i}^{q})}\m{X}_{i}^{q} + \tr{(\m{A}^{l})}\m{X}^{l} &=~~\ve{b} \\[5pt]
\m{X}_{j}^{s} &\in~\mathcal{S}^{s_{j}}~~\forall~j \\[5pt]
\m{X}_{i}^{q} &\in ~\mathcal{Q}^{q_{i}}~~\forall~i \\[5pt]
\end{array}
\end{array}
\label{eq:con2}
\end{equation}

Here, \emph{svec} takes the upper triangular elements of a matrix (including the diagonal) in a column-wise fashion and vectorizes them. In general for an $n \times p$ matrix $\m{X}=[x_{ij}]$, $svec(\m{X})$ will have the following form $\tr{[x_{11},x_{12},x_{22},x_{13},...,x_{(n-1)p},x_{np}]}$. Recall that matrices in $\mathcal{S}$ are symmetric, so it is sufficient to constrain only the upper triangular elements of the matrix $\m{X}^{s}$. For this formulation, $\m{A}_{j}^{s}$, $\m{A}_{i}^{q}$ and $\m{A}^{l}$ are the constraint matrices of the appropriate size.

Some important problems in statistics can be formulated to fit this form of the optimization problem. 

\subsubsection{The nearest correlation matrix}

First addressed by \cite{higham2002computing} 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 \cite{higham2002computing} as

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & ||\m{R} - \m{X}||_{F}\\
\text{subject to} & \\
& \begin{array}{rl}
diag(\m{X}) &=~~\ve{1} \\
\m{X} &\in~~\mathcal{S}^{n}
\end{array}
\end{array}
\]

\noindent where $\m{R}$ is the approximate correlation matrix and $||\cdot||_{F}$ denotes the Frobenius norm. Unfortunately, the Frobenius norm in the objective function prevents the problem being formatted as a conic linear optimization problem.

Since the matrix $\m{X}$ is constrained to have unit diagonal and be symmetric, and the matrix $\m{R}$ is an approximate correlation matrix, meaning it will also have unit diagonal and be symmetric, we can re-write the objective function as

\[
||\m{R}-\m{X}||_{F} = 2*||svec(\m{R}) - svec(\m{X})|| = 2*||\ve{e}||
\]

Now, introduce a variable $e_{0}$ such that $e_{0} \geq ||\ve{e}||$, and define $\ve{e}^{\ast} = [e_{0};~\ve{e}]$. The vector $\ve{e}^{\ast}$ is now restricted to be in the quadratic cone $\mathcal{Q}^{n(n+1)/2+1}$. This work leads to the formulation of \cite{toh1999sdpt3}

\[
\begin{array}{ll}
\underset{\ve{e}^{\ast},~\m{X}}{\text{minimize}} & e_{0} \\
\text{subject to} & \\
& 
\begin{array}{rl}
svec(\m{R}) - svec(\m{X}) &=~~ [\ve{0}, \m{I}_{n(n+1)/2}]~\ve{e}^{\ast} \\
diag(\m{X}) &=~~\ve{1} \\
\m{X} &\in~~\mathcal{S}^{n} \\
\ve{e}^{\ast} &\in~~\mathcal{Q}^{n(n+1)/2+1}
\end{array}
\end{array}
\]

Here, $[\m{X},\m{Y}]$ denotes column binding of the two matrices $\m{X}_{n\times p}$ and $\m{Y}_{n \times m}$ to form a matrix of size $n \times (p+m)$. By minimizing $e_{0}$, we indirectly minimize $\ve{e} = svec(\m{R}) - svec(\m{X})$, since recall we have $e_{0} \geq ||\ve{e}||$, which is the goal of the original objective function.

To see this as a conic linear optimization problem, notice that $e_{0}$ can be written as $\langle \m{C}^{q},~\m{X}^{q} \rangle$ by letting $\m{C}^{q} = [1;\ve{0}_{n(n+1)/2}]$ and $\m{X}^{q} = \ve{e}^{\ast}$. Since the matrix $\m{X}$ (i.e. $\m{X}^{s}$) does not appear in the objective function, the matrix $\m{C}^{s}$ is an $n \times n$ matrix of zeros.

Re-writing the first constraint as

\[
svec(\m{X}) + [\ve{0}, \m{I}_{n(n+1)/2}]~\ve{e}^{\ast} =~~  svec(\m{R})
\]

\noindent we can easily define the constraint matrices and right hand side of the first constraint as

\[
\begin{array}{rl}
\m{A}_{1}^{s} &= \m{I}_{n(n+1)/2} \\
\m{A}_{1}^{q} &= [\ve{0}, \m{I}_{n(n+1)/2}] \\
\ve{b}_{1} &= svec(\m{R})
\end{array}
\]

The second constraint is identical to the constraint from the Max-Cut problem, where each diagonal element of $\m{X}$ is constrained to be equal to 1. Define $\ve{b}_{2} = \ve{1}$, and for the $k^{th}$ diagonal element of $\m{X}$, define the matrix $\m{A}_{k}$ as

\[
\m{A}_{k} = [a_{ij}] = 
\begin{cases}
1, & i = j = k \\
0, & \text{otherwise}
\end{cases}
\]

\noindent yielding $\langle \m{A}_{k},\m{X} \rangle = x_{kk}$. To write this as $\tr{(\m{A}_{2}^{s})}\m{X}^{s}$, define 

\[
\m{A}_{2}^{s} = [svec(\m{A}_{1}),...,svec(\m{A}_{n})]
\]

Since $\ve{e}^{\ast}$ does not appear in the second constraint, $\m{A}_{2}^{q} = \ve{0}_{n(n+1)/2+1}$. 

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 $\m{A}^{s} = [\m{A}_{1}^{s},~\m{A}_{2}^{s}]$, $\m{A}^{q} = [\m{A}_{1}^{q},~\m{A}_{2}^{q}]$. We also concatenate both right hand side vectors to form a single vector by defining $\ve{b} = [\ve{b}_{1};~\ve{b}_{2}]$. Here, the notation $[\m{X};\m{Y}]$ is used to denote two matrices $\m{X}_{p \times m}$ and $\m{Y}_{q \times m}$ bound vertically to form a matrix of size $(p+q) \times m$. With this, the nearest correlation matrix problem is written as a conic linear optimization.

\newpage

\subsection{Semidefinite quadratic linear programming}\label{sec:gen3}

While Equation \ref{eq:con2} 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 \emph{semidefinite quadratic linear programming} (SQLP) problem.

The first generality afforded by an SQLP is the addition of an unconstrained variable $\m{X}^{u}$, which, as the name suggests, is not bound to a cone, but instead, it is ``constrained'' to the reals in the appropriate dimension. The second generalization is to allow for what are known as \emph{log-barrier} terms to exist in the objective function. In general, a barrier function in an optimization problem is a term that approaches infinity as the point approaches the boundary of the feasible region. As we will see, these log-barrier terms appear as log terms in the objective function.

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

\subsubsection{The primal problem}

The primal formulation of the SQLP problem is

\begin{equation}
\begin{array}{ll}
\underset{\m{X}_{j}^{s},\m{X}_{i}^{q},\m{X}^{l},\m{X}^{u}}{\text{minimize}} &  \sum_{j=1}^{n_{s}} [\langle \m{C}_{j}^{s},~ \m{X}_{j}^{s} \rangle - \ve{v}_{j}^{s} ~ log~det~\m{X}_{j}^{s}] ~ + ~ \sum_{i=1}^{n_{q}} [\langle \m{C}_{i}^{q},~ \m{X}_{i}^{q} \rangle - \ve{v}_{i}^{q} ~ log ~\gamma(\m{X}_{i}^{q})] \\
& \quad \quad \quad + ~ \langle \m{C}^{l},~ \m{X}^{l} \rangle ~ - ~ \sum_{k=1}^{n_{l}}\ve{v}_{k}^{l}~log~\m{X}_{k}^{l} ~ + ~ \langle \m{C}^{u},~ \m{X}^{u} \rangle \\
\text{subject to} & \\
& \begin{array}{rll}
\sum_{j=1}^{n_{s}} \m{A}_{j}^{s}(\m{X}_{j}^{s}) + \sum_{i=1}^{n_{q}} \m{A}_{i}^{q}\m{X}_{i}^{q} + \m{A}^{l}\m{X}^{l} + \m{A}^{u}\m{X}^{u} & = ~ \ve{b} & \\[5pt]
\m{X}_{j}^{s} & \in ~ \mathcal{S}^{s_{j}} & \forall~j \\[5pt]
\m{X}_{i}^{q} & \in ~ \mathcal{Q}^{q_{i}} & \forall~i \\[5pt] 
\m{X}^{l} & \in ~ \mathcal{L}^{n_{l}} & \\[5pt]
\m{X}^{u} & \in ~ \mathcal{R}^{n_{u}} &
\end{array}
\end{array}
\label{eq:Primal}
\end{equation}

For each $j$, $\m{C}_{j}^{s}$ and $\m{X}_{j}^{s}$ are symmetric matrices of dimension $s_{j}$, restricted to the cone of positive semidefinite matrices of the same dimension. Similarly, for all $i$, $\m{C}_{i}^{q}$ and $\m{X}_{i}^{q}$ are real vectors of dimension $q_{i}$, restricted to the the quadratic cone of dimension $q_{i}$. For a vector $\ve{u} = [u_{0}; \tilde{\ve{u}}]$ in a second order cone, define $\gamma(u) = \sqrt{u_{0}^{2} - \tr{\tilde{\ve{u}}}\tilde{\ve{u}}}$. Finally, $\m{C}^{l}$ and $\m{X}^{l}$ are vectors of dimension $n_{l}$, restricted to linear cone of the same dimension, and $\m{C}^{u}$ and $\m{X}^{u}$ are unrestricted real vectors of dimension $n_{u}$.

As before, the matrices $\m{A}_{i}^{q}$, $\m{A}^{l}$, and $\m{A}^{u}$ are constraint matrices in $q_{i}$, $n_{l}$, and $n_{u}$ dimensions respectively, each corresponding to their respective quadratic, linear, or unrestricted block. $\m{A}_{j}^{s}$ is defined to be a linear map from $\mathcal{S}^{s_{j}}$ to $\mathcal{R}^{m}$ defined by

\[
\m{A}_{j}^{s_{j}}(\m{X}_{j}^{s}) = [\langle \m{A}_{j,1}^{s}, \m{X}_{j}^{s} \rangle; \dots ; \langle \m{A}_{j,m}^{s}, \m{X}_{j}^{s} \rangle]
\]

where $\m{A}_{j,1}^{s} \dots \m{A}_{j,m}^{s} \in \mathcal{S}^{s_{j}}$ are constraint matrices associated with the $j^{th}$ semidefinite variable $\m{X}_{j}^{s}$. 

\subsubsection{The dual problem}

The dual problem associated with the semidefinite quadratic linear programming formulation is

\begin{equation}
\begin{array}{lll}
\underset{\m{Z}^{s}_{j},\m{Z}^{q}_{i},\m{Z}^{l},\ve{y}}{\text{maximize}} & \tr{\ve{b}}\ve{y} &+~\sum_{j=1}^{n_{s}}[\ve{v}_{j}^{s}~log~det~\m{Z}_{j}^{s} ~ + ~ s_{j}~\ve{v}_{j}^{s}~(1-log~\ve{v}_{j}^{s})] \\ 
& & +~\sum_{i=1}^{n_{q}} [\ve{v}_{i}^{q}~log~\gamma(\m{Z}_{i}^{q}) ~+~ \ve{v}_{i}^{q}~(1-log~\ve{v}_{i}^{q})] \\[10pt]
& & +~\sum_{k=1}^{n_{l}}[\ve{v}_{k}^{l}~log~\m{Z}_{k}^{l} ~ + ~ \ve{v}_{k}^{l}~(1-log~\ve{v}_{k}^{l})] \\
\text{subject to} & & \\
& & \begin{array}{llllll}
\tr{(\m{A}_{j}^{s})}\ve{y} & +~~ \m{Z}_{j}^{s} & = ~ \m{C}_{j}^{s}, & \m{Z}_{j}^{s} & \in ~ \mathcal{S}^{s_{j}}, & j = 1,\dots, n_{s} \\[5pt]
\tr{(\m{A}_{i}^{q})}\ve{y} & +~~ \m{Z}_{i}^{q} & = ~ \m{C}_{i}^{q}, & \m{Z}_{i}^{q} & \in ~ \mathcal{Q}^{q_{i}}, &  i = 1,\dots, n_{q} \\[5pt]
\tr{(\m{A}^{l})}\ve{y} & +~~ \m{Z}^{l} & = ~ \m{C}^{l}, & \m{Z}^{l} & \in ~ \mathcal{L}^{n_{l}} & \\[5pt]
\tr{(\m{A}^{u})}\ve{y} & & = ~ \m{C}^{u}, & \ve{y} & \in ~ \mathcal{R}^{m} & \\
\end{array}
\end{array}
\label{eq:Dual}
\end{equation}

\noindent where $(\m{A}_{j}^{s})^{T}$ is defined to be the adjoint operator of $\m{A}_{j}^{s}$, where $\tr{(\m{A}_{j}^{s})}\ve{y} = \sum_{k=1}^{m}\ve{y}_{k}\m{A}_{j,k}^{s}$. Equations \ref{eq:Primal} and \ref{eq:Dual} represent the most general form of the linear conic optimization problem that can be solved using \pkg{sdpt3r}.

\subsubsection{Optimal design of experiments}

Consider the problem of estimating a vector $\ve{x}$ from measurements $\ve{y}$ given by the relationship

\[
\ve{y} = \m{A}\ve{x} + \sv{\epsilon}, \quad \sv{\epsilon} \sim \mathcal{N}(0,1). 
\]

The variance-covariance matrix of such an estimator is proportional to $(\tr{\m{A}}\ve{A})^{-1}$. A reasonable goal during the design phase of an experiment would therefore be to minimize $(\tr{\m{A}}\ve{A})^{-1}$ in some way. 

There are many different ways in which $(\tr{\m{A}}\ve{A})^{-1}$ might be made minimal. For example, minimization of the trace of $(\tr{\m{A}}\ve{A})^{-1}$ (A-Optimality), minimization of the maximum eigenvalue of $(\tr{\m{A}}\ve{A})^{-1}$ (E-Optimality), minimization of the determinant of $(\tr{\m{A}}\ve{A})^{-1}$ (D-Optimilaity), and maximization of the trace of $\tr{\m{A}}\ve{A}$ (T-Optimality) all have their merits. 

Perhaps the most commonly used of these optimality criteria is D-Optimality, which is equivalent to maximizing the determinant of $\tr{\m{A}}\ve{A}$. Typically, the rows of $\m{A} = [\ve{a}_{1},...,\ve{a}_{q}]^{T}$ are chosen from $M$ possible test vectors $\ve{u}_{i} \in \mathcal{R}^{p}, ~ i=1,...M$, which are known in advance. That is,

\[
\ve{a}_{i} \in \{\ve{u}_{1},...,\ve{u}_{M}\}, \quad i = 1,...,q
\]

Given that the matrix $\m{A}$ is made up of these test vectors $\ve{u}_{i}$, \cite{vandenberghe1998determinant} write the matrix $\tr{\m{A}}\ve{A}$ as 

\begin{equation}
\tr{\m{A}}\ve{A} = q \sum_{i=1}^{M}\lambda_{i}\ve{u}_{i}\tr{\ve{u}_{i}}
\label{eq:dopt}
\end{equation}

\noindent where $\lambda_{i}$ is the fraction of rows in $\m{A}$ that are equal to the vector $\ve{u}_{i}$.  Then, \cite{vandenberghe1998determinant} write the D-optimal experimental design problem as a minimum determinant problem

\[
\begin{array}{ll}
\underset{\sv{\lambda}}{\text{minimize}} & \text{log}~\text{det}~(\sum_{i=1}^{M}\lambda_{i}\ve{u}_{i}\tr{\ve{u}_{i}})^{-1} \\
\text{subject to} & \\
& 
\begin{array}{rl}
\lambda_{i} &\geq~~0, \quad i = 1,...,m \\
\sum_{i=1}^{M}\lambda_{i} &=~~1
\end{array}
\end{array}
\]

Due to the inequality constraint, this primal formulation cannot be interpreted as an SQLP of the form of Equation \ref{eq:Primal}. By defining $\m{Z} = \ve{u}~diag(\sv{\lambda})~\tr{\ve{u}}$, the dual problem is (\cite{toh1999sdpt3}) 

\[
\begin{array}{ll}
\underset{\m{Z},~\ve{z}^{l},~\sv{\lambda}}{\text{maximize}} & \text{log det}~(\m{Z})\\
\text{subject to} &\\
& 
\begin{array}{rll}
-\sum_{i=1}^{p} \lambda_{i}(\ve{u}_{i}\tr{\ve{u}_{i}}) + \m{Z} &=~~0, & \m{Z} \in \mathcal{S}^{n} \\
-\sv{\lambda} + \ve{z}^{l} & =~~0, & \ve{z}^{l} \in \mathcal{R}^{p}_{+} \\
\ve{1}^{T}\sv{\lambda} & =~~1, & \sv{\lambda} \in \mathcal{R}^{p} \\
\end{array}
\end{array}
\]

Keeping in mind that this is a dual configuration, and thus follows Equation \ref{eq:Dual}, 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 $\m{Z}$, which is the log-barrier. This indicates that the variable $v^{s}$ in Equation \ref{eq:Dual} is equal to 1 in this formulation, while $v^{q}$ and $v^{l}$ are both zero. Since $\sv{\lambda}$ does not appear in the objective function, the vector $\ve{b}$ is equal to $\ve{0}$.

The constraint matrices $\m{A}$ are easy to define in the case of the dual formulation, as they multiply the vector $\ve{y}$ in Equation \ref{eq:Dual}, so therefore multiply $\sv{\lambda}$ in our case. In the first constraint, each $\lambda_{i}$ is multiplied by the matrix formed by $-\ve{u}_{i}\tr{\ve{u}}_{i}$, so define $\m{A}_{i}$ to be

\[
\m{A}_{i} = -\ve{u}_{i}\tr{\ve{u}}_{i},~~ i = 1,...,p.
\]

Then, the constraint matrix is $\m{A}^{s} = [svec(\m{A}_{1}),...,svec(\m{A}_{p})]$. In the second constraint containing the linear variable $\ve{z}^{l}$, the constraint matrix is $\m{A}^{l} = -\m{I}_{p}$, and in the third constraint containing only the unconstrained variable $\sv{\lambda}$, the constraint matrix is $\m{A}^{u} = \tr{\ve{1}}$. Since there is no quadratic variable, $A^{q} = \ve{0}$.

Finally, define the right hand side of each constraint
\[
\begin{array}{rl}
\m{C}^{s} &= \m{0}_{n\times n} \\
\m{C}^{l} &= \ve{0}_{p\times 1} \\
\m{C}^{u} &= 1
\end{array}
\]

\noindent which fully specifies the D-Optimal design problem as an SQLP. 

In the next section, we will demonstrate using \proglang{R} how these definitions can be translated for use in the main function of \pkg{sdpt3} so an SQLP problem can be solved.

\section[Examples using R]{Examples using \proglang{R}} \label{sec:usage}

In this section, we introduce the main function call available in \pkg{sdpt3r} and use the work done in Section \ref{sec:math} to demonstrate how an SQLP can be converted to a form solvable by \pkg{sdpt3r}. Using three numerical examples, we will also demonstrate the output provided by \pkg{sdpt3r}.

For \pkg{sdpt3r}, each optimization variable will be referred to as a \emph{block} in the space in which it is restricted. For instance, if we have an optimization variable $\m{X} \in \mathcal{S}^{n}$, we will refer to this as a semidefinite block of size $n$. It is important to note that it is possible to have multiple blocks from the same space, that is, it is possible to have both $\m{X} \in \mathcal{S}^{n}$ as well as $\m{Y} \in \mathcal{S}^{m}$ in the same problem.

\subsection{Input variables}

The main function call in \pkg{sdpt3r} is \code{sqlp}, which takes the following input variables

\begin{table}[h!]
\centering
\begin{tabular}{rl}
\code{blk} & A matrix object describing the block structure of the optimization variables. \\
\code{At} & A matrix object containing constraint matrices $\m{A}^{s}$, $\m{A}^{q}$, $\m{A}^{l}$, and $\m{A}^{u}$ \\
 & for the primal-dual problem.  \\
\code{b} & A vector containing the right hand side of the equality constraints, $\ve{b}$, \\
  & in the primal problem, or equivalently the constant vector in the dual. \\
\code{C} & A matrix object containing the constant $\m{C}$ matrices in the primal objective \\
  & function or equivalently the corresponding right hand side of the equality \\
  & constraints in the dual problem. \\
\code{X0}, \code{y0}, \code{Z0} & Matrix objects containing an initial iterate for the $\m{X}$, $\ve{y}$, and $\m{Z}$ variables for \\
   & the SQLP problem. If not provided, an initial iterate is computed internally. \\
\code{OPTIONS} & A list object providing additional parameters for use in \code{sqlpr}. \\
 & If not provided, default values are used.
\end{tabular}
\end{table}

The input variable \code{blk} describes the block structure of the problem. Letting \code{L} be the total number of semidefinite, quadratic, linear, and unrestricted blocks in the SQLP problem, define \code{blk} to be an $L \times 2$ matrix object, with the first column describing the type of block, and the second denoting the size of the optimization variable, summarized in Table \ref{tab:blk}.

\begin{table}[h!]
\centering
\begin{tabular}{|l|c|c|}
\hline
Block Type & Column 1 & Column 2 \\ \hline
semidefinite & s & $s_{j}$ \\
quadratic & q & $q_{i}$ \\
linear & l & $n_{l}$ \\
unrestricted & u & $n_{u}$ \\ \hline
\end{tabular}
\caption{\emph{Structure of \code{blk}.}}
\label{tab:blk}
\end{table}

The input variable \code{At} corresponds to the constraint matrices in Equation \ref{eq:Primal}, and \code{C} the constant matrices in the objective function. The size of these input variables depends on the block they are representing, summarized in Table \ref{tab:size1} for each block type.

\begin{table}[hbtp]
\centering
\begin{tabular}{|r|cccc|}
\hline
 & \multicolumn{4}{c|}{Block Type} \\
 & semidefinite & quadratic & linear & unrestricted \\ \hline
\code{At} & $\bar{s}_{j} \times m$ & $q_{j} \times m$ & $n_{l} \times m$ & $n_{u} \times m$ \\
\code{C} & $s_{j} \times s_{j}$ & $q_{j} \times 1$ & $n_{l} \times 1$ & $n_{u} \times 1$ \\ \hline
\end{tabular}
\caption{\emph{Size of \code{At} and \code{C} for each block type.}}
\label{tab:size1}
\end{table}

Note that in Table \ref{tab:size1}, $\bar{s}_{j} = s_{j}(s_{j} + 1)/2$. The size of \code{At} in the semidefinite block reflects the upper-triangular input format that has been discussed previously. In a semidefinite block, the optimization variable $\m{X}$ is necessarily symmetric and positive semidefinite, it is therefore more efficient to consider only the upper-triangular portion of the corresponding constraint matrix.

It is important to note that both input variables \code{At} and \code{C} are \emph{matrices of matrices}, a constraint matrix and a constant matrix for each optimization variable. While \code{blk} does not have this same requirement, we will nonetheless use the same structure to initialize \code{blk} as we do \code{At} and \code{C} for the sake of consistency.

In general, the user need not supply initial iterates \code{X0}, \code{y0}, and \code{Z0} for a solution to be found using \code{sqlp}. The infeasible starting point generated internally by \code{sqlp} tends to be sufficient to find a solution. If the user wishes to provide a starting point however, the size parameters in Table \ref{tab:init} must be met for each block.

\begin{table}[hbtp]
\centering
\begin{tabular}{|r|cccc|}
\hline
 & \multicolumn{4}{c|}{Block Type} \\
 & semidefinite & quadratic & linear & unrestricted \\ \hline
\code{X0} & $s_{j} \times s_{j}$ & $q_{j} \times 1$ & $n_{l} \times 1$ & $n_{u} \times 1$ \\
\code{y0} & $s_{j} \times 1$ & $q_{j} \times 1$ & $n_{l} \times 1$ & $n_{u} \times 1$ \\
\code{Z0} & $s_{j} \times s_{j}$ & $q_{j} \times 1$ & $n_{l} \times 1$ & $n_{u} \times 1$ \\ \hline
\end{tabular}
\caption{\emph{Required size for initial iterates \code{X0}, \code{y0}, and \code{Z0}.}}
\label{tab:init}
\end{table}

The user may choose to depart from the default values of several parameters which could affect the optimization by specifying alternative values in the \code{OPTIONS} list. A complete list of all parameters that can be altered can be found in Appendix \ref{sec:A}.

An important example is the specification of the \code{parbarrier} parameter in \code{OPTIONS}, which specifies the presence of a log-barrier in the objective function. The default case in \code{OPTIONS} assumes that the parameters $\ve{v}_{j}^{s},~\ve{v}_{i}^{q},~\ve{v}_{k}^{l}$ in Equation \ref{eq:Primal} are all $\ve{0}$. If this, however, is not the case, then the user must specify an $L \times 1$ matrix object in \code{OPTIONS$parbarrier} to store the values of these parameters (including zeros). If the $j^{th}$ block is a semidefinite block containing $p$ variables, $parbarrier_{j} = [v_{j1}^{s},...,v_{jn}^{s}]$. If the $j^{th}$ block is a quadratic block containing $p$ variables, $parbarrier_{j} = [v_{j1}^{q},...,v_{jn}^{q}]$. If the $j^{th}$ block is a linear block $parbarrier_{j} = [v_{1}^{l},...,v_{n_{l}}^{l}]$. Finally, if the $j^{th}$ block is the unrestricted block, then $parbarrier_{j} = [0,...,0]$, where 0 is repeated $n_{u}$ times. Section \ref{sec:dopt} contains an example where \code{OPTIONS$parbarrier} is specified.

When executed, \code{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. Both solutions are accessible through the following list of output variables

\newpage

\begin{table}[hbtp]
\centering
\begin{tabular}{rl}
\code{pobj} & the value of the primary objective function \\
\code{dobj} & the value of the dual objective function \\
\code{X} & A matrix object containing the optimal matrix $\m{X}$ for the primary problem \\
\code{y} & A vector object containing the optimal vector $\ve{y}$ for the dual problem \\
\code{Z} & A matrix object containing the optimal matrix $\m{Z}$ for the dual problem \\
\end{tabular}
\end{table}

The examples in subsequent subsections will demonstrate the output provided by \code{sqlp}.

\subsection{The Max-Cut problem} \label{sec:maxcut}

Recall that the maximum cut of a graph $\m{G}$ with adjacency matrix $\m{B}$ can be found as the solution to

\[
\begin{array}{ll}
\text{Minimize} & \langle \m{C},\m{X} \rangle \\
\text{subject to} & \\
& 
\begin{array}{rl}
diag(\m{X}) &=~~ \ve{1} \\
\m{X} & \in~~\mathcal{S}^{n} \\
\end{array}
\end{array}
\]

\noindent where $\m{C} = -(diag(\m{B}\ve{1}) - \m{B})/4$. In Section \ref{sec:math}, we wrote this in the form of an SQLP

\[
\begin{array}{ll}
\text{Minimize} & \langle \m{C},\m{X} \rangle \\
\text{subject to} & \\
& 
\begin{array}{rl}
\langle \m{A}_{k},\m{X} \rangle &=~~ 1, \quad k~=~1,\dots,n \\
\m{X} & \in~~\mathcal{S}^{n} \\
\end{array}
\end{array}
\]

where we defined $\m{A}_{k}$ as

\[
\m{A}_{k} = [a_{ij}] = 
\begin{cases}
1, & i = j = k \\
0, & \text{otherwise}
\end{cases}
\]

To convert this to a form usable by \code{sqlp}, we begin by noting that we have one optimization variable, $\m{X}$. Therefore, with $L=1$, we initialize the required input variables as follows

\begin{CodeChunk}
\begin{CodeInput}
R> blk <- matrix(list(),nrow=1,ncol=2)
R> At <- matrix(list(),nrow=1,ncol=1)
R> C <- matrix(list(),nrow=1,ncol=1)
\end{CodeInput}
\end{CodeChunk}

This initializes \code{blk}, \code{At}, and \code{C} as a matrix of matrices. While not required mathematically in this problem, this is the format required for all problems solved using \code{sqlp}, and will be required for any problem with more than one optimization variable. 

This initialization has the advantage of allowing \code{blk} to contain character and numerical values, without using the overhead of a data frame. Having $\m{X}$ constrained to the space of semidefinite matrices of size $n$, we specify \code{blk} as

\begin{CodeChunk}
\begin{CodeInput}
R> blk[[1,1]] <- "s"
R> blk[[1,2]] <- n
\end{CodeInput}
\end{CodeChunk}

With the objective function in the form $\langle \m{C}, \m{X} \rangle$, we define the input \code{C} as

\begin{CodeChunk}
\begin{CodeInput}
R> one <- matrix(1,nrow=n,ncol=1)
R> C[[1,1]] <- -(diag(B %*% one) - B)/4
\end{CodeInput}
\end{CodeChunk}

\noindent where $\m{B}$ is the adjacency matrix for a graph on which we would like to find the maximum cut, such as the one in Figure \ref{fig:maxcut}.

\begin{figure}[hbtp]
\begin{minipage}{.49\textwidth}
\includegraphics[width=.8\textwidth]{MaxCutGraph}
\end{minipage}
%
\begin{minipage}{.49\textwidth}
$\m{B} = \left[
\begin{array}{cccccccccc}
0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 1\\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\
1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\
1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 1\\
1 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1\\
0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 & 1 & 0\\
\end{array} \right]$
\end{minipage}
\caption{\emph{A graph object and associated adjacency matrix for which we would like to find the maximum cut.}}
\label{fig:maxcut}
\end{figure}

The matrix \code{At} is constructed using the upper triangular portion of the $\m{A}_{k}$ matrices. To do this in \proglang{R}, the function \code{svec} is made available in \pkg{sdpt3r}.

\begin{CodeChunk}
\begin{CodeInput}
R> #Construct Ak matrices
R> A <- matrix(list(),nrow=1,ncol=n)
R> for(k in 1:n){
R>   A[[k]] <- matrix(0,nrow=n,ncol=n)
R>   diag(A[[k]])[k] <- 1
R> }

R> #Combine to form At
R> At[[1,1]] <- svec(blk[1,],A)
\end{CodeInput}
\end{CodeChunk}

Having each of the diagonal elements of $\m{X}$ constrained to be 1, \code{b} is a $n \times 1$ matrix of ones

\begin{CodeChunk}
\begin{CodeInput}
R> b <- matrix(1,nrow=n,ncol=1)
\end{CodeInput}
\end{CodeChunk}

With all the input variables now defined, we can now call \code{sqlp} to solve the Max-Cut problem

\begin{CodeChunk}
\begin{CodeInput}
R> sqlp(blk, At, C, b)
\end{CodeInput}
\end{CodeChunk}

\subsubsection[A numerical example and the maxcut function]{A numerical example and the \code{maxcut} function}

The built-in function \code{maxcut} takes as input a (weighted) adjacency matrix \code{B} and returns all of the input variables required for use in \code{sqlp}. If we wish to find to the maximum cut of the graph in Figure \ref{fig:maxcut}, given the adjacency matrix $\m{B}$ we can compute the input variables for \code{sqlp} using \code{maxcut}

\begin{CodeChunk}
\begin{CodeInput}
R> out <- maxcut(B)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\begin{CodeOutput}
$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
\end{CodeOutput}
\end{CodeChunk}

Note that the value of the primary objective function is negative as we have defined $\m{C} = -(diag(\m{B}\ve{1}) - \m{B})/4$ since we require the primal formulation to be a minimization problem. The original formulation given in \cite{goemans1995improved} frames the Max-Cut problem as a maximization problem with $\m{C} = (diag(\m{B}\ve{1}) - \m{B})/4$. Therefore, the approximate value of the maximum cut for the graph in Figure \ref{fig:maxcut} is 14.68 (recall we are solving a relaxation). 

As an interesting aside, we can show that the matrix $\m{X}$ is actually a correlation matrix by considering its eigenvalues - we can see it clearly is symmetric, with unit diagonal and all elements in [-1,1].

\begin{CodeChunk}
\begin{CodeInput}
R> eigen(X)
\end{CodeInput}
\begin{CodeOutput}
$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
\end{CodeOutput}
\end{CodeChunk}

The fact that $\m{X}$ is indeed a correlation matrix comes as no surprise. \cite{goemans1995improved} show that the set of feasible solutions for the Max-Cut problem is in fact the set of correlation matrices. So while we may not be interested in $\m{X}$ as an output for solving the Max-Cut problem, it is nonetheless interesting to see that it is in fact in the set of feasible solutions.

\subsection{Nearest correlation matrix}

Recall that the nearest correlation matrix is found as the solution to

\[
\begin{array}{ll}
\underset{\ve{e}^{\ast},~\m{X}}{\text{minimize}} & e_{0} \\
\text{subject to} & \\
& 
\begin{array}{rl}
svec(\m{R}) - svec(\m{X}) &=~~ [\ve{0}, \m{I}_{n(n+1)/2}]~\ve{e}^{\ast} \\
diag(\m{X}) &=~~\ve{1} \\
\m{X} &\in~~\mathcal{S}^{n} \\
\ve{e}^{\ast} &\in~~\mathcal{Q}^{n(n+1)/2+1}
\end{array}
\end{array}
\]

In Section \ref{sec:gen2} we wrote this as the following SQLP

\[
\begin{array}{ll}
\underset{\ve{e}^{\ast},~\m{X}}{\text{minimize}} & \langle \m{C}, \ve{e}^{\ast} \rangle \\
\text{subject to} & \\
& 
\begin{array}{rl}
\tr{(\m{A}^{s})}svec(\m{X}) + \tr{(\m{A}^{q})}\ve{e}^{\ast} &=~~ \ve{b} \\
\m{X} &\in~~\mathcal{S}^{n} \\
\ve{e}^{\ast} &\in~~\mathcal{Q}^{n(n+1)/2+1}
\end{array}
\end{array}
\]

\noindent for $\m{C} = [1,\ve{0}_{n(n+1)/2}]$, and \\[10pt]
\begin{minipage}{\textwidth}
\begin{minipage}[t]{.49\textwidth}
\begin{center}
$
\begin{array}{rl}
\m{A}^{s} &=~~ [\m{A}_{1}^{s},~\m{A}_{2}^{s}]~~~~~~~~~~~~~~~~~~~~ \\
\m{A}^{q} &=~~ [\m{A}_{1}^{q},~\m{A}_{2}^{q}]~~~~~~~~~~~~~~~~~~~~ \\
\end{array}
$
\end{center}
\end{minipage}
%
\begin{minipage}[t]{.49\textwidth}
\begin{center}
$
\begin{array}{ll}
\ve{b} &=~ [\ve{b}_{1};~\ve{b}_{2}] \\
\end{array}
$
\end{center}
\end{minipage}
\end{minipage}
\\[10pt]
\noindent where
\\[10pt]
\begin{minipage}{\textwidth}
\begin{minipage}[t]{.49\textwidth}
\begin{center}
$
\begin{array}{rl}
\m{A}_{1}^{s} &=~~ \m{I}_{n_{2}} \\
\m{A}_{1}^{q} &=~~ [\ve{0}, \m{I}_{n_{2}}] \\
& \\
\m{A}_{2}^{s} &=~~ [svec(\m{A}_{1}),\dots,svec(\m{A}_{n})] \\
\m{A}_{2}^{q} &=~~ \m{0}_{n_{2}} \\
\end{array}
$
\end{center}
\end{minipage}
%
\begin{minipage}[t]{.49\textwidth}
\begin{center}
$
\begin{array}{ll}
\ve{b}_{1} &=~ svec(\m{R}) \\
\ve{b}_{2} &=~ \tr{\ve{1}} \\
\end{array}
$
\end{center}
\end{minipage}
\end{minipage}
\\[10pt]
\noindent and $\m{A}_{1},\dots,\m{A}_{n}$ are given by

\[
\m{A}_{k} = [a_{ij}] = 
\begin{cases}
1, & i = j = k \\
0, & \text{otherwise}
\end{cases}
\]

To be solved using \code{sqlp}, we first define \code{blk}. There are two optimization variables in the formulation of the nearest correlation matrix - $\m{X}$ is an $n \times n$ matrix constrained to be in a semidefinite cone, and $\ve{y}$ is an $n(n+1)/2 + 1$ length vector constrained to be in a quadratic cone, so 

\newpage

\begin{CodeChunk}
\begin{CodeInput}
R> blk <- matrix(list(),nrow=2,ncol=2)

R> blk[[1,1]] <- "s"
R> blk[[1,2]] <- n

R> blk[[2,1]] <- "q"
R> blk[[2,2]] <- n * (n + 1) / 2 + 1
\end{CodeInput}
\end{CodeChunk}

Note that $\m{X}$ does not appear in the objective function, so the \code{C} entry corresponding to the block variable $\m{X}$ is an $n \times n$ matrix of zeros, which defines \code{C} as

\begin{CodeChunk}
\begin{CodeInput}
R> C <- matrix(list(),nrow=2,ncol=1)
R> C[[1,1]] <- matrix(0,nrow=n,ncol=n)
R> C[[2,1]] <- rbind(1, matrix(0,nrow=n2,ncol=1))
\end{CodeInput}
\end{CodeChunk}

Next comes the constraint matrix for $\m{X}$

\begin{CodeChunk}
\begin{CodeInput}
R> At <- matrix(list(),nrow=2,ncol=1)
R>
R> #Constraint Matrix for Upper Triangular Elements of X
R> A1s <- diag(1,nrow=n2,ncol=n2)
R>
R> #Construct Ak matrices
R> Aks <- matrix(list(),nrow=1,ncol=n)
R> for(k in 1:n){
R>   Aks[[k]] <- matrix(0,nrow=n,ncol=n)
R>   diag(Aks[[k]])[k] <- 1
R> }
R>
R> A2s <- svec(blk[1,],Aks)
R> 
R> #Combined Constraint Matrix for X
R> At[[1,1]] <- cbind(A1s,A2s)
\end{CodeInput}
\end{CodeChunk}

\noindent then the constraint matrix for $\ve{e}^{\ast}$.

\begin{CodeChunk}
\begin{CodeInput}
R> A1q<- matrix(0,nrow=n,ncol=n2+1)
R>
R> A2q1 <- matrix(0,nrow=n2,ncol=1)
R> A2q2 <- diag(1,nrow=n2,ncol=n2)
R> A2q <- cbind(A211, A212)
R>
R> At[[2,1]] <- rbind(A1q, A2q)
\end{CodeInput}
\end{CodeChunk}

\noindent and the right hand side vector \code{b}

\begin{CodeChunk}
\begin{CodeInput}
R> b <- rbind(svec(blk[1,],R),matrix(1,n,1))
\end{CodeInput}
\end{CodeChunk}

The nearest correlation matrix problem is now solved by

\begin{CodeChunk}
\begin{CodeInput}
R> sqlp(blk, At, C, b)
\end{CodeInput}
\end{CodeChunk}

\subsubsection[A numerical example and the nearcorr function]{A numerical example and the \code{nearcorr} function}

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 \proglang{R} package \pkg{quantmod} (\cite{quantmod}) to access stock data from five tech firms (Microsoft, Apple, Amazon, Alphabet/Google, and IBM) beginning in 2007.

\begin{CodeChunk}
\begin{CodeInput}
R> library(quantmod)

R> getSymbols(c("MSFT", "AAPL", "AMZN", "GOOGL", "IBM"))
R> stock.close <- as.xts(merge(MSFT, AAPL, AMZN, GOOGL,IBM))[,c(4,10,16,22,28)]
\end{CodeInput}
\end{CodeChunk}

The correlation matrix for these five stocks is

\begin{CodeChunk}
\begin{CodeInput}
R> stock.corr <- cor(stock.close)
R> stock.corr
\end{CodeInput}
\begin{CodeOutput}
            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
\end{CodeOutput}
\end{CodeChunk}

Next, consider the effect of having a positive correlation between Microsoft and Apple

\begin{CodeChunk}
\begin{CodeInput}
R> stock.corr[1,2] <- -1 * stock.corr[1,2]
R> stock.corr[2,1] <- stock.corr[1,2]
R> stock.corr
\end{CodeInput}
\begin{CodeOutput}
            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
\end{CodeOutput}
\end{CodeChunk}

Unfortunately, this correlation matrix is not positive semidefinite

\begin{CodeChunk}
\begin{CodeInput}
R> eigen(stock.corr)$values
\end{CodeInput}
\begin{CodeOutput}
[1]  2.8850790  1.4306393  0.4902211  0.3294150 -0.1353544
\end{CodeOutput}
\end{CodeChunk}

Given the approximate correlation matrix \code{stock.corr}, the built-in function \code{nearcorr} provides the input variables required for the problem to be solved using \code{sqlp}

\begin{CodeChunk}
\begin{CodeInput}
R> out <- nearcorr(stock.corr)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> foo <- sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

Since this is a minimization problem, and thus a primal formulation of the SQLP, the output \code{X} from \code{sqlp} will provide the optimal solution to the problem - that is, \code{X} will be the nearest correlation matrix to \code{stock.corr}.

\begin{CodeChunk}
\begin{CodeOutput}
foo$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
\end{CodeOutput}
\end{CodeChunk}

The matrix above is symmetric with unit diagonal and all entries in $[-1,1]$. By checking the eigenvalues,

\begin{CodeChunk}
\begin{CodeInput}
eigen(X)
\end{CodeInput}
\begin{CodeOutput}
$values

[1] 2.846016e+00 1.384062e+00 4.570408e-01 3.128807e-01 9.680507e-11
\end{CodeOutput}
\end{CodeChunk}

\noindent we can see that \code{X} is indeed a correlation matrix. 

\subsection{D-optimal experimental design} \label{sec:dopt}

Recall from Section \ref{sec:gen3} that the D-Optimal experimental design problem was stated as the following dual SQLP

\[
\begin{array}{ll}
\underset{\m{Z},~\ve{z}^{l},~\sv{\lambda}}{\text{maximize}} & \text{log det}~(\m{Z})\\
\text{subject to} &\\
& 
\begin{array}{rll}
-\sum_{i=1}^{p} \lambda_{i}(\ve{u}_{i}\tr{\ve{u}_{i}}) + \m{Z} &=~~0, & \m{Z} \in \mathcal{S}^{n} \\
-\sv{\lambda} + \ve{z}^{l} & =~~0, & \ve{z}^{l} \in \mathcal{R}^{p}_{+} \\
\ve{1}^{T}\sv{\lambda} & =~~1, & \sv{\lambda} \in \mathcal{R}^{p} \\
\end{array}
\end{array}
\]

\noindent which we wrote as 

\[
\begin{array}{ll}
\underset{\m{Z},~\ve{z}^{l},~\sv{\lambda}}{\text{maximize}} & \text{log det}~(\m{Z})\\
\text{subject to} &\\
& 
\begin{array}{rll}
\tr{(\m{A}^{s})} \sv{\lambda} &+~~ \m{Z} &=~~\m{C}^{s}, \quad \m{Z} \in \mathcal{S}^{n} \\
\tr{(\m{A}^{l})} \sv{\lambda} &+~~ \ve{z}^{l} & =~~\m{C}^{q}, \quad \ve{z}^{l} \in \mathcal{R}^{p}_{+} \\
\tr{(\m{A}^{u})}\sv{\lambda} & & =~~\m{C}^{u}, \quad \sv{\lambda} \in \mathcal{R}^{p} \\
\end{array}
\end{array}
\]

\noindent where $\ve{b} = \ve{0}$, and
\\[10pt]
\begin{minipage}{\textwidth}
\begin{minipage}{.49\textwidth}
\begin{center}
$
\begin{array}{rl}
\m{A}^{s} &= -[svec(\m{A}_{1}),\dots,svec(\m{A}_{p})] \\
\m{A}^{l} &= -\m{I}_{p} \\
\m{A}^{u} &= \tr{\ve{1}}
\end{array}
$
\end{center}
\end{minipage}
%
\begin{minipage}{.49\textwidth}
\begin{center}
$
\begin{array}{rl}
\m{C}^{s} &= \m{0}_{n\times n} \\
\m{C}^{l} &= \ve{0}_{p\times 1} \\
\m{C}^{u} &= 1
\end{array}
$
\end{center}
\end{minipage}
\end{minipage}
\\[10pt]
Here, $\m{A}_{1},\dots,\m{A}_{p}$ are given by 

\[
\m{A}_{i} = \ve{u}_{i}\tr{\ve{u}_{i}}, \quad i = 1,\dots,p
\]

To convert this to a form usable by \pkg{sdpt3r}, we initialize our input variables by noting we have three blocks - $\m{X}$, $\ve{z}^{l}$, and $\sv{\lambda}$

\begin{CodeChunk}
\begin{CodeInput}
R> blk <- matrix(list(),nrow=3,ncol=2)
R> At <- matrix(list(),nrow=3,ncol=1)
R> C <- matrix(list(),nrow=3,ncol=1)
\end{CodeInput}
\end{CodeChunk}

As before, we declare the three blocks in \code{blk}. The first block is semidefinite containing the matrix $\m{Z}$, the second a linear block containing $\ve{z}^{l}$, and the third an unrestricted block containing $\sv{\lambda}$

\begin{CodeChunk}
\begin{CodeInput}
R> blk[[1,1]] <- "s"
R> blk[[1,2]] <- n

R> blk[[2,1]] <- "l"
R> blk[[2,2]] <- p

R> blk[[3,1]] <- "u"
R> blk[[3,2]] <- 1
\end{CodeInput}
\end{CodeChunk}

Next, by noting the variable $\sv{\lambda}$ does not appear in the objective function, we specify \code{b} as a vector of zeros

\begin{CodeChunk}
\begin{CodeInput}
R> b <- matrix(0,nrow=p,ncol=1)
\end{CodeInput}
\end{CodeChunk}

Next, looking at the right-hand side of the constraints, we define the matrices \code{C}

\begin{CodeChunk}
\begin{CodeInput}
R> C[[1,1]] <- matrix(0,nrow=n,ncol=n)
R> C[[2,1]] <- matrix(0,nrow=p,ncol=1)
R> C[[3,1]] <- 1
\end{CodeInput}
\end{CodeChunk}

Finally, we construct \code{At} for each variable

\begin{CodeChunk}
\begin{CodeInput}
R> A <- matrix(list(),nrow=p,ncol=1)
$>
R> for(k in 1:p){
R>   A[[k]] <- -uk %*% t(uk)
R> }
R>
R> At[[1,1]] <- svec(blk[1,], A)
R> At[[2,1]] <- diag(-1,nrow=p,ncol=p)
R> At[[3,1]] <- matrix(1,nrow=1,ncol=p)
\end{CodeInput}
\end{CodeChunk}

The final hurdle necessary to address in this problem is the existence of the log-barrier. Recall that it is assumed that $v^{s}, v^{q}$, and $v^{l}$ in Equation \ref{eq:Dual} are all zero in \code{OPTIONS}. In this case, we can see that is not true, as we have a log term containing $\m{Z}$ in the objective function, meaning $v^{s}$ is equal to one. To pass this to \code{sqlp}, we define the \code{OPTIONS$parbarrier} variable as

\begin{CodeChunk}
\begin{CodeInput}
R> OPTIONS$parbarrier <- matrix(list(),nrow=3,ncol=1)
R> OPTIONS$parbarrier[[1]] <- 1       #for vs
R> OPTIONS$parbarrier[[2]] <- 0       #for vq
R> OPTIONS$parbarrier[[3]] <- 0       #for vl
\end{CodeInput}
\end{CodeChunk}

The D-Optimal experimental design problem can now be solved using \code{sqlp}

\begin{CodeChunk}
\begin{CodeInput}
R> sqlp(blk, At, C, b, OPTIONS)
\end{CodeInput}
\end{CodeChunk}

\subsubsection[A numerical example and the doptimal function]{A numerical example and the \code{doptimal} function}

To demonstrate the output generated from a D-optimal experimental design problem, we consider a simple $3 \times 25$ matrix containing the known test vectors $\ve{u}_{1},...,\ve{u}_{25}$ (the data is available in the \code{sqlp} package). To generate the required input for \code{sqlp}, we use the function \code{doptimal}, which takes as input an $n \times p$ matrix $\m{U}$ containing the known test vectors, and returns the input necessary for \code{sqlp}. The output we are interested in is \code{y}, corresponding to $\sv{\lambda}$ in our formulation, the percentage of each $\ve{u}_{i}$ necessary to achieve maximum information in the experiment.

\begin{CodeChunk}
\begin{CodeInput}
R> data(DoptDesign)

R> out <- doptimal(DoptDesign)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b
R> OPTIONS <- out$OPTIONS

R> sqlp(blk,At,C,b,OPTIONS)
\end{CodeInput}
\begin{CodeOutput}
$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
\end{CodeOutput}
\end{CodeChunk}

The information matrix $\tr{\m{A}}\m{A}$ is a linear combination of the test vectors $\ve{u}_{i}$, weighted by the optimal vector \code{y} above. 

\section{Additional problems} \label{sec:morex}

The \pkg{sdpt3r} package considerably broadens the set of optimization problems that can be solved in \proglang{R}. In addition to those problems presented in detail in Section \ref{sec:usage}, 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 \ref{eq:Primal} or \ref{eq:Dual}. Accompanying each problem will be an \proglang{R} helper function, which will produce the input variables \code{blk}, \code{At}, \code{C}, and \code{b} so that the problem can be solved using \code{sqlp}. Each helper function is made available to the user in the \pkg{sdpt3r} package.

\subsection{Minimum volume ellipsoids}

The problem of finding the ellipsoid of minimum volume containing a set of points $\ve{v}_{1},...,\ve{v}_{n}$ is stated as the following optimization problem (\cite{vandenberghe1998determinant})

\[
\begin{array}{ll}
\underset{\m{B},~\ve{d}}{\text{maximize}} & log~det(\m{B}) \\
\text{subject to} & \\
& \begin{array}{rl}
||\m{B}\ve{x} + \ve{d}|| &\leq~~ 1, \quad \forall~]ve{x}~\in~[\ve{v}_{1},...,\ve{v}_{n}]
\end{array}
\end{array}
\]

The function \code{minelips} takes as input an $n \times p$ matrix $\m{V}$ containing the points around which we would like to find the minimum volume ellipsoid, and returns the input variables necessary to solve the problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- minelips(V)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b
R> OPTIONS <- out$OPTIONS

R> sqlp(blk,At,C,b,OPTIONS)
\end{CodeInput}
\end{CodeChunk}

\subsection{Distance weighted discrimination}

Given two sets of points in a matrix $\m{X} \in \mathcal{R}^{n}$ with associated class variables [-1,1] in $\m{Y} = diag(\ve{y})$, distance weighted discrimination (\cite{marron2007distance}) seeks to classify the points into two distinct subsets by finding a hyperplane between the two sets of points. Mathematically, the distance weighted discrimination problem seeks a hyperplane defined by a normal vector, $\sv{\omega}$, and position, $\beta$, such that each element in the residual vector $\bar{\ve{r}} = \m{Y}\tr{\m{X}}\sv{\omega} + \beta \ve{y}$ is positive and large. Since the class labels are either 1 or -1, having the residuals be positive is equivalent to having the points on the proper side of the hyperplane.

Of course, it may be impossible to have a perfect separation of points using a linear hyperplane, so an error term $\xi$ is introduced. Thus, the perturbed residuals are defined to be

\[
\ve{r} =  \m{Y}\tr{\m{X}}\sv{\omega} + \beta \ve{y} + \sv{\xi}
\]

Distance Weighted Discrimination (\cite{marron2007distance}) solves the following optimization problem to find the optimal hyperplane.

\[
\begin{array}{ll}
\underset{\ve{r},~\sv{\omega},~\beta,~\sv{\xi}}{\text{minimize}} & \sum_{i=1}^{n}(1/r_{i}) + C\tr{\ve{1}}\sv{\xi} \\
\text{subject to} & \\
& \begin{array}{rl}
\ve{r} &=~~ \m{Y}\tr{\m{X}}\sv{\omega} + \beta \ve{y} + \sv{\xi} \\
\tr{\sv{\omega}}\sv{\omega} &\leq~~ 1 \\
\ve{r} &\geq~~ \ve{0} \\
\sv{\xi} &\geq~~ \ve{0}
\end{array}
\end{array}
\]

\noindent where $C > 0$ is a penalty parameter to be chosen. 

The function \code{dwd} takes as input two $n \times p$ matrices \code{X1} and \code{X2} containing the points to be separated, as well as a penalty term \code{C} $\geq 0$ penalizing the movement of a point on the wrong side of the hyperplane to the proper side, and returns the input variables necessary for \code{sqlp} to solve the distance weighted discrimination problem. 

\begin{CodeChunk}
\begin{CodeInput}
R> out <- dwd(X1,X2,C)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{Max-kCut}

Similar to the Max-Cut problem, the Max-kCut problem asks, given a graph $\m{G} = (\m{V},\m{E})$ and an integer $k$, does a cut exist of at least size $k$. For a given (weighted) adjacency matrix $\m{B}$ and integer $k$, the Max-kCut problem is formulated as the following primal problem

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & \langle \m{C},~\m{X} \rangle \\
\text{subject to} & \\
 & \begin{array}{rl}
diag(\m{X}) &=~~\ve{1} \\
X_{ij} &\geq~~ 1/(k-1) \quad \forall~i \neq j\\
\m{X} &\in ~~\mathcal{S}_{n}
\end{array}
\end{array}
\]

Here, $\m{C} = -(1-1/k)/2 * (diag(\m{B}\ve{1})-\m{B})$. The Max-kCut problem is slightly more complex than the Max-Cut problem due to the inequality constraint. In order to turn this into a standard SQLP, we must replace the inequality constraints with equality constraints, which we do by introducing a slack variable $\ve{x}^{l}$, allowing the problem to be restated as

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}}  & \langle \m{C},~\m{X} \rangle \\
\text{subject to} & \\
 & \begin{array}{rl}
diag(\m{X}) &=~~\ve{1} \\
X_{ij} - x^{l} &=~~ 1/(k-1) \quad \forall~i \neq j\\
\m{X} &\in ~~\mathcal{S}^{n} \\
\ve{x}^{l} &\in ~~ \mathcal{L}^{n(n+1)/2} 
\end{array}
\end{array}
\]

The function \code{maxkcut} takes as input an adjacency matrix \code{B} and an integer \code{k}, and returns the input variables necessary for the problem to be solved using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- maxkcut(B,k)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{Graph partitioning problem}

The graph partitioning problem can be formulated as the following primal optimization problem

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & tr(\m{C}\m{X}) \\
\text{subject to} & \\
 & \begin{array}{rl}
tr(\ve{1}\tr{\ve{1}}\m{X}) &=~~ \alpha \\
diag(\m{X}) &=~~ \ve{1}
\end{array}
\end{array}
\]

Here, $\m{C} = -(diag(\m{B}\ve{1}) - \m{B})$, for an adjacency matrix $\m{B}$, and $\alpha$ is any real number. 

The function \code{gpp}, takes as input a weighted adjacency matrix \code{B} and a real number \code{alpha} and returns the input necessary to solve the problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- gpp(B,alpha)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{The Lovasz number}

The Lovasz Number of a graph $\m{G}$, denoted $\vartheta(\m{G})$, is the upper bound on the Shannon capacity of the graph. For an adjacency matrix $\m{B} = [B_{ij}]$ the problem of finding the Lovasz number is given by the following primal SQLP problem

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & tr(\m{C}\m{X}) \\
\text{subject to} & \\
& \begin{array}{rl}
tr(\m{X}) & = ~~1 \\
X_{ij} &=~~ 0 \quad \text{if $B_{ij}$ = 1} \\
\m{X} &\in~~ \mathcal{S}^{n}
\end{array}
\end{array}
\]

The function \code{lovasz} takes as input an adjacency matrix $\m{B}$, and returns the input variables necessary for the Lovasz number to be found using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- lovasz(B)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{Toeplitz approximation}

Given a symmetric matrix $\m{F}$, the Toeplitz approximation problem seeks to find the nearest symmetric positive definite Toeplitz matrix. In general, a Toeplitz matrix is one with constant descending diagonals, i.e.

\[
\m{T} = \left[
\begin{array}{ccccc}
a & b & c & d & e \\
f & a & b & c & d \\
g & f & a & b & c \\
h & g & f & a & b \\
i & h & g & f & a
\end{array}
\right]
\]

is a general Toeplitz matrix. The problem is formulated as the following optimization problem 

\[
\begin{array}{ll}
\underset{\m{X}}{\text{maximize}} & -y_{n+1} \\
\text{subject to} & \\
&
\begin{array}{rl}
\left[\begin{array}{cc}
\m{I} & \m{0} \\
\m{0} & -\beta
\end{array}\right] ~~+~~ \sum_{k=1}^{n}y_{k}
\left[\begin{array}{cc}
\m{0} & \gamma_{k}\ve{e}_{k} \\
\gamma_{k}\ve{e}_{k}^{T} & -2q_{k}
\end{array}\right] ~~+~~y_{n+1}\m{B} &\geq~~ \m{0} \\
\tr{[y_{1},...,y_{n}]} + y_{n+1}\m{B} & \geq~~ \ve{0}
\end{array}
\end{array}
\]

\noindent where $\m{B}$ is an $(n+1)\times(n+1)$ matrix of zeros, and $\m{B}_{(n+1)(n+1)} = 1$, $q_{1} = -tr(\m{F})$, $q_{k} =$ sum of $k^{th}$ diagonal upper and lower triangular matrix, $\gamma_{1} = \sqrt{n}$, $\gamma_{k} = \sqrt{2*(n-k+1)},~k=2,...,n$, and $\beta = ||\m{F}||_{F}^{2}$. 

The function \code{toep} takes as input a symmetric matrix \code{F} for which we would like to find the nearest Toeplitz matrix, and returns the input variables required to solve the problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- toep(F)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{The educational testing problem}

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 (\cite{fletcher1981nonlinear}). 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 $\m{S}$ such that the resulting matrix remains positive semidefinite (\cite{chu1995educational}).

The Educational Testing Problem (ETP) is formulated as the following dual problem

\[
\begin{array}{rl}
\underset{\ve{d}}{\text{maximize}} & \tr{\ve{1}}\ve{d} \\
\text{subject to} & \\
& \begin{array}{rl}
\m{A} - diag(\ve{d}) &\succeq~~ \ve{0} \\
\ve{d} &\geq~~ \ve{0}
\end{array}
\end{array}
\]

\noindent where $\ve{d} = [d_{1},~d_{2},...,~d_{n}]$ is a vector of size $n$ and $diag(\ve{d})$ denotes the corresponding $n \times n$ diagonal matrix. In the second constraint, having each element in $\ve{d}$ be greater than or equal to 0 is equivalent to having $diag(\ve{d}) \succeq 0$. 

The corresponding primal problem is

\[
\begin{array}{rl}
\underset{\m{X}}{\text{minimize}} & tr(\m{A}\m{X}) \\
\text{subject to} & \\
& \begin{array}{rl}
diag(\m{X}) &\geq~~ \ve{1} \\
\m{X} &\succeq~~ 0
\end{array}
\end{array}
\]

The function \code{etp} takes as input an $n \times n$ positive definite matrix \code{A}, and returns the input variables required to solve the educational testing problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- etp(A)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{Logarithmic Chebyshev approximation}

For a $p \times n$ ($p > n$) matrix $\m{B}$ and $p \times 1$ vector $\ve{f}$, the Logarithmic Chebyshev Approximation problem is stated as the following optimization problem (\cite{vandenberghe1998determinant})

\[
\begin{array}{ll}
\underset{\ve{x},~t}{\text{minimize}} & t \\
\text{subject to} & \\
 & \begin{array}{rl}
1/t & \leq~~(\tr{\m{x}}\m{B}_{i\cdot})/\ve{f}_{i}~~\leq~~t, \quad i = 1,...,p
\end{array}
\end{array}
\]

\noindent where $\m{B}_{i\cdot}$ denotes the $i^{th}$ row of the matrix $\m{B}$. Note that we require each element of $\m{B}_{\cdot j}/\ve{f}$ to be greater than or equal to 0 for all $j$. 

The function \code{logcheby} takes as input a matrix \code{B} and vector \code{f}, and returns the input variables necessary to solve the Logarithmic Chebyshev Approximation problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- logcheby(B,f)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\subsection{Linear matrix inequality problems}

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 $n\times p$ matrix $\m{B}$ known in advance

\[
\begin{array}{ll}
\underset{\eta,~\m{Y}}{\text{maximize}} & -\eta \\
\text{subject to} & \\
 & \begin{array}{rl}
\m{B}\m{Y} + \m{Y}\tr{\m{B}} &\preceq~~ 0 \\
-\m{Y} &\preceq~~ -\m{I} \\
\m{Y} - \eta \m{I} &\preceq~~ 0 \\
Y_{11} &=~~ 1, \quad \m{Y} \in \mathcal{S}^{n}
\end{array}
\end{array}
\]

The function \code{lmi1} takes as input a matrix $\m{B}$, and returns the input variables \code{blk}, \code{At}, \code{C}, and \code{b} for \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- lmi1(B)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

The second linear matrix inequality problem is

\[
\begin{array}{ll}
\underset{\m{P},~\ve{d}}{\text{maximize}} & -tr(\m{P}) \\
\text{subject to} & \\
& \begin{array}{rl}
\m{A}_{1}\m{P} + \m{P}\tr{\m{A}_{1}} + \m{B}*diag(\ve{d})*\tr{\m{B}} &\preceq~~ 0 \\
\m{A}_{2}\m{P} + \m{P}\tr{\m{A}_{2}} + \m{B}*diag(\ve{d})*\tr{\m{B}} &\preceq~~ 0 \\
-\ve{d} &\preceq~~ 0 \\
\sum_{i}^{p} d_{i}  &=~~ 1
\end{array}
\end{array}
\]

Here, the matrices $\m{B}$, $\m{A}_{1}$, and $\m{A}_{2}$ are known in advance. 

The function \code{lmi2} takes the matrices \code{A1}, \code{A2}, and \code{B} as input, and returns the input variables necessary for \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- lmi2(A1,A2,B)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

The final linear matrix inequality problem originates from a problem in control theory (\cite{boyd1994linear}) and requires three matrices be known in advance, $\m{A}$, $\m{B}$, and $\m{G}$

\[
\begin{array}{ll}
\underset{\eta,~\m{P}}{\text{maximize}} & \eta \\
\text{subject to} & \\
 & \left[\begin{array}{cc}
\m{A}\m{P} + \m{P}\tr{\m{A}} & \m{0} \\
\m{B}\m{P} & \m{0}
\end{array} \right]
+ \eta
\left[\begin{array}{cc}
\m{0} & \m{0} \\
\m{0} & \m{I}
\end{array} \right]
\preceq
\left[\begin{array}{cc}
-\m{G} & \m{0} \\
\m{0} & \m{0}
\end{array} \right]
\end{array}
\]

The function \code{lmi3} takes as input the matrices \code{A}, \code{B}, and \code{G}, and returns the input variables necessary to solve the problem using \code{sqlp}.

\begin{CodeChunk}
\begin{CodeInput}
R> out <- lmi3(A,B,G)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> sqlp(blk,At,C,b)
\end{CodeInput}
\end{CodeChunk}

\section{Concluding remarks}\label{sec:conc}

In Section \ref{sec:math}, 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 \ref{sec:gen3}.

In Section \ref{sec:usage}, we introduced the \proglang{R} package \pkg{sdpt3r}, and the main function call available in the package, \code{sqlp}. The specifics of the necessary input variables, the optional input variables, and the output variables provided by \code{sqlp} were presented. Using the examples from Section \ref{sec:math}, we showed how a problem written as a semidefinite quadratic linear program could be solved in \proglang{R} using \pkg{sdpt3r}.

Finally, in Section \ref{sec:morex}, we presented a number of additional problems that can be solved using the \pkg{sdpt3r} package, and presented the helper functions available so these problems could be easily solved using \code{sqlp}. 

The \pkg{sdpt3r} package significantly broadens the range of problems that can be solved using \proglang{R}. Here, we discussed a number of problems that can be solved using \pkg{sdpt3r}, including important problems in the statistical sciences, graph theory, classification, control theory, and general matrix theory. The \code{sqlp} function in \pkg{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 \ref{eq:Primal} or \ref{eq:Dual} by specifying the input variables \code{blk}, \code{At}, \code{C}, and \code{b} for their particular problem. 

\newpage

\appendix

\section{OPTIONS} \label{sec:A}

\begin{table}[hbtp]
\centering
\begin{tabular}{rl}
\code{vers} & specifies the search direction \\
 & 0, HKM if semidefinite blocks present, NT otherwise (default) \\
 & 1, HKM direction \\
 & 2, NT direction \\
\code{predcorr} & TRUE, use Mehrotra prediction-correction (default) \\
 & FALSE, otherwise \\
\code{gam} & step-length (default 0)\\
\code{expon} & exponent used to decrease sigma (default 1) \\
\code{gaptol} & tolerance for duality gap as a fraction of the objective function (default $1e-8$) \\
\code{inftol} & tolerance for stopping due to infeasibility (default 1e-8)\\
\code{steptol} & tolerance for stopping due to small steps (default 1e-6)\\
\code{maxit} & maximum number of iterations (default 100) \\
\code{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 \\
\code{scale_data} & TRUE, scale data prior to solving \\
 & FALSE, otherwise (default) \\
\code{rmdepconstr} & TRUE, remove nearly dependent constraints\\
 & FALSE, otherwise (default) \\
 \code{parbarrier} & declare the existence of a log barrier term \\
  & default value is 0 (i.e. no log barrier)
\end{tabular}
\end{table}

\bibliography{sdpt3r}
\end{document}