The source files for all examples can be found in /examples.
Markowitz Portfolio Optimization, Three Ways
This example shows how to solve a Markowitz portfolio optimization problem using both the GenericSolver and the QPSolver.
Specifically, we want to solve the problem
\[\begin{array}{ll} \text{minimize} & (1/2)\gamma x^T\Sigma x - \mu^Tx \\ \text{subject to} & \mathbf{1}^Tx = 1 \\ & x \geq 0, \end{array}\]
where $\Sigma$ is the covariance matrix of the returns of the assets, $\mu$ is the expected return of each asset, and $\gamma$ is a risk adversion parameter. The variable $x$ represents the fraction of the total wealth invested in each asset.
using GeNIOS
using Random, LinearAlgebra, SparseArraysGenerating the problem data
Note that we generate the covariance matrix $\Sigma$ as a diagonal plus low-rank matrix, which is a common model for financial data and referred to as a 'factor model'
Random.seed!(1)
k = 5
n = 100k
# Σ = F*F' + diag(d)
F = sprandn(n, k, 0.5)
d = rand(n) * sqrt(k)
μ = randn(n)
γ = 1;QPSolver interface
To use the QPSolver, we need to specify $P$, $q$, $M$, $l$, and $u$. The Markowitz portfolio optimization problem is equivalent to the following 'standard form' QP:
\[\begin{array}{ll} \text{minimize} & (1/2)x^T(\gamma \Sigma) x + (-\mu)Tx \\ \text{subject to} & \begin{bmatrix}1 \\ 0\end{bmatrix} \leq \begin{bmatrix} I \\ \mathbf{1}^T \end{bmatrix} x \leq \begin{bmatrix}\infty \\ 1\end{bmatrix} \end{array}\]
P = γ*(F*F' + Diagonal(d))
q = -μ
M = vcat(I, ones(1, n))
l = vcat(zeros(n), ones(1))
u = vcat(Inf*ones(n), ones(1))
solver = GeNIOS.QPSolver(P, q, M, l, u)
res = solve!(solver; options=GeNIOS.SolverOptions(eps_abs=1e-6))
println("Optimal value: $(round(solver.obj_val, digits=4))")Starting setup...
Setup in 0.000s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 0.000e+00 Inf Inf 1.000e+00 0.000
1 0.000e+00 1.000e+00 3.201e+01 1.000e+00 0.223
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 0.266
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.317
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.373
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.434
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.494
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.555
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.620
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.687
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.755
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.826
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.887
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.948
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 1.008
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.069
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.130
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.190
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.248
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.305
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 1.361
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 1.418
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.475
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.532
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.589
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.646
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.698
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.748
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.760
SOLVED in 1.760s, 525 iterations
Total time: 1.761s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4289Performance improvements
We can also define custom operators for $P$ and $M$ to speed up the computation. For this, we will use LinearMaps.jl, which allows us to build these operators without ever explicitly constructing the final matrices.
If we don't implement size for all objects passed to the solver, we need set check_dims to false.
using LinearMaps
# P = γ*(F*F' + Diagonal(d))
F_lm = LinearMap(F)
P = γ*(F_lm*F_lm' + Diagonal(d))
# M = vcat(I, ones(1, n))
M = vcat(LinearMap(I, n), ones(1, n))
solver = GeNIOS.QPSolver(P, q, M, l, u; check_dims=false);
res = solve!(solver; options=GeNIOS.SolverOptions(eps_abs=1e-6));
println("Optimal value: $(round(solver.obj_val, digits=4))")Starting setup...
Setup in 0.410s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 0.000e+00 Inf Inf 1.000e+00 0.000
1 0.000e+00 1.000e+00 3.201e+01 1.000e+00 0.882
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 1.189
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 1.190
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 1.192
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 1.194
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 1.196
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 1.198
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 1.201
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 1.203
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 1.206
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 1.221
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 1.223
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 1.224
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 1.226
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.228
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.229
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.231
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.232
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.234
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 1.236
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 1.237
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.239
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.240
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.242
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.243
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.245
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.246
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.246
SOLVED in 1.246s, 525 iterations
Total time: 1.657s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4289Custom P and M
Of course, we can also define each needed operations ourselves, which may give an additional speedup for some matrices. Here, we avoid allocations for the intermediate product when multiplying by $F^T$ and then $F$.
struct FastP
F
d
γ
vk
end
function LinearAlgebra.mul!(y::AbstractVector, P::FastP, x::AbstractVector)
mul!(P.vk, P.F', x)
mul!(y, P.F, P.vk)
@. y += P.d*x
@. y *= P.γ
return nothing
end
Pfast = FastP(F, d, γ, zeros(k))
struct FastM
n::Int
end
Base.size(M::FastM) = (M.n+1, M.n)
Base.size(M::FastM, d::Int) = d <= 2 ? size(M)[d] : 1
function LinearAlgebra.mul!(y::AbstractVector, M::FastM, x::AbstractVector)
y[1:M.n] .= x
y[end] = sum(x)
return nothing
end
LinearAlgebra.adjoint(M::FastM) = Adjoint{Float64, FastM}(M)
function LinearAlgebra.mul!(x::AbstractVector{T}, M::Adjoint{T, FastM}, y::AbstractVector{T}) where T <: Number
@. x = y[1:M.parent.n] + y[end]
return nothing
end
function LinearAlgebra.mul!(x::AbstractVector{T}, M::Adjoint{T, FastM}, y::AbstractVector{T}, α::T, β::T) where T <: Number
@. x = α * ( y[1:M.parent.n] + y[end] ) + β * x
return nothing
end
P = FastP(F, d, γ, zeros(k))
M = FastM(n)
solver = GeNIOS.QPSolver(P, q, M, l, u; check_dims=false);
res = solve!(solver; options=GeNIOS.SolverOptions(eps_abs=1e-6));
println("Optimal value: $(round(solver.obj_val, digits=4))")Starting setup...
Setup in 0.179s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 0.000e+00 Inf Inf 1.000e+00 0.000
1 0.000e+00 1.000e+00 3.201e+01 1.000e+00 0.308
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 0.309
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.310
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.312
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.313
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.315
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.316
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.318
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.319
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.321
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.323
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.324
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.326
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 0.327
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 0.329
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 0.339
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 0.340
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 0.342
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 0.343
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 0.345
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 0.346
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 0.347
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 0.349
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 0.350
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 0.351
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 0.353
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 0.354
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 0.354
SOLVED in 0.354s, 525 iterations
Total time: 0.533s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4289GenericSolver interface
The GenericSolver interface is more flexible and allows for some speedups via an alternative problem splitting. We will solve the equivalent problem
\[\begin{array}{ll} \text{minimize} & (1/2)\gamma x^T\Sigma x - \mu^Tx + I(z) \\ \text{subject to} & -x + z = 0 \end{array}\]
where $I(z)$ is the indicator function for the set $\{z \mid z \ge 0 \text{ and } \mathbf{1}^Tz = 1\}$. The gradient and Hessian of $f(x) = (1/2)\gamma x^T\Sigma x$ are easy to compute. The proximal operator of $g(z) = I(z)$ is simply the projection on this set, which can be solved via a one-dimensional root-finding problem (see the appendix of our paper).
params = (; F=F, d=d, μ=μ, γ=γ, tmp=zeros(k))
# f(x) = γ/2 xᵀ(FFᵀ + D)x - μᵀx
function f(x, p)
F, d, μ, γ, tmp = p.F, p.d, p.μ, p.γ, p.tmp
mul!(tmp, F', x)
qf = sum(abs2, tmp)
qf += sum(i->d[i]*x[i]^2, 1:length(x))
return γ/2 * qf - dot(μ, x)
end
# ∇f(x) = γ(FFᵀ + D)x - μ
function grad_f!(g, x, p)
F, d, μ, γ, tmp = p.F, p.d, p.μ, p.γ, p.tmp
mul!(tmp, F', x)
mul!(g, F, tmp)
@. g += d*x
@. g *= γ
@. g -= μ
return nothing
end
# ∇²f(x) = γ(FFᵀ + D)
struct HessianMarkowitz{T, S <: AbstractMatrix{T}} <: HessianOperator
F::S
d::Vector{T}
vk::Vector{T}
end
function LinearAlgebra.mul!(y, H::HessianMarkowitz, x)
mul!(H.vk, H.F', x)
mul!(y, H.F, H.vk)
@. y += d*x
return nothing
end
Hf = HessianMarkowitz(F, d, zeros(k))
# g(z) = I(z)
function g(z, p)
T = eltype(z)
return all(z .>= zero(T)) && abs(sum(z) - one(T)) < 1e-6 ? zero(T) : Inf
end
function prox_g!(v, z, ρ, p)
z_max = maximum(w->abs(w), z)
l = -z_max - 1
u = z_max
# bisection search to find zero of F
while u - l > 1e-8
m = (l + u) / 2
if sum(w->max(w - m, zero(eltype(z))), z) - 1 > 0
l = m
else
u = m
end
end
ν = (l + u) / 2
@. v = max(z - ν, zero(eltype(z)))
return nothing
end
solver = GeNIOS.GenericSolver(
f, grad_f!, Hf, # f(x)
g, prox_g!, # g(z)
I, zeros(n); # M, c: Mx - z = c
params=params
)
res = solve!(solver)
println("Optimal value: $(round(solver.obj_val, digits=4))")Starting setup...
Setup in 0.156s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 Inf Inf Inf 1.000e+00 0.000
1 0.000e+00 4.472e-02 2.292e+01 1.000e+00 0.230
20 -2.496e+00 2.936e-02 2.738e-03 1.000e+00 0.231
40 -2.457e+00 1.970e-02 1.993e-03 1.000e+00 0.232
60 -2.447e+00 1.151e-02 4.535e-03 2.000e+00 0.234
80 -2.438e+00 6.178e-03 1.995e-03 2.000e+00 0.235
100 -2.433e+00 4.025e-03 9.519e-04 2.000e+00 0.236
120 -2.431e+00 2.971e-03 5.204e-04 2.000e+00 0.237
140 -2.429e+00 2.248e-03 3.124e-04 2.000e+00 0.238
SOLVED in 0.238s, 140 iterations
Total time: 0.394s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4295This page was generated using Literate.jl.