Vertex and barycentric coordinate enumeration for a polytope

Problem 1

Enumerate all the vertices of the polytope defined by: ${\bf Ax}\leq {\bf b}$ where ${\bf A} \in R^{m \times n}$, ${\bf x} \in R^n$, ${\bf b} \in R^m$ and each element of ${\bf x} = \{x_1, x_2, \cdots, x_n \}$ satisfies $x_i \geq 0$. For my problem, $n = 4$, and $m = n + 3$. The first $n$ constraints are always a diagonal matrix (upper bounds on $x$). The structure of the last 3 constraints varies, but there are always two non-zero entries in each row (of the last three rows of ${\bf A}$). For example:

$\left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ a & c & 0 & 0 \\ 0 & d & 0 & e \\ 0 & f & g & 0 \end{array}\right] \left[ \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array}\right] \leq \left[ \begin{array}{c} b_1 \\ b_2 \\ b_3 \\ b_4 \\ b_5 \\ b_6 \\ b_7 \end{array}\right]$

Note: I believe one approach is to use the Simplex algorithm with a constant cost function. I would appreciate any insight whether that is a correct approach and/or how to go about that. Ultimately, I would prefer a (faster) algorithm that exploits the structure of this problem.

Problem 2

Define the set of $n + 1$ evenly-spaced points from 0 to 1

$Q = \{0, ~~1/n, ~~2/n, ~~\cdots, ~~(n-1)/n, ~~1 \}$

Define a sequence of $p$ points where $w_i \in Q$ that sum to unity

$w = \{w_1,w_{2},\cdots, w_{p}\}$ and $\sum_{i=1}^p w_i = 1$.

Suppose repeated points are admissible (i.e. if $p = 4$ the sequences $\{(n-2)/n,~~1/n,~~1/n,~~0\}$ or $\{0,~~ 0,~~ 1,~~ 0\}$ are admissible). The problem is to enumerate all possible sequences $w$.

Note: One approach might be to first compute all unique combinations that sum to unity, and then consider permuations of these. I'm not sure how to do this other than brute-force and I am interested in alternatives.

Motivation for Problem 1 + 2

It might be worth mentioning my motivation, perhaps someone will have a similar problem or point to a method not involving Problem 1,2:

The goal is write an algorithm that uniformly samples a large number of points from the interior of the polytope in Problem 1. I've tried a rejection sampling method but that is too slow. I've avoided hit and run sampling because of its random nature (I need repeatable results). My current approach is :

1) Enumerate all $p$ vertices (Problem 1)

2) Form the barycentric coordinate system: Let $v_j \in R^n$ denote the coordinates of the $j$-th vertex where $j \in \{1, 2, \cdots, p\}$. Any point ${\bf p}$ inside the polytope can be expressed using barycentric coordinates as: ${\bf p} = \sum_{j=1}^p w_j v_j$ where the weights $\sum_{j=1}^p w_j= 1$.

3) Use the set of weights from Problem 2 to obtain a set of points contained in the polytope. Of course they won't be "uniformly" distributed but I think they will be sufficient representative of the polytope for my particular problem.

For posterity...

Problem 1)

I ended up using a vertex enumeration library written by Komei Fukuda: https://www.inf.ethz.ch/personal/fukudak/cdd_home/

Problem 2)

This problem is actually equivalent to finding all the sequences of length $m$ consisting of positive integers (and zero) that sum to $n$. I found the number of such sequences to be $\left( \begin{array}{c} n + m - 1 \\ m - 1\end{array}\right)$ which turns out to be (not surprisingly) a "polytopic" number. For example for $n = 5$ and $m=3$ there are 21 sequences and 21 is a triangular number. The sequences are:

            0        0   5.0000
0   1.0000   4.0000
0   2.0000   3.0000
0   3.0000   2.0000
0   4.0000   1.0000
0   5.0000        0
1.0000        0   4.0000
1.0000   1.0000   3.0000
1.0000   2.0000   2.0000
1.0000   3.0000   1.0000
1.0000   4.0000        0
2.0000        0   3.0000
2.0000   1.0000   2.0000
2.0000   2.0000   1.0000
2.0000   3.0000        0
3.0000        0   2.0000
3.0000   1.0000   1.0000
3.0000   2.0000        0
4.0000        0   1.0000
4.0000   1.0000        0
5.0000        0        0


I wrote a recursive C++ algorithm for generating this set. I'm sure there is a more efficient way to write this, but here it is ...

  // third-party linear algebra library

// returns factorial n!
int factorial(int n){
int nfact = 1;
for (int i = 0; i < n; i++){
nfact = nfact*(i+1);
}
return nfact;
}

// returns bionomial "n choose k"
int binomial(int n, int k){
// binomial(n, k) = n!/[k!(n-k)!]
return (factorial(n) / (factorial(k)*factorial(n-k)));
}

// returns matrix with rows of all permutations of integer pairs that sum to n
// i.e. {(n,0), (n-1,1) , ... , (0,n)}
arma::mat constantSumPairs(int n){
arma::mat out(n+1,2);
for (int i = 0; i < n+1; i++){
out(i,0) = i;
out(i,1) = n-i;
}
return out;
}

// main algorithm: recursively populates the matrix with the desired sequences
void constantSumSubsetsRecursive(int n, int &m, int &curDepth, int &curRow,
arma::mat &seqMat){
if (curDepth == m - 2){ // stopping condition
// last two cols of seqMat will now be populated
arma::mat lastCols = constantSumPairs(n);
// iterate through the list of constant sum pairs
for (int i = 0; i < lastCols.n_rows; i++){
// at each iteration insert a row into seqMat
if (i != 0){ // once curRow is incremented need to copy previous values
for (int j = 0; j < curDepth; j++){ // populate cols 0 to curDepth
seqMat(curRow, j) = seqMat(curRow - 1, j);
}
}
// insert new values
seqMat(curRow, curDepth) = lastCols(i,0); // second-to-last col
seqMat(curRow, curDepth + 1) = lastCols(i,1); // last col
// increment row counter
curRow++;
}
return;
}
else {
// generate a list of constant sum pairs from n
arma::mat listPairs = constantSumPairs(n);
// iterate through this list
for (int i = 0; i < listPairs.n_rows; i++){
if (i != 0){ // for all but the first pair
// copy values between cols 0 and curDepth from previous row to curRow
for (int j = 0; j < curDepth; j++){
seqMat(curRow, j) = seqMat(curRow - 1, j);
}
}
// insert first value of listPair at the current depth
seqMat(curRow, curDepth) = listPairs(i,0);
curDepth++; // prepare to insert the second value
// if the second value is not the last col, then we need to continue
// expanding the curDepth col recurisvely
constantSumSubsetsRecursive(listPairs(i,1), m, curDepth, curRow, seqMat);
}
}
}

// returns a matrix with rows that contain all possible sequences of m integers
// (including zero) that sum to n
arma::mat constantSumSubsets(int n, int m){
int numPerms = binomial(n+m-1, m-1); // get number of sequences
arma::mat output = arma::ones<arma::mat>(numPerms, m); // initialize output
int curDepth = 0; // col of output matrix being populated
int curRow = 0; // row of output matrix being populated
constantSumSubsetsRecursive(n, m, curDepth, curRow, output); // start alg.
return output;
}

// the program
int main() {
// n defines the n+1 length set: N = {0, 1, 2, ... , n}
int n = 5;
// m is the length of the sequence that should sum to m
// i.e. we want all the sequences {s1, s2, ... , m} with elements si in N
// that have the property : s1 + s2 + ... + sm = 1
int m = 3;
// runs the algorithm to generate an armadillo matrix with the result
arma::mat out = constantSumSubsets(n, m);
// prints the result
out.print("out");
return 0;
}


(Note: the code uses the armadillo linear algebra library for the arma::mat type but this could easily be replaced with an array)