在Python中解决线性整数方程组

Question

我正在寻找在Python中解决线性方程组的方法。特别地，我正在寻找比所有零向量更大且能够解决给定方程的最小整数向量。例如，我有以下方程：

想要解决以下问题： 。

在此情况下，解决此方程的最小整数向量为 。

然而，我该如何自动确定此解？如果我使用scipy.optimize.nnls，例如：

A = np.array([[1,-1,0],[0,2,-1],[2,0,-1]])
b = np.array([0,0,0])
nnls(A,b)

结果是(array([ 0., 0., 0.]), 0.0)。这也是正确的，但不是期望的解决方案...

编辑：对于某些方面的不准确，我深表歉意。 如果有人对细节感兴趣，问题来自论文《数字信号处理中同步数据流程序的静态调度》(Static Scheduling of Synchronous Data Flow Programs for Digital Signal Processing)，作者为 Edward A. Lee 和 David G. Messerschmitt，IEEE 计算机学报，Vol. C-36，No. 1，pp. 24-35，1987 年 1 月。

定理 2 表示：

对于具有拓扑矩阵 A、秩(A)=s-2 和 s 个节点的连通 SDF 图，我们可以找到一个正整数向量 b != 0，使得 Ab = 0，其中 0 是零向量。

在定理 2 的证明之后，他们说：

可能希望解决零空间中最小的正整数向量。为此，将 u' 中的每个有理数条目约简，使其分子和分母互质。欧几里得算法可用于此。

- Apoptose

---

“the smallest integer vector”是什么意思？ - Code-Apprentice

“smallest” 是如何定义的，(1, 1, 2) 怎么比 (0, 0, 0) 更小？ - user2357112

一般来说，这可能非常困难，我认为NumPy或SciPy都没有提供任何工具来完成它。 - user2357112

你已经编辑了要排除零向量的事实，但你仍然没有定义“最小”。(2, 0, 0)比(0, 1, 1)小吗？(1, 0, 0)比(-1, 0, 0)小吗？“最小”是含糊不清的。 - user2357112

A是秩不足的，你提出的x其实是它的零空间。这里是否存在一个“XY问题”？ - Ahmed Fasih

Answer 1

为了找到您想要的确切解决方案，numpy和scipy可能不是最好的工具。它们的算法通常使用浮点数，不能保证给出确切的答案。

您可以使用sympy来获得此问题的确切答案。在sympy中，Matrix类提供了nullspace()方法，返回零空间的基向量列表。以下是一个示例：

In [20]: from sympy import Matrix, lcm

In [21]: A = Matrix([[1, -1, 0], [0, 2, -1], [2, 0, -1]])

获取零空间中的向量。 nullspace() 返回一个列表；此代码假设 A 的秩为2，因此该列表长度为1：

In [22]: v = A.nullspace()[0]

In [23]: v
Out[23]: 
Matrix([
[1/2],
[1/2],
[  1]])

找出 v 中所有数的分母的最小公倍数，以便我们可以将结果缩放为最小的整数：

In [24]: m = lcm([val.q for val in v])

In [25]: m
Out[25]: 2

x保存整数答案：

In [26]: x = m*v

In [27]: x
Out[27]: 
Matrix([
[1],
[1],
[2]])

要将该结果转换为一个整数的numpy数组，您可以执行以下操作：

In [52]: np.array([int(val) for val in x])
Out[52]: array([1, 1, 2])

- Warren Weckesser

---

sympy 似乎从未返回过元素共享因子（例如 (3, 3, 6)）的零空间向量，这可能会使您的算法出错。我想知道这是一个保证还是实现细节。无论如何，当有更安全的解决方案可以轻松获得时，我不建议依赖它。 - Paul Panzer

Answer 2

实际上，这只是基本的线性代数。

>>> A = np.array([[1,-1,0], [0,2,-1],[2,0,-1]])    

让我们计算这个矩阵的特征值和特征向量。

>>> 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])

四舍五入后，我们发现0是矩阵A的一个特征值。因此，特征向量s对于0特征值的方程组是一个好的候选。

>>> s = e[1][:,1]
>>> s
array([-0.40824829, -0.40824829, -0.81649658])    

将特征向量与相关矩阵相乘，会得到特征向量本身，乘以相关特征值的比例。因此，在上述情况下，我们可以看到：As = 0s = 0

>>> np.round(A.dot(s), 10)
array([ 0.,  0.,  0.])

既然我们对整数解感兴趣，那么我们就需要对解向量进行缩放。

>>> x = s / s[1]
>>> x
array([ 1.,  1.,  2.])

希望这个答案能解决你的问题。

- Berkan Erol

---

设置 x = s / s [1]，并仔细查看 x [0]。我得到的是 0.999999999999999，而不是 1。但这可能已经足够解决问题了。这取决于Apoptose接下来要对结果做什么。 - Warren Weckesser

但是如果x的组成部分例如是三个不同的质数，你该如何重构它们呢？ - Paul Panzer

Answer 3

以下是使用ortools的解决方案:

s = Solver('')
A = [[1, -1, 0], [0, 2, -1], [2, 0, -1]]
vars = [s.IntVar(1, 2**32) for _ in A[0]]
for row in A:
    s.Add(s.ScalProd(vars, row) == 0)
s.NewSearch(s.Phase(vars, s.CHOOSE_FIRST_UNBOUND, s.ASSIGN_MIN_VALUE))
s.NextSolution()
print([var.Value() for var in vars])

- Björn Lindqvist

Answer 4

这种类型的方程有一个求解器，可以在pypi找到。它显然可以计算矩阵的Hermite normal form，从而解决您的问题。

该软件包的版本为0.1。

Sage似乎也支持Hermite正规形式。

同齐次系统的特殊情况，即b=0的情况相比较而言更容易一些。这里有一个求解器，适用于矩阵仅具有秩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

样例输出：

>>> 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 - >

- Paul Panzer

---

返回的解向量往往是具有最小范数的向量。 - sascha

Answer 5

这也许是一个大胆的想法，但听起来好像你正在寻求构建一个约束求解器。

minizinc 是一个通用的约束求解器。也许可以以这样的方式表达您的约束条件，使minizinc能够解决它？

然后似乎有一个Python库与之接口：https://pypi.python.org/pypi/pymzn/

- Irmen de Jong

---

CP在目标函数方面确实很困难。但是如果没有更正式的“smallest”描述，很难回答。 - sascha

Answer 6

这个问题在评论中看起来相当非正式。如果不知道您对最小值的定义（例如：l1-范数，l2-范数），很难回答您的具体问题。

一般问题与解决丢番图方程组有关，但这些问题很难解决（在一般情况下），并且没有太多的软件。

一种自然的方法是使用整数规划，它是NP-hard的，但商业和一些开源求解器非常强大。

在numpy/scipy中没有内置的方法可以解决您的问题，而不需要进行巨大的修改（例如：基于numpy/scipy实现一些自己的算法）。遗憾的是，在numpy/scipy中也没有IP求解器。

让我们假设：

• smallest 是变量之和（大多数情况下这不是正确的方法；l2-范数更受欢迎，但它更难以表述，并且需要更强大的求解器）
• 您想禁止零向量
• x为非负数

这里是使用pulp和numpy的一些简单的基于IP的实现。我不太喜欢pulp，但它很容易在各种流行的系统上安装（pip install pulp）。

代码

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])

输出

Optimal
[1.0, 1.0, 2.0]

- sascha