我为解决 Ax = b
而实现的最速下降显示出一些奇怪的行为:对于任何足够大的矩阵 (~10 x 10
,只测试了方阵,所以far),返回的 x
包含所有巨大的值(按 1x10^10
的顺序)。
def steepestDescent(A, b, numIter=100, x=None):
"""Solves Ax = b using steepest descent method"""
warnings.filterwarnings(action="error",category=RuntimeWarning)
# Reshape b in case it has shape (nL,)
b = b.reshape(len(b), 1)
exes = []
res = []
# Make a guess for x if none is provided
if x==None:
x = np.zeros((len(A[0]), 1))
exes.append(x)
for i in range(numIter):
# Re-calculate r(i) using r(i) = b - Ax(i) every five iterations
# to prevent roundoff error. Also calculates initial direction
# of steepest descent.
if (numIter % 5)==0:
r = b - np.dot(A, x)
# Otherwise use r(i+1) = r(i) - step * Ar(i)
else:
r = r - step * np.dot(A, r)
res.append(r)
# Calculate step size. Catching the runtime warning allows the function
# to stop and return before all iterations are completed. This is
# necessary because once the solution x has been found, r = 0, so the
# calculation below divides by 0, turning step into "nan", which then
# goes on to overwrite the correct answer in x with "nan"s
try:
step = np.dot(r.T, r) / np.dot( np.dot(r.T, A), r )
except RuntimeWarning:
warnings.resetwarnings()
return x
# Update x
x = x + step * r
exes.append(x)
warnings.resetwarnings()
return x, exes, res
(返回exes
和res
用于调试)
我认为问题一定出在计算 r
或 step
(或一些更深层次的问题)上,但我无法弄清楚它是什么。
最佳答案
代码似乎是正确的。例如,以下测试对我有效(linalg.solve 和 steepestDescent 大部分时间都给出了接近的答案):
import numpy as np
n = 100
A = np.random.random(size=(n,n)) + 10 * np.eye(n)
print(np.linalg.eig(A)[0])
b = np.random.random(size=(n,1))
x, xs, r = steepestDescent(A,b, numIter=50)
print(x - np.linalg.solve(A,b))
问题出在数学上。如果 A 是正定矩阵,则该算法保证收敛到正确的解。通过将 10 * 单位矩阵添加到随机矩阵,我们增加了所有特征值都为正的概率
如果您使用大型随机矩阵(例如 A = random.random(size=(n,n))
进行测试,您几乎可以肯定具有负特征值,并且算法不会收敛。
关于python - 最速下降吐出不合理的大值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/38594289/