# Uniform lattice sample inside a particular convex polytope

[update]:

hardmath suggests using tools from linear programming. This looks like a good idea indeed. I can now tell that my feasible set is described by:

$Set = \{d \in \mathbb{N}^c, -B.d\leqslant x_0\}, x_0 \in \mathbb{N}^a, B \in {\cal M}_{a,c}(\mathbb{N})$

.. which is a convex polytope of $\mathbb{R}^c$, or rather its intersection with $\mathbb{N}^c$. From the underlying process of my matrix, I know it is finite and not empty. The question is, then: how do I uniformly sample from $Set$?

[update]:

Even better than a toy example, here is the very structure of the matrix B defining my particular polytope: $B = \left(\begin{matrix}C(l,r)\\Identity(l.r)\end{matrix}\right)$, with $C(l,r)$ a patterned matrix entirely defined by only $(l,r) \in \mathbb{N}^2$. You will understand the pattern with this small example:

$C(3,4)=\left(\begin{matrix} -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & -1 & -1 & -1 \\ -1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -1 \\ \end{matrix}\right)$

As far as $x_0$ is concerned, it only contains $l+r$ values taken in $\mathbb{N^*}$ then $l.r$ zeroes, for a small example $x_0=(6,11,7,5,2,6,8,0,0,0,0,0,0,0,0,0,0,0,0)$.

This is the reason why I can tell that zero is a vertex of the polytope, that each axis of $\mathbb{R}^{c=l.r}$ contains an edge of the polytope, and that each direction is bounded between zero and the minimum of only two planes, whose value I can calculate from the current value of $d$ in $B.d \leqslant x_0$. (This is very light from a computational point of view since I don't even need to use $B$ to get the bounds: the pattern helps me a lot and I have got an efficient $upperbound=f(d{\tiny(current state)},i{\tiny(direction)},l,r)$)

Here are my final questions then, with total loss of generality ;)

• Let $\cal{D}$ be the uniform distribution of $d$ in $Set$: $d \hookrightarrow \cal D$. Is it true that each full conditional distribution of $\cal D$ is uniform on $\{0,...,upperbound\}$ ? $d_{i|_{d_{j\neq i}}} \hookrightarrow {\cal U}(\{0,\dots,upperbound(d_{j\neq i},i,l,r)\}$

• If it is, is it true that starting from $d = 0$ then Gibbs-sampling successively in each successive direction will produce a Markov Chain whose stationnary distribution is $\cal D$?

• If it is not, how do I uniformly sample from $Set$?

[original question]:

I need to generate random values of $x \in \mathbb{N}^a$ (column), and they must respect the constraint

$f=A^Tx\quad$

, with $f \in \mathbb{N}^b$ and $A \in \cal{M}_{a,b}$

I don't want to generate naive $x$'s then reject the incorrect values. All the possible values of $x$ should be uniformly sampled.

One interesting first step is to consider the kernel of $A^T$, $ker(A^T) = \left\{v,A^Tv = 0 \right\}$

I have a basis of this space, which I call $B \in \cal{M}_{a,c}$, and I have a $x_0$ verifying the constraint. I can then generate, for any $d \in \mathbb{N}^c$, $B.d = v \in ker(A^T)$ then $x = x_0 + v$ verifying the constraint since $A^Tx=A^T(x_0 + v) = A^Tx_0 = f$.

But I cannot safely sample $d$ from $\{0,\dots,constant\}^c$ since there is no guarantee that each element of $x = x_0 + v$ will be positive.

How can I find the correct 'uniform' degrees of freedom $\nu \in \mathbb{N}^c$ and the corresponding function $g: \nu \mapsto d$ that would correctly and smoothly describe all possible values of x, then sampling from them?

• I'm stumped by the constraint. What does $^T Ax$ mean? I suspect something is being transposed, but I don't know this notation. – hardmath Mar 16 '15 at 17:00
• @hardmath Wops, sorry. I used ${}^T\!A$ for the transposition of $A$. The way we build $A$ is just not the same as the way we write the constraint. One should read the constraint $f=({}^T\!A).x$ – iago-lito Mar 16 '15 at 17:14
• Would it not be better to write $A^T$? Or at least equvalent? – hardmath Mar 16 '15 at 17:19
• @hardmath Maybe it would. I learned the other way round and I thought this just depended on your religion, but.. well, I'll change it. What matters most is that we should agree on the meaning, right? ;) – iago-lito Mar 16 '15 at 17:58
• Right! You may already have considered this, but I would try using linear programming to identify vertices of the feasible region, perhaps with a way of eliminating some unknowns in terms of "free" variables. If the feasible region of free variables can be bounded, then we could sample them uniformly and test for the dependent variables to be positive. – hardmath Mar 16 '15 at 20:20

The topic of uniformly sampling from (integer) lattice points within a rational convex polytope is often considered in the literature in connection with the problem of counting such points (known to be NP-hard, as the corresponding decision problem is NP-complete).

See the paper "Sampling lattice points," by Kannan and Vempala (1997), which proposes a nearly uniform sampling method when the polytope contains a ball of sufficient radius. In some special cases their sampling is exactly uniform. See also the lengthy survey paper "An Algorithmic Theory of Lattice Points in Polyhedra," by Barnivok and Pommersheim (1999).

You seem to have made the right first steps toward solving, or at least formulating the problem. Using exact integer row and column operations, the matrix $A^T$ can be reduced to "diagonal" Smith normal form:

$$S A^T T = D = \begin{pmatrix} d_1 & 0 & 0 & \ldots & 0 \\ 0 & d_2 & 0 & \ldots & 0 \\ 0 & 0 & d_3 & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 \\ 0 & 0 & 0 & \ldots & 0 \end{pmatrix}$$

where $S,T$ are invertible $b\times b$ and $a\times a$ matrices respectively, and the diagonal entries of the product $S A^T T$ satisfy $d_i \mid d_{i+1}$ for $1 \le i \lt b$, with trailing entries allowed to be zeros. Note that $D$ is $b\times a$ and a matrix equivalent to $A^T$.

The system $A^T x = f$ can then be written in terms of $x = Ty$ where $y$ solves $D y = S f$. In this form the system if consistent (feasible) if and only if each $i$th component of $S f$ is a multiple of $d_i$ respectively. Corresponding components of $y$ are determined when $d_i$ is nonzero and otherwise free. A basis for the free portion of the solution space, previously discussed in terms of a basis for the kernel of $A^T$, can be constructed.

The equality constraints' integer solutions are thus easily parameterized (in exact integer arithmetic), but finding a nonnegative integer solution remains potentially quite difficult (even though nonnegative (feasible) rational solutions can be found using typical linear programming techniques). We are fortunate that here the existence of a feasible integer solution is guaranteed by the underlying problem.

The data $-B d \le x_0$ for the specific problem (see Question revision) can be interpreted in terms of a nonnegative integer $\ell \times r$ matrix $M$:

$$M = \begin{pmatrix} d_{11} & \ldots & d_{1r} \\ \vdots & & \vdots \\ d_{\ell 1} &\ldots & d_{\ell r} \end{pmatrix}$$

whose row sums are bounded by the first $\ell$ entries of $x_0$, and whose column sums are bounded by the next $r$ entries of $x_0$. Taking $M$ as the zero matrix affords one feasible integer solution.

That is, let $u$ be a positive integer column vector of length $\ell$ and $v$ a positive integer column vector of length $r$ such that:

$$x_0^T = \left( \begin{array}{c|c|c} u^T & v^T & 0 \end{array} \right)$$

By identifying the integer column vector $d$ with the unrolled entries of $M$, we have $\ell$ "row sum" constraints:

$$d_{i1} + d_{i2} + \ldots + d_{ir} \le u_i \; , \; 1 \le i \le \ell$$

with $r$ "column sum" constraints:

$$d_{1j} + d_{2j} + \ldots + d_{\ell j} \le v_j \; , \; 1 \le j \le r$$

and $\ell\cdot r$ nonnegativity constraints $d_{ij} \ge 0$ for all $1 \le i \le \ell$ and $1 \le j \le r$, corresponding to the zero entries in the tail of $x_0$.

• Thank you for this update! I'll have a glance at those papers today. As far as my particular problem is concerned, I know that 1) my polytope is positive and that 2) every axis of $\mathbb{R}^c$ contains an edge of my polytope (so 0 is a vertex) and that 3) for any given dimension $i \in {1,\dots,c}$, if $d$ stands inside the polytope, it can move along dimension $i$ from zero to a "ceiling" which is the minimum of only 2 planes, whose expression I know from $(B,d,x_0)$. Is there anything wrong then using Gibbs sampling in this case? (I can sample from every full conditional distribution) – iago-lito Mar 23 '15 at 7:38
• Perhaps you mean that your $x_0$ is a vertex and the polytope edges at $x_0$ are parallel to the coordinate axes, proceeding in the positive direction? Is there a "toy" version of $A^Tx = f$? – hardmath Mar 23 '15 at 14:14
• Yes there is. See the update in a minut :) – iago-lito Mar 23 '15 at 15:25