#16714 closed enhancement (fixed)
Add a matrix of constraints in a LP
Reported by:  ncohen  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  linear programming  Keywords:  
Cc:  vbraun, dimpase, ingolfured, novoselt  Merged in:  
Authors:  Volker Braun  Reviewers:  Dima Pasechnik 
Report Upstream:  N/A  Work issues:  
Branch:  d3d0d10 (Commits)  Commit:  
Dependencies:  Stopgaps: 
Description
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
Change History (59)
comment:1 Changed 6 years ago by
 Cc ingolfured added
comment:2 Changed 6 years ago by
Looks like a duplicate of #7293... (which has no code and no comments)
comment:3 followups: ↓ 4 ↓ 5 Changed 6 years ago by
One technical problem is that matrices in Sage can always be multiplied, so you can't create matrices over linear functions.
comment:4 in reply to: ↑ 3 Changed 6 years ago by
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
comment:5 in reply to: ↑ 3 Changed 6 years ago by
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...
comment:6 Changed 6 years ago by
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)
.
comment:7 Changed 6 years ago by
 Branch set to u/vbraun/add_a_matrix_of_constraints_in_a_lp
comment:8 followups: ↓ 9 ↓ 12 Changed 6 years ago by
 Commit set to 846cba629c68e46f274de8d6cbbaabd3cdbb0a99
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

comment:9 in reply to: ↑ 8 ; followup: ↓ 10 Changed 6 years ago by
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
comment:10 in reply to: ↑ 9 ; followup: ↓ 11 Changed 6 years ago by
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.
comment:11 in reply to: ↑ 10 ; followup: ↓ 13 Changed 6 years ago by
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
comment:12 in reply to: ↑ 8 Changed 6 years ago by
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_{n1})  x_0^2 >= sum_{j=1}^{n1} x_j^2}
, known as "icecream cone") you have a partial order <_K
on R^n
defined by x<_K y
iff yk 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
).
comment:13 in reply to: ↑ 11 ; followup: ↓ 14 Changed 6 years ago by
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...
comment:14 in reply to: ↑ 13 ; followup: ↓ 15 Changed 6 years ago by
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
comment:15 in reply to: ↑ 14 ; followup: ↓ 16 Changed 6 years ago by
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.
comment:16 in reply to: ↑ 15 ; followup: ↓ 17 Changed 6 years ago by
Come on, it is implemented e.g. here: networkx
No, that is Edmond's algorith for general graphs (it is not bipartitespecific) and it is also the default algorithm called by Graph.matching
in case. Plus it's pure Python.
Nathann
comment:17 in reply to: ↑ 16 Changed 6 years ago by
Replying to ncohen:
Come on, it is implemented e.g. here: networkx
No, that is Edmond's algorith for general graphs (it is not bipartitespecific) 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 LPbased methods hopelessly lose here. I don't make baseless claims here, I know well what I am talking about...
comment:18 Changed 6 years ago by
 Commit changed from 846cba629c68e46f274de8d6cbbaabd3cdbb0a99 to 04f9ce174c2c7f9b205dfb1d5ef67710d9553f15
Branch pushed to git repo; I updated commit sha1. New commits:
04f9ce1  Add a parent for MIPVariables

comment:19 Changed 6 years ago by
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...
comment:20 Changed 6 years ago by
 Commit changed from 04f9ce174c2c7f9b205dfb1d5ef67710d9553f15 to 0a09269b45ec12f7dbe2567bb403db7eee6debdb
comment:21 Changed 6 years ago by
 Commit changed from 0a09269b45ec12f7dbe2567bb403db7eee6debdb to c09be57f5a5b80d6e9b45735cdf15502f0c2182b
Branch pushed to git repo; I updated commit sha1. New commits:
c09be57  Match x_i print order

comment:22 Changed 6 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:23 Changed 6 years ago by
 Commit changed from c09be57f5a5b80d6e9b45735cdf15502f0c2182b to a94e77e87c4acd07252f7dfdf31095e55cc61309
Branch pushed to git repo; I updated commit sha1. New commits:
a94e77e  Symbolic (in)equality involving linear tensors

comment:24 followup: ↓ 40 Changed 6 years ago by
Inequalities with vector/matrixvalued 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]
comment:25 followup: ↓ 26 Changed 6 years ago by
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? ;)
comment:26 in reply to: ↑ 25 ; followup: ↓ 28 Changed 6 years ago by
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
comment:27 Changed 6 years ago by
 Commit changed from a94e77e87c4acd07252f7dfdf31095e55cc61309 to fd10630a8e4884d3a44f55cc6e1d4c83ea7ee154
comment:28 in reply to: ↑ 26 ; followup: ↓ 30 Changed 6 years ago by
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.
comment:29 Changed 6 years ago by
This:
cpdef add_linear_constraint(self, constraints, lower_bound, upper_bound, name=*) cpdef add_linear_constraints(self, int number, lower_bound, upper_bound, names=*)
comment:30 in reply to: ↑ 28 Changed 6 years ago by
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
comment:31 followup: ↓ 32 Changed 6 years ago by
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?
comment:32 in reply to: ↑ 31 ; followup: ↓ 39 Changed 6 years ago by
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
comment:33 Changed 6 years ago by
 Commit changed from fd10630a8e4884d3a44f55cc6e1d4c83ea7ee154 to 08bc41aec5233abd6c0f215b6875c38bcd80aca5
Branch pushed to git repo; I updated commit sha1. New commits:
08bc41a  Allow specifying constraints via vectorvalued linear functions

comment:34 Changed 6 years ago by
 Commit changed from 08bc41aec5233abd6c0f215b6875c38bcd80aca5 to e7d48ef775732ce5e3f72ea976116c19888e7d26
Branch pushed to git repo; I updated commit sha1. New commits:
e7d48ef  Also document matrix * MIPVariable notation

comment:35 Changed 6 years ago by
I'm proposing a new interface for vectorvalued 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.
comment:36 Changed 6 years ago by
 Commit changed from e7d48ef775732ce5e3f72ea976116c19888e7d26 to 0090edc953f5930e0a10b0c981e99c317db6c27c
Branch pushed to git repo; I updated commit sha1. New commits:
0090edc  Improve documentation

comment:37 Changed 6 years ago by
 Status changed from new to needs_review
This is all that should go into this ticket, please review.
comment:38 Changed 6 years ago by
typo:
implinicly > implicitly (?) in linear_tensor.py
and there is also at least one instanciate* instead of instantiate* in the diff...
comment:39 in reply to: ↑ 32 Changed 6 years ago by
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.
comment:40 in reply to: ↑ 24 Changed 6 years ago by
Replying to vbraun:
Inequalities with vector/matrixvalued 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?
comment:41 followup: ↓ 42 Changed 6 years ago by
Right now you can only add scalar and vectorvalued linear constraints. I intended the matrix notation for semidefinite programming. But you can't do anything with them yet.
comment:42 in reply to: ↑ 41 Changed 6 years ago by
Replying to vbraun:
Right now you can only add scalar and vectorvalued 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.
comment:43 Changed 6 years ago by
 Commit changed from 0090edc953f5930e0a10b0c981e99c317db6c27c to 039a835cd4b2a47fb3fe3ef0670a687f8b7ee3e1
Branch pushed to git repo; I updated commit sha1. New commits:
039a835  Spelling and docstring fixes

comment:44 Changed 6 years ago by
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.
comment:45 followup: ↓ 46 Changed 6 years ago by
In the previous comment I was talking about vectorvalued linear functions.
For matrixvalued 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]
comment:46 in reply to: ↑ 45 ; followup: ↓ 49 Changed 6 years ago by
Replying to vbraun:
In the previous comment I was talking about vectorvalued linear functions.
For matrixvalued 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.
comment:47 Changed 6 years ago by
 Cc novoselt added
comment:48 Changed 6 years ago by
 Commit changed from 039a835cd4b2a47fb3fe3ef0670a687f8b7ee3e1 to e2ddf8c74a71662934b9d095c04b58c33999740a
Branch pushed to git repo; I updated commit sha1. New commits:
e2ddf8c  Simplify and document how to extract coefficients

comment:49 in reply to: ↑ 46 Changed 6 years ago by
comment:50 Changed 6 years ago by
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()
comment:51 Changed 6 years ago by
 Commit changed from e2ddf8c74a71662934b9d095c04b58c33999740a to d3d0d108142add4ef9aa3c441079715327aa77ab
Branch pushed to git repo; I updated commit sha1. New commits:
d3d0d10  Document generator syntax some more

comment:52 Changed 6 years ago by
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
.
comment:53 Changed 6 years ago by
 Reviewers set to Dima Pasechnik
 Status changed from needs_review to positive_review
OK!
comment:54 Changed 6 years ago by
 Branch changed from u/vbraun/add_a_matrix_of_constraints_in_a_lp to d3d0d108142add4ef9aa3c441079715327aa77ab
 Resolution set to fixed
 Status changed from positive_review to closed
comment:55 followup: ↓ 56 Changed 6 years ago by
 Commit d3d0d108142add4ef9aa3c441079715327aa77ab deleted
Hello !
I was just trying to see how this new code works (Thanks Volker !!!) and see what we can earn by implementing solverspecific "add_constraints". And while playing, I found that :
sage: p.<x> = MixedIntegerLinearProgram()  TypeError Traceback (most recent call last) <ipythoninput14300a6e09d28b> in <module>() > 1 p = MixedIntegerLinearProgram(names=('x',)); (x,) = p._first_ngens(1) /home/ncohen/.Sage/local/lib/python2.7/sitepackages/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
comment:56 in reply to: ↑ 55 Changed 6 years ago by
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 solverspecific "add_constraints". And while playing, I found that :
sage: p.<x> = MixedIntegerLinearProgram()  TypeError Traceback (most recent call last) <ipythoninput14300a6e09d28b> in <module>() > 1 p = MixedIntegerLinearProgram(names=('x',)); (x,) = p._first_ngens(1) /home/ncohen/.Sage/local/lib/python2.7/sitepackages/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...
comment:57 Changed 6 years ago by
Oh right. a touch numerical/*
and sage b
solved it ! Thanks !
Nathann
comment:58 Changed 6 years ago by
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.
comment:59 Changed 6 years ago by
Yo Dima ! If you happen to write a LPrelated 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
This would also be very natural syntax to specify polyhedra.