如何求解非线性方程组

Question

我有一组可以表示为矩阵的方程，我还要求对于解中的所有变量 x_i，(x_i)³ - x_i = 0.例如，

A = [0 1 0 0]
    [0 0 1 0]
    [1 0 0 1]

我也有，不是所有的变量 = 0。

这意味着

x₂ = 0
x₃ = 0
x₁ + x₄ = 0
(x₁)³ - x₁ = 0
(x₂)³ - x₂ = 0
(x₃)³ - x₃ = 0
(x₄)³ - x₄ = 0

一个简单的解决方案是 x₁ = 1 和 x₄ = -1。

您如何求解此类方程组的小实例？最好我想要一个至少可以从 python 调用的解决方案。

我目前解决问题的方法是尝试所有 3ⁿ 个不同的向量，其值从 -1,0,1 开始。

for v in itertools.product([-1,0,1], repeat = n):
    vector = np.asarray(v)
    if (not np.dot(M,v).any()):
        print "Solution found!"
        break

编辑

这应该是对@alko 答案的评论，但它太长了。让我通过一个例子来完成这个方法。

A = np.matrix([[0, 1, 1, 1, 1, 0, 1], [1, 0, 1, 1, 1, 1, 0], [0, 1, 0, 1, 1, 1, 1], [1, 0, 1, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1, 1]])
p,l,u=scipy.linalg.lu(A)
print u
[[ 1.  0.  1.  1.  1.  1.  0.]
 [ 0.  1.  1.  1.  1.  0.  1.]
 [ 0.  0. -1.  0.  0.  1.  0.]
 [ 0.  0.  0. -1.  0.  0.  1.]
 [ 0.  0.  0.  0. -1.  0.  0.]]

我不清楚下一步会是什么......？

Answer 1

既然你放了 SymPy 标签，我会指出 SymPy 可以很容易地通过符号方式解决这个问题

In [6]: x1, x2, x3, x4 = symbols('x1:5')

In [7]: solve([x2, x3, x1 + x4, x1**3 - x1, x2**3 - x2, x3**3 - x3, x4**3 - x4], [x1, x2, x3, x4], dict=True)
Out[7]: [{x₁: -1, x₂: 0, x₃: 0, x₄: 1}, {x₁: 0, x₂: 0, x₃: 0, x₄: 0}, {x₁: 1, x₂: 0, x₃: 0, x₄: -1}]

如果您的解决方案不是整数，事情会变得更加棘手，因为根式的解决方案可能不存在，或者至少 SymPy 可能无法找到它们。如果是这种情况，并且您只对数字解决方案感兴趣，您应该坚持使用像 numpy 或 scipy 这样的数字库，但由于您的解决方案都是 -1、0 或 1，因此这不是问题。

编辑：

如果你有一个矩阵，说：

A = Matrix([[0, 1, 0, 0],
            [0, 0, 1, 0],
            [1, 0, 0, 1]])

然后将其转换为系统很容易。只需将其乘以包含符号的向量（为方便起见，我已在此处切换为基于 0 的索引）：

In [13]: syms = symbols("x:4")

In [14]: s = Matrix(syms)

In [15]: constraints = [xi**3 - xi for xi in syms]

In [16]: A*s
Out[16]:
⎡  x₁   ⎤
⎢       ⎥
⎢  x₂   ⎥
⎢       ⎥
⎣x₀ + x₃⎦

In [17]: solve(list(A*s) + constraints, syms, dict=True)
Out[17]: [{x₀: -1, x₁: 0, x₂: 0, x₃: 1}, {x₀: 0, x₁: 0, x₂: 0, x₃: 0}, {x₀: 1, x₁: 0, x₂: 0, x₃: -1}]

这是您更大系统的解决方案：

In [35]: A = np.matrix([[0, 1, 1, 1, 1, 0, 1], [1, 0, 1, 1, 1, 1, 0], [0, 1, 0, 1, 1, 1, 1], [1, 0, 1, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1, 1]])

In [36]: A = Matrix(A).applyfunc(int)

In [37]: syms = symbols("x:7")

In [38]: s = Matrix(syms)

In [39]: constraints = [xi**3 - xi for xi in syms]

In [40]: solve(list(A*s) + constraints, syms, dict=True)
Out[40]:
[{x₀: -1, x₁: 1, x₂: 1, x₃: -1, x₄: 0, x₅: 1, x₆: -1}, {x₀: 0, x₁: 0, x₂: 0, x₃: 0, x₄: 0, x₅: 0, x₆: 0}, {x₀: 1, x₁: -1, x₂: -1, x₃: 1, x₄: 0, x₅: -1, x₆:
1}]

两个音符：

• 您不需要在 LU 中获取它。 SymPy 的解决方案会为您解决这个问题（如果您的计算时间过长，您可能会开始担心这些事情）。

• In [36] 将矩阵条目转换为整数（默认情况下它们是浮点数）。这不是必需的，但总的来说，当你知道它们是精确的时，SymPy 会更好地处理精确的数字，特别是因为你知道你的解决方案无论如何都是整数。如果您从一开始就使用 SymPy Matrix，则无需担心这一点。

Answer 2

这种问题称为约束规划。有一些python库可以解决这个问题。例如下面的代码使用or-tools

from constraint_solver import pywrapcp as cs
import numpy as np

A = np.array(
    [[0, 1, 0, 0],
     [0, 0, 1, 0],
     [1, 0, 0, 1]], np.bool)

#A = np.array( [[0, 1, 1, 1, 1, 0, 1],
#               [1, 0, 1, 1, 1, 1, 0], 
#               [0, 1, 0, 1, 1, 1, 1], 
#               [1, 0, 1, 0, 1, 1, 1], 
#               [0, 1, 0, 1, 0, 1, 1]], np.bool)

values = [-1, 0, 1]
solver = cs.Solver("name")
X = np.array([solver.IntVar(values) for i in range(A.shape[1])])
Xl = X.tolist()

for row in A:
    x = X[row].tolist()
    solver.Add(solver.Sum(x) == 0)

db = solver.Phase(Xl, 
                  solver.INT_VAR_DEFAULT, 
                  solver.INT_VALUE_DEFAULT)

solver.NewSearch(db)
while solver.NextSolution():
    solution = [x.Value() for x in Xl]
    print solution

输出：

[-1, 0, 0, 1]
[0, 0, 0, 0]
[1, 0, 0, -1]

Answer 3

你可以将A分解为lu decomposition：

import numpy as np
import scipy.linalg
A = np.matrix([[0, 1, 1, 1, 1, 0, 1],
               [1, 0, 1, 1, 1, 1, 0], 
               [0, 1, 0, 1, 1, 1, 1], 
               [1, 0, 1, 0, 1, 1, 1], 
               [0, 1, 0, 1, 0, 1, 1]])
p, l, u = scipy.linalg.lu(A)

为简单起见，让我们假设 A 有最大排名，即它的一些适当的未成年人是可逆的。如果不是这样，您可以使用与下面相同的 LU 分解将矩阵 A 减少为更小的 A'，其中 A' 满足此属性并且在方程方面是等效的（u 的最后一行将为零，删除它们并继续）。你的方程Ax=0等于(plu)x=0，并且p和l是可逆矩阵ux=0。 u，依次为上三角

print u
# [[ 1.  0.  1.  1.  1.  1.  0.]
#  [ 0.  1.  1.  1.  1.  0.  1.]
#  [ 0.  0. -1.  0.  0.  1.  0.]
#  [ 0.  0.  0. -1.  0.  0.  1.]
#  [ 0.  0.  0.  0. -1.  0.  0.]]

您只对最后（两）列感兴趣，它根据最后一个变量值定义所有变量。探测所有 x4 可能的非零值，1 和 -1，您可以检查您的方程组是否有解：

from itertools import product
idx = u.shape[0] - u.shape[1] 
m, v = u[:, :idx], u[:,idx:]
for spare in product({1, 0, -1}, repeat=-idx):
    if any(spare): # to exclude all zeros
        sol = scipy.linalg.solve(m, -np.dot(v, spare))
        if not set(sol).difference({0,1,-1}):
            print 'solution: ', sol, spare
# solution:  [-1.  1.  1. -1. -0.] (1, -1)
# solution:  [ 1. -1. -1.  1. -0.] (-1, 1)