Discrete dynamical systems with Haskell

A discrete dynamical system is given by xₖ₊₁ = Axₖ for k = 0, 1, 2, … with a matrix A and an initial state vector x₀. Often times it’s useful to observe the long-term behavior of a system by finding its steady state vector—that is, a vector q for which q = Aq. We can do this in Haskell by declaring a linear transformation x ↦ Ax and finding its fixed point given the initial state. In other words, we will repeatedly apply the linear transformation, starting with the initial state, until we converge to a single vector (conventionally, one would just solve the equation q = Aq). Conveniently enough, Haskell defines some data structures and functions for computing fixed points, and I’ll cheat a bit by using a linear algebra library, although implementing one from scratch might be left to a future post.

For simplicity, we’ll work with 3 × 3 matrices and 3-vectors. First, we must define a data structure for state vectors and systems.

{-# LANGUAGE RecordWildCards #-}
module Main where

import Data.Ratio
import Linear.Matrix
import Linear.V3
import Linear.Epsilon
import Control.Monad.Fix

type State = V3 Double

data System = System {
    x0 :: State,
    a :: M33 Double
}

Now for a brief interlude on fixed points. The module Control.Monad.Fix defines a function fix that infinitely applies a function with a starting parameter, like so: f(f(…f(x)…)).

fix :: (a -> a) -> a
fix f = let x = f x in x

We can “fix” recursive functions by redefining them with an auxiliary parameter rec that is used for the recursive call instead of the function itself. For example, the factorial function is normally defined as the following.

f n = if n == 0 then 1 else n * f $ n - 1

Our “fixed” version is:

f rec n = if n == 0 then 1 else n * rec $ n - 1

This essentially converts a recursive function to a non-recursive one, because fix will keep passing said function as rec. As a result, fix enables support for recursion in primitive languages like the lambda calculus. Furthermore, notice that we always need to have a base case/terminating condition, so that fix doesn’t recur forever.

So how do we get from fix to a system’s steady state vector? First, we need to define a terminating condition: we want our linear transformation t to terminate when q = Aq. However, we can’t check for strict equality because (1) we’re dealing with floating-point numbers and (2) we’re checking for convergence with this method, and not actual equality. With that in mind, let’s define an operator === that returns whether or not two vectors are close enough to each other to be considered equal.

(===) :: State -> State -> Bool
a === b = nearZero $ a - b

Now we’re ready to define a function steadyState that takes a system and returns its steady state vector. Thus, we will pass fix a function t that will return the result of computing Ax if it === x, otherwise it will recur on Ax. Once we pass fix the initial state vector as well, it will magically converge to the steady state of the system.

steadyState :: System -> State
steadyState System{..} = fix t x0 where
    t rec xk = let xk1 = a !* xk in
        if xk1 === xk then
            xk1 -- We found it, so terminate
        else
            rec xk1 -- Keep going

Now that we’ve defined everything, let’s test it on a discrete-time Markov chain, with an initial probability vector and stochastic matrix.

dtmc = System {
    x0 = V3 1 0 0,
    a = V3 (V3 0.5 0.2 0.3) -- Row-major representation
           (V3 0.3 0.8 0.3)
           (V3 0.2 0.0 0.4)
}

Let’s load up all of this code in ghci and run it.

$ ghci steadystate.hs
...
*Main> steadyState dtmc
V3 0.30000038146975283 0.5999988555908204 0.10000076293942693

That’s reasonably close to the vector (.3, .6, .1), which is the right answer according to my textbook. The source for this project is on github.