Add a matrix of constraints in a LP

[ncohen]

Some users reported two problems with Sage's interface with LP solvers.

1) It is slow to add constraints one by one to the solver's data structure

2) It should be possible to add a matrix of constraints all at once

With a new syntax like that (thanks Dima), both problems can be circumvented:

p = MixedIntegerLinearProgram()
x = p.new_variable()
p.add_constraint(A_matrix*x <= a_vector)

This would consider x as the vector x[0],x[1],...,x[size_of_the_matrix - 1].

In the add_constraint function, we can then call the backend's function to add 10 000 constraints at the same time (probably better for the data structure), and then fill the entries.

Nathann

[dimpase]

This would also be very natural syntax to specify polyhedra.

[jdemeyer]

Looks like a duplicate of #7293... (which has no code and no comments)

[vbraun]

One technical problem is that matrices in Sage can always be multiplied, so you can't create matrices over linear functions.

[ncohen]

One technical problem is that matrices in Sage can always be multiplied, so you can't create matrices over linear functions.

I don't understand O_o

We do not mind if two matrices are multiplied, we mind when two LP variables are multiplied together. We just have to scream when two variables are multiplied, don't we ? So if a guy multiplies two matrices containing variables it will try to multiply the variables at some point so we are safe, aren't we ?

Nathann

[dimpase]

Replying to vbraun:

One technical problem is that matrices in Sage can always be multiplied, so you can't create matrices over linear functions.

but why would you want this? All that is needed is to pull the matrix apart, row by row, and create a linear function from each of them...

[dimpase]

what might help, conceptually, if one can do linear maps just as well as linear functions. Then, when, say, A and B are matrices, one could imagine to do

x = p.new_variable()
y = p.new_variable()
p.add_constraint(A*x+B*y <= c)

so that A*x+B*y is a linear map from the space of dimension A.ncol()+B.ncol() to the space of dimension len(c).

[vbraun]

I've started with tensors (over the base ring) of linear functions and free modules. This implements:

sage: m = matrix([[1,2],[3,4]]);  m
[1 2]
[3 4]
sage: p = MixedIntegerLinearProgram()
sage: v = p.new_variable(nonnegative=True)
sage: m.row(0) * v[0]              # vector * linear function
(1.0, 2.0)*x_0
sage: m * v[0]                     # matrix * linear function
[x_0   2*x_0]
[3*x_0 4*x_0]
sage: v * m                        # MIPVariable * matrix
(3.0, 4.0)*x_1 + (1.0, 2.0)*x_0

The last multiplication including MIPvariables is incomplete: You can't write m * v since MIPvariables do not fit into the parent/element framework, so they don't really participate in coercion. This needs to be fixed, probably best by making MIPvariables the elements of the parent MixedIntegerLinearProgram. Once that is done it'll be easy to add (in)equalities of the new tensor elements.

---

New commits:

​846cba6

Implement tensor product of linear functions and free modules

[ncohen]

Hellooooo !

This needs to be fixed, probably best by making MIPvariables the elements of the parent MixedIntegerLinearProgram.

Can it result in any slowdown for the usual usage of those variables ?

Nathann

[dimpase]

Replying to ncohen:

Hellooooo !

This needs to be fixed, probably best by making MIPvariables the elements of the parent MixedIntegerLinearProgram.

Can it result in any slowdown for the usual usage of those variables ?

I'd be surprised. The real bottleneck is in interacting with the solver and solving itself, anyway.

Besides, the current prevalent use of MIP frontend is very inefficient, as you collect the linear system into a matrix, internally, and the backend is called each time a new constraint is added, right? It's certainly quicker to pass the matrix directly and call the backend once.

[ncohen]

Replying to dimpase:

Replying to ncohen:

Hellooooo !

This needs to be fixed, probably best by making MIPvariables the elements of the parent MixedIntegerLinearProgram.

Can it result in any slowdown for the usual usage of those variables ?

I'd be surprised. The real bottleneck is in interacting with the solver and solving itself, anyway.

Well, I need to compute a max matching on a large bipartite graph for some designs I build these days, and ...

sage: g=graphs.CompleteBipartiteGraph(100,200)
sage: %time _=g.matching(algorithm="LP", verbose=2)
...

Root node processing (before b&c):
  Real time             =    0.17 sec. (104.83 ticks)
Parallel b&c, 4 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.17 sec. (104.83 ticks)
CPU times: user 1.25 s, sys: 24 ms, total: 1.28 s
Wall time: 1.06 s

Sage reports more than 1second and the solver report 0.17. I don't really know where the computations go to be honest at %prun does not say much. And I suspect that it does not say much because the code that takes time is written in Cython, i.e. the LP variables. But I haven't found a proof yet.

Besides, the current prevalent use of MIP frontend is very inefficient, as you collect the linear system into a matrix, internally, and the backend is called each time a new constraint is added, right? It's certainly quicker to pass the matrix directly and call the backend once.

It may be faster, but the two must be possible. You don't want to have to build a matrix for any of these problems, the code would be ugly as hell. Being able to index variables with "anything" certainly adds something, and not only for readability.

​http://steinertriples.fr/ncohen/tut/LP_examples/

Nathann

[dimpase]

Replying to vbraun:

I've started with tensors (over the base ring) of linear functions and free modules. This implements:

...
sage: v * m                        # MIPVariable * matrix
(3.0, 4.0)*x_1 + (1.0, 2.0)*x_0

this is very nice; I hope this can be a basis for implementing an interface to "conic programming". That is, whenever you have a convex cone K in R^n (e.g. the positive ortant, or {(x_0,...,x_{n-1}) | x_0^2 >= sum_{j=1}^{n-1} x_j^2}, known as "icecream cone") you have a partial order <_K on R^n defined by x<_K y iff y-k is in K.

then one can do linear optimisation on the intersection of K with an affine subspace. For K being the positive ortant this is just the usual LP; for K the icecream cone this is a kind of norm minimisation, etc. CVXOPT has implementations for several different cones like this (by the way, semidefinite programming is yet another example, K being the cone of psd matrices in R^nxn).

[dimpase]

Replying to ncohen:

Well, I need to compute a max matching on a large bipartite graph for some designs I build these days, and ...

BTW, using LP for this is a very inefficient thing. The classical combinatorial algorithms (see e.g. ​Hungarian method) will be much faster...

[ncohen]

BTW, using LP for this is a very inefficient thing.

How do you think you can prove this assertion ? Especially when the polytope of perfect matchings in a bipartite graph is integer ?..

The classical combinatorial algorithms (see e.g. ​Hungarian method) will be much faster...

Claim without proof. Also, this isn't implemented, lest of all efficiently implemented.

Nathann

[dimpase]

Replying to ncohen:

BTW, using LP for this is a very inefficient thing.

How do you think you can prove this assertion ? Especially when the polytope of perfect matchings in a bipartite graph is integer ?..

The classical combinatorial algorithms (see e.g. ​Hungarian method) will be much faster...

Claim without proof. Also, this isn't implemented, lest of all efficiently implemented.

Come on, it is implemented e.g. here: ​networkx

Convert your graph into networkx one and call it... Although indeed this might not be the fastest implementation around, sure.

[ncohen]

Come on, it is implemented e.g. here: ​networkx

No, that is Edmond's algorith for general graphs (it is not bipartite-specific) and it is also the default algorithm called by Graph.matching in case. Plus it's pure Python.

Nathann

[dimpase]

Replying to ncohen:

Come on, it is implemented e.g. here: ​networkx

No, that is Edmond's algorith for general graphs (it is not bipartite-specific) and it is also the default algorithm called by Graph.matching in case. Plus it's pure Python.

Edmonds' on bipartite graphs will do more or less the same amount of work as the Hungarian on bipartite graphs. Trust me, I taught these things few times, I used to know details, and even implemented Hungarian method in C and in Python :-)

Although Hungarian isn't the fastest guy in town.

PS. It is common (to the combinatorial optimisation community) knowledge that LP-based methods hopelessly lose here. I don't make baseless claims here, I know well what I am talking about...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​04f9ce1

Add a parent for MIPVariables

[vbraun]

I've added a Parent for MIPVariables, this allows for matrix multiplication to work from both sides:

            sage: p = MixedIntegerLinearProgram()
            sage: v = p.new_variable()
            sage: m = matrix([[1,2], [3,4]])
            sage: v * m
            (3.0, 4.0)*x_1 + (1.0, 2.0)*x_0
            sage: m * v
            (2.0, 4.0)*x_1 + (1.0, 3.0)*x_0

The parent must refer back to the mip so there is a reference cycle (#12616). Deal with it...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​4c64eb7	Make the parent of MIPVariables a singleton
​0a09269	Enable p.<x> = MixedIntegerLinearProgram() syntax

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c09be57

Match x_i print order

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a94e77e

Symbolic (in)equality involving linear tensors

[vbraun]

Inequalities with vector/matrix-valued linear functions:

sage: mip.<x> = MixedIntegerLinearProgram()
sage: x[0] * matrix([[0,0,1],[0,1,0],[1,0,0]]) + x[1] * identity_matrix(3) >= 0
[0 0 0]    [x_1 0         x_0]
[0 0 0] <= [0   x_0 + x_1 0  ]
[0 0 0]    [x_0 0         x_1]

[vbraun]

Whats the point of the add_linear_constraints method in backends? Add multiple trivial constraints without a way to hand it any coefficient data? Really? ;-)

[ncohen]

Whats the point of the add_linear_constraints method in backends? Add multiple trivial constraints without a way to hand it any coefficient data? Really? ;-)

They also expose the solver's corresponding functions, which do the very same. And that's probably where we will get a speedup, for it seems that calling add_constraint n times instead of calling add_constraints(n) does make a difference.

Nathann

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1c5c024	Full doctest coverage
​fd10630	Clean up add_constraint()

[vbraun]

Replying to ncohen:

They also expose the solver's corresponding functions, which do the very same. And that's probably where we will get a speedup, for it seems that calling add_constraint n times instead of calling add_constraints(n) does make a difference.

Obviously it is very fast to solve trivial constraints. Read my question again.

[vbraun]

This:

    cpdef add_linear_constraint(self, constraints, lower_bound, upper_bound, name=*)
    cpdef add_linear_constraints(self, int number, lower_bound, upper_bound, names=*)

[ncohen]

Obviously it is very fast to solve trivial constraints. Read my question again.

Volker you make mistakes too, and if there is something I don't stand for long it is people giving me orders. If there is a misunderstanding, we will clear it.

What I said is that people (some I met, or myself in the past) noticed that it took time to add constraints to a LP. I am not talking of solving the LP itself. Just imagine that the matrix is stored as a dense contiguous matrix, and that adding a row means having to copy everything.

I don't know if this is how it is implemented, but on some occasions adding constraints one by one (each with the same number of nonzero coefficients and looking the same) takes a time which is not constant per call, i.e. calling add_constraint on a matrix that is already big is longer than on an empty matrix.

Which means that calling add_linear_constraints may make a difference. Which would also explain why the solvers contain this in their API, even though it seems useless at first.

I just tried it again by adding the same constraint "b[3]<=8" many many times but in this case the time seems to be constant. Perhaps there is something else missing.

Nathann

[vbraun]

I'm not trying to give you orders, I just want to know whether add_linear_constraints is supposed to work, and if not how you'd envision it to work. Create empty rows and then mutate later? Or is the constraint data just missing from the argspec? or what?

[ncohen]

I'm not trying to give you orders, I just want to know whether add_linear_constraints is supposed to work, and if not how you'd envision it to work. Create empty rows and then mutate later? Or is the constraint data just missing from the argspec? or what?

I guess that it is meant to allocate space for new rows, before changing the coefficient. As I said above I implemented this function only to expose those that are available in the solvers, and they do the very same thing (if not less).

GLPK: takes only a number of rows as an argument. Nothing else

glp_add_rows(self.lp, number)

Cplex: takes a number of rows and the corresponding upper/lower bound. Plus names

CPXnewrows(self.env, self.lp, number, bound, sense, rng, c_names if names else NULL) #

Looks like the others don't have this function implemented.

Nathann

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​08bc41a

Allow specifying constraints via vector-valued linear functions

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e7d48ef

Also document matrix * MIPVariable notation

[vbraun]

I'm proposing a new interface for vector-valued linear functions. This should be the most natural way to hand over multiple constraints in one call. By default, this just falls back to adding scalar constraints individually. In another ticket we can add special backend support for that.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0090edc

Improve documentation

[vbraun]

This is all that should go into this ticket, please review.

[dimpase]

typo:

implinicly -> implicitly (?) in linear_tensor.py

and there is also at least one instanciate* instead of instantiate* in the diff...

[dimpase]

Replying to ncohen:

add_linear_constraints is supposed to work, and if not how you'd envision it to work. Create empty rows and then mutate later? Or is the constraint data just missing from the argspec? or what?

I guess that it is meant to allocate space for new rows, before changing the coefficient. As I said above I implemented this function only to expose those that are available in the solvers, and they do the very same thing (if not `less).

GUROBI does have ​`GRBaddconstrs` which a full thing, in the sense you can specify a bunch of constraints, specified by a sparse matrix, at once.

GLPK: takes only a number of rows as an argument. Nothing else

glp_add_rows(self.lp, number)

Cplex: takes a number of rows and the corresponding upper/lower bound. Plus names

CPXnewrows(self.env, self.lp, number, bound, sense, rng, c_names if names else NULL) #

there is ​`CPXaddrows` in CPLEX, similar to the one in GUROBI.

[dimpase]

Replying to vbraun:

Inequalities with vector/matrix-valued linear functions:

sage: mip.<x> = MixedIntegerLinearProgram()
sage: x[0] * matrix([[0,0,1],[0,1,0],[1,0,0]]) + x[1] * identity_matrix(3) >= 0
[0 0 0]    [x_1 0         x_0]
[0 0 0] <= [0   x_0 + x_1 0  ]
[0 0 0]    [x_0 0         x_1]

is it expected that mip.show() shows some kind of nonsense after such inequalities are entered?

[vbraun]

Right now you can only add scalar and vector-valued linear constraints. I intended the matrix notation for semidefinite programming. But you can't do anything with them yet.

[dimpase]

Replying to vbraun:

Right now you can only add scalar and vector-valued linear constraints. I intended the matrix notation for semidefinite programming.

Great!

How does one recover the original matrices x[i]'s are multiplied with? I see that now each constraint gets .rhs which is a matrix of linear forms. For SDP backends one would typically want the original matrices. (Of course they can still be found...)

Well, this is probably for another ticket.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​039a835

Spelling and docstring fixes

[vbraun]

IMHO if you want to pass a dense matrix you should just implement a separate method for that. Since the point of the interface is to simplify the variable ordering/indexing I'm cutting the matrix into columns and only store the columns. The solver backend should then put them back into a matrix (or whatever storage the backend uses) in the backend's order.

[vbraun]

In the previous comment I was talking about vector-valued linear functions.

For matrix-valued linear functions you get the coefficient matrices in the usual way:

sage: m = x[0] * matrix([[0,0,1],[0,1,0],[1,0,0]])
sage: m.dict()[0]
[0.0 0.0 1.0]
[0.0 1.0 0.0]
[1.0 0.0 0.0]

[dimpase]

Replying to vbraun:

In the previous comment I was talking about vector-valued linear functions.

For matrix-valued linear functions you get the coefficient matrices in the usual way:

sage: m = x[0] * matrix([[0,0,1],[0,1,0],[1,0,0]])
sage: m.dict()[0]
[0.0 0.0 1.0]
[0.0 1.0 0.0]
[1.0 0.0 0.0]

can this be made explicit in a docstring? Otherwise looks good to me.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e2ddf8c

Simplify and document how to extract coefficients

[vbraun]

Replying to dimpase:

can this be made explicit in a docstring?

done

[dimpase]

is the syntax

blah.<foo>=MixedIntegerLinearProgram()

new? I can't recall seeing this before. It should be documented in mip.pyx, as an alternative to

blah=MixedIntegerLinearProgram()
foo=blah.new_variable()

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d3d0d10

Document generator syntax some more

[vbraun]

IMHO accepting generator syntax not really something that needs to be spelled out, but is something that should be expected to work. I've used it in most new doctests so its not like it is hidden away. In any case, added another note to new_variable.

[dimpase]

OK!

[ncohen]

Hello !

I was just trying to see how this new code works (Thanks Volker !!!) and see what we can earn by implementing solver-specific "add_constraints". And while playing, I found that :

sage: p.<x>  = MixedIntegerLinearProgram()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-300a6e09d28b> in <module>()
----> 1 p = MixedIntegerLinearProgram(names=('x',)); (x,) = p._first_ngens(1)

/home/ncohen/.Sage/local/lib/python2.7/site-packages/sage/numerical/mip.so in sage.numerical.mip.MixedIntegerLinearProgram.__init__ (build/cythonized/sage/numerical/mip.c:1740)()

TypeError: __init__() got an unexpected keyword argument 'names'

Nathann

[dimpase]

Replying to ncohen:

Hello !

I was just trying to see how this new code works (Thanks Volker !!!) and see what we can earn by implementing solver-specific "add_constraints". And while playing, I found that :

sage: p.<x>  = MixedIntegerLinearProgram()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-300a6e09d28b> in <module>()
----> 1 p = MixedIntegerLinearProgram(names=('x',)); (x,) = p._first_ngens(1)

/home/ncohen/.Sage/local/lib/python2.7/site-packages/sage/numerical/mip.so in sage.numerical.mip.MixedIntegerLinearProgram.__init__ (build/cythonized/sage/numerical/mip.c:1740)()

TypeError: __init__() got an unexpected keyword argument 'names'

hmm, it seems you must have something weird in the codebase in your local Sage install...

[ncohen]

Oh right. a touch numerical/* and sage -b solved it ! Thanks !

Nathann

[dimpase]

Here is the work of my MSc student on an interface for SDP solver, that builds upon this ticket: ​https://github.com/ingolfured/sageproject/

it has a backend to the SDP solver in CVXOPT.

[ncohen]

Yo Dima ! If you happen to write a LP-related ticket, could you do that ?

+++ b/src/sage/numerical/mip.pyx
@@ -1183,7 +1183,7 @@ cdef class MixedIntegerLinearProgram(SageObject):
             if b.is_variable_integer(i):
                 var_type = 'an integer'
             elif b.is_variable_binary(i):
-                var_type = 'a boolean variable'
+                var_type = 'a boolean'
             else:
                 var_type = 'a continuous'
             if varid_name[i] == str(self.gen(i)):

Otherwise we get stuff like that :

sage: p.show()
Maximization:
 
Constraints:
  x_0 <= 1.0
Variables:
  x_0 is a boolean variable variable (min=0.0, max=1.0)

And it does not seem worth a ticket of its own...

Nathann