# coding: utf-8 # # Power Iteration and its Variants # In[2]: import numpy as np import numpy.linalg as la # Let's prepare a matrix with some random or deliberately chosen eigenvalues: # In[6]: n = 6 if 1: np.random.seed(70) eigvecs = np.random.randn(n, n) eigvals = np.sort(np.random.randn(n)) # Uncomment for near-duplicate largest-magnitude eigenvalue # eigvals[1] = eigvals[0] + 1e-3 A = eigvecs.dot(np.diag(eigvals)).dot(la.inv(eigvecs)) print(eigvals) else: # Complex eigenvalues np.random.seed(40) A = np.random.randn(n, n) print(la.eig(A)[0]) # Let's also pick an initial vector: # In[4]: x0 = np.random.randn(n) x0 # ### Power iteration # In[5]: x = x0 # Now implement plain power iteration. # # Run the below cell in-place (Ctrl-Enter) many times. # In[6]: x = np.dot(A, x) x # * What's the problem with this method? # * Does anything useful come of this? # * How do we fix it? # ### Normalized power iteration # Back to the beginning: Reset to the initial vector. # In[7]: x = x0/la.norm(x0) # Implement normalized power iteration. # # Run this cell in-place (Ctrl-Enter) many times. # In[8]: x = np.dot(A, x) nrm = la.norm(x) x = x/nrm print(nrm) print(x) # * What do you observe about the norm? # * What about the sign? # * What is the vector \$x\$ now? # # Extensions: # # * Now try the matrix variants above. # * Suggest a better way of estimating the eigenvalue. [Hint](https://en.wikipedia.org/wiki/Rayleigh_quotient) # ------ # What if we want the smallest eigenvalue (by magnitude)? # # Once again, reset to the beginning. # In[9]: x = x0/la.norm(x0) # Run the cell below in-place many times. # In[10]: x = la.solve(A, x) nrm = la.norm(x) x = x/nrm print(1/nrm) print(x) # * What's the computational cost per iteration? # * Can we make this method search for a specific eigenvalue? # * What is this [method](https://en.wikipedia.org/wiki/Inverse_iteration) called? # -------------- # Can we feed an estimate of the current approximate eigenvalue back into the calculation? (Hint: Rayleigh quotient) # # Reset once more. # In[11]: x = x0/la.norm(x0) # Run this cell in-place (Ctrl-Enter) many times. # In[12]: sigma = np.dot(x, np.dot(A, x))/np.dot(x, x) x = la.solve(A-sigma*np.eye(n), x) x = x/la.norm(x) print(sigma) print(x) # * What's this [method](https://en.wikipedia.org/wiki/Rayleigh_quotient_iteration) called? # * What's a reasonable stopping criterion? # * Computational downside of this iteration?