Solve system of linear integer equations in Python
I am looking for a way to solve a system of linear equations in Python. In particular, I am looking for the smallest integer vector that is greater than all zeros and solves the given equation. For example, I have the following equation:
and want to decide
.
In this case, the smallest integer vector that solves this equation is
.
However, how can I detect this decision automatically? If I use scipy.optimize.nnls
for example
A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]]) b = np.array([0,0,0]) nnls(A,b)
result (array([ 0., 0., 0.]), 0.0)
. Which is also correct, but not the desired solution ...
Edit: I apologize for being inaccurate in certain aspects. If anyone is interested in the details, the problem comes from "Static Scheduling of Synchronous Data Flow Programs for Digital Signal Processing", Edward A. Lee and David G. Messerschmitt, IEEE Transactions on Computers, Vol. C-36, No. 1, pp. 24-35, January 1987
Theorem 2 says that
For a connected SDF graph with s-nodes and a matrix of topology A and rank (A) = s-2, one can find an integer vector b! = 0 such that Ab = 0, where 0 is a zero vector.
Immediately after the proof of Theorem 2, it is said that
It may be desirable to solve for the smallest positive integer vector in zero space. To do this, decrease each rational entry in u 'so that its numerator and denominator are coprime. For this, the Euclidean algorithm will work.
source to share
To find the exact solution you want, numpy and scipy are probably not the best tools. Their algorithms usually work in floating point and do not guarantee an exact answer.
You can use sympy
to get the exact answer to this problem. Class Matrix
c sympy
provides a method nullspace()
that returns a list of base vectors for nullspace. Here's an example:
In [20]: from sympy import Matrix, lcm
In [21]: A = Matrix([[1, -1, 0], [0, 2, -1], [2, 0, -1]])
Get a vector in nullspace. nullspace()
returns a list; this code assumes rank of A is 2, so the list is one length:
In [22]: v = A.nullspace()[0]
In [23]: v
Out[23]:
Matrix([
[1/2],
[1/2],
[ 1]])
Find the smallest common multiple of the denominators of the values in v
so we can scale the result to the smallest integers:
In [24]: m = lcm([val.q for val in v])
In [25]: m
Out[25]: 2
x
contains an integer answer:
In [26]: x = m*v
In [27]: x
Out[27]:
Matrix([
[1],
[1],
[2]])
To convert this result to a numpy array of integers, you can do something like:
In [52]: np.array([int(val) for val in x])
Out[52]: array([1, 1, 2])
source to share
Actually, it's just basic linear algebra.
>>> A = np.array([[1,-1,0], [0,2,-1],[2,0,-1]])
Let the eigenvalues and the eigenvector be calculated for this matrix.
>>> e = np.linalg.eig(A)
>>> e
(array([ -4.14213562e-01, -1.05471187e-15, 2.41421356e+00]), array([[-0.26120387, -0.40824829, 0.54691816],
[-0.36939806, -0.40824829, -0.77345908],
[-0.89180581, -0.81649658, 0.32037724]]))
>>>> np.round(e[0], 10)
array([-0.41421356, -0. , 2.41421356])
After rounding off, we see that 0 is the eigenvalue of our matrix A. So the eigenvector s for the 0-eigenvalue is a good candidate for your system of equations.
>>> s = e[1][:,1]
>>> s
array([-0.40824829, -0.40824829, -0.81649658])
Multiplying an eigenvector by an associated matrix results in the eigenvector itself, scaled by the corresponding eigenvalue. So in the above case we see: As = 0s = 0
>>> np.round(A.dot(s), 10)
array([ 0., 0., 0.])
Since we are interested in an integer solution, we must scale the solution vector.
>>> x = s / s[1]
>>> x
array([ 1., 1., 2.])
Hope this answer solves your problem.
source to share
This might be a wild idea, but it sounds like you want to create a constraint solver.
minizinc is a general purpose constraint solver. Perhaps you can express your constraint in such a way that minizinc can solve it?
Then the Python library appears to interact with it: https://pypi.python.org/pypi/pymzn/
source to share
This question is rather informal as you can see from the comments. Without knowing your smallest definition (examples: l1-norm, l2-norm), it is difficult to answer your specific problem.
A general problem is related to solving a system of Diophantine equations, but they are difficult to solve (in general) and there is little software.
The natural approach is to use Integer programming , which is NP-hard, but commercial and some open source solutions are very powerful.
There is no built-in method in numpy / scipy that solves your problem without huge changes (eg: implementing your own algorithm based on numpy / scipy). Unfortunately, the IP solver is not within numpy / scipy either.
Let's assume:
- the smallest sum is the sum of the variables (this is the wrong approach most of the time; l2-norm is more popular, but a little more difficult to formulate and needs more powerful solvers)
- you want to ban the null vector
-
x is nonnegative
Here is some simple IP based implementation using cellulose and numpy. I don't really like cellulose, but it's easy to install ( pip install pulp
) on all popular systems.
Code
from pulp import *
import numpy as np
EPS = 1e-3
""" Input """
A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])
""" MIP """
# Variables
x = np.empty(b.shape[0], dtype=object)
for i in range(b.shape[0]):
x[i] = LpVariable("x" + str(i), lowBound=0, upBound=None, cat='Integer')
# Problem
prob = LpProblem("prob", LpMinimize)
# Objective
prob += np.sum(x)
# Constraints
for row in range(A.shape[0]):
prob += np.dot(A[row], x) == b[row]
prob += np.sum(x) >= EPS # forbid zero-vector
# Solve
status = prob.solve()
print(LpStatus[status])
print([value(x_) for x_ in x])
Output
Optimal
[1.0, 1.0, 2.0]
source to share
There is a solution for such an equation in pypi . Apparently it can compute the Errite normal form of a matrix, which in turn can be used to solve your problem.
The package is version 0.1.
Sage also appears to support Hermite's normal form.
A special case of a homogeneous system, i.e. b = 0, a little simpler, here is a solver for the simplest possible case of a matrix of rank n-1
import sympy
import numpy as np
def create_rd(n, defect=1, range=10):
while True:
res = sympy.Matrix((np.random.randint(-range+1,range,(n, n-defect))
@np.random.randint(0,2,(n-defect, n)))
.astype(object))
if res.rank() == n-defect:
break
return res
def solve(M):
ns = M.nullspace()
ns = [n / sympy.gcd(list(n)) for n in ns]
nsnp = np.array([[int(k) for k in n] for n in ns])
if len(ns) == 1:
return ns[0], nsnp[0]
else:
raise NotImplementedError
Output example:
>>> M = create_rd(4) # creates a rank-deficient matirx
>>> ns, nn = solve(M) # finds the 1d nullspace and a minimal integer basis vector
>>> M
Matrix([
[-7, -7, -7, -12],
[ 0, 6, 0, 6],
[ 4, 1, 4, -3],
[-4, -7, -4, -9]])
>>> ns
Matrix([
[-1],
[ 0],
[ 1],
[ 0]])
>>> M*ns
Matrix([
[0],
[0],
[0],
[0]])
>>> M = create_rd(40) # we can go to higher dimensions
>>> ns, nn = solve(M) # but solutions quickly become unwieldy
>>> ns
Matrix([
[ 8620150337],
[-48574455644],
[ -6216916999],
[-14703127270],
< - snip - >
source to share