Recently my fellow students and I have been spending our free time using Professor B------'s remarkable calculating engine to experiment with cellular automata, being mathematical contrivances that might be thought of as crude models of the lives of those most humble of creatures; amoebas. In their simplest form they are unending lines of boxes, some of which contain a living cell that at each generation will live, die or reproduce according to the contents of its neighbouring boxes. For example, we might say that each cell divides and its two offspring migrate to the left and right, dying if they encounter another cell's progeny.

Over the last few months we have seen how we can efficiently implement the Householder transformations and shifted Givens rotations used by Francis's algorithm to diagonalise a real symmetric matrix M, yielding its eigensystem in a matrix V whose columns are its eigenvectors and a diagonal matrix Λ whose diagonal elements are their associated eigenvalues, which satisfy

    M = V × Λ × V^T

and together are known as the spectral decomposition of M.
In this post, we shall add it to the ak library using the householder and givens functions that we have put so much effort into optimising.

Greetings Sir R-----. I trust that I find you in good spirits this evening? Will you take a glass of this excellent porter and join me in a little sport?

Splendid!

I propose a game that is popular amongst Antipodean opal scavengers as a means to improve their skill at guesswork.
Opals, as any reputable botanist will confirm, are the seeds of the majestic opal tree which grows in some abundance atop the vast monoliths of that region. Its mouth-watering fruits are greatly enjoyed by the Titans on those occasions when, attracted by its entirely confused seasons, they choose to winter thereabouts.

We have recently been looking at how we can use a special case of Francis's QR transformation to reduce a real symmetric matrix M to a diagonal matrix Λ by first applying Householder transformations to put it in tridiagonal form and then using shifted Givens rotations to zero out the off diagonal elements.
The columns of the matrix of transformations V and the elements on the leading diagonal of Λ are the unit eigenvectors and eigenvalues of M respectively and they consequently satisfy

    M × V = V × Λ

and, since the product of V and its transpose is the identity matrix

    M = V × Λ × V^T

which is known as the spectral decomposition of M.
Last time we saw how we could efficiently apply the Householder transformations in-place, replacing the elements of M with those of the matrix of accumulated transformations Q and creating a pair of arrays to represent the leading and off diagonal elements of the tridiagonal matrix. This time we shall see how we can similarly improve the implementation of the Givens rotations.

The Baron's most recent wager with Sir R----- set him the challenge of being the last to remove a horizontally, vertically or diagonally adjacent pair of draughts from a five by five square of them, with the Baron first taking a single draught and Sir R----- and he thereafter taking turns to remove such pairs.

When I heard these rules I was reminded of the game of Cram and could see that, just like it, the key to figuring the outcome is to recognise that the Baron could always have kept the remaining draughts in a state of symmetry, thereby ensuring that however Sir R----- had chosen he shall subsequently have been free to make a symmetrically opposing choice.

Over the last few months we have seen how we can use a sequence of Householder transformations followed by a sequence of shifted Givens rotations to efficiently find the spectral decomposition of a symmetric real matrix M, formed from a matrix V and a diagonal matrix Λ satisfying

    M × V = V × Λ

implying that the columns of V are the unit eigenvectors of M and their associated elements on the diagonal of Λ are their eigenvalues so that

    V × V^T = I

where I is the identity matrix, and therefore

    M = V × Λ × V^T

From a mathematical perspective the combination of Householder transformations and shifted Givens rotations is particularly appealing, converging on the spectral decomposition after relatively few matrix multiplications, but from an implementation perspective using ak.matrix multiplication operations is less than satisfactory since it wastefully creates new ak.matrix objects at each step and so in this post we shall start to see how we can do better.

Over the course of the year, my fellow students and I have been experimenting with an ethereal orrery which models the motion of heavenly bodies using nought but Sir N-----'s laws of gravitation and motion. Whilst the consequences of those laws are not generally subject to solution by mathematical reckoning, we were able to approximate them with a scheme that admitted errors of the order of the sixth power of the steps in time by which we advanced the positions of those bodies.
We have thus far employed it to model the solar system itself, uniformly distributed bodies of matter and the accretion of bodies that are close to Earth's orbit about the Sun. Whilst we were most satisfied by its behaviour, I should now like to report upon an altogether more surprising consequence of its engine's action.

Last time we saw how we could use a sequence of Householder transformations to reduce a symmetric real matrix M to a symmetric tridiagonal matrix, having zeros everywhere other than upon the leading, upper and lower diagonals, which we could then further reduce to a diagonal matrix Λ using a sequence of Givens rotations to iteratively transform the elements upon the upper and lower diagonals to zero so that the columns of the accumulated transformations V were the unit eigenvectors of M and the elements on the leading diagonal of the result were their associated eigenvalues, satisfying

    M × V = V × Λ

and, since the transpose of V is its own inverse

    M = V × Λ × V^T

which is known as the spectral decomposition of M.
Unfortunately, the way that we used Givens rotations to diagonalise tridiagonal symmetric matrices wasn't particularly efficient and I concluded by stating that it could be significantly improved with a relatively minor change. In this post we shall see what it is and why it works.

Hello there Sir R-----! Come join me by the hearth for a dram of warming spirits! I trust that this cold spell has not chilled your desire for a wager?

Good man! Good man!

I must say that the contrast between the warmth of this fire and the frost outside brings most vividly to my mind an occasion during my tenure as the Empress's ambassador to the land of Oz; specifically the time that I attended King Quadling Rex's winter masked ball during which his southern palace was overrun by an infestation of Snobbles!

Some years ago we saw how we could use the Jacobi algorithm to find the eigensystem of a real valued symmetric matrix M, which is defined as the set of pairs of non-zero vectors v_i and scalars λ_i that satisfy

    M × v_i = λ_i × v_i

known as the eigenvectors and the eigenvalues respectively, with the vectors typically restricted to those of unit length in which case we can define its spectral decomposition as the product

    M = V × Λ × V^T

where the columns of V are the unit eigenvectors, Λ is a diagonal matrix whose i^th diagonal element is the eigenvalue associated with the i^th column of V and the T superscript denotes the transpose, in which the rows and columns of the matrix are swapped.
You may recall that this is a particularly convenient representation of the matrix since we can use it to generalise any scalar function to it with

    f(M) = V × f(Λ) × V^T

where f(Λ) is the diagonal matrix whose i^th diagonal element is the result of applying f to the i^th diagonal element of Λ.
You may also recall that I suggested that there's a more efficient way to find eigensystems and I think that it's high time that we took a look at it.