使用 scipy.optimize.root 求解稀疏非线性方程组

Question

我想求解以下非线性方程组。

备注

a_k和x之间的dot表示dot product.
第一个方程中的0表示0 vector，第二个方程中的0表示scaler 0
如果重要的话，所有矩阵都是稀疏的。

已知

K 是一个 n x n（正定）矩阵
每个 A_k 是一个已知的（对称）矩阵
每个 a_k 是一个已知的 n x 1 向量
N 是已知的（假设 N = 50）。但是我需要一种可以轻松更改 N 的方法。

未知（试图求解）

x 是一个 n x 1 向量。
每个 alpha_k 用于 1 <= k <= N 一个定标器

我的想法。

我正在考虑使用 scipy root 来查找 x 和每个 alpha_k。我们基本上有来自第一个方程的每一行的 n 方程和来自约束方程的另一个 N 方程来求解我们的 n + N 变量。因此，我们有所需数量的方程来求解。

我对 x 和 alpha_k's 也有一个可靠的初始猜测。

玩具示例。

n = 4
N = 2
K = np.matrix([[0.5, 0, 0, 0], [0, 1, 0, 0],[0,0,1,0], [0,0,0,0.5]])
A_1 = np.matrix([[0.98,0,0.46,0.80],[0,0,0.56,0],[0.93,0.82,0,0.27],[0,0,0,0.23]])
A_2 = np.matrix([[0.23, 0,0,0],[0.03,0.01,0,0],[0,0.32,0,0],[0.62,0,0,0.45]])
a_1 = np.matrix(scipy.rand(4,1))
a_2 = np.matrix(scipy.rand(4,1))

我们正在努力解决

 x = [x1, x2, x3, x4] and alpha_1, alpha_2

问题：

我实际上可以暴力破解这个玩具问题并将其提供给求解器。但是我该如何解决这个玩具问题，以便我可以轻松地将它扩展到我说 n=50 和 N=50
我可能必须显式计算更大矩阵的雅可比矩阵？？。

谁能给我指点一下？

Answer 1

我认为 scipy.optimize.root 方法站得住脚，但避开平凡的解决方案可能是这个方程组的真正挑战。

无论如何，此函数使用root求解方程组。

def solver(x0, alpha0, K, A, a):
'''
x0     - nx1 numpy array. Initial guess on x.
alpha0 - nx1 numpy array. Initial guess on alpha.
K      - nxn numpy.array.
A      - Length N List of nxn numpy.arrays.
a      - Length N list of nx1 numpy.arrays.
'''

# Establish the function that produces the rhs of the system of equations.
n = K.shape[0]
N = len(A)
def lhs(x_alpha):
    '''
    x_alpha is a concatenation of x and alpha.
    '''

    x = np.ravel(x_alpha[:n])
    alpha = np.ravel(x_alpha[n:])
    lhs_top = np.ravel(K.dot(x))
    for k in xrange(N):
        lhs_top += alpha[k]*(np.ravel(np.dot(A[k], x)) + np.ravel(a[k]))

    lhs_bottom = [0.5*x.dot(np.ravel(A[k].dot(x))) + np.ravel(a[k]).dot(x)
                  for k in xrange(N)]

    lhs = np.array(lhs_top.tolist() + lhs_bottom)

    return lhs

# Solve the system of equations.
x0.shape = (n, 1)
alpha0.shape = (N, 1)
x_alpha_0 = np.vstack((x0, alpha0))
sol = root(lhs, x_alpha_0)
x_alpha_root = sol['x']

# Compute norm of residual.
res = sol['fun']
res_norm = np.linalg.norm(res)

# Break out the x and alpha components.
x_root = x_alpha_root[:n]
alpha_root = x_alpha_root[n:]


return x_root, alpha_root, res_norm

运行然而，在玩具示例中，只产生了平凡的解决方案。

# Toy example.
n = 4
N = 2
K = np.matrix([[0.5, 0, 0, 0], [0, 1, 0, 0],[0,0,1,0], [0,0,0,0.5]])
A_1 = np.matrix([[0.98,0,0.46,0.80],[0,0,0.56,0],[0.93,0.82,0,0.27],      
                [0,0,0,0.23]])
A_2 = np.matrix([[0.23, 0,0,0],[0.03,0.01,0,0],[0,0.32,0,0],
      [0.62,0,0,0.45]])
a_1 = np.matrix(scipy.rand(4,1))
a_2 = np.matrix(scipy.rand(4,1))
A = [A_1, A_2]
a = [a_1, a_2]
x0 = scipy.rand(n, 1)
alpha0 = scipy.rand(N, 1)

print 'x0 =', x0
print 'alpha0 =', alpha0

x_root, alpha_root, res_norm = solver(x0, alpha0, K, A, a)

print 'x_root =', x_root
print 'alpha_root =', alpha_root
print 'res_norm =', res_norm

输出是

x0 = [[ 0.00764503]
 [ 0.08058471]
 [ 0.88300129]
 [ 0.85299622]]
alpha0 = [[ 0.67872815]
 [ 0.69693346]]
x_root = [  9.88131292e-324  -4.94065646e-324   0.00000000e+000        
          0.00000000e+000]
alpha_root = [ -4.94065646e-324   0.00000000e+000]
res_norm = 0.0

使用 scipy.optimize.root 求解稀疏非线性方程组

solving a sparse non linear system of equations using scipy.optimize.root

numpy

scipy

sparse-matrix

mathematical-optimization

equation-solving

备注

已知

未知（试图求解）

我的想法。

玩具示例。

问题：