add exercise to eigen and fix cross references

Author: HumphreyYang

lectures/eigen.md renamed to lectures/eigen_I.md

@@ -51,6 +51,10 @@ We will use the following imports:
 ```{code-cell} ipython3
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib.lines import Line2D
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib.patches import FancyArrowPatch
+from mpl_toolkits.mplot3d import proj3d
 ```
 
 (matrices_as_transformation)=
@@ -587,6 +591,8 @@ $$
 Let's first see examples of a sequence of iterates $(A^k v)_{k \geq 0}$ under
 different maps $A$.
 
+(plot_series)=
+
 ```{code-cell} ipython3
 from numpy.linalg import matrix_power
 from matplotlib import cm
@@ -638,7 +644,7 @@ def plot_series(B, v, n):
 B = np.array([[sqrt(3) + 1, -2],
               [1, sqrt(3) - 1]])
 B = (1/(2*sqrt(2))) * B
-v = (-3,-3)
+v = (-3, -3)
 n = 12
 
 plot_series(B, v, n)
@@ -652,7 +658,7 @@ In this case, repeatedly multiplying a vector by $A$ makes the vector "spiral in
 B = np.array([[sqrt(3) + 1, -2],
               [1, sqrt(3) - 1]])
 B = (1/2) * B
-v = (2.5,0)
+v = (2.5, 0)
 n = 12
 
 plot_series(B, v, n)
@@ -667,7 +673,7 @@ an ellipse".
 B = np.array([[sqrt(3) + 1, -2],
               [1, sqrt(3) - 1]])
 B = (1/sqrt(2)) * B
-v = (-1,-0.25)
+v = (-1, -0.25)
 n = 6
 
 plot_series(B, v, n)
@@ -834,4 +840,337 @@ to one.
 
 The eigenvectors and eigenvalues of a map $A$ determine how a vector $v$ is transformed when we repeatedly multiply by $A$.
 
-This is discussed further later.
+This is discussed further later.
+
+
+```{exercise}
+:label: eig1_ex1
+
+Power iteration is a method for finding the largest absolute eigenvalue of a diagnalizable matrix.
+
+The method starts with a random vector $b_0$ and repeatedly applies the matrix $A$ to it
+
+$$
+b_{k+1}=\frac{A b_k}{\left\|A b_k\right\|}
+$$
+
+A thorough discussion of the method can be found [here](https://pythonnumericalmethods.berkeley.edu/notebooks/chapter15.02-The-Power-Method.html).
+
+In this exercise, implement the power iteration method and use it to find the largest eigenvalue of the matrix.
+
+Visualize your results by plotting the eigenvalue as a function of the number of iterations.
+```
+
+```{solution-start} eig1_ex1
+:class: dropdown
+```
+
+```{code-cell} ipython3
+# Define a matrix A
+A = np.array([[1, 0, 3], 
+              [0, 2, 0], 
+              [3, 0, 1]])
+
+# Define a number of iterations
+num_iters = 20
+
+# Define a random starting vector b
+b = np.random.rand(A.shape[1])
+
+# Initialize a list to store the eigenvector approximations
+res = []
+
+# Power iteration loop
+for i in range(num_iters):
+    # Multiply b by A
+    b = A @ b
+    # Normalize b
+    b = b / np.linalg.norm(b)
+    # Append b to the list of eigenvector approximations
+    dis = np.linalg.norm(np.array(b) - np.linalg.eig(A)[1][:, 0])
+    res.append(dis)
+    
+# Plot the eigenvector approximations for each iteration
+plt.figure(figsize=(10, 6))
+plt.xlabel('Iterations')
+plt.ylabel('L2 Norm')
+plt.title('Power Iteration and Eigenvector Approximations')
+_ = plt.plot(res)
+```
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+# Define a matrix A
+A = np.array([[1, 0, 3], 
+              [0, 2, 0], 
+              [3, 0, 1]])
+
+# Define a number of iterations
+num_iters = 20
+
+# Define a random starting vector b
+b = np.array([0.5, 1, 0.5])
+
+# Initialize a list to store the eigenvector approximations
+res = [b]
+
+# Power iteration loop
+for i in range(num_iters):
+    # Multiply b by A
+    b = A @ b
+    # Normalize b
+    b = b / np.linalg.norm(b)
+    # Append b to the list of eigenvector approximations
+    res.append(b)
+
+# Get the actual eigenvectors of matrix A
+eigenvector = np.linalg.eig(A)[1][:, 0]
+# Set up the figure and axis for 3D plot
+fig = plt.figure()
+ax = fig.add_subplot(111, projection='3d')
+
+# Plot the actual eigenvectors
+
+ax.scatter(eigenvector[0], eigenvector[1], eigenvector[2], color='r', s = 80)
+
+# Plot the approximated eigenvectors (b) at each iteration
+for i, vec in enumerate(res):
+    ax.scatter(vec[0], vec[1], vec[2], color='b', alpha=(i + 1) / (num_iters+1), s = 80)
+
+ax.set_xlabel('X')
+ax.set_ylabel('Y')
+ax.set_zlabel('Z')
+ax.set_title('Power iteration eigenvector approximations')
+points = [plt.Line2D([0], [0], linestyle='none', c=i, marker='o') for i in ['r', 'b']]
+ax.legend(points, ['Actual eigenvectors', 'Approximated eigenvectors (b)'], numpoints=1)
+ax.set_box_aspect(aspect=None, zoom=0.8)
+
+# Show the plot
+plt.show()
+```
+
+```{solution-end}
+```
+
+```{exercise}
+:label: eig1_ex2
+
+We have discussed the trajectory of the vector $v$ after being transformed by $A$.
+
+Consider the matrix $A = \begin{bmatrix} 1 & 2 \\ 1 & 1 \end{bmatrix}$ and the vector $v = \begin{bmatrix} 2 \\ -2 \end{bmatrix}$.
+
+Try to compute the trajectory of $v$ after being transformed by $A$ for $n=6$ iterations and plot the result.
+
+```
+
+
+```{solution-start} eig1_ex2
+:class: dropdown
+```
+
+```{code-cell} ipython3
+A = np.array([[1, 2], 
+              [1, 1]])
+v = (0.4, -0.4)
+n = 11
+
+# Compute right eigenvectors and eigenvalues
+eigenvalues, eigenvectors = np.linalg.eig(A)
+
+print(f"eigenvalues:\n {eigenvalues}")
+print(f"eigenvectors:\n {eigenvectors}")
+
+plot_series(A, v, n)
+```
+
+We find the trajectory of the vector $v$ after being transformed by $A$ for $n=6$ iterations and plot the result seems to converge to the eigenvector of $A$ with the largest eigenvalue.
+
+Let's use a vector field to visualize the transformation brought by A.
+
+```{code-cell} ipython3
+# Create a grid of points (vector field)
+x, y = np.meshgrid(np.linspace(-5, 5, 15), np.linspace(-5, 5, 20))
+
+# Apply the matrix A to each point in the vector field
+vec_field = np.stack([x, y])
+u, v = np.tensordot(A, vec_field, axes=1)
+
+
+# Plot the transformed vector field
+c = plt.streamplot(x, y, u - x, v - y, density=1, linewidth=None, color='#A23BEC')
+c.lines.set_alpha(0.5)
+c.arrows.set_alpha(0.5)
+
+# Plot the eigenvectors as long blue and green arrows
+origin = np.zeros((2, len(eigenvectors)))
+plt.quiver(*origin, eigenvectors[0], eigenvectors[1], color=['b', 'g'], angles='xy', scale_units='xy', scale=0.1, width=0.01)
+plt.quiver(*origin, - eigenvectors[0], - eigenvectors[1], color=['b', 'g'], angles='xy', scale_units='xy', scale=0.1, width=0.01)
+
+colors = ['b', 'g']
+lines = [Line2D([0], [0], color=c, linewidth=3) for c in colors]
+labels = ["2.4 eigenspace", "0.4 eigenspace"]
+plt.legend(lines, labels,loc='center left', bbox_to_anchor=(1, 0.5))
+
+plt.title("Convergence/Divergence towards Eigenvectors")
+plt.xlabel("x-axis")
+plt.ylabel("y-axis")
+plt.grid()
+plt.gca().set_aspect('equal', adjustable='box')
+plt.show()
+```
+
+Note that the vector field converges to the eigenvector of $A$ with the largest eigenvalue and diverges from the eigenvector of $A$ with the smallest eigenvalue.
+
+In fact, the eigenvectors are also the directions in which the matrix $A$ stretches or shrinks the space.
+
+Specifically, the eigenvector with the largest eigenvalue is the direction in which the matrix $A$ stretches the space the most.
+
+We will see more intriguing examples of eigenvectors in the following exercise.
+
+```{solution-end}
+```
+
+```{exercise}
+:label: eig1_ex3
+
+{ref}`Previously <plot_series>`, we demonstrated the trajectory of the vector $v$ after being transformed by $A$ for three different matrices.
+
+Use the visualization in the previous exercise to explain why the trajectory of the vector $v$ after being transformed by $A$ for the three different matrices.
+
+```
+
+
+```{solution-start} eig1_ex3
+:class: dropdown
+```
+
+Here is one solution
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+
+figure, ax = plt.subplots(1,3, figsize = (15,5))
+A = np.array([[sqrt(3) + 1, -2],
+              [1, sqrt(3) - 1]])
+A = (1/(2*sqrt(2))) * A
+
+B = np.array([[sqrt(3) + 1, -2],
+              [1, sqrt(3) - 1]])
+B = (1/2) * B
+
+C = np.array([[sqrt(3) + 1, -2],
+              [1, sqrt(3) - 1]])
+C = (1/sqrt(2)) * C
+
+examples = [A, B, C]
+
+for i, example in enumerate(examples):
+    M = example
+
+    # Compute right eigenvectors and eigenvalues
+    eigenvalues, eigenvectors = np.linalg.eig(M)
+    print(f'Example {i+1}:\n')
+    print(f'eigenvalues:\n {eigenvalues}')
+    print(f'eigenvectors:\n {eigenvectors}\n')
+
+    eigenvalues_real = eigenvalues.real
+    eigenvectors_real = eigenvectors.real
+
+    # Create a grid of points (vector field)
+    x, y = np.meshgrid(np.linspace(-20, 20, 15), np.linspace(-20, 20, 20))
+
+    # Apply the matrix A to each point in the vector field
+    vec_field = np.stack([x, y])
+    u, v = np.tensordot(M, vec_field, axes=1)
+
+    # Plot the transformed vector field
+    c = ax[i].streamplot(x, y, u - x, v - y, density=1, linewidth=None, color='#A23BEC')
+    c.lines.set_alpha(0.5)
+    c.arrows.set_alpha(0.5)
+
+    # Plot the eigenvectors as long blue and green arrows
+    parameters = {'color':['b', 'g'], 'angles':'xy', 'scale_units':'xy', 'scale':1, 'width':0.01, 'alpha':0.5}
+    origin = np.zeros((2, len(eigenvectors)))
+    ax[i].quiver(*origin, eigenvectors_real[0], eigenvectors_real[1], **parameters)
+    ax[i].quiver(*origin, - eigenvectors_real[0], - eigenvectors_real[1], **parameters)
+
+    ax[i].set_xlabel("x-axis")
+    ax[i].set_ylabel("y-axis")
+    ax[i].grid()
+    ax[i].set_aspect('equal', adjustable='box')
+
+plt.show()
+```
+
+The pattern demonstrated here is because we have complex eigenvalues and eigenvectors.
+
+It is important to acknowledge that there is a complex plane.
+
+If we add the complex axis for the plot, the plot will be more complicated.
+
+Here we used the real part of the eigenvalues and eigenvectors.
+
+We can try to plot the complex plane for one of the matrix using `Arrow3D` class retrieved from [stackoverflow](https://stackoverflow.com/questions/22867620/putting-arrowheads-on-vectors-in-matplotlibs-3d-plot).
+
+```{code-cell} ipython3
+class Arrow3D(FancyArrowPatch):
+    def __init__(self, xs, ys, zs, *args, **kwargs):
+        super().__init__((0,0), (0,0), *args, **kwargs)
+        self._verts3d = xs, ys, zs
+
+    def do_3d_projection(self, renderer=None):
+        xs3d, ys3d, zs3d = self._verts3d
+        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M)
+        self.set_positions((0.1*xs[0],0.1*ys[0]),(0.1*xs[1],0.1*ys[1]))
+
+        return np.min(zs)
+        
+
+# Define matrix A with complex eigenvalues
+A = np.array([[sqrt(3) + 1, -2],
+              [1, sqrt(3) - 1]])
+A = (1/(2*sqrt(2))) * A
+
+# Find eigenvalues and eigenvectors
+eigenvalues, eigenvectors = np.linalg.eig(A)
+
+# Create meshgrid for vector field
+x, y = np.meshgrid(np.linspace(-2, 2, 10), np.linspace(-2, 2, 10))
+
+# Calculate vector field (real and imaginary parts)
+u_real = A[0][0] * x + A[0][1] * y
+v_real = A[1][0] * x + A[1][1] * y
+u_imag = np.zeros_like(x)
+v_imag = np.zeros_like(y)
+
+# Create 3D figure
+fig = plt.figure()
+ax = fig.add_subplot(111, projection='3d')
+vlength = np.linalg.norm(eigenvectors)
+ax.quiver(x, y, u_imag, u_real-x, v_real-y, v_imag-u_imag, colors = 'b', alpha=0.3, length = .2, arrow_length_ratio = 0.01)
+
+arrow_prop_dict = dict(mutation_scale=2, arrowstyle='-|>', shrinkA=0, shrinkB=0)
+
+# Plot 3D eigenvectors
+for c, i in zip(['b', 'g'], [0, 1]):
+    a = Arrow3D([0, eigenvectors[0][i].real], [0, eigenvectors[1][i].real], 
+            [0, eigenvectors[1][i].imag], color=c, **arrow_prop_dict)
+    ax.add_artist(a)
+
+# Set axis labels and title
+ax.set_xlabel('X')
+ax.set_ylabel('Y')
+ax.set_zlabel('Im')
+ax.set_box_aspect(aspect=None, zoom=0.8)
+
+plt.draw()
+plt.show()
+```
+
+```{solution-end}
+```