Vertex and barycentric coordinate enumeration for a polytope

Question

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.

Answer 1

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
  #include<armadillo>

  // 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);
        curDepth--; // return to previous col
      } 
    }
  }

  // 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)