Plotting an ellipse with eigenvectors using matplotlib and numpy

Question

I ran into a deadlock with the following code, and I hope you can help me out. It seems like a deadlock situation because I need to understand the math in order to write the code, but I also need the code to understand the math. I have spent half a day trying to figure out what's wrong, but with no success. It could be possible that everything is correct, but who knows?

Let me first show you when things go well. My goal is to plot an ellipsoid (an ellipse) in a 2D space.

Theoretical background

Ellipsoid

For this, I take the standard definition of the ellipsoid which tells you that an ellipsoid is set of the form (x-z)^T * D * (x-z) <= 1 where D is a positive definite matrix, z the center of an ellipsoid and x is a point in the 2D space. For simplicity reasons we assume that z=0 hence an ellipsoid is just x^T*D*x <= 1 for x in R^2.

An ellipsoid is defined for any dimension. An ellipsoid is 2D is just an (filled) ellipse.

On the other hand, an ellipse in 2D can be described using the formula (x^2 / a^2) + (y^2 / b^2) = 1 where a and b are the lengths of the semi-axes of the ellipsoid. If the ellipsoid (or ellipse) is centered at the origin and oriented such that its principal axes align with the coordinate axes, then the semi-axes are the distances from the origin to the points where the ellipsoid (or ellipse) intersects the coordinate axes.

Eigenvalues and eigenvectors

Let the positive definite matrix defining the ellipsoid be as follows:

D = 1 0
    0 4

To find the eigenvalues of the matrix, we solve the equation det(D - λI) = 0 which yields λ1 = 1 and λ2 = 4. We then take the square roots of the reciprocals of these eigenvalues to obtain the lengths of the semi-axes of the ellipsoid. Specifically, a = sqrt(1/1) = 1 is the length of the first axis, and b = sqrt(1/4) = 1/2 is the length of the second axis.

To determine the eigenvectors of the matrix that define the directions in which the space is going to be stretched or compressed, we solve the systems D - λ1*I = 0 and D - λ2*I = 0. This gives us the eigenspace [u 0] for λ1 and the eigenspace [0 w] for λ2. We set u and w to 1 and conclude that [1 0] is an eigenvector for the eigenvalue 1, while [0 1] is an eigenvector for the eigenvalue 4.

Parameter Transformation and Scaling

Next, we will plot the ellipse or the border of the ellipsoid. An ellipse can be described by the following equation: (x^2 / a^2) + (y^2 / b^2) = 1, where a and b are the lengths of the semi-axes of the ellipsoid. To plot the ellipse using matplotlib, we use the parameter transformation. We set x = a * cos(t) and y = b * sin(t) for 0 <= t <= 2pi.

Finally, we will plot the eigenvectors that define the axes of the ellipsoid. We take the desired length of the semi-axis, multiply it by the corresponding eigenvector, and scale it by the norm of this vector.

Plot

Plotting everything together we get

If I'm correct, the plot shows two semi-axes with lengths 1 and 1/2 that align with the x and y-axis, respectively, in line with the eigenvalues. Additionally, we can shift or rotate the ellipsoid using the rotation matrix:

Implementation

We start by creating an array of points and get eigenvalues and eigenvectors from the built-in functions:

# Generate points on the ellipse.
theta = np.linspace(0, 2 * np.pi, 10000)
eigenvalues, eigenvectors = np.linalg.eig(positive_definite_matrix)

As far as I understand, np.linalg.eig() returns a 2D array as the second output. It is called "eigenvectors", but these are not eigenvectors, since the first array contains the x coordinates of the vectors and the second array contains the y coordinates of the vectors. Therefore, we need to transpose it:

# Transpose the result to get eigenvectors we can calculate with.
eigenvectors = eigenvectors.T

Now we can get the lengths of the semi-axes as described in the theory section:

a = np.sqrt(1/eigenvalues[0])
b = np.sqrt(1/eigenvalues[1])

We get the ellipse points:

ellipse_points = a * np.cos(theta)[:, np.newaxis] * eigenvectors[:, 0] + b * np.sin(theta)[:, np.newaxis] * eigenvectors[:, 1]

If an angle for the rotation was given, we can rotate the ellipsoid using the rotation matrix:

rotation_matrix = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)], [np.sin(rotation_angle), np.cos(rotation_angle)]])
rotated_points = np.dot(rotation_matrix, ellipse_points.T).T

We shift all these points to a new center:

rotated_points += center

We plot the ellipse and the center:

ax.plot(rotated_points[:, 0], rotated_points[:, 1], 'b-')
ax.scatter(center[0], center[1], c='b', s=100, marker='o', linewidths=1)

Finally, we need to plot the eigenvectors, which are scaled by the length of the corresponding semi-axes:

# Rotate eigenvectors
rotated_eigenvectors = np.dot(rotation_matrix, eigenvectors).T
# Scale the eigenvectors according to the axis
rotated_eigenvectors[0] = a * rotated_eigenvectors[0] / np.linalg.norm(rotated_eigenvectors[0])
rotated_eigenvectors[1] = b * rotated_eigenvectors[1] / np.linalg.norm(rotated_eigenvectors[1])
plot_vector(ax, center, rotated_eigenvectors[0] + center[0], **BASIS_VECTOR_PROPERTIES)
plot_vector(ax, center, rotated_eigenvectors[1] + center[1], **BASIS_VECTOR_PROPERTIES)

The entire code along with the plotting properties is shown below:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import collections as mc
from pylab import *

# Plotting area properties. Normally are in a separate module.
PLOTTING_BOUND = 2
PLOT_BOTTOM_LEFT_CORNER = np.array([-PLOTTING_BOUND, -PLOTTING_BOUND])
PLOT_TOP_RIGHT_CORNER = np.array([PLOTTING_BOUND, PLOTTING_BOUND])
BASIS_VECTOR_PROPERTIES = {
   'scale_units' : 'xy', 
   'angles' : 'xy',
   'scale' : 1,
   'color' : 'b'
}

def plot_ellipsoid_in_work(ax, center, positive_definite_matrix, rotation_angle, want_to_plot_eigenvectors):
    # Generate points on the ellipse.
    theta = np.linspace(0, 2 * np.pi, 10000)
    eigenvalues, eigenvectors = np.linalg.eig(positive_definite_matrix)
    
    # Transpose the result to get eigenvectors in one array.
    eigenvectors = eigenvectors.T
    print("Eigenvalues: " + str(eigenvalues))
    print("Eigenvectors: " + str(eigenvectors))
    a = np.sqrt(1/eigenvalues[0])
    b = np.sqrt(1/eigenvalues[1])
    print("Semi-axis: " + str(a) + ", " + str(b))
    
    # Get ellipse points.
    ellipse_points = a * np.cos(theta)[:, np.newaxis] * eigenvectors[:, 0] + b * np.sin(theta)[:, np.newaxis] * eigenvectors[:, 1]
    
    # Rotate the points.
    rotation_matrix = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)], [np.sin(rotation_angle), np.cos(rotation_angle)]])
    rotated_points = np.dot(rotation_matrix, ellipse_points.T).T
    
    # Shift by center.
    rotated_points += center
    
    # Plot the ellipse
    ax.plot(rotated_points[:, 0], rotated_points[:, 1], 'b-')
            
    ax.scatter(center[0], center[1], c='b', s=100, marker='o', linewidths=1)

    # Show eigenvectors
    if want_to_plot_eigenvectors:
        # Rotate eigenvectors
        rotated_eigenvectors = np.dot(rotation_matrix, eigenvectors).T
        # Scale the eigenvectors according to the axis
        rotated_eigenvectors[0] = a * rotated_eigenvectors[0] / np.linalg.norm(rotated_eigenvectors[0])
        rotated_eigenvectors[1] = b * rotated_eigenvectors[1] / np.linalg.norm(rotated_eigenvectors[1])
        plot_vector(ax, center, rotated_eigenvectors[0] + center[0], **BASIS_VECTOR_PROPERTIES)
        plot_vector(ax, center, rotated_eigenvectors[1] + center[1], **BASIS_VECTOR_PROPERTIES)

def plot_vector(ax: plt.Axes, start_point: np.ndarray, vector: np.ndarray, **properties) -> None:
    """
    Plot a vector given its starting point and displacement vector.

    Args:
        ax (matplotlib.axes.Axes): The matplotlib axes on which to plot the vector.
        start_point (np.ndarray): A 2D numpy array representing the starting point of the vector.
        vector (np.ndarray): A 2D numpy array representing the displacement vector.
        **properties: Additional properties to be passed to the quiver function.

    Returns:
        None
    """
    displacement = vector - start_point
    ax.quiver(start_point[0], start_point[1], displacement[0], displacement[1], **properties)
    
def setup_plot(ax: plt.Axes, ldown: np.ndarray, rup: np.ndarray) -> None:
    """
    Setup the given Matplotlib axis for plotting.

    Args:
        ax (matplotlib.axes.Axes): The Matplotlib axis to be setup.
        ldown (np.ndarray): A 2D numpy array representing the lower left corner of the rectangular region.
        rup (np.ndarray): A 2D numpy array representing the upper right corner of the rectangular region.

    Returns:
        None
    """
    ax.set_aspect('equal')
    ax.grid(True, which='both')
    ax.set_xlim(ldown[0], rup[0])
    ax.set_ylim(ldown[1], rup[1])

    plt.axhline(0, color='black')
    plt.axvline(0, color='black')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')

    mngr = plt.get_current_fig_manager()
    mngr.resize(960, 1080)

# Create plot
fig, ax = plt.subplots()
setup_plot(ax, PLOT_BOTTOM_LEFT_CORNER, PLOT_TOP_RIGHT_CORNER)

# Plot the ellipsoid
plot_ellipsoid_in_work(ax,np.array([1,1]),np.array([[1,0],[0,4]]),np.pi/4,True)

# Show the plot
plt.show()

Problem

UPD: center problem solved. The error was in the wrong indexing. Even though I went carefully through every step, I just do not understand what is happening if I vary the inputs. Putting the center at [-1,1], we get

I am just shifting everything by the center in the lines

rotated_points += center
...
plot_vector(ax, center, rotated_eigenvectors[0] + center[0], **BASIS_VECTOR_PROPERTIES)
plot_vector(ax, center, rotated_eigenvectors[1] + center[1], **BASIS_VECTOR_PROPERTIES)

Why it makes the eigenvectors to point to some strange direction?

UPD2: This problem is not solved. I add a plot and description at the end of the post in bold Putting D =

1 3
1 4

and calling

plot_ellipsoid_in_work(ax,np.array([0,0]),np.array([[1,3],[1,4]]),0,True)

results in

Here, I am wondering why the ellipsoid gets rotated, even though the angle is zero. I have picked some points and checked if they satisfy the ellipsoid equation, and they do. So, I assume that the ellipsoid is plotted correctly. However, I am curious why the eigenvectors do not align with the axis of the ellipsoid and are not scaled correctly. The left one is too short, and the center one is too long.

I calculate the first eigenvalue by hand and get λ1 = (5 + sqrt(21)) / 2, which is approximately 4.7912878474779195. I check the output in the terminal and see that this corresponds to one eigenvalue:

Eigenvalues: [0.20871215 4.79128785]

Hence, the eigenvalue is computed correctly. Then I take a = sqrt(1/λ1) = sqrt(2 / (5 + sqrt(21))). This must be the length of this axis.

Then I compute the eigenspace by hand and get a * [1,(3 + sqrt(21))/6]. I set a=1 and plot the scaled vector [1,(3 + sqrt(21))/6]:

a = np.sqrt(2 / (5 + np.sqrt(21)))
plot_vector(ax, np.array([0,0]), a * eigenvector / np.linalg.norm(eigenvector), color='r', scale=1, scale_units='xy')

Now this new vector points somewhere into eternity, and I am giving up.

I am struggling to understand why my code is not working as expected. I started by trying to understand how the matrix D of the ellipsoid defines the ellipsoid. I computed the eigenvalues and they seemed to be correct. However, when I computed the eigenvectors, they were different from those returned by the np.linalg.eig() function. I am now questioning my understanding of math.

I wrote the code, and it seems to work fine, with the eigenvectors always aligned with the ellipsoid axes. However, when I used the matrix [1 3; 1 4] as described before, the axes were not aligned anymore. I am questioning whether I have made a mistake somewhere or if there is a problem with the library functions.

I am frustrated and confused, and my question is essentially "What is wrong with this code?". I am not sure if my understanding of math, code, or the Python library is flawed.

I hope you can help me figure out what I am missing, because at this point, I am just lost.

UPD3: here is the problem in detail. If we use the same matrix D, we get the following output:

You can see the output of the program in blue. In red I added the manually calculated eigenvector and the axis of the ellipsoid. I consulted myself regarding the math part, and there must be no error there. However, we can see that both program and my manually calculated eigenvectors are out of the axes. Mathematically they must be on the same axis, and we can assume that both the blue and red eigenvectors are calculated correctly. And, mathematically we can assume that the eigenvectors must be on the ellipse axes. The question is, why do they get plotted with wrong length and angle. Is there an issue with the quiver function? What is the reason?

UPD4: Example of a symmetric matrix

12 3
3 15

for which the plot is still not aligned: