在 Python 中求解线性整数方程组答案

【问题标题】：Solve system of linear integer equations in Python在 Python 中求解线性整数方程组
【发布时间】：2017-07-27 01:05:55
【问题描述】：

我正在寻找一种在 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)。这也是正确的，但不是理想的解决方案...

编辑：我很抱歉在某些方面不够精确。如果有人对细节感兴趣，问题来自论文 “用于数字信号处理的同步数据流程序的静态调度”，Edward A. Lee 和 David G. Messerschmitt， IEEE 计算机汇刊，卷。 C-36，第 1 期，第 24-35 页，1987 年 1 月。

定理 2 说

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

在证明定理 2 之后他们说

可能需要求解零空间中的最小正整数向量。为此，减少 u' 中的每个有理项，使其分子和分母互质。欧几里得算法适用于此。

【问题讨论】：

“最小整数向量”是什么意思？
“最小”是如何定义的，(1, 1, 2) 是如何小于(0, 0, 0) 的？
一般来说，这可能是really difficult，我认为 NumPy 或 SciPy 都没有提供任何工具来做到这一点。
您已经编辑了您特别想排除零向量的事实，但您仍然没有定义“最小”。 (2, 0, 0) 是否小于 (0, 1, 1)？ (1, 0, 0) 是否小于 (-1, 0, 0)？ “最小”是模棱两可的。
A 等级不足，您建议的理想 x 实际上是空空间。我们这里有XY problem 吗？

标签： python numpy scipy linear-algebra

【解决方案1】：

这也许是一个疯狂的想法，但听起来您正在寻找构建一个约束求解器。

minizinc 是一个通用的约束求解器。也许可以用 minizinc 可以解决的方式来表达你的约束？

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

【讨论】：

CP 确实在与目标函数作斗争。但是如果没有对smallest 的更正式的描述，就很难回答了。

【解决方案2】：

pypi 有一个求解此类方程的求解器。它显然可以计算矩阵的Hermite normal form，而矩阵又可以用来解决您的问题。

软件包版本为 0.1。

Sage 也出现在 support 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 - >

【讨论】：

The returned solution vector will tend to be one with the smallest norms.

【解决方案3】：

正如在 cmets 中看到的那样，这个问题非常不正式。在不知道您对最小的定义（例如：l1-norm、l2-norm）的情况下，很难回答您的具体问题。

一般问题与求解丢番图方程组有关，但这些问题很难求解（在一般情况下），而且软件不多。

一种自然的方法是使用 整数编程，这是 NP 难的，但商业和一些开源求解器非常强大。

numpy/scipy 中没有内置方法可以在不进行大量修改的情况下解决您的问题（例如：基于 numpy/scipy 实现一些自己的算法）。可悲的是，numpy/scipy 中也没有 IP-solver。

假设：

smallest 是变量的总和（大多数时候这不是正确的方法；l2-norm 更流行，但它有点难以制定，它需要更强大的求解器）
你想禁止零向量
x is nonnegative

这里是一些使用纸浆和 numpy 的简单的基于 IP 的实现。我不太喜欢纸浆，但它很容易在各种流行的系统上安装 (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]

【讨论】：

【解决方案4】：

要找到您想要的确切解决方案，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])

【讨论】：

sympy 似乎永远不会返回其元素共享一个因子的空空间向量（例如 (3, 3, 6)，这会使您的算法出错）。我想知道这是保证还是实现细节。在任何情况下，我都不建议您依赖它，因为它可以以几乎免费的方式获得更安全的解决方案。

【解决方案5】：

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

>>> 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 的特征值。因此，0 特征值的特征向量 s 非常适合您的方程系统。

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

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

【讨论】：

设置x = s / s[1]，仔细看看x[0]。我得到0.999999999999999，而不是1。但这可能足以解决问题。这取决于 Apoptose 下一步将如何处理结果。
但是如果 x 的分量是例如三个不同的素数呢？你将如何重建它们？

【解决方案6】：

这是一个使用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])

【讨论】：