Sir R-----! Come join me for a glass of chilled wine! I have a notion that you're in the mood for a wager. What say you?

I knew it!

I have in mind a game of dice that reminds me of my time as the Russian military attaché to the city state of Coruscant and its territories during the traitorous popular uprising fomented by the blasphemous teachings of a fundamentalist religious sect known as the Jedi.

Last time we took a look at basis function interpolation with which we approximate functions from their values at given sets of arguments, known as nodes, using weighted sums of distinct functions, known as basis functions. We began by constructing approximations using polynomials before moving on to using bell shaped curves, such as the normal probability density function, centred at the nodes. The latter are particularly useful for approximating multi-dimensional functions, as we saw by using multivariate normal PDFs.
An easy way to create rotationally symmetric functions, known as radial basis functions, is to apply univariate functions that are symmetric about zero to the distance between the interpolation's argument and their associated nodes. PDFs are a rich source of such functions and, in fact, the second bell shaped curve that we considered is related to that of the Cauchy distribution, which has some rather interesting properties.

Recall that the Baron’s game consisted of guessing under which of a pair of cups was to be found a token for a stake of four cents and a prize, if correct, of one. Upon success, Sir R----- could have elected to play again with three cups for the same stake and double the prize. Success at this and subsequent rounds gave him the opportunity to play another round for the same stake again with one more cup than the previous round and a prize equal to that of the previous round multiplied by its number of cups.

A few years ago we saw how we could approximate a function f between pairs of points (x_i, f(x_i)) and (x_i+1, f(x_i+1)) by linear and cubic spline interpolation which connect them with straight lines and cubic polynomials respectively, the latter of which yield smooth curves at the cost of somewhat arbitrary choices about their exact shapes.
An alternative approach is to construct a single function that passes through all of the points and, given that n^th order polynomials are uniquely defined by n+1 values at distinct x_i, it's tempting to use them.

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.