How to use this notebook

This solves a Linear Programming problem of the form
min  c x
          A x= b
          x ≥ 0
via a primal-dual method, i.e. Newton iterations of the complementarity equations.

Routines

showiteration[xx_,lastpoints_]:=Block[{xxx},
  xxx=Take[xx,-Min[lastpoints,Length[xx]]];
 iterlines0=Graphics[{RGBColor[1,0.,0.],Line[xxx]},
           PlotRange->  box[xxx]];
  iterpoints0=Graphics[Table[{PointSize[0.005+0.015/ i],Point[xxx[[i]]]},
         {i,Length[xxx]}],
         PlotRange-> box[xxx]];
  ax=Graphics[{RGBColor[0.6,0.2,0.4]},Axes->Automatic,
            AxesStyle->{RGBColor[0,0.7,1]},
             PlotRange-> box[xxx]];
  Show[ax,iterlines0,iterpoints0,PlotRange->All];

        
]

CompositeMatrix[M_]:=Block[{n,m},
    n=Length[M[[1]]];
    m=Length[M];
    Apply[Join,Table[Transpose[
    Apply[Join,Table[Transpose[M[[i,j]]]
       ,{j,n}]]],{i,m}]]
]

zmat[n_, m_] := Table[0, {n}, {m}]

zeromatrix[m_,n_]:=Table[0,{m},{n}]

box[points_]:=Block[
{pp,xmax,xmin,ymax,ymin,xsize,ysize,maxratio,center},
     pp=Transpose[points];
     xmax=Max[pp[[1]]];
     xmin=Min[pp[[1]]];
     ymax=Max[pp[[2]]];
     ymin=Min[pp[[2]]];
     xsize=xmax-xmin;
     ysize=ymax-ymin;
     maxratio=2;
     If[xsize > maxratio ysize, 
        ysize = xsize/maxratio];
     If[ysize > maxratio xsize, 
        xsize = ysize/maxratio];
     center=N[{xmax+xmin,ymax+ymin}/2];
     center+{{-1,1}xsize,{-1,1}ysize}/2 1.1
];

Problem Data

Instance #1

A={{2,1,3,1},{0,1,2,1}}

c={1,-2,3,-1};

b={7,4};

n=Length[A[[1]]]

m=Length[A]

Instance #2  Assignment

This is a random nn x nn weighted assignment instance

nn=10;

A=Table[0,{i,2 nn},{j,nn^2}];

Do[A[[i,(j-1)nn+i]]=1;A[[nn+j,(j-1)nn+i]]=1,
{i,nn},{j,nn}]; A=Drop[A,-1]

A

b=Table[1,{2 nn-1}]

SeedRandom[10]

c=Table[Random[Integer,{1,9}],{nn^2}]

nn

Perturbation

c = 
      c  + Table[Random[], {nn^2}]/10

n=Length[A[[1]]];
m=Length[A];

Do[If[Length[A[[j]]]!=n,Print["Check row ",j," of A !"]],
    {j,m}];
If[Length[b]!=m,Print["check b !"]];
If[Length[c]!=n, Print["check c !"]];

Manual iteration

Initialization

set initial values for x and s (must be strictly positive !) and y (no constraints for y)

x = Table[1, {n}]; s = x;
y = Table[0, {m}];

set an initial (fairly large) value for tau

tau = 1024;

this is the Jacobian

CompositeMatrix[{
{Table[0,{n},{n}],Transpose[A],IdentityMatrix[n]},
{A,Table[0,{m},{m}],Table[0,{m},{n}]},
{DiagonalMatrix[s],Table[0,{n},{m}],DiagonalMatrix[x]}}]

fulliter={Join[{tau},x,y,s,{0}]};partiter={Join[{tau}, x,y,s ]}

Iteration

computation of Newton direction

sol=LinearSolve[CompositeMatrix[{
{Table[0,{n},{n}],Transpose[A],IdentityMatrix[n]},
{A,Table[0,{m},{m}],Table[0,{m},{n}]},
{DiagonalMatrix[s],Table[0,{n},{m}],DiagonalMatrix[x]}}],
N[Join[c-y.A-s,b-A.x,-x s +tau]]]

new values

xx=x+Take[sol,n]

yy=y+Take[sol,{n+1,n+m}]

ss=s+Take[sol,-n]

-x s +tau

gap=Max[Abs[-xx ss+tau]]

AppendTo[fulliter,Join[{tau},xx,yy,ss,{gap}]];

x=xx;s=ss;y=yy;

If gap is sufficiently small and x and s are non negative halve tau and redo previous computations

AppendTo[partiter, Join[{tau}, xx, yy, ss ]];

tau=0.5 tau

Automatic iteration

x=Table[1,{n}];s=x;
tau=1024;y=Table[0,{m}];
iter={Join[{"*"},x,y,s,
{Max[Max[Abs[A.x-b]],Max[Abs[y.A+s-c]]]},
{Max[Abs[-x s +tau]]}]};
tol=0.0001;
While[tau > 0.001,
   gap=1;cv=1;
   While[Max[gap,cv] > tol,
         sol=LinearSolve[CompositeMatrix[{
         {Table[0,{n},{n}],Transpose[A],IdentityMatrix[n]},
         {A,Table[0,{m},{m}],Table[0,{m},{n}]},
         {DiagonalMatrix[s],Table[0,{n},{m}],DiagonalMatrix[x]}}],
         N[Join[c-y.A-s,b-A.x,-x s +tau]]];
         xx=x+Take[sol,n];
         yy=y+Take[sol,{n+1,n+m}];
         ss=s+Take[sol,-n];
         gap=Max[Abs[-xx ss+tau]]/xx.ss;
         cv=Max[Max[Abs[A.x-b]],Max[Abs[y.A+s-c]]];
         curr=Join[{tau},xx,yy,ss,{cv},{gap}];
         AppendTo[iter,curr];
         x=xx;s=ss;y=yy;
        ];
      tau=0.5 tau; Print["tau= ",tau];
]

Chopping small elements

delta = 0.01

xch = Chop[x, delta]

sch = Chop[s, delta]

Selecting positive elements

xnz = 
      xch[[Complement[Range[Length[x]], 
            Flatten[Position[xch, 0]]]]]

Rationalize[xnz, 0.06]

snz = 
      sch[[Complement[Range[Length[s]], 
            Flatten[Position[sch, 0]]]]]

Rationalize[snz, 0.1]

Recovering solution  via least squares

Dual system

Anz = 
      Transpose[Transpose[
            A][[Complement[
              Range[Length[s]], 
       Flatten[
                Position[sch, 0]]]]]]

cnz = 
      c[[Complement[Range[Length[s]], 
            Flatten[Position[sch, 0]]]]]

Az = 
      Transpose[Transpose[A][[Flatten[Position[sch, 0]]]]]

cz = c[[Flatten[Position[sch, 0]]]]

n1 = Length[Anz]

n2 = Length[Az[[1]]]

n3 = Length[Anz[[1]]]

specmat = 
      CompositeMatrix[{{zmat[n1, n1], Az, 
            Anz}, 
{Transpose[Az], zmat[n2, n2], 
            zmat[n2, n3]}, 
{Transpose[Anz], zmat[n3, n2], 
            IdentityMatrix[n3]}}]

ssoll = LinearSolve[specmat, Join[Anz . snz, cz, cnz]]

finaly = Take[ssoll, {1, n1}]

finals = cnz - finaly . Anz

raty = Rationalize[finaly,  .03]

raty . Az

cz

rats = cnz - raty . Anz

posvecs = 
      Complement[Range[Length[s]], Flatten[Position[sch, 0]]]

trues = Table[0, {Length[s]}];
Do[trues[[posvecs[[i]]]] = 
        rats[[i]]
 , {i, Length[posvecs]}]

Primal system

Anzx = 
      Transpose[Transpose[
            A][[Complement[
              Range[Length[x]], 
       Flatten[
                Position[xch, 0]]]]]]

z = LinearSolve[Anzx . Transpose[Anzx], Anzx . xnz - b]

xxx = xnz - Transpose[Anzx] . z

posvec = 
      Complement[Range[Length[x]], Flatten[Position[xch, 0]]]

truex = Table[0, {Length[x]}];
Do[truex[[posvec[[i]]]] = 
        xxx[[i]]
 , {i, Length[posvec]}]

truex

A . truex - b

trues . truex

raty . A + trues - c

Iteration Display

Select two primal variable indeces (between 1 and n)

n

firstindex = 1

secondindex = 2

Iteration history for the selected variables

L11 are all computed values of first selected variable

L11 = Transpose[fulliter][[1 + firstindex]]

ListPlot[L11, PlotJoined -> True]

L12 are  the values  of first selected variable before halving tau

L12 = Transpose[partiter][[1 + firstindex]]

ListPlot[L12, PlotJoined -> True]

same for second selected variable

L21 = Transpose[fulliter][[1 + secondindex]];

L22 = Transpose[partiter][[1 + secondindex]];

ListPlot[L21, PlotJoined -> True]

ListPlot[L22, PlotJoined -> True]

bidimensional display of selected variables

itt1 = Transpose[{L11, L21}]

itt2 = Transpose[{L12, L22}]

the second parameter in showiteration is the number of points displayed starting from the end

showiteration[itt1, 30]

showiteration[itt2, 30]