使用scipy.optimize.root解决稀疏的非线性方程组 [英] solving a sparse non linear system of equations using scipy.optimize.root

问题描述

我想解决以下非线性方程组.

I want to solve the following non-linear system of equations.

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

• the dot between a_k and x represents dot product.
• the 0 in the first equation represents 0 vector and 0 in the second equation is scaler 0
• all the matrices are sparse if that matters.

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

• K is an n x n (positive definite) matrix
• each A_k is a known (symmetric) matrix
• each a_k is a known n x 1 vector
• N is known (let's say N = 50). But I need a method where I can easily change N.

• x是n x 1一个向量.
• 1 <= k <= N缩放器的每个alpha_k

• x is an n x 1 a vector.
• each alpha_k for 1 <= k <= N a scaler

我正在考虑使用 scipy根查找x和每个alpha_k.本质上，第一个方程式的每一行都有n方程式，约束方程式还有另一个N方程式来求解我们的n + N变量.因此，我们拥有所需数量的方程式才具有解决方案.

I am thinking of using scipy root to find x and each alpha_k. We essentially have n equations from each row of the first equation and another N equations from the constraint equations to solve for our n + N variables. Therefore we have the required number of equations to have a solution.

对于x和alpha_k's我也有一个可靠的初步猜测.

I also have a reliable initial guess for x and the 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. 对于较大的矩阵，我可能必须显式地计算雅可比行列式?

1. I can actually brute force this toy problem and feed it to the solver. But how do I do I solve this toy problem in such a way that I can extend it easily to the case when I have let's say n=50 and N=50
2. I will probably have to explicitly compute the Jacobian for larger matrices??.

有人可以给我任何指示吗?

Can anyone give me any pointers?

推荐答案

我认为scipy.optimize.root方法固若金汤，但是对于这个方程组来说，避免琐碎的解决方案可能是真正的挑战.

I think the scipy.optimize.root approach holds water, but steering clear of the trivial solution might be the real challenge for this system of equations.

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

In any event, this function uses root to solve the system of equations.

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

但是，以玩具示例为例，只会产生平凡的解决方案.

Running on the toy example, however, only produces the trivial solution.

# 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解决稀疏的非线性方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！