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.006
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 0.101
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.159
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.222
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.289
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.357
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.424
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.496
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.569
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.655
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.735
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.802
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.869
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 0.936
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.003
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.070
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.135
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.198
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.261
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 1.323
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 1.385
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.448
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.511
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.573
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.644
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.704
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.759
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.772
SOLVED in 1.773s, 525 iterations
Total time: 1.773s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.250s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.497
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 1.084
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 1.085
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 1.088
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 1.090
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 1.093
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 1.096
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 1.099
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 1.102
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 1.106
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 1.109
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 1.112
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 1.115
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 1.117
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.120
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.123
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.126
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.136
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.138
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 1.140
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 1.143
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.145
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.147
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.149
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.151
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.153
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.154
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.155
SOLVED in 1.155s, 525 iterations
Total time: 1.405s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.154s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.095
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 0.143
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.145
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.148
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.150
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.153
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.156
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.159
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.163
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.166
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.170
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.173
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.177
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 0.180
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 0.183
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 0.186
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 0.189
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 0.193
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 0.196
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 0.206
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 0.209
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 0.211
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 0.214
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 0.216
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 0.219
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 0.221
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 0.223
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 0.224
SOLVED in 0.224s, 525 iterations
Total time: 0.377s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.144s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.148
20 -2.496e+00 2.936e-02 2.738e-03 1.000e+00 0.149
40 -2.457e+00 1.970e-02 1.993e-03 1.000e+00 0.151
60 -2.447e+00 1.151e-02 4.535e-03 2.000e+00 0.152
80 -2.438e+00 6.178e-03 1.995e-03 2.000e+00 0.154
100 -2.433e+00 4.025e-03 9.519e-04 2.000e+00 0.155
120 -2.431e+00 2.971e-03 5.204e-04 2.000e+00 0.157
140 -2.429e+00 2.248e-03 3.124e-04 2.000e+00 0.158
SOLVED in 0.158s, 140 iterations
Total time: 0.302s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4295This page was generated using Literate.jl.