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, SparseArrays
Generating 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.4289
Performance 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.4289
Custom 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.4289
GenericSolver 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.4295
This page was generated using Literate.jl.