用 numpy 求解整数上的线性系统答案

【问题标题】：Solving linear system over integers with numpy用 numpy 求解整数上的线性系统
【发布时间】：2012-12-03 14:36:36
【问题描述】：

我正在尝试用 numpy.目前，我正在做这样的事情（作为一个简单的例子）：

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

print np.linalg.lstsq(a,b)[0]

[1. 1.]

这可行，但使用浮点数。有没有办法只解决整数系统？我已经尝试过类似

print map(int, np.linalg.lstsq(a,b)[0])

[0, 1]

为了将解决方案转换为整数数组，期待[1, 1]，但显然我遗漏了一些东西。谁能指出我正确的方向？

【问题讨论】：

你真的应该澄清这个问题。文本状态超定，所以也许我解释错了，但这表明你不知道有一个精确的解决方案，我也看不出最佳/精确解决方案应该是整数的原因（如果你知道其中一个是一个巨大的缺失信息）。在这种情况下，这个例子过于简单了，因为它有一个最好的（甚至是精确的）整数解。

标签： python numpy linear-algebra

【解决方案1】：

您应该使用专门的整数问题求解器（请注意，整数问题甚至都不容易解决）。 openopt 是一个包，例如应该为整数二次优化提供良好的包装器，例如您正在做的。尝试使用线性代数根本不会直接为您提供正确的解决方案。

您的问题可以写成quadratic program，但它是integer，所以使用openopt 或其他模块。由于它是一种非常简单、不受约束的方法，也许还有其他方法。但对于初学者来说，这并不是一开始看起来那么简单的问题，openopt等中有程序准备好有效地解决这种事情。

【讨论】：

【解决方案2】：

您正在查看一个线性系统diophantine equations。一个快速的谷歌搜索出现了 Systems of Linear Diophantine Equations 费利克斯·拉泽布尼克。在那篇论文中，作者考虑了以下问题：

给定一个线性方程组 Ax = b，其中 A = a(i,j) 是一个 m × n 矩阵具有整数条目，并且 b 是具有整数分量的 m × 1 列向量，系统有一个整数解，即具有整数分量的 n × 1 解向量 x？

【讨论】：

【解决方案3】：

当您转换为int 时，元素的小数部分会被截断，因此会向下舍入。

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

x = np.linalg.lstsq(a,b)[0]

结果：

>>> x
array([ 1.,  1.])
>>> x[0]
0.99999999999999967
>>> x[1]
1.0000000000000002
>>> x.astype(int)
array([0, 1])
>>> map(int, x)
[0, 1]
>>> np.array([1.,1.]).astype(int) # works fine here
array([1, 1])

【讨论】：

我认为它是这样的 - 那么有没有办法在计算中只使用整数来避免这些浮点数不准确？
我认为np.round 应该这样做。
对不起，但这恰好在这里工作，因为对于这个例子，精确的解决方案恰好是整数。舍入/截断没有理由为您提供最佳整数解决方案。

【解决方案4】：

我可能误解了你的问题，但我认为你只需要 round 和 astype(int) 的组合？

例如

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

x = np.linalg.lstsq(a,b)[0]
print x.round().astype(int)

【讨论】：

谢谢，这也很有帮助。所以 numpy 没有专门用于解决整数上的线性系统的功能吗？我本来打算在万不得已的情况下使用这样的舍入方法，但它应该不会有太大的不同。
@A.R.S.：您希望使用整数算术求解线性代数，还是希望结果为整数数据类型？我认为前者是不可能的。对于后者，我认为这里给出的建议应该足够了。
我同意 tiango，从问题中我会说这是不正确的，因为四舍五入并不能保证找到最佳整数解决方案（除非它以整数开头），但系统在手，这似乎不太可能。尽管我承认，这个例子恰好有一个精确的解决方案。
@A.R.S. - 不，据我所知，numpy 没有任何方法可以求解保证产生整数解的方程组。在内部，numpy 只是从LAPACK 等调用例程，所以像lstsq 这样的函数本质上是浮点数。我认为在 NPE 确定是否存在解决方案的答案和 seberg 尝试非线性求解器的建议之间，您可能可以将一些东西放在一起。如果没有别的，一个舍入或截断的浮点解可能会为非线性求解器提供一个很好的开始猜测。

【解决方案5】：

+1 到 seberg，这里有一个反例来说明你不应该映射圆形：
（抱歉，是matlab风格，但是你很容易pythonize）

a =
     3     0
     0     3
     1     1
b = 
    2.71
   11.7
    0.5
x = a\b =
    0.5121
    3.5088
round(x) =
    1
    4
norm(a*round(x)-b) = 4.5193
norm(a*[0;4]-b) = 4.4367
norm(a*[1;3]-b) = 4.4299

【讨论】：

【解决方案6】：

我需要这样做并最终将 Keith Matthews 编写的 PHP 程序（您可以在 http://www.numbertheory.org/php/php.html 上找到）移植到 Python 中。我最初使用 Numpy 数组，但遇到整数溢出问题，因此切换到使用任意精度数值表示的 Sympy 矩阵。

代码在 GitHub 上发布，地址为 https://github.com/tclose/Diophantine，遵循 MIT 许可证，因此请随时使用它，如果您遇到任何问题，请告诉我（抱歉，它没有更好的文档记录）。 master 分支使用 Sympy，但您可以在“numpy”分支中访问原始 Numpy 实现，这对于相当稀疏的系统似乎可以正常工作。

如果您最终将它用于科学出版物，请引用 Keith 的论文（并可能添加指向 GitHub 存储库的链接）。

【讨论】：

【解决方案7】：

我的方法是先找到非整数解，然后放大到整数一

from fractions import Fraction, gcd
from functools import reduce

def lcm(a, b):
    return a * b // gcd(a, b)

def common_int(*numbers):
    fractions = [Fraction(n).limit_denominator() for n in numbers[0]]
    multiple = reduce(lcm, [f.denominator for f in fractions])
    ints = [f * multiple for f in fractions]
    divisor = reduce(gcd, ints)

    return [int(n / divisor) for n in ints]

sol = np.linalg.solve(np.array([[1, 2, 3], [2, 1, 0], [2, 1, 4]]), np.array([1., 1., 1.]))  # system of equation
# [0.3333333, 0.3333333, 0.]

common_int(sol)
# [1., 1., 0.]

【讨论】：

【解决方案8】：

有一种方法叫做block lanczos. 这可以让你在有限域上回答。您可以找到针对此特定问题的块 lanczos 求解器。

【讨论】：