Lower bound greater than upper bound

[johnnybonney]

library(ivmte)

# error doesn't reproduce if I use all years
dt <- AE[AE$yob < 55,]

args <- list(data = dt,
             outcome = "worked",
             propensity = morekids ~ factor(samesex)*factor(yob),
             m0 = ~factor(yob)*u + factor(yob):I(u^2),
             m1 = ~factor(yob)*u + factor(yob):I(u^2),
             target = 'late',
             late.to = list(samesex = 1),
             late.from = list(samesex = 0),
             solver = "gurobi",
             direct = "lp",
             initgrid.nu = 100,
             audit.nu = 200)  
est <- do.call(ivmte, args)
print(est)

Bounds on the target parameter: [-0.08306298, -0.0857325]
Audit terminated successfully after 2 rounds

Just to verify that I'm not going crazy...

> est$bounds[["upper"]] - est$bounds[["lower"]]
[1] -0.00266952

[a-torgovitsky]

Minimum criterion is not zero, right?

Is it sensitive to the audit grid?

[johnnybonney]

The minimum criterion is 2.681385e-06, so not zero.

Is it sensitive to the audit grid?

I tried a few different u-grid values. For initgrid.nu values less than 200, I get these flipped bounds. For initgrid.nu values of 200 or greater, I get

Error:  The minimization problem was proven to be infeasible or unbounded. Since a minimum criterion was found, the model is most likely feasible but unbounded. This can happen if the initial grid is too small. Try increasing the parameters 'initgrid.nx' and 'initgrid.nu'. If the model is indeed infeasible, consider exporting the model and passing it to the solver for more details.

(As an aside, this isn't the first time I've run into cases where increasing the grid size leads the problem to go from bounded to unbounded...)

[a-torgovitsky]

Do you mean initgrid.nu or audit.nu? It's the audit grid that should matter.

Also I assume that the default value of audit.nx is large enough here that changing it doesn't affect anything?

[johnnybonney]

Also I assume that the default value of audit.nx is large enough here that changing it doesn't affect anything?

Yes, that's right. There are only 11 covariate values here, and the default audit.nx is 2500. (And I've verified that playing around with audit.nx and initgrid.nx doesn't change anything, as we would expect.)

Do you mean initgrid.nu or audit.nu?

I mean initgrid.nu. Keeping audit.nu = 200 fixed, here is what I get for a few different values of initgrid.nu:

• initgrid.nu = 20: bounds are [-0.08487877, -0.0857325], criterion is 2.681385e-06
• initgrid.nu = 50: minimization problem is infeasible/unbounded, criterion is 2.681385e-06
• initgrid.nu = 100: bounds are [-0.08306298, -0.0857325], criterion is 2.681385e-06
• initgrid.nu = 150: bounds are [-0.08400145, -0.0857325], criterion is 2.681385e-06
• initgrid.nu = 200: minimization problem is infeasible/unbounded, criterion is 2.681385e-06

So, something is going on with the lower bound. The upper bound is stable, as we would hope both bounds would be (given a fixed audit grid). Changing audit.nu does not seem to change any of the above results.

[a-torgovitsky]

I would think that making audit.nu finer is going to make the lower/upper bounds more stable... Anyway, the results you get already look pretty stable when they exist.

I'm more concerned with what's going on with initgrid.nu = 50 and initgrid.nu = 200.

We know the 2nd step problems will be feasible, at least up to any numerical issues. So if they are unbounded that means that the audit grid is not sufficiently constraining the MTRs to prevent getting an unbounded target parameter. But your MTR specification is not that flexible, and target parameter not that exotic, so that doesn't really make sense.

@johnnybonney can you start by determining whether these problems are actually infeasible or unbounded?

Actually, maybe I'm wrong about this not being exotic. Your MTR specification is a quadratic for each covariate value. And your target parameter involves covariate-specific LATEs. So all you need is for one LATE to be over a region that is not constrained by the audit grid. Which could happen if the x-specific first stage is quite small for one covariate value?

[johnnybonney]

So all you need is for one LATE to be over a region that is not constrained by the audit grid. Which could happen if the x-specific first stage is quite small for one covariate value?

This is something that I have tried to be mindful about. The smallest first stage here is 0.034. Given the selection of u grid points via the Halton sequence, 200 audit points should be more than enough to constrain the problem. (In fact, 50 points is enough to constrain any first stages larger than 0.031.) To wit:

library(data.table)
setDT(dt)
fs <- dt[, .(mean(morekids[samesex == 1]) - mean(morekids[samesex == 0])),
         by = yob]

ugrid <- sort(ivmte::rhalton(50))
max_udiff <- max(ugrid[2:length(ugrid)] - ugrid[1:length(ugrid)-1])

message(paste("Minimum first stage is", min(fs$V1)))
message(paste("Initial u-grid fine enough to cover differences exceeding", max_udiff))

Minimum first stage is 0.0341548985110104
Initial u-grid fine enough to cover differences exceeding 0.03125

can you start by determining whether these problems are actually infeasible or unbounded?

What would be the right way to go about this?

My first thought is, the problem was feasible and bounded with an initial grid of 20 u-points and an audit grid of 200 u-points. (This problem required an audit iteration.) The problem with 50 initial u-points (which include those 20 initial points) must then be bounded, right?

And it is feasible and bounded with an initial grid of 100 points (solved with an audit iteration). So it must be feasible with an initial grid of 50 points, right? Or am I missing something?

[a-torgovitsky]

What would be the right way to go about this?

I think just run it with the Gurobi output on and see what Gurobi says.

My first thought is, the problem was feasible and bounded with an initial grid of 20 u-points and an audit grid of 200 u-points. (This problem required an audit iteration.) The problem with 50 initial u-points (which include those 20 initial points) must then be bounded, right? And it is feasible and bounded with an initial grid of 100 points (solved with an audit iteration). So it must be feasible with an initial grid of 50 points, right? Or am I missing something?

No this is the right logic...I am thinking this is a bug.

Here's another test you can do that would help Instead of the target parameter being the unconditional LATE, make it just the conditional LATE for a single x covariate. Then loop through all x covariate values and compute bounds. There should be at least one value for which you get something unbounded.

[johnnybonney]

OK, so if I run it with debug = T, the Gurobi output is

Lower bound optimization statistics:
------------------------------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 2355 rows, 198 columns and 11598 nonzeros
Model fingerprint: 0x8e66f7a4
Coefficient statistics:
  Matrix range     [1e-06, 1e+00]
  Objective range  [1e-03, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e-06, 1e+00]
Presolve removed 1146 rows and 0 columns
Presolve time: 0.00s
Presolved: 1209 rows, 1340 columns, 7290 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s

Solved in 311 iterations and 0.01 seconds
Infeasible or unbounded model

ivmte doesn't return any further details about the problem. I did find that if I run it with solver.options = list(presolve = 0), we get bounds: [-0.08487877, -0.08573251]. The lower bound Gurobi output for the first iteration of that audit is

Lower bound optimization statistics:
------------------------------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (mac64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 2355 rows, 198 columns and 11598 nonzeros
Model fingerprint: 0xf033d455
Coefficient statistics:
  Matrix range     [1e-06, 1e+00]
  Objective range  [1e-03, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e-06, 1e+00]
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s
     511   -8.4878781e-02   0.000000e+00   0.000000e+00      0s

Solved in 511 iterations and 0.02 seconds
Optimal objective -8.487878089e-02

At first blush, and as a non-expert of Gurobi, it would seem to me that presolving may have something to do with it. But what I find below suggests the problem arises somewhere else.

---

Instead of the target parameter being the unconditional LATE, make it just the conditional LATE for a single x covariate. Then loop through all x covariate values and compute bounds.

If I set late.X = list(yob = 54), then I get an an unbounded problem. All other values of yob produce bounds. yob = 54 is the least common yob in the subsetted data, at 2.5% of the sample, but it still has over 5000 observations. This also happens to be the group with the smallest first stage, of 0.034, that I mentioned earlier...

Also, presolve = 0 doesn't save us here---this problem remains unbounded.

---

During the above exercise, I wanted to try comparing three different estimates for a given yob value: (1) the bounds I get using the full yob < 55 subsample, using late.X = list(yob = X); (2) the bounds I get using only the subsample of data that is limited to yob == X; and (3) the linear IV estimate for the subsample.

For yob = 54 (the problematic yob):

• full data: unbounded (as mentioned above), with min. criterion of 2.681385e-06
• subset to just yob == 54: bounds of [0.5297372, 0.5297985], with a criterion of 0.0001145109
• IV estimate: 0.951389

For yob = 53 (as an example):

• full data: [-0.2190124, -0.2614573] (criterion of 2.681385e-06)
• subset to yob == 53: [-0.2565897, -0.2565897] (criterion of 0)
• IV estimate: -0.2565897

One last piece of evidence: if I re-run target = "late" with the yob < 54 subset (that is, removing the problematic yob), I get [-0.1035615, -0.1035615] and a minimum criterion of 0. This exactly reproduces the average LATE that I get by estimating things manually.

---

Taken together, this suggests that there are two different things going on here.

1. There is something weird with yob = 54 that is leading to problems being unbounded when they should be bounded, and that leads to "point estimates" for LATE(yob = 54) that differ from their true point estimates.
2. Some nonseparability is introduced between the different covariate values. The MTRs are specified such that they should be fully interacted with yob dummies, but problems with yob = 54 are reverberating through to cause problems with, e.g., LATE(yob = 53). These problems include lower bounds that exceed upper bounds (the initial issue). My understanding is that this "nonseparability" is from the criterion.

[a-torgovitsky]

Ok good that's progress. Let's focus on this:

For yob = 54 (the problematic yob):

full data: unbounded (as mentioned above), with min. criterion of 2.681385e-06
subset to just yob == 54: bounds of [0.5297372, 0.5297985], with a criterion of 0.0001145109
IV estimate: 0.951389

If $E[Y(1) - Y(0) \vert G = \text{cp}] = .95$, then my guess is that either $E[Y(1) \vert G = \text{cp}] > 1$ or $E[Y(0) \vert G = \text{cp}] < 0$. What do you get when you estimate those? If my guess is right, that means we are rejecting the LATE model here (basically a Kitagawa 2015 test). However that really shouldn't be a problem, as the method in ivmte is capable of handling all sorts of misspecification.

---

So really I think the issue here is with presolve. @jkcshea do you know anything about that? Have we played around with presolve before?

[johnnybonney]

Estimating the complier levels:

dt <- AE[AE$yob == 54,]
ey_gz1 <- mean(dt[dt$samesex == 1,]$worked)
ey_gz0 <- mean(dt[dt$samesex == 0,]$worked)
ey_at <- mean(dt[dt$morekids == 1 & dt$samesex == 0,]$worked)
ey_nt <- mean(dt[dt$morekids == 0 & dt$samesex == 1,]$worked)
p1 <- mean(dt[dt$samesex == 1,]$morekids)
p0 <- mean(dt[dt$samesex == 0,]$morekids)

ey1_cp <- (ey_gz1 - ey_nt*(1 - p1) - ey_at*p0) / (p1 - p0)
ey0_cp <- (ey_gz0 - ey_nt*(1 - p1) - ey_at*p0) / (p1 - p0)

message(paste("E[Y(1)|G=cp] =", round(ey1_cp, 4)))
message(paste("E[Y(0)|G=cp] =", round(ey0_cp, 4)))

E[Y(1)|G=cp] = 0.6736
E[Y(0)|G=cp] = -0.2778

Sure enough, $E[Y(0)|G = cp] < 0$.

Also, as I mentioned, even setting presolve = 0 doesn't fix the problem when we estimate late.X = list(yob = 54). That problem is unbounded with or without presolving.

[a-torgovitsky]

What happens if you specify the target parameter manually? I worry that late.X = list(yob = 54) may not be doing what we think it is...

[johnnybonney]

If I specify the target parameter (the conditional on yob = 54 LATE) manually, then I still get an unbounded problem. This is true whether or not presolve = 0.

If I switch to the yob = 53 LATE, I get an upper bound that differs from what I got using late.X = list(yob = 53):

dt <- AE[AE$yob < 55,]

# custom weights, knots to reflect LATE | yob = 53
p1_x53 <- mean(dt[dt$samesex == 1 & dt$yob == 53,]$morekids)
p0_x53 <- mean(dt[dt$samesex == 0 & dt$yob == 53,]$morekids)
Pr_x53 <- mean(dt$yob == 53)

late_wt_m1 <- function(yob) as.integer(yob == 53) / ((p1_x53 - p0_x53) * Pr_x53)
late_wt_m0 <- function(yob) -1 * late_wt_m1(yob)

args <- list(data = dt,
             outcome = "worked",
             target.weight0 = c(0, late_wt_m0, 0),
             target.weight1 = c(0, late_wt_m1, 0),
             target.knots0 = c(p0_x53, p1_x53),
             target.knots1 = c(p0_x53, p1_x53),
             propensity = morekids ~ factor(samesex) * factor(yob),
             m1 = ~factor(yob)*u + factor(yob):I(u^2),
             m0 = ~factor(yob)*u + factor(yob):I(u^2),
             solver = "gurobi",
             direct = "lp",
             initgrid.nu = 50,
             audit.nu = 200)  

do.call(ivmte, args)

Bounds on the target parameter: [-0.2190124, -0.3006759]
Audit terminated successfully after 2 rounds 

Recall that the late.X = list(yob = 53) bounds were [-0.2190124, -0.2614573]. So my manual "lower bound" is the same, but my upper bound is different here.

[johnnybonney]

To confirm that $E[Y(0) | G = cp] < 0$ is causing the problem---I find that if I remove the upper and lower bounds on m1 and m0, then all these problems vanish. That is, I get bounds, I get point identification, and point estimates are not sensitive to the way I have specified m1 or m0 (issue #224).

edit: and the estimates aren't sensitive to whether I specify the target parameter manually or using late.X.

[a-torgovitsky]

Yeah if you relax the bounds on m0 and m1 then you are effectively making the model correctly specified and allowing the MTRs to fit the moments exactly. That's not something we can do in general, but good to know that pinpointed the issue.

---

On this:

If I switch to the yob = 53 LATE, I get an upper bound that differs from what I got using late.X = list(yob = 53):

Can you add that to #220 ? I suspect it's the same problem just conditional on covariates, and maybe @jkcshea has fixed it for conditional LATEs too (it hasn't been pushed yet), but if not making a post there will be a reminder.

---

I'm having trouble reproducing this:

If I specify the target parameter (the conditional on yob = 54 LATE) manually, then I still get an unbounded problem. This is true whether or not presolve = 0.

Here's my code, maybe I'm missing something:

library(ivmte)

dt <- AE[AE$yob == 54,]
args <- list(data = dt,
             outcome = "worked",
             propensity = morekids ~ factor(samesex),
             m0 = ~ u + I(u^2),
             m1 = ~ u + I(u^2),
             target = 'late',
             late.to = list(samesex = 1),
             late.from = list(samesex = 0),
             solver = "gurobi",
             direct = "lp",
             initgrid.nu = 100,
             audit.nu = 200)
est <- do.call(ivmte, args)
print(est)

for me gives

Bounds on the target parameter: [0.5297372, 0.5297985]
Audit terminated successfully after 1 round 

[johnnybonney]

I was unclear when I said this:

If I specify the target parameter (the conditional on yob = 54 LATE) manually, then I still get an unbounded problem. This is true whether or not presolve = 0.

What I meant was, when I use the full sample but use knots/weights to manually specify LATE(yob = 54). That is, take my example from my previous comment and just replace 54 with 53 everywhere. Then the problem is unbounded.

This is what I meant:

dt <- AE[AE$yob < 55,]

# custom weights, knots to reflect LATE | yob = 54
p1_x54 <- mean(dt[dt$samesex == 1 & dt$yob == 54,]$morekids)
p0_x54 <- mean(dt[dt$samesex == 0 & dt$yob == 54,]$morekids)
Pr_x54 <- mean(dt$yob == 54)

late_wt_m1 <- function(yob) as.integer(yob == 54) / ((p1_x54 - p0_x54) * Pr_x54)
late_wt_m0 <- function(yob) -1 * late_wt_m1(yob)

args <- list(data = dt,
             outcome = "worked",
             target.weight0 = c(0, late_wt_m0, 0),
             target.weight1 = c(0, late_wt_m1, 0),
             target.knots0 = c(p0_x54, p1_x54),
             target.knots1 = c(p0_x54, p1_x54),
             propensity = morekids ~ factor(samesex) * factor(yob),
             m1 = ~factor(yob)*u + factor(yob):I(u^2),
             m0 = ~factor(yob)*u + factor(yob):I(u^2),
             solver = "gurobi",
             direct = "lp",
             initgrid.nu = 50,
             audit.nu = 200)  

est <- do.call(ivmte, args)

When I subset the data itself on yob = 54, as in your example, then I get bounds just like you.

[johnnybonney]

Can you add that to 220 ? I suspect it's the same problem just conditional on covariates

I posted it, but as soon as I did, I realized it probably isn't the same problem, since turning presolve off leads to the same bounds for package and manual estimation.

[a-torgovitsky]

Ok I see, is there a way to strip this down and get a more minimal example though? There's a lot going on here and it's a bit difficult to isolate.

I tried increasing criterion.tol and still get "infeasible or unbounded." So I'm guessing it's an unbounded problem. That doesn't make sense with the initial grid and MTR specification as it is... But it's sufficiently complicated that it's hard to see exactly what's going on.

[johnnybonney]

Here is a more minimal example.

library(ivmte)

dt <- AE[AE$yob %in% 53:54,]
dt$X <- as.integer(dt$yob == 54)

args <- list(data = dt,
             outcome = "worked",
             propensity = morekids ~ samesex + X + X:samesex,
             m1 = ~u + I(u^2) + X + X:u + X:I(u^2),
             m0 = ~u + I(u^2) + X + X:u + X:I(u^2),
             target = 'late',
             late.to = list(samesex = 1),
             late.from = list(samesex = 0),
             solver = "gurobi",
             direct = "lp")  

# increasing initial grid
est1 <- do.call(ivmte, c(args, list(initgrid.nu = 10, audit.nu = 100)))
est2 <- do.call(ivmte, c(args, list(initgrid.nu = 12, audit.nu = 100)))
est3 <- do.call(ivmte, c(args, list(initgrid.nu = 15, audit.nu = 100)))

# increasing criterion.tol
est4 <- do.call(ivmte, c(args, list(initgrid.nu = 12, audit.nu = 100, criterion.tol = 0.05)))

For initgrid.nu = 10, we have

> print(est1)

Bounds on the target parameter: [0.1005246, 0.100737]
Audit terminated successfully after 2 rounds 

For initgrid.nu = 12, we get the error message

Error:  The minimization problem was proven to be infeasible or unbounded. Since a minimum criterion was found, the model is most likely feasible but unbounded. This can happen if the initial grid is too small. Try increasing the parameters 'initgrid.nx' and 'initgrid.nu'. If the model is indeed infeasible, consider exporting the model and passing it to the solver for more details.

For initgrid.nu = 15, we have

> print(est3)

Bounds on the target parameter: [0.1005246, 0.100737]
Audit terminated successfully after 2 rounds 

If we increase criterion.tol high enough, we can bound the initgrid.nu = 12 problem:

> print(est4)

Bounds on the target parameter: [0.09586238, 0.1045862]
Audit terminated successfully after 2 rounds 

This suggests that the problem is infeasible, though it should be feasible.

Edit: turning presolve off also allows us to get bounds for the initgrid.nu = 12 problem, without changing criterion.tol:

do.call(ivmte, c(args, list(initgrid.nu = 12, audit.nu = 100,
                            solver.options = list(presolve = 0))))

Bounds on the target parameter: [0.1005337, 0.100737]
Audit terminated successfully after 2 rounds 

[a-torgovitsky]

Excellent, thanks. I will play around with this on Monday.

[jkcshea]

I think these errors may have something to do with the feasibility tolerance. Gurobi's default feasibility tolerance is 1e-6. The bounds are no longer flipped in original example once I pass the option solver.options = list(FeasibilityTol = 1e-07). All the errors in the most recent examples also disappear once I reduce the tolerance.

In all these examples, the minimum criteria are very small. e.g. @johnnybonney's original example had a minimum criterion of 2.681385e-06, very close to Gurobi's feasibility tolerance. So the estimated bounds could be sensitive to violations of the criterion constraints within the tolerance level.

(This is related to my earlier comment here about Gurobi's feasibility tolerance. But maybe my understanding of the math is wrong/incomplete.)

[a-torgovitsky]

An interesting twist here is that the upper bound problem appears to solve correctly. But that should have exactly the same constraints as the lower bound problem...

library(ivmte)

dt <- AE[AE$yob %in% 53:54,]
dt$X <- as.integer(dt$yob == 54)

args <- list(data = dt,
             outcome = "worked",
             propensity = morekids ~ samesex + X + X:samesex,
             m1 = ~u + I(u^2) + X + X:u + X:I(u^2),
             m0 = ~u + I(u^2) + X + X:u + X:I(u^2),
             target = 'late',
             late.to = list(samesex = 1),
             late.from = list(samesex = 0),
             solver = "gurobi",
             noisy = T,
             debug = T,
             direct = "lp")
est2 <- do.call(ivmte, c(args, list(initgrid.nu = 12, audit.nu = 100)))

and the output looks like this:

Obtaining propensity scores...

Generating target moments...
    Integrating terms for control group...
    Integrating terms for treated group...

Performing direct MTR regression...
    MTR is not point identified.

Performing audit procedure...
    Solver: Gurobi ('gurobi')
    Generating initial constraint grid...

    Audit count: 1

Minimum criterion optimization statistics:
------------------------------------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 124 rows, 36 columns and 576 nonzeros
Model fingerprint: 0x00aef5dd
Coefficient statistics:
  Matrix range     [2e-05, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e-04, 1e+00]
Presolve removed 58 rows and 0 columns
Presolve time: 0.00s
Presolved: 66 rows, 90 columns, 388 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s
      34    4.3432674e-05   0.000000e+00   0.000000e+00      0s

Solved in 34 iterations and 0.00 seconds
Optimal objective  4.343267353e-05

    Minimum criterion: 4.343267e-05
    Obtaining bounds...

Lower bound optimization statistics:
------------------------------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 125 rows, 36 columns and 600 nonzeros
Model fingerprint: 0x9f77e0a0
Coefficient statistics:
  Matrix range     [2e-05, 1e+00]
  Objective range  [2e-02, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e-05, 1e+00]
Presolve removed 58 rows and 0 columns
Presolve time: 0.00s
Presolved: 67 rows, 90 columns, 412 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s

Solved in 39 iterations and 0.00 seconds
Infeasible or unbounded model

Upper bound optimization statistics:
------------------------------------
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 125 rows, 36 columns and 600 nonzeros
Model fingerprint: 0x2ff30c46
Coefficient statistics:
  Matrix range     [2e-05, 1e+00]
  Objective range  [2e-02, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e-05, 1e+00]
Presolve removed 58 rows and 0 columns
Presolve time: 0.00s
Presolved: 67 rows, 90 columns, 412 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s
      43    1.0076018e-01   0.000000e+00   0.000000e+00      0s

Solved in 43 iterations and 0.00 seconds
Optimal objective  1.007601756e-01

    Model was infeasible or unbounded.
Error:  The minimization problem was proven to be infeasible or unbounded. Since a minimum criterion was found, the model is most likely feasible but unbounded. This can happen if the initial grid is too small. Try increasing the parameters 'initgrid.nx' and 'initgrid.nu'. If the model is indeed infeasible, consider exporting the model and passing it to the solver for more details.

@jkcshea could you try flipping the order of the solves (max first, then min) and report what you get for this example? Another thing to try would be: min, max, then min again. (I wonder if Gurobi is warm-starting these...)

[jkcshea]

I tried both patterns of minimizing/maximizing and got the same results. So it doesn't look like Gurobi is warm-starting... I'll keep playing with this to see if I can figure out what's going on.

[a-torgovitsky]

Odd, very odd. Is it possible to reproduce this by writing out the .mps and trying to solve directly with Gurobi. If all that is different is the max vs. min, I would say that's an issue for the Gurobi support board.

Another sanity check would be to change $min\quad c'x$ to $-1 * max\quad-c'x$ Surely that should give the same thing...

[jkcshea]

Is it possible to reproduce this by writing out the .mps and trying to solve directly with Gurobi.

Yep, it is possible to reproduce this directly with Gurobi.

Another sanity check would be to change... Surely that should give the same thing...

Yep, it does. Gurobi returns the same upper bound, but is unable to give us a lower bound.

If all that is different is the max vs. min, I would say that's an issue for the Gurobi support board.

That is indeed the only difference. In the package, I just switch model$modelsense from min to max to get both bounds. The issue has been posted on Gurobi's support portal, here.

[a-torgovitsky]

Ok let's see what they say. Of course this behavior suggests that the model is in fact unbounded not infeasible. But @johnnybonney 's example of raising criterion.tol suggests the opposite. Hmm...

[jkcshea]

Gurobi was quick to respond! The link above has the full response, but here is what the support staff said:

I think the issue here is largely due to the numeric instabilities. If you tighten the tolerances by one more order of magnitude as compared to the defaults, namely

intFeasTol = 1e-7 feasibilityTol = 1e-8

then both runs show as infeasible on my end. The differences you have observed are thus likely to the solver taking a different algorithmic pass depending on solving a maximization or minimization version of the model, and come down to a different round of errors etc.

They also suggested I look into DualReductions parameter of Gurobi. It is related to presolve, and turning it off allows Gurobi to determine whether the model is unbounded or infeasible. When I turn off dual reductions in the minimization problem, Gurobi reports the problem is infeasible.

[a-torgovitsky]

Can you repost the link? The one above doesn't seem to work on my end

[jkcshea]

Hm, the link I would've posted would be identical to the one above. I think you cannot access the link because I posted the question using their support portal for academic licensees. So my request isn't public. That may also explain why Gurobi responded almost immediately. I'll forward you the email with their response instead.

[a-torgovitsky]

Got the email, thanks.

I'm not really sure what the right solution is here.

One possibility is that this is related to #218, and introducing additional stabilizing fixes will deal with this problem too.

But barring that, it seems like the solution is to just twiddle with the Gurobi knobs a bit. Which is unfortunate since it's going to be hard for an average user to do that. One thing we could potentially do is loop through some common settings (e.g. presolve) until we get a clean solve before returning an error. It's not very elegant though.