While running experiments on a simple LP model [1], I wanted to see what Simplex solvers returned when they hit an iteration limit early on. Depending on the iteration limit, the model status could be
- intermediate infeasible,
- intermediate nonoptimal,
- or optimal.
GAMS Model
$ontext LAD (Least Absolute Deviation) Regression LP model using variable splitting $offtext
options lp=xpress, seed=1234, iterlim=10;
*--------------------------------------------------------------------- * data *---------------------------------------------------------------------
set i 'observations' /observation1*observation30/;
parameter data(i,*);
data(i,'x') = uniform(0,10); data(i,'y') = 10+0.25*data(i,'x')+normal(0,3); * for a few observations, generate more wild y data data(i,'y')$(uniform(0,1)<0.1) = 10+0.25*data(i,'x')+normal(0,10); option data:8; display data;
*--------------------------------------------------------------------- * derived data *---------------------------------------------------------------------
parameters x(i) 'independent variable' y(i) 'dependent variable' ; x(i) = data(i,'x'); y(i) = data(i,'y');
*--------------------------------------------------------------------- * LAD regression LP model *---------------------------------------------------------------------
set pm 'plus/minus' /'+','-'/;
variables a 'intercept' b 'slope' z 'sum of absolute deviations' ; positive variables r(i,pm) 'residuals' ;
equations objective 'objective function' fit(i) 'model to be fitted' ;
objective.. z =e= sum((i,pm), r(i,pm)); fit(i).. y(i) =e= a + b * x(i) + r(i,'+') - r(i,'-');
model lad /all/; solve lad using lp minimizing z; display a.l, b.l, r.l, z.l; |
S O L V E S U M M A R Y
MODEL lad OBJECTIVE z
TYPE LP DIRECTION MINIMIZE
SOLVER XPRESS FROM LINE 64
**** SOLVER STATUS 2 Iteration Interrupt
**** MODEL STATUS 14 No Solution Returned
**** OBJECTIVE VALUE 0.0000
RESOURCE USAGE, LIMIT 0.032 10000000000.000
ITERATION COUNT, LIMIT 10 10
No solution returned
---- 65 VARIABLE a.L = 0.000 intercept
VARIABLE b.L = 0.000 slope
---- 65 VARIABLE r.L residuals
( ALL 0.000 )
---- 65 VARIABLE z.L = 0.000 sum of absolute deviations
From the log we have:
Starting concurrent solve with dual (1 thread) Concurrent-Solve, 0s Dual objective dual inf stopped ----- interrupted ------ Concurrent statistics: Dual: 10 simplex iterations, 0.00s Maximum number of iterations reached ... XPRESS stop status: iteration limit (4) ... XPRESS sol status: not available (0) ... Status: Iteration Limit (2) / No Solution (14) ... Solved in 0.032s, 10 iterations
Obviously, the solver stopped nicely after 10 iterations, but unfortunately did not give us a solution we could study. The GAMS people said this was an Xpress problem; nothing they can do about it. So, that piqued my interest. Why would Xpress withhold evidence? Is this something we can replicate?
Xpress/Python replication
My hypothesis is that the GAMS link is using problem.getSolution() which fails when there is no optimal solution. I.e., using:
prob.setControl("LPITERLIMIT", 10) solvestatus, solstatus = prob.optimize() try: sol = prob.getSolution() unpackSolution() except: print("Get lost! I don't want to give you the current solution.")
will show: Get lost! I don't want to give you the current solution.
Complete Python Code
import xpress as xp from io import StringIO import pandas as pd iterlimit = 10 #---------------------------------------------- # data #---------------------------------------------- df = pd.read_table(StringIO(""" i x y observation1 0.44210132 6.44750203 observation2 3.82274024 10.01983898 observation3 6.23060668 7.78926079 observation4 5.72676411 9.70417381 observation5 4.60680179 12.42981755 observation6 8.77921471 11.61833643 observation7 9.27713358 14.67683786 observation8 5.93381615 28.95514017 observation9 2.72686678 6.06008592 observation10 8.38980378 15.84345556 observation11 7.66500308 11.44049644 observation12 5.79829196 8.44430907 observation13 9.09576650 7.25589802 observation14 7.77304183 11.28781064 observation15 4.45932581 7.59624812 observation16 1.59348806 6.75018412 observation17 3.95136687 7.49105201 observation18 2.77230871 11.99033842 observation19 1.75187971 11.37407422 observation20 4.91071592 11.55804245 observation21 5.61335506 12.91631670 observation22 6.16796837 11.60374255 observation23 7.55066563 11.72576766 observation24 5.05663151 6.75010567 observation25 2.54691604 6.26558134 observation26 8.15682045 11.05210464 observation27 2.39781739 8.13913190 observation28 5.16244536 10.73444211 observation29 8.78656221 12.71992827 observation30 2.94928071 11.67903419 """), sep="\\s+") x = df['x'].to_numpy() y = df['y'].to_numpy() #---------------------------------------------- # LP Model #---------------------------------------------- n = len(x) p = xp.problem() rplus = p.addVariables(n, name='rplus') rmin = p.addVariables(n, name='rmin') a = p.addVariable('a', lb=-xp.infinity) b = p.addVariable('b', lb=-xp.infinity) p.addConstraint(rplus[i] - rmin[i] == y[i] - (a + b * x[i]) for i in range(n)) p.setObjective(sum([rplus[i] + rmin[i] for i in range(n)]), sense=xp.ObjSense.MINIMIZE) #---------------------------------------------- # Solve and reporting #---------------------------------------------- p.setControl("LPITERLIMIT", iterlimit) solvestatus, solstatus = p.optimize() print(f"{solvestatus=}, {solstatus=}") print(f"Num infeas: {p.attributes.primalinfeas}") print(f"Num nonopt: {p.attributes.dualinfeas}") print(f"Objective value: {p.attributes.lpobjval:.3f}") try: sol = p.getSolution() print(f"Estimates: a = {sol[2*n]:.3f}, b = {sol[2*n+1]:.3f}") print("Residuals (rplus, rmin):") dfsol = pd.DataFrame({'rplus': sol[0:n], 'rmin': sol[n:2*n]}) print(dfsol.round(3).astype(str).replace(['0.00','0.0','0'], '')) except: print("Get lost! I don't want to give you the current solution.")
The output looks like:
Starting concurrent solve with dual (1 thread) Concurrent-Solve, 0s Dual objective dual inf stopped ----- interrupted ------ Concurrent statistics: Dual: 10 simplex iterations, 0.00s Maximum number of iterations reached solvestatus=<SolveStatus.STOPPED: 1>, solstatus=<SolStatus.NOTFOUND: 0> Num infeas: 20 Num nonopt: 0 Objective value: 36.198 Get lost! I don't want to give you the current solution.
The Fix
prob.setControl("LPITERLIMIT", 10) solvestatus, solstatus = prob.optimize() sol = prob.getLpSol()[0]
Complete Python Code
import xpress as xp from io import StringIO import pandas as pd iterlimit = 10 #---------------------------------------------- # data #---------------------------------------------- df = pd.read_table(StringIO(""" i x y observation1 0.44210132 6.44750203 observation2 3.82274024 10.01983898 observation3 6.23060668 7.78926079 observation4 5.72676411 9.70417381 observation5 4.60680179 12.42981755 observation6 8.77921471 11.61833643 observation7 9.27713358 14.67683786 observation8 5.93381615 28.95514017 observation9 2.72686678 6.06008592 observation10 8.38980378 15.84345556 observation11 7.66500308 11.44049644 observation12 5.79829196 8.44430907 observation13 9.09576650 7.25589802 observation14 7.77304183 11.28781064 observation15 4.45932581 7.59624812 observation16 1.59348806 6.75018412 observation17 3.95136687 7.49105201 observation18 2.77230871 11.99033842 observation19 1.75187971 11.37407422 observation20 4.91071592 11.55804245 observation21 5.61335506 12.91631670 observation22 6.16796837 11.60374255 observation23 7.55066563 11.72576766 observation24 5.05663151 6.75010567 observation25 2.54691604 6.26558134 observation26 8.15682045 11.05210464 observation27 2.39781739 8.13913190 observation28 5.16244536 10.73444211 observation29 8.78656221 12.71992827 observation30 2.94928071 11.67903419 """), sep="\\s+") x = df['x'].to_numpy() y = df['y'].to_numpy() #---------------------------------------------- # LP Model #---------------------------------------------- n = len(x) p = xp.problem() rplus = p.addVariables(n, name='rplus') rmin = p.addVariables(n, name='rmin') a = p.addVariable('a', lb=-xp.infinity) b = p.addVariable('b', lb=-xp.infinity) p.addConstraint(rplus[i] - rmin[i] == y[i] - (a + b * x[i]) for i in range(n)) p.setObjective(sum([rplus[i] + rmin[i] for i in range(n)]), sense=xp.ObjSense.MINIMIZE) #---------------------------------------------- # Solve and reporting #---------------------------------------------- p.setControl("LPITERLIMIT", iterlimit) solvestatus, solstatus = p.optimize() print(f"{solvestatus=}, {solstatus=}") print(f"Num infeas: {p.attributes.primalinfeas}") print(f"Num nonopt: {p.attributes.dualinfeas}") print(f"lpstatus: {p.attributes.lpstatus.name}") sol = p.getLpSol()[0] dfsol = pd.DataFrame({'rplus': sol[0:n], 'rmin': sol[n:2*n]}) print(f"Objective value: {p.attributes.lpobjval:.3f}") print(f"Estimates: a = {sol[2*n]:.3f}, b = {sol[2*n+1]:.3f}") print("Residuals (rplus, rmin):") print(dfsol.round(3).astype(str).replace(['0.00','0.0','0'], ''))
The output looks like:
Starting concurrent solve with dual (1 thread) Concurrent-Solve, 0s Dual objective dual inf stopped ----- interrupted ------ Concurrent statistics: Dual: 10 simplex iterations, 0.00s Maximum number of iterations reached solvestatus=<SolveStatus.STOPPED: 1>, solstatus=<SolStatus.NOTFOUND: 0> Num infeas: 20 Num nonopt: 0 lpstatus: UNFINISHED Objective value: 36.198 Estimates: a = 5.812, b = 1.436 Residuals (rplus, rmin): rplus rmin 0 1 2 3 4 5 6 7 14.619 8 3.669 9 10 11 12 13 5.69 14 15 1.351 16 3.997 17 2.196 18 3.045 19 20 21 22 23 24 25 26 27 28 29 1.63
Conclusion
References
- A note on variable splitting: complementarity. https://yetanothermathprogrammingconsultant.blogspot.com/2026/11/a-note-on-variable-splitting.html
No comments:
Post a Comment