ROI: An Extensible R Optimization Infrastructure

Optimization plays an important role in many methods routinely used in statistics, machine learning and data science. Often, implementations of these methods rely on highly specialized optimization algorithms, designed to be only applicable within a specific application. However, in many instances recent advances, in particular in the field of convex optimization, make it possible to conveniently and straightforwardly use modern solvers instead with the advantage of enabling broader usage scenarios and thus promoting reusability. This paper introduces the R optimization infrastructure ROI which provides an extensible infrastructure to model linear, quadratic, conic and general nonlinear optimization problems in a consistent way. Furthermore, the infrastructure administers many different solvers, reformulations, problem collections and functions to read and write optimization problems in various formats.


Introduction
Optimization is at the core of inference in modern statistics since solving statistical inference problems goes hand in hand with solving optimization problems (OPs). As such statisticians, data scientists, and others who regularly employ computational methods ranging from various types of regression (e.g., constrained least squares, regularized least squares, nonlinear least squares), and classification (e.g., support vector machines, convex clustering) to covariance estimation and low rank approximations (e.g., multidimensional scaling, non-negative matrix factorization) benefit from advances in optimization, in particular in mixed integer and convex optimization. For example, Bertsimas, King, and Mazumder (2016) show that, thanks to a striking speedup factor of 450 billion in mixed integer optimization in the period of 1991-2015, the NP-hard best subset problem (Miller 2002) can now be solved reasonably fast (number of observations in the 100s and number of variables in the 1000s is solved within minutes). O'Donoghue, Chu, Parikh, and Boyd (2016) introduce the SCS solver for convex optimization problems, which can be used to solve among others (logistic) regression with l {1,2} regularization, support vector machines, convex clustering, non-negative matrix factorization and graphical lasso.
For R (R Core Team 2020), being a general-purpose tool for scientific computing and data science, optimization and access to highly efficient solvers play an important role. The field of optimization already has many resources to offer, like software for modeling, solving and randomly generating optimization problems, as well as optimization problem collections used to benchmark optimization solvers. In order to exploit the available resources more conveniently, over the years many modeling tools have emerged. One of the first systems used to model linear optimization problems is the so-called mathematical programming system (MPS) format (see Kallrath 2004). Developed in the 1960's, the MPS format today seems rather archaic but it is still widely used to store and exchange linear problems and is supported by most of the linear optimization solvers. Later, algebraic modeling languages (AMLs; e.g., GAMS, Bisschop andMeeraus 1982, andAMPL, Fourer, Gay, andKernighan 1989) became available. AMLs are domain specific languages (DSLs) dedicated to optimization. Today modern optimization systems are typically implemented in high-level programming languages like Julia (Bezanson, Edelman, Karpinski, and Shah 2017), MATLAB (The MathWorks Inc. 2019), Python (Python Software Foundation 2017) or R. Among the modern optimization systems, many are DSLs specially suited for convex optimization, such as YALMIP (Löfberg 2004) and CVX (Grant and Boyd 2014) in MATLAB, CVXPY (Diamond and Boyd 2016) and CVXOPT (Andersen, Dahl, and Vandenberghe 2016) in Python, Convex.jl (Udell, Mohan, Zeng, Hong, Diamond, and Boyd 2014) in Julia and CVXR (Fu, Narasimhan, and Boyd 2020) in R. JuMP (Lubin and Dunning 2015) is a DSL implemented in Julia designed for mixed-integer programming. pyOpt (Perez, Jansen, and Martins 2012) is a Python package for nonlinear constrained optimization.
Despite R having access to many modern optimization solvers which are capable of solving a wide class of optimization problems (see, e.g., the CRAN Optimization and Mathematical Programming Task View by Theußl, Borchers, and Schwendinger 2020), it is still commonplace to develop highly sophisticated special purpose code (SPC) for many statistical problems. The reasons are many. To name but a few: 1) availability, i.e., many solvers have not been easily available in R, 2) capability, i.e., problems could not be solved due to a lack of adequate solvers, and 3) efficiency, i.e., SPC tends to be faster. This paper introduces an extensible object oriented R optimization infrastructure (ROI) promoting the usage of optimization in R and R as a tool for optimization. In doing so it strives to enable users to formulate problems and experiment with different solvers in a straightforward way, help researchers to find the appropriate solver for their particular problem, or assist package developers to streamline their package dependencies. The framework is composed of package ROI (Hornik, Meyer, Schwendinger, and Theußl 2020) and its (at the time of this writing) 23 companion packages. Package ROI is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=ROI.
In contrast to DSLs, the ROI package does not aim to create a new language but provides a modeling mechanism borrowing its strength from the rich language features R has to offer. Optimization problems are constructed in a consistent way and stored in a single object. This makes it possible that problems are easily altered (reused) and shared before they are passed to a unified solve function. Such problems are then formulated and manipulated by using the provided R functions instead of special syntax from DSLs for which highly specialized knowledge would be required. Moreover, we believe that this approach makes it more attractive to add new solvers to the R solver landscape, e.g., to take advantage of recent advances in conic optimization (increase availability). Another key feature of ROI is that it is designed to be extensible. Companion packages equip ROI with state of the art optimization solvers, benchmark collections and functions to read and write optimization problems in various formats (increase capability). Furthermore, allowing package developers to plug-in new solvers quite effortlessly not only makes it easy to use their highly efficient code for a given problem but possibly also in many other applications (eliminate efficiency detriments). Currently ROI can be used to model and solve linear, quadratic, second order cone, semidefinite, exponential cone, power cone and general nonlinear optimization problems as well as mixed integer problems. This covers many optimization problems encountered in statistics, machine learning and data science (see, e.g., Koenker and Mizera 2014, for a survey of convex problems in statistics).
The remainder of this paper is organized as follows: In Section 2 we discuss the basic optimization problem classes, with a special focus on the newer developments in convex optimization. A survey of available R packages concerned with solving these problem classes is given in Section 3. Sections 4 and 5 show, respectively, how to formulate and solve optimization problems with the ROI package. Based on the tools presented in the previous sections, Section 6 provides basic examples. Section 7 is dedicated to the extension of ROI. Applications in the field of statistics are presented in Section 8. Section 9 concludes this paper.

Problem classes
Optimization is the process of allocating scarce resources to a feasible set of alternative solutions in order to minimize (or maximize) the overall outcome. Given a function f 0 : R n → R and a set C ⊆ R n we are interested in finding an x * ∈ R n that solves minimize f 0 (x) subject to x ∈ C. (1) The function f 0 is called the objective function. A point x is said to be feasible if it satisfies every constraint given by the set C of all feasible points defining the feasible region. If C is empty, then we say that the optimization problem is infeasible. Since maximization problems can be expressed as minimization problems by just changing the sign in the objective function, we subsequently will mainly deal with minimization problems.
An OP is called bounded if there exists a finite constant u such that f 0 (x) > u, ∀x ∈ C, if such constant does not exist, the problem is unbounded. Thus, a problem like in Equation 1 may or may not have a solution. Given a feasible and bounded problem a vector x * ∈ C that satisfies is commonly referred to as a solution of the OP.
Since any feasible set C can be expressed by the combination of constraint functions, the OP from Equation 1 can be written as: where b ∈ R m is the so-called right-hand-side (Nocedal and Wright 2006). The constraints f i , i = 1, . . . , m are sometimes referred to as functional constraints (Ben-Tal and Nemirovski 2001;Nesterov 2004). Since any equality constraint can be expressed by two inequality constraints and vice versa any inequality constraint can be expressed as an equality constraint by adding additional variables (also called slack variables), it is common practice to define OPs only in terms of either equality, less than or equal or greater than or equal constraints, to avoid redundancies.
Equation 2 is also sometimes referred to as the primal problem, which highlights the fact that there exists an alternative problem formulation, the dual problem. The dual problem is typically defined via the Lagrangian function (Lagrange duality; Nocedal and Wright 2006).
Several interconnected characteristics exist which determine how efficiently a given OP can be solved, namely convexity, the functional form of the objective, the functional form of the constraints and if the variable x is binary, integer, or continuous. An OP as displayed in Equation 1 is convex, if f 0 is convex and the set C is convex. Whereas modern solvers can efficiently solve a wide range of convex OPs and verify that a global solution (i.e., one as good or better than all other feasible solutions) was obtained, the same is mostly not true for non-convex problems (several local optima may exist). More information about convex programming can be found in, e.g., Boyd and Vandenberghe (2004); Ben-Tal and Nemirovski (2019).
Based on the functional form of the objective function and of the constraints, OPs can be divided into linear and nonlinear OPs. Furthermore, the class of nonlinear OPs can be further subdivided into conic, quadratic and general nonlinear OPs. In the following we give a formal definition of the different classes of OPs and overview their properties.

Linear programming
A linear program (LP) is an OP where all f i (i = 0, . . . , m) in Equation 2 are linear. Thus an LP can be defined as: where x is the vector of objective variables which has to be optimized. The coefficients of the objective function are represented by a 0 ∈ R n . A ∈ R m×n is a matrix of coefficients representing the constraints of the LP. Hence, in accordance with Equation 2, Ax ≤ b could also be written as a i x ≤ b i , i = 1, . . . , m (here a i refers to the ith row of the coefficient matrix A). All LPs are convex and usually solved via interior-point or simplex methods. For more information about the origination and mathematical properties of these methods we refer the reader to the book of Nocedal and Wright (2006).
A typical statistical problem which falls into this problem class is solving the least absolute deviations (LAD) or L 1 regression problem. Following, Wagner (1959) the objective function (4)

Quadratic programming
A quadratic program (QP) is a generalization of the standard LP shown in Equation 3, where the objective function contains a quadratic part in addition to the linear term. The quadratic part is typically represented by a matrix Q 0 ∈ R n×n . Therefore QPs can be expressed in the following manner: Unlike LPs, not all QPs are convex. A QP is convex if and only if Q 0 is positive semidefinite. A generalization of the QP is the quadratically constrained quadratic program (QCQP): A QCQP is convex if and only if all Q i (i = 0, . . . , m) are positive semidefinite (Lobo, Vandenberghe, Boyd, and Lebret 1998). Whereas convex QPs or even QCQPs are commonly solved by reformulations (transformations) to second-order cone programming (SOCP) or semidefinite programming (SDP; see Section 2.3), the question how to obtain a reliable global solution for a non-convex QCQP is still an active field of research. Details on the necessary transformations to cast a convex QCQP into an SOCP or SDP can be found in, e.g., Lobo et al. (1998); Alizadeh and Goldfarb (2003); Bao, Sahinidis, and Tawarmalani (2011).

Conic programming
Conic programming refers to a class of problems designed to model convex OPs. The most prominent members of this class are LP, SOCP and SDP. We follow the common practice to define a conic program (CP) as: where the set K is a nonempty closed convex cone.
The standard form of a CP as given in Equation 7 minimizes a linear objective over a convex cone (b − Ax = s ∈ K). As Nemirovski (2006) points out, representing CPs in this form has two main advantages. First, this formulation has strong unifying abilities which means only a few cones allow modeling of many different types of OPs. Additionally, the nonlinearities are no longer represented by general nonlinear objective and constraint functions but vectors and matrices which allows the algorithms to utilize the structure present in the convex OPs. Second, the convexity is built-in into the definition of CPs. At the same time, theoretically, any convex OP can be reformulated into the form given in Equation 7. Thereby nonlinear objective functions are expressed in epigraph form (see, e.g., Boyd and Vandenberghe 2004): Practically the number of OPs which can be solved via CP is limited by the number of cones supported by a given optimization solver. State of the art solvers distinguish between up to eight different types of cones. Following the definitions in Diamond and Boyd (2015) and O'Donoghue et al. (2016), a convex cone K is typically a Cartesian product from simple convex cones of the following types.

Zero cone and free cone
The zero and free cones are, respectively, given by where for a cone K we write K * = {y|x y ≥ 0 for all x ∈ K} for the dual cone (see, e.g., Boyd and Vandenberghe 2004 for more information about dual cones). From Equation 7 it can be immediately seen, that in the case of linear equality constraints s i has to be zero, i.e.,

Linear cone (non-negative orthant)
The linear cone is given by This cone is used to represent linear inequality (less than or equal) constraints, by requiring s i to be non-negative, i.e., From the definition of the free cone and non-negative cone, it is apparent that any LP can be written as a CP where K is a product of free and non-negative cones.

Second-order cone
The second-order cone is given by This cone is used to model sums of norms as well as convex QPs and QCQPs (Lobo et al. 1998;Alizadeh and Goldfarb 2003). CPs where K is a product of free, non-negative and second-order cones are commonly referred to as SOCPs.

Positive semidefinite cone
The positive semidefinite (PSD) cone is given by Here S n refers to the space of real-symmetric n × n matrices. CPs restricted to the positive semidefinite cone are referred to as SDPs. They are commonly used for solving combinatorial problems (e.g., maximum cut problem) and for solving convex QPs and QCQPs (Vandenberghe and Boyd 1996;Helmberg 2000;Freund 2009;Bao et al. 2011). Lobo et al. (1998) show that each SOCP can be rewritten as an SDP.

Exponential cone
The primal exponential cone is defined as Its dual is given by As can be inferred from Equation 12, the exponential cone can be used to model exponential functions and logarithms. More details about the exponential cone and functions representable by the exponential cone can be found in Chares (2009) andSerrano (2015).

Power cone
The 3-dimensional primal power cone has already been investigated in Koecher (1957) and is defined as Its dual is given by The power cone can be used to model powers and p-norms. For more information about the power cone and its modeling capabilities we refer to Chares (2009).
Putting the hierarchies described above all together we get the following ordering among OPs LP ⊂ convex QP ⊂ convex QCQP ⊂ SOCP ⊂ SDP ⊂ CP.

Nonlinear optimization
The most general problem class is nonlinear optimization or nonlinear programming (NLP). This is the problem where at least one f i , i = 0, . . . , m in Equation 2 is not linear. NLPs are not required to be convex, which makes it in general hard to obtain a reliable global solution.
Contrary to the convex case, in a non-convex setting most optimization algorithms only find the extremum of f 0 in the neighborhood of the starting value (local optimum).

Mixed integer programming
A mixed integer program (MIP) adds the additional requirement to the optimization problem that some of the objective variables can only take integer values. Considering Equation 2, a problem is called a mixed integer problem if the (type) constraint x k ∈ Z for at least one k is added. In the case where all n objective variables are integral we speak of a pure integer programming (IP) problem. An IP where all variables are bounded between zero and one, i.e., x ∈ {0, 1} n , is called a binary (integer) program.
Since MIPs are non-convex, even mixed integer linear programs (MILP) can already be hard to solve. Nevertheless an increase in quantity and quality of free and nonfree solvers was observed in the last decades (Linderoth and Ralphs 2005;Bixby 2012). Typically solvers use branch-and-bound (Land and Doig 1960) and the cutting plane (Gomory 1960) algorithms or a combination of both. The algorithms avoid solving the problem directly, but instead solve multiple relaxations where the integer constraint is dropped.

Software
Recently, an increase of the available packages handling many different OPs in R has been observed. The CRAN task view Optimization and Mathematical Programming  currently lists around 100 different optimization related packages. The capability these packages provide range from solvers which can solve a wide range of optimization problems (e.g., optimx, Nash and Varadhan 2011; Nash 2014a) to very specialized solvers which are created to solve a specific problem type very fast (e.g., nonlinear regression solvers). This section provides an overview of the solver landscape in R. The insights gained in this section will be used to derive a consistent solver infrastructure. First, we investigate the available (open source) solvers, splitting these into linear solvers, quadratic solvers, conic solvers and general purpose solvers. We then discuss commercial solvers (i.e., any solver developed for sale) and the NEOS server.

Overview
As pointed out in Section 2, in the field of optimization we are typically facing different problem classes. The possibly three most important distinctions are between linear versus nonlinear problems, integer versus continuous and convex versus non-convex problems.
Ordered based on increasing complexity, an objective function might be of type linear, (convex) quadratic, conic (i.e., any objective expressible as a CP) or functional (i.e., any objective expressible as a function). Similarly constraints are typically of type box, linear, (convex) quadratic, conic or functional. Box constraints (or variable bounds) are a special type of linear constraints which enforce lower and upper bounds on the objective variables.
The terms conic objective/constraints are used in a general way and refer to any linear and nonlinear objective/constraints that can be reformulated as a conic problem. Therefore this also includes problems with linear and convex quadratic objective/constraints. Note that, solvers that take as input values a linear objective and conic constraints are also applicable to OPs with conic objective and conic constraints by making use of the epigraph form transformation. The most general form are functional objective/constraints which includes all linear and nonlinear objective/constraints. Table 1 gives an overview on optimization packages available on CRAN (https://CRAN. R-project.org) with a focus on open source solvers. The position of a particular package in the table indicates its ability to solve a given problem. Each problem class to the left and above of the current position can be handled by the package including its current position. For instance, the ECOSolveR package (Fu and Narasimhan 2019) which provides an interface to the ECOS library (Domahidi, Chu, and Boyd 2013) can solve conic problems restricted Table 1: Selected R packages displayed based on the types of optimization problems they are applicable to. Here * indicates that the solver is restricted to convex problems and + indicates that the solver can model integer constraints.
to combinations of the zero, non-negative, second-order and primal exponential cone. Since ECOS is equipped with a branch-and-bound algorithm, it can also be used to solve mixed integer conic problems.

The R solver landscape
The solver landscape can be split into two parts. First, solvers where the functional form is fixed and only the coefficients are provided, which includes all LP, QP, QCQP and CP solvers currently available in R. Second, solvers which can optimize any functional form expressible as an R function. This includes most NLP solvers, sometimes summarized as general purpose solvers.

Linear solvers
Interfaces to several open source LP and MILP solvers are available in R. Most of these packages provide a high-level access to the solver, those explicitly designed to provide a lowlevel access are commonly marked with the suffix API. lp_solve (Berkelaar, Eikland, and Notebaert 2016) uses the simplex algorithm combined with branch-and-bound to solve LPs and MILPs. It furthermore allows modeling of semicontinuous and special ordered sets problems. Packages lpSolve (Berkelaar 2020) and lp-SolveAPI (Konis and Schwendinger 2020) provide access to the lp_solve solver in R.
Additionally the function lpcdd() from package rcdd (Geyer, Meeden, and Fukuda 2019) and the function simplex() from package boot (Canty and Ripley 2020) can be used to solve LPs via the simplex algorithm.
By taking a closer look at the elements needed by packages capable of solving LPs and MILPs 1 we can conclude that the following elements should be present in a consistent and convenient optimization infrastructure for modeling LPs and MILPs.
objective: A numeric vector giving the coefficients of the linear objective. constraints: • Includes a constraint matrix A (see Equation 3), • a vector giving the direction of the constraints (i.e., ==, <= or >=), and • a vector giving the right-hand-side b (see Equation 3).
bounds: Two vectors giving the lower and upper bounds.
types: A vector storing the type information, i.e., binary, integer and numeric.
maximum: A Boolean indicating if the objective function should be maximized or minimized.
Note that the elements bounds and maximum, as well as the constraint directions and the binary types are not strictly necessary. Their inclusion is motivated by the fact that they are supported by many solvers and simplify the problem specification.

Quadratic solvers
As Table 1 shows, most of the quadratic solvers are designed to solve convex quadratic problems with linear constraints. The quadprog (Turlach and Weingessel 2019) package uses the dual method described in Goldfarb and Idnani (1983). LowRankQP (Ormerod and Wand 2020) is based on an interior-point algorithm described in Fine and Scheinberg (2001). Dykstra (Helwig 2018) implements Dykstra's cyclic projection algorithm (Dykstra 1983), coneproj (Meyer and Liao 2018) transforms the original QP problem into a cone projection problem (Liao and Meyer 2014) and osqp (Stellato, Banjac, Goulart, and Boyd 2019) uses the alternating direction method of multipliers described in Stellato, Banjac, Goulart, Bemporad, and Boyd (2017).
Additionally, the package ROI.plugin.qpoases (Schwendinger 2020) and the ipop() function from kernlab (Karatzoglou, Smola, Hornik, and Zeileis 2004;Karatzoglou, Smola, and Hornik 2019) can be used to obtain solutions for non-convex quadratic problems with linear constraints. However, for the non-convex case there is no guarantee that the returned solution is a global optimum. ROI.plugin.qpoases is an interface to the qpOASES (Ferreau, Kirches, Potschka, Bock, and Diehl 2014;Ferreau, Potschka, and Kirches 2018) library, which uses an online active set strategy to solve quadratic optimization problems.
QP solvers generally take the same arguments as LP solvers plus an additional matrix parameter storing the coefficients of the quadratic term Q 0 .

Conic solvers
Most of the conic solvers use a standard form similar to Equation 7, where the objective function is assumed to be linear and the vector b − Ax is restricted to a certain cone K. Nevertheless, in Table 1 they are shown to have a conic objective function and conic constraints to express that they are able to solve any LP and convex NLP expressible by a CP. Therefore, which types of NLPs a given solver can solve, depends on the types of cones the solver can model. Table 2 shows the conic solvers available in R and the types of cones they support.
Package CLSOCP (Rudy 2011) is specialized in solving SOCPs, it is a pure R implementation of the one-step smoothing Newton method based on the algorithm described in Tang, He, Dong, and Fang (2012

General purpose solvers
Solvers capable of handling nonlinear objective functions without further restrictions are called general purpose solvers (GPSs). These solvers can minimize (or maximize) any functional form representable as an R function with -depending on solver capabilities -different types of constraints, where again the most general form of constraint is the functional constraint (i.e., an R function). The generality of GPSs comes at the price of performance and that there is usually no guarantee that a global optimum is reached.
Important properties of GPSs are whether they are designed to search for a local or global optimum, if gradient information has to be provided or the method is gradient free and which type of constraints can be set. Table 3 shows the number of GPS methods grouped by these properties (the counts are based on Table 5 where additional details can be found) and reveals some interesting details about the R GPS landscape. Around 60 percent of the GPS are designed for local optimization, even though most of the local solvers utilize gradient information, only four of the global solvers use gradient information. The difference in distribution of gradient based and gradient free optimization algorithms between global and local GPSs can be explained by the fact that in global optimization, metaheuristics like evolutionary methods or particle swarm optimization are commonly used. In a recent study Mullen (2014a) surveys the continuous global optimization packages available in R and compares their performance on a set of tests bundled in the globalOptTests (Mullen 2014b) package. Table 5 gives an extensive listing of which methods are designed to search for a global solution. For more information about the methods we refer to Mullen (2014a).
Based on the type of constraints, the GPSs can be divided into no constraints, box constraints, linear constraints, quadratic constraints and functional constraints. As Table 3 shows, most of the GPSs support no constraint or box constraints. Fortunately, package optimx (Nash and Varadhan 2011) provides a unified interface to many of these solvers, consolidating methods from packages stats, ucminf (Nielsen and Mortensen 2016), minqa (Bates, Mullen, Nash, and Varadhan 2014), Rcgmin (Nash 2014b), Rvmmin (Nash 2018) and BB (Varadhan and Gilbert 2009). It was designed as a possible successor of optim which is part of the stats package and can be used to solve OPs with box constraints. Another package which incorporates many different algorithms is nloptr (Ypma and Johnson 2020). It is an R interface to the NLopt (Johnson 2016) library, which bundles several global and local optimization algorithms. Depending on the algorithm it can solve NLPs with box constraints or functional constraints.
Most of the GPSs able to handle functional constraints allow to specify functional equality and/or functional inequality constraints.
To model functional equality constraints the following two forms are most commonly used  Ghalanos and Theußl 2015), where h is a function and b ∈ R k gives the right-hand-side. Similarly, functional inequality constraints are commonly given in one of the following three forms: • g j (x) ≤ 0, j = k + 1, . . . , m (e.g., DEoptimR, nloptr::nloptr, NlcOptim, Chen and Yin 2019; Rnlminb2, csr::snomadr, Racine and Nie 2019), • g j (x) ≥ 0, j = k+1, . . . , m (e.g., alabama, nloptr::auglag, nloptr::cobyla, nloptr:: ires, nloptr::mma, neldermead, Bihorel and Baudin 2018; and nloptr::slsqp), where g is a function and l ∈ R m−k , u ∈ R m−k are the lower and upper bounds of the constraints. A general optimization infrastructure should be designed in a way that the functional form employed can be transformed into the commonly used forms shown above.
An analysis of the above solver spectrum reveals that the critical arguments to GPSs are:

start:
The initial values for the (numeric) parameter vector.
objective: The function to be optimized.
constraints: Depending on the GPS the constraints can be none, linear, quadratic, functional equality or functional inequality constraints. To model functional constraints consistently with linear and quadratic constraints the following elements are needed: • a function representing the constraints; • a vector giving the direction of the constraints; and • a vector giving the right-hand-side.
bounds: Variable bounds, commonly given as lower and upper bounds.
Additionally some GPSs make use of the gradient and/or Hessian of the objective and the Jacobian of the constraints. The optional elements can be summarized by: gradient: A function that evaluates the gradient of the argument objective.
hessian: A function that evaluates the Hessian of the argument objective.
jacobian: A function that evaluates the Jacobian of the argument constraints.
maximum: A Boolean indicating whether the objective function should be maximized or minimized.
control: Further control arguments specific to the solver.
value/objective: The value of the objective function evaluated at the "solution".
convergence, status: An integer information about the convergence and exit status of the optimization task.

gradient:
The gradient evaluated at the solution found.
hessian: The Hessian evaluated at the solution found.
message: A text message giving additional information about the optimization / exit status.
iterations/evaluations: The number of iterations and / or evaluations.

Commercial solvers
Since commercial solver packages often bundle a variety of solvers, it is often not possible to assign them to a certain problem class. At the time of this writing R interfaces are available to the commercial solver software CPLEX (IBM ILOG 2019), MOSEK (MOSEK ApS 2017), Gurobi (Gurobi Optimization, Inc. 2016), Lindo (Lindo Systems 2003) and localsolver (Benoist, Estellon, Gardi, Megel, and Nouioua 2011). To cover also the commercial side, the interface packages Rcplex (Theußl and Bravo 2016) and Rmosek (MOSEK ApS 2019) were included in defining the requirements for a consistent solver infrastructure.

Network-enabled optimization system (NEOS)
The NEOS (Czyzyk, Mesnier, and Moré 1998;Dolan 2001;Gropp and Moré 1997) server (https://neos-server.org) provides free access to more than 60 numerical optimization solvers. It can be accessed via the internet by submitting OPs via the homepage, email, the XML-RPC application programming interface or Kestrel. Depending on the solver the OPs have to be formulated in different ways, overall most solvers support input from AMPL and GAMS. The R packages rneos (Pfaff 2020) and ROI.plugin.neos (Hochreiter and Schwendinger 2020) use the XML-RPC API to communicate with the NEOS server. In order to upload OPs with rneos, the OPs have to be formulated in the input format supported by the solver (e.g., AMPL, GAMS, MPS). In contrast to that, ROI.plugin.neos supports OPs generated with ROI, internally each OP is transformed to GAMS before they are submitted to the server. The result is again converted back into a suitable solution format.

Data structures
After reviewing the optimization resources available in R, it is apparent that the main function of a general optimization infrastructure package should take at least three arguments: problem representing an object containing the description of the corresponding OP, solver specifying the solver to be used (e.g., "glpk", "nlminb", "scs"), control containing a list of additional control arguments to the corresponding solver.
The arguments solver and control are easily understood, since from the available solver spectrum we only have to choose those which are capable to handle the corresponding OP and (optionally) supply appropriate control parameters. However, building the object for the problem argument, in a general and intuitive way, is a challenging task which leads to several design issues.
Based on the review in Sections 2 and 3 it seems natural to instantiate OPs based on an objective function, one or several constraints, types and bounds of the objective variables, as well as the direction of optimization (whether a minimum or a maximum is sought).
In the remainder of this section we discuss the conceptual design of ROI and how to use it to formulate optimization problems. For illustrative purposes we already use functionality of package ROI.

Objective function
The survey of optimization solvers in Section 3 reveals that the way the objective function is stored depends primarily on its functional form. If the objective function is linear (L), i.e., a 0 x, then it is common practice to only supply a coefficient vector a 0 ∈ R n . For quadratic objective functions (Q) of the form 1 2 x Q 0 x + a 0 x most solvers take a vector a 0 ∈ R n and a matrix Q 0 ∈ R n×n as input. General nonlinear objective functions (F, i.e., nonlinear functions which cannot be represented as an QP or CP), are represented as an R function which takes the vector of objective variables as argument and returns the objective value. Depending on the type of the objective function, i.e., F, Q, or L, only a subset of the solver spectrum can be used.
Objective function types and corresponding constructors implemented in ROI are: F The most general form of an objective function is created with the F_objective(F, n, G, H, names) constructor by simply supplying F, an R function representing f 0 (x), and n the length of x. Optionally, information about the gradient and the Hessian can be provided via the arguments G and H. If no gradient is provided it will be calculated numerically if needed. The optional names argument is propagated to the solution object to make the solution more readable.
Q Objective functions representing a quadratic form as outlined above can be easily created with the Q_objective(Q, L, names) constructor taking Q, the quadratic part Q 0 , and optionally L, the linear part a 0 , as arguments. The names argument is again optional.
L If the objective to be optimized is a linear function then one should use the L_objective(L, names) constructor supplying L (the coefficients of the objective variables) as a numeric vector. The names argument is again optional.
All three constructors return an object inheriting from class 'objective'.

Constraints
To model all the problem classes introduced in Section 2, four different types of constraints are sufficient. Thereby arguments with the same name have the same functionality irrespective of the constraint type and will therefore be explained only once.
F The most general form of constraints can express any constraint representable by an R function. They are created via F_constraint(F, dir, rhs, J, names). Here F is either a function or a list of functions, dir is a character vector giving the direction of the constraint and rhs is a numeric vector giving the right-hand-side of the constraint. The optional arguments J and names can be used to provide the Jacobian and the variable names of the constraints.
C Conic constraints are constructed via the function C_constraint(L, cones, rhs, names), where L can be either a numeric vector of length n or a matrix of dimension m × n.
In accordance with Equation 7 the cones argument imposes a restriction on the slack variables s. A conic constraint can be comprised of several cones, where each cone type can occur multiple times. The cone constructors all start with K_ followed by a short cut of the cone name, as defined in Section 2.3. Currently ROI implements constructors for the cones K_zero, K_lin, K_soc, K_psd, K_expp, K_expd, K_powp, and K_powd. To combine different cones the generic combine c() can be used.
Q Quadratic constraints as defined in Equation 6 can be easily created with the constructor Q_constraint(Q, L, dir, rhs, names). The quadratic constraints Q are given as a list of length m where the entries are either n × n matrices or NULL.
L Linear constraints are constructed via the function L_constraint(L, dir, rhs, names).
All constructors return an object inheriting from class 'constraint'. Since in many situations it is desirable to optimize a given objective function subject to composite constraints of different kinds, ROI can combine multiple constraints into a single constraint using the generic functions c() or rbind(). Since the matrices L and Q can become huge but are typically sparse we use the simple triplet matrix format from the slam (Hornik, Meyer, and Buchta 2019) package to store them internally. In the simple triplet matrix format (sometimes referred to as coordinate list format) only the row indices, the column indices and the values of the non-zero elements are stored in a list. We choose the slam package not only for efficiency reasons but also due to the fact that many solver APIs demand such a format.

Objective variable types
As it is common practice in mixed-integer solvers to distinguish between the variable types continuous, integer and binary we encode the variable choice with the following characters: "C" for continuous, "I" for integer and "B" for binary. By default all the variables are assumed to be of continuous type.

Bounds
A variable bound is a special type of constraint used to restrict an objective variable between a real lower and upper bound. Therefore, variable bounds are often also called "box bounds" or "box constraints". Although variable bounds could be easily modeled as linear constraints, many GPSs only support variable bounds (see Table 3). Furthermore, most solvers that support any type of constraint allow to specify variable bounds directly. Thus, it is reasonable but also convenient to consider them separately.

ROI Optimization Problem:
Maximize a linear objective function of length 3 with -3 continuous objective variables, subject to -3 constraints of type linear.
-1 lower and 1 upper non-standard variable bound.
The setter functions make it easy to alter previously created OPs. The getter function objective() returns the objective as function, which can be directly used to evaluate parameters. The number of parameters required can be obtained by the generic function length().

$names NULL
For all the other elements the corresponding getter function returns directly the underlying data representation.
Function OP() always returns an S3 object of class 'OP' which stores the entire OP. Storing the OP in a single R object has many advantages, among others: • The OP can be checked for consistency during the creation of the problem.
• The different elements of the OP can be easily accessed after the creation of the problem.
• The OP can be easily altered, e.g., a constraint can be added, bounds can be changed, without the need to redefining the entire OP.
The consistency checks verify that the dimensions of the arguments fit together.

Functionality
The R optimization infrastructure is structured into the package ROI and its accompanying extensions (plug-ins and models). Package ROI provides all the necessary classes and methods and manages the extensions (i.e., automatically loads plug-ins and manages meta data about the plug-ins). The extension packages add optimization solvers, read/write functions and additional resources (e.g., model collections). The plug-in extensions play a special role, hence all plug-ins are loaded automatically when ROI is loaded. When a plug-in is loaded it provides data about its capabilities. This data is stored in an in-memory database and includes information about to which problems the plug-in is applicable, which formats it can read/write, the control arguments available for the solver and how the solver specific control arguments relate to arguments commonly used.
This mechanism makes it possible that ROI is aware of all the installed plug-ins, without the need to change ROI when a new plug-in is added. To make the automatic loading possible the plug-ins have to follow the name convention ROI.plugin.<name>, where <name> is typically the name of an optimization solver (e.g., ROI.plugin.glpk; Theußl 2017). The prefix ROI.models (e.g., ROI.models.netlib; Schwendinger 2019a) is used for data packages with predefined OPs. In Section 5.6 we give an overview on the data packages available in the ROI format.

Solving optimization problems
After formulating an OP as described in Section 4, it can be solved by calling the function ROI_solve(x, solver, control, ...). This function takes an R object of class 'OP' containing the formulation of the OP, the name of the solver to be used and a list containing solver specific parameters as arguments. The solver and control arguments are optional, if no solver argument is provided ROI will choose an applicable solver automatically (see Section 5.7). Alternatively the solver specific parameters can be specified via the dots arguments.

Status code
Solver status codes are used to inform the user about the exit status of the solver. Despite the common usage of status codes in optimization solvers there is no widely used standard. Nevertheless, we believe it is desirable to provide unified status codes. The status codes used in ROI_solve are simple and consistent with the common practice, to return 0 on success (if a "solution" meeting the solver specific requirements was found) and 1 otherwise. For optimization solvers specialized on solving convex LPs, QPs and CPs it can be assumed that a global solution was found, if the status code is 0. Unfortunately, the same is not true for non-convex QPs and GPS, here status codes are less informative. Some packages like optimx and alabama check the Karush-Kuhn-Tucker (KKT) conditions to verify that the solution found meets the criteria of a local minimum.

Solution
To make the solutions of the various solvers easy to understand, all the solutions are canonicalized. After the canonicalization each solution contains the following components: solution the solution of the OP, objval the optimal objective value, status the canonicalized status code, message the original solver message and a meta attribute containing the solver name and additional optional arguments.
To obtain the (primal) "solution" the generic function solution(x, type) should be used:

R> solution(lp_sol)
[1] 0.000000 9.238806 -1.835821 Some OPs have multiple solutions, in the case of BLP (MILP) some solvers can retrieve all (multiple) solutions. For MILPs it is in general not possible to obtain all the solutions but only multiple solutions, since even this simple MILP has an infinite number of solutions. The following example is based on Fischetti and Salvagnin (2010) and will be used later to illustrate how multiple solutions can be obtained.

R> solution(blp_sol)
[ alternatively it is also possible to set nsol_max to Inf. This is advantageous if an upper bound on the number of solutions is hard to guess and all the solutions should be retrieved.
If the status code is not zero, solution will return NA to prevent the user from using solutions with a status code different from 0.
R> lp_inaccurate_sol <-ROI_solve(lp, solver = "scs", tol = 1e-32) R> solution(lp_inaccurate_sol) [1] NA NA NA However in a few situations it can be desirable to obtain solutions even if the solver signals no success. In these cases ROI can be forced to return the solution provided by the solver regardless of the status code.

Reformulations
Reformulations are often used to transform a problem of class A into a problem of class B, where the solution of the original problem can be derived from the solution of the reformulation (which is typically easier to solve). Although reformulation techniques are commonly used in optimization, the functions performing these reformulations are generally hidden within the optimization software. To facilitate the comparison of different reformulation algorithms, ROI provides functions for managing reformulations. The available reformulations are listed by calling function ROI_registered_reformulations() and ROI_reformulate(x, to, method) performs the reformulation. Following Boros and Hammer (2002) we illustrate the transformation of a binary QP into a MILP. The code for the reformulation is based on the implementation in the relations (Meyer and Hornik 2019) package.
The binary QP is given by: In the following the ROI code is given to define the binary QP, transform the problem into a MILP and solve it: R> Q <-rbind(c(0, 3, 0), c(0, 0, 1), c(0, 0, 0)) R> bqp <-OP(Q_objective(Q = Q + t(Q), L = c(-1, -4, -1)), + types = rep("B", 3)) R> glpk_signature <-ROI_solver_signature("glpk") R> head(glpk_signature, 3) Optimal solution found. The objective value is: -4.000000e+00 Here ROI selects the applicable reformulations based on the provided signatures. A method is considered to be applicable if it can transform the given OP into a new OP, where the signature of the new OP is a subset of the signature provided in the argument to. Since it is possible that several methods are applicable, the argument method can be used to select a specific reformulation method.

ROI solvers
ROI can currently make use of more than thirty different solver methods, applicable to a wide range of OPs. Inspired by R's available.packages() function, ROI can return a listing of the solver plug-ins available at CRAN, R-Forge (https://R-Forge.R-project.org/, Theußl and Zeileis 2009) and GitHub (https://github.com/). ROI_available_solvers() without an argument lists all the available solvers. If an OP is provided as argument, only the available solvers which are also applicable will be returned.

com/FlorianSchwendinger
A listing of all the available plug-ins on CRAN and R-Forge could be easily compiled by just using the available.packages() function. But to be able to find all the solvers available and applicable to a given OP also the solver signature is needed. Therefore a database (an rds file on R-Forge) containing the solver signatures and the information provided by available.packages() was compiled and is queried whenever ROI_available_solvers() is called.
Both return values are based on the solver registry, which stores the solver method and information about the solver registered by the plug-ins. The solver registry is an in-memory database based on the registry (Meyer 2019) package.
ROI_installed_solvers() gives a listing of all the installed plug-ins (not necessarily loaded) delivered directly with ROI and found by searching for the prefix 'ROI.plugin' in the R library trees.

R> head(ROI_installed_solvers(), 3)
nlminb alabama cbc "ROI.plugin.nlminb" "ROI.plugin.alabama" "ROI.plugin.cbc" An overview on the currently available solver plug-ins based on the problem types is given in Table 4. Please note that the functionality provided in a plug-in does not necessarily have to be the same as the functionality of the solver, e.g., ROI.plugin.nlminb can take functional constraints, while nlminb can only take box constraints.
Furthermore we want to emphasize that ROI was built to be extended, as shown in Section 7.

ROI models
Test problem collections are commonly used in optimization to evaluate and compare the performance of solvers. As each class of optimization problems has its own test sets stored in various formats, ROI currently provides access to the NETLIB-LP and MIPLIB collections and the globalOptTests package. The NETLIB-LP collection (Gay 1985) is a collection of linear programming problems, which, even though the main part was created more than 30 years ago is still used today. Mixed integer optimization problems are commonly evaluated using the MIPLIB collection (Koch et al. 2011), an extensive collection of academic and industrial MILP applications. The globalOptTests package contains 50 box constrained nonlinear global OPs for benchmarking purposes.
The problems contained in these collections and packages were transformed into ROI optimization problems and can be accessed through the ROI.models.netlib, ROI.models.miplib (Schwendinger and Theußl 2019) and ROI.models.globalOptTests (Schwendinger 2017) packages. Since MIPLIB provides no license file, the OPs are not included in the package but can be easily obtained with the miplib_download_*() functions. Since the problems are stored as objects of class 'OP', they can be directly entered into ROI_solve.
R> library("ROI.models.netlib") R> agg <-netlib("agg") R> ROI_solve(agg, "glpk") Optimal solution found. The objective value is: -3.599177e+07 ROI makes these data collections (test problem sets) available in a common format, so users can easily compare the different solvers and developers interested in creating optimization software can use them to test their packages. Furthermore, the intuitive structure of ROI objects and its use of sparse data structures make it possible to directly derive a new format for the exchange of linear, quadratic and conic optimization problems.

R> library("jsonlite") R> nested_unclass <-function
x <-lapply(x, nested_unclass) + x + } R> agg_json <-toJSON(nested_unclass(agg)) R> tmp_file <-tempfile() R> writeLines(agg_json, tmp_file) The resulting text file can be easily imported into any programming language supporting JavaScript Object Notation (JSON). JSON is an open-standard file format that can be parsed by almost all programming languages. For historic reasons, OP collections are commonly provided in flat file formats (e.g., MPS, QPS). We believe, that today it would be advantageous to store them in general data exchange formats like JSON or XML (Extensible Markup Language).

ROI settings
Many general and/or solver-related settings can be set or modified using the ROI_options() function. Apart from that ROI recognizes some environment variables.

Gradient and Jacobian
When creating a plug-in, the function G should be used to derive the gradient and J should be used to derive the Jacobian. In ROI the gradient and Jacobian are derived analytically for linear and quadratic terms. For the derivation of nonlinear terms, by default, the numDeriv (Gilbert and Varadhan 2019) package with the Richardson extrapolation is used. However, this can be easily changed by providing customized functions to derive the gradient or Jacobian function.

Solver selection
If no solver is provided in ROI_solve, the default solver set in ROI_options will be used.

R> ROI_options("default_solver")
[1] "auto" By default the option "default_solver" is set to "auto" which enables automatic solver selection, if any other solver name (e.g., "glpk") is provided the automatic solver selection is discarded in favor of the specified solver.

Load plug-ins
The plug-ins are loaded automatically. However, in some situations it is desirable to deactivate the automatic loading and require plug-in packages one at a time. This can be accomplished by setting the environment variable ROI_LOAD_PLUGINS to FALSE.

R> Sys.setenv(ROI_LOAD_PLUGINS = FALSE)
Afterwards the default load behavior of ROI is altered and only the "nlminb" solver (which is included in ROI) gets registered when library("ROI") is called. Therefore, all the other plug-ins have to be loaded manually if needed (e.g., library("ROI.plugin.glpk")).

Examples
In this section we provide small examples to introduce the reader into the modeling capabilities of ROI.

Linear cone
The linear cone is used to model linear less than equal constraints, since Ax ≤ b is equivalent to Ax + s = b, s ∈ K lin . The linear constraint 5x 1 + x 2 ≤ 6 ⇐⇒ 6 − (5 1)x = s ∈ K lin can be expressed as follows: R> cpleq <-C_constraint(c(5, 1), K_lin(1), 6) By combining the zero cone constraint with the linear cone constraint the LP stated in Equation 20 can also be formulated as a CP.

R> solution(zlcp_sol)
[1] 1 1 Note that since in the definition of K zero and K lin only one variable is involved, a single zero cone or linear cone constraint can be expressed by a single row of the constraint matrix A. For all the other cones at least two rows of the constraint matrix A will be needed to express a single conic constraint.

Second-order cone
The second-order cone is used to express constraints of the type u 2 ≤ w (see Equation 10), where the variables u ∈ R n−1 and w ∈ R are expressed by b − Ax ∈ K n soc . Specifically, w is expressed by b 1 − a 1 x and u i is expressed by b i − a i x, i = 2, . . . , n.

R> solution(soc1_sol)
[1] 19.055671 6.300041 9.000000 Consider the SOCP minimize by making use of the epigraph form the OP can be rewritten into the standard form.

Power cone
The primal power cone is used to express constraints of the type u α v 1−α ≥ |w| (see Equation 14), here u ∈ R, v ∈ R, w ∈ R, u, v ≥ 0 and α ∈ [0, 1]. Since three scalar variables are involved, one primal exponential cone adds three rows to the constraint matrix A. The variables u, v and w are again expressed by the corresponding elements of b − Ax. Specifically Consider the CP minimize y 3y 1 + 5y 2 subject to 5 + y 1 ≥ (2 + y 2 ) 4 y 1 ≥ 0, y 2 ≥ 2.

Positive semidefinite cone
The positive semidefinite cone is used to express constraints of the type X ∈ S n and z Xz ≥ 0 for all z ∈ R n (see Equation 11). Positive semidefinite constraints are expressed by the constraint matrix A and right-hand-side b. To express the linear matrix inequality n i=1 in terms of b − Ax ∈ K psd the symmetric matrices F i ∈ R d×d are transformed into vectors by a half-vectorization. Half-vectorization is a special kind of matrix vectorization for symmetric matrices, which transforms a symmetric matrix R> (A <-matrix(c(1, 2, 3, 2, 4, 5, 3, 5, 6), nrow = 3)) subject to x 1 10 3 3 10 + x 2 6 −4 −4 10 + x 3 8 1 1 6 16 −13 −13 60 x 1 , x 2 , x 3 ≥ 0.

General nonlinear optimization problems
The following example from Rosenbrock (1960) is known as Rosenbrock's post office problem.
The following code can be used to solve this problem with ROI.

ROI Optimization Problem:
Maximize a nonlinear objective function of length 3 with -3 continuous objective variables, subject to -1 constraint of type nonlinear. -0 lower and 3 upper non-standard variable bounds.

ROI Optimization Problem:
Maximize a nonlinear objective function of length 3 with -3 continuous objective variables, subject to -1 constraint of type linear. -0 lower and 3 upper non-standard variable bounds.
using L and Q constraints rather than F_constraint has the advantage that for L and Q constraints the Jacobian is derived analytically if needed and thus the analytical Jacobian does not need to be provided.

))
No optimal solution found. The solver message was: Convergence due to lack of progress in parameter updates. The objective value is: 1.314286e+308

R> solution(nlp_1_sol_2)
[1] NA NA NA There are several possibilities to help the algorithm to find a good solution. In almost all practical applications it is possible to specify lower and upper bounds. Carefully chosen bounds can improve the quality of the solution and decrease the runtime of the algorithm. The runtime of the algorithm also strongly depends on the tolerances set; if the tolerances are set too small the algorithm will reach the maximum number of iterations before convergence.

R>
Optimal solution found. The objective value is: 3.456000e+03 R> solution(nlp_1_sol_4) [1] 23.99997 11.99991 12.00011 In general it is often hard to obtain a global optimum for non-convex optimization problems.
Often the only option is to try different algorithms, parameters and starting values and hope that one of the solutions is a global optimum. ROI lowers the burden to compare different algorithms and therefore can assist in finding a global solution.

Extending ROI
To stay abreast of changes and to further the availability of different solvers in the ROI ecosystem, ROI allows developers to integrate their own extensions, so called plug-ins. This can be seen as a key feature, since it allows the use of new solvers with no or minimal code changes.
Extending ROI with a new solver method can be split into three parts. First, a function to be called by ROI has to be written. Second, the function plus information about the function are added into the ROI solver registry. Third, a mapping from the solver specific arguments and the status codes to their ROI counterpart has to be provided.

Signatures
In order to establish a connection between the OP and the solvers provided via plug-ins, both are equipped with a signature. The signature captures all the information necessary to determine which solver is applicable to a given problem.

R> OP_signature(lp)
objective constraints bounds cones maximum New signatures are created by the function ROI_plugin_make_signature(). The following shows how to create the signature for the "glpk" solver, R> glpk_signature <-ROI_plugin_make_signature(objective = "L", + constraints = "L", types = c("C", "I", "B", "CI", "CB", "IB", "CIB"), + bounds = c("X", "V"), maximum = c(TRUE, FALSE)) where the objective and the constraints have to be linear. Furthermore, this signature indicates that, the variable types are allowed to be binary ("B"), integer ("I"), continuous ("C") or any combinations thereof. The bounds have to be variable bounds ("V") or no bounds at all, encoded by "X". The last argument maximum specifies that GLPK can find both maxima and minima.

Writing a new solver method
Any function supposed to add a solver to ROI has to take the arguments x and control, where x is of class 'OP' and control a list containing the additional arguments. Furthermore, the solution has to be canonicalized before it is returned. The following shows the code from ROI.plugin.glpk for solving linear problems.

Registering solver methods
Registering a solver method can be seen as telling ROI which function it should use when ROI_solve() with argument solver set to the name of the plug-in (e.g., "glpk") is called. In order to avoid ambiguity, each plug-in should at most provide one method for each problem type. Solver methods are registered via the ROI_plugin_register_solver_method() function, which takes as arguments the problem types (as signatures), the solver name and a wrapper function, ROI_solve() is dispatched to.
The following code registers the solver, but is not executed as the entry already exists in the registry.

ROI_plugin_register_solver_method(glpk_signature, "glpk", glpk_solve_OP)
After the solver registration the name of the solver will appear among the registered solvers.

Adding additional information
To be able to provide a consistent interface, each plug-in has to define a mapping between the solver specific status codes and the status codes used by ROI, as well as a mapping between solver specific control variables and ROI control variables.

Status codes
Status codes can be added via the function ROI_plugin_add_status_code_to_db(). The following code is not executed as the entry already exists in the registry.

Control variables
Plug-ins are contracted to provide a mapping between the names of the control variables used by ROI and the names of the control variables used by the plug-in. The following maps the "glpk" argument tm_limit to the ROI equivalent max_time. The following code is not executed as the entry already exists in the registry.

Registering reformulations
While in Section 5.3 we showed how to use reformulations, here we explain how new reformulations can be added through plug-ins. Again, the signature is used to define which transformations can be performed by a given method.

Registering reader/writer
Plug-ins can also add new read and write functions. Any method to be registered as read function has to take as arguments file, the file name, and ..., for optional additional arguments. The write functions need the additional argument x, which is the OP to be written out.

ROI Optimization Problem:
Maximize a linear objective function of length 2 with -2 continuous objective variables, subject to -2 constraints of type linear. -2 lower and 2 upper non-standard variable bounds.

ROI tests
Writing tests is an important task in software development. The ROI.tests (Schwendinger 2019b) package provides a collection of tests which should be applied to any ROI plug-in during development. Since ROI knows the signature of each solver, ROI.tests can select the appropriate tests based on the solver name.

Applications
In the following we demonstrate how ROI can be used to solve selected problems from statistics, namely L 1 regression, best subset selection, relative risk regression, sum-of-norms clustering, and graphical lasso.

L 1 regression
The linear programming formulation of the L 1 regression problem as shown in Section 2.1 can be constructed using ROI methods via the following R function.

Best subset selection
Recently Bertsimas et al. (2016) reported a bewildering 450 billion factor speedup from 1991 to 2015 for solving MIP, which is partly due to algorithmic improvements and partly because of hardware speedups. They show how this speed gain can be utilized to solve the best subset selection problem in regression (see, for example, Miller 2002), which is an NP-hard combinatorial OP. The best subset selection problem is a variable selection scheme which extends linear least-squares by adding a constraint on the number of predictor variables.
As Equation 34 suggests, the best subset selection is in spirit similar to ridge regression and lasso. However instead of the l 2 and l 1 norm best subset selection uses the l 0 norm which makes it non-convex and therefore hard to solve. The leaps (Lumley 2020) package implements an efficient branch-and-bound algorithm which is significantly faster than exhaustive search. Using an optimization solver has the additional advantage that it is possible to impose additional restrictions, e.g., if the quadratic term of a covariate is selected to be in the equation the linear term has also to be selected. In ROI best subset selection can be either implemented as mixed integer quadratic problem or as mixed integer second order cone problem. This problem can be solved with a mixed integer QP solver or a mixed integer SCOP solver. An implementation of the second order cone version can be found in the supplementary material.

Relative risk regression
Generalized linear models (GLMs) are often the statisticians' first choice for regression analysis of binary response data. The most prominent model of the GLM family is logistic regression. 3 A GLM which is mainly used in epidemiology but closely related to logistic regression is relative risk regression. The main difference between these two models is the fact that relative risk regression uses the log link and therefore estimates relative risks instead of odds ratios. Lumley, Kronmal, and Ma (2006) reviews the algorithms proposed in the literature to perform relative risk regression. Luo, Zhang, and Sun (2014) suggest to use a log-binomial model with Here X i * refers to the ith row of the data matrix X, the constraint Xβ ≤ 0 ensures that the estimated probabilities are in the interval [0, 1]. The log-binomial regression model can be formulated as a general nonlinear optimization problem: R> logbin_gps <-function(y, X) { + loglikelihood <-function(beta) { + xb <-drop(X %*% beta) + if (any(xb > 0)) NaN else sum(y * xb + (1 -y) * log(1 -exp(xb))) + } + + gradient <-function(beta) { + exb <-exp(drop(X %*% beta)) + drop(crossprod(X, (y -exb) / (1 -exb))) + } + + OP(F_objective(loglikelihood, n = ncol(X), G = gradient), + L_constraint(L = X, dir = leq(nrow(X)), rhs = double(nrow(X))), + bounds = V_bound(ld = -Inf, nobj = ncol(X)), maximum = TRUE) + } which can be solved by any GPS which allows to specify linear constraints (e.g., using package alabama). Alternatively the log-binomial regression model can be solved by any CP solver which supports the linear and the primal exponential cone (e.g., using package scs). The latter formulation has attractive theoretical and numerical convergence guarantees without the need to find suitable starting values.

Sum-of-norms clustering
Borrowing ideas from regularization, sum-of-norms (SON) clustering (convex clustering) is an interesting alternative to established clustering approaches like hierarchical or k-means clustering, which has attracted a lot of research in recent years (Pelckmans, De Brabanter, De Moor, and Suykens 2005;Lindsten, Ohlsson, and Ljung 2011;Hocking, Joulin, Bach, and Vert 2011;Zhu, Xu, Leng, and Yan 2014;Chi and Lange 2015;Tan, Witten et al. 2015). Pelckmans et al. (2005) and Hocking et al. (2011) describe SON clustering as a convexification of hierarchical clustering and Lindsten et al. (2011) establish that SON clustering can be seen as a convex relaxation of k-means clustering. Due to its convexity, SON clustering is not dependent on the starting values, which is a clear advantage over the non-convex k-means and hierarchical clustering.
SON clustering solves the following convex OP, where q ∈ {1, 2, ∞}, X ∈ R m×n is the data matrix and M i * the ith row of the optimization variable M . The regularization term λ i<j M i * − M j * q induces equal rows M i * . For λ = 0 all rows are unique and M is equal to X, when λ increases the number of unique rows of M will decrease. This gives a clustering where all equal rows belong to the same cluster. By solving Equation 36 for different λ i , where λ 1 < λ 2 < · · · < λ k−1 < λ k , one can obtain a hierarchical clustering tree (Pelckmans et al. 2005).
At least two implementations of SON clustering exist in R, Hocking et al. (2011) provide their code on R-Forge (Hocking 2015) and Chi and Lange (2015) provide a fast implementation of SON clustering on CRAN (Chi and Lange 2014). A ROI formulation as SOCP of SON clustering can be found in the supplementary material.

Graphical lasso
Obtaining good estimates of the covariance matrix Σ is important in modern statistics. Often Σ is not estimated directly but its inverse, the precision matrix Θ = Σ −1 (e.g., Meinshausen and Bühlmann 2006). Estimating the precision matrix instead of the covariance matrix has the advantage that there is a direct connection between the precision matrix and Gaussian graphical models, in the sense that the precision matrix defines the structure of the Gaussian graphical model. The elements of the precision matrix are the partial correlations, Θ ij is zero if and only if i and j are conditionally independent. Translated to Gaussian graphical models, two edges A and B are only connected if the corresponding entry in the precision matrix is non-zero (Lauritzen 1996).
Several authors proposed an algorithm connected to the lasso (Tibshirani 1996), to obtain a sparse estimate of the precision matrix, the so-called graphical lasso (glasso; Friedman, Hastie, and Tibshirani 2008). The glasso solves the following convex OP, minimize Θ 0 − log(det(Θ)) + tr(S Θ) + λ X 1 , where the data matrix X ∈ R n×p is assumed to be generated from a p-dimensional multivariate normal distribution N p (µ, Σ) and S is the sample covariance matrix of X. Making use of the exponential and semidefinite cone, this can be brought into the CP standard form and solved by ROI using SCS. The corresponding R code can be found in the supplementary material.

Conclusions
In this paper we presented the ROI package and its extensions. ROI provides a consistent way to model OPs in R. ROI makes strong use of R's generic functions, such that users already familiar with R are not obliged to learn a new language. The plug-in packages equip ROI with optimization solvers and predefined optimization models. ROI is currently applicable to linear, quadratic, conic and general nonlinear OPs and provides access to nineteen solvers and three model plug-ins.
We illustrated how ROI can be used to solve OPs from many different problem classes. Furthermore, we have shown how ROI can be used to solve challenging statistical problems like best subset selection, convex clustering and glasso. The plug-in package ROI.plugin.msbinlp (Hornik, Meyer, and Schwendinger 2017b) serves as an example to highlight the benefit of the development of new packages based on ROI. As package ROI.plugin.msbinlp shows, the main benefit is that there is no need to have a dependency on a specific solver which could be also interesting for the implementation of nonlinear optimization algorithms (e.g., sequential quadratic programming). Another benefit is that package authors can reuse test cases from other packages based on ROI and the plug-in package ROI.tests provides a standardized way to test new solver plug-ins.
Although ROI already is able to cope with a wide range of optimization problems, there are still many possibilities for extensions. These include extending the supported functional forms of the objective functions and constraints as well as adding additional solvers or model collections through plug-ins. The following gives an overview of the planned extensions and possible future work: • Currently, ROI is lacking a plug-in capable of solving general nonlinear optimization problems with mixed integer constraints. Within the COIN-OR project the Couenne (Belotti, Lee, Liberti, Margot, and Wächter 2009) and the Bonmin (Bonami and Lee 2013) solvers are designed to try to obtain a global solution for non-convex MINLPs. Therefore both solvers would be a valuable extension of ROI.
• Another popular solver from the COIN-OR project is Ipopt (Wächter and Biegler 2006), which aims to solve general nonlinear optimization problems. As mentioned before, there exists already the ipoptr package, but since there are many steps needed for the installation we assume it is currently not accessible to many R users. Therefore a ROI plugin installable directly from CRAN would be preferable.
• Add a reader for the QPLIB file format (Furini et al. 2017).
• Explore possibilities of supervised solver recommendation.
We are working on extending the amount of solvers and resources available through ROI.