问题描述
我想解决以下非线性方程组.
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
betweena_k
andx
representsdot product
. - the
0
in the first equation represents0 vector
and0
in the second equation isscaler 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 ann 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 ann x 1
a vector.- each
alpha_k
for1 <= 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
问题:
- 我实际上可以蛮力解决这个玩具问题,并将其输入求解器.但是我该如何解决这个玩具问题,以便可以轻松地将其扩展到我说
n=50
和N=50
的情况. - 对于较大的矩阵,我可能必须显式地计算雅可比行列式?
- 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
andN=50
- 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解决稀疏的非线性方程组的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!