<jupyter_start><jupyter_text>Arnoldi Iteration with Complex Eigenvalues<jupyter_code>import numpy as np
import numpy.linalg as la

import matplotlib.pyplot as pt<jupyter_output><empty_output><jupyter_text>Let us make a matrix with a defined set of eigenvalues and eigenvectors, given by `eigvals` and `eigvecs`.<jupyter_code>np.random.seed(40)

# Generate matrix with eigenvalues 1...25
nhalf = 12
n = 2*nhalf
eigvals = np.zeros(n,dtype=np.complex)
eigvals[0:nhalf] = np.linspace(1., nhalf, nhalf)+1j*np.linspace(1., nhalf, nhalf)
eigvals[nhalf:] = eigvals[0:nhalf].conj()
eigvecs = np.random.randn(n, n)
print(eigvals)

A = la.solve(eigvecs, np.dot(np.diag(eigvals), eigvecs))
print(la.eig(A)[0])<jupyter_output>[  1. +1.j   2. +2.j   3. +3.j   4. +4.j   5. +5.j   6. +6.j   7. +7.j
   8. +8.j   9. +9.j  10.+10.j  11.+11.j  12.+12.j   1. -1.j   2. -2.j
   3. -3.j   4. -4.j   5. -5.j   6. -6.j   7. -7.j   8. -8.j   9. -9.j
  10.-10.j  11.-11.j  12.-12.j]
[ 12.-12.j  12.+12.j  11.-11.j  11.+11.j  10.-10.j  10.+10.j   9. +9.j
   9. -9.j   8. -8.j   8. +8.j   7. +7.j   7. -7.j   6. +6.j   5. +5.j
   4. +4.j   3. +3.j   2. +2.j   1. +1.j   6. -6.j   5. -5.j   1. -1.j
   4. -4.j   3. -3.j   2. -2.j]<jupyter_text>Initialization Set up $Q$ and $H$:<jupyter_code>Q = np.zeros((n, n),dtype=np.complex)
H = np.zeros((n, n),dtype=np.complex)

k = 0<jupyter_output><empty_output><jupyter_text>Pick a starting vector, normalize it<jupyter_code>x0 = np.random.randn(n)#+1j*np.random.randn(n)
x0 = x0/la.norm(x0)

# Poke it into the first column of Q
Q[:, k] = x0

del x0<jupyter_output><empty_output><jupyter_text>Make a list to save arrays of Ritz values, including the maximum Ritz value:<jupyter_code>ritz_values = []
ritz_max = []<jupyter_output><empty_output><jupyter_text>Algorithm Carry out one iteration of Arnoldi iteration.Run this cell in-place (Ctrl-Enter) until H is filled.<jupyter_code>print(k)

u = A @ Q[:, k]

# Carry out Gram-Schmidt on u against Q
for j in range(k+1):
    qj = Q[:, j]
    H[j,k] = np.inner(qj.conj(), u)
    u = u - H[j,k]*qj

if k+1 < n:
    H[k+1, k] = la.norm(u)
    Q[:, k+1] = u/H[k+1, k]

k += 1

pt.spy(H)
D = la.eig(H)[0]
max_ritz = D[np.argmax(np.abs(D))]
print(max_ritz)
ritz_max.append(max_ritz)
ritz_values.append(D)<jupyter_output>23
(12+12j)<jupyter_text>Check that $Q^T A Q =H$:<jupyter_code>la.norm(Q[:,:k-1].T.conj() @ A @ Q[:,:k-1] - H[:k-1,:k-1])/ la.norm(A)<jupyter_output><empty_output><jupyter_text>Check that Q is orthogonal:<jupyter_code>la.norm((Q.T.conj() @ Q)[:k-1,:k-1] - np.eye(k-1))<jupyter_output><empty_output><jupyter_text>Look at convergence of largest Ritz value to dominant eigenvalue<jupyter_code>max_eigval = eigvals[np.argmax(np.abs(eigvals))]
pt.plot(np.abs(max_eigval-ritz_max)/np.abs(max_eigval))<jupyter_output><empty_output><jupyter_text>Plot convergence of Ritz values Enable the Ritz value collection above to make this work.<jupyter_code>for i, rv in enumerate(ritz_values[1:]):
    pt.plot([i] * len(rv), rv, "x")<jupyter_output>/home/edgar/environments/py35env/lib/python3.5/site-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
{'fraction_non_alphanumeric': 0.18063314711359404, 'fraction_numerical': 0.05276225946617008, 'mean_word_length': 3.6914119359534205, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '644833', 'license': 'MIT', 'n_tokens_mistral': 1558, 'n_tokens_neox': 1504, 'n_words': 316, 'path': 'eigenvalue/Arnoldi Iteration with Complex Eigenvalues.ipynb', 'repo_name': 'JiaheXu/MATH', 'size': '3222'}