Sierpinski (<1K) Primal Simplex Algorithm
    "Fast cars, fast women, fast algorithms... what more could a man want?" -- Joe Mattis.

The Basic Algorithm
This is a technique used to solve the linear programming problem:
Maximise a linear objective function q.x subject to the linear constraints   Ax = b over x > 0 where q and x are n dimensional vectors, b is an m dimensional vector, and A is an mn matrix.

If m=F, n=F+3 and q has the form (q1,q2,q3,0,0,..,0) then the geometrical interpretation of the problem is that we are trying to find the point inside an F-faced convex hull closest to an exterior plane having normal (q1,q2,q3).

The primal simplex algorithm operates by hoping from vertex to adjacent vertex, always moving closer to the desired plane; calculating the vertex positions from the given set of planes ("facets.").   

We will illustrate its use using the same numeric example discussed in "Algorithms" (Robert Sedgewick, Addison-Wesley 1984) but approached in a slightly different way. The problem as presented by Sedgewick is:
Maximise x1 + x2x3 subject to the constraints
-x1+x2 £ 5
x1+4x2 £ 45
2x1+x2 £ 27
3x1-4x2 £ 24
x3 £ 4
over x1,x2,x3 ³ 0.
We reexpress this in standard form by introducing  slack variables x4,x5,x6,x7,x8 thus
Maximise M=x1 + x2x3 subject to the constraints
-x1+ x2+x4=5
x1+4x2+x5=45
2x1+ x2+x6=27
3x1-4x2+x7=24
x3+ x8=4
over x1,x2,x3,x4,x5,x6,x7,x8 ³ 0.
This is the  standard form of the problem with q=(1,1,1,0,0,0,0,0)T, b=(5,45,27,24,4)T and
A = æ -11010000 ö
   ç 1 4001000 ÷
   ç 2 1000100 ÷
   ç 3-4000010 ÷
   è 0 0100001 ø

Viewed geometrically, the inequalities for x1,x2 and x3 (including x1 ³ 0,x2 ³ 0,x3 ³ 0) describe a 3D convex hull having the origin (0,0,0) as one vertex. The slack variable x4 gives the distance of the point (x1,x2,x3) from the plane -x1 + x2 = 5, x5 gives the distance from x1 + 4x2 = 45, and so on. x1,x2 and x3 may themselves be thought of as describing the distances of the point (x1,x2,x3) from the planes x1=0, x2=0, and x3=0 respectively

The first step of the Primal Simplex Algorithm is to find a feasible solution to the problem, that is, find an x>0 satisfying Ax=b. If any solutions exist, ie. if the problem is solvable, then there will be an optimal solution having at most m non-zero terms in x (and thus at least nm zero terms), so we look for a feasible solution of this form.

Geometrically, this is saying that there will always be an optimal point for the objective function that is a vertex of the hull.

An obvious candidate solution (used in Sedgewick) is the "origin"   x=(0,0,0,5,45,27,24,4) but for the sake of generality we will start with x=(8,0,0,13,37,11,0,4).

We call the m non-zero variables in x (x1,x4,x5,x6 and x8)  basis variables and need to rearrange our equations into  canonical form with respect to the basis variables. This means that each basis variable can appear in only one equation and within that equation must have coefficient 1.
When the equations are arranged in canonical form, it is clear what value x takes since the basis variables each appear in only one equation and the non-basis variables are known to be zero.
(If we had taken x4,x5,x6,x7, and x8 as our basis variables as Sedgewick does then the equations would already be canonical form)
We can rearrange the fourth constraint as x1 - 1.33x2+ 0.33x7 = 8 and then replace x1 in the other equations with 8+1.33x2-0.33x7 thus:

Maximise 2.33x2 + x3 -0.33x7 = M-8
subject to the constraints
 -0.33x2+x4+0.33x7=13
  5.33x2+x5-0.33x7=37
  3.67x2+x6-0.67x7=11
      x1-1.33x2+0.33x7=8
        x3+x8=4
over x1,x2,x3,x4,x5,x6,x7,x8 ³ 0.
Geometrically, the non-basis variables correspond to the planes that meet at the vertex we are currently at.

We can express these equations more succinctly using an (m+1)(n+1) matrix thus:
B= æ 0.00-2.33-1.00 0.00 0.00 0.00 0.33 0.00  8.00 ö
  ç 0.00-0.33 0.00 1.00 0.00 0.00 0.33 0.00 13.00 ÷
  ç 0.00 5.33 0.00 0.00 1.00 0.00-0.33 0.00 37.00 ÷   ç 0.00 3.67 0.00 0.00 0.00 1.00-0.67 0.00 11.00 ÷
  ç 1.00-1.33 0.00 0.00 0.00 0.00 0.33 0.00  8.00 ÷
  è 0.00 0.00 1.00 0.00 0.00 0.00 0.00 1.00  4.00 ø

The top (0th) row of B is the negation of the objective function; the rightmost entry of this row (B0n) contains the current value of the objective function. The remainder of the rightmost column contains the values of the basis variables.
Note that we will be counting rows from zero, but columns from one.

Suppose we wish to change our basis variables from x1,x4,x5,x6,x8 to x1,x2,x4,x5,x8 , ie. to replace x6 with x2.
To rearrange the equations into canonical form with respect to the new basis we rewrite the third constraint as
x2 + 0.27x6 - 0.18x7 = 3
and then plug it into the other equations to obtain
æ 0.000.00-1.000.000.00 0.64-0.090.0015.00 ö
ç 0.000.00 0.001.000.00 0.09 0.270.0014.00 ÷
ç 0.000.00 0.000.001.00-1.45 0.640.0021.00 ÷
ç 0.001.00 0.000.000.00 0.27-0.180.00 3.00 ÷
ç 1.000.00 0.000.000.00 0.36 0.090.0012.00 ÷
è 0.000.00 1.000.000.00 0.00 0.001.00 4.00 ø

which exemplifies the feasible solution (12,3,0,14,21,0,0,4). This is equivalent to multiplying the third row by 1/3.67 and then adding appropriate multiples of it to the other rows in order to zero their second entries - an operation very similar to that used in Gaussian Elimination.

If we next wanted to replace x8 with x3 in our basis we would use the fifth row (the row which defines x8) to zero the third entry of each of the other rows to move to the solution (12,3,4,14,21,0,0,0).
This operation is known as "pivoting about" B53.

Pivoting about Bij >0 corresponds to bringing xj into the basis instead of xk, where k is the unique value in {1,2,..,n} for which Bik ¹ 0. We force xk to be zero and allow xj to be nonzero. Geometrically this corresponds to moving (x1,x2,x3) away from the (j+4)th face and into the (k+4)th face.

This pivoting procedure gives an easy way of obtaining successive   potential solutions to the problem. If, after pivoting, the rightmost column is non-negative then the new solution is feasible but as yet we have no way of guaranteeing that it is "better" than the last solution. This is done by careful choice of pivot point with reference to the top row.
If xi is any basis variable then B0i = 0. Thus if B0j < 0 then making xj positive nonzero at the expense of xi will increase the objective function.

The rightmost column will transform as Bk(n+1)' = Bk(n+1)  BkjBi(n+1)/Bij so we can ensure it remains nonnegative by choosing i to minimise Bi(n+1)/Bij over those i such that Bij >0.
Bij £ 0 for all i Î {1,2,..,n} can only occur if the constaint A x=b is unbounded, ie. permits some xi to be infinite.

At each stage, therefore, we first choose a column j satisfying B0j <0. If there are more than one such then we can choose the one with the smallest entry in row zero (greatest increment); the one that would give the greatest increase in the objective function (steepest descent); or select at random. Having selected j we then choose i to minimise Bi(n+1)/Bij over those i such that Bij >0. If there are two or more such, then we select at random or choose the one that removes the variable with lowest index from the basis (ie. minimises k where k is uniquely defined by Bik¹0). It is important to choose randomly rather than, say, the first or last ones found in order to avoid cycling, an infinite sequence of pivots that do not increase the objective function.

Once there are no negative entries in the top row then we have an optimal solution.

The cost of each pivot is (n-m+1)D+(m(n-m+1))M., however all the D's are by the same value and one is into 1 so we can replace them with 1R+(nm)M giving 1R+((n-m)(m+1)+m)M.

The primal simplex algorithm as presented above needs a feasible solution having at most m nonzero values in order to begin. If such a solution cannot be found by inspection then we can use the primal simplex algorithm to find one by introducing n new variables xn+1,xn+2,..,xn+m, (adding one to each of the m constraint equations) and replacing the objective function with åi=n+1n+m xi.

If we can achieve an objective value of 0 for this new problem then (possibly after a few further pivots to remove any remaining "new" variables from the basis) we will have a feasible solution to the original constraints.

Application to 3D Convex Hulls
    The standard form for the primal simplex algorithm is not quite the form of the problem that arises when trying to find the vertex of a convex hull H =(H,d,V) having centre c and  orientation A nearest the plane r.q=g since we do not necessarily have (0,0,0) as a vertex.
We take variables x1=r1,x2=r2,x3=r3, and slack variables x4,x5,..,xF+3.
Our equations are
Maximise q1x1+q2x2+q3x3
Subject to  (Ani)1x1+(Ani)2x2+(Ani)3x3x3+i  £ eii=1,2,..F.
where ei = di+(Ani).c ; i=1,2,..F.
We do not impose x1,x2,x3 ³ 0 and do not necessarily have 0 Î H. We can assume that d ³ 0 however.
The matrix for these equations is the (F+1)(F+4) matrix
B= æ -q1-q2-q3000..00 ö
  ç   100..0e1 ÷
  ç AH 010..0e2 ÷
  ç   ........ ÷
  è   000..1eF ø
Our initial solution is r=Avh where vh Î V is one of the vertices (preferably one likely to be only a few edges away from the optimal vertex) so that
x1=(Avh1, , )
x2=(Avh2, , )
x3=(Avh3, , )
x3+i=ei-ni.vh i=12..F
Since vh is a vertex, it must lie in at least three of the planes and so at least three of the slack variables x3+i are zero (non-basis). Let these be x3+s,x3+t, and x3+u.
The B presented above is in canonical form with respect to the basis x4,x5,...,x3+F. We need to rearrange it so that it is canonical form with respect to the basis formed by replacing x3+s,x3+t, and x3+u with x1,x2, and x3. The sth,tth, and uth equations taken together enable us to do this by expressing x1,x2, and x3 in terms of x3+s,x3+t, and x3+u. The three equations are
(Ans)1x1 + (Ans)2x2 + (Ans)3x3 + x3+s = es
(Ant)1x1 + (Ant)2x2 + (Ant)3x3 + x3+t = et
(Anu)1x1 + (Anu)2x2 + (Anu)3x3 + x3+u = eu
which we can rewrite as
(AN)T æ x1 ö + æ x3+t ö = æ es ö
         ç x2 ÷ ç x3+t ÷ ç et ÷
         è x3 ø è x3+u ø è eu ø
where N=(ns,nt,nu) is bound to be invertible since no two of the normals can be parallel (or the planes would not have met at vh).
Let M=((AN)T)-1=(NTAT)-1=AN-T. We have
æ x1 ö +M æ x3+t ö =M æ es ö
ç x2 ÷ ç x3+t ÷ ç et ÷
è x3 ø è x3+u ø è eu ø
and so we can replace the sth,tth,and suh rows of B with the rows
123...s t u F+4
(100...m110..0m120..0m130..0 es'=m11es+m12et+m13eu)
(010...m210..0m220..0m230..0 et'=m21es+m22et+m23eu)
and(001...m310..0m320..0m330..0 eu'=m31es+m32et+m33eu)
(Cost 9M)
Having done this we subtract multiples of these three rows from all the other rows. For i¹0, s,t,u we replace the ith row
123 ...s ...t ...i ...u ...F+4
(bi1bi2bi30...00...00...10...00...ei)
with
(0000...bis'0...bit'0...10...biu'0... ei')
where
bis' = -bi1m11-bi2m21-bi3m31
bit' = -bi1m12-bi2m22-bi3m32
biu' = -bi1m13-bi2m23-bi3m33
ei'  = ei-bi1es'-bi2et'-bi3eu'
(Cost 12M for each of F-3 rows).
We replace the top row
123 ...s ...t ...i ...u ...F+4
(-q1-q2-q30...00...00...00...00...0)
with
(0000...qs'0...qt'0...10...qu'0... Q)
qs' = q1m11+q2m21+q3m31
qt' = q1m12+q2m22+q3m33
qu' = q1m13+q2m23+q3m33
Q  = q1es'+q2et'+q3eu'
(Cost 12M).
Having performed this "simultaneous triple-pivot", we have B in canonical form with respect to a vertex of the hull and we can proceed with the primal simplex algorithm.

We cannot use the algorithm exactly as given above because the assumption x ³ 0 does not hold for the first three components of x. We have "lost" the three planes corresponding to the first three variables (x1=0, x2=0, and x3=0) and so must be more careful in our choice of pivot. We get round this by keeping x1,x2, and x3 in the basis at all times. Recalling that pivoting about Bij >0 corresponds to bringing xj into the basis at the expense of xk, where k is the unique value in {1,2,..,n} for which Bik ¹; 0, we see that we will never want to pivot about Bij if j=1,2, or 3 and must avoid pivoting about Bij if any of Bi1,Bi2, or Bi3 is nonzero.
This means that we cannot, as before, for a given choice of j choose i to minimise Bi(n+1)/Bij over those i such that Bij >0 so we cannot guarantee that the rightmost column will remain nonnegative.
However, if we choose i to minimise Bi(n+1)/Bij over those i such that Bij >0 and Bi1=Bi2=Bi3=0 then the only entries in the rightmost column that can become negative are the first second and third which is okay since this corresponds merely to vertices having negative coordinates.
On each iteration we will have a choice of three pivot columns, the cost of choosing the pivot row is (F-3)D or less, and the cost of the pivot operation is 4D+4FM simplifying to 1R+(4F+3)M.


Generating the Vertices of a Convex Hull

The Primal simplex algorithm's ability to moves from vertex to vertex of a convex hull is potentially useful as a means to generate the vertices V from H and d. Might AV be so derived from AH and d faster than by explicitly rotationg each vertex? A "Vertex Generation Algorithm" has no need for an objective function so we drop the zeroth row of B. At each stage we will have a free choice   of three pivot columns, however (apart from on the first iteration) one of these would take us back to the vertex we came from.

Topologically, we need to visit each node of a trinary graph, a   depth first recursive traversal will serve but requires B to be placed on a stack every bifurcation. Each vertex is the intersection of three planes (the i,j, and kth say).  On travelling from one vertex <i,j,k> to another <i,j,l> we flag <i,j> as a "traversed edge-pair". We begin our first iteration with no traversed pairs. At each subsequent iteration we will be at a vertex <i,j,k> considering the i,j, and kth columns as possible pivot columns.
If <i,j> is a traversed pair we know not to choose the kth column because <i,j,k> must have been already visited. Likewise we do not choose the jth column if <i,k> has been traversed, and reject the ith column if <j,k> has been traversed. If we reject all three pivot columns we exit from the recursive procedure; if only one column is acceptable we choose it; if two (or all three) are acceptable we must choose both (or all three) and impliment them one after another by preserving B and then calling the traversal procedure recursively.

Though this method will generate V it has two drawbacks. Firstly, it is recursive and requires the large matrix B to be preserved over the recursion. Secondly, it will generate any vertex at which more than three planes intersect (the peak of a square based pyramid for example) more than once.
It is not easy to solve either of these problems if all that is known about H is H and d (possibly the case for the hull-compiler), but if the connectivity of the hull is known (likely for the hull-rotater) then we can eliminate both these problems if we provide the algorithm with  an "edge-walk" {i1,i2,...iP}, a sequence of pivot columns that specify a non- recursive traversal of the hull that traverses all the edges of the hull at the expense of possibly visting some vertices more than once.
We can save the algorithm a further (F-3)D per iteration by providing the pivot row as well, effectively providing a "vertex-walk" {(i1,j1),(i2,j2),..(iP,jP)}.
Each iteration will take 1R+(4F-1)M since we no longer have the top row of B to maintain. This compares badly with the 9M (or 6M+3T) required to rotate a v Î V by A.

We conclude therefore that it is faster to rotate hull vertices than to recreate them from pre-rotated planes.


Glossary   Contents   Author
Copyright (c) Ian C G Bell 1998
Web Source: www.iancgbell.clara.net/maths or www.bigfoot.com/~iancgbell/maths
18 Nov 2006.