如何求解非负整数上的线性方程组？

Question

如何求解非负整数上的线性方程组？

pythonmatrixlinear-algebrasympyinteger-programming

9

给定一个线性系统 Ax = b，其中矩阵 A 和向量 b 均具有整数值，我希望找到所有解决此方程的 非负整数 向量 x。

到目前为止，我已经发现了一些技巧，例如使用矩阵的 Smith正规形式或 Hermite正规形式来寻找整数解，然后可以使用线性求解器来找到非负解。是否存在库可以更轻松地实现这一点？

Python 解决方案是理想的，但如果其他语言中存在库，我也想知道。

- while1fork

你尝试过使用Numpy或SciPy吗？ - Julio Daniel Reyes

@Julio Daniel Reyes 这两个工具可以使用浮点数解决线性系统并进行优化，但据我所知，它们不能解决我在问题中询问的特定事情。 - while1fork

1

如果A是可逆的，那么只有一个解，你可以找到并检查是否为全整数。但我猜你的矩阵是不可逆的（可能甚至不是方阵）。 - Tom Karzes

1

相关：获取线性方程的所有正整数解，使用numpy/sympy找到线性系统的整数解（丢番图方程）。 - user6655984

实际上，你的问题涉及到线性代数中的一个问题。如果一个矩阵方程有解，那么它要么是唯一解，要么有无限多个解。因此，在第二种情况下，将所有解表示为向量而不是公式是有问题的。 - MaximTitarenko

4个回答

3

如果至少存在一个非负整数解Ax=b，那么整数规划应该可以找到它。所有非负整数解的集合可以通过A的零空间找到。

例子

使用 Erwin 的答案中的 A 和 b:

>>> from sympy import *
>>> A = Matrix([[ 1, 2, 1],
                [ 3, 1,-1]])
>>> b = Matrix([20,12])

计算零空间：

>>> A.nullspace()
[Matrix([
[ 3/5],
[-4/5],
[   1]])]

空间的零空间是一维的，其基向量为有理数。使用Erwin找到的特定解(2,8,2)，所有实数解的集合是一个线，参数为(2,8,2) + t (3,-4,5)。由于我们只对非负整数解感兴趣，因此我们将该线与非负八卦面相交，然后再与整数晶格相交，这会产生t ∈ {0,1,2}。因此，如果：

t=1，我们得到解(5,4,7)。
t=2，我们得到解(8,0,12)。

请注意，这些是Erwin使用Z3找到的3个解决方案。

然而，当零空间的维度高于1时，情况就变得更加困难。看一下:

Jesús De Loera, 计算多面体中点格点数的许多方面, 2005年。

- Rodrigo de Azevedo

非常好。谢谢。 - Erwin Kalvelagen

2

这是一个Z3模型。 Z3 是微软的定理证明器。该模型与之前提出的MIP模型非常相似。

在MIP中枚举整数解并不完全是微不足道的（可以通过一些努力 [链接] 或使用高级MIP求解器中的“解池”技术来完成）。使用Z3会更容易一些。甚至更容易的方法是使用约束编程（CP）求解器：它们可以自动枚举解决方案。

我们开始吧：

from z3 import *
A = [[1, 2, 1], [3, 1, -1]]
b = [20, 12]
n = len(A[0]) # number of variables
m = len(b)    # number of constraints

X = [ Int('x%d' % i) for i in range(n) ]
s = Solver()
s.add(And([ X[i] >= 0 for i in range(n) ]))
for i in range(m):
    s.add( Sum([ A[i][j]*X[j] for j in range(n) ]) == b[i])

while s.check() == sat:
    print(s.model())
    sol = [s.model().evaluate(X[i]) for i in range(n)]
    forbid = Or([X[i] != sol[i] for i in range(n)])
    s.add(forbid)

它解决了什么问题。

x0+2x1+x2 = 20
3x0+x1-x2 = 12

解决方案如下：

[x2 = 2, x0 = 2, x1 = 8]
[x2 = 7, x1 = 4, x0 = 5]
[x2 = 12, x1 = 0, x0 = 8]

我们可以打印最终模型，以便查看如何添加“禁止”约束：

[And(x0 >= 0, x1 >= 0, x2 >= 0),
 1*x0 + 2*x1 + 1*x2 == 20,
 3*x0 + 1*x1 + -1*x2 == 12,
 Or(x0 != 2, x1 != 8, x2 != 2),
 Or(x0 != 5, x1 != 4, x2 != 7),
 Or(x0 != 8, x1 != 0, x2 != 12)]

- Erwin Kalvelagen

如果A的零空间包含非负半轴上的有理数射线，会怎么样呢？在这种情况下，如果存在一个整数非负解，那么就会存在无限多个整数非负解。Z3如何处理这种情况？它会超时吗？ - Rodrigo de Azevedo

1

它会表现出人们所预期的行为：循环将不会终止（好吧，在某些时候你会耗尽内存），并且不断产生新的解决方案。只需添加一个计数器，并在看到N个解决方案时停止即可。 - Erwin Kalvelagen

0

请查看这里的4ti2包中的zsolve方法。

- user124697

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- josliber · Accepted Answer

您可以使用整数规划来实现这一点，为您拥有的每个x值定义一个非负整数决策变量。然后，您可以使用约束条件Ax=b和目标函数0来解决问题，该函数搜索可行的整数解来解决您的方程组。

这可以很容易地在Python的pulp包中实现。例如，考虑解决以下系统:

x+2y+z = 5
3x+y-z = 5

你可以使用pulp来解决这个问题：

import pulp
A = [[1, 2, 1], [3, 1, -1]]
b = [5, 5]
mod = pulp.LpProblem('prob')
vars = pulp.LpVariable.dicts('x', range(len(A[0])), lowBound=0, cat='Integer')
for row, rhs in zip(A, b):
    mod += sum([row[i]*vars[i] for i in range(len(row))]) == rhs
mod.solve()
[vars[i].value() for i in range(len(A[0]))]
# [1.0, 2.0, 0.0]

这表明存在一组可行的整数解，x=1，y=2，z=0。