[Bug] `get_polytope_samples` returns unfeasible samples

[yoavnavon]

🐛 Bug

When generating samples with get_polytope_samples for a constrained problem that is determined in some dimensions (e.g. there is only one solution on that dimension), the function returns samples that don't satisfy the constraints.

To reproduce

>>> from botorch.utils.sampling import get_polytope_samples
>>> import torch

>>> bounds = torch.tensor([[0,0,0],[10,10,10]]).double() # 3-dimensional space

# This doesn't satisfy the constraint
>>> constraint = (torch.tensor([0,1]).double(),torch.tensor([1,1]),20) # dim 0,1 have to be 10 to satisfy the constraint
>>> get_polytope_samples(1,bounds,inequality_constraints=[constraint])
tensor([[8.2054, 2.7486, 4.8757]], dtype=torch.float64)

# This works fine
>>> constraint = (torch.tensor([0,1]).double(),torch.tensor([1,1]),19.9) # There is some freedom
>>> get_polytope_samples(1,bounds,inequality_constraints=[constraint])
tensor([[9.1134, 9.9935, 0.4924]], dtype=torch.float64)

Stack trace/error message

I'm getting this warning the first time I run get_polytope_samples:

.venv/lib/python3.8/site-packages/botorch/utils/sampling.py:471: OptimizeWarning: Solving system with option 'cholesky':True failed. It is normal for this to happen occasionally, especially as the solution is approached. However, if you see this frequently, consider setting option 'cholesky' to False.
  result = scipy.optimize.linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
.venv/lib/python3.8/site-packages/botorch/utils/sampling.py:471: OptimizeWarning: Solving system with option 'sym_pos':True failed. It is normal for this to happen occasionally, especially as the solution is approached. However, if you see this frequently, consider setting option 'sym_pos' to False.
  result = scipy.optimize.linprog(c=c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
.venv/lib/python3.8/site-packages/scipy/optimize/_linprog_ip.py:117: LinAlgWarning: Ill-conditioned matrix (rcond=1.9303e-18): result may not be accurate.
  return sp.linalg.solve(M, r, sym_pos=sym_pos)

Expected Behavior

The constraint is feasible so a point that satisfies the constraint should be returned.

System information

Please complete the following information:

• Botorch: 0.6
• GPyTorch: 1.6.0
• PyTorch:1.10.1
• Computer OS: MacOS 12.0.1

[Balandat]

Hmm interesting, this clearly seems to be due to the constraint being tight. We'll have to look in more detail what exactly is happening, but my hunch is that while the approach of the sampler is to work in the lower-dimension subspace described by A_eq @ x = b_eq, we don't do any preprocessing of the input in oder to detect binding constraints that would reduce the problem dimension further, and so we end up with this issue.

This might be rather hard to handle in full generality, but there might be some things we can do (e.g. by solving some auxiliary LP to detect this).

[Balandat]

cc @jelena-markovic

[yoavnavon]

I think I found part of the root of the issue, in the hit-and-run algorithm here b - A @ x is supposed to be always >= 0, but due to numerical imprecisions the value can be a small negative number (e.g. -1e-12), this obviously happens when the constraints are really tight.

I tried clamping to zero, but for some reason I'm getting always the same value in the dimensions that are not constrained.

[Balandat]

due to numerical imprecisions the value can be a small negative number

Hmm it seems that this should be handled by this piece of code just below?

I tried clamping to zero, but for some reason I'm getting always the same value in the dimensions that are not constrained.

What is your problem specification? Do you have a repro for this? Note that this sampler is meant for bounded polytopes and I don't think it would work for unbounded ones.

[yoavnavon]

What is your problem specification? Do you have a repro for this? Note that this sampler is meant for bounded polytopes and I don't think it would work for unbounded ones.

Sorry I meant to say not constrained except for the bounds, that are obviously also a constraint.

Hmm it seems that this should be handled by this piece of code just below?

I don't think so, looking at the code there:

# given x, the next point in the chain is x+alpha*r
# it also satisfies A(x+alpha*r)<=b which implies A*alpha*r<=b-Ax
# so alpha<=(b-Ax)/ar for ar>0, and alpha>=(b-Ax)/ar for ar<0.
w = (b - A @ x).squeeze() / ar  # b - A @ x is always >= 0
pos = w >= 0
alpha_max = w[pos].min()
# important to include equality here in cases x is at the boundary
# of the polytope
neg = w <= 0
alpha_min = w[neg].max()

The algorithm selects the entries in w that are positive assuming that they are the ones where ar was positive, but if b - A @ x results to have a negative sign, we are actually selecting the entries where ar was negative.

What I think is correct is to do:

w = (torch.clamp(b - A @ x, min=0)).squeeze() / ar 

But this results in that alpha_min = alpha_max = 0, so basically the only option is to stay where we are in the polytope, producing all equal points.

I realize this is probably a limitation of the hit-and-run algorithm for polytope sampling, and probably the way to go is to reduce the space via lp somewhere as you mentioned.

[Balandat]

But this results in that alpha_min = alpha_max = 0, so basically the only option is to stay where we are in the polytope, producing all equal points.

Is this happening in your specific setting? I don't think that's generally true since ar may have both positive and negative entries, in which case it generally won't be that alpha_min = alpha_max = 0.

So I think your fix is correct and needed to deal with numerical imprecisions more generally, but in this specific case where the constraints reduce the effective dimension of the problem we need to do some more work/preprocessing.

[yoavnavon]

Is this happening in your specific setting? I don't think that's generally true since ar may have both positive and negative entries, in which case it generally won't be that alpha_min = alpha_max = 0.

Oh yeah you are right, actually it also depends on the random sampling of the r values as well of my particular setting.

[saitcakmak]

@Balandat was this resolved by #1052?

[Balandat]

Not quite, I think - This improved some cases in which this can happen, but if the (non-equality) constraints are such that they indeed result in a zero-volume polytope then I believe we need to do something smarter here. This doesn't seem like a quick fix though. Let's leave this open for now to track this.