7

我正在尝试用 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]，但显然我错过了一些东西。谁能指出我正确的方向？

4

8 回答 8

7

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

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

于 2012-12-16T10:45:28.180 回答
4

您正在查看一个线性丢番图方程组。快速的谷歌搜索出现了 Felix Lazebnik的线性丢番图方程系统。在那篇论文中，作者考虑了以下问题：

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

于 2012-12-16T12:43:36.857 回答
2

当您转换为 时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])
于 2012-12-16T03:18:02.177 回答
1

我可能误解了你的问题，但我认为你只需要roundand then的组合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)
于 2012-12-16T03:20:01.617 回答
1

+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
于 2012-12-16T17:00:41.473 回答
1

我需要这样做并最终将 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 存储库的链接）。

于 2015-07-19T05:38:21.897 回答
1

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

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.]
于 2021-04-04T01:31:43.437 回答
0

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

于 2018-06-10T23:46:02.573 回答