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.