Solving Markowitz style portfolio optimization with Rmosek

I am trying to understand how mosek works and formulated a formulated a simple min xTSigmax

I use the example data from the Rmosek codebook and tried to formulate the problem in conic form.

The quadratic term is substituted by t following the guides on the mosek website. t is then minimized.

# 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 ) )`

library("Rmosek")

# Sample covariance matrix
Sigma <- t(GT)%*%GT

# Number of assets
n <- ncol(GT)

# Problem definition
prob <- list(sense="min")

# Objective: We're minimizing t, which represents portfolio variance
prob$c <- c(rep(0, n), 1)

# Conic constraint: t >= x^T Sigma x
prob$F <- rbind(
  cbind(Matrix(0.0,nrow=1,ncol=n), 1), 
  cbind(GT, Matrix(0, n, 1))
)
prob$g <- c(rep(0, n),1)

prob$cones <- matrix(list("QUAD", 1+n, NULL), nrow=3, ncol=1)
rownames(prob$cones) <- c("type","dim","conepar")

# Linear constraint: Sum of x[i] = 1
prob$A <- Matrix(c(rep(1, n), 0), 1)

prob$bc <- rbind(blc=1,
                 buc=1)

# Bounds: 0 <= x[i] <= 1, 0 <= t
prob$bx <- rbind(blx=c(rep(0, n + 1)), 
                 bux=c(rep(1, n), Inf))

# Solve the problem
r <- mosek(prob,list(verbose=10))

xTestMosekSimpleMinVarConic <- r$sol$itr$xx[1:n]

stddevsol <- r$sol$itr$xx[n+1]

stddevcalc <- (t(xTestMosekSimpleMinVarConic)%*%Sigma%*%xTestMosekSimpleMinVarConic)^0.5
expret <- drop(mu %*% xTestMosekSimpleMinVarConic)

print(sprintf("xTestMosekSimpleMinVarConic"))
print(sprintf("Expected return: %.4e Std. deviation solver: %.4e Std. deviation calc: %.4e", round(expret,5), round(stddevsol,5), stddevcalc))
print(round(xTestMosekSimpleMinVarConic,5))

When I solve this the result is:

[1] "Expected return: 1.7350e-01 Std. deviation solver: 1.0214e+00 Std. deviation calc: 2.0803e-01"
[1] 0.13469 0.14807 0.29008 0.23863 0.00000 0.11060 0.07793 0.00000

However when I solve the same problem as a quadratic problem i get:

[1] "Expected return: 1.6625e-01 Std. deviation: 2.0373e-01 "
[1] 0.11332 0.11437 0.30221 0.18154 0.00000 0.05632 0.04519 0.18706

the code:

library("Rmosek")

Sigma <- t(GT)%*%GT
# Number of assets
n <- length(mu)

# Define the QP problem
prob <- list(sense = "min")

# Convert the lower triangular part of the covariance matrix Sigma 
# into the required list format for qobj
i_indices <- c()
j_indices <- c()
values <- c()

for (i in 1:n) {
  for (j in 1:i) {
    if (Sigma[i, j] != 0) {
      i_indices <- append(i_indices, i)
      j_indices <- append(j_indices, j)
      values <- append(values, Sigma[i, j])
    }
  }
}

prob$qobj <- list(i = i_indices, j = j_indices, v = values)

prob$c <-c(rep(0,n))

# Equality constraint for sum of weights: 1^T * x = 1
prob$A <- matrix(rep(1, n), 1)
prob$bc <- rbind(blc = 1, buc = 1)

# Non-negativity constraints on asset weights: x >= 0
prob$bx <- rbind(blx = rep(0, n), bux = rep(Inf, n))

# Solve the problem
solution <- mosek(prob,list(verbose=1))

xTestMosekSimpleQuad <- solution$sol$itr$xx
stddev <- (t(xTestMosekSimpleQuad)%*%Sigma%*%xTestMosekSimpleQuad)^0.5
expret <- drop(mu %*% xTestMosekSimpleQuad)

print(sprintf("xTestMosekSimpleQuad"))
print(sprintf("Expected return: %.4e Std. deviation: %.4e ", round(expret,5), round(stddev,5)))
print(round(xTestMosekSimpleQuad,5))

The same result I get when I solve the problem with CVXR using the standard solver but also with the mosek-solver

CXVR code:

library("CVXR")

Sigma <- t(GT)%*%GT

## Form problem
w <- Variable(n)
ret <- t(mu) %*% w
risk <- quad_form(w, Sigma)
constraints <- list(w >= 0, sum(w) == 1)
objective <- risk
prob <- Problem(Minimize(objective), constraints)
result <- solve(prob, solver='MOSEK', verbose = TRUE)
risk_data <- result$getValue(sqrt(risk))
ret_data <- result$getValue(ret)
w_data <- result$getValue(w)

CVXR is handling it as a conic problem:

Problem Name :
Objective sense : minimize
Type : CONIC (conic optimization problem) Constraints : 19
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 1
Scalar variables : 19
Matrix variables : 0
Integer variables : 0

Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - primal attempts : 1 successes : 1
Lin. dep. - dual attempts : 0 successes : 0
Lin. dep. - primal deps. : 0 dual deps. : 0
Presolve terminated. Time: 0.00
Optimizer - threads : 12
Optimizer - solved problem : the primal
Optimizer - Constraints : 10
Optimizer - Cones : 1
Optimizer - Scalar variables : 18 conic : 10
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00
Factor - dense det. time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 54 after factor : 54
Factor - dense dim. : 0 flops : 1.20e+03
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.0e+00 6.0e-01 1.4e+00 0.00e+00 0.000000000e+00 -3.973982284e-01 1.0e+00 0.00
1 4.1e-01 2.5e-01 1.1e-01 7.52e-01 4.724270646e-02 -3.621647720e-01 4.1e-01 0.00
2 1.5e-01 8.9e-02 2.3e-02 2.20e+00 3.208365717e-02 -4.069524762e-02 1.5e-01 0.00
3 2.0e-02 1.2e-02 1.1e-03 1.39e+00 4.349220651e-02 3.528925879e-02 2.0e-02 0.00
4 4.1e-03 2.4e-03 9.7e-05 1.05e+00 4.314371279e-02 4.153013503e-02 4.1e-03 0.00
5 6.2e-04 3.7e-04 7.0e-06 1.01e+00 4.148713867e-02 4.126481240e-02 6.2e-04 0.00
6 7.6e-05 4.6e-05 3.0e-07 1.00e+00 4.150471852e-02 4.147717911e-02 7.6e-05 0.00
7 1.0e-05 6.2e-06 1.5e-08 1.00e+00 4.150531100e-02 4.150163293e-02 1.0e-05 0.00
8 2.5e-07 1.5e-07 5.7e-11 1.00e+00 4.150536024e-02 4.150527034e-02 2.5e-07 0.00
9 5.1e-10 3.0e-10 5.0e-15 1.00e+00 4.150536306e-02 4.150536288e-02 5.0e-10 0.00
Optimizer terminated. Time: 0.00

Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: 4.1505363060e-02 nrm: 1e+00 Viol. con: 2e-10 var: 0e+00 cones: 0e+00
Dual. obj: 4.1505362883e-02 nrm: 2e-01 Viol. con: 0e+00 var: 1e-10 cones: 0e+00
Optimizer summary Optimizer - time: 0.00
Interior-point - iterations : 9 time: 0.00
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00

Where is my mistake in the Mosek code?

Solution

First thing I noticed is that your comment says

# Conic constraint: t >= x^T Sigma x

where x^T Sigma x = x^T GT^T GT x = || GT x ||_2^2, but then you use the quadratic cone and essentially get

t >= || GT x ||_2

which corresponds to t >= sqrt(x^T Sigma x). This difference is fine if t only appears in the objective (in fact, it is beneficial, since it gives you better scaling), but you need to be aware of the difference when you interpret the result.

Second thing I noticed is that you have

prob$g <- c(rep(0, n),1)

but this adds a constant 1 to last entry of the conic constraint essentially giving you t >= || GT x + e_n ||_2 where e_n is a unit vector. Try using

prob$g <- c(rep(0, n+1))

and see if that helps.