1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.