Пример Wiki для итерации Арнольди работает только для реальных матриц? - PullRequest
Новая
       15

Пример Wiki для итерации Арнольди работает только для реальных матриц?

0 голосов
/ 29 октября 2018

Запись в Википедии для метода Арнольди предоставляет пример Python, который создает основу подпространства Крылова матрицы A. Предположительно, если A является эрмитовым (т.е. если A == A.conj().T), то матрица Гессенберга h, сгенерированная этим алгоритмом, является трехдиагональной ( source ). Однако, когда я использую код Википедии на реальной эрмитовой матрице, матрица Хессенберга вовсе не является трехдиагональной. Когда я выполняю вычисления на действительной части A (так что A == A.T), тогда я получаю трехдиагональную матрицу Гессенберга, поэтому, похоже, существует проблема с мнимыми компонентами 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()

Результат:

enter image description here

1 Ответ

0 голосов
/ 29 октября 2018

После просмотра слайдов презентации на конференции я понял, что в какой-то момент 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
...