Animating Linear Transformations with Quiver

scientist inevitably means working on multiple layers of abstraction, foremost abstractions of code and math. This is great, because this allows you to get astonishing results quickly. But sometimes it is well-advised to pause for a second and ponder what actually happens behind a neat interface. This process of pondering is often assisted by visualizations. In this article, I want to present how animated quiver plots can help to ponder about linear transformations, which often toil away reliably in the obscurity of machine learning algorithms and associated interfaces. In the end, we will be able to visualize concepts like Singular Value Decomposition with our quiver plot.

Plotting Static Quiver Plots

A quiver plot from the matplotlib python package allows us to plot arrows (which in our case represent vectors). Let’s first take a look at a static quiver plot:

We can directly derive the transformation matrix from the image by looking at the target positions of the two base vectors. The first base vector is starting at position (1, 0) and landing on (1, 1), while the second base vector travels from (0, 1) to (-1, 1). Therefore the matrix, that describes this transformation is:

\[
\begin{pmatrix}
1 & -1 \\
1 & 1 \\
\end{pmatrix}
\]

Visually this corresponds to an anti-clockwise rotation by 45 degrees (or \(\pi/4\) in radian) and a slight stretch (by the factor \(\sqrt{2}\)).

With this information, let’s look at how this is implemented with quiver (Note that I omit some boilerplate code like scaling of axis):

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

def quiver_plot_base_vectors(transformation_matrix: np.ndarray):
    # Define vectors
    basis_i = np.array([1, 0])
    basis_j = np.array([0, 1])
    i_transformed = transformation_matrix[:, 0]
    j_transformed = transformation_matrix[:, 1]
    
    # plot vectors with quiver-function
    cmap = cm.inferno
    fig, ax = plt.subplots()
    ax.quiver(0, 0, basis_i[0], basis_i[1], 
        color=cmap(0.2), 
        scale=1, 
        angles="xy", 
        scale_units="xy", 
        label="i", 
        alpha=0.3)
    ax.quiver(0, 0, i_transformed[0], i_transformed[1], 
        color=cmap(0.2), 
        scale=1,   
        angles="xy",  
        scale_units="xy", 
        label="i_transformed")
    ax.quiver(0, 0, basis_j[0], basis_j[1], 
        color=cmap(0.5), 
        scale=1, 
        angles="xy", 
        scale_units="xy", 
        label="j", 
        alpha=0.3)
    ax.quiver(0, 0, j_transformed[0], j_transformed[1], 
        color=cmap(0.5), 
        scale=1, 
        angles="xy", 
        scale_units="xy", 
        label="j_transformed")

if __name__ == "__main__":
    matrix = np.array([
        [1, -1],
        [1, 1]  
    ])
    quiver_plot_base_vectors(matrix)

As you can see we defined one quiver plot per vector. This is just for illustrative purposes. If we look at the signature of the quiver function – quiver([X, Y], U, V, [C], /, **kwargs) – we can observe that U and V take numpy arrays as input, which is better than providing scalar values. Let’s refactor this function to using only one quiver invocation. Additionally, let’s add a vector v = (1.5, -0.5) to see the transformation applied on it.

def quiver_plot(transformation_matrix: np.ndarray, vector: np.ndarray):
    # Define vectors
    basis_i = np.array([1, 0])
    basis_j = np.array([0, 1])
    i_transformed = transformation_matrix[:, 0]
    j_transformed = transformation_matrix[:, 1]
    vector_transformed = transformation_matrix @ vector
    U, V = np.stack(
        [
            basis_i, i_transformed,
            basis_j, j_transformed,
            vector, vector_transformed,
        ],
        axis=1)

    # Draw vectors
    color = np.array([.2, .2, .5, .5, .8, .8])
    alpha = np.array([.3, 1.0, .3, 1.0, .3, 1.0])
    cmap = cm.inferno
    fig, ax = plt.subplots()
    ax.quiver(np.zeros(6), np.zeros(6), U, V,
        color=cmap(color),
        alpha=alpha,
        scale=1,
        angles="xy",
        scale_units="xy",
    )

if __name__ == "__main__":
    matrix = np.sqrt(2) * np.array([
        [np.cos(np.pi / 4), np.cos(3 * np.pi / 4)],
        [np.sin(np.pi / 4), np.sin(3 * np.pi / 4)]
    ])
    vector = np.array([1.5, -0.5])
    quiver_plot(matrix, vector)

This is much shorter and convenient than the first example. What we did here was to stack each vector horizontally producing the following array:

The first row corresponds to the U-parameter of quiver and the second to V. While the columns hold our vectors, where \(\vec{i}\) is the first base vector, \(\vec{j}\) is the second and \(\vec{v}\) is our custom vector. The indices, b and a, stand for before and after (i.e. whether the linear transformation is applied or not). Let’s look at the output:

Taking a second look at the code it might be confusing what happened to our neat and simple transformation matrix, which was restated to:

\[
{\scriptsize
M=\begin{pmatrix}
{1}&{-1}\\
{1}&{1}\\
\end{pmatrix}={\sqrt{2}}
\begin{pmatrix}
{\cos\left(\frac{1}{4}\pi\right)}&{\cos\left(\frac{3}{4}\pi\right)}\\
{\sin\left(\frac{1}{4}\pi\right)}&{\sin\left(\frac{3}{4}\pi\right)}\\
\end{pmatrix}
}
\]

The reason is, as we move on by adding animations, this representation will come in handy. The scalar multiplication by the square root of two represents how much our vectors get stretched, while the elements of the matrix are rewritten in trigonometric notation to depict the rotation in the unit circle.

Let’s animate

Reasons to add animations may include cleaner plots as we can get rid of the ghost vectors and create a more engaging experience for presentations. In order to enhance our plot with animations we can stay in the matplotlib ecosystem by utilizing the FuncAnimation() function from matplotlib.animation. The function takes the following arguments:

• a matplotlib.figure.Figure object
• an update function
• the number of frames

For each frame the update function gets invoked producing an updated version of the initial quiver plot. For more details check the official documentation from matplotlib.

With this information in mind our task is to define the logic to implement in the update function. Let’s start simple with only three frames and our base vectors. On frame 0 we are in the initial state. While on the last frame (frame 2) we need to arrive at the restated matrix M. Therefore we would expect to be half way there on frame 1. Because the arguments of \(\cos\) and \(\sin\) in M show the radians (i.e. how far we have traveled on the unit circle), we can divide them by two in order to get our desired rotation. (The second vector gets a negative \(\cos\) , because we are now in the second quadrant). Similarly we need to account for the stretch, represented by the scalar factor. We do this by computing the change in magnitude, which is \(\sqrt{2}-1\), and adding half of that change to the initial scaling.

\[
{\scriptsize
\begin{aligned}
\text{Frame 0:} \quad &
\begin{pmatrix}
\cos(0) & \cos\left(\frac{\pi}{2}\right) \\
\sin(0) & \sin\left(\frac{\pi}{2}\right)
\end{pmatrix}
\\[1em]
\text{Frame 1:} \quad &
s \cdot \begin{pmatrix}
\cos\left(\frac{1}{2} \cdot \frac{\pi}{4}\right) & -\cos\left(\frac{1}{2} \cdot \frac{3\pi}{4}\right) \\
\sin\left(\frac{1}{2} \cdot \frac{\pi}{4}\right) & \sin\left(\frac{1}{2} \cdot \frac{3\pi}{4}\right)
\end{pmatrix}, \quad \text{with } s = 1 + \frac{\sqrt{2} – 1}{2}
\\[1em]
\text{Frame 2:} \quad &
\sqrt{2} \cdot \begin{pmatrix}
\cos\left(\frac{\pi}{4}\right) & \cos\left(\frac{3\pi}{4}\right) \\
\sin\left(\frac{\pi}{4}\right) & \sin\left(\frac{3\pi}{4}\right)
\end{pmatrix}
\end{aligned}
}
\]

One Caveat to the explanation above: It serves the purpose to give intuition to the implementation idea and holds true for the base vectors. However the actual implementation contains some more steps, e.g. some transformations with \(\arctan\) to get the desired behavior for all vectors in the two-dimensional space.

So let’s inspect the main parts of the implementation. The full code can be found on my github.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import cm

class AnimationPlotter:
[...]
def animate(self, filename='output/mat_transform.gif'):
        self.initialize_plot()
        anim = animation.FuncAnimation(
            self.fig,
            self.update_quiver,
            frames=self.frames + 1,
            init_func=self.init_quiver,
            blit=True,
        )
        anim.save(filename, writer='ffmpeg', fps=self.frames/2)
        plt.close()
   
if __name__ == "__main__":
    matrix = np.sqrt(2) * np.array([
        [np.cos(np.pi / 4), np.cos(3 * np.pi / 4)],
        [np.sin(np.pi / 4), np.sin(3 * np.pi / 4)]
       
    ])
    vector = np.array([1.5, -0.5]).reshape(2, 1)
    transformer = Transformer(matrix)
    animation_plotter = AnimationPlotter(transformer, vector)
    animation_plotter.animate()

The animate() method belongs to a custom class, which is called AnimationPlotter. It does what we already learned with the inputs as provided above. The second class on the scene is a custom class Transformer, which takes care of computing the linear transformations and intermediate vectors for each frame. The main logic lies within the AnimationPlotter.update_quiver() and Transformer.get_intermediate_vectors() methods, and looks as follows.

class AnimationPlotter:
    [...]
    def update_quiver(self, frame: int):
        incremented_vectors = self.transformer.get_intermediate_vectors(
            frame, self.frames
        )
        u = incremented_vectors[0]
        v = incremented_vectors[1]
        self.quiver_base.set_UVC(u, v)
        return self.quiver_base,

class Transformer:
    [...]
    def get_intermediate_vectors(self, frame: int, total_frames: int) -> np.ndarray:
         change_in_direction = self.transformed_directions - self.start_directions
         change_in_direction = np.arctan2(np.sin(change_in_direction), np.cos(change_in_direction))
         increment_direction = self.start_directions + change_in_direction * frame / total_frames
         increment_magnitude = self.start_magnitudes + (self.transformed_magnitudes - self.start_magnitudes) * frame / total_frames
         incremented_vectors = np.vstack([np.cos(increment_direction), np.sin(increment_direction)]) @ np.diag(increment_magnitude)
         return incremented_vectors

What happens here is that for each frame the intermediate vectors get computed. This is done by taking the difference between the end and start directions (which represent vector angles). The change in direction/angle is then normalized to the range \([-\pi, \pi]\) and added to the initial direction by a ratio. The ratio is determined by the current and total frames. The magnitude is determined as already described. Finally, the incremented vector gets computed based on the direction and magnitude and this is what we see at each frame in the animation. Increasing the frames to say 30 or 60 makes the animation smooth.

Animating Singular Value Decomposition (SVD)

Finally I want to showcase how the introductory animation was created. It shows how 4 vectors (each for every quadrant) get transformed consecutively three times. Indeed, the three transformations applied correspond to our meanwhile well-known transformation matrix M from above, but decomposed via Singular Value Decomposition (SVD). You can gain or refresh your knowledge about SVD in this great and intuitive tds article. Or take a look here if you prefer a more math-focused read. However, with numpy.linalg.svd() it is straightforward to compute the SVD of our matrix M. Doing so results in the following decomposition:

\[
{\scriptsize
\begin{align}
A \vec{v} &= U\Sigma V^T\vec{v} \\[1em]
\sqrt{2} \cdot \begin{pmatrix}
\cos\left(\frac{\pi}{4}\right) & \cos\left(\frac{3\pi}{4}\right) \\
\sin\left(\frac{\pi}{4}\right) & \sin\left(\frac{3\pi}{4}\right)
\end{pmatrix} \vec{v} &=
\begin{pmatrix}
\cos\left(\frac{3\pi}{4}\right) & \cos\left(\frac{3\pi}{4}\right) \\
\sin\left(\frac{-\pi}{4}\right) & \sin\left(\frac{\pi}{4}\right)
\end{pmatrix}
\begin{pmatrix}
\sqrt{2} & 0 \\
0 & \sqrt{2}
\end{pmatrix}
\begin{pmatrix}
-1 & 0 \\
0 & 1
\end{pmatrix} \vec{v}
\end{align}
}
\]

Note how the stretching by the square root gets distilled by the middle matrix. The following animation shows how this looks in action (or motion) for v = (1.5, -0.5).

In the end the purple vector \(\vec{v}\) arrives at its determined position in both cases.

Conclusion

To wrap it up, we can use quiver() to display vectors in 2D space and, with the help of matplotlib.animation.FuncAnimation(), add appealing animations on top. This results in clear visualizations of linear transformations that you can use, for instance, to demonstrate the underlying mechanics of your machine learning algorithms. Feel free to fork my repository and implement your own visualizations. I hope you enjoyed the read!