使用 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 变量。因此，我们有所需数量的方程来获得解。

我对@9​​87654345@ 和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

问题：

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

任何人都可以给我任何指示吗？

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