Skip to content

Latest commit

 

History

History
879 lines (762 loc) · 30 KB

kalman-folding-04-008-streams-observables.org

File metadata and controls

879 lines (762 loc) · 30 KB

Kalman Folding 4: Streams and Observables (Review Draft)

1 Abstract

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.

2 Kalman Folding in the Wolfram Language

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 $\mathbold{x}$ are $\left[\left[\mathbold{x}\right]\right] \stackrel{\text{\tiny def}}{=} (\mathcal{X}, n×{1})$ and of $\mathbold{z}$ are $\left[\left[\mathbold{z}\right]\right] \stackrel{\text{\tiny def}}{=} (\mathcal{Z}, b×{1})$, then

Dimensional arguments, regarding both matrix dimensions and physical dimensions, are invaluable for checking code and derivations in this topic at-large.

2.1 A Test Example

In the following example, the observations $\mathbold{z}$ are $1×{1}$ matrices, equivalent to scalars, so $b=1$.

The function in equation \ref{eqn:kalman-cume-definition} lambda-lifts/lmlf $\mathbold{Z}$, meaning that it is necessary to call /kalmanStatic with a constant $\mathbold{Z}$ to get the actual accumulator function used in folds. This is desirable to reduce coupling between the accumulator function and its calling environment.

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 $\mathbold{P}$. For details about this test case, see the first paper in the series, Kalman Folding, Part 1.klfl

Below, we reproduce these values exactly, to the bit level, by running kalmanStatic over lazy streams and asynchronous observables.

3 Types for Kalman Folding

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:

  1. another function (called the accumulator function)
  2. an initial instance of type $S$
  3. 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 $\text{Sequence}\left[\,\text{Observation}\,\right]$

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.

4 Over Lazy Streams and Asynchronous Observables

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.

4.1 Folding Over Lazy Streams

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 $n$:

Calling, say, integersFrom[42] produces {42, integersFrom[42 + 1]&}, a pair of an integer, $42$, and another stream, 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 $2×{2}$ identity matrix, returns a lazy stream of matrices full of Fibonacci numbers:

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.

4.1.1 Disperse :: List $→$ Stream

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.

4.1.2 Reify :: Stream $→$ List

We need to go the other way, too; don’t call this on a stream of infinite length:

4.1.3 foldStream

Our equivalent for Wolfram’s FoldList is foldStream.uncl Its type is similar

Here is an implementation:

4.1.4 Test

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.

4.2 Folding Over an Asynchronous Observable

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

4.2.1 Subscribe :: Observable $→$ Observer $→$ Null

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:

4.2.2 Dispense :: List $→$ Observable

The following is a specification of a task to run. Nothing happens till you subscribe something to it.

4.2.3 Harvest :: Observable $→$ List

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.

4.2.4 foldObservable

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.

4.2.5 Test

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}.

5 Concluding Remarks

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 $×$ symbols to clarify matrix multiplication. They have no other significance and we can always write matrix multiplication just by juxtaposing the matrices. clos https://en.wikipedia.org/wiki/Closure_(computer_programming) cold This convention only models so-called cold observables, but it’s enough to demonstrate Kalman’s working over them. cons This is quite similar to the standard — not Wolfram’s — definition of a list as a pair of a value and of another list. cova We use the terms covariance for matrices and variance for scalars. csoc https://en.wikipedia.org/wiki/Separation_of_concerns ctsc https://en.wikipedia.org/wiki/Catastrophic_cancellation dstr http://tinyurl.com/ze6qfb3 elib Brookner, Eli. Tracking and Kalman Filtering Made Easy, New York: Wiley, 1998. http://tinyurl.com/h8see8k fldl http://tinyurl.com/jmxsevr fwik https://en.wikipedia.org/wiki/Fold_%28higher-order_function%29 gama https://en.wikipedia.org/wiki/Gauss%E2%80%93Markov_theorem intr http://introtorx.com/ jplg JPL Geodynamics Program http://www.jpl.nasa.gov/report/1981.pdf just justified by the fact that $\mathbold{D}$ is a diagonal matrix that commutes with all other products, therefore its left and right inverses are equal and can be written as a reciprocal; in fact, $\mathbold{D}$ is a $1×{1}$ matrix — effectively a scalar — in all examples in this paper klde B. Beckman, Kalman Folding 3: Derivations, http://vixra.org/abs/1607.0059. klf1 B. Beckman, Kalman Folding, Part 1, http://vixra.org/abs/1606.0328. klf2 B. Beckman, Kalman Folding 2: Tracking and System Dynamics, http://vixra.org/abs/1606.0348. klf3 B. Beckman, Kalman Folding 3: Derivations, http://vixra.org/abs/1607.0059. klf4 B. Beckman, Kalman Folding 4: Streams and Observables, to appear. klfl B. Beckman, Kalman Folding, Part 1, http://vixra.org/abs/1606.0328. layi https://en.wikipedia.org/wiki/Fundamental_theorem_of_software_engineering lmbd Many languages use the keyword lambda for such expressions; Wolfram uses the name Function. lmlf https://en.wikipedia.org/wiki/Lambda_lifting lsqo LINQ’s Standard Query Operators lssq https://en.wikipedia.org/wiki/Least_squares ltis http://tinyurl.com/hhhcgca matt https://www.cs.kent.ac.uk/people/staff/dat/miranda/whyfp90.pdf mcmc https://en.wikipedia.org/wiki/Particle_filter mond https://en.wikipedia.org/wiki/Monad musc http://www1.cs.dartmouth.edu/~doug/music.ps.gz ndim https://en.wikipedia.org/wiki/Nondimensionalization patt http://tinyurl.com/j5jzy69 pseu http://tinyurl.com/j8gvlug rasp http://www.wolfram.com/raspberry-pi/ rcrn https://en.wikipedia.org/wiki/Recurrence_relation rsfr http://rosettacode.org/wiki/Loops/Foreach rxbk http://www.introtorx.com/content/v1.0.10621.0/07_Aggregation.html scan and of Haskell’s scans and folds, and Rx’s scans and folds, etc. scla http://tinyurl.com/hhdot36 scnd A state-space form containing a position and derivative is commonplace in second-order dynamics like Newton’s Second Law. We usually employ state-space form to reduce \(n\)-th-order differential equations to first-order differential equations by stacking the dependent variable on $n-1$ of its derivatives in the state vector. scnl http://learnyouahaskell.com/higher-order-functions stsp https://en.wikipedia.org/wiki/State-space_representation uncl The initial uncial (lower-case) letter signifies that we wrote this function; it wasn’t supplied by Wolfram. wfld http://reference.wolfram.com/language/ref/FoldList.html?q=FoldList wlf1 http://tinyurl.com/nfz9fyo wlf2 http://rebcabin.github.io/blog/2013/02/04/welfords-better-formula/ wolf http://reference.wolfram.com/language/ zarc Zarchan and Musoff, Fundamentals of Kalman Filtering, A Practical Approach, Fourth Edition, Ch. 4