dot v w = sum $ zipWith (*) v w
matrixFun f vs = dot vs <$> f
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
echelon = \case
[] -> []
xs@([]:_) -> xs
xs | (zs, nzs) <- partition ((== 0) . head) xs -> case nzs of
[] -> map (0:) $ echelon $ tail <$> zs
(pivot:pivotTail):nzt -> (1:pt) : rest
where
pt = map (/pivot) pivotTail
rest = map (0:) $ echelon $ tail <$> zs <|> elim <$> nzt
elim (h:t) = zipWith (-) t $ map (h*) pt
reduce = reverse . reduceReverse . reverse . echelon
reduceReverse [] = []
reduceReverse (row:rest) = row : (go <$> reduceReverse rest) where
k = length $ takeWhile (0 ==) row
go other
| k < length other = zipWith (-) other $ map ((other!!k)*) row
| otherwise = other
echelon' dust = \case
[] -> []
xs@([]:_) -> xs
xs | (zs, nzs) <- partition ((== 0) . head) (map dust <$> xs) -> case nzs of
[] -> map (0:) $ echelon' dust $ tail <$> zs
(pivot:pivotTail):nzt -> (1:pt) : rest
where
pt = map (/pivot) pivotTail
rest = map (0:) $ echelon' dust $ tail <$> zs <|> elim <$> nzt
elim (h:t) = zipWith (-) t $ map (h*) pt
reduce' dust = reverse . reduceReverse . reverse . echelon' dust
reduceReverse [] = []
reduceReverse (row:rest) = row : (go <$> reduceReverse rest) where
k = length $ takeWhile (0 ==) row
go other
| k < length other = zipWith (-) other $ map ((other!!k)*) row
| otherwise = other
forwardElim' [] = []
forwardElim' xs = p : rest where
(zs, p@((pivot:pivotTail), _):nzs) = partition ((== 0) . head . fst) xs
rest = first (0:) <$> forwardElim' (elim <$> (zs ++ nzs))
elim (h:t, (idx, ks)) = let k = h*(1/pivot) in
(zipWith (-) t $ map (k*) pivotTail, (idx, k:ks))
buildLinv acc mults = row:acc where
idRow = replicate (length mults) 0 ++ [1]
row = foldl (zipWith (-)) idRow $ (++ repeat 0) <$> scaledRows
scaledRows = zipWith map (map (*) mults) acc
scale a = map (a*)
backElim u inv = revrev $ go (revrev u) (revrev inv) where
revrev = reverse . map reverse
go = \cases
[] [] -> []
-- `one` might not actually be 1 due to floating point.
((one:_):rest) (row:rows) -> row : go (tail <$> rest) rows' where
rows' = zipWith (zipWith (-)) rows
$ map (`scale` row) (head <$> rest)
matrixInvert m = matrixMul (backElim (applyD u) (applyD lInv)) p where
(u, info) = unzip $ forwardElim' $ zip m $ (, []) <$> [0..]
(perm, lRaw) = unzip info
n = length u
p = take n <$> zipWith (++) (map (`replicate` 0) perm) (repeat (1 : repeat 0))
lInv = reverse $ map (take n . (++ repeat 0)) $ foldl buildLinv [] lRaw
applyD = zipWith scale $ (1/) <$> leads u
leads = \case
[] -> []
row:rest -> head row : leads (tail <$> rest)
instance (Eq a, Ring a, Field a) => Field (Matrix a) where
recip = \case
Matrix m -> Matrix $ matrixInvert m
Scalar a -> Scalar $ recip a
norm = sqrt . sum . map (^2)
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
\[ \begin{align} x_{n+1} &= x_n + y_n \\ y_{n+1} &= x_n \end{align} \]
We express this with 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.
We apply \(T\) ten times to the initial state:
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 \(A^n\) is trivial to write down because \(A\) is a diagonal matrix, so:
\[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}\]
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\).