User:Raul Miller/Simplex

Today, I had occasion to be interested in linear programming, and -- since I usually find that I get a lot more out of J implementations than random descriptions, I performed the obvious google search. That search should find this page, before too long, but all I got out of the search was that I should be searching for "simplex" rather than for "linear programming".

After some more searching, I found a variety of relevant and irrelevant pages:

• Useful: Stories/RichardBrown
• Interesting: User:Devon McCormick/convexHull

Different concept, same word: http://www.jsoftware.com/help/jforc/loopless_code_iii--adverbs__.htm and http://www.jsoftware.com/help/phrases/geometry.htm

Also, some other files from J's math/misc package use the word simplex, and some of them seem relevant:

system/packages/math/simplex.ijs	Lots of comment, with a short implementation.	uses pivot from '~system/packages/math/matutil.ijs'
system/packages/math/simplexnr.ijs	Lots of comments, with a bulky implementation	has lots of test cases
system/packages/math/gamesolver.ijs	Example using simplexnr implementation

I initially decided that Richard Brown's treatment was most useful to me, but that I wanted to simplify his approach. Note, also that Richard Brown's example is the same example used in simplex.ijs.

However, since other pages might change, and since his initial statement differed slightly from his re-expressions of that example, here's a re-statement of that example:

maximize $z=2x+3y$

subject to $4x+y\leq 12$ and $2x+5y\leq 15$

with $x\geq 0$ and $y\geq 0$ .

Using brute force, I expect to find a case as least as good as the case where $x=2,y=2$ and $z=10$

z=:2&*@[ + 3&*@[
c1=: 12 >: 4&*@[ + ]
c2=: 15 >: 2&*@[ + 5&*@]

   (#"1~[:(=>./)z/)(#"1~c2/)(#"1~c1/),:~0.1*i:200
2.1
2.1
   2.1 z 2.1
10.5

But, of course, that's not a very good way to go about solving this kind of problem, and I really should not assume that x and y are equal.

One approach to representing this kind of problem uses a system of equations: $\%\%betteralignmentneeded{\begin{aligned}z-2x-3y&=0\\4x+y+s&=12\\2x+5y+t&=15\end{aligned}}$

  z - 2x - 3y         = 0
      4x +  y + s     = 12
      2x + 5y     + t = 15

where the first equation is a re-arranged version of our original "maximize this quantity" equation, and the remaining equations introduce unknowns which must remain positive (s and t, in this example). In J, we can represent the particulars of this kind of system of equations using a table of coefficients:

M=: 1 _2 _3 0 0 0,0 4 1 1 0 12,:0 2 5 0 1 15

1	_2	_3	0	0	0
0	4	1	1	0	12
0	2	5	0	1	15

Following Richard Brown's narrative, we can pivot this matrix using an expression of the form:

M2=: M - (C-B)*/(R%P)

  NB. where
  R=: 0 2 5 0 1 15 NB. pivot row
  P=: 5            NB. pivot value
  C=: _3 1 5       NB. pivot column
  B=: 0 0 1        NB. pivot row selector

This, of course, just begs to be expressed in J.

NB. define pivotat such that  M2 -: M pivotat 2 2
   Mpat=. [
   Rpat=. {~ {.
   Ppat=. {~ <
   Cpat=. {"1~ {:
   Bpat=. {.@] = i.@#@[
   pivotat=: (Mpat - (Cpat-Bpat) */ Rpat % Ppat) f.
   assert M2 -: M pivotat 2 2

Note that this expression works like pivot~ defined from '~system/packages/math/matutil.ijs. (pivot has a swapped order for its arguments, which is why I have used a different name, here.)

   pivotat
[ - (({"1~ {:) - {:@] = i.@#@[) */ ({~ {.) % ({~ <)
   M pivotat 2 2
1 _0.8 0 0  0.6 9
0  3.6 0 1 _0.2 9
0  0.4 1 0  0.2 3

I have elected to retain Richard Brown's argument order, in my pivotat, rather than the matutil.ijs argument order because Richard Brown's argument order lends itself nicely to a heavy use of indexing forks.

Note, however, that I have started to deviate from Richard Brown's treatment of this subject. His essay was about teaching students the simplex method, and he goes into quite a bit of detail on the concept of picking the row/column indices. He also uses indices where 1 selects the first element, which is not how J works (in J, 0 picks the first element).

Anyways, you pick the pivot column by excluding the last column of the matrix and then finding the column with the lowest value in the first row. You then find the pivot row by finding the row with the largest value in that pivot column. In other words:

pivM=. }:"1
pivCol=. (i. <./)@{.
pivRow=. [: (i. >./) {"1~
pivotinds=: ([: ((pivRow, ]) pivCol) pivM) f.
assert 2 2 -: pivotinds M

And that, it seems to me, is enough to define simplex:

splex=: (pivotat pivotinds)^:_
   splex M
1 0 0  0.222222   0.555556  11
0 1 0  0.277778 _0.0555556 2.5
0 0 1 _0.111111   0.222222   2

Note that I did not take any special steps to detect when I had found the optimum values. It seems to me that that is adequately handled by the mechanics of ^:. Logically speaking, however, this state has occurred when the value from the first row in the pivot column is zero.

As Richard Brown noted on his page, this means that

 z      + 0.22s + 0.55t = 11
   x    + 0.27s - 0.55t = 2.5
      y - 0.11s + 0.22t = 2

Here, the last column represents the solution, as we set s and t both to zero (remember that originally we let them take on any non-negative values). 11 is the value for "z", 2.5 is the value for "x" and 2 is the value for "y".

   2.5 z 2
11
   2.5 c1 2
1
   2.5 c2 2
1

When I test my code against the examples from simplexnr.ijs, I get the same result

   simplexnr 2 3;(4 1,:2 5);10 15;_1 _1
+-+-----+--+
|0|2.5 2|11|
+-+-----+--+

further links

You might have better luck searching the J Forums directly, as in linear program, linear programming, and simplex. For example, I was able to find a simplex script on on Keith Smillie's J page.

And a quote of historical interest:

The third page had an illustration that, in a few short lines, described George Dantzig’s simplex algorithm simply and precisely.

That was the overwhelming, crucial experience.

-- Michael Montalbano, A Personal History of APL, 1982, quoted in the Ken Iverson anecdote page.

---