If m=F, n=F+3 and q has the form (q_{1},q_{2},q_{3},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 (q_{1},q_{2},q_{3}).
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 x_{1} + x_{2} + x_{3} subject to the constraints
-x_{1} | + | x_{2} | £ 5 |
x_{1} | + | 4x_{2} | £ 45 |
2x_{1} | + | x_{2} | £ 27 |
3x_{1} | - | 4x_{2} | £ 24 |
x_{3} | £ 4 |
-x_{1} | + | x_{2} | + | x_{4} | = | 5 |
x_{1} | + | 4x_{2} | + | x_{5} | = | 45 |
2x_{1} | + | x_{2} | + | x_{6} | = | 27 |
3x_{1} | - | 4x_{2} | + | x_{7} | = | 24 |
x_{3} | + | x_{8} | = | 4 |
A = | æ | -1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | ö |
ç | 1 | 4 | 0 | 0 | 1 | 0 | 0 | 0 | ÷ | |
ç | 2 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | ÷ | |
ç | 3 | -4 | 0 | 0 | 0 | 0 | 1 | 0 | ÷ | |
è | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | ø |
Viewed geometrically, the inequalities for x_{1},x_{2} and x_{3} (including x_{1} ³ 0,x_{2} ³ 0,x_{3} ³ 0) describe a 3D convex hull having the origin (0,0,0) as one vertex. The slack variable x_{4} gives the distance of the point (x_{1},x_{2},x_{3}) from the plane -x_{1} + x_{2} = 5, x_{5} gives the distance from x_{1} + 4x_{2} = 45, and so on. x_{1},x_{2} and x_{3} may themselves be thought of as describing the distances of the point (x_{1},x_{2},x_{3}) from the planes x_{1}=0, x_{2}=0, and x_{3}=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
(x_{1},x_{4},x_{5},x_{6} and x_{8})
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 x_{4},x_{5},x_{6},x_{7}, and x_{8} as our basis variables as Sedgewick
does then the equations would already be canonical form)
We can rearrange the fourth constraint as x_{1} - 1.33x_{2}+ 0.33x_{7} = 8
and then replace x_{1} in the other equations with 8+1.33x_{2}-0.33x_{7}
thus:
Maximise 2.33x_{2} + x_{3 }-0.33x_{7} = M-8_{
}
subject to the constraints
-0.33x_{2} | + | x_{4} | + | 0.33x_{7} | =13 |
5.33x_{2} | + | x_{5} | - | 0.33x_{7} | =37 |
3.67x_{2} | + | x_{6} | - | 0.67x_{7} | =11 |
x_{1} | - | 1.33x_{2} | + | 0.33x_{7} | =8 |
x_{3} | + | x_{8} | =4 |
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 (0^{th}) row of B is the negation of the objective function; the
rightmost entry of this row (B_{0n}) 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
x_{1},x_{4},x_{5},x_{6},x_{8}
to
x_{1},x_{2},x_{4},x_{5},x_{8}
, ie. to replace x_{6} with x_{2}.
To rearrange the equations into canonical form with respect to the new basis
we rewrite the third constraint as
x_{2} + 0.27x_{6} - 0.18x_{7} = 3
and then plug it into the other equations to obtain
æ | 0.00 | 0.00 | -1.00 | 0.00 | 0.00 | 0.64 | -0.09 | 0.00 | 15.00 | ö |
ç | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.09 | 0.27 | 0.00 | 14.00 | ÷ |
ç | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | -1.45 | 0.64 | 0.00 | 21.00 | ÷ |
ç | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.27 | -0.18 | 0.00 | 3.00 | ÷ |
ç | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.36 | 0.09 | 0.00 | 12.00 | ÷ |
è | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.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 x_{8} with x_{3}
in our basis we would use the fifth row (the row which defines x_{8})
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" B_{53}.
Pivoting about B_{ij} >0 corresponds to bringing x_{j}
into the basis instead of x_{k}, where k is the unique value in {1,2,..,n} for
which B_{ik} ¹ 0.
We force x_{k} to be zero and allow x_{j} to be nonzero.
Geometrically this corresponds to moving (x_{1},x_{2},x_{3}) 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 x_{i} is any basis variable then B_{0i} = 0.
Thus if B_{0j} < 0 then making x_{j}
positive nonzero at the expense of x_{i} will increase the objective
function.
The rightmost column will transform as
B_{k(n+1)}' = B_{k(n+1)} B_{kj}B_{i(n+1)}/B_{ij} so
we can ensure it remains nonnegative by choosing i to minimise B_{i(n+1)}/B_{ij}
over those i such that B_{ij} >0.
B_{ij} £ 0 for all i Î {1,2,..,n} can only occur if the
constaint A x=b is unbounded, ie. permits some x_{i} to be infinite.
At each stage, therefore, we first choose a column j satisfying B_{0j} <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 B_{i(n+1)}/B_{ij}
over those i such that B_{ij} >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
B_{ik}¹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 x_{n+}1,x_{n+}2,..,x_{n+m}, (adding one to each of the m constraint equations) and replacing the objective function with å_{i=n+1}^{n+m} x_{i}.
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 x_{1}=r_{1},x_{2}=r_{2},x_{3}=r_{3}, and slack variables x_{4},x_{5},..,x_{F+}3.
Our equations are
Maximise q_{1}x_{1}+q_{2}x_{2}+q_{3}x_{3}
Subject to
(An_{i})_{1}x_{1}+(An_{i})_{2}x_{2}+(An_{i})_{3}x_{3}+ x_{3+i} £ e_{i} ; i=1,2,..F.
where e_{i} _{}= d_{i}+(An_{i}).c ; i=1,2,..F.
We do not impose x_{1},x_{2},x_{3} ³ 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= | æ | -q_{1} | -q_{2} | -q_{3} | 0 | 0 | 0..0 | 0 | ö |
ç | 1 | 0 | 0..0 | e_{1} | ÷ | ||||
ç | AH | 0 | 1 | 0..0 | e_{2} | ÷ | |||
ç | . | . | .... | .. | ÷ | ||||
è | 0 | 0 | 0..1 | e_{F} | ø |
x_{1} | = | (Av_{h}_{1}, , ) | |||||
x_{2} | = | (Av_{h}_{2}, , ) | |||||
x_{3} | = | (Av_{h}_{3}, , ) | |||||
x_{3+i} | = | e_{i}-n_{i}.v_{h} | i=1 | 2 | .. | F |
(AN)^{T} | æ | x_{1} | ö | + | æ | x_{3} | + | t | ö | = | æ | e_{s} | ö |
ç | x_{2} | ÷ | ç | x_{3} | + | t | ÷ | ç | e_{t} | ÷ | |||
è | x_{3} | ø | è | x_{3} | + | u | ø | è | e_{u} | ø |
æ | x_{1} | ö | + | M | æ | x_{3} | + | t | ö | = | M | æ | e_{s} | ö |
ç | x_{2} | ÷ | ç | x_{3} | + | t | ÷ | ç | e_{t} | ÷ | ||||
è | x_{3} | ø | è | x_{3} | + | u | ø | è | e_{u} | ø |
1 | 2 | 3 | ... | s | t | u | F+4 | ||||||
( | 1 | 0 | 0 | ... | m_{11} | 0..0 | m_{12} | 0..0 | m_{13} | 0..0 | e_{s}'=m_{11}e_{s}+m_{12}e_{t}+m_{13}e_{u} | ) | |
( | 0 | 1 | 0 | ... | m_{21} | 0..0 | m_{22} | 0..0 | m_{23} | 0..0 | e_{t}'=m_{21}e_{s}+m_{22}e_{t}+m_{23}e_{u} | ) | |
and | ( | 0 | 0 | 1 | ... | m_{31} | 0..0 | m_{32} | 0..0 | m_{33} | 0..0 | e_{u}'=m_{31}e_{s}+m_{32}e_{t}+m_{33}e_{u} | ) |
1 | 2 | 3 | ... | s | ... | t | ... | i | ... | u | ... | F+4 | |||||||
( | b_{i1} | b_{i2} | b_{i3} | 0 | ... | 0 | 0 | ... | 0 | 0 | ... | 1 | 0 | ... | 0 | 0 | ... | e_{i} | ) |
with | |||||||||||||||||||
( | 0 | 0 | 0 | 0 | ... | b_{is}' | 0 | ... | b_{it}' | 0 | ... | 1 | 0 | ... | b_{iu}' | 0 | ... | e_{i}' | ) |
1 | 2 | 3 | ... | s | ... | t | ... | i | ... | u | ... | F+4 | |||||||
( | -q_{1} | -q_{2} | -q_{3} | 0 | ... | 0 | 0 | ... | 0 | 0 | ... | 0 | 0 | ... | 0 | 0 | ... | 0 | ) |
with | |||||||||||||||||||
( | 0 | 0 | 0 | 0 | ... | q_{s}' | 0 | ... | q_{t}' | 0 | ... | 1 | 0 | ... | q_{u}' | 0 | ... | Q | ) |
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 (x_{1}=0, x_{2}=0, and x_{3}=0)
and so must be more careful in our choice of pivot.
We get round this by keeping
x_{1},x_{2}, and x_{3} in the basis at all times.
Recalling that pivoting about B_{ij}
>0 corresponds to bringing x_{j} into the basis at the expense of x_{k}, where k is the unique value in {1,2,..,n} for which B_{ik} ¹; 0, we see that we will
never want to pivot about B_{ij} if j=1,2, or 3 and must avoid pivoting about B_{ij}
if any of B_{i1},B_{i2}, or B_{i3} is nonzero.
This means that we cannot, as before, for a given choice of j choose i
to minimise B_{i(n+1)}/B_{ij} over those i such that B_{ij} >0 so we cannot guarantee
that the rightmost column will remain nonnegative.
However, if we choose i to minimise B_{i(n+1)}/B_{ij} over those i such that
B_{ij} >0 and B_{i1}=B_{i2}=B_{i3}=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" {i_{}1,i_{2},...i_{P}}, 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"
{(i_{1},j_{1}),(i_{2},j_{2}),..(i_{P},j_{P})}.
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.