Arnoldi 迭代的 Wiki 示例仅适用于实矩阵？

Question

Wikipedia entry for the Arnoldi method provides a Python example that produces basis of the Krylov subspace of a matrix A. Supposedly, if A is Hermitian (i.e. if A == A.conj().T) then the Hessenberg matrix h generated by this algorithm is tridiagonal (source)。但是，当我在现实世界的 Hermitian 矩阵上使用 Wikipedia 代码时，Hessenberg 矩阵根本不是三对角矩阵。当我对 A 的实部进行计算时（因此 A == A.T），我确实得到了一个三对角 Hessenberg 矩阵，因此 A 的虚部似乎有问题.有人知道为什么维基百科代码没有产生预期的结果吗？

工作示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import circulant


def arnoldi_iteration(A, b, n):
    m = A.shape[0]

    h = np.zeros((n + 1, n), dtype=np.complex)
    Q = np.zeros((m, n + 1), dtype=np.complex)

    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 + 1):  # Subtract the projections on previous vectors
            h[j, k] = np.dot(Q[:, j], v)
            v = v - h[j, k] * Q[:, j]

        h[k + 1, k] = np.linalg.norm(v)
        eps = 1e-12  # 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 Q, h
    return Q, h


# Construct matrix A
N = 2**4
I = np.eye(N)
k = np.fft.fftfreq(N, 1.0 / N) + 0.5
alpha = np.linspace(0.1, 1.0, N)*2e2
c = np.fft.fft(alpha) / N
C = circulant(c)
A = np.einsum("i, ij, j->ij", k, C, k)

# Show that A is Hermitian
print(np.allclose(A, A.conj().T))

# Arbitrary (random) initial vector
np.random.seed(0)
v = np.random.rand(N)
# Perform Arnoldi iteration with complex A
_, h = arnoldi_iteration(A, v, N)
# Perform Arnoldi iteration with real A
_, h2 = arnoldi_iteration(np.real(A), v, N)

# Plot results
plt.subplot(121)
plt.imshow(np.abs(h))
plt.title("Complex A")
plt.subplot(122)
plt.imshow(np.abs(h2))
plt.title("Real A")
plt.tight_layout()
plt.show()

结果：

Answer 1

浏览了一些会议演示幻灯片后，我意识到在某些时候 Q 必须在 A 是复数的情况下共轭。正确的算法贴在下面供参考，并标出代码更改（注意，此更正也已提交到维基百科词条）：

import numpy as np

def arnoldi_iteration(A, b, n):
    m = A.shape[0]

    h = np.zeros((n + 1, n), dtype=np.complex)
    Q = np.zeros((m, n + 1), dtype=np.complex)

    q = b / np.linalg.norm(b)
    Q[:, 0] = q

    for k in range(n):
        v = A.dot(q)
        for j in range(k + 1):
            h[j, k] = np.dot(Q[:, j].conj(), v)  # <-- Q needs conjugation!
            v = v - h[j, k] * Q[:, j]

        h[k + 1, k] = np.linalg.norm(v)
        eps = 1e-12
        if h[k + 1, k] > eps:
            q = v / h[k + 1, k]
            Q[:, k + 1] = q
        else:
            return Q, h
    return Q, h

Arnoldi 迭代的 Wiki 示例仅适用于实矩阵？

Wiki example for Arnoldi iteration only works for real matrices?

python

matrix

wikipedia

complex-numbers