Simplex Algorithm

[Albert Chan]

(11-12-2023 12:38 AM)Albert Chan Wrote:  Fastest is remove excess variables, turning it into slack variables.

Note that we scale pivot as 1 before clearing other cells.
This does not affect other constraints, not even cost functions.
In other words, row operations are safe.

Time (including pivot operations) = 0.0075 s
Production code need pivots selection, but should not cost much.

I don't know a good name for this algorithm, so I just use remove_eq_c(a, n)
It remove top n "=" constraints. (n variables, if constraints are independent)
In other words, only "≤" constraints remain.

simplex_le(a, ±n) ≡ simplex_le(remove_eq_c(a,n), ±∞)

simplex_le(a, dir) only had 1 change. Instead of more constraints, it now do less variables.

< IF b<len(a) THEN a := extend(-a[1 .. b],a) END
> IF b<len(a) THEN a := remove_eq_c(a, b) END

Code:

pivot_index(a) :=
BEGIN    
    LOCAL b := contains((a:=abs(a)), 1.0);
    IF b THEN RETURN b END
    IF (b:=max(a)) THEN b := contains(a,b) END
    RETURN b;
END;

remove_eq_c(a, n) :=
BEGIN
  LOCAL p, r, c := len(a[1]);
  FOR r FROM n DOWNTO 1 DO
    p := pivot_index(suppress(a[r],c));
    IF p==0 and a[r][c] THEN error(a) END
    a := p ? pivot((a[r]/=a[r][p]),r,p) : delrows(a,r);
  END
  return a;
END;

simplex_le(a, dir) :=
BEGIN
 LOCAL b := abs(dir);
 IF b<len(a) THEN a := remove_eq_c(a, b) END
 b := dim(a);
 a := transpose(extend(
    [col(a,1 .. b(2)-1)],
    [col(identity(b(1)),1 .. b(1)-1)],
    [col(a,b(2))]
 ));
 IF dir>0 THEN a[0] := -a[0] END /* maximize */
 a := simplex2(a,range(b(2),sum(b)-1),sum(b)-2,0,0);
 IF dir>0 THEN a[1] := -a[1] END /* undo flip */
 return a;
END;

Benchmark example:

> m := [[3,2,1,10], [2,5,3,15], [-2,-3,-4,0]]
> time(r := simplex_le(m, -2))                     → 0.0077
> r

\([\displaystyle
\frac{-130}{7},[\frac{15}{7},0,\frac{25}{7},\frac{15}{7},0],\left(\begin{array}{cccccc}
1 & \frac{1}{7} & 0 & 1 & \frac{1}{7} & \frac{15}{7} \\
0 & \frac{11}{7} & 1 & 0 & \frac{11}{7} & \frac{25}{7} \\
0 & \frac{25}{7} & 0 & 0 & \frac{25}{7} & \frac{130}{7}
\end{array}\right) ]\)

Update Nov 12, 2023: simplex_le directly call simplex2, instead of wrapper simplex.
see https://www.hpmuseum.org/forum/thread-20...#pid180152