Thanks for trying this @andrewdelong.
I'm not really sure what to do here. Do you agree that the code accurately reflects the algorithm as described in reference [4]? Does it not seem like an inherent issue in the algorithm?

@mdhaber I'm a little puzzled too.

Regarding the original bug, that's a case where the problem is unambiguously unbounded, so I think it should definitely return unbounded in that case. The bug only shows up with tol=1e-8 (the default), so it's somewhat pathological and could be an unlucky accumulation of numerical errors (b.T @ y - tol is only around 1.5e-8). I am not sure.

For the failing test, that's the case where the system is unbounded wrt one variable and infeasible wrt another. In this case there doesn't seem to be a consistent relation between b.T @ y and -c.T @ x, so I am not sure how to detect infeasibility separately from unboundedness.

For example, the system below returns "unbounded", both for the original code and proposed changes:

c = [-1, -1]
A_eq = [[0, 1], [0, 2]]
b_eq = [1, 1]
linprog(c, None, None, A_eq, b_eq, options={"presolve": False})

because at termination it has homogeneous quantities

 b.T @ y = -0.520742401565883 < tol
-c.T @ x = 1.3264191994902952
 x = [1.32641920e+00 1.95971331e-11]
 y = [-0.1735808 -0.3471616]
 z = [1.39314638e-13 8.67904003e-01]
 tau = 2.4461587693169447e-11
 kappa = 0.8056767979120855

but merely changing to A_eq = [[0, 1], [0, 3]] returns "infeasible" because it terminates with completely different values:

 b.T @ y = 1.088743805719052e+37 > tol
-c.T @ x = 2.0511091275047628e+36
 x = [2.05110913e+36 2.93509133e-52]
 y = [ 8.46344192e+41 -8.46333304e+41]
 z = [6.56465468e-54 1.69266213e+42]
 tau = 3.34693300342076e-57
 kappa = 2.0523495828590013e+36

Negative coefficients like A_eq = [[0, 1], [0, -2]] (like those in test_enzo_example_c_with_infeasibility test) consistently return "infeasible", probably because of the default bounds:

 b.T @ y = 48.1936206703016 > tol
-c.T @ x = 0.9116045616117406
 x = [9.11604562e-01 1.78034291e-12]
 y = [31.80237933 16.39124134]
 z = [1.69922691e-12 9.80103352e-01]
 tau = 1.8596743291897465e-12
 kappa = 1.1076188140147998

In any of these cases, the [x, tau] and [z, kappa] vectors at termination are strictly complementary, with tau near zero and kappa very large, so the conditions of Theorem 8.3 seem satisfied. There also don't seem to be any relative or absolute tolerances being satisfied between b.T @ y - c.T @ x and kappa in these bad cases.

Still, if you run the above examples with presolve=True then linprog will return "trivially unbounded" due to the column of zeros. This isn't consistent with what test_enzo_example_c_with_infeasibility seems to want: it uses presolve=False and expects "infeasible" to be returned, despite "unbounded" being returned for m <= 3.

so it's somewhat pathological and could be an unlucky accumulation of numerical errors

I'm tempted to leave this alone unless there is a really sound way to resolve it. I think this part of the code faithfully corresponds with the algorithm, and we can't guarantee that the code will work on every problem with the default tolerances. Does that make sense? I appreciate your thoughts on this. I don't really have any other ideas, though.

I think it's totally reasonable not to make the change, so no worries closing it out. I'll just add one last comment on why I lean slightly toward the change. There are two cases under consideration:

1. An "obviously unbounded" case where, though rare, the current code incorrectly returns "infeasible."
2. A "unbounded in some variables, infeasible in others" case where the current code is biased towards "infeasible" but isn't consistent (can return "unbounded", depending on coefficients).

The PR fixes (1) but alters the already-inconsistent behaviour of (2) by removing the "bias" and thereby the pretence of consistency in scenario (2). There exists a unit test concerning (2), so someone out there considers "infeasible" to be [always?] correct, matching simplex, but the current code doesn't achieve that consistently for interior-point.

As for fidelity to the MOSEK paper, I didn't encounter any specifics in [4] about how to test b.T @ y > 0 versus c.T @ x < 0 in practice, so it's hard for me to say. The current code seems to place more "trust" in testing b.T @ y > tol than in testing c.T @ x < -tol, and I guess this was an example where that trust breaks down.

@Kai-Striega I'd appreciate a third opinion on this. I wouldn't ask you to delve into the specifics, just the broad strokes - linprog interior point is currently incorrect for a very simple problem (without presolve), this PR fixes simple problem but breaks another problem in the test suite. There are reasonable arguments for and against the change. Can you weigh in?

@andrewdelong One thing that would help convince me would be benchmarking on infeasible and unbounded problems. I've wanted to add the infeasible problems from the Netlib LP test suite to our benchmarks. Would you be interested in taking a shot at that? (If you need help getting started, I think I can point you in the right direction.)

@mdhaber I could but not in the coming weeks. Is your goal to evaluate speed (scipy/benchmarks) or correctness (scipy/optimize/tests)? And in a separate PR I assume. I see you've converted NETLIB instances to (c, A_ub, b_ub, A_eq, b_eq, bound) format and baked those arrays into npz, so if you have a consistent/scripted way of doing that conversion let me know.

Of course, if scipy provided MPS file reading, like cvxopt's fromfile, then users (and tests) could load NETLIB and other files directly to scipy, but that's a whole other feature and not much advantage.

The goal is to evaluate correctness in classifying them as infeasible or unbounded, not to measure speed.

Yes, I did do those conversions using my own MPS reader a few years ago. (I'm not sure if cvxopt's was available, or maybe I just didn't find it.) I just put it up on a repo for you after converting to Python 3 : )
https://github.com/mdhaber/mpsreader

I tested using the script at the bottom that it loads uncompressed MPS files correctly (unless they have ranges - that's not implemented) and saves and loads .npz files. I think you can pretty easily extend the script to run a batch job.

Note that the files you download from Netlib are compressed MPS files. They can be uncompressed using emps.

Whenever you get around to it, thanks!

First, welcome @andrewdelong and thank you taking the time and effort, it's made understanding the issue and PR much easier for me. Second my apologies for the delayed response - I've been quite ill over the last month and haven't focused too much on FOSS.

I am slightly for the change, but I would prefer to see the results of the benchmarks before we commit to a change. Such a simple problem shouldn't be failing. Especially if the alternative test catches it reliably and the existing code is inconsistent to begin with.

I'm tempted to leave this alone unless there is a really sound way to resolve it. I think this part of the code faithfully corresponds with the algorithm, and we can't guarantee that the code will work on every problem with the default tolerances.

If I understand it correctly, @mdhaber's main concern with changing it isn't the inconsistent behaviour of the existing code, but that this issue essentially stems from numerical errors which could always occur if the problem is picked correctly. But the proposed changes look solid and follow the paper (or at least I couldn't find a reference on which method to use).

One more thing to consider, what if someone currently relies on having a specific tolerance and being able to change it via the tol parameter?

If I understand it correctly, @mdhaber's main concern with changing...

Maybe my concern can be addressed by the answer to this question:

If we get to this comparison in the code:
if b.transpose().dot(y) > -c.transpose().dot(x):
Is it guaranteed that at least one of:
b.transpose().dot(y) > 0
or
c.transpose().dot(x) < 0
is true?

@Kai-Striega Thanks! I don't know when I could contribute to the benchmarks. This PR can go on the backlog if everyone's comfortable.

@mdhaber Good question but I don't know the answer. It's empirically true for the handful of instances I inspected. I suppose it's true to the degree that the gap constraint b.T @ y - c.T @ x = kappa is satisfied at the moment inf1 or inf2 becomes true, since the algorithm at least seems to successfully maintain kappa >= 0 throughout.

Alternative fix that is probably a bad idea. I noticed that reference [4] uses different internal tolerances than linprog (see Table 8.1). I tried using these different tolerances (while also reverting back to b.T @ y > tol) and this fixes bug #11617 and allows test_enzo_example_c_with_infeasibility to pass. See andrewdelong#1 for the small change. However, it increases the effective default precision of linprog and so some tests that rely on fast solvers (cholesky, sym_pos) now fail. Even if those failing tests could be fixed, I do not think the "use tolerances from [4]" approach is a good fix at this point:

• it has a broad effect on the default behaviour of linprog even for feasible/bounded cases;
• it might slow down linprog due to extra iterations;
• there's no reason the thresholds from [4] would fundamentally rule out bugs like linprog returns 'infeasible' for simple unbounded case #11617, it may just be a serendipitous fix.

I only mention the Table 8.1 tolerances to emphasize a meaningful difference between linprog and what Anderson & Anderson implemented.

I only mention the Table 8.1 tolerances to emphasize a meaningful difference between linprog and what Anderson & Anderson implemented.

Yes, I remember reviewing the table when I wrote the method, and I believe we agreed to leave separating and exposing all the tolerances to the user in a future PR.

It's OK with me for this PR to go on the backburner. Would you like me to check in with you in a few months?

👍 for this going on the backlog (for now).

Maybe, once the lapack wrappers are done, I can take a look at setting up the infeasible benchmarks if no one else has gotten around to them by that point.

Update. The previous version of this PR (test b.T @ y > -c.T @ x) failed to fix several new cases I added to the test. The crux of the original bug seems to be b.T @ y not satisfying any absolute tolerance (mathematically, not float error). I've updated the PR with a fix that seems more justified.

Old fix fails. The original code (b.T @ y > tol) and the old fix both fail if the LP is scaled up by 1e9:

s = 1e9  # or 1e4, ... 1e9
c = [s, 0., -s]
A_eq = [[s, s, -s]]
b_eq = [-s]
linprog(c, None, None, A_eq, b_eq)

terminates due to inf2 (not inf1) after two iterations with:

 b.T @ y = 0.9997708945929584   <-- b.T @ y > -c.T @ x would return infeasible (WRONG)
-c.T @ x = -0.0002696514129638672
 x = [9.9999999999986522e-01 3.5023982751083223e-09 1.0000000000001348e+00]
 y = [-9.997708945929584e-10]
 z = [3.4997708950917548e-09 9.9977089809734287e-01 3.4997708940888678e-09]
 tau = 2.4999999999987055e-09
 kappa = 1.0000399904359316

The above also indicates that b.T @ y doesn't satisfy any absolute tolerance, which makes sense.

New fix. When b.T @ y > 0, the proof of infeasibility (see Ye et al. 1994, Theorem 3(ii)) relies on complementary slackness, so I tried using x.T @ z as the tolerance (related to mu). This worked in every case that terminates due to inf1 (including the original buggy LP where s=1 and tol didn't work), but there were still cases with large coefficients that terminated due to inf2 and failed this condition (including the example above). So, I added condition A.T @ y <= 0 based on Farkas' lemma as described in Lemma 2.1 of Anderson 2001. All current tests now pass.

Thoughts on alternatives. I'm not sure the new conditions are perfect, but they're somewhat justified and work so far. I can't think of a legitimate role for tol because all the quantities it guarantees at termination (rho_p, rho_d, rho_mu, etc) are just measurements of relative reductions of residuals (A @ x - b * tau etc) from the blind-start residuals; I am not sure those relative reductions in residuals tell us anything about the absolute or relative accuracy of b.T @ y or other relevant terms (A.T @ y, etc), as compared to their unknown optimal values.

@andrewdelong I'll try to take a look at your most recent work this week
@Kai-Striega I'm working on the infeasible benchmarks right now. If you find unbounded benchmarks, let me know. (I haven't found any yet. I supposed the duals of the infeasible ones could be unbounded, but we'd need ground truth on that.)

Update: currently linprog correctly reports infeasible for all the Netlib infeasible benchmarks except cpex1, cplex2, gosh, and greanbea with dense arithmetic. greanbea doesn't seem to load properly, so I couldn't test it. gosh took more than a few minutes and was running into numerical issues, so I stopped it. Incorrect status on cplex1 (status 3 - unbounded) and cplex2 (status 0 - success) are not so surprising based on the descriptions:

CPLEX1, CPLEX2:  medium and large problems respectively.  CPLEX1
referred to as CPLEX problem in Chinneck [1993], and is remarkably
non-volatile, showing a single small IIS regardless of the IIS algorithm
applied.  CPLEX2 is an almost-feasible problem. Contributor:  Ed Klotz,
CPLEX Optimization Inc.

With sparse arithmetic, linprog correctly reports that gosh is infeasible. It has the same trouble with the cplex problems and has some new issues: we get numerical difficulties / status 4 with pang, and unbounded (status 1 - iteration limit reached) on pilot4i. The issues with pang seems to be caused by presolve; without presolve, linprog correctly reports status 2.

@andrewdelong Yes, linprog doesn't consider problem scaling very carefully.

I see that the new approach has justification. I'm curious how it does on the benchmarks I just added in gh-11729. If we can also find a few unbounded benchmarks and test on those, too, I'd be happy to merge this. Do you have access to Matlab? I thought we could check the duals of these infeasible benchmarks for unboundedness, and use the ones that are unbounded.

@mdhaber Thanks for #11729, as this PR was indeed breaking a couple of the tests you added. The only way I could get it to pass all tests (unit tests + NETLIB tests) was a two-tiered logic: first check for for well-behaved infeasible/unbounded cases (the first if, elif), then use A.T @ y <= 0 to catch some poorly-behaved infeasible cases (the second elif). I couldn't find something simpler that worked, unfortunately.

Separate infeasibility detection issue. Trivially infeasible instances with small coefficients cause linprog to raise an uncaught ValueError, claiming the instance was easy to solve:

ValueError: The algorithm terminated successfully and determined that the problem is infeasible.

This happens for simplex and interior-point, caused by tol checks in _linprog_util._check_result. This error came up when I scaled coefficients of NETLIB tests ex72a and ex73a by 1e-5, but can be triggered by the example test below. Not sure what the intended behaviour is.

    def test_bug_postprocess_raises_trivial_infeas(self):
        # postprocess raised an uncaught ValueError for a trivially infeasible
        # instance with small coefficients
        c = [1e-5]
        A_eq = [[1e-5]]
        b_eq = [0]
        bounds = [(1e-5, None)]
        res = linprog(c, None, None, A_eq, b_eq, bounds, method=self.method)
        _assert_infeasible(res)

I can't remember how we ended up with this ValueErrors instead of a graceful termination with the appropriate status code. @Kai-Striega do you remember?

@andrewdelong just checking - are you suggesting that this is feasible within certain tolerances, or are you mostly concerned about the ValueError?

@mdhaber My concern was mainly that raising a ValueError was not the intention, based on the code comment:

@andrewdelong Yeah, and your example problem was assigned status=2 in presolve, not during simplex. Hopefully @Kai-Striega remembers why it's raising a ValueError instead of properly returning the status in the OptimizeResult. If not, or if there's not a really good reason, that will need to be changed.

I think that I set it as an error before we had merged the two solvers, then didn't want to break backwards compatibility. Looking at it now, I think it will be OK to return the same status and set the message to something similar to comment above.

I think it will be OK to return the same status and set the message to something similar to comment above.

@Kai-Striega You prefer I try to tack that fix/test into this branch, or open a separate issue? (I don't have a test that triggers the scenario that originally motivated the comment.)

You prefer I try to tack that fix/test into this branch, or open a separate issue?

I would prefer this as its own PR. I don't think there's a need for an issue, just open a PR and reference this conversation.

I don't have a test that triggers the scenario that originally motivated the comment.

These occur for the Simplex method on problems 5140, 6139, 7237. There are quite a few combinations of the tests so you will need to find the exact combination yourself (maybe just change the code and see which ones fail)