祝大家新年快乐。我目前正在尝试理解Arnoldi算法。作为第一步,我尝试实现Wikipedia implementation of the Arnoldi Iteration。为了检查我的理解,我决定对照numpy linalg eig函数检查结果。为此,我计算计算出的H_n矩阵的特征值与eig给出的特征值之间的差值
不幸的是,只有一个特征值似乎是相同的。我对阿诺迪有什么误解?(注意:我知道代码在当前状态下效率很低。)
这是我的代码:
import numpy as np
from numpy.linalg import eig
def arnoldi_iteration(A, b, n: int):
"""Computes a basis of the (n + 1)-Krylov subspace of A: the space
spanned by {b, Ab, ..., A^n b}.
Arguments
A: m × m array
b: initial vector (length m)
n: dimension of Krylov subspace, must be >= 1
Returns
Q: m x (n + 1) array, the columns are an orthonormal basis of the
Krylov subspace.
h: (n + 1) x n array, A on basis Q. It is upper Hessenberg.
"""
m = A.shape[0]
h = np.zeros((n + 1, n))
Q = np.zeros((m, n + 1))
q = b / np.linalg.norm(b) # Normalize the input vector
Q[:, 0] = q # Use it as the first Krylov vector
for k in range(n):
v = A.dot(q) # Generate a new candidate vector
for j in range(k): # Subtract the projections on previous vectors
h[j, k] = np.dot(Q[:, j].conj(), v)
v = v - h[j, k] * Q[:, j]
h[k + 1, k] = np.linalg.norm(v)
eps = 1e-13 # If v is shorter than this threshold it is the zero vector
if h[k + 1, k] > eps: # Add the produced vector to the list, unless
q = v / h[k + 1, k] # the zero vector is produced.
Q[:, k + 1] = q
else: # If that happens, stop iterating.
return h
return h
n = 50 #specify number of dimensions
A = np.random.rand(n,n) #create Matrix with random values
b = np.random.rand(n) #create initial vector with random values
#print difference in eigenvalues in arnoldi and in numpy method
print(eig(arnoldi_iteration(A, b, n)[:-1])[0]-eig(A)[0])
目前没有回答
相关问题 更多 >
编程相关推荐