A Simple Class Assignment Exemplifies Computational Linear Algebra

By Paul Davis

The John von Neumann Lecture is among SIAM’s oldest and most prestigious prizes, recognizing “outstanding and distinguished contributions to the field of applied mathematics.” In that sense, the lecture is similar to a mathematical version of highlights from the Winter Olympics — a collage of victorious moments over mathematical challenges as varied as slalom runs, half-pipes, ski jumps, and ice rinks.

Charles Van Loan’s prize lecture at the 2018 SIAM Annual Meeting, which took place this July in Portland, Ore., had an entirely different feel. The Cornell University professor’s audience was treated to a view from Van Loan’s helmet camera as he navigated down an expert slope, beginning with a deceptively simple problem statement and culminating in a rich set of discoveries and conjectures about untangling random polygons. In that single run, he exhibited the insights, skills, and accomplishments that comprise the recognition in his prize citation: “pioneering contributions to research in numerical linear algebra and to the exposition of the subject…a brilliant communicator…(whose) book Matrix Computations…has shaped and influenced the way we think about matrix computations."

The crux of Van Loan’s exemplary lecture originated 10 years ago in a burst of inspiration on the eve of a first-year undergraduate class called “Introduction to Computing.” He prepared the following assignment: Display a sequence of polygons where each polygon is obtained from its predecessor by connecting the midpoints of its sides. Let the original polygon be random. Figure 1 illustrates the averaging process for five vertices and reveals that every polygon shares the same centroid.

After seeing these polygons shrink around their common centroid, Van Loan added an additional requirement: that the new arrays of averaged \(x\) and \(y\) coordinates be renormalized to unit length: \(\textrm{x=x/norm(x)}\), \(\textrm{y=y/norm(y)}\). He also fixed the centroid at the origin: \(\textrm{x=x-mean(x)}\), \(\textrm{y=y-mean(y)}\).

The normalized version of the problem produces a limiting polygon inscribed in an ellipse tilted at a 45-degree angle (see Figure 2). The vertices appear to crawl around the ellipse as the iteration proceeds. Van Loan asked his students to provide explanations of these observations, and posed questions drawn from a computational mathematician’s repertoire of reflexive queries: How long does it take the iterations to converge? Will they always converge? What is the inverse of this process?

The original just-in-time problem formulation debuted in a 2008 freshman-level MATLAB course at Cornell and resurfaced in its normalized form in a 2009 upper-level course on matrix computations. Van Loan and then-undergraduate Adam Elmachtoub answered those aforementioned questions and posed others in a 2010 SIAM Review article [1]. This entire adventure in computational linear algebra came full circle as the subject of this year’s John von Neumann Lecture.

When speaking to his SIAM audience, Van Loan sought to “explain it all using matrix computations.” The problem he posed to his first-year students in 2008 was part of a strategy to help beginning computer scientists “see \(n\) things as one” and get them “thinking at the vector level before they headed over to the math department,” he said. Van Loan aimed to stimulate a hunger for mathematics in his students that could explain their observations, opening the door to one of the core experiences of applied mathematics. At a basic computational level, he led them to recognize matrix-vector operations lurking implicitly within algorithms.

Establishing the matrix-vector formulation of the vertex update was the first important turn in Van Loan’s run. The second was recognizing that the iterative process of averaging polygons amounts to the power method for the update matrix \(M_n\), whose five-vertex version appears in Figure 3; the \(k\)th \(x\)-vector is the product \(M_n^k {\cdot} (\textrm{initial}\) \(x\)-\(\textrm{vector})\), and the \(k\)th \(y\)-vector is obtained similarly.

An analysis of the eigenvalues and eigenvectors of \(M_n\) now lay ahead; the linear algebra was picking up speed. Van Loan noted that one can write the update matrix as \(M_n=(I_n+S_n)/2\), where \(S_n\) is the upshift matrix, i.e., the identity matrix with its last column shifted to the first position and all other columns displaced one spot to the right. Its eigenvalues and eigenvectors are completely known. With the exception of the details—a Schur decomposition—needed to actually cross the finish line, this race was over.

As an example, Figure 4 shows the structure of the Schur decomposition of \(M_5\). To follow the effects of powers of the update matrix on an array of \(x\) or \(y\) coordinates of an initial polygon, one can expand those arrays in terms of a real Schur basis, as in Figure 5.

Because the polygon’s centroid is at the origin, each coordinate array averages to zero. Hence, the principle (black) eigenvector drops out of the orthonormal expansion in Figure 5. Applying powers of \(M_5\) to a vertex coordinate array amounts to applying those powers to the red and blue basis vectors. Since both pairs of complex conjugate eigenvalues have magnitude less than one, the successive averaged polygons collapse to the origin.

If each iteration renormalizes the vertex coordinate arrays to unit vectors, only the red component survives because the complex conjugate pair of red eigenvalues has larger magnitude. At that point, we might as well have started with coordinate arrays that are norm-one linear combinations of the red vectors in Figure 5’s Schur basis expansion. With just a touch more analysis via a singular-value decomposition, the limiting ellipse emerges — tilted 45 degrees from the coordinate axes.

There is still more before the finish. The vertices appear to move around the limiting ellipse because they alternate between even and odd steps. Iterating in reverse yields starburst polygons instead of the smoother, ellipsoidal polygons of forward iteration; the invariant subspace associated with the complex eigenvalues of smallest magnitude dominates, and the components of its basis vectors tend to flip-flop in sign. Computational evidence also points to a conjectural version of Kepler’s second law: the centered “pizza slice” triangles lying on the limiting ellipse have equal area.

Having quickly negotiated this trail through the highpoints of computational linear algebra, Van Loan crossed the finish line with a graceful verbal embrace rather than the familiar arms-raised victory salute. He was explicit in his appreciation of the opportunity to speak and the presence of so many students in the audience. Recalling his first SIAM meeting—at Newport News, Va., in 1973—Van Loan spoke with the conviction of personal experience. “SIAM does [student involvement] better than any other professional society,” he said.

Van Loan’s lecture is available from SIAM either as slides with synchronized audio or a PDF of slides only.

Since 1960, SIAM has annually recognized a John von Neumann lecturer “for outstanding and distinguished contributions to the field of applied mathematics, and for the effective communication of these ideas to the community.” The award honors John von Neumann (1903-1957), one of the most prolific and articulate practitioners of applied mathematics in the 20th century. It includes an honorarium of $5,000 from a fund established by IBM and maintained by SIAM.

The figures in this article were provided by Charles Van Loan.

References
[1] Elmachtoub, A.N., & Van Loan, C.F. (2010). From Random Polygon to Ellipse: An Eigenanalysis. SIAM Rev., 52(1), 151-170.

Paul Davis is professor emeritus of mathematical sciences at Worcester Polytechnic Institute.