Constrained System / Combinatorial Problem

Question

Let $x\in \mathbb{R}^{n}$, $Y\in \mathbb{R}^{mxn}$.

We can then define:

$row_{i}(Y)=$ $i^{th}$ row of $Y$
$column_{i}(Y)=$ $i^{th}$ column of $Y$
$x_{i}=i^{th}$ element of $x$
$sum(x)=$ sum of the element of $x$
$card(x)=$ number of non-zero element of $x$.

My Goal is to devise a method to determine if it is possible to construct a matrix $k\in \mathbb{R}_{\geq 0}^{mxn}$ given some $a\in \mathbb{R}^{n}$, $b, c\in\mathbb{R}^{m}$ for some $m, n \in \mathbb{N}_{>0}$ where $sum(a)=sum(b)$

s.t. for all $ i, j \in \mathbb{N}_{\geq0}$ with $0<i<m$ and $0<j<n$:

$\bullet$ $sum(row_{i}(K))=b_{i}$
$\bullet$ $sum(column_{j}(K))=a_{j}$
$\bullet$ $card(row_{i}(K))=c_{i}$

Any suggestions?

Answer 1

So you are essentially asking how you can constrain the cardinality of a vector. The answer below is based on a standard approach called big-M modelling.

Let $x$ be a vector for which you want exactly $c$ non-zeros. Introduce a new binary vector $y$ which indicates if an element is zero or non-zero. Let $M$ be an upper bound on the possible value of an element in $x$ at optimality (from domain knowledge)

If $y_i=0$ then $x_i$ should be $0$:

$x \leq My$

This also ensures $y_i=1$ when $x_i$ is non-zero.

If $y_i=1$ then $x_i$ non-zero. This is a bit tricky, since you have to define what non-zero is in the numerical precision you are working with:

$ x \geq 0.00001y$

It also ensures $y_i=0$ if $x_i=0$

Exactly $c$ elements are non-zero

$\sum y = c$

All remaining constraints are standard linear programming constraints, and the logic above is applied on every row of $K$.

The code below implements an example in the MATLAB Toolbox YALMIP (disclaimer:developed by me)

m = 5;
n = 10;
% create a known solution
K0 = sprand(m,n,.25)
a = sum(K0,2);
b = sum(K0,1);
c = sum(K0>0,2);

% Define decision variables
K = sdpvar(m,n,'full')
Y = binvar(m,n,'full');
M = max(max(a),max(b));

Model = [sum(K,1)==b, sum(K,2)==a, K>=0];
Model = [Model, K <= M*Y, K >= 0.00001*Y,sum(Y,2)==c];

% Let us find a solution with smallest largest element
Objective = max(max(K));
solvesdp(Model,Objective);