CVXvsMOSEK.py

import time

import numpy as np
import cvxpy as cp

from mosek.fusion import Model, Domain, Expr, ObjectiveSense


def cvx_explicit(cov_factor, residual_var, F, mu, gamma, Lmax):
    n, m = F.shape
    w = cp.Variable(n)
    f = cp.Variable(m)
    ret = mu.T @ w
    risk = cp.sum_squares(
        np.transpose(np.linalg.cholesky(cov_factor)) @ f
    ) + cp.sum_squares(np.sqrt(residual_var) @ w)
    problem = cp.Problem(
        cp.Maximize(ret - gamma * risk),
        [cp.sum(w) == 0, f == F.T @ w, cp.norm(w, 1) <= Lmax],
    )
    problem.solve()
    return w.value


def cvx_implicit(cov_factor, residual_var, F, mu, gamma, Lmax):
    n, m = F.shape
    w = cp.Variable(n)
    ret = mu.T @ w
    risk = cp.sum_squares(
        np.transpose(np.linalg.cholesky(cov_factor)) @ (F.T @ w)
    ) + cp.sum_squares(np.sqrt(residual_var) @ w)
    problem = cp.Problem(
        cp.Maximize(ret - gamma * risk),
        [cp.sum(w) == 0, cp.norm(w, 1) <= Lmax],
    )
    problem.solve()
    return w.value


# a case for monkey-patching?
def sum_squares(model, vector):
    x = model.variable()
    model.constraint(Expr.vstack(x, 0.5, vector), Domain.inRotatedQCone())
    return x


# a case for monkey-patching
def one_norm(model, vector):
    p = model.variable(n, Domain.unbounded())
    model.constraint(Expr.hstack(p, vector), Domain.inQCone())
    return Expr.sum(p)


def mosek_implicit(cov_factor, residual_var, F, mu, gamma, Lmax):
    n, m = F.shape

    with Model("cash-neutral-factor") as M:
        w = M.variable("w", n, Domain.unbounded())

        # position in factor space
        f = Expr.mul(F.T, w)

        expected_return = Expr.dot(mu, w)

        # cash-neutral
        M.constraint("cash-neutral", Expr.sum(w), Domain.equalsTo(0.0))

        # Leverage
        M.constraint("leverage", one_norm(M, w), Domain.lessThan(Lmax))

        # Residual variances
        res_var = sum_squares(M, Expr.dot(np.sqrt(residual_var), w))

        # Variance in factor space
        factor_var = sum_squares(
            M, Expr.mul(np.transpose(np.linalg.cholesky(cov_factor)), f)
        )

        risk = Expr.add(factor_var, res_var)
        M.objective(
            "obj",
            ObjectiveSense.Maximize,
            Expr.sub(expected_return, Expr.mul(gamma, risk)),
        )
        M.solve()
        return w.level()


if __name__ == "__main__":
    n = 2000
    m = 200

    mu = np.random.randn(n)

    # construct a factor covariance matrix
    Sigma_tilde = np.random.randn(m, m)
    Sigma_tilde = Sigma_tilde.T.dot(Sigma_tilde)

    # construct idiosyncratic variances
    D = np.random.uniform(0, 0.9, size=n)

    # construct the factor loadings
    F = np.random.randn(n, m)

    Lmax = 2
    gamma = 0.1

    # solve the problem with the explicit construction of the f variable
    t = time.time()
    for i in range(5):
        w = cvx_explicit(
            cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
        )
        print(mu.T @ w)
    print(time.time() - t)

    # solve the problem but without the explicit construction of the f variable
    t = time.time()
    for i in range(5):
        w = cvx_implicit(
            cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
        )
        print(mu.T @ w)
    print(time.time() - t)

    # solve the same problem with Mosek fusion
    t = time.time()
    for i in range(5):
        w = mosek_implicit(
            cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
        )
        print(mu.T @ w)
    print(time.time() - t)

Careful please. This is work in progress. We are solving a standard Markowitz problem with realistic dimensions, e.g. 2000 assets and 200 factors. The core idea is to project the weights (in asset space) down into the lower dimensional factor space.