Discrete constrained optimization problem

Question

What would be the best way to solve -- either analytically or algorithmically (in this case preferably using Python) -- a discrete constrained optimization problem of the form

$$ \vec x^\star = \arg\min_{\vec x}(\vec x^T C \ \vec x) \qquad\text{subject to } \sum_{i=1}^n x_i = p $$ where $x \in \{0,1\}^n$, $C \in \mathbb R^{n \times n}$, for some $p \in \mathbb N$.

Update

In response to Brian's answer, $C$ is a correlation matrix restricted to $-1 \leq C_{ij} \leq 1$ and in particular positive semi-definite while $n \approx 500$.

Answer 1

This is a 0-1 quadratic integer programming problem. Like all 0-1 integer programming problems with polynomial objective functions, it can be transformed into a 0-1 linear integer programming problem and solved by branch-and-bound or cutting plane methods. There are a number of open-source and proprietary commercial software packages that can solve the transformed problem for you (CBC, SCIP, CPLEX, GuRoBi, etc.)

The linearization technique is as follows. For each pair of variables $x_{i}$ and $x_{j}$ (We only need to consider the case $i<j$ since #$x_{i}x_{j}=x_{j}x_{i}$), add a 0-1 variable $y_{i,j}$, where we will ensure that $y_{i,j}=x_{i}x_{j}$. Then replace every term $C_{i,j}x_{i}x_{j}$ in the objective function with $C_{i,j}y_{i,j}$. Also, since $x_{i}^{2}=x_{i}$, replace every objective function term $C_{i,i}x_{i}^{2}$ with $C_{i,i}x_{i}$. Now your objective function will be linear in the $x_{i}$ and $y_{i,j}$ variables.

To ensure that $y_{i,j}=x_{i}x_{j}$, add the linear constraints:

$x_{i}+x_{j} - y_{i,j} \leq 1$

$ -x_{i}+y_{i,j} \leq 0$

$ -x_{j} +y_{i,j} \leq 0$.

You can easily check that for every 0-1 combination of $x_{i}$ and $x_{j}$ values, these constraints ensure that $y_{i,j}=x_{i}x_{j}$.

The major problem with this technique is that it adds a large number of 0-1 variables ($n(n-1)/2$ additional variables) and constraints ($3n(n-1)/2$ additional constraints) to the problem. When $n$ is large this becomes impractical.

If the matrix $C$ is positive semidefinite then it may be more efficient to solve the problem using branch and bound on the quadratic formulation. This is supported, for example, by GuRoBi.

How big is your $n$? Is $C$ positive semidefinite?

In practice, a lot will depend on the size of the problem. If $n$ is small (say 10 to 100), then there's a good chance you can solve this problem in a reasonable amount of time using the linearization approach. If $n$ is larger but $C$ is PSD, then it is possible to solve the problem by applying branch and bound directly to the quadratic formulation. If $n$ is large (thousands or millions), but $C$ is not PSD, then it may be practically impossible to solve your problem to optimality.

The OP has updated the question to explain that $C$ is a PSD covariance matrix and $n \approx 500$. For problems like this, which commonly arise in mathematical finance, branch and bound on the QP formulation is generally the best approach. I'd look at CPLEX, GuRoBi, Mosek, or XPRESS for such problems. I'm not aware of any good open-source solver for them.

Answer 2

If you do a linearization with $y_{ij}=x_i x_j$ for $i < j$, you might consider adding an RLT cut based on your cardinality constraint: $$\sum_{i=1}^{j-1} y_{ij} + \sum_{i=j+1}^n y_{ji} = (p - 1) x_j.$$

Answer 3

Here's the code I ended up writing to solve this problem to optimality using Gurobi's Python API:

import pandas as pd
from gurobipy import GRB, Model, quicksum

# Required to help Gurobi finds its license file. An academic license can be obtained
# for free at https://www.gurobi.com/downloads/end-user-license-agreement-academic.
# Only works in Jupyter notebooks. Else use os.environ["GRB_LICENSE_FILE"] = "path/to/gurobi.lic".
%env GRB_LICENSE_FILE path/to/gurobi.lic

# Create a model for solving the quadratic optimization problem of selecting n_select
# out of n_vars candidates with the least pairwise correlation according to the matrix corr_mat.
model = Model("quadratic_problem")
# model.params.LogFile = "results/gurobi.log"

n_vars, n_select = len(candidates), 20

# Create decision variables.
dvars = model.addVars(n_vars, vtype=GRB.BINARY).values()

# Define the model objective to minimize the sum of pairwise correlations.
obj = corr_mat.dot(dvars).dot(dvars)
model.setObjective(obj, GRB.MINIMIZE)

# Add L1 constraint on dvars so that the optimization returns at least n_select formulas.
constr = model.addConstr(quicksum(dvars) >= n_select, "l1_norm")

model.optimize()

# Save selected materials to dataframe and CSV file.
gurobi_candidates = candidates.iloc[[bool(var.x) for var in dvars]]
gurobi_candidates.to_csv("gurobi_candidates.csv")