<jupyter_start><jupyter_text>Семинар 5. Степенной метод. QR алгоритм  Экспериментально посмотрим на сходимость степенного метода$$x_{k+1} = \frac{Ax_k}{\|Ax_k\|} $$- Теретическая скорость сходимости **линейная** с фактором $\frac{|\lambda_2|}{|\lambda_1|}$<jupyter_code>import numpy as np

def power_method(A, x0, max_iter, eps):
    x = x0.copy()
    conv = [x]
    for i in range(max_iter):
        x = A @ x
        x = x / np.linalg.norm(x)
        conv.append(x)
        eigval = x @ (A @ x)
        res = A @ x - eigval * x
        if np.linalg.norm(res) < eps:
            break
    return x, eigval, conv<jupyter_output><empty_output><jupyter_text>Симметричная положительно определённая матрица<jupyter_code>np.random.seed(0)
n = 10
A = np.random.randn(n, n)
A = A.T @ A

true_eigvals, true_eigvec = np.linalg.eigh(A)
print(true_eigvals)
print(true_eigvec)
x0 = np.random.randn(n)
max_eigvec, max_eigval, conv = power_method(A, x0, 50, 1e-6)
print(max_eigvec)
print(max_eigval)
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)

power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")
power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, -1] * true_eigvec[:, -1]) for x in conv])
plt.plot(power_method_eigvec_conv, linewidth=4, label="Eigvec conv")
plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)
print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])
print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])
print(true_eigvals[-2] / true_eigvals[-1])<jupyter_output><empty_output><jupyter_text>Симметричная матрица- Покажем, что собственные значений вещественные    - $Ax = \lambda x$    - Умножим слева на $x^*$: $x^* A x = \lambda \|x\|^2_2$    - Возьмём сопряжение от обеих частей: $x^* A^* x = \lambda^* \|x\|_2^2$    - В силу эрмитовости: $x^* A x = \lambda^* \|x\|_2^2$    - Сравнивая со вторым равенством, получаем $\lambda^* = \lambda$    - Значит собственные значения действительные<jupyter_code>A = np.random.randn(n, n)
A = A + A.T
true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)
print(true_eigvec)
x0 = np.random.randn(n)
max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")

idx = np.argmax(np.abs(true_eigvals))

power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])
plt.plot(power_method_eigvec_conv, linewidth=4, label="Eigvec conv")
plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)
print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])
print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])
sorted_abs_eigvals = np.sort(np.abs(true_eigvals))
print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])<jupyter_output><empty_output><jupyter_text>Несимметричная матрица<jupyter_code>A = np.random.randn(n, n)

true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)
print(np.sort(np.abs(true_eigvals)))
print(true_eigvec)
x0 = np.random.randn(n)
max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, label="Eigval conv")

idx = np.argmax(np.abs(true_eigvals))

power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])
plt.plot(power_method_eigvec_conv, label="Eigvec conv")
plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)
print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])
print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])
sorted_abs_eigvals = np.sort(np.abs(true_eigvals))
print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])
print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])
print(power_method_eigvec_conv[1:] / power_method_eigvec_conv[:-1])
sorted_abs_eigvals = np.sort(np.abs(true_eigvals))
print(sorted_abs_eigvals[-2] / sorted_abs_eigvals[-1])<jupyter_output><empty_output><jupyter_text>Если $\lambda^*_2 = \lambda_1$, то степенной метод осциллирует и не сходится!   Другие проблемные случаи1) $\lambda_1 = \lambda_2$2) $\lambda_1 = -\lambda_2$  Если $\lambda_1 = \lambda_2$<jupyter_code>A = np.random.randn(n, n)
Q, _ = np.linalg.qr(A)
A = Q @ np.diagflat([2, 2] + list(np.random.rand(n - 2))) @ Q.T
true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)
print(np.abs(true_eigvals))
print(true_eigvec)
x0 = np.random.randn(n)
max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")

idx = np.argmax(np.abs(true_eigvals))

power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])
plt.plot(power_method_eigvec_conv, linewidth=4, label="Eigvec conv")
plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)
alpha = np.linalg.lstsq(true_eigvec[:, n-2:], max_eigvec, rcond=None)[0]
print("Accuracy = {}".format(np.linalg.norm(true_eigvec[:, n-2:] @ alpha - max_eigvec)))<jupyter_output><empty_output><jupyter_text>Если $\lambda_1 = -\lambda_2$<jupyter_code>A = np.random.randn(n, n)
Q, _ = np.linalg.qr(A)
A = Q @ np.diagflat([2, -2] + list(np.random.rand(n - 2))) @ Q.T
true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)
print(np.abs(true_eigvals))
print(true_eigvec)
x0 = np.random.randn(n)
max_eigvec, max_eigval, conv = power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - x @ (A @ x) * x) for x in conv])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")

idx = np.argmax(np.abs(true_eigvals))

power_method_eigvec_conv = np.array([np.linalg.norm(x - x[0] / true_eigvec[0, idx] * true_eigvec[:, idx]) for x in conv])
plt.plot(power_method_eigvec_conv, linewidth=4, label="Eigvec conv")
plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Блочный степенной метод- Ищем $k$ cтарших собственных значения и вектора- Этапы аналогичны степенному методу для максимального по модулю собственного значения    - Умножение матрицы на набор векторов    - Аналог нормализации для матриц?<jupyter_code>def block_power_method(A, x0, max_iter, eps):
    x = x0.copy()
    conv_vec = [x]
    conv_val = [np.einsum('ik,ij,jk->k',x, A, x)]
    for i in range(max_iter):
        x = A @ x
        x, _ = np.linalg.qr(x)
        conv_vec.append(x)
        eigval = np.einsum('ik,ij,jk->k', x, A, x)
        conv_val.append(eigval)
        res = A @ x - (eigval[:, np.newaxis] * x.T).T
        if np.linalg.norm(res) < eps:
            break
    return x, eigval, conv_vec, conv_val<jupyter_output><empty_output><jupyter_text>Немного про einsum<jupyter_code>A = np.random.randn(n, n)
print(np.trace(A), np.einsum('ii', A))
print(np.diag(A) - np.einsum("ii->i", A))
x = np.random.randn(n)
print(A @ x - np.einsum('ij,j', A, x))
B = np.random.randn(n, n)
print(np.diag(A @ B) - np.einsum('ij,ji->i', A, B))<jupyter_output><empty_output><jupyter_text>Более подробное об этой операции см. [тут](https://obilaniu6266h16.wordpress.com/2016/02/04/einstein-summation-in-numpy/)  Вернёмся к блочному степенному методу<jupyter_code>block_size = 3
A = np.random.randn(n, n)
A = A.T @ A
true_eigvals, true_eigvec = np.linalg.eigh(A)
print(true_eigvals)
print(true_eigvec)
x0 = np.random.randn(n, block_size)
max_eigvec, max_eigval, conv_vec, conv_val = block_power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
print(true_eigvec[:, n-block_size:n][:, ::-1])
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - lam * x) for x, lam in zip(conv_vec, conv_val)])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")

idx = np.argsort(np.abs(true_eigvals))

plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)
print(power_method_eigval_conv[1:] / power_method_eigval_conv[:-1])
print(true_eigvals[-3] / true_eigvals[-2])
print(true_eigvals[-1], true_eigvals[-2], true_eigvals[-3])<jupyter_output><empty_output><jupyter_text>Несимметричный случай<jupyter_code>A = np.random.randn(n, n)
true_eigvals, true_eigvec = np.linalg.eig(A)
block_size = 5
print(true_eigvals)
print(np.sort(np.abs(true_eigvals)))
print(true_eigvec.round(5))
x0 = np.random.randn(n, block_size)
max_eigvec, max_eigval, conv_vec, conv_val = block_power_method(A, x0, 500, 1e-6)
print(max_eigvec)
print(max_eigval)
print(true_eigvec[:, n-block_size:n][:, ::-1])
power_method_eigval_conv = np.array([np.linalg.norm(A @ x - lam * x) for x, lam in zip(conv_vec, conv_val)])
plt.figure(figsize=(10, 6))
plt.plot(power_method_eigval_conv, linewidth=4, label="Eigval conv")

idx = np.argsort(np.abs(true_eigvals))

plt.legend(fontsize=20)
plt.yscale("log")
plt.ylabel("Residual norm", fontsize=20)
plt.yticks(fontsize=16)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=16)
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Какие собственные значения удаётся восстановить, а какие нет?  QR алгоритм- Приведение матрицы к форме Шура с помощью унитарных преобразований$$ A = UTU^*, $$$T$ – верхнетреугольная матрица.<jupyter_code>def qr_algorithm(A, num_iter, eps):
    T = A.copy()
    U = np.eye(A.shape[0])
    conv = [(T, U)]
    for i in range(num_iter):
        Q, R = np.linalg.qr(T)
        T = R @ Q
        U = U @ Q
        conv.append((T, U))
        if np.sum(np.abs(np.tril(T, k=-1))) < eps:
            break
    return T, U, conv[1:]
n = 7
A = np.random.randn(n, n)
# A = A.T @ A
A = A + A.T
true_eigvals, true_eigvec = np.linalg.eig(A)
print(true_eigvals)
T, U, conv = qr_algorithm(A, 2000, 1e-6)
print(np.linalg.norm(A - U @ T @ U.T))
plt.spy(T, markersize=10, precision=1e-6)
plt.figure(figsize=(10, 7))
conv_qr = np.array([np.sum(np.abs(np.tril(T, k=-1))) for T, U in conv])
plt.plot(conv_qr)
plt.yscale("log")
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("Sum of lower triangular part in modulos", fontsize=20)
plt.grid(True)
U.round(4), true_eigvec.round(4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1651727438865199, 'fraction_numerical': 0.024105821947611173, 'mean_word_length': 5.481669486745629, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '172991', 'license': 'MIT', 'n_tokens_mistral': 5371, 'n_tokens_neox': 5242, 'n_words': 864, 'path': 'seminars/seminar5/seminar5.ipynb', 'repo_name': 'anastasiya-zhiboedova/nla2019_ozon', 'size': '11491'}