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.106
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.166
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.242
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.310
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.378
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.446
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.518
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.593
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.668
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.750
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.818
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.886
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 0.954
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.022
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.090
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.157
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.233
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.297
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.425
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.489
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.552
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.616
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.680
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.738
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.794
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.808
SOLVED in 1.808s, 525 iterations
Total time: 1.808s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.244s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.524
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 1.120
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 1.122
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 1.124
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 1.127
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 1.129
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 1.132
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 1.135
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 1.138
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 1.141
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 1.144
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 1.147
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 1.150
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 1.152
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 1.155
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 1.158
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 1.161
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 1.163
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 1.166
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 1.169
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 1.171
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 1.174
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 1.177
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 1.179
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 1.182
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 1.184
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 1.187
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 1.187
SOLVED in 1.187s, 525 iterations
Total time: 1.431s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.167s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.086
20 -2.347e+01 5.416e-01 2.331e-01 1.000e+00 0.136
40 -1.095e+01 1.979e-01 4.927e-02 1.000e+00 0.138
60 -7.166e+00 1.048e-01 1.925e-02 1.000e+00 0.140
80 -5.457e+00 6.878e-02 9.222e-03 1.000e+00 0.143
100 -4.660e+00 5.029e-02 5.210e-03 1.000e+00 0.145
120 -4.140e+00 3.793e-02 3.837e-03 1.000e+00 0.148
140 -3.742e+00 2.874e-02 2.895e-03 1.000e+00 0.151
160 -3.435e+00 2.181e-02 2.191e-03 1.000e+00 0.153
180 -3.198e+00 1.657e-02 1.660e-03 1.000e+00 0.156
200 -3.017e+00 1.260e-02 1.259e-03 1.000e+00 0.159
220 -2.802e+00 8.214e-03 2.483e-03 2.000e+00 0.162
240 -2.686e+00 6.028e-03 1.527e-03 2.000e+00 0.164
260 -2.615e+00 4.589e-03 1.027e-03 2.000e+00 0.167
280 -2.571e+00 3.580e-03 7.681e-04 2.000e+00 0.169
300 -2.537e+00 2.818e-03 5.827e-04 2.000e+00 0.172
320 -2.512e+00 2.233e-03 4.470e-04 2.000e+00 0.174
340 -2.488e+00 1.662e-03 3.209e-04 2.000e+00 0.177
360 -2.472e+00 1.245e-03 2.334e-04 2.000e+00 0.179
380 -2.460e+00 9.374e-04 1.717e-04 2.000e+00 0.181
400 -2.450e+00 1.097e-03 6.081e-04 2.000e+00 0.184
420 -2.443e+00 7.788e-04 4.377e-04 2.000e+00 0.186
440 -2.438e+00 5.548e-04 3.146e-04 2.000e+00 0.189
460 -2.434e+00 3.963e-04 2.261e-04 2.000e+00 0.191
480 -2.432e+00 2.837e-04 1.625e-04 2.000e+00 0.193
500 -2.430e+00 2.036e-04 1.168e-04 2.000e+00 0.195
520 -2.429e+00 1.464e-04 8.414e-05 2.000e+00 0.198
525 -2.429e+00 1.349e-04 7.753e-05 2.000e+00 0.198
SOLVED in 0.198s, 525 iterations
Total time: 0.366s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.163s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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.142
20 -2.496e+00 2.936e-02 2.738e-03 1.000e+00 0.144
40 -2.457e+00 1.970e-02 1.993e-03 1.000e+00 0.145
60 -2.447e+00 1.151e-02 4.535e-03 2.000e+00 0.146
80 -2.438e+00 6.178e-03 1.995e-03 2.000e+00 0.148
100 -2.433e+00 4.025e-03 9.519e-04 2.000e+00 0.149
120 -2.431e+00 2.971e-03 5.204e-04 2.000e+00 0.151
140 -2.429e+00 2.248e-03 3.124e-04 2.000e+00 0.152
SOLVED in 0.152s, 140 iterations
Total time: 0.316s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Optimal value: -2.4295This page was generated using Literate.jl.