In Kalman Folding, Part 1,klfl we present basic, static Kalman filtering as a functional fold, highlighting the unique advantages of this form for deploying test-hardened code verbatim in harsh, mission-critical environments. In that paper, all examples folded over arrays in memory for convenience and repeatability. That is an example of developing filters in a friendly environment.
Here, we prototype a couple of less friendly environments and demonstrate exactly the same Kalman accumulator function at work. These less friendly environments are
- lazy streams, where new observations are computed on demand but never fully realized in memory, thus not available for inspection in a debugger
- asynchronous observables, where new observations are delivered at arbitrary times from an external source, thus not available for replay once consumed by the filter
Streams are a natural fit for integration of differential equations, which often arise in applications. As such, they enable unique modularization for all kinds of filters, including non-linear Extended Kalman Filters.
The fact that the Kalman accumulator function gives bit-for-bit identical results in all cases gives us high confidence that code developed in friendly environments will behave as intended in unfriendly environments. This level of repeatability is available only because of functional decomposition, which minimizes the coupling between the accumulator function and the environment and makes it possible to deploy exactly the same code, without even recompilation, in all environments.
In this series of papers, we use the Wolfram languagewolf because it excels at concise expression of mathematical code. All examples in these papers can be directly transcribed to any modern mainstream language that supports closures. For example, it is easy to write them in C++11 and beyond, Python, any modern Lisp, not to mention Haskell, Scala, Erlang, and OCaml. Many can be written without full closures; function pointers will suffice, so they are easy to write in C. It’s also not difficult to add extra arguments to simulate just enough closure-like support in C to write the rest of the examples in that language.
In Kalman Folding,klfl we found the following elegant formulation for the accumulator function of a fold that implements the static Kalman filter:
\noindent where
\noindent and all quantities are matrices:
-
$\mathbold{z}$ is a${b}×{1}$ column vector containing one multidimensional observation -
$\mathbold{x}$ is an${n}×{1}$ column vector of model states -
$\mathbold{Z}$ is a${b}×{b}$ matrix, the covariance of observation noise -
$\mathbold{P}$ is an${n}×{n}$ matrix, the theoretical covariance of$\mathbold{x}$ -
$\mathbold{A}$ is a${b}×{n}$ matrix, the observation partials -
$\mathbold{D}$ is a${b}×{b}$ matrix, the Kalman denominator -
$\mathbold{K}$ is an${n}×{b}$ matrix, the Kalman gain
In physical or engineering applications, these quantities carry physical
dimensions of units of measure in addition to their matrix dimensions as numbers
of rows and columns.
If the physical and matrix dimensions of
Dimensional arguments, regarding both matrix dimensions and physical dimensions, are invaluable for checking code and derivations in this topic at-large.
In the following example, the observations
The function in equation \ref{eqn:kalman-cume-definition}
lambda-lifts/lmlf
In Wolfram, this function is
We test it on a small case
\noindent expecting results within one or two standard deviations of the ground
truth $\aleph=\begin{bmatrix}-3& 9& -4& -5\end{bmatrix}^\intercal$, where the
standard deviations can be found as square roots of the diagonal
elements of
Below, we reproduce these values exactly, to the bit level, by running kalmanStatic over lazy streams and asynchronous observables.
Kalman and all its variants are examples of statistical function inversion. We have models that predict outcomes from inputs; we observe outcomes and want estimates of the inputs. Structurally, all such incremental model inversions take a pair of a state estimate (with uncertainty) and an observation, and produce a new state estimate (with uncertainty). Such an inverted model has signature, using a type notation similar to that of Haskell or Scala
\noindent where the return type is on the far right and the other types that appear before arrows are the types of input arguments. This function signature is exactly that required for the first argument of a functional fold (more precisely, a left fold). The signature of fold is as follows:
Read this, abstractly, as follows
\emph{Fold} over types
$S$ and$T$ is a function that takes three arguments:
- another function (called the accumulator function)
- an initial instance of type
$S$ - a sequence of instances of type
$T$ and produces an instance of type
$S$ . The accumulator function, in turn, is a binary function that takes an$S$ and a$T$ and produces an$S$ .
More concretely, In the context of Kalman filtering:
\noindent where the types Accumulation and Observation are arbitrary.
It’s the job of Fold to pass the elements of the input sequence to the
accumulator function one observation at a time, and to maintain and ultimately
return the final accumulation. The second argument to Fold is
the desired, initial value of the accumulation. The
third and final argument to Fold is the sequence of observations, of type
Fold looks like a trinary function of an accumulator function, an initial accumulation, and a sequence, yielding an accumulation. Folds thus have the following type:
\noindent where Sequence can be List, Stream, Observable, or any type that can be accessed sequentially.
The accumulator function knows nothing about the source of the observations. If we can figure out how to implement Fold and FoldList for things other than List, we will have Kalman filtering over those sources, too.
The following are research-grade sketches of implementations of Fold over lazy streamsmusc and asynchronous observables.intr They provide just enough to support the Kalman-folding examples.
Represent a lazy stream as a pair of a value and a thunk (function of no arguments).cons The thunk must produce another lazy stream when called. Such a stream can be infinite in abstract length because the elements of the stream are only concretized in memory when demanded by calling thunks.
Streams are a natural fit for integrals of differential equations. We see in other papers of this series how we can use them to deeply modularize filters over rich non-linear models. In this paper, we show only how to fold a linear Kalman filter over a stream.
By convention, a finite stream has a Null thunk at the end. Thus, the empty
stream, obtained by invoking such a thunk, is Null[], with square brackets
denoting invocation with no arguments.
One of Wolfram’s notations for a literal thunk is an expression with an
ampersand in postfix position. An ampersand turns the expression to its left
into a thunk. For instance, here’s a function that returns an infinite stream of
natural numbers starting at
Calling, say, integersFrom[42] produces {42, integersFrom[42 + 1]&}, a pair
of an integer, integersFrom[42+1]&. We get the
stream by extracting the second part of the pair via Wolfram’s double-bracket notation
\noindent and then call it with empty brackets (it’s a thunk, and takes no arguments):
\noindent and so on. We can get a few more by repeating the process
\noindent but the best way to extract values from streams is to write recursive functions to demand any number of elements from the head. The variety of such functions, which include map, select, fold, is well known, large, and identical across lists, streams, observables, and, in fact, any collection that can support a next operator. A good, contemporary full-service library for collection types is LINQ’s Standard Query Operators (SQO),lsqo. If building up a library from the present prototype level into something of product grade, presentable to intolerant users, the SQO are an excellent framework to emulate.
As another example, the following function, when called with an appropriate
input, say the
Here is an explicit invocation a few values down:
\noindent the point being that lazy streams are versatile.
We now write bi-directional conversions between streams and lists so we can test an example, then we write foldStream.
We’ll need a way to convert a list into a stream. There are three cases: an empty list, a singleton list, and the inductive or recursive case.
We need to go the other way, too; don’t call this on a stream of infinite length:
Our equivalent for Wolfram’s FoldList is foldStream.uncl Its type is similar
Here is an implementation:
Test it over the dispersion of the example data:
The only changes to the earlier fold over lists is the initial call of disperse to
convert the test case into a stream, and the final postfix call // reify
to turn the result back into a list for display. The final results are identical
to those in equation \ref{eqn:kalman-filter-results}, but we see all the
intermediate results as well, confirming that Kalman folds over observations one
at a time. We would have seen exactly the same output had we called FoldList
instead of Fold over lists above.
Just as FoldList produces a list from a list, and foldStream produces a stream from a stream, foldObservable produces an observable from an observable. Its full signature is
Lists provide data elements distributed in space (memory). Lazy streams provide data in constant memory, but distributed in a kind of virtual time, delivered when demanded, the way a debugger fakes time. Observables provide data elements distributed asynchronously in real time. To consume elements of an observable, subscribe an observer to it. An observer has a callback function, and the observable will invoke the callback for each observation, asynchronously, as the observation arrives. The callback function takes a single argument that receives the observation.
We do not develop observables fully, here. For that, see a reference like Campbell’s Intro to Rx.intr Instead, we content ourselves with just enough to demonstrate Kalman folding over them and, as with lazy streams, a way to get back and forth from lists.
We model observables as stateful thunks that produce new values every time they’re invoked, then invoke the thunks inside asynchronous Wolfram tasks that start at the moment some observer subscribes.cold
Wolfram supplies a primitive, RunScheduledTask, for evaluating expressions asynchronously, once per second by default. The expression that we pass to RunScheduledTask, just calls the observer on the evaluated observable:
The following is a specification of a task to run. Nothing happens till you subscribe something to it.
Set up a conventional, external variable, r$, so that we can interactively
look at the results in a Wolfram Dynamic[r$] form. Our harvest
subscribes an observer that appends observations to a list held in r$.
Semicolon-separated expressions are sequenced, as with Scheme’s begin or
Lisp’s progn.
We must eventually clean up the tasks and the external variable.
The concrete type of foldObservable is obvious: just replace Stream with Observable in a copy of the type of foldStream.
One might ask about the appropriate generalization of higher-order types like this, where we could go up a level, parameterize on types like Stream and Observable, and make the concrete types of foldStream and foldObservable instances of that higher, parameterized type. This is a sensible question, and the answer leads to category theory and monads,mond out of scope for this paper.
This implementation isn’t hygienic: it uses global variables (suffixed with $
signs). It’s just enough to test Kalman folding over observables.
The following call has the same shape as our call of foldStream above, except calling dispense instead of disperse and harvest instead of reify.
The results are exactly the same as in equation \ref{eql:full-big-results}.
With prototypes for foldStream and foldObservable, we have demonstrated Kalman folding with exactly the same accumulator function over wildly different data-delivery environments. This demonstrates the primary thesis of this series of papers: that writing filters as functional folds enables verbatim deployment of code in both friendly, synchronous environments with all data in memory, and unfriendly asynchronous environments using only constant memory. Verbatim means with no changes at all, not even recompilation.
We have tested these prototypes against bigger examples like the tracking exampleklf2 and the accelerometer example,klfl and there are no surprises.
affn https://en.wikipedia.org/wiki/Affine_transformation
bars Bar-Shalom, Yaakov, et al. Estimation with applications to tracking and navigation. New York: Wiley, 2001.
bier http://tinyurl.com/h3jh4kt
bssl https://en.wikipedia.org/wiki/Bessel’s_correction
busi https://en.wikipedia.org/wiki/Business_logic
cdot We sometimes use the center dot or the