Visualizing the Normalized Power Iteration and Inverse Iteration in Python 3

In this problem, we seek to show how the Power Iteration and Inverse Iteration with a shift of 0 act on the norm balls in R² that we saw in the norms homework. Before beginning the homework, recall what both methods do to a random vector. What happens to the vector? Does what happens depend on certain properties of the matrix? Once you think you have a reasonable answer to these questions, begin.

Instructions

In this homework we will visualize two numerical methods. One is the Power Iteration and the other is the Inverse Power Iteration. Your job is to create four functions, power_iteration inverse_iteration, visualize_power_iteration and visualize_inverse_iteration.

For power_iteration the function must do the following: (Do the same for inverse_iteration except use the Inverse Power Iteration with a shift of 0.)

First generate the norm ball with respect to the p norm using the given function get_normball(p).

Run the Power Iteration with the matrix A on all generated vectors in the norm ball for 4 iterations. Be sure to normalize the vectors using the p-norm. (How might you be able to np.linalg.norm to normalize a set of vectors?) Save the set of vectors at every step, starting with the original norm ball. Once finished you will have a list of 5 numpy arrays which represent different stages of the Power Iteration.

Return this list of numpy arrays, which represent a sequence of sets of points.

For visualize_power_iteration the function must do the following: (Do the same for visualize_inverse_iteration except use the Iverse Power Iteration.)

Calculate the eigenvectors of A . Normalize these using the specified p p norm then rescale them by their respective eigenvalue. Feel free to use numpy.linalg to compute the eigenvectors and eigenvalues.

Run your power_iteration function and use the provided plot_seq function to plot the sequence of points and the rescaled eigenvectors.

Return a list of rescaled eigenvectors. Make sure that the eigenvector with the largest eigenvalue is the first element in the list.

NOTE: It may be useful to verify for yourself that get_normball(p) will produce the ball with respect to thep norm by plotting the results of the function using plot_seq or by plotting it yourself. The demo function implements some of these functionalities.

Think About It

Run your newly made functions to see how Power Iteration and Inverse Iteration act on the 1,2, balls. Then respond to the following questions. Write your response as a string and assign it to the variable thoughts. This part is only graded for completion NOT for correctness.

Is there a relationship between your iteration results and the eigenvectors of A? How does the rate of convergence of the vectors change with the direction the vectors are pointing in? Will every vector always converge? Why or why not?

INPUT:

A: The matrix functions will be tested with.

OUTPUT:

power_iteration: Function performing Normalized Power Iteration the p norm ball.

inverse_iteration: Function performing Normalized Inverse Iteration with a shift of 0 on the p norm ball.

visualize_power_iteration: Function plotting the sequence of the set of vectors during Power Iteration.

visualize_inverse_iteration: Function plotting the sequence of the set of vectors during Inverse Iteration.

thoughts: Answers to the questions posed in the 'Think About It' section. Should be a string

Problem set-up code (click to view)

This is the code that is used to set up the problem and produce test data for your submission. It is reproduced here for your information, or if you would like to run your submission outside of . You should not copy/paste this code into the code box below. This code is run automatically 'behind the scenes' before your submitted code.

import numpy as np import numpy.linalg as la import matplotlib.pyplot as plt from scipy import linalg A = np.array([ [1.1, 0.3], [1.1, -0.1]]) def get_normball(p): """ :param p: The norm to use. :type p: Real number greater than .5. :returns r: A (2, 50) numpy array, representing a set of vectors in R^2 :type r: np.ndarray """ if p < .5: print("'p' must be >= .5") return t = np.linspace(0, 2*np.pi, 50) r = np.array([np.cos(t), -np.sin(t)]) r = r/la.norm(r, p, axis=0) return r def plot_seq(pts, vecs=None): """ :param pts: list of numpy array objects. Where each column is a point in R^2 The first row being the x coordinate and the second row being the y coordinate :type pts: list of np.ndarray :param vecs: A list of two 2-D numpy arrays to be plotted as arrows If this is not given, then the pts are only plotted :param vecs: list of np.ndarray :returns: None """ N = len(pts) if N < 0: print("'pts' is too short.") if N > 20: print("'pts' is too long.") return if vecs is not None: if type(vecs) != list: print("error in plot_seq(pts,vecs), 'vecs' should be a list of numpy arrays") return if len(vecs) != 2: print("'vecs' should be of length 2") return if type(pts[0]) != np.ndarray: print("'pts' should be a list of numpy arrays") return fig = plt.figure(figsize=(N*2,2)) plt.title("Power/Inverse Iteration in Normed spaces") for i in range(N): R = pts[i] ax = fig.add_subplot("1"+str(N)+str(i+1)) ax.set_aspect("equal") ax.set_ylim([-1.1, 1.1]) ax.set_xlim([-1.1, 1.1]) ax.grid() ax.set_ylabel("y") ax.set_xlabel("x") ax.scatter(R[0, :], R[1, :], c="b") if vecs: e1, e2 = vecs ax.arrow(0,0, e1[0], e1[1]) ax.arrow(0,0, e2[0], e2[1]) plt.show() def not_allowed(*args, **kwargs): raise RuntimeError("Calling this function is not allowed") la.inv = not_allowed linalg.inv = not_allowed 

This is the skeletal structure of what your code should look like. Fill in definitions with the missing functions.

import numpy as np import numpy.linalg as la def power_iteration(A, p): """ :param A: The matrix to apply Power Iteration to :type A: np.ndarray :param p: The specified norm to use for normalizing the vectors :type p: A real number greater than or equal to 1 :returns: A list of numpy arrays. In each array, each column is a point in R^2. The first row is the x coordinates and the second row is the y coordinates :type ret: list of np.ndarray """ pass def inverse_iteration(A, p): """ :param A: The matrix to apply Inverse Iteration to :type A: np.ndarray :param p: The specified norm to use for normalizing the vectors :type p: A real number greater than or equal to 1 :returns: A list of numpy arrays. In each array, each column is a point in R^2. The first row is the x coordinates and the second row is the y coordinates :type ret: list of np.ndarrays """ pass def visualize_power_iteration(A, p): """ :param A: The matrix to apply Power Iteration to :type A: np.ndarray :param p: The specified norm to use for normalizing the vectors :type p: A real number greater than or equal to 1 :returns: A list of two 2-D numpy arrays representing the two rescaled eigenvectors of A. The first element will be the eigenvector with the larger eigenvector. type ret: list of np.ndarrays """ pass def visualize_inverse_iteration(A, p): """ :param A: The matrix to apply Inverse Iteration to :type A: np.ndarray :param p: The specified norm to use for normalizing the vectors :type p: A real number greater than or equal to 1 :returns: A list of two 2-D numpy arrays representing the two rescaled eigenvectors of A. The first element will be the eigenvector with the larger eigenvector. type ret: list of np.ndarrays """ pass # This function demonstrates how to use the provided functions def demo(): rs = [] for p in [1, 2, np.inf]: R = get_normball(p) rs.append(R) # Two vectors that will be plotted by plot_seq vecs = [np.array([1,1]), np.array([-1,1])] # plots a sequence of a set of points. # plots the vectors in vecs as arrows plot_seq(rs, vecs) demo() visualize_power_iteration(A, 2) visualize_inverse_iteration(A, 2)