Sunday, October 18, 2020

PuLP mystery

PuLP is a popular Python-based modeling tool for LP and MIP models. In [1], a user asked a (somewhat trivial) question about PuLP syntax. But there was an interesting wrinkle in the model that was unfamiliar to me. The basic issue is that sometimes PuLP accepts a NumPy array in expressions:
import pulp as p
import numpy as np



prob = p.LpProblem("test",p.LpMaximize)
prob += x+(a-y)

This fragment creates an objective. But we have a strange term here. a-y is funny because a is a NumPy array (vector) and y is a scalar PuLP variable. Usually, adding a scalar to a vector is interpreted as elementwise addition: \[\begin{pmatrix} a_0 - y \\ a_1 - y \\ a_2-y \end{pmatrix}\] How one would interpret an objective like: \[ x + \begin{pmatrix} a_0 - y \\ a_1 - y \\ a_2-y \end{pmatrix}\] is not clear to me. I would say this is an error. However, PuLP accepts this, and prints the model as:

1*x + -3*y + 6
x free Continuous
y free Continuous

So, apparently the interpretation is: \[ x + \sum_i (a_i-y)=x + \sum_i a_i - n\cdot y\] where \(n\) is the length of vector \(a\). But then again, we would expect PuLP to accept
  prob += (a-y)

as objective. However, this produces the error message: 

TypeError                                 Traceback (most recent call last)
<ipython-input-6-bf3a3a41e9da> in <module>()
      9 prob = p.LpProblem("test",p.LpMaximize)
---> 10 prob += (a-y)
     11 prob

/usr/local/lib/python3.6/dist-packages/pulp/ in __iadd__(self, other)
   1528    = name
   1529         else:
-> 1530             raise TypeError("Can only add LpConstraintVar, LpConstraint, LpAffineExpression or True objects")
   1531         return self

TypeError: Can only add LpConstraintVar, LpConstraint, LpAffineExpression or True objects

I have no idea what is going on here. It could be all just a bug, but even that can not easily explain the behavior of sometimes accepting (a-y) and sometimes refusing it.

The underlying problem can actually be traced back to operator overloading as shown by [2]. Who owns the + or -? The following experiment demonstrates the issue in stark terms:

>>> a+x
array([1*x + 1, 1*x + 2, 1*x + 3], dtype=object)
>>> x+a
1*x + 6

In the first case, NumPy handles the +, leading to the interpretation: \[\begin{pmatrix} a_0  \\ a_1  \\ a_2 \end{pmatrix} + x = \begin{pmatrix} x+a_0  \\ x+a_1  \\x+a_2  \end{pmatrix}\] In the second case, the + is dealt with by PuLP. This gives \[x+\begin{pmatrix} a_0  \\ a_1  \\ a_2 \end{pmatrix}=x+\sum_i a_i\]  


This is mind-bending (a.k.a. crazy). It would be better if things were a bit more predictable ("orthogonal" is the term used in programming language design). Preferably, we don't want to spend time on figuring out how to interpret a + or a -

Advise: it is better not to mix PuLP with NumPy.



  1. When you add an Iterable to an LpVariable, it adds each element of that Iterable separately. I wrote an explanation here: