This R session will introduce different solvers in R that can be used for portfolio optimization.

For an overwhelming amount of information on the many solvers for R, see the Task View on Optimization and Mathematical Programming.

John Nash, the author of many of the optimization packages (not to be confused with the main character of the movie A Beautiful Mind), has an interesting and insightful blog: https://nashjc.wordpress.com. In particular, the entry on choosing which method to use for “optimizing” functions is recommended as a guide on how to navigate the overwhelming amount of solvers in R.

General solvers

General solvers can deal with arbitrary nonlinear optimization problems, at the expense of perhaps slow convergence and risk of numerical issues.

Package stats

Package stats (base R package installed by default) offers several general purpose optimization routines.

  • optimize(): For one-dimensional unconstrained function optimization over an interval (for one-dimensional root finding use uniroot()).

  • optim(): General-purpose optimization with six different optimization methods (see a conversation with the author):
    • Nelder-Mead: relatively robust method (default) that does not require derivatives (function nmkb() of package dfoptim is probably better)
    • CG: low-memory optimizer for unconstrained problems with high dimensions (albeit generally not satisfactory, better use package Rcgmin)
    • BFGS: simple unconstrained quasi-Newton method (package Rvmmin is probably better)
    • L-BFGS-B: modest-memory optimizer for bounds-constrained problems
    • SANN: simulated annealing method
    • Brent: for one-dimensional problems (it actually calls optimize())

This example does a least squares fitting: \[ \begin{array}{ll} \underset{a,b}{\textsf{minimize}} & \sum_{i=1}^{m} \left(y_i - (a+bx_i)\right)^2 \end{array} \]

The next example illustrates the use of the gradient with the famous Rosenbrock banana function: \[ f(\mathbf{x}) = 100(x_2-x_1^2)^2 + (1-x_1)^2 \] with gradient \[ \nabla f(\mathbf{x}) = \left( \begin{array}{c} -400x_1(x_2-x_1^2) -2(1-x_1)\\ 200(x_2-x_1^2) \end{array}\right) \] with the unconstrained minimization problem \[ \begin{array}{ll} \underset{\mathbf{x}\in \mathbf{R}^2}{\textsf{minimize}} & 100(x_2-x_1^2)^2 + (1-x_1)^2 \end{array} \]

The following example uses box constraints: \[ \begin{array}{ll} \underset{\mathbf{x}\in \mathbf{R}^{25}}{\textsf{minimize}} & [1, 4, 4, \ldots, 4] \times \left[ \begin{array}{c} (x_1-1^2)^2\\ (x_2-x_1^2)^2\\ (x_3-x_2^2)^2\\ \cdots\\ (x_{25}-x_{24}^2)^2\end{array}\right]\\ {\textsf{subject to}} & 2 \le x_i \le 4, \qquad i=1,\ldots,25 \end{array} \]

This example uses the simulated annealing method (for global optimization):

  • constrOptim(): Minimize a function subject to linear inequality constraints using an adaptive barrier algorithm (it calls optim() under the hood).
  • nlm(): This function carries out a minimization of the objective function using a Newton-type algorithm.
  • nlminb(): Unconstrained and box-constrained optimization using PORT routines.

Package optimr

The author of the base function optim() (which he claims is obsolete now) has also created the package optimr as a convenient wrapper to many other solvers so that they can be easily used and compared:

The package optimrx is similar to optimr but accepts more methods. It is not available in CRAN, instead check R-Forge and to install do: install.packages("optimrx", repos="http://R-Forge.R-project.org").

Package gloptim

Package gloptim is a wrapper to many global optimization packages (similar in spirit to what optimr is lo local optimization). Note that global optimization is totally different in philosophy to local optimization (global optimization solvers are usually called stochastic solvers and attempt to escape local optimum points). To install, do: install.packages("gloptim", repos="http://R-Forge.R-project.org").

Solvers for specific classes of problems

If the problem to be solved falls within a particular class of problems, such as LS, LP, MILP, QP, SOCP, or SDP, then it is much better to use a solver specific for that class.

Least Squares (LS)

A linear least-squares (LS) problem minimizes the \(\ell_2\)-norm \(\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_2\) possibly with bounds or linear constraints. The function qr.solve(A, b) from base R solves over- and underdetermined linear systems in the least-squares sense. The package colf performs least squares constrained optimization on a linear objective function and contains a number of algorithms to choose from.

Linear Programming (LP)

The package linprog contains the function solveLP() to conveniently solve LPs of the form: \[ \begin{array}{ll} \underset{\mathbf{x}\in \mathbf{R}^n}{\textsf{minimize}} & \mathbf{c}^T\mathbf{x}\\ {\textsf{subject to}} & \begin{array}[t]{l} \mathbf{A}\mathbf{x} \le \mathbf{b}\\ \mathbf{x} \ge \mathbf{0}\end{array}\\ \end{array} \]

Quadratic Programming (QP)

The package quadprog contains the function solve.QP() to conveniently solve QPs of the form: \[ \begin{array}{ll} \underset{\mathbf{x}\in \mathbf{R}^n}{\textsf{minimize}} & -\mathbf{d}^T\mathbf{x} + \frac{1}{2}\mathbf{x}^T\mathbf{D}\mathbf{x}\\ {\textsf{subject to}} & \mathbf{A}^T\mathbf{x} =/\ge \mathbf{b}\\ \end{array} \]

The package quadprogXT extends the quadprog to solve quadratic programs with absolute value constraints and absolute values in the objective function.

Second-Order Cone Programming (SOCP)

There are several packages:

  • Package ECOSolveR provides an interface to the Embedded COnic Solver (ECOS), a well-known, efficient, and robust C library for convex problems. Conic and equality constraints can be specified in addition to integer and boolean variable constraints for mixed-integer problems.

  • The CLSOCP package provides an implementation of a one-step smoothing Newton method for the solution of SOCP problems.

Cone Programming and Semi-Definite Programming (SDP)

There are several good packages:

  • Package scs applies operator splitting to solve linear programs, cone programs (SOCP), and semidefinite programs; cones can be second-order, exponential, power cones, or any combination of these.

  • Package sdpt3r solves general semidefinite programs using an R implementation of SDPT3, a MATLAB software for semidefinite quadratic-linear programming.

  • Package Rcsdp is an interface to the CSDP semidefinite programming library of routines that implement a primal-dual barrier method for solving semidefinite programming problems.

  • Package Rmosek provides an interface to the (commercial) MOSEK optimization library for large-scale LP, QP, and MIP problems, with emphasis on (nonlinear) conic, semidefinite, and convex tasks; academic licenses are available.

Optimization infrastructure packages

We have already seen two packages that are wrappers to many other solvers, namely: optimr for local optimization and gloptim for global optimization.

It is also worth mentioning the NEOS Server for optimization that provides online access to state-of-the-art optimization problem solvers. It is a free internet-based service for solving numerical optimization problems and provides access to more than 60 state-of-the-art solvers. The R package rneos enables the user to pass optimization problems to NEOS and retrieve results within R.

The following two packages are extremely useful and powerful: ROI and CVXR.

ROI for Convex Problems, MIP, and Nonconvex Problems

The ROI package for R Optimization Infrastructure provides a framework for handling optimization problems in R. It uses an object-oriented approach to define and solve various optimization tasks in R which can be from different problem classes (e.g., linear, quadratic, non-linear programming problems). This makes optimization transparent for the R user as the corresponding workflow is completely abstracted from the underlying solver. The approach allows for easy switching between solvers and thus enhances comparability. See ROI paper for details. Let’s see now several examples.

LP – Consider the LP: \[ \begin{array}{ll} \underset{\mathbf{x}}{\textsf{maximize}} & 3x_1 + 7x_2 - 12x_3\\ \textsf{subject to} & \begin{array}[t]{l} 5x_1 + 7x_2 + 2x_3 &\le 61\\ 3x_1 + 2x_2 - 9x_3 &\le 35\\ x_1 + 3x_2 + x_3 &\le 31\\ x_1,x_2 \ge 0, \quad x_3 \in &[-10,10] \end{array} \end{array} \]

library(ROI)  #install.packages("ROI")
#> ROI: R Optimization Infrastructure
#> Registered solver plugins: nlminb, ecos, lpsolve, scs.
#> Default solver: auto.

prob <- OP(objective = L_objective(c(3, 7, -12)),
           constraints = L_constraint(L = rbind(c(5, 7,  2), 
                                                c(3, 2, -9), 
                                                c(1, 3,  1)),
                                      dir = c("<=", "<=", "<="), 
                                      rhs = c(61, 35, 31)),
           bounds = V_bound(li = 3, ui = 3, lb = -10, ub = 10, nobj = 3),
           maximum = TRUE)
prob
#> 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.

# let's look at available solvers
ROI_available_solvers(prob)$Package  # solvers available for the particular problem on the internet
#>  [1] "ROI.plugin.alabama"  "ROI.plugin.ecos"     "ROI.plugin.lpsolve"  "ROI.plugin.neos"     "ROI.plugin.alabama" 
#>  [6] "ROI.plugin.cplex"    "ROI.plugin.ecos"     "ROI.plugin.glpk"     "ROI.plugin.gurobi"   "ROI.plugin.ipop"    
#> [11] "ROI.plugin.lpsolve"  "ROI.plugin.mosek"    "ROI.plugin.neos"     "ROI.plugin.nloptr"   "ROI.plugin.qpoases" 
#> [16] "ROI.plugin.scs"      "ROI.plugin.symphony" "ROI.plugin.cbc"      "ROI.plugin.clp"      "ROI.plugin.gurobi"  
#> [21] "ROI.plugin.mosek"
library(ROI.plugin.ecos)  #install.packages("ROI.plugin.ecos")
library(ROI.plugin.lpsolve)  #install.packages("ROI.plugin.lpsolve")
ROI_installed_solvers()  # solvers installed (not necessarily loaded)
#>               nlminb                 ecos              lpsolve                  scs 
#>  "ROI.plugin.nlminb"    "ROI.plugin.ecos" "ROI.plugin.lpsolve"     "ROI.plugin.scs"
ROI_registered_solvers()  # solvers installed and loaded (registered) in system
#>               nlminb                 ecos              lpsolve                  scs 
#>  "ROI.plugin.nlminb"    "ROI.plugin.ecos" "ROI.plugin.lpsolve"     "ROI.plugin.scs"
ROI_applicable_solvers(prob) # solvers installed and loaded in system applicable to the particular problem
#> [1] "ecos"    "lpsolve" "scs"

# solve it
res <- ROI_solve(prob)
res
#> Optimal solution found.
#> The objective value is: 8.670149e+01
solution(res)            # also: res$solution
#> [1]  1.941150e-10  9.238806e+00 -1.835821e+00
solution(res, "objval")  # also: res$objval
#> [1] 86.70149
solution(res, "status")
#> $code
#> [1] 0
#> 
#> $msg
#>   solver ecos
#>     code 0
#>   symbol ECOS_OPTIMAL
#>  message Optimal solution found.
#> roi_code 0

MILP – Consider the previous LP and make it into an MILP by adding the constraints \(x_2,x_3 \in \mathbb{Z}\).

BLP – Consider the binary linear program (BLP): \[ \begin{array}{ll} \underset{\mathbf{x}}{\textsf{minimize}} & -x_1-x_2-x_3-x_4-99x_5\\ \textsf{subject to} & \begin{array}[t]{l} x_1 + x_2 & \le 1\\ x_3 + x_4 & \le 1\\ x_4 + x_5 & \le 1\\ x_i\in \{0,1\} \end{array} \end{array} \]

SOCP – Consider the SOCP: \[ \begin{array}{ll} \underset{\mathbf{x}}{\textsf{maximize}} & x + y + z\\ \textsf{subject to} & \begin{array}[t]{l} \sqrt{x^2 + z^2} \le \sqrt{2}\\ x,y,z \ge 0, \quad y \le 5\end{array} \end{array} \] and note that the SOC constraint \(\sqrt{x^2 + z^2} \le \sqrt{2}\) can be written as \(\|(x,z)\|_2 \le \sqrt{2}\) or \((\sqrt{2}, (x, z)) \in \mathcal{K}_\textsf{soc(3)}\) which is achieved in the code as: \(\textsf{rhs} - \mathbf{L}\mathbf{x} = \left[\begin{array}{c} \sqrt{2}\\ x\\ z\end{array} \right] \in \mathcal{K}_\textsf{soc(3)}\).

SDP – Consider the SDP: \[ \begin{array}{ll} \underset{\mathbf{x}}{\textsf{minimize}} & x_1 + x_2 - x_3\\ \textsf{subject to} & \begin{array}[t]{l} x_1\left(\begin{array}{cc} 10 & 3\\ 3 & 10\end{array} \right) + x_2\left(\begin{array}{cc} 6 & -4\\ -4 & 10\end{array} \right) + x_3\left(\begin{array}{cc} 8 & 1\\ 1 & 6\end{array} \right) \preceq \left(\begin{array}{cc} 16 &-13\\-13 & 60\end{array} \right)\\ x_1,x_2,x_3 \ge 0\end{array} \end{array} \] and note that the SDP constraint \(\sum_{i=1}^m x_i\mathbf{F}_i \preceq \mathbf{F}_0\) can be written as \(\mathbf{F}_0 - \sum_{i=1}^m x_i\mathbf{F}_i \in \mathcal{K}_\textsf{sdp(3)}\) (size 3 is because in our problem the matrices are \(2\times2\) but vech() extracts 3 independent variables since the matrices are symmetric).

NLP – Consider the nonlinear program (NLP): \[ \begin{array}{ll} \underset{\mathbf{x}}{\textsf{maximize}} & x_1 x_2 x_3\\ \textsf{subject to} & \begin{array}[t]{l} x_1 + 2x_2 + 2x_3 \le 72\\ x_1,x_2,x_3 \in [0,42] \end{array} \end{array} \]

CVXR for Convex Problems

CVXR is an R package that provides an object-oriented modeling language for convex optimization, similar to CVX, CVXPY, YALMIP, and Convex.jl. It allows the user to formulate convex optimization problems in a natural mathematical syntax rather than the restrictive standard form required by most solvers. The user specifies an objective and set of constraints by combining constants, variables, and parameters using a library of functions with known mathematical properties. CVXR then applies signed disciplined convex programming (DCP) to verify the problem’s convexity. Once verified, the problem is converted into standard conic form using graph implementations and passed to a cone solver such as ECOS or SCS. See CVX manual for details. Let’s see now several examples.

Least squares – Let’s start with a simple LS example: \[ \begin{array}{ll} \underset{\boldsymbol{\beta}}{\textsf{minimize}} & \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2\\ \end{array} \] We can of course the R base linear model fitting function lm():

Now let’s do it with CVXR:

We can now easily add a contraint to solve a nonnegative LS:

We can add other constraints:

Robust Huber regression – Let’s consider another simple example of robust regression: \[ \begin{array}{ll} \underset{\boldsymbol{\beta}}{\textsf{minimize}} & \sum_{i=1}^m \phi(y_i - \mathbf{x}_i^T\boldsymbol{\beta}) \end{array} \] where \[ \phi(u) = \begin{cases} u^2 & \text{if }|u| \le M \\ 2Mu-M^2 & \text{if }|u| > M \end{cases} \]

To see a list of functions available in the CVXR package, see https://cvxr.rbind.io/cvxr_functions/.

Elastic net regularization – We will now solve the problem: \[ \begin{array}{ll} \underset{\boldsymbol{\beta}}{\textsf{minimize}} & \frac{1}{2m}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \lambda \left( (1-\alpha)\|\boldsymbol{\beta}\|_2^2 + \alpha\| \boldsymbol{\beta} \|_1 \right) \end{array} \]

Sparse inverse covariance estimation – Consider the matrix-valued convex problem: \[ \begin{array}{ll} \underset{\mathbf{X}}{\textsf{maximize}} & \textsf{logdet}(\mathbf{X}) - \textsf{Tr}(\mathbf{X}\mathbf{S})\\ \textsf{subject to} & \mathbf{X} \succeq \mathbf{0}, \qquad \sum_{i=1}^n\sum_{j=1}^n\left|X_{ij}\right| \le \alpha \end{array} \]

Worst case covariance – Consider the matrix-valued convex problem: \[ \begin{array}{ll} \underset{\boldsymbol{\Sigma}}{\textsf{maximize}} & \mathbf{w}^T \boldsymbol{\Sigma} \mathbf{w}\\ \textsf{subject to} & \boldsymbol{\Sigma} \succeq \mathbf{0}, \qquad L_{ij} \le \Sigma_{ij} \le U_{ij} \end{array} \]

Portfolio optimization – Consider a Markowitz portfolio design: \[ \begin{array}{ll} \underset{\mathbf{w}}{\textsf{maximize}} & \boldsymbol{\mu}^T\mathbf{w} - \gamma \mathbf{w}^T \boldsymbol{\Sigma} \mathbf{w}\\ \textsf{subject to} & \mathbf{w} \ge \mathbf{0}, \quad \mathbf{1}^T\mathbf{w}=1 \end{array} \]

Conclusion

The amount of available solvers in R is overwhelming and it may be difficult to navigate. The following steps are suggested:

  • If the problem is convex, then start with CVXR for initial tests.
  • If the speed is not fast enough, then move on to the next level and use ROI at the expense of more effort to use.
  • If faster speed is still required, then if the problem belongs to one of the well defined classes, such as LP, MILP, QP, SOCP, or SDP, use a solver specific for that class (for example, for LPs I recommend lpSolve and for QPs quadprog).
  • However, if the problem does not fit any class, then one has to use a general solver for nonlinear optimization. In that sense, if a local solution is enough, then a good starting point is either ROI or the package optimr, which is a wrapper for many solvers. If a global solver is required, then the package gloptim is a good choice, which is a wrapper for many global solvers.