\documentstyle[12pt]{article}
\topmargin -0.6in
\textwidth 6.0in
\textheight 9.0in
\oddsidemargin 0.1in
\title{Efficient Parallel Implementation of \\ Target Tracking Kalman Filter} 
\author{Radhakisan S. Baheti and David R. O'Hallaron\\
GE Research and Development Center, Schenectady, NY 12301}
\begin{document}
\bibliographystyle{ieeetr}
\maketitle
\begin{abstract} 
A parallel algorithm for solving an $n$-state Kalman filter on an
$(n+1)$-cell linear array is described.  The approach is to update the
columns of the filter covariance matrix in parallel to balance the
computations and minimize the  communications.
 The algorithm is the basis for
the mapping of an extended-Kalman filter on the Warp linear array processor.
The Warp implementation of the filter is written in a high-level language and
achieves a measured speedup of 400 over the same filter running
on a Sun workstation.  The filter mapping is very efficient because
8-cell Warp implementation (parallel algorithm) achieves a measured
speedup of 7 over 1-cell implementation (sequential algorithm).
Linear arrays have traditionally been applied to regular and uniform
computations such as matrix-matrix multiplications.  The problem of
interest is the mapping of irregular and
non-uniform computations required by the extended-Kalman filter onto
a linear processor.
With the availability of high-performance multi-computer chips
in the near future, the results are directly applicable to 
on-board target tracking.
\end{abstract}
\section{Introduction}
Kalman filtering is a general and popular technique for recursively
estimating the state variables of a dynamic system from noisy
measurements. The solution of the Kalman equations is expensive, in
general requiring $O(n^3)$ floating-point computations. This
time-complexity has motivated the design of a number of parallel
algorithms, \cite{andrews81,chen86,jover86,shaffer87b,sung87}, 
all of which assume 2-D arrays of computational cells. Such arrays
require numerous external (and expensive) high-speed data transmission
channels, and cannot always be synchronized with a single global clock
\cite{fisher83}. In general, 2-D arrays of cells are more costly and
difficult to construct than linear arrays of cells. 

This paper describes a parallel algorithm for solving an $n$-state
Kalman filter using an $(n+1)$-cell linear array.  The algorithm is
the basis for the mapping of a 9-state target tracking filter on the
Warp computer \cite{annaratone87e}, a linear array of 10 or more
identical, powerful, and programmable cells connected to a
general-purpose host.  The Warp implementation is written in a
high-level language and achieves a measured sample time of 125 $\mu
s$, which is a speedup of 400 over the same filter running on a
Sun 3/60 workstation.

Section 2 describes a Kalman filter for target tracking applications.
A parallel filter algorithm for mapping on a linear array of computational
cells is given in Section 3. In Section 4, the
mapping of the algorithm on Warp and its performance is described.

\section{Target Tracking Kalman Filter}

In this section we describe a Kalman filter for target tracking applications.
Let $n$ denote the number of elements in the state vector and $m$ denote the
number of observations in the measurement vector.  In this paper, we assume
a target tracking filter with $n=9$ states and $m=3$ measurements.
The filter is used to estimate from noisy measurements the position,
velocity and acceleration of a maneuvering target.  The target may be 
surface, subsonic aerial, or supersonic aerial, such as a ship, an 
aircraft, or a missile.  The radar, sonar, or optical instruments measure
the target range and azimuth and elevation angles at high data rates.  

The tracking filter is based on the laws of motion and a stochastic
acceleration model.  It is assumed that the target moves at constant
velocity if there is no maneuver or atmospheric turbulence.  Evasive
maneuvers and other disturbances are modelled as perturbations upon the
constant-velocity trajectory.

The target acceleration is correlated in time.  For example, a lazy turn
will use a highly correlated acceleration.  On the other hand, evasive
maneuvers will be less correlated in time.  The correlation function
$ \gamma ( \tau ) $ for the target acceleration $ a(t) $ can be 
represented by 
\begin{displaymath}
\gamma ( \tau ) = E\{ a(t)a(t+ \tau ) \} = \sigma _ m  ^ 2 e ^{ - \alpha \mid \tau \mid }
\end{displaymath}
where $E$ denotes the expected value, $ \sigma _ m  ^ 2 $ denotes the variance
of the target acceleration, and $ 1 / \alpha $ is the maneuver time constant.
The target dynamic model in one dimension  has the form
\begin{displaymath}
d/dt  \left [
       \begin{array}{c}
        x \\
        \dot x \\
        \ddot x
       \end{array} \right ] =
\left [ \begin{array}{c c c}
      0 & 1 & 0 \\
      0 & 0 & 1 \\
      0 & 0 & - \alpha
     \end{array}
      \right ]  \left [ \begin{array}{c}
         x \\
         \dot x \\
         \ddot x \end{array}
    \right ]
 + \left [ \begin{array}{c}
    0 \\
    0 \\
    1 
 \end{array}
 \right ]   
 w(t)    
\end{displaymath}
where $ x,  \dot x, $ and $ \ddot x $ denote the target position, velocity,
and acceleration at time $ t $, and  $ \{ w(t) \} $ is a white Gaussian
noise process with zero mean and variance $ \sigma _ m ^ 2 $.  Note that
the uncertainty in the target acceleration is represented by a first-order
Gauss-Markov process.

The filtering may be performed in polar coordinates with range, azimuth, 
and elevation angles, and their derivatives as the state variables.  In 
polar coordinates the measurement system model is linear and uncoupled.
However, even for a target moving in a straight line at constant speed,
the dynamic model is nonlinear, and the linearization leads to
large errors. In Cartesian coordinates, a constant-speed target does obey
a a linear law, but the measurement equations are nonlinear.  Both 
coordinate systems have been successfully used in many target tracking
applications.  Target state estimation is normally part of a major goal
to predict the target position and improve the accuracy of a fire control 
system.
A discrete-time equation of target motion in Cartesian coordinates
can be expressed by a state variable model given by \cite{baheti86}. 
\begin{eqnarray}
{\bf x}(t+1) & = & \Phi ( \Delta t, {\bf \Lambda })  {\bf x}(t) + {\bf w}(t) \nonumber \\
{\bf y}(t)   & = & {\bf h(x}(t)) + {\bf v}(t), ~~t=0, 1, \cdots \label{model}
\end{eqnarray}
where $t$ is the normalized discrete time, $ \Delta t $ denotes the
sampling, $ {\bf x} (t) $ is a nine-dimensional state vector with position,
velocity, and acceleration in each of the Cartesian coordinates axes $x,y,z$.
The first-order Markov parameters of discretized stochastic acceleration models in
$ x,y,$ and $ z $ axes are denoted by $ ( 3 \times 3 ) $ matrix
$ {\Lambda } $ = diag. $[ \lambda _ x , \lambda  _ y , \lambda _ z ]$,
where $ \lambda  = (1 - e ^ {- \alpha \Delta t}).  $
The matrix $ \Phi ( \Delta t, {\Lambda } )$ is given by
\begin{equation}
\Phi ( \Delta t, {\Lambda } ) = \left [
\begin{array}{c c c}
    \Phi _x & 0 & 0  \\
    0 &  \Phi _y       & 0 \\
    0 &  0          &  {\Phi _z }
\end{array}
\right ]
\label{eqnphi}
\end{equation}
where $ \Phi _x $ is $ (3 \times 3) $ matrix given by
\begin{displaymath}
\Phi _x = \left [
\begin{array} {c c c}
1  & \Delta t  & 1/2 (\Delta t) ^ 2 \\
0  &  1        & \Delta t \\
0  &  0        & \lambda _x
\end{array}
\right ]
\end{displaymath}
The symbols $ \Phi _x,\Phi _y, \Phi _z $ denote the state transition matrices
for $x,y$ and $z$ axes respectively.

The state vector
$ {\bf x}(t) $ is represented by

\begin{displaymath}
{\bf x}(t)= [ x, \dot x,\ddot x,y,\dot y,\ddot y ,z,\dot z,
\ddot z ] ^ T
\end{displaymath}
where $ T $ is a transpose operation.  The vector $ {\bf h(x}(t)) $
is the transformation from the Cartesian coordinates to the polar
coordinates defined by   
\begin{eqnarray}
 r & = & \sqrt{x ^ 2 + y ^ 2 + z ^ 2} \nonumber \\
 \theta _ T & = & \arctan (y/x) \nonumber \\
 \theta _ E & = & \arcsin (z/r) \nonumber
\end{eqnarray}
where $ r, \theta _ T, $ and $ \theta _ E $ denote the target range,
azimuth angle, and elevation angle, respectively.  The noisy 
measurements of range, azimuth, and elevation angles are denoted by
$ {\bf y}(t) $.
\begin{displaymath}
{\bf y}(t) = [r, \theta _ T, \theta _ E ] ^ T.
\end{displaymath}
It is assumed that the state and measurement noise sequences $ \{ {\bf w}(t) \}$ and $\{ {\bf v}(t) \} $ are zero mean gaussian random vector sequences
with covariance matrices $ Q (t) $ and $ R (t) $ respectively.

The  measurement errors in the target range, azimuth angle, and
elevation angle are assumed uncorrelated and may depend on the target
range.  Let $ \sigma ^ 2 _ R , \sigma ^ 2 _ {\theta T }, $ and
$ \sigma ^ 2 _ {\theta E} $ denote the variance of the measurement
noise in range, azimuth angle, and elevation angle, respectively.  The
measurement noise covariance matrix $ R(t) $ is given by
\begin{displaymath}
R(t)= diag [ \sigma ^ 2 _ R, \sigma ^ 2 _ {\theta T}, \sigma ^ 2 _ {\theta E} ].
\end{displaymath}
The dynamics and measurement models presented above are commonly used in
many tracking and navigation applications.  An extended Kalman filter
can be used to recursively estimate the state vector $ {\bf x}(t) $ 
using the dynamics and the measurement model. 
The  extended Kalman filter uses linearization
of nonlinear equations about the current state estimates.  The structure
of the covariance and gain computation equations in the extended Kalman
filter is similar to the structure of the Kalman filter.  However, it is 
important to note that the extended Kalman filter depends on the current state
estimate and thus cannot be computed off-line. 

A block diagram of a target tracking extended Kalman filter is given
in Figure \ref{kfblock}.  The measurements ${\bf y}(t) $ include the effect of sensor
noise, reference input, disturbance, and the target dynamics.  The 
filter consists of time update, Jacobian, and measurement update equations
summarized in the following.

\begin{figure}
\setlength{\unitlength}{0.155in}
\begin{center}
\begin{picture}(19,16)(0,0)
\scriptsize
\put(7,11){\framebox(5,2){\shortstack{Target\\Dynamics}}}
\put(2,7){\framebox(5,2){\shortstack{Time\\Update}}}
\put(12,7){\framebox(5,2){\shortstack{Measurement\\Update}}}
\put(12,4){\framebox(5,1){Jacobian}}
\put(8,6){\framebox(3,1){delay}}
\put(8,1){\framebox(3,1){delay}}
\put(14.5,12){\circle{1}}
\put(14.5,12){\makebox(0,0){+}}
% main loop
\put(1,12){\vector(1,0){6}}
\put(4,12.2){\makebox(0,0)[b]{Reference Input}}
\put(12,12){\vector(1,0){2}}
\put(15,12){\line(1,0){3}}
\put(17.8,10.5){\makebox(0,0)[r]{Measurement}}
\put(18,12){\line(0,-1){3.5}}
\put(18,8.5){\vector(-1,0){1}}
\put(12,8){\vector(-1,0){5}}
\put(9.5,8.2){\makebox(0,0)[b]{\shortstack{$\hat{x}(t|t)$\\$P(t|t)$}}}
\put(2,8){\line(-1,0){1}}
\put(1,8){\line(0,-1){6.5}}
\put(1,1.5){\vector(1,0){7}}
\put(4.5,1.7){\makebox(0,0)[b]{\shortstack{$\hat{x}(t+1|t)$\\$P(t+1|t)$}}}
\put(11,1.5){\line(1,0){7}}
\put(14.5,1.7){\makebox(0,0)[b]{\shortstack{$\hat{x}(t|t-1)$\\$P(t|t-1)$}}}
\put(18,1.5){\line(0,1){6}}
\put(18,7.5){\vector(-1,0){1}}
%jacobian loop
\put(9.5,8){\vector(0,-1){1}}
\put(9.5,6){\line(0,-1){1.5}}
\put(9.3,5){\makebox(0,0)[r]{$\hat{x}(t-1|t-1)$}}
\put(9.5,4.5){\vector(1,0){2.5}}
\put(14.5,5){\vector(0,1){2}}
\put(15,6){\makebox(0,0)[l]{$H$}}
%disturbance and sensor noise
\put(9.5,14){\vector(0,-1){1}}
\put(9.5,14.2){\makebox(0,0)[b]{Disturbance}}
\put(14.5,14){\vector(0,-1){1.5}}
\put(14.5,14.2){\makebox(0,0)[b]{Sensor Noise}}
\end{picture}
\end{center}
\caption{Tracking Filter}
\label{kfblock}
\end{figure}

\subsection{Time Update}

The time-update predicts the state and estimation error covariance matrix using the
following equations:
\begin{eqnarray}
{\bf {\hat{x}}}(t+1|t) & = & \Phi {\bf \hat{x}}(t|t) \label{eqntux} \\
P(t+1|t) & = & \Phi P(t|t) {\Phi} ^ T + Q(t) \label{eqntup}
\end{eqnarray}
where $ {\bf \hat{x}}(t+1|t) $ denotes one-step state prediction
at time $ (t+1) $
based on  the measurements at time $t$. The state estimate at time $t$
denoted by ${\bf {\hat{x}}}(t|t) $ represents the estimate based on the
measurements at time $t$. The state transition matrix $ \Phi $ 
is defined in Equation (\ref{eqnphi}) . The state estimation error covariance
matrix is given by
$ P(t|t) = E \{ {{\bf \hat{x}}(t|t)}{{\bf \hat{x}} ^ T (t|t) } \} $.
The state prediction covariance matrix is defined as $ P(t+1|t) = 
E \{ {{\bf \hat{x}}(t+1|t)}{{\bf \hat{x}} ^ T (t+1|t)} \} $.  The 
time update of the covariance matrix $P(t+1|t)$ accounts
for modeling errors and disturbances by a suitable choice of $Q(t)$
in Equation (\ref{eqntup}) . 


\subsection{Jacobian Update}

The measurement model in Equation (\ref{model}) is nonlinear and is 
linearized  at the optimal state estimate $ {\bf \hat{x}}(t-1|t-1) $ at
time $t-1$ so that it can be used at $t$ in the measurement update. 
Let $H$ denote the Jacobian matrix of ${\bf h(x)}$ evaluated at 
${\bf \hat{x}}(t-1|t-1)$.
\begin{eqnarray}
H & = & \frac { \partial {\bf h(x)}}  { \partial {\bf x}} 
|  {{\bf{x}}  = {{\bf \hat{x}}(t-1|t-1)}} \label{eqnjac}
\end{eqnarray}

\subsection{Measurement Update}

The  measurement update corrects the predicted state and error
covariance using the measurement ${\bf y}(t)$ and the measurement noise
covariance matrix $R(t)$ as follows:
The prediction error known as {\em innovations}
is denoted by ${\bf e}(t)$ and defined as
\begin{displaymath}
{\bf e}(t) = {\bf y}(t)-{\bf h}({\bf \hat{x}}(t|t-1)
\end{displaymath}
The state estimate is given as
\begin{equation}
{\bf \hat{x}}(t|t)={\bf \hat{x}}(t|t-1)+K(t){\bf e}(t) \label{eqnmux}
\end{equation}
The covariance of the innovations is denoted by $R_e(t)$ and defined by
\begin{equation}
R_e(t)=HP(t|t-1)H ^ T +R(t) \label{eqninnovcov}
\end{equation}
The Kalman gain matrix $K(t)$ and the covariance matrix $P(t|t)$ are
given by
\begin{eqnarray}
K(t) & = & P(t|t-1) H^T  { R_e(t)} ^ {-1} \label{eqngain} \\
P(t|t) & = & P(t|t-1) - K(t)HP(t|t-1) \label{eqnmup}
\end{eqnarray}
The initial conditions to start the recursive computations are given by
$ {\bf \hat{x}}(0|0) = {\bf x}(0) $ and $P(0|0)=P(0)$, the {\em a priori}
estimates.  In practice, the first few measurements are used to determine
the initial conditions. 


\begin{figure}
\setlength{\unitlength}{0.14in}
\begin{center}
\begin{picture}(25,10)(1,0)
\scriptsize
\put(4,6){\framebox(5,3){\shortstack{Measurement\\Update\\(range)}}}
\put(11,6){\framebox(5,3){\shortstack{Measurement\\Update\\(azimuth)}}}
\put(18,6){\framebox(5,3){\shortstack{Measurement\\Update\\(elevation)}}}
\put(11,2){\framebox(5,2){Jacobian}}
% x
\multiput(2,8)(7,0){4}{\vector(1,0){2}}
\put(2.3,8.3){\makebox(0,0)[b]{$\hat{x}(t|t-1)$}}
\put(10,8.3){\makebox(0,0)[b]{$\hat{x}^{(1)}$}}
\put(17,8.3){\makebox(0,0)[b]{$\hat{x}^{(2)}$}}
\put(24,8.3){\makebox(0,0)[b]{$\hat{x}(t|t)$}}
% p
\multiput(2,7)(7,0){4}{\vector(1,0){2}}
\put(2.3,6.4){\makebox(0,0){$P(t|t-1)$}}
\put(10,6.4){\makebox(0,0){$P^{(1)}$}}
\put(17,6.4){\makebox(0,0){$P^{(2)}$}}
\put(24,6.4){\makebox(0,0){$P(t|t)$}}
% jac
\put(6.5,3){\vector(0,1){3}}
\put(13.5,4){\vector(0,1){2}}
\put(20.5,3){\vector(0,1){3}}
\put(6.5,3){\line(1,0){4.5}}
\put(16,3){\line(1,0){4.5}}
\put(13.5,1){\vector(0,1){1}}
\put(13.5,0){\makebox(0,0)[b]{$\hat{x}(t-1|t-1)$}}
\put(7,5){\makebox(0,0)[l]{$H_1$}}
\put(14,5){\makebox(0,0)[l]{$H_2$}}
\put(21,5){\makebox(0,0)[l]{$H_3$}}
\end{picture}
\end{center}
\caption{Measurement Update}
\label{measblock}
\end{figure}


\subsection{Scalar Measurement Processing}
In the measurement update, Equation (\ref{eqngain}) requires the inverse
of an $(m \times m)$  matrix that can be avoided by processing 
each measurement individually \cite{bierman77}.  Let $y_i(t)$ and
$h_i({\bf \hat{x}})$ denote the $ith$ element of vectors ${\bf {y}}(t)$ and
${\bf h(\hat{x})}$ respectively. Let $H_i$ denote the $ith$ row of
the measurement partial derivative matrix $H$. Then Equation (\ref{eqnmux})
can be written as
\begin{eqnarray}
{\bf \hat{x}}^{(0)} & = & {\bf \hat{x}}(t|t-1) \nonumber \\
e ^{(i)} &      = & y_{i}(t) - h_{i}({\bf \hat{x}}(t|t-1)) \nonumber \\
{\bf \hat{x}}^{(i)} & = & {\bf \hat{x}}^{(i-1)} + K^{(i)}e^{(i)}, \label{eqnmuxk}\\
{\bf \hat{x}}(t|t) & = & {\bf \hat{x}}^{(m)} \nonumber
\end{eqnarray}
Equation (\ref{eqnmup}) can be written as 
\begin{eqnarray}
P^{(0)} & = & P(t|t-1) \nonumber \\ 
P^{(i)} & = & P^{(i-1)} - K^{(i)}H_{i}P^{(i-1)}, ~i = 1,...,m \label{eqnmupk}\\
P(t|t) & = & P^{(m)} \nonumber
\end{eqnarray}
Equation (\ref{eqninnovcov}) can be written as 
\begin{displaymath}
R_e^{(i)} = H_iP^{(i-1)}H_i^T+R_i(t)
\end{displaymath}
And Equation (\ref{eqngain}) can be written as 
\begin{eqnarray}
K^{(i)} &=& P^{(i-1)}H_i^T (R_e^{(i)})^{-1} \label{eqngaink} 
\end{eqnarray}			
The computations for the range, azimuth angle, and elevation angle
measurements using the scalar measurement algorithm for target
tracking are summarized in Figure \ref{measblock}.  A block diagram of
the $ith$ scalar update is given in Figure \ref{scalarblock}.
Note that $K^{(i)}$ in equation (\ref{eqngaink}) is not the $ i-th $
column of matrix $K$ in equation (\ref{eqngain})
\begin{figure}
\setlength{\unitlength}{0.155in}
\begin{center}
\begin{picture}(20,13)(1,0)
\scriptsize
\put(1,8){\framebox(5,3){\shortstack{Jacobian\\(eq. \ref{eqnjac})}}}
\put(8,8){\framebox(5,3){\shortstack{Innovation\\Covariance\\(eq. \ref{eqngaink})}}}
\put(15,8){\framebox(5,3){\shortstack{Innovation\\Error\\(eq. \ref{eqnmuxk})}}}
\put(1,3){\framebox(5,3){\shortstack{Covariance\\Update\\(eq. \ref{eqnmupk})}}}
\put(8,3){\framebox(5,3){\shortstack{Gain\\Computation\\(eq. \ref{eqngaink})}}}
\put(15,3){\framebox(5,3){\shortstack{State\\Update\\(eq. \ref{eqnmuxk})}}}
%top
\put(3.5,12){\vector(0,-1){1}}
\put(3.5,12.2){\makebox(0,0)[b]{$\hat{x}(t-1|t-1)$}}
\put(9,12){\vector(0,-1){1}}
\put(9,12.2){\makebox(0,0)[b]{$y_{i}(t)$}}
\put(12,12){\vector(0,-1){1}}
\put(12,12.2){\makebox(0,0)[b]{$P^{(i-1)}$}}
\put(17.5,12){\vector(0,-1){1}}
\put(17.5,12.2){\makebox(0,0)[b]{$\hat{x}^{(i-1)}$}}
%middle
\put(6,9.5){\vector(1,0){2}}
\put(7,9.7){\makebox(0,0)[b]{$H_i$}}
\put(10.5,8){\vector(0,-1){2}}
\put(11,7){\makebox(0,0)[l]{$R_e^{(i)}$}}
\put(17.5,8){\vector(0,-1){2}}
\put(18,7){\makebox(0,0)[l]{$e^{(i)}$}}
\put(8,4.5){\vector(-1,0){2}}
\put(7,4.7){\makebox(0,0)[b]{$K^{(i)}$}}
\put(13,4.5){\vector(1,0){2}}
\put(14,4.7){\makebox(0,0)[b]{$K^{(i)}$}}
%bottom
\put(3.5,3){\vector(0,-1){1}}
\put(3.5,1.8){\makebox(0,0)[t]{$P^{(i)}$}}
\put(17.5,3){\vector(0,-1){1}}
\put(17.5,1.8){\makebox(0,0)[t]{$\hat{x}^{(i)}$}}
\end{picture}
\end{center}
\caption{$ith$ Scalar Measurement Update}
\label{scalarblock}
\end{figure}


\section{Mapping on Linear Array}
The Kalman equations can be characterized at a course level by the
{\em computation graph} in Figure \ref{compgraph}.  For our purposes, a
computation graph is a directed acyclic graph consisting of a starting
vertex, $s$, which has no input edges, an ending vertex, $e$, which
has no output edges, and a collection of computation vertices, which
represent a partial ordering of computations.  Vertices are labeled
with types of computations; edges are labeled with input data items.
An edge $(u,v)$ labeled with $A$ indicates that $u$ produces $A$ and
$v$ consumes $A$. Edge $(u,v)$ also represents a weak precedence
relation between $u$ and $v$ in that $v$ is guaranteed not to complete
its computation until $u$ has produced all the elements of $A$.

\begin{figure}
\setlength{\unitlength}{0.15in}
\begin{center}
\begin{picture}(25,45)(0,0)
\scriptsize
%nodes
\put(15,45){\makebox(0,0){$\bullet$}}
\put(14.55,46){\fbox{s}}

\put(10,40){\makebox(0,0){$\bullet$}}
\put(7.5,39){\fbox{\em tux}}

\put(15,40){\makebox(0,0){$\bullet$}}
\put(12.75,40){\fbox{\em jac}}

\put(20,40){\makebox(0,0){$\bullet$}}
\put(20.5,39){\fbox{\em tup}}

\put(5,35){\makebox(0,0){$\bullet$}}
\put(5.5,35){\fbox{\em trans}}

\put(15,35){\makebox(0,0){$\bullet$}}
\put(12,35){\fbox{\em gain}}

\put(5,30){\makebox(0,0){$\bullet$}}
\put(3.75,30){\fbox{s}}

\put(10,30){\makebox(0,0){$\bullet$}}
\put(7.25,29){\fbox{\em mux}}

\put(20,30){\makebox(0,0){$\bullet$}}
\put(20.5,29){\fbox{\em mup}}
\put(5,25){\makebox(0,0){$\bullet$}}
\put(5.5,25){\fbox{\em trans}}

\put(15,25){\makebox(0,0){$\bullet$}}
\put(12,25){\fbox{\em gain}}

\put(5,20){\makebox(0,0){$\bullet$}}
\put(3.75,20){\fbox{s}}

\put(10,20){\makebox(0,0){$\bullet$}}
\put(7.25,19){\fbox{\em mux}}

\put(20,20){\makebox(0,0){$\bullet$}}
\put(20.5,19){\fbox{\em mup}}

\put(5,15){\makebox(0,0){$\bullet$}}
\put(5.5,15){\fbox{\em trans}}

\put(15,15){\makebox(0,0){$\bullet$}}
\put(12,15){\fbox{\em gain}}

\put(5,10){\makebox(0,0){$\bullet$}}
\put(3.75,10){\fbox{s}}

\put(10,10){\makebox(0,0){$\bullet$}}
\put(7.25,9){\fbox{\em mux}}

\put(20,10){\makebox(0,0){$\bullet$}}
\put(20.5,9){\fbox{\em mup}}

\put(15,5){\makebox(0,0){$\bullet$}}
\put(13.5,3.75){\fbox{\em e(xout)}}
%circles
\put(5,40){\circle{.75}}
\put(5,40){\makebox(0,0){3}}

\put(13,38){\circle{.75}}
\put(13,38){\makebox(0,0){2}}

\put(17,38){\circle{.75}}
\put(17,38){\makebox(0,0){1}}

\put(5,27.5){\circle{.75}}
\put(5,27.5){\makebox(0,0){3}}

\put(15,30){\circle{.75}}
\put(15,30){\makebox(0,0){1}}

\put(5,17.5){\circle{.75}}
\put(5,17.5){\makebox(0,0){3}}

\put(15,20){\circle{.75}}
\put(15,20){\makebox(0,0){2}}

%labeled edges
\put(15,45){\vector(-1,-1){5}}
\put(10.25,43){\makebox(0,0){$\hat{x}(t-1|t-1)$}}

\put(15,45){\vector(0,-1){5}}
\put(15.5,42){\makebox(0,0){$\hat{x}$}}

\thicklines
\put(15,45){\vector(1,-1){5}}
\thicklines
\put(15,45){\vector(1,-1){5}}
\thinlines
\put(19.75,43){\makebox(0,0){$P(t-1|t-1)$}}

\put(10,40){\line(-1,0){4.5}}
\put(7.5,40.5){\makebox(0,0){$\hat{x}^{(0)}$}}

\put(5,39.5){\vector(0,-1){4.5}}
\put(15,40){\vector(-1,-1){1.75}}

\put(15,40){\vector(0,-1){5}}
\put(14.5,37){\makebox(0,0){$H_1$}}

\put(15,40){\vector(1,-2){5}}
\put(18.5,34.5){\makebox(0,0){$H_1$}}
\put(15,40){\vector(1,-1){1.75}}

\multiput(20,40)(0,-10){3}{\vector(-1,-1){5}}
\thicklines
\put(20,10){\vector(-1,-1){5}}
\thinlines
\put(18,37){\makebox(0,0){$P^{(0)}$}}
\put(18,27){\makebox(0,0){$P^{(1)}$}}
\put(18,17){\makebox(0,0){$P^{(2)}$}}
\put(18.35,7){\makebox(0,0){$P(t|t)$}}

\thicklines
\multiput(20,40)(0,-10){3}{\vector(0,-1){10}}
\thinlines
\put(21,35){\makebox(0,0){$P^{(0)}$}}
\put(21,25){\makebox(0,0){$P^{(1)}$}}
\put(21,15){\makebox(0,0){$P^{(2)}$}}
\multiput(10,40)(0,-10){3}{\vector(0,-1){10}}
\put(10.85,35){\makebox(0,0){$\hat{x}^{(0)}$}}
\put(10.85,25){\makebox(0,0){$\hat{x}^{(1)}$}}
\put(10.85,15){\makebox(0,0){$\hat{x}^{(2)}$}}

\multiput(5,35)(0,-10){3}{\vector(1,-1){5}}
\put(8,33){\makebox(0,0){$h_1$}}
\put(8,23){\makebox(0,0){$h_2$}}
\put(8,13){\makebox(0,0){$h_3$}}

\multiput(5,30)(0,-10){3}{\vector(1,0){5}}
\put(7.5,30.5){\makebox(0,0){$y_1$}}
\put(7.5,20.5){\makebox(0,0){$y_2$}}
\put(7.5,10.5){\makebox(0,0){$y_3$}}

\multiput(15,35)(0,-10){3}{\vector(-1,-1){5}}
\put(13,32){\makebox(0,0){$K^{(1)}$}}
\put(13,22){\makebox(0,0){$K^{(2)}$}}
\put(13,12){\makebox(0,0){$K^{(3)}$}}

\multiput(15,35)(0,-10){3}{\vector(1,-1){5}}
\put(16.75,32){\makebox(0,0){$K^{(1)}$}}
\put(16.75,22){\makebox(0,0){$K^{(2)}$}}
\put(16.75,12){\makebox(0,0){$K^{(3)}$}}

\multiput(5,27)(0,-10){2}{\vector(0,-1){2}}

\multiput(15,29.5)(0,-10){2}{\vector(0,-1){4.6}}
\put(14,27.5){\makebox(0,0){$H_2$}}
\put(14,17.5){\makebox(0,0){$H_3$}}

\multiput(15.15,29.65)(0,-10){2}{\vector(1,-2){4.85}}
\put(18.5,24.25){\makebox(0,0){$H_2$}}
\put(18.5,14.25){\makebox(0,0){$H_3$}}

\put(10,10){\vector(1,-1){5}}
\put(11,7){\makebox(0,0){$\hat{x}(t|t)$}}

\end{picture}

\scriptsize
\begin{tabular}{|l|l|l|l|}
\hline
Vertex & Description & Eqn \# & Complexity \\
\hline
\hline
gain & Kalman gain computation & \ref{eqngaink} & $O(n)$ \\
jac & Jacobian transform & \ref{eqnjac} & $O(m)$ \\
mup & scalar measurement update of $P$ & \ref{eqnmupk} & $O(n^2)$ \\
mux & scalar measurement update of $\hat{x}$ & \ref{eqnmuxk} & $O(n)$ \\
trans & coordinate transform & \ref{eqnmuxk} & $O(m)$ \\
tup & time update of $P$ & \ref{eqntup} & $O(n^3)$ \\
tux & time update of $\hat{x}$ & \ref{eqntux} & $O(n^{2})$ \\
\hline
\end{tabular}
\vspace{.25in}
\end{center}
\caption{Kalman Filter Computation Graph}
\label{compgraph}
\end{figure}

The problem to be solved is the design of a fast parallel Kalman
filtering algorithm for a linear array of computational cells.
Inspection of Figure \ref{compgraph} shows that the Kalman filter
computations are an eclectic mix of matrix operations ranging in
complexity from linear to cubic in the number of states.  But linear
arrays have traditionally been applied to regular and uniform
computations.  Thus, in a more general sense, the problem to be solved
is the mapping of an irregular and non-uniform computation onto a
linear array.
 
It is easy to see that the covariance updates 
(Equations \ref{eqntup}) and (\ref{eqnmupk})) require most
of the computations.  This motivates the first and rather obvious
design goal:
\begin{quote}
{\em Focus on covariance matrix updates.}
\end{quote}
Along with their large computational requirements, the covariance
updates in Figure \ref{compgraph} have the largest input and output
requirements. This motivates the second design goal:
\begin{quote}
{\em Minimize data transfers involving covariance.}
\end{quote}
On systems like the Warp computer, the data transfer rate between the
external host and the array of cells is an order of magnitude less
than the transfer rate between cells.  This fact, along with the fact
that many of the inputs in Figure \ref{compgraph} are small vectors
and matrices motivates the third design goal:
\begin{quote}
{\em Minimize data transfers between host and array.}
\end{quote}
The next section describes a parallel algorithm for the Kalman
filter that achieves these goals.


A high level diagram of a linear array with $n$ computational cells is given 
in Figure \ref{arch}. Boxes represent computational cells.
Arrows represent unidirectional buffered communication channels. Data
flows from left to right along the upper channel and from right to
left along the lower channel.
A parallel algorithm to implement an $n$-state Kalman
filter on an $(n+1)$-cell linear array is given in the following.   The
9-state target tracking filter is used as an example to illustrate
the concepts.

\paragraph{Covariance Time Update:}

The time update of the covariance matrix given by equation (\ref{eqntup})
is accomplished by updating one column of the covariance matrix on a 
single computational cell.  Thus, for an $(n \times n)$ covariance
matrix update $n$ cells are used.  The time update requires $O(n^3)$
computations and by using $n$ processing elements the computations
can be done in $O(n^2)$.

Equation (\ref{eqntup}) can be written as
\begin{equation}
p_{ij}(t+1|t) = \sum_{k=1}^{n} \sum_{\zeta = 1}^{n} p_{k\zeta}(t|t)
\phi_{ik}\phi_{j\zeta}+q_{ij}(t)
\label{eqnmap2}
\end{equation}
where $p_{ij},\phi_{ij}$ and $q_{ij}$ denote the $ijth$ elements in
$P,\Phi,Q$ respectively.

\begin{figure}
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(20,5)(1,0)
\scriptsize
\put(5,1){\framebox(2,2){\shortstack{{\em cell}\\{\em 1}}}}
\put(8,1){\framebox(2,2){\shortstack{{\em cell}\\{\em 2}}}}
\put(11,1){\framebox(2,2){\shortstack{{\em cell}\\{\em 3}}}}
\put(14,1){\makebox(2,2){$\cdots$}}
\put(17,1){\framebox(2,2){\shortstack{{\em cell}\\{\em n}}}}
\put(3,2.5){\vector(1,0){2}}
\put(5,1.5){\vector(-1,0){2}}
%\put(4,3.1){\makebox(0,0){$y(t)$}}
%\put(4,1.1){\makebox(0,0){$\hat{x}(t|t)$}}
\multiput(7,2.5)(3,0){4}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){4}{\vector(-1,0){1}}
\end{picture}
\end{center}
\caption{Linear Array of Computational Cells}
\label{arch}
\end{figure}

For the target tracking filter, the state transition matrix given by
equation (\ref{eqnphi}) has 18 non-zero elements.
The computations can be minimized by simplifications.  For example,
$p_{12}(t+1|t)$ is given by
\begin{eqnarray}
p_{12} &=& p_{12} + p_{22}\delta + p_{32}\sigma + p_{13}\delta + 
p_{23}\delta^2 + p_{33}\delta\sigma  \nonumber \\
& & +q_{12} \nonumber
\end{eqnarray}
where $\delta = \Delta t$ and $\sigma = 1/2(\Delta t)^2$.  


The parallel algorithm for the time update is summarized in the
following.  Before
the time update, the $kth$ cell contains the $kth$ column of $P$.
\begin{center}
\begin{tabular}{rccccccccc}
cell:&1&2&3&4&5&6&7&8&9\\ 
Start:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$\\ 
\end{tabular}
\end{center}
First, the columns  $P_2,P_3,P_5,P_6,P_8,P_9$ are shifted to the left cells and
the results are stored in appropriate cells.
\begin{center}
\begin{tabular}{rccccccccc}
cell:&1&2&3&4&5&6&7&8&9\\ 
Start:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
Shift-left:&$P_2$&$P_3$&     &$P_5$&$P_6$&     &$P_8$&$P_9$&      \\
\end{tabular}
\end{center}
Next, the columns  $P_3,P_6,P_9$ are again shifted to the left cells and
the results are stored.
\begin{center}
\begin{tabular}{rccccccccc}
cell:&1&2&3&4&5&6&7&8&9\\ 
Start:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
Shift-left:&$P_2$&$P_3$&     &$P_5$&$P_6$&     &$P_8$&$P_9$&      \\
Shift-left:&$P_3$&     &     &$P_6$ &    &    &$P_9$  &    &      \\
\end{tabular}
\end{center}

No further shift operations are required due to the sparse structure
of the state transition matrix.  
With the appropriate columns in place, each cell updates its column of
$P$ independently to compute the time update of the covariance matrix.


Compare the algorithm in this section with a known systolic algorithm
for performing matrix multiplication on $n$ linearly-connected
cells\cite{annaratone87e}. The algorithm in \cite{annaratone87e} works
for arbitrary matrices, requires $2n^{3} = 1458$ scalar
multiplications to perform the pair of matrix multiplications in
Equation (\ref{eqntup}), and does not exploit the constant and sparse
nature of $\Phi$.
The algorithm described above, on the other hand, is appropriate for
filters with sparse, state transition
matrices.  For the tracking filter,
the time update requires only 288 scalar multiplications. Further,
the elements of $P$ that need to be transferred to different cells is
minimal.

\paragraph{Covariance Scalar Measurement Update:}
The $ith$ scalar measurement update of $P$ (Equation \ref{eqnmup}))
can be rewritten in terms of the $kth$ column of $P$.
\begin{equation}
P_k^{(i)}  =  P_k^{(i-1)} - K^{(i)}H_{i}P_k^{(i-1)}, ~~k = 1,...,n
\label{eqnmupcol}
\end{equation}
Since $H_i$ is a vector, $H_{i}P$ and $PH_{i}^{T}$ are vectors, and since
$P$ is symmetric,
\begin{equation}
(H_{i}P)^T = PH_{i}^T.
\label{eqnmupequal}
\end{equation}
Since a vector stored in computer memory is identical to its
transpose, Equation (\ref{eqnmupequal}) implies that we can compute
$H_{i}P$ once and then use the resulting vector wherever $H_{i}P$ and
$PH_{i}^T$ appear in Equations (\ref{eqnmup}) and (\ref{eqngaink}).
It is easy to see that 
\begin{equation}
H_{i}P = \left[ H_{i}P_{1} ~ H_{i}P_{2} ~\cdots~ H_{i}P_{n} \right] 
\label{eqnhp}
\end{equation}
Equations (\ref{eqnmupcol}), (\ref{eqnmupequal}), and (\ref{eqnhp}) suggest
the parallel algorithm for computing the $ith$ scalar measurement
update of $P$.  Before the $ith$ scalar measurement update of $P$, the
$kth$ cell contains the $kth$ column of $P$:
\begin{center}
\begin{tabular}{rccccccccc}
cell:&1&2&3&4&5&6&7&8&9\\ 
Start:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
\end{tabular}
\end{center}
The elements of $H_i$ are sent from the left to the array of cells.
\begin{displaymath}
\cdots ~H_{i2} ~H_{i1} \rightarrow \framebox{1} \rightarrow \framebox{2} 
\rightarrow \cdots \rightarrow \framebox{9}
\end{displaymath}
As $H_i$ flows past the $kth$  cell, the cell accumulates
the $kth$ scalar element of $H_i P$. As scalar elements of $H_i
P$ are computed, they are passed back to the left and out of the array.
\begin{displaymath}
H_{i}P_1 ~H_{i}P_2 ~\cdots \leftarrow \framebox{1} \leftarrow
\framebox{2}
\leftarrow \cdots \leftarrow \framebox{9}
\end{displaymath}
When the elements of $H_i P$ have been computed, the array of cells
is in the following state.
\begin{center}
\begin{tabular}{rccccccccc}
cell:&1&2&3&4&5&6&7&8&9\\ 
Start:&$P_1$&$P_2$&$P_3$&$P_4$&$P_5$&$P_6$&$P_7$&$P_8$&$P_9$ \\
Compute:&$a_1$&$a_2$&$a_3$&$a_4$&$a_5$&$a_6$&$a_7$&$a_8$&$a_9$ \\
\end{tabular}
\end{center}
where $a_k = H_i P_k$. Next, $H_i P$ is used to
compute the gain vector $K^{(i)}$ (Equation (\ref{eqngaink})) As each
element of $K^{(i)}$ is computed it is sent from the left to the array of 
cells.
\begin{displaymath}
\cdots ~K_2^{(i)} ~K_1^{(i)} \rightarrow \framebox{1} \rightarrow \framebox{2} 
\rightarrow \cdots \rightarrow \framebox{9}
\end{displaymath}
As $K^{(i)}$ flows through the array, the $kth$ cell uses
$K^{(i)}$ and $a_k = H_{i} P_k$ to perform the scalar measurement
update of $P_k$ according to Equation (\ref{eqnmupcol}). 
Note that the $ith$ scalar measurement update requires $2n$ scalar multiplications on each cell. 

It is seen that $n$ cells are used to update the time and measurement
covariance matrix. The Jacobian transformation,
state update and coordinate transformations used in the target
tracking filter are computed sequentially. A separate cell on the
left that communicates with the $n$-cells can be used for these computations.
\cite{ohallaron88d}. 
\section{Mapping on Warp} 
The Warp computer \cite{annaratone87e} (or simply Warp) is a linear
array of  identical and programmable cells connected to a
general-purpose host. Figure \ref{warp} shows the architecture of
a 10-cell Warp.  A Sun-3 workstation called the {\em master}
is connected via a VME repeater to a pair of {\em clusters}. Each
cluster consists of an MC68020 processor(P) with an MC68881 floating
point coprocessor, 3 Mbytes of data memory(M), 1 Mbyte of program
memory, and a switch(S), all connected by a local VSB bus. The
clusters are known collectively as the external host, so called
because they are external to the Master.

\begin{figure}
\setlength{\unitlength}{0.15in}
\begin{center}
\begin{picture}(30,20)(3,2)
\small
%host
%\multiput(4,25)(6,0){4}{\framebox(4,2){Sun-3}}
\put(13,19){\framebox(4,2){\shortstack{Sun-3\\Master}}}
%\put(6,23){\line(1,0){18}}
%\multiput(6,23)(6,0){4}{\line(0,1){2}}
%\put(15,21){\line(0,1){2}}
\put(15,16){\line(0,1){3}}
%external host
\put(9,16){\line(1,0){12}}
\put(9,13){\framebox(1,2){P}}
\put(11,13){\framebox(1,2){M}}
\put(13,11){\framebox(1,4){S}}
\put(20,13){\framebox(1,2){P}}
\put(18,13){\framebox(1,2){M}}
\put(16,11){\framebox(1,4){S}}
\multiput(9.5,15)(2,0){3}{\line(0,1){1}}
\multiput(16.5,15)(2,0){3}{\line(0,1){1}}
\multiput(9.5,12)(2,0){2}{\line(0,1){1}}
\multiput(18.5,12)(2,0){2}{\line(0,1){1}}
\put(9,12){\line(1,0){4}}
\put(17,12){\line(1,0){4}}
\put(13.5,9){\line(0,1){2}}
\put(16.5,9){\vector(0,1){2}}

%iu and array
%\put(12,8){\framebox(6,1){IU}}
\put(10,4){\framebox(2,2){\shortstack{\em Cell\\0}}}
\put(13,4){\framebox(2,2){\shortstack{\em Cell\\1}}}
\put(15.5,5){$\ldots$}
\put(18,4){\framebox(2,2){\shortstack{\em Cell\\9}}}
%from IU to array
\put(9,9){\line(1,0){4.5}}
\put(21,9){\line(-1,0){4.5}}
\put(9,5.25){\line(0,1){3.75}}
\put(21,5.25){\line(0,1){3.75}}
%X channel X
\put(9,5.25){\vector(1,0){1}}
\put(12,5.25){\vector(1,0){1}}
\scriptsize
\put(12.25,5.5){X}
\small
\put(17,5.25){\vector(1,0){1}}
\put(20,5.25){\line(1,0){1}}
%Y channel Y
\put(12,4.75){\vector(1,0){1}}
\put(13,4.75){\vector(-1,0){1}}
\scriptsize
\put(12.25,4.0){Y}
\small
\put(17,4.75){\vector(1,0){1}}
\put(18,4.75){\vector(-1,0){1}}

%dashed boxes labeling clusters and array
\put(9.5,3.5){\dashbox{.1}(11,3)}
\put(8.5,10.5){\dashbox{.1}(6,6)}
\put(15.5,10.5){\dashbox{.1}(6,6)}
\put(11.5,17){\makebox(0,0){Input Cluster}}
\put(18.5,17){\makebox(0,0){Output Cluster}}
\put(15.3,2.5){\makebox(0,0){Warp Array}}
\end{picture}
\end{center}
\caption{The Warp Computer}
\label{warp}
\end{figure}


A Warp cell
contains a pipelined floating-point adder and multiplier, 8K
instruction words, and 32K data words. The clock cycle of each cell is
200ns. Each functional unit can emit one result per cycle, for a
maximum performance of 5 MFLOPS per functional unit, 10 MFLOPS per
cell, and 100 MFLOPS for the 10-cell Warp array. Each cell is
connected to its nearest neighbor by two 32-bit buffered data channels
(X and Y).  Data flows from left to right along the X channel; the Y
channel can be statically reconfigured to pass data in either
direction. Each cell can transfer up to 20 million 32-bit data words
to its neighboring cells each second.

Warp is programmed in W2, a high-level Algol-like language with {\em
send} and {\em receive} statements for cell/cell and host/array data
transfers.  The W2 compiler inputs a W2 source file, generates
microcode for the Warp array, and generates C code for the
clusters.  Warp is accessed with a programming environment, called
WPE, which provides a programmable shell for interactively running and
debugging W2 programs from any Sun connected via Ethernet to the
Master, as well as a set of lower-level routines that allow
user-written C and Lisp programs to access Warp.

The mapping of the 9-state tracking filter on Warp is 500 lines of W2
code. As shown in Figure \ref{mapping}, each {\em p-cell}, {\em
k-cell}, and {\em x-cell} is assigned a Warp cell.
A {\em p-cell} computes a column of the covariance matrix, {\em k-cell}
computes the gain vector and {\em x-cell} computes the state update.
During the initialization, {\em p-cell} 9 stores the 9th column of the
covariance matrix.  Similiarly {\em p-cell} 1 stores the first column.
Since the current
Warp has only 10 cells, we combine the {\em k-cell} and the {\em
x-cell} on Cell 9.

Data flows on the X channel from left to right; data flows on the Y
channel from right to left. In Figure \ref{arch} the {\em x-cell}
handles all data transfers between the host and the array. Since data
may enter the Warp array only from the left hand side and may exit the
array only from the right hand side, the mapping uses Cell
0 (corresponding to {\em p-cell}~9) to input the measurement vector
${\bf y}(t)$ after the time update of $P$ and to stream it along the X
channel through the array to Cell 9 (corresponding to the {\em
x-cell} and the {\em k-cell}). The updated state vector ${\bf \hat{x}}$
exits the array along the X channel.
\begin{figure}
\setlength{\unitlength}{0.175in}
\begin{center}
\begin{picture}(20,5)(1,0)
\scriptsize
\put(5,1){\framebox(2,2){\shortstack{\em Cell\\0}}}
\put(6,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{{\em p-cell}\\9}}$}}
\put(8,1){\makebox(2,2){$\cdots$}}
\put(11,1){\framebox(2,2){\shortstack{\em Cell\\8}}}
\put(12,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{{\em p-cell}\\1}}$}}
\put(14,1){\framebox(2,2){\shortstack{\em Cell\\9}}}
\put(15,4.25){\makebox(0,0){$\overbrace{}^{\shortstack{{\em k-cell}\\{\em x-cell}}}$}}
\put(3,2.5){\vector(1,0){2}}
\put(4,3.1){\makebox(0,0){$y(t)$}}
\put(16,2.5){\vector(1,0){2}}
\put(17,3.1){\makebox(0,0){$\hat{x}(t|t)$}}
\multiput(7,2.5)(3,0){3}{\vector(1,0){1}}
\multiput(8,1.5)(3,0){3}{\vector(-1,0){1}}
\put(7.5,2.8){X}
\put(7.5,0.9){Y}
\end{picture}
\end{center}
\caption{Mapping on Warp}
\label{mapping}
\end{figure}






The key performance measure for the Kalman filter is {\em sample
time}, that is, the real time required to process one measurement
sample.  To compare the sample time of the Warp filter with a
conventional sequential machine, we used the first 300 samples from a
typical maneuvering target trajectory.  
We ran a C version of the same Kalman filter,
with compiler optimization enabled on a Sun-3/60 with 68881
floating-point coprocessor. 


The C program does not access disk; all data are loaded with the
program and are available when the program begins execution.  The
sample time for the Sun was obtained using the {\em time} command of
the Unix shell. Total user CPU seconds was divided by a factor of 300
(the number of samples) to arrive at the sample time.  The sample time
for Warp was obtained using a built-in timing facility on the external
host. We noted the time that the input cluster started executing
before processing the first sample and we noted the time that the
output cluster finished executing after processing the last sample.
We divided the difference by 300 to arrive at the sample time. The
Warp sample time includes the time required to transfer data between
the external host and the Warp array, as well as startup times for the
external host and Warp array.

In the final implementation only 8 cells of the Warp have been used
by combining multiple columns of the covariance matrix. 
The sample times for Sun and Warp with 1-cell (sequential processing)
and Warp with 8-cells (parallel processing) 
are shown in Figure \ref{performance}.
A speedup is the ratio of 1-cell sample time to 8-cell sample
time.  It is seen that the target tracking extended-Kalman filter 
can be efficiently mapped on 8-cell Warp to achieve a speedup of 7.
There are a number of reasons for the two orders of magnitude speedup
of the Warp filter over the Sun filter. First, each of the Warp cells
is a fast scalar processor. Second, the Warp filter partitions the
covariance matrix updates over $n$ cells and performs these updates
in parallel. Third, the updates of the columns of the
covariance matrix are performed in-place, and except during the first
phase of the covariance time update, elements of the covariance matrix
are not transferred between cells. Finally, by performing all of the
Kalman filter computations on the Warp array, interaction between
the host and the array is minimized.

\begin{figure}
\begin{center}
\begin{tabular}{|c|r|}
\hline
\hline
 & Filter Sample Time\\
\hline
\hline
Sun-3/60 & 50000 $\mu$sec. \\
\hline
Warp (1-cell) & 875 $\mu$sec. \\
\hline
Warp (8-cell) & 125 $\mu$sec. \\
\hline
speedup       & 7 \\
\hline
\hline
\end{tabular}
\end{center}
\caption{Kalman Filter Performance}
\label{performance}
\end{figure}

\section{Conclusions}
We have developed a parallel algorithm for solving an $n$-state Kalman
filter on $n+1$ linearly connected computational cells and we have
implemented the algorithm on the Warp computer.  The algorithm
distributes the expensive covariance matrix updates across $n$ cells
and minimizes data flow through the array and interaction between the
external host and the array.
The Warp implementation
achieves a measured sample time of 125 $\mu$s, which represents a
speedup of 400 over the same Kalman filter running on a Sun
3/60 workstation. 
We have also implemented a square-root filtering algorithm \cite{bierman77}, 
where a factorized form of the covariance matrix is
updated to provide better numerical accuracy. 
However, the speedup
achieved was not as dramatic as the speedup reported here. 


In the current Warp design, each cell is a board. For on-board
applications, iWarp has been designed by Carnegie Mellon
University and Intel Corporation, where each cell computation
and communication function is integrated into a single VLSI
chip \cite{borkar88}.	Since the algorithms developed here
are applicable to the iWarp architecture, they open new solutions
to on-board target tracking applications.

\section{Acknowledgements}
This work was supported in part by Office of Naval Research Contract
N00014-86-K-0607 under CMU subcontract 30511-51189. Thanks to H.T.
Kung for graciously allowing us unlimited use of the Warp computer at
Carnegie Mellon, and to Chuck Fiduccia for a discussion that led to
the algorithm for the covariance time update.  Howard Itzkowitz
contributed very significantly in implementing the W2 code.
\bibliography{/usr/droh//bib/defs,/usr/droh/bib/refs,/usr/droh/bib/me}


\end{document}


