15 List of examples

List of examples shipped in the distribution of Rmosek package:

Table 15.1 List of distributed examples

File

Description

affco1.R

A simple problem using affine conic constraints

affco2.R

A simple problem using affine conic constraints

ceo1.R

A simple conic exponential problem

cqo1.R

A simple conic quadratic problem

gp1.R

A simple geometric program (GP) in conic form

helloworld.R

A Hello World example

lo1.R

A simple linear problem

logistic.R

Implements logistic regression and simple log-sum-exp (CEO)

mico1.R

A simple mixed-integer conic problem

milo1.R

A simple mixed-integer linear problem

mioinitsol.R

A simple mixed-integer linear problem with an initial guess

normex.R

Demonstrates least squares and other norm minimization problems

parameters.R

Shows how to set optimizer parameters and read information items

portfolio_1_basic.R

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.R

Portfolio optimization - efficient frontier

portfolio_3_impact.R

Portfolio optimization - market impact costs

portfolio_4_transcost.R

Portfolio optimization - transaction costs

portfolio_5_card.R

Portfolio optimization - cardinality constraints

portfolio_6_factor.R

Portfolio optimization - factor model

pow1.R

A simple power cone problem

qo1.R

A simple quadratic problem

reoptimization.R

Demonstrate how to modify and re-optimize a linear problem

response.R

Demonstrates proper response handling

sdo1.R

A simple semidefinite problem with one matrix variable and a quadratic cone

sdo2.R

A simple semidefinite problem with two matrix variables

sdo_lmi.R

A simple semidefinite problem with an LMI using the SVEC domain.

simple.R

A simple I/O example: read problem from a file, solve and write solutions

solutionquality.R

Demonstrates how to examine the quality of a solution

Additional examples can be found on the MOSEK website and in other MOSEK publications.

affco1.R

Listing 15.1 affco1.R Click here to download.
##
#  File : affco1.R
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  Description :
#    Implements a basic tutorial example with affine conic constraints:
#
#    maximize x_1^(1/3) + (x_1+x_2+0.1)^(1/4)
#    st.      (x_1-0.5)^2 + (x_2-0.6)^2 <= 1
#             0 <= x_1 <= x_2 + 1
# 
##
library("Rmosek")

affco1 <- function()
{

    prob <- list(sense="max")

    # Variables [x1; x2; t1; t2]
    prob$c <- c(0, 0, 1, 1)

    # Linear inequality x_1 - x_2 <= 1
    prob$A <- Matrix(c(1, -1, 0, 0), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=-Inf, buc=1)
    prob$bx <- rbind(blx=rep(-Inf,4), bux=rep(Inf,4))

    # The quadratic cone
    FQ <- rbind(c(0,0,0,0), c(1,0,0,0), c(0,1,0,0))
    gQ <- c(1, -0.5, -0.6)
    cQ <- matrix(list("QUAD", 3, NULL), nrow=3, ncol=1)

    # The power cone for (x_1, 1, t_1) \in POW3^(1/3,2/3)
    FP1 <- rbind(c(1,0,0,0), c(0,0,0,0), c(0,0,1,0))
    gP1 <- c(0, 1, 0)
    cP1 <- matrix(list("PPOW", 3, c(1/3, 2/3)), nrow=3, ncol=1)

    # The power cone for (x_1+x_2+0.1, 1, t_2) \in POW3^(1/4,3/4)
    FP2 <- rbind(c(1,1,0,0), c(0,0,0,0), c(0,0,0,1))
    gP2 <- c(0.1, 1, 0)
    cP2 <- matrix(list("PPOW", 3, c(1.0, 3.0)), nrow=3, ncol=1)

    # All cones
    prob$F <- rbind(FQ, FP1, FP2)
    prob$g <- cbind(gQ, gP1, gP2)
    prob$cones <- cbind(cQ, cP1, cP2)
    rownames(prob$cones) <- c("type","dim","conepar")

    r <- mosek(prob, list(soldetail=1))
    stopifnot(identical(r$response$code, 0))

    print(r$sol$itr$pobjval)
    print(r$sol$itr$xx[1:2])
}

affco1()

affco2.R

Listing 15.2 affco2.R Click here to download.
##
#  File : affco2.R
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  Description :
#    Implements a basic tutorial example with affine conic constraints:
#
#    minimize t
#    st.      (d, z1*y1,... zn*yn) \in Q^{n+1}
#             (yi, 1, ai*t)        \in EXP, i=1,\ldots,n
#
#    with input ai<0, zi, d.
#
#    See also https://docs.mosek.com/modeling-cookbook/expo.html#hitting-time-of-a-linear-system
# 
##
library("Rmosek")

firstHittingTime <- function(n, z, a, d)
{
    prob <- list(sense="min")
    # Variables [t, y1, ..., yn]
    prob$A <- Matrix(nrow=0, ncol=n+1)
    prob$bx<- rbind(rep(-Inf,n+1), rep(Inf,n+1))
    prob$c <- c(1, rep(0, n))

    # Quadratic cone
    FQ <- Diagonal(n+1, c(0, z))
    gQ <- c(d, rep(0, n))

    # All exponential cones
    FE <- sparseMatrix(  c(seq(1,3*n,by=3), seq(3,3*n,by=3)),
                         c(2:(n+1),         rep(1,n)),
                       x=c(rep(1,n),        t(a)))
    gE <- rep(c(0, 1, 0), n)

    # Assemble input data
    prob$F <- rbind(FQ, FE)
    prob$g <- c(gQ, gE)
    prob$cones <- cbind(matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1),
                        matrix(list("PEXP", 3, NULL), nrow=3, ncol=n))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(soldetail=1))
    stopifnot(identical(r$response$code, 0))

    r$sol$itr$xx[1]
}

n = 3
z = c(3.3, 2.2, 1.3)
a = c(-0.08, -0.3, -0.06)
d = 0.5

t <- firstHittingTime(n, z, a, d)
print(sprintf("t = %.4e", t))

ceo1.R

Listing 15.3 ceo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      ceo1.R
#
#  Purpose :   To demonstrate how to solve a small conic exponential
#              optimization problem using the Rmosek package.
##
library("Rmosek")

ceo1 <- function()
{
    # Specify the non-conic part of the problem.
    prob <- list(sense="min")
    prob$c  <- c(1, 1, 0)
    prob$A  <- Matrix(c(1, 1, 1), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=1, 
                     buc=1)
    prob$bx <- rbind(blx=rep(-Inf,3), 
                     bux=rep( Inf,3))
    
    # Specify the affine conic constraints.
    prob$F <- rbind(c(1,0,0),c(0,1,0),c(0,0,1))
    prob$g <- c(0, 0, 0)
    prob$cones <- matrix(list("PEXP", 3, NULL), nrow=3, ncol=1)
    rownames(prob$cones) <- c("type","dim", "conepar")

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

ceo1()

cqo1.R

Listing 15.4 cqo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      cqo1.R
#
#  Purpose :   To demonstrate how to solve a small conic quadratic
#              optimization problem using the Rmosek package.
##
library("Rmosek")

cqo1 <- function()
{
    # Specify the non-conic part of the problem.
    prob <- list(sense="min")
    prob$c  <- c(0, 0, 0, 1, 1, 1)
    prob$A  <- Matrix(c(1, 1, 2, 0, 0, 0), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=1, 
                     buc=1)
    prob$bx <- rbind(blx=c(rep(0,3), rep(-Inf,3)), 
                     bux=rep(Inf,6))
    
    # Specify the affine conic constraints.
    # NOTE: The F matrix is internally stored in the sparse
    #       triplet form. Use 'giveCsparse' or 'repr' option 
    #       in the sparseMatrix() call to construct the F 
    #       matrix directly in the sparse triplet form. 
    prob$F     <- sparseMatrix(i=c(1, 2, 3, 4, 5, 6),
                               j=c(4, 1, 2, 5, 6, 3), 
                               x=c(1, 1, 1, 1, 1, 1),
                               dims = c(6,6))
    prob$g     <- c(1:6)*0
    prob$cones <- matrix(list(), nrow=3, ncol=2)
    rownames(prob$cones) <- c("type","dim","conepar")

    prob$cones[,1] <- list("QUAD", 3, NULL)
    prob$cones[,2] <- list("RQUAD",3, NULL)
    
    #
    # Use cbind to extend this chunk of ACCs, if needed:
    #
    #    oldcones <- prob$cones
    #    prob$cones <- cbind(oldcones, newcones)
    #

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

cqo1()

gp1.R

Listing 15.5 gp1.R Click here to download.
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      gp1.R
#
#  Purpose:    Demonstrates how to solve a simple Geometric Program (GP)
#              cast into conic form with exponential cones and log-sum-exp.
#
#              Example from
#                https://gpkit.readthedocs.io/en/latest/examples.html//maximizing-the-volume-of-a-box
#

#
# maximize     h*w*d
# subjecto to  2*(h*w + h*d) <= Awall
#                       w*d <= Afloor
#              alpha <= h/w <= beta
#              gamma <= d/w <= delta
#
# Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
#
# maximize     x+y+z
# subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
#                              y+z <= log(Afloor)  
#              log( alpha ) <= x-y <= log( beta )
#              log( gamma ) <= z-y <= log( delta )
#
#
# log(sum(exp(A*x+b))) <= 0 is equivalent to sum(u) == 1, (ui,1,ai*x+bi) in Kexp,
# so we have two exp-cones and two auxilliary variables u1,u2.
#
# We order variables as (x,y,z,u1,u2),

library("Rmosek")

gp1 <- function()
{
    # Input data
    Awall  <- 200.0
    Afloor <- 50.0
    alpha  <- 2.0
    beta   <- 10
    gamma  <- 2
    delta  <- 10

    # Objective
    prob <- list(sense="max")
    prob$c <- c(1, 1, 1, 0, 0)

    # Linear constraints:
    # [ 0  0  0  1  1 ]                    == 1
    # [ 0  1  1  0  0 ]                    <= log(Afloor)     
    # [ 1 -1  0  0  0 ]                    in [log(alpha), log(beta)]
    # [ 0 -1  1  0  0 ]                    in [log(gamma), log(delta)]
    #
    prob$A <- rbind(c(0, 0, 0, 1, 1),
                    c(0, 1, 1, 0, 0),
                    c(1,-1, 0, 0, 0),
                    c(0,-1, 1, 0, 0))
               
    prob$bc <- rbind(c(1, -Inf,        log(alpha), log(gamma)),
                     c(1, log(Afloor), log(beta),  log(delta)))

    prob$bx <- rbind(c(-Inf, -Inf, -Inf, -Inf, -Inf),
                     c( Inf,  Inf,  Inf,  Inf,  Inf))

    # The conic part FX+g \in Kexp x Kexp
    #   x  y  z  u  v
    # [ 0  0  0  1  0 ]    0 
    # [ 0  0  0  0  0 ]    1               in Kexp
    # [ 1  1  0  0  0 ]    log(2/Awall)
    #
    # [ 0  0  0  0  1 ]    0
    # [ 0  0  0  0  0 ]    1               in Kexp
    # [ 1  0  1  0  0 ] +  log(2/Awall)
    #
    # NOTE: The F matrix is internally stored in the sparse
    #       triplet form. Use 'giveCsparse' or 'repr' option 
    #       in the sparseMatrix() call to construct the F 
    #       matrix directly in the sparse triplet form. 
    prob$F <- sparseMatrix(i = c(1, 3, 3, 4, 6, 6),
                           j = c(4, 1, 2, 5, 1, 3), 
                           x = c(1, 1, 1, 1, 1, 1),
                           dims = c(6,5))
                
    prob$g <- c(0, 1, log(2/Awall), 0, 1, log(2/Awall))

    prob$cones <- matrix(rep(list("PEXP", 3, NULL), 2), ncol=2)
    rownames(prob$cones) <- c("type","dim","conepar")

    # Optimize and print results
    r <- mosek(prob, list(soldetail=1))
    print(exp(r$sol$itr$xx[1:3]))
}

gp1()

helloworld.R

Listing 15.6 helloworld.R Click here to download.
##
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      helloworld.R
#
#  The most basic example of how to get started with MOSEK.

library("Rmosek")

prob <- list(sense="min")          # Minimization problem
prob$A <- Matrix(nrow=0, ncol=1)   # 0 constraints, 1 variable
prob$bx <- rbind(blx=2.0, bux=3.0) # Bounds on the only variable
prob$c <- c(1.0)                   # The objective coefficient

# Optimize
r <- mosek(prob)

# Print answer
r$sol$itr$xx

lo1.R

Listing 15.7 lo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      lo1.R
#
#  Purpose :   To demonstrate how to solve a small linear
#              optimization problem using the Rmosek package.
##
library("Rmosek")

lo1 <- function()
{
    prob <- list()

    # Objective sense (maximize or minimize)
    prob$sense <- "max"

    # Objective coefficients
    prob$c <- c(3, 1, 5, 1)

    # Specify matrix 'A' in sparse format.
    asubi  <- c(1, 1, 1, 2, 2, 2, 2, 3, 3)
    asubj  <- c(1, 2, 3, 1, 2, 3, 4, 2, 4)
    aval   <- c(3, 1, 2, 2, 1, 3, 1, 2, 3)

    prob$A <- sparseMatrix(asubi,asubj,x=aval)

    # Bound values for constraints
    prob$bc <- rbind(blc=c(30,  15, -Inf), 
                     buc=c(30, Inf,   25))

    # Bound values for variables
    prob$bx <- rbind(blx=rep(0,4), 
                     bux=c(Inf, 10, Inf, Inf))

    # Solve the problem
    r <- mosek(prob)
    
    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

lo1()

logistic.R

Listing 15.8 logistic.R Click here to download.
##
#  File : logistic.R
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  Description : Implements logistic regression with regulatization.
#
#          Demonstrates using the exponential cone and log-sum-exp.
#
#          Plots an example for 2D datasets.
# 
##
library("Rmosek")

logisticRegression <- function(X, y, lamb)
{
    prob <- list(sense="min")
    n <- dim(X)[1];
    d <- dim(X)[2];
    
    # Variables: r, theta(d), t(n), z1(n), z2(n)
    prob$c <- c(lamb, rep(0,d), rep(1, n), rep(0,n), rep(0,n));
    prob$bx <-rbind(rep(-Inf,1+d+3*n), rep(Inf,1+d+3*n));

    # z1 + z2 <= 1
    prob$A <- sparseMatrix( rep(1:n, 2), 
                            c((1:n)+1+d+n, (1:n)+1+d+2*n),
                            x = rep(1, 2*n));
    prob$bc <- rbind(rep(-Inf, n), rep(1, n));

    # (r, theta) \in \Q
    FQ <- cbind(diag(rep(1, d+1)), matrix(0, d+1, 3*n));
    gQ <- rep(0, 1+d);

    # (z1(i), 1, -t(i)) \in \EXP, 
    # (z2(i), 1, (1-2y(i))*X(i,) - t(i)) \in \EXP
    FE <- Matrix(nrow=0, ncol = 1+d+3*n);
    for(i in 1:n) {
        FE <- rbind(FE,
                    sparseMatrix( c(1, 3, 4, rep(6, d), 6),
                                  c(1+d+n+i, 1+d+i, 1+d+2*n+i, 2:(d+1), 1+d+i),
                                  x = c(1, -1, 1, (1-2*y[i])*X[i,], -1),
                                  dims = c(6, 1+d+3*n) ) );
    }
    gE <- rep(c(0, 1, 0, 0, 1, 0), n);

    prob$F <- rbind(FQ, FE)
    prob$g <- c(gQ, gE)
    prob$cones <- cbind(matrix(list("QUAD", 1+d, NULL), nrow=3, ncol=1),
                        matrix(list("PEXP", 3, NULL), nrow=3, ncol=2*n));
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve, no error handling!
    r <- mosek(prob, list(soldetail=1))

    # Return theta
    r$sol$itr$xx[2:(d+1)]
}

mapFeature <- function(x, d)
{
    g = as.matrix(expand.grid(0:d,0:d));
    g = g[g[,1]+g[,2]<=d,];
    apply(g, 1, function(gel) x[1]^gel[1]*x[2]^gel[2]);
}

mapFeatures <- function(X, d)
{
    t(apply(X, 1, function(x) mapFeature(x, d)));
}

runExample <-function(x1, x2, y, deg, lamb, pngname)
{
    X = mapFeatures(cbind(x1,x2),deg);
    theta = logisticRegression(X, y, lamb);
    print(theta);

    xplot <- seq(-1, 1, 0.05);
    yplot <- seq(-1, 1, 0.05);
    nplot <- length(xplot);
    zplot <- matrix(mapFeatures(as.matrix(expand.grid(xplot, yplot)),deg) %*% theta, ncol=nplot, nrow=nplot);

    # Uncomment to plot
    #png(pngname);
    #plot(x1, x2, pch=19, col=y+1);
    #contour(xplot, xplot, zplot, levels=c(0), add=TRUE);    
}

example1 <- function()
{
    # Example from documentation is contained in
    # https://datascienceplus.com/wp-content/uploads/2017/02/ex2data2.txt
    ex <- read.csv("ex2data2.txt");
    x1 <- ex[,1];
    x2 <- ex[,2];
    y <- ex[,3];
    deg <- 6;

    for(lamb in c(1e-0,1e-2,1e-4))
    {
        runExample(x1, x2, y, deg, lamb, sprintf("ex2-%f.png", lamb));
    }
}

example2 <- function()
{
    # Example 2: discover a circle
    n <- 300;
    deg <- 2;
    x1 <- runif(n, -1, 1);
    x2 <- runif(n, -1, 1);
    y  <- ifelse(x1^2 + x2^2 - 0.69 < 0, 0, 1);
    x1 <- x1 + runif(n, -0.2, 0.2);
    x2 <- x2 + runif(n, -0.2, 0.2);

    runExample(x1, x2, y, deg, 0.1, "circle.png");
}

#example1();
example2();

mico1.R

Listing 15.9 mico1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      mico1.R
#
#  Purpose :   To demonstrate how to solve a small mixed integer conic
#              optimization problem using the Rmosek package.
#
#              minimize    x^2 + y^2
#              subject to  x >= e^y + 3.8
#                          x, y - integer
##
library("Rmosek")

mico1 <- function()
{
    # Specify the continuous part of the problem.
    # Variables are [t; x; y]
    prob <- list(sense="min")
    prob$c <- c(1, 0, 0)
    prob$A <- Matrix(nrow=0, ncol=3)   # 0 constraints, 3 variables
    prob$bx <- rbind(blx=rep(-Inf,3), bux=rep(Inf,3))

    # Conic part of the problem
    prob$F <- rbind(diag(3), c(0,1,0), c(0,0,0), c(0,0,1))
    prob$g <- c(0, 0, 0, -3.8, 1, 0)
    prob$cones <- cbind(matrix(list("QUAD", 3, NULL), nrow=3, ncol=1),
                        matrix(list("PEXP", 3, NULL), nrow=3, ncol=1))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Specify the integer constraints
    prob$intsub <- c(2 ,3)

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    print(r$sol$int$xx[2:3])
}

mico1()

milo1.R

Listing 15.10 milo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      milo1.R
#
#  Purpose :   To demonstrate how to solve a small mixed integer linear
#              optimization problem using the Rmosek package.
##
library("Rmosek")

milo1 <- function()
{
    # Specify the continuous part of the problem.
    prob <- list(sense="max")
    prob$c  <- c(1, 0.64)
    prob$A  <- Matrix(rbind(c(50, 31),
                            c( 3, -2)), sparse=TRUE)
    prob$bc <- rbind(blc=c(-Inf,  -4),
                     buc=c( 250, Inf))
    prob$bx <- rbind(blx=c(  0,   0),
                     bux=c(Inf, Inf))

    # Specify the integer constraints
    prob$intsub <- c(1 ,2)

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

milo1()

mioinitsol.R

Listing 15.11 mioinitsol.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      mioinitsol.R
#
#  Purpose :   To demonstrate how to solve a linear mixed-integer
#              problem with a start guess.
##
library("Rmosek")

mioinitsol <- function()
{
    # Specify the problem.
    prob <- list(sense="max")
    prob$c  <- c(7, 10, 1, 5)
    prob$A  <- Matrix(c(1, 1, 1, 1), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=-Inf,
                     buc=2.5)
    prob$bx <- rbind(blx=rep(0,4),
                     bux=rep(Inf,4))
    prob$intsub <- c(1 ,2, 3)

    # Specify start guess for the integer variables.
    prob$sol$int$xx <- c(1, 1, 0, NaN)

    # Request constructing the solution from integer variable values
    prob$iparam <- list(MIO_CONSTRUCT_SOL=1)

    # Solve the problem
    r <- mosek(prob, list(getinfo=TRUE))

    # The solution
    stopifnot(identical(r$response$code, 0))
    print(r$sol$int$xx)

    # Was the initial solution used ?
    print(r$iinfo$MIO_CONSTRUCT_SOLUTION)
    print(r$dinfo$MIO_CONSTRUCT_SOLUTION_OBJ)
}

mioinitsol()

normex.R

Listing 15.12 normex.R Click here to download.
##
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      normex.R
#
#  Purpose:   Demonstrates various norm minimization problems
#
#             * least squares regression
#             * ridge regularization
#             * lasso regularization
#             * p-norm minimization
#
##
library("Rmosek")

# This function runs a few examples
normex <- function()
{
    # Create some random data
    # We will attempt to fit various polynomials to
    # a noisy degree 3 polynomial p(t)
    N  <- 100
    x1 <- sort(runif(N))
    p  <- function(t) { 1.5*t^3-2*t^2+0.5*t^1-1 }
    y1 <- p(x1) + 0.01*runif(N);

    # Least squares regression
    deg<- 5
    xp <- outer(x1, 0:deg, '^')
    x <- norm_lse(xp, y1, matrix(0, 0, deg+1), numeric(0))
    print(sprintf("LSE: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))

    # With ridge regularization
    x <- norm_lse_reg(xp, y1, matrix(0, 0, deg+1), numeric(0), 0.2)
    print(sprintf("Ridge: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))

    # With squared version of ridge regularization
    x <- norm_lse_reg_quad(xp, y1, matrix(0, 0, deg+1), numeric(0), 0.2)
    print(sprintf("QRidge: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", x[deg+1], deg, x[deg], deg-1, x[1]))

    # Completely random large data for lasso example
    N <- 10
    K <- 3000
    F <- matrix(rnorm(K*N), K, N)
    g <- runif(K)
    print(sprintf("Lasso regularization"))
    for (gamma in c(0.01, 0.1, 0.3, 0.6, 0.9, 1.3))
    {
        x <- norm_lse_lasso(F, g, matrix(0,0,N), numeric(0), gamma)
        z <- data.frame(x)
        print(sprintf('Gamma %.4f  density %.0f  |Fx-g|_2: %.4f', gamma, sum(abs(z)>1e-6)/N*100, norm(F%*%x-g, type="2")))
    }

    # Example with the p-norm cone for various p
    # We add a far outlier to the first example
    x12 <- c(x1, 0.73)
    y12 <- c(y1, -0.99)
    xp2 <- outer(x12, 0:deg, '^')
    for (p in c(1.1, 2.0, 3.0, 6.0))
    {
       x <- norm_p_norm(xp2, y12, matrix(0, 0, deg+1), numeric(0), p)
       print(sprintf("p = %.1f: best fit %.2f x^%d + %.2f x^%d + ... + %.2f", p, x[deg+1], deg, x[deg], deg-1, x[1]))
    }
}

# Least squares regression
# minimize \|Fx-g\|_2
norm_lse <- function(F,g,A,b)
{
    n <- dim(F)[2]
    k <- length(g)
    m <- dim(A)[1]

    # Linear constraints in [x; t]
    prob <- list(sense="min")
    prob$A <- cbind(A, rep(0, m))
    prob$bx <- rbind(rep(-Inf, n+1), rep(Inf, n+1))
    prob$bc <- rbind(b, b)
    prob$c <- c(rep(0, n), 1)

    # Affine conic constraint
    prob$F <- rbind( c(rep(0,n), 1),
                     cbind(F, rep(0, k)) )
    prob$g <- c(0, -g)
    prob$cones <- matrix(list("QUAD", k+1, NULL))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$response$code, 0))
    r$sol$itr$xx[1:n]
}

# Least squares regression with regularization
# minimize \|Fx-g\|_2 + gamma*\|x\|_2
norm_lse_reg <- function(F,g,A,b,gamma)
{
    n <- dim(F)[2]
    k <- length(g)
    m <- dim(A)[1]

    # Linear constraints in [x; t1; t2]
    prob <- list(sense="min")
    prob$A <- cbind(A, rep(0, m), rep(0, m))
    prob$bx <- rbind(rep(-Inf, n+2), rep(Inf, n+2))
    prob$bc <- rbind(b, b)
    prob$c <- c(rep(0, n), 1, gamma)

    # Affine conic constraint
    prob$F <- rbind(     c(rep(0,n), 1,         0),
                     cbind(F,        rep(0, k), rep(0, k)),
                         c(rep(0,n), 0,         1),
                     cbind(diag(n),  rep(0, n), rep(0, n)))
    prob$g <- c(0, -g, rep(0, n+1))
    prob$cones <- cbind(matrix(list("QUAD", k+1, NULL)), 
                        matrix(list("QUAD", n+1, NULL)))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$response$code, 0))
    r$sol$itr$xx[1:n]
}

# Least squares regression with regularization
# The "classical" quadratic version
# minimize \|Fx-g\|_2^2 + gamma*\|x\|_2^2
norm_lse_reg_quad <- function(F,g,A,b,gamma)
{
    n <- dim(F)[2]
    k <- length(g)
    m <- dim(A)[1]

    # Linear constraints in [x; t1; t2]
    prob <- list(sense="min")
    prob$A <- cbind(A, rep(0, m), rep(0, m))
    prob$bx <- rbind(rep(-Inf, n+2), rep(Inf, n+2))
    prob$bc <- rbind(b, b)
    prob$c <- c(rep(0, n), 1, gamma)

    # Affine conic constraint
    prob$F <- rbind(     c(rep(0,n), 1,         0),
                         c(rep(0, n+2)),
                     cbind(F,        rep(0, k), rep(0, k)),
                         c(rep(0,n), 0,         1),
                         c(rep(0, n+2)),
                     cbind(diag(n),  rep(0, n), rep(0, n)))
    prob$g <- c(0, 0.5, -g, 0, 0.5, rep(0, n))
    prob$cones <- cbind(matrix(list("RQUAD", k+2, NULL)), 
                        matrix(list("RQUAD", n+2, NULL)))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$response$code, 0))
    r$sol$itr$xx[1:n]
}

# Least squares regression with lasso regularization
# minimize \|Fx-g\|_2 + gamma*\|x\|_1
norm_lse_lasso <- function(F,g,A,b,gamma)
{
    n <- dim(F)[2]
    k <- length(g)
    m <- dim(A)[1]

    # Linear constraints in [x; u; t1; t2]
    prob <- list(sense="min")
    prob$A <- rbind(cbind(A,    matrix(0, m, n+2)),
                    cbind(diag(n), diag(n), matrix(0, n, 2)),
                    cbind(-diag(n), diag(n), matrix(0, n, 2)),
                    c(rep(0, n), rep(-1, n), 0, 1))
    prob$bx <- rbind(rep(-Inf, 2*n+2), rep(Inf, 2*n+2))
    prob$bc <- rbind(c(b, rep(0, 2*n+1)), c(b, rep(Inf, 2*n+1)))
    prob$c <- c(rep(0, 2*n), 1, gamma)

    # Affine conic constraint
    prob$F <- rbind(c(rep(0, 2*n), 1, 0),
                    cbind(F, matrix(0, k, n+2)))
    prob$g <- c(0, -g)
    prob$cones <- matrix(list("QUAD", k+1, NULL))
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$response$code, 0))
    r$sol$itr$xx[1:n]
}

# P-norm minimization
# minimize \|Fx-g\|_p
norm_p_norm <- function(F,g,A,b,p)
{
    n <- dim(F)[2]
    k <- length(g)
    m <- dim(A)[1]

    # Linear constraints in [x; r; t]
    prob <- list(sense="min")
    prob$A <- rbind(cbind(A,    matrix(0, m, k+1)),
                    c(rep(0, n), rep(1, k), -1))
    prob$bx <- rbind(rep(-Inf, n+k+1), rep(Inf, n+k+1))
    prob$bc <- rbind(c(b, 0), c(b, 0))
    prob$c <- c(rep(0, n+k), 1)

    # Permutation matrix which picks triples (r_i, t, F_ix-g_i)
    M <- cbind(sparseMatrix(seq(1,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)),
               sparseMatrix(seq(2,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)),
               sparseMatrix(seq(3,3*k,3), seq(1,k), x=rep(1,k), dims=c(3*k,k)))

    # Affine conic constraint
    prob$F <- M %*% rbind(cbind(matrix(0, k, n), diag(k), rep(0, k)),
                          cbind(matrix(0, k, n+k), rep(1, k)),
                          cbind(F, matrix(0, k, k+1)))
    prob$g <- as.numeric(M %*% c(rep(0, 2*k), -g))
    prob$cones <- matrix(list("PPOW", 3, c(1, p-1)), nrow=3, ncol=k)
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$sol$itr$prosta, "PRIMAL_AND_DUAL_FEASIBLE"))
    r$sol$itr$xx[1:n]
}

normex()

parameters.R

Listing 15.13 parameters.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      parameters.R
#
#  Purpose :   To demonstrate how to set parameters 
#              using the Rmosek package.
##
library("Rmosek")

parameters <- function()
{
    # Specify the optimization problem
    prob <- list(sense="max")
    prob$c <- c(3, 1, 5, 1)
    prob$A <- sparseMatrix(i=c(1, 1, 1, 2, 2, 2, 2, 3, 3),
                           j=c(1, 2, 3, 1, 2, 3, 4, 2, 4),
                           x=c(3, 1, 2, 2, 1, 3, 1, 2, 3))
    prob$bc <- rbind(blc=c(30,  15, -Inf), 
                     buc=c(30, Inf,   25))
    prob$bx <- rbind(blx=rep(0,4), 
                     bux=c(Inf, 10, Inf, Inf))

    # Specify the parameters

    # Set log level (integer parameter)
    # and select interior-point optimizer... (integer parameter)    
    # ... without basis identification (integer parameter)    
    prob$iparam <- list(LOG=0, OPTIMIZER="OPTIMIZER_INTPNT", INTPNT_BASIS="BI_NEVER")

    # Set relative gap tolerance (double parameter)
    prob$dparam <- list(INTPNT_CO_TOL_REL_GAP=1.0e-7)    

    # Solve the problem
    r <- mosek(prob)

    # Demonstrate retrieving information items

    opts <- list(getinfo=TRUE)
    r <- mosek(prob, opts)
    
    print(r$dinfo$OPTIMIZER_TIME)
    print(r$iinfo$INTPNT_ITER)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

parameters()

portfolio_1_basic.R

Listing 15.14 portfolio_1_basic.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_1_basic.R
#
#  Purpose :   To implement a basic Markowitz optimization model computing 
#              the optimal expected return and portfolio for a given risk.
##
library("Rmosek")

BasicMarkowitz <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    gamma)      # Maximum risk (=std. dev) accepted
{
    prob <- list(sense="max")
    prob$c <- mu
    prob$A <- Matrix(1.0, ncol=n)
    prob$bc <- rbind(blc=w+sum(x0), 
                     buc=w+sum(x0))
    prob$bx <- rbind(blx=rep(0.0,n),
                     bux=rep(Inf,n))

    # Specify the affine conic constraints.
    NUMCONES <- 1
    prob$F <- rbind(
        Matrix(0.0,ncol=n), 
        GT
    )
    prob$g <- c(gamma,rep(0,n))
    prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
    rownames(prob$cones) <- c("type","dim","conepar")

    prob$cones[-3,1] <- list("QUAD", n+1)

    # Solve the problem
    r <- mosek(prob,list(verbose=1))
    stopifnot(identical(r$response$code, 0))

    # Return the solution
    x <- r$sol$itr$xx
    list(expret=drop(mu %*% x), stddev=gamma, x=x)
}

# Example of input
n      <- 8
w      <- 59.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0)
gamma  <- 36.0
GT     <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
                 c(0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
                 c(0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
                 c(0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
                 c(0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ) )

r <- BasicMarkowitz(n,mu,GT,x0,w,gamma)
print(sprintf('Expected return: %.4e   Std. deviation: %.4e', r$expret, r$stddev))

portfolio_2_frontier.R

Listing 15.15 portfolio_2_frontier.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_2_frontier.R
#
#  Purpose :   To implements a basic portfolio optimization model.
#              Computes points on the efficient frontier.
##
library("Rmosek")

EfficientFrontier <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    alphas)     # List of risk penalties (we maximize expected return - alpha * variance)
{
    prob <- list(sense="max")
    prob$A <- cbind(Matrix(1.0, ncol=n), 0.0)
    prob$bc <- rbind(blc=w+sum(x0), 
                     buc=w+sum(x0))
    prob$bx <- rbind(blx=c(rep(0.0,n), -Inf),
                     bux=rep(Inf,n+1))

    # Specify the affine conic constraints.
    NUMCONES <- 1
    prob$F <- rbind(
        cbind(Matrix(0.0,ncol=n), 1.0),
        rep(0, n+1),
        cbind(GT                , 0.0)
    )
    prob$g <- c(0, 0.5, rep(0, n))
    prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
    rownames(prob$cones) <- c("type","dim","conepar")

    prob$cones[-3,1] <- list("RQUAD", n+2)

    frontier <- matrix(NaN, ncol=3, nrow=length(alphas))
    colnames(frontier) <- c("alpha","exp.ret.","variance")

    for (i in seq_along(alphas))
    {
        prob$c <- c(mu, -alphas[i])

        r <- mosek(prob, list(verbose=1))
        stopifnot(identical(r$response$code, 0))

        x     <- r$sol$itr$xx[1:n]
        gamma <- r$sol$itr$xx[n+1]
        
        frontier[i,] <- c(alphas[i], drop(mu%*%x), gamma)
    }

    frontier
}

# Example of input
n      <- 8
w      <- 1.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT     <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
                 c(0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
                 c(0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
                 c(0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
                 c(0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ) )
alphas <- c(0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0)

frontier <- EfficientFrontier(n,mu,GT,x0,w,alphas)
print(frontier)

portfolio_3_impact.R

Listing 15.16 portfolio_3_impact.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_3_impact.R
#
#  Purpose :   Implements a basic portfolio optimization model
#              with x^(3/2) market impact costs.
##
library("Rmosek")

MarkowitzWithMarketImpact <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    gamma,      # Maximum risk (=std. dev) accepted
    m)          # Market impacts (we use m_j|x_j-x0_j|^3/2 for j'th asset)
{
    prob <- list(sense="max")
    prob$c <- c(mu, rep(0,n))
    prob$A <- cbind(Matrix(1.0,ncol=n), t(m))
    prob$bc <- rbind(blc=w+sum(x0),
                     buc=w+sum(x0))
    prob$bx <- rbind(blx=rep(0.0,2*n),
                     bux=rep(Inf,2*n))

    # Specify the affine conic constraints.
    # 1) Risk
    Fr <- rbind(
        Matrix(0.0,nrow=1,ncol=2*n), 
        cbind(GT, Matrix(0.0,nrow=n,ncol=n))
    )
    gr <- c(gamma,rep(0,n))
    Kr <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)

    # 2) Market impact (t_j >= |x_j-x0_j|^3/2)
    # [    t_j     ]
    # [     1      ] \in PPOW(2,1)
    # [ x_j - x0_j ]
    Fm <- sparseMatrix(
                 i=c(seq(from=1,by=3,len=n), seq(from=3,by=3,len=n)),
                 j=c(seq(from=n+1,len=n),    seq(from=1,len=n)),
                 x=c(rep(1.0,n),             rep(1.0,n)),
                 dims=c(3*n, 2*n))
    gm <- rep(c(0,1,0), n)
    gm[seq(from=3,by=3,len=n)] <- -x0
    Km <- matrix(list("PPOW", 3, c(2,1)), nrow=3, ncol=n)
    
    prob$F <- rbind(Fr, Fm)
    prob$g <- c(gr, gm)
    prob$cones <- cbind(Kr, Km)
    rownames(prob$cones) <- c("type","dim","conepar")

    # Solve the problem
    r <- mosek(prob,list(verbose=1))
    stopifnot(identical(r$response$code, 0))

    # Return the solution
    x <- r$sol$itr$xx[1:n]
    tx <- r$sol$itr$xx[(n+1):(2*n)]
    list(expret=drop(mu %*% x), stddev=gamma, cost=drop(m %*% tx), x=x)
}

# Example of input
n      <- 8
w      <- 1.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT     <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
                 c(0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
                 c(0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
                 c(0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
                 c(0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ) )
gamma  <- 0.36
m      <- rep(0.01, n)

r <- MarkowitzWithMarketImpact(n,mu,GT,x0,w,gamma,m)
print(sprintf('Expected return: %.4e   Std. deviation: %.4e  Market impact cost: %.4e', r$expret, r$stddev, r$cost))

portfolio_4_transcost.R

Listing 15.17 portfolio_4_transcost.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_4_transcost.R
#
#  Purpose :   Implements a basic portfolio optimization model
#              with fixed setup costs and transaction costs
#              as a mixed-integer problem.
##
library("Rmosek")

MarkowitzWithTransactionCosts <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    gamma,      # Maximum risk (=std. dev) accepted
    f,          # Fixed transaction cost
    g)          # Linear part of transaction cost
{

    # Upper bound on the traded amount
    u <- w+sum(x0)

    prob <- list(sense="max")
    prob$c <- c(mu, rep(0,2*n))

    # Specify linear constraints
    # [ e'  g'  f' ]   [ x ]  =   w + e'*x0
    # [ I  -I   0  ] * [ z ]  <=  x0
    # [ I   I   0  ]   [ y ]  >=  x0
    # [ 0   I  -U  ]          <=  0
    prob$A <- rbind(cbind(Matrix(1.0,ncol=n), t(g),              t(f)),
                    cbind(Diagonal(n, 1.0),   -Diagonal(n, 1.0), Matrix(0,n,n)),
                    cbind(Diagonal(n, 1.0),   Diagonal(n, 1.0),  Matrix(0,n,n)),
                    cbind(Matrix(0,n,n),      Diagonal(n, 1.0),  Diagonal(n, -u)))
    prob$bc <- rbind(blc=c(w+sum(x0), rep(-Inf,n), x0, rep(-Inf,n)),
                     buc=c(w+sum(x0), x0, rep(Inf,n), rep(0.0,n)))
    # No shortselling and the linear bound 0 <= y <= 1     
    prob$bx <- rbind(blx=c(rep(0.0,n), rep(-Inf,n), rep(0.0,n)),
                     bux=c(rep(Inf,n), rep(Inf, n), rep(1.0,n)))

    # Specify the affine conic constraints for risk
    prob$F <- rbind(
        Matrix(0.0,nrow=1,ncol=3*n), 
        cbind(GT, Matrix(0.0,nrow=n,ncol=2*n))
    )
    prob$g <- c(gamma,rep(0,n))
    prob$cones <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
    rownames(prob$cones) <- c("type","dim","conepar")

    # Demand y to be integer (hence binary)
    prob$intsub <- (2*n+1):(3*n);

    # Solve the problem
    r <- mosek(prob,list(verbose=1))
    stopifnot(identical(r$response$code, 0))

    # Return the solution
    x <- r$sol$int$xx[1:n]
    z <- r$sol$int$xx[(n+1):(2*n)]
    y <- r$sol$int$xx[(2*n+1):(3*n)]
    list(expret=drop(mu %*% x), stddev=gamma, cost = drop(f %*% y)+drop(g %*% z), x=x)
}

# Example of input
n      <- 8
w      <- 1.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT     <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
                 c(0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
                 c(0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
                 c(0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
                 c(0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ) )
gamma  <- 0.36
f      <- rep(0.01,n)
g      <- rep(0.001,n)

r <- MarkowitzWithTransactionCosts(n,mu,GT,x0,w,gamma,f,g)
print(sprintf("Markowitz portfolio optimization with transactions cost"))
print(sprintf("Expected return: %.4e Std. deviation: %.4e Transactions cost: %.4e", r$expret, r$stddev, r$cost))

portfolio_5_card.R

Listing 15.18 portfolio_5_card.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_5_card.R
#
#  Description :  Implements a basic portfolio optimization model
#                 with cardinality constraints on number of assets traded.
##
library("Rmosek")

MarkowitzWithCardinality <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    gamma,      # Maximum risk (=std. dev) accepted
    k)          # Cardinality bound
{

    # Upper bound on the traded amount
    u <- w+sum(x0)

    prob <- list(sense="max")
    prob$c <- c(mu, rep(0,2*n))

    # Specify linear constraints
    # [ e'  0   0  ]           =   w + e'*x0
    # [ I  -I   0  ]   [ x ]  <=  x0
    # [ I   I   0  ] * [ z ]  >=  x0
    # [ 0   I  -U  ]   [ y ]  <=  0
    # [ 0   0   e' ]          <=  k
    prob$A <- rbind(cbind(Matrix(1.0,ncol=n), Matrix(0.0,ncol=2*n)),
                    cbind(Diagonal(n, 1.0),   -Diagonal(n, 1.0), Matrix(0,n,n)),
                    cbind(Diagonal(n, 1.0),   Diagonal(n, 1.0),  Matrix(0,n,n)),
                    cbind(Matrix(0,n,n),      Diagonal(n, 1.0),  Diagonal(n, -u)),
                    cbind(Matrix(0.0,ncol=2*n), Matrix(1.0,ncol=n)))
    prob$bc <- rbind(blc=c(w+sum(x0), rep(-Inf,n), x0, rep(-Inf,n), 0.0),
                     buc=c(w+sum(x0), x0, rep(Inf,n), rep(0.0,n), k))
    # No shortselling and the linear bound 0 <= y <= 1     
    prob$bx <- rbind(blx=c(rep(0.0,n), rep(-Inf,n), rep(0.0,n)),
                     bux=c(rep(Inf,n), rep(Inf, n), rep(1.0,n)))

    # Specify the affine conic constraints for risk
    prob$F <- rbind(
        Matrix(0.0,nrow=1,ncol=3*n), 
        cbind(GT, Matrix(0.0,nrow=n,ncol=2*n))
    )
    prob$g <- c(gamma,rep(0,n))
    prob$cones <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
    rownames(prob$cones) <- c("type","dim","conepar")

    # Demand y to be integer (hence binary)
    prob$intsub <- (2*n+1):(3*n);

    # Solve the problem
    r <- mosek(prob,list(verbose=1))
    stopifnot(identical(r$response$code, 0))

    # Return the solution
    x <- r$sol$int$xx[1:n]
    list(card=k, expret=drop(mu %*% x), x=x)
}

# Example of input
n      <- 8
w      <- 1.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
GT     <- rbind( c(0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638),
                 c(0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506),
                 c(0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914),
                 c(0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742),
                 c(0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 ),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327),
                 c(0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 ) )
gamma  <- 0.25

print(sprintf("Markowitz portfolio optimization with cardinality constraints"))
for (k in 1:n) {
    r <- MarkowitzWithCardinality(n,mu,GT,x0,w,gamma,k)
    print(sprintf("Bound: %d   Expected return: %.4e  Solution:", r$card, r$expret))
    print(r$x)

}

portfolio_6_factor.R

Listing 15.19 portfolio_6_factor.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      portfolio_6_factor.R
#
#  Purpose :   Implements a basic Markowitz portfolio model of factor type
##
library("Rmosek")

BasicMarkowitz <- function(
    n,          # Number of assets
    mu,         # An n-dimmensional vector of expected returns
    GT,         # A matrix with n columns so (GT')*GT  = covariance matrix
    x0,         # Initial holdings 
    w,          # Initial cash holding
    gamma)      # Maximum risk (=std. dev) accepted
{
    k <- dim(GT)[1]
    prob <- list(sense="max")
    prob$c <- mu
    prob$A <- Matrix(1.0, ncol=n)
    prob$bc <- rbind(blc=w+sum(x0), 
                     buc=w+sum(x0))
    prob$bx <- rbind(blx=rep(0.0,n),
                     bux=rep(Inf,n))

    # Specify the affine conic constraints.
    NUMCONES <- 1
    prob$F <- rbind(
        Matrix(0.0,ncol=n), 
        GT
    )
    prob$g <- c(gamma,rep(0,k))
    prob$cones <- matrix(list(), nrow=3, ncol=NUMCONES)
    rownames(prob$cones) <- c("type","dim","conepar")

    prob$cones[-3,1] <- list("QUAD", k+1)

    # Solve the problem
    r <- mosek(prob,list(verbose=1))
    stopifnot(identical(r$response$code, 0))

    # Return the solution
    x <- r$sol$itr$xx
    list(expret=drop(mu %*% x), stddev=gamma, x=x)
}

# Example of input
n      <- 8
w      <- 1.0
mu     <- c(0.07197349, 0.15518171, 0.17535435, 0.0898094, 0.42895777, 0.39291844, 0.32170722, 0.18378628)
x0     <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
gammas <- c(0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48)

B    <- rbind(c(0.4256, 0.1869),
              c(0.2413, 0.3877),
              c(0.2235, 0.3697),
              c(0.1503, 0.4612),
              c(1.5325, -0.2633),
              c(1.2741, -0.2613),
              c(0.6939, 0.2372),
              c(0.5425, 0.2116) )

S_F <- rbind(c(0.0620, 0.0577),
             c(0.0577, 0.0908))

S_theta <- diag(c(0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459))

# Compute the small factorization required for the model
P  <- t(chol(S_F))
G  <- cbind(B %*% P,   sqrt(S_theta))
GT <- t(G)

for (gamma in gammas)
{
    # We call the same model as in portfolio_1_basic.R
    # but now we use the GT matrix obtained from the factor model.
    r <- BasicMarkowitz(n,mu,GT,x0,w,gamma)
    print(sprintf('Gamma: %.4e  Expected return: %.4e   Std. deviation: %.4e', gamma, r$expret, r$stddev))
}

pow1.R

Listing 15.20 pow1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      pow1.R
#
#  Purpose :   To demonstrate how to solve a small power cone
#              optimization problem using the Rmosek package.
##
library("Rmosek")

pow1 <- function()
{
    # Specify the non-conic part of the problem.
    prob <- list(sense="max")
    prob$c  <- c(-1, 0, 0, 1, 1)
    prob$A  <- Matrix(c(1, 1, 0.5, 0, 0), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=2, 
                     buc=2)
    prob$bx <- rbind(blx=c(rep(-Inf,5)), 
                     bux=c(rep( Inf,5)))
    
    # Specify the affine conic constraints.
    # NOTE: The F matrix is internally stored in the sparse
    #       triplet form. Use 'giveCsparse' or 'repr' option 
    #       in the sparseMatrix() call to construct the F 
    #       matrix directly in the sparse triplet form. 
    prob$F     <- sparseMatrix(i=c(1, 2, 3, 4, 6),
                               j=c(1, 2, 4, 3, 5), 
                               x=c(1, 1, 1, 1, 1), 
                               dims = c(6,5))
    prob$g     <- c(0, 0, 0, 0, 1, 0)
    prob$cones <- matrix(list(), nrow=3, ncol=2)
    rownames(prob$cones) <- c("type","dim","conepar")

    prob$cones[,1] <- list("PPOW", 3, c(0.2, 0.8))
    prob$cones[,2] <- list("PPOW", 3, c(0.4, 0.6))

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

pow1()

qo1.R

Listing 15.21 qo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      qo1.R
#
#  Purpose :   To demonstrate how to solve a small quadratic
#              optimization problem using the Rmosek package.
##
library("Rmosek")

qo1 <- function()
{
    # Specify the non-quadratic part of the problem.
    prob <- list(sense="min")
    prob$c <- c(0, -1, 0)
    prob$A <- Matrix(c(1, 1, 1), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=1, 
                     buc=Inf)
    prob$bx <- rbind(blx=rep(0,3), 
                     bux=rep(Inf,3))

    # Specify the quadratic objective matrix in triplet form.
    prob$qobj$i <- c(1,  3,   2,  3)
    prob$qobj$j <- c(1,  1,   2,  3)
    prob$qobj$v <- c(2, -1, 0.2,  2)

    # Solve the problem
    r <- mosek(prob)

    # Return the solution
    stopifnot(identical(r$response$code, 0))
    r$sol
}

qo1()

reoptimization.R

Listing 15.22 reoptimization.R Click here to download.
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      reoptimization.R
#
#  Purpose : Demonstrates how to modify and re-optimize a linear problem
#
library("Rmosek")

reoptimization <- function()
{
    # Specify the c vector.
    prob   <- list(sense="max")
    prob$c <- c(1.5, 2.5, 3.0)

    # Specify a in sparse format.
    subi   <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
    subj   <- c(1, 2, 3, 1, 2, 3, 1, 2, 3)
    valij  <- c(2, 4, 3, 3, 2, 3, 2, 3, 2)

    prob$A <- sparseMatrix(subi,subj,x=valij);

    # Specify bounds of the constraints.
    prob$bc<- rbind(blc=rep(-Inf, 3),
                    buc=c(100000, 50000, 60000))

    # Specify bounds of the variables.
    prob$bx<- rbind(blx=rep(0,3),
                    bux=rep(Inf,3))

    # Perform the optimization.
    prob$iparam <- list(OPTIMIZER="OPTIMIZER_DUAL_SIMPLEX")
    r <- mosek(prob)

    # Show the optimal x solution.
    print(r$sol$bas$xx)

    #####################################################
    # Change a coefficient
    prob$A[1,1] <- 3.0

    # Reoptimize with changed coefficient
    # Use previous solution to perform very simple hot-start.
    # This part can be skipped, but then the optimizer will start
    # from scratch on the new problem, i.e. without any hot-start.
    prob$sol <- list(bas=r$sol$bas)     
    r <- mosek(prob)
    print(r$sol$bas$xx)

    #####################################################
    # Add a variable
    prob$c       <- c(prob$c, 1.0)
    prob$A       <- cbind(prob$A, c(4.0, 0.0, 1.0))
    prob$bx      <- cbind(prob$bx, c(0.0,Inf))

    # Reoptimize with a new variable and hot-start
    # All parts of the solution must be extended to the new dimensions.
    prob$sol <- list(bas=r$sol$bas)     
    prob$sol$bas$xx  <- c(prob$sol$bas$xx, 0.0)
    prob$sol$bas$slx <- c(prob$sol$bas$slx, 0.0)
    prob$sol$bas$sux <- c(prob$sol$bas$sux, 0.0)
    prob$sol$bas$skx <- c(prob$sol$bas$skx, "UN")
    r <- mosek(prob)
    print(r$sol$bas$xx)

    #####################################################
    # Add a constraint
    prob$A       <- rbind(prob$A, c(1.0, 2.0, 1.0, 1.0)) 
    prob$bc      <- cbind(prob$bc, c(-Inf, 30000.0)) 

    # Reoptimization with a new constraint
    # All parts of the solution must be extended to the new dimensions.
    prob$sol <- list(bas=r$sol$bas)     
    prob$sol$bas$xc  <- c(prob$sol$bas$xc, 0.0)
    prob$sol$bas$slc <- c(prob$sol$bas$slc, 0.0)
    prob$sol$bas$suc <- c(prob$sol$bas$suc, 0.0)
    prob$sol$bas$skc <- c(prob$sol$bas$skc, "UN")
    r <- mosek(prob) 
    print(r$sol$bas$xx)

    #####################################################
    # Change constraint bounds
    prob$bc<- rbind(blc=rep(-Inf, 4),
                    buc=c(80000, 40000, 50000, 22000))
    r <- mosek(prob) 
    print(r$sol$bas$xx)

}

reoptimization()

response.R

Listing 15.23 response.R Click here to download.
##
#  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File:      response.R
#
#  Purpose:   This example demonstrates proper response handling
#             for problems solved with the interior-point optimizers.
#
library("Rmosek")

# Set up a simple conic problem for test purposes
setupProblem <- function()
{
    prob <- list(sense="min")
    prob$c  <- c(0, 0, 0, 1, 1, 1)
    prob$A  <- Matrix(c(1, 1, 2, 0, 0, 0), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=1, buc=1)
    prob$bx <- rbind(blx=c(rep(0,3), rep(-Inf,3)), bux=rep(Inf,6))

    prob$F  <- sparseMatrix(i=c(1, 2, 3, 4, 5, 6),
                            j=c(4, 1, 2, 5, 6, 3), 
                            x=c(1, 1, 1, 1, 1, 1), 
                            dims = c(6,6))
    prob$g  <- c(1:6)*0
    prob$cones <- matrix(list(), nrow=3, ncol=2)
    rownames(prob$cones) <- c("type","dim","conepar")
    prob$cones[,1] <- list("QUAD",  3, NULL)
    prob$cones[,2] <- list("RQUAD", 3, NULL)
    prob
}

response <- function()
{
    # In this example we set up a simple problem 
    prob <- setupProblem()

    # (Optionally) Uncomment the next line to get solution status Unknown
    # prob$iparam <- list(INTPNT_MAX_ITERATIONS=1)

    # Perform the optimization.
    r <- mosek(prob, list(verbose=0))
    # Use the line below instead to get more log output
    #r <- mosek(prob, list(verbose=10))

    # Expected result: The solution status of the interior-point solution is optimal.

    # Check if there was a fatal error
    if (r$response$code != 0 && startsWith(r$response$msg, "MSK_RES_ERR"))
    {
        print(sprintf("Optimization error: %s (%d),", r$response$msg, r$response$code))
    }
    else
    {
        if (r$sol$itr$solsta == 'OPTIMAL')
        {
            print('An optimal interior-point solution is located:')
            print(r$sol$itr$xx)
        }
        else if (r$sol$itr$solsta == 'DUAL_INFEASIBLE_CER')
        {
            print('Dual infeasibility certificate found.')
        }
        else if (r$sol$itr$solsta == 'PRIMAL_INFEASIBLE_CER')
        {
            print('Primal infeasibility certificate found.')
        }
        else if (r$sol$itr$solsta == 'UNKNOWN')
        { 
            # The solutions status is unknown. The termination code 
            # indicates why the optimizer terminated prematurely. 
            print('The solution status is unknown.')
            print(sprintf('Termination code: %s (%d).', r$response$msg, r$response$code))
        }
        else
        {
            printf('An unexpected solution status is obtained.')
        }
    }
}

response()

sdo1.R

Listing 15.24 sdo1.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo1.R
#
#  Purpose :   To solve the mixed semidefinite and conic quadratic optimization problem
#
#                 minimize    Tr [2, 1, 0; 1, 2, 1; 0, 1, 2]*X + x(1)
#
#                 subject to  Tr [1, 0, 0; 0, 1, 0; 0, 0, 1]*X + x(1)               = 1
#                             Tr [1, 1, 1; 1, 1, 1; 1, 1, 1]*X        + x(2) + x(3) = 0.5
#                             X>=0,  x(1) >= sqrt(x(2)^2 + x(3)^2)
##
library("Rmosek")

getbarvarMatrix <- function(barvar, bardim)
{
    N <- as.integer(bardim)
    new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}

sdo1 <- function()
{
    # Specify the non-matrix variable part of the problem.
    prob <- list(sense="min")
    prob$c     <- c(1, 0, 0)
    prob$A <- sparseMatrix(i=c(1, 2, 2),
                           j=c(1, 2, 3),
                           x=c(1, 1, 1), dims=c(2, 3))
    prob$bc    <- rbind(blc=c(1, 0.5), 
                        buc=c(1, 0.5))
    prob$bx    <- rbind(blx=rep(-Inf,3), 
                        bux=rep( Inf,3))

    # NOTE: The F matrix is internally stored in the sparse
    #       triplet form. Use 'giveCsparse' or 'repr' option 
    #       in the sparseMatrix() call to construct the F 
    #       matrix directly in the sparse triplet form. 
    prob$F     <- sparseMatrix(i=c(1,2,3), 
                               j=c(1,2,3), 
                               x=c(1,1,1), 
                               dims = c(3,3))
    prob$g     <- c(1:3)*0
    prob$cones <- cbind(list("QUAD", 3, NULL))

    # Specify semidefinite matrix variables (one 3x3 block)
    prob$bardim <- c(3)

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barc':
    prob$barc$j <- c(1, 1, 1, 1, 1)
    prob$barc$k <- c(1, 2, 3, 2, 3)
    prob$barc$l <- c(1, 2, 3, 1, 2)
    prob$barc$v <- c(2, 2, 2, 1, 1)

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barA':
    prob$barA$i <- c(1, 1, 1, 2, 2, 2, 2, 2, 2)
    prob$barA$j <- c(1, 1, 1, 1, 1, 1, 1, 1, 1)
    prob$barA$k <- c(1, 2, 3, 1, 2, 3, 2, 3, 3)
    prob$barA$l <- c(1, 2, 3, 1, 2, 3, 1, 1, 2)
    prob$barA$v <- c(1, 1, 1, 1, 1, 1, 1, 1, 1)

    # Solve the problem
    r <- mosek(prob)

    # Print matrix variable and return the solution
    stopifnot(identical(r$response$code, 0))
    print( list(barx=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1])) )
    r$sol
}

sdo1()

sdo2.R

Listing 15.25 sdo2.R Click here to download.
#
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo2.R
#
#  Purpose :   Solves the semidefinite problem with two symmetric variables:
#
#                 min   <C1,X1> + <C2,X2>
#                 st.   <A1,X1> + <A2,X2> = b
#                             (X2)_{1,2} <= k
#                
#                 where X1, X2 are symmetric positive semidefinite,
#
#                 C1, C2, A1, A2 are assumed to be constant symmetric matrices,
#                 and b, k are constants.
library("Rmosek")

getbarvarMatrix <- function(barvar, bardim)
{
    N <- as.integer(bardim)
    new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}

sdo2 <- function()
{
    # Sample data in sparse lower-triangular triplet form
    C1_k <- c(1, 3);
    C1_l <- c(1, 3);
    C1_v <- c(1, 6);
    A1_k <- c(1, 3, 3);
    A1_l <- c(1, 1, 3);
    A1_v <- c(1, 1, 2);
    C2_k <- c(1, 2, 2, 3);
    C2_l <- c(1, 1, 2, 3);
    C2_v <- c(1, -3, 2, 1);
    A2_k <- c(2, 2, 4);
    A2_l <- c(1, 2, 4);
    A2_v <- c(1, -1, -3);
    b <- 23.0;
    k <- -3.0;
 
    # Specify the dimensions of the problem
    prob <- list(sense="min");
    # Two constraints
    prob$A <- Matrix(nrow=2, ncol=0);                   # 2 constraints
    prob$c <- numeric(0);
    prob$bx <- rbind(blc=numeric(0), 
                     buc=numeric(0));

    # Dimensions of semidefinite matrix variables 
    prob$bardim <- c(3, 4);
    # Constraint bounds    
    prob$bc <- rbind(blc=c(b, -Inf), 
                     buc=c(b, k));

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barc':
    prob$barc$j <- c(rep(1, length(C1_v)), 
                     rep(2, length(C2_v)));         # Which PSD variable (j)
    prob$barc$k <- c(C1_k, C2_k);                   # Entries: (k,l)->v
    prob$barc$l <- c(C1_l, C2_l);
    prob$barc$v <- c(C1_v, C2_v);

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barA':
    prob$barA$i <- c(rep(1, length(A1_v)+length(A2_v)), 2);   # Which constraint (i)
    prob$barA$j <- c(rep(1, length(A1_v)),
                     rep(2, length(A2_v)),
                     2);                                      # Which PSD variable (j)
    prob$barA$k <- c(A1_k, A2_k, 2);                          # Entries: (k,l)->v
    prob$barA$l <- c(A1_l, A2_l, 1);
    prob$barA$v <- c(A1_v, A2_v, 0.5);

    # Solve the problem
    r <- mosek(prob);

    # Print matrix variable and return the solution
    stopifnot(identical(r$response$code, 0));
    print( list(X1=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1])) );
    print( list(X2=getbarvarMatrix(r$sol$itr$barx[[2]], prob$bardim[2])) );

}

sdo2();

sdo_lmi.R

Listing 15.26 sdo_lmi.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      sdo_lmi.R
#
#  Purpose :   To solve a problem with an LMI and an affine conic constrained problem with a PSD term
#
#                 minimize    Tr [1, 0; 0, 1]*X + x(1) + x(2) + 1
#
#                 subject to  Tr [0, 1; 1, 0]*X - x(1) - x(2) >= 0
#                             x(1) [0, 1; 1, 3] + x(2) [3, 1; 1, 0] - [1, 0; 0, 1] >> 0
#                             X >> 0
##
library("Rmosek")

getbarvarMatrix <- function(barvar, bardim)
{
    N <- as.integer(bardim)
    new("dspMatrix", x=barvar, uplo="L", Dim=c(N,N))
}

sdo_lmi <- function()
{
    # Specify the non-matrix variable part of the problem.
    prob       <- list(sense="min")
    prob$c     <- c(1, 1)
    prob$c0    <- 1

    # Specify variable bounds
    prob$bx    <- rbind(blx=rep(-Inf,2), bux=rep( Inf,2))
    # The following two entries must always be defined, even if set to zero.
    prob$A  <- Matrix(c(0, 0), nrow=1, sparse=TRUE)
    prob$bc <- rbind(blc=rep(-Inf,1),buc=rep(Inf, 1))

    prob$F     <- sparseMatrix(i=c(1, 1, 2, 3, 3, 4),
                               j=c(1, 2, 2, 1, 2, 1),
                               x=c(-1, -1, 3, sqrt(2), sqrt(2), 3), dims=c(4, 2))
    prob$g     <- c(0, -1, 0, -1)
    prob$cones <- matrix(list(), nrow=3, ncol=2)
    prob$cones[,1] <- list("RPLUS", 1, NULL)
    prob$cones[,2] <- list("SVEC_PSD_CONE", 3, NULL)

    # Specify semidefinite matrix variables (one 2x2 block)
    prob$bardim <- c(2)

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barc':
    prob$barc$j <- c(1, 1, 1)
    prob$barc$k <- c(1, 2, 2)
    prob$barc$l <- c(1, 1, 2)
    prob$barc$v <- c(1, 0, 1)

    # Block triplet format specifying the lower triangular part 
    # of the symmetric coefficient matrix 'barF' for the ACC:
    prob$barF$i <- c(1, 1)
    prob$barF$j <- c(1, 1)
    prob$barF$k <- c(1, 2)
    prob$barF$l <- c(1, 1)
    prob$barF$v <- c(0, 1)

    # Solve the problem
    r <- mosek(prob)

    # Print matrix variable and return the solution
    stopifnot(identical(r$response$code, 0))
    print( list(barx=getbarvarMatrix(r$sol$itr$barx[[1]], prob$bardim[1]), xx=r$sol$itr$xx) )
    r$sol
}

sdo_lmi()

simple.R

Listing 15.27 simple.R Click here to download.
##
#  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
#  File :      simple.R
#
#  Purpose :   To demonstrate a very simple example using the Rmosek package
#              to read a problem file, solve the optimization problem and
#              write the solution.
##
library("Rmosek")

simple <- function(filename)
{
    # Read problem
    r <- mosek_read(filename, list(verbose=1, usesol=FALSE, useparam=TRUE))
    stopifnot(identical(r$response$code, 0))
    prob <- r$prob

    # Solve problem
    r <- mosek(prob, list(verbose=1))
    stopifnot(identical(r$response$code, 0))
    sol <- r$sol

    # Print solution
    print(sol)
}

# Run only if called directly from the command-line
if( sys.nframe() == 0L )
{
    args <- commandArgs(trailingOnly=TRUE)
    if( !file.exists(args[1]) ) {
        stop("Expected problem file as input argument.")
    }

    simple(args[1])
}

solutionquality.R

Listing 15.28 solutionquality.R Click here to download.
###
##  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
##
##  File :      solutionquality.R
##
##  Purpose :   To demonstrate how to examine the quality of a solution. 
###
library("Rmosek")

solutionquality <- function()
{
    # Specify problem.
    prob <- list(sense="max")
    prob$c <- c(3, 1, 5, 1)
    prob$A <- sparseMatrix(i=c(1, 1, 1, 2, 2, 2, 2, 3, 3),
                           j=c(1, 2, 3, 1, 2, 3, 4, 2, 4),
                           x=c(3, 1, 2, 2, 1, 3, 1, 2, 3))
    prob$bc <- rbind(blc=c(30,  15, -Inf), 
                     buc=c(30, Inf,   25))
    prob$bx <- rbind(blx=rep(0,4), 
                     bux=c(Inf, 10, Inf, Inf))

    # Solve the problem.
    r <- mosek(prob, list(verbose=1, soldetail=2))
    stopifnot(identical(r$response$code, 0))

    # Examine interior-point solution quality.
    isol <- r$sol$itr
    if( is.element(isol$solsta, list("OPTIMAL")) )
    {
      pobj <- isol$pobjval
      dobj <- isol$dobjval

      abs_obj_gap <- abs(dobj-pobj)
      rel_obj_gap <- abs_obj_gap/(1.0 + min(abs(pobj),abs(dobj)))

      max_primal_viol <- max( isol$maxinfeas$pcon, 
                              isol$maxinfeas$pbound,
                              isol$maxinfeas$pbarvar,
                              isol$maxinfeas$pcone )

      max_dual_viol <- max( isol$maxinfeas$dcon,
                            isol$maxinfeas$dbound,
                            isol$maxinfeas$dbarvar,
                            isol$maxinfeas$dcone )

      ## Assume the application needs the solution to be within
      ## 1e-6 of optimality in an absolute sense. Another approach
      ## would be looking at the relative objective gap 

      cat("Customized solution information.\n")
      cat("  Absolute objective gap: ",abs_obj_gap,"\n")
      cat("  Relative objective gap: ",rel_obj_gap,"\n")
      cat("  Max primal violation  : ",max_primal_viol,"\n")
      cat("  Max dual violation    : ",max_dual_viol,"\n")

      accepted <- TRUE

      if( rel_obj_gap > 1e-6 ) {
        print("Warning: The relative objective gap is LARGE.")
        accepted <- FALSE
      }

      ## We will accept a primal and dual infeasibility of 1e-6. 
      ## These numbers should chosen problem dependent.

      if( max_primal_viol > 1e-6 ) {
        print("Warning: Primal violation is too LARGE")
        accepted <- FALSE
      }

      if( max_dual_viol > 1e-6 ) {
        print("Warning: Dual violation is too LARGE.")
        accepted <- FALSE
      }

      if( accepted ) {
        print("Optimal primal interior-point solution")
        print(isol$xx)
      }

    } else if( is.element(isol$solsta, list("DUAL_INFEASIBLE_CER","PRIMAL_INFEASIBLE_CER")) ) {
        cat("Primal or dual infeasibility certificate found.\n")

    } else {
        cat("Other solution status.\n")
    }
}

solutionquality()