dot v w = sum $ zipWith (*) v w
matrixAdd = zipWith $ zipWith (+)
matrixScalarMul a = map $ map (a*)
matrixMul a b = map (\row -> dot row <$> transpose b) a
onDiag f a = go <$> zipWith splitAt [0..] a where
go (as, b:bs) = as ++ f b : bs
data Matrix a = Matrix [[a]] | Scalar a deriving Eq
instance Ring a => Ring (Matrix a) where
fromInteger z = Scalar $ fromInteger z
(+) = \cases
(Scalar a) (Scalar b) -> Scalar $ a + b
(Scalar a) (Matrix b) -> Matrix $ onDiag (a+) b
(Matrix a) (Scalar b) -> Matrix $ onDiag (b+) a
(Matrix a) (Matrix b) -> Matrix $ matrixAdd a b
a - b = a + negate b
(*) = \cases
(Scalar a) (Scalar b) -> Scalar $ a * b
(Scalar a) (Matrix b) -> Matrix $ matrixScalarMul a b
(Matrix a) (Scalar b) -> Matrix $ matrixScalarMul b a
(Matrix a) (Matrix b) -> Matrix $ matrixMul a b
negate = \case
(Scalar a) -> Scalar $ negate a
(Matrix a) -> Matrix $ map negate <$> a
instance Show a => Show (Matrix a) where
show = \case
Scalar a -> show a
Matrix a -> unlines $ unwords . map show <$> a
Rabbits
Suppose there are \(x_0\) pairs of adult rabbits and \(y_0\) pairs of kits. Each year, each kit grows up and becomes an adult, and each pair of adult rabbits gives birth to a pair of kits.
Then we can compute the number of pairs of rabbits and kits after a given number of years by iterating \(x_{n+1} = x_n + y_n\) and \(y_{n+1} = x_n\), equations which we write using matrices:
\[ \begin{bmatrix} x_{n+1} \\ y_{n+1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} x_n \\ y_n \end{bmatrix} \]
This is an example of a discrete linear dynamical system. The transition matrix is:
\[ T = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \]
and the \(n\)th state is:
\[ \begin{bmatrix} x_n \\ y_n \end{bmatrix} = T^n \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} \]
Let’s simulate a decade of this model from the initial state: \( \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \)
We throw together a barebones matrix library. A matrix is represented with a list of lists, and there are no checks that dimensions match.
Then start from the initial state, and apply \(T\) to it 10 times:
rabbit = Matrix [[1, 1], [1, 0]] putStr $ unlines $ map show $ take 11 $ iterate (rabbit *) $ Matrix [[0], [1]]
"Then a Miracle Occurs"
What is the long-term behaviour of this system?
We answer this with the help of a few magic numbers. Let:
\[\begin{align}
\phi &= \frac{1+\sqrt{5}}{2} \\
\psi &= \frac{1-\sqrt{5}}{2} = 1 - \phi = -\phi^{-1} \\
d &= \sqrt{5} = \phi -\psi \\
A &= \begin{bmatrix} \phi & 0 \\ 0 & \psi \end{bmatrix} \\
P &= \frac{1}{d} \begin{bmatrix} 1 & -\psi \\ -1 & \phi \end{bmatrix} \\
\end{align}\]
phi = 1.618034 psi = 1 - phi d = phi - psi a = Matrix [[phi, 0], [0, psi]] p = Matrix $ map (map (*recip d)) [[1, -psi], [-1, phi]]
We can check:
\[ P^{-1} = \begin{bmatrix} \phi & \psi \\ 1 & 1 \end{bmatrix} \]
p1 = Matrix [[phi, psi], [1, 1]] p1 * p
We can also check:
\[ T = P^{-1} A P \]
p1 * a * p
Then:
\[T^n = (P^{-1} A P)^n = P^{-1} A^n P = \frac{1}{d}
\begin{bmatrix}
\phi^{n+1}-\psi^{n+1} & \phi^n-\psi^n \\
\phi^n-\psi^n & \phi^{n-1}-\psi^{n-1}
\end{bmatrix}\]
where \(A^n\) is trivial to write down because \(A\) is a diagonal matrix. Therefore:
\[ \begin{align} x_n &= ((\phi^{n+1} - \psi^{n+1}) x_0 + (\phi^{n} - \psi^{n}) y_0) / d \\ y_n &= ((\phi^{n} - \psi^{n}) x_0 + (\phi^{n-1} - \psi^{n-1}) y_0) / d \end{align} \]
When \(x_0 = 0, y_0 = 1\). we have:
\[ x_n = \frac{1}{\sqrt{5}} \left( \left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n \right) \]
This is the famous Binet’s formula for the Fibonacci numbers.
binet n = recip d * (phi^n - psi^n) binet <$> [0..10]
But using magic numbers that happen to fit perfectly is like pulling a rabbit out of a hat. How were they found?
I first learned to derive Binet’s formula via generating functions (a rabbit hole that is worth going down). Over the next few sections, we’ll see that linear algebra provides another path, and one that is more general when the recurrences are linear: instead of a sequence of single numbers, we can handle a sequence of vectors of numbers. We shall learn that \(\phi\) and \(\psi\) are the eigenvalues of \(T\), and \(A\) is the Jordan normal form of \(T\).