我无法从维基百科实现python Arnoldi方法来查找特征值。我做错了什么?

2024-06-01 13:26:19 发布

您现在位置:Python中文网/ 问答频道 /正文

祝大家新年快乐。我目前正在尝试理解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])

Tags: oftheinnumpyisnpbasisrandom