Real-time Linear Operator Construction and State Estimation with The Kalman Filter

1 Introduction

A quick tool for noise reduction and short-term prediction is important for areas such as weather forecasting and adjusting a scanning probe microscope (SPM). In weather forecasting, engineers require a fast denoising method in order to use the results for instantaneous forecasting. Since SPMs require adjustment of some parameters for use, a method for real-time adjustment is helpful in carrying out experiments.

A number of noise reduction methods have been proposed [2], [7]

. In recent years, denoising methods using deep learning, such as Noise2Noise and Noise2Void, have achieved great success

[21], [16], [14]. However, these methods cannot predict future states.

For predicting future states, a combination of convolutional neural networks (CNNs) and recurrent neural networks (RNNs) has been proposed since ConvLSTM

[18]. These methods use no model constraint and are applied to image sequence data [1], [17], [19], [23]. They, however, only focus on future forecasting, and we cannot interpret the time-series model.

Focusing on the state space model, we can avoid a lack of interpretation. The Kalman filter (KF) is a powerful tool for estimate states of the linear Gaussian state space model in a parameter-given scenario [12]. Combined with an expectation maximization (EM) algorithm, we can estimate the parameters and the states, i.e., achieve noise reduction as refined states and short-term forecasting from the constructed model [9]. However, the algorithm requires a great deal of calculation time and is unsuitable for real-time situation.

We therefore propose a method by which to estimate the states and the state transition of the model in real time. Using three ideas, namely, assumption of an observation transition, a time-invariant interval, and an online learning framework, the method works well for some experiments. We refer to this method as linear operator construction with the Kalman filter (LOCK).

In addition, by introducing localization and spatial uniformity to the proposed method, referred to respectively as local LOCK (LLOCK) and spatially uniform LOCK (SLOCK), we can apply the proposed method to spatio-temporal data, such as image sequence data and grid sequence data. Applied to synthetic data, the proposed methods were found to be superior to the existing method and to be efficient in terms of calculation and memory cost.

We introduce the proposed methods in the following section. In Section 3, we state the experimental results obtained using a damped oscillation model, object moving, global flow, and local stationary flow data¹¹1Our code is available at https://github.com/ZoneTsuyoshi/lock. Finally, we conclude the present paper in Section 4.

2 Methods

In this section, we first explain the notation, followed by state Kalman filter as a benchmark method. We then propose new methods in order to estimate the states of the system and the transition matrices in real time.

2.1 Notation

We use bold small letters, e.g., $x$ , to denote the variable as vector. The set of natural numbers is denoted as $N$ , and we denote as $N_{n}$ as the set of natural numbers that are less than or equal to $n$ , i.e., $N_{n} = {1, 2, \dots, n}$ . For a matrix $M \in R^{d \times d}$ and a vector $v \in R^{d}$ that depends on time $t$ , we use $M_{t}$ and $v_{t}$ for the time dependency. For a matrix $M \in R^{d \times d}$ and a vector $v \in R^{d}$ , we use $M_{i j}$ and $v_{i}$ to denote the $(i, j)$ -th element of the matrix and the $i$ -th element of the vector, respectively. Moreover, for a matrix $M \in R^{d \times d}$ and a vector $v \in N^{n} (n \leq d)$ , $M_{v, v}$ denotes the matrix

⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} M_{v_{1} v_{1}} & \dots & M_{v_{1} v_{n}} ⋮ & ⋱ & ⋮ M_{v_{n} v_{1}} & \dots & M_{v_{n} v_{n}} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠ \in R^{n \times n} .

Similarly, for vectors $v \in R^{d}$ and $u \in N^{n} (n \leq d)$ , $v_{u}$ denotes the vector $(v_{u_{1}}, \dots, v_{u_{n}})^{T} \in R^{n}$ . Moreover, for a set $A \subset R$ that satisfies $| A | = n$ , $v e c (A)$ denotes the vector, the $i$ -th element of which corresponds to $i$ -th smallest value of the matrix.

2.2 Kalman Filter

Consider the state space model

	$x_{t}$	$= f_{t} (x_{t - 1}) + v_{t},$		(1)
	$y_{t}$	$= h_{t} (x_{t}) + w_{t},$		(2)

where $x_{t} \in R^{m}$ , $y_{t} \in R^{l}$ , $v_{t} \in R^{m}$ , and $w_{t} \in R^{l}$ represent the state of the system, the observation, the system noise, and the observation noise, respectively, at a given time $t$ . In the linear scenario, using

	$f_{t} (x_{t - 1})$	$= F_{t} x_{t - 1},$		(3)
	$h_{t} (x_{t})$	$= H_{t} x_{t},$		(4)

we can obtain

	$x_{t}$	$= F_{t} x_{t - 1} + v_{t},$		(5)
	$y_{t}$	$= H_{t} x_{t} + w_{t},$		(6)

from equations (1) and (2

). Suppose the system noise and the observation noise follow a Gaussian distribution, the mean of which is 0, i.e.,

v_{t} \sim N (0, Q_{t})

and

w_{t} \sim N (0, R_{t})

. Under this assumption, we can apply the Kalman filter (KF) to the estimation of the state vectors.

The filter uses the following update formula

$x_{t \| t - 1}$	$= F_{t} x_{t - 1 \| t - 1},$	(7)
$V_{t \| t - 1}$	$= F_{t} V_{t - 1 \| t - 1} F_{t}^{T} + Q_{t},$	(8)
$K_{t}$	$= V_{t \| t - 1} H_{t}^{T} (H_{t} V_{t \| t - 1} H_{t}^{T} + R_{t})^{- 1},$	(9)
$x_{t \| t}$	$= x_{t \| t - 1} + K_{t} (y_{t} - H_{t} x_{t \| t - 1}),$	(10)
$V_{t \| t}$	$= V_{t \| t - 1} - K_{t} H_{t} V_{t \| t - 1},$	(11)

where $x_{t | s}$ and $V_{t | s}$ represent the estimated value and the covariance, respectively, of the state at time $t$ given observation ${y_{1}, \dots, y_{s}}$ .

The expectation maximization (EM) algorithm is one of the most powerful tools for estimation of the transition matrix. However, this algorithm is unsuitable for a real-time situation due to the requirement for a sequence of smoothed values of ${x_{t | T}}$ and ${V_{t | T}}$ , where $T$ is final time-step of time-series data. We can consider an ad hoc strategy, namely, expectation maximization KF (EMKF), which applies the algorithm every $τ$ step. However, in the next section, we show that the proposed method only works for a limited situation. Thus, the following sections propose new methods for estimation of the transition matrices and the states in real time.

2.3 Linear Operator Construction with Kalman Filter (LOCK)

In order to overcome the real-time estimation problem, we propose a new method referred to as linear operator construction with the Kalman filter (LOCK) by introducing three ideas: assumption of an observation transition, a time-invariant interval, and an online learning framework.

First, we introduce the assumption of a linear and Gaussian observation transition

y_{t} = G_{t} y_{t - 1} + u_{t},

(12)

where $u_{t} \in R^{l}$ is the noise of the observation transition at time $t$ and follows Gaussian distribution $N (0, S_{t})$ . If the observation matrices are regular, this assumption is satisfied by

	$G_{t}$	$= H_{t} F_{t} H_{t - 1}^{- 1},$		(13)
	$S_{t}$	$= H_{t - 1}^{- T} F_{t}^{T} H_{t}^{T} R_{t - 1} H_{t} F_{t} H_{t - 1}^{- 1} + H_{t}^{T} Q_{t} H_{t} + R_{t} .$		(14)

Figure 1 shows these relationships.

Using equations (5), (6), and (12), we can obtain

y_{t} = H_{t} (F_{t} x_{t - 1} + v_{t}) + w_{t} = G_{t} (H_{t - 1} x_{t - 1} + w_{t - 1}) + u_{t} .

(15)

Then, by taking the expectation, the following relation is satisfied:

(H_{t} F_{t} - G_{t} H_{t - 1}) E [x_{t - 1}] = 0.

(16)

Therefore, we can obtain an unbiased estimator

{^F}_{t} = H_{t}^{-} G_{t} H_{t - 1},

(17)

where $H_{t}^{-}$ and ${^F}_{t}$ indicate the Moore-Penrose type pseudo-inverse matrix of $H_{t}$ and the estimated value of $F_{t}$ , respectively.

In addition, suppose these matrices are time-invariant in interval $τ$ , then, we can obtain

^G = Y_{t} Y_{t - 1}^{-},

(18)

where $Y_{t} = (y_{t - τ + 1}, \dots y_{t})$ is a matrix composed of the observation vectors from time $t - τ + 1$ to time $t$ . Therefore, we can obtain the estimation algorithm

^F = H^{-} Y_{t} Y_{t - 1}^{-} H .

(19)

Moreover, based on an idea incorporated from online learning[4]

, we introduce parameters in order to treat outliers in the observation

	$F_{t} = F_{t - τ} - η c r o p (F_{t - τ} -^F, - c, c),$		(20)
	${(c r o p (F, a, b))}_{i j} = min {max {F_{i j}, a}, b},$		(21)

where $η$ and $c$ are the learning rate and the cutoff distance, respectively, which control the maximum amount of difference between the old estimate and the new estimate. Algorithm 1 summarizes the LOCK method.

Figure 1: Schematic representation of the relationships among (5), (6), and (12). The state at time $t - 1$ decides the observation at time $t - 1$ and the state at time $t$ by $H_{t - 1}$ and $F_{t}$ , respectively. Moreover, the observation at time $t - 1$ decides the observation at time $t$ by $G_{t}$ .

x_{0 | 0}, V_{0 | 0}, η, c, τ

Set initial time

t = 0

while observation at a time

t

y_{t}

exists do

Calculate

x_{t | t - 1}, V_{t | t - 1}

by equations (7) and (8)

Calculate

K_{t}, x_{t | t}, V_{t | t}

by equations (9),(10), and (11)

t \equiv 0 (mod τ)

then

Update

F

by the equation (20)

end if

t \leftarrow t + 1

end while

Algorithm 1 Algorithm of LOCK

2.4 Local LOCK (LLOCK)

We introduced the LOCK method in the previous section. However, this method cannot function properly for the case in which $τ < l$ , where $τ$ and $l$ are the update interval and the dimension of the observation, respectively. In particular, applying LOCK to movie data of large dimension, due to the requirement for a large $τ$ , is unrealistic for a real-time situation. Hence, we propose the local LOCK (LLOCK) method, which adopts localization in order to overcome this problem[11].

Consider a lattice point in spatio-temporal data. The variable of the point at time $t$ is affected only by the neighborhood of the point at time $t - 1$ . This localization can reduce the effective observation dimension. In other words, in order to update the $(i, j)$ -th element of the transition matrix, we can use only variables of the observation in the vicinity of points $i$ and $j$ .

We define localization matrix $L = A + I$ , where $A$ is adjacency matrix

A_{i j} = {\begin{matrix} 1, & (i th variable is close to j % th variable), 0, & (e l s e) . \end{matrix}

(22)

Using this matrix, we can obtain

	$^Gij={(Yt)k,k((Yt−1)k,k)−}lilj,$		(23)
	$k = v e c ({k \in N_{l} ∣ L_{i k} + L_{j k} > 0}),$		(24)
	$l_{i} = \| {k \in N_{i} ∣ L_{i, k} + L_{j, k} > 0} \|,$		(25)

where $l_{i}$ and $l_{j}$ are local indices of points $i$ and $j$ , respectively. Algorithm 2 summarizes the local calculation of the observation transition matrix $G$ .

L, {y_{s}}_{s = t - τ}^{t}

for

i \in {1, \dots, l}

for

j \in {1, \dots, l}

L_{i, j} = 1

then

Construct the vector

k

Calculate the local indices

l_{i}

and

l_{j}

corresponding to the global indices

i

and

j

Localize

{y_{s}}_{s = t - τ}^{t}

{y_{s}^{l o c a l}}_{s = t - τ}^{t}

, where

y_{s}^{l o c a l} = (y_{s})_{k}

Y_{t - 1}^{l o c a l} \leftarrow (y_{t - τ}^{l o c a l}, \dots, y_{t - 1}^{l o c a l})

Y_{t}^{l o c a l} \leftarrow (y_{t - τ + 1}^{l o c a l}, \dots, y_{t}^{l o c a l})

{^G}^{l o c a l} \leftarrow Y_{t}^{l o c a l} (Y_{t - 1}^{l o c a l})^{-}

{^G}_{i j} \leftarrow {^G}_{l_{i} l_{j}}^{l o c a l}

end if

end for

Algorithm 2 Local calculation of the observation transition matrix

2.5 Spatially Uniform LOCK (SLOCK)

The LOCK and LLOCK methods consider the elements of the transition matrix as independent parameters, as a result, these methods treat many parameters and require a long update interval in order to avoid $τ < l$ . Thus, both methods are unsuitable for rapid change of the transition matrix. Therefore, we introduce the spatially uniform LOCK (SLOCK) method to solve this problem. This method considers multiple elements of the transition matrix as a common parameter. In other words, we assume that parameters are spatially uniform in given vector field. Since the number of parameters is smaller than that of LOCK or LLOCK, the algorithm of this method can update with a shorter interval.

First, we design a matrix $P \in ({0} \cup N)^{l} \times ({0} \cup N)^{l}$ , each element of which represents a parameter number, in which if two elements are same, then the corresponding two elements of $G$ are same, i.e.,

P_{i_{1} j_{1}} = P_{i_{2} j_{2}} \Rightarrow G_{i_{1} j_{1}} = G_{i_{2} j_{2}} .

(26)

In addition, similar to the localization matrix, the zero elements of $P$ correspond to the zero elements of $G$ , i.e.,

P_{i j} = 0 \Rightarrow G_{i j} = 0.

(27)

If we use the numbers in ascending order, then we have $α \leq l^{2}$ unique values, excluding zero, which are ${1, \dots, α}$ . Next, we use the vector $θ \in R^{α}$ to denote the values of $G$

G_{i j} = {\begin{matrix} θ_{P_{i j}}, & (P_{i j} > 0), 0, & (P_{i j} = 0) . \end{matrix}

(28)

Thanks to this notation, we can obtain

	$y_{t} = Ξ_{t - 1} θ,$		(29)
	$(Ξ_{t})_{i j} = \sum k \in B_{i j} (y_{t})_{k}, Ξ \in R^{l} \times R^{α}$		(30)
	$B_{i, j} = {k \in N_{l} ∣ P_{i k} = j} (i \in N_{l}, j \in N_{α}),$		(31)

where $θ$ is the parameter vector. Using this equation, we can obtain the following update equations:

$^θ$	$= Ξ_{t - 1}^{-} y_{t},$	(32)
${(^G)}_{i j}$	$= {(^θ)}_{P_{i j}},$	(33)
$^F$	$= H_{t}^{-}^G H_{t - 1} .$	(34)

Algorithm 3 summarizes this method.

P, {y_{s}}_{s = t - 1}^{t}

for

i \in {1, \dots, l}

for

j \in {1, \dots, l}

Construct the set

B_{i j}

Calculate the

(i, j)

-th element of the matrix

Ξ_{t - 1}

end for

Calculate

^θ,^G,^F

by equations (32)-(34)

Algorithm 3 Estimation algorithm of the transition matrix used by SLOCK

3 Results

We applied EMKF and LOCK to the damped oscillation model, SLOCK to object moving data and global flow data, and LLOCK to global flow data and local stationary flow data.

3.1 Damped Oscillation Model

3.1.1 Model Equation

Damped oscillation model

	$\frac{d x}{d t}$	$= v,$		(35)
	$m \frac{d v}{d t}$	$= - k x - r v$		(36)

represents the behavior of a damped oscillator, which is subject to a damping force $r v$ , where $x$ is the position of the object, $v$ is the velocity, $m$ is the mass, $k$ is the oscillation constant, and $r$ is the damping constant. Using the Euler forward method, we can obtain

	$x_{t}$	$= x_{t - 1} + Δ t v_{t - 1},$		(37)
	$v_{t}$	$= - \frac{k}{m} Δ t x_{t - 1} + (1 - \frac{r}{m} Δ t) v_{t - 1},$		(38)

where $Δ t$ is the step size. By combining $x_{t}$ and $v_{t}$ into the state vector $x_{t}$ , we can obtain

	$x_{t} = F^{t r u e} x_{t - 1},$		(39)
	$x_{t} = (\begin{matrix} x_{t} v_{t} \end{matrix}), F^{t r u e} = (\begin{matrix} 1 & Δ t - k Δ t / m & 1 - r Δ t / m \end{matrix}),$		(40)

where $F^{t r u e}$ is the true transition matrix of this model. Appling this model to the linear state space model, we can use the EMKF and LOCK methods.

3.1.2 Experimental Conditions

We conducted five experiments. In the first and second experiments, the parameters of the model are time-invariant. In the other experiments, the parameters of the model are time-variant. The five experiments use the same setting, excluding the model parameters, as shown in Table 1.

In experiment 1, we first simulated the true model from $x_{0} = (5.0, 0.0)^{T}$ for $T = 100$ adding system noise $v_{t} \sim N (0, {0.01}^{2} \times I)$ (referred to herein as “true” states). Then, adding observation noise $w_{t} \sim N (0, {0.2}^{2} \times I)$ , we can obtain “pseudo-observation” data. Then, we applied EMKF and LOCK to the observation and obtained the experimental results. Here, the initial transition matrix $F_{0}$ is the true one, and the initial state and the transition covariance are different from the condition for generating the true data.

Experiment 2 differs from the previous experiment with respect to the initial transition matrix and the number of simulations. The initial transition matrix is not the true one but rather follows isotropic Gaussian distribution $I N (F^{t r u e}, {1.0}^{2} \times I)$ , and each element $(i, j)$ of the matrix independently follows Gaussian distribution $N (F_{i j}, {1.0}^{2})$ . Regarding the number of simulations, considering randomness, we executed 100 simulations and obtained the experimental results.

The main difference between experiment 3 and experiments 1 and 2 is whether parameters are time-invariant or time-variant. Experiments 3 and 4 use the same settings as experiments 1 and 2, respectively, excluding $k$ , $r$ , and $η$ (see Table 1). Experiments 4 and 5 have different distributions of the initial transition matrices, which are $I N (F_{0}^{t r u e}, {1.0}^{2} \times I)$ and $I N (F_{0}^{t r u e}, {0.01}^{2} \times I)$ for experiments 4 and 5, respectively.

damped oscillation model	$Δ t = 1.0$ , $m = 1.0$ , final time-step $T = 100$
state space model	$x_{0} = (\begin{matrix} 60 \end{matrix})$ , $H = V_{0} = I$ , $R = {0.2}^{2} \times I$ , $Q = {0.01}^{2} (\begin{matrix} 0 & 0 0 & (Δ t / m)^{2} \end{matrix})$
LOCK	$τ = 4$ , $c = 0.5$
EMKF	$τ = 4$ , $c = 0.5$ , #(iteration) $= 5$
experiment 1	$k = 0.5$ , $r = 0.52$ , $η = 0.6$ , $F_{0} = F^{t r u e}$
experiment 2	$k = 0.5$ , $r = 0.52$ , $η = 0.6$ , $F_{0} \sim I N (F^{t r u e}, {1.0}^{2} \times I)$
experiment 3	$k (t) = 0.65 t / T + 0.35 (1 - t / T)$ , $r (t) = 0.37 t / T + 0.67 (1 - t / T)$ ,
experiment 3	$η = 0.8$ , $F_{0} = F_{0}^{t r u e}$
experiment 4	$k (t) = 0.65 t / T + 0.35 (1 - t / T)$ , $r (t) = 0.37 t / T + 0.67 (1 - t / T)$ ,
experiment 4	$η = 0.8$ , $F_{0} \sim I N (F_{0}^{t r u e}, {1.0}^{2} \times I)$
experiment 5	$k (t) = 0.65 t / T + 0.35 (1 - t / T)$ , $r (t) = 0.37 t / T + 0.67 (1 - t / T)$ ,
experiment 5	$η = 0.8$ , $F_{0} \sim I N (F_{0}^{t r u e}, {0.01}^{2} \times I)$

Table 1: Parameters of the experiments using the damped oscillation model

3.1.3 Experimental Results and Considerations

Figure 2 illustrates the time transition of the true states, the observation, and the estimated results used by LOCK and EMKF in experiment 1. From this figure, LOCK and EMKF function well in terms of the estimation of the states.

Figure 3 shows similar results to those for experiment 3. Although the estimated results for LOCK sometimes extend away from the true states, both methods work well overall for state estimation of the time-variant model.

In addition, Figures 4 and 5 show the time transition of the mean squared error (MSE) between the true matrices and the estimated matrices used by LOCK and EMKF, respectively. The four panels correspond to the four elements of the matrices; e.g., the upper-left panel corresponds to $F_{11}$ . From Figure 4, the estimated matrices of LOCK are similar to the true matrix, even if the initial transition matrices are far from the true matrix. On the other hand, as shown in Figure 5, EMKF works only when the initial transition matrix is similar to the true matrix. In the real world, since the initial transition matrix is often unknown, LOCK is expected to be better than EMKF, from the perspective of approaching the true matrix.

Similar results for experiments 2 and 5 are presented in Appendix A.

Figure 2: Time transition of the true states, the observations, and the estimated results used by LOCK and EMKF in experiment 1.

Figure 3: Time transition of the true states, the observations, and the estimated results used by LOCK and EMKF in experiment 3.

Figure 4: Time transition of the true transition matrices (red) and the estimated transition matrices by LOCK (blue) in experiment 4. The four panels correspond to the four elements of the matrices; e.g., the upper-left panel corresponds to $F_{11}$ .

Figure 5: Time transition of the true transition matrices (red) and the estimated transition matrices by EMKF (blue) in experiment 4. The four panels correspond to the four elements of the matrices; e.g., the upper-left panel corresponds to $F_{11}$ .

3.2 Object Moving

3.2.1 Experimental condition

We generated movie data in which an object moves various directions for $T = 100$ , referred to herein as “object moving” data. Details are presented in Appendix B.1. After generating the true data, we added Gaussian noise $N (0, 20^{2})$ . Figure 6 shows the pseudo-observation at $t = 0, 5, 10, 15, 20$ . The directions of transition often change and summarized in Figure 7. In addition, we used the “adjacency distance” as the neighbors considered during the localization phase, as shown in Figure 8. We set the parameters of SLOCK as $τ = 1$ , $η = c = 1.0$ , and $d = 1$ and the state space model given by $F_{0} = H = V_{0} = I$ and $Q = R = {0.2}^{2} \times I$ .

Figure 6: Observation of object moving data at time $t = 0, 5, 10, 15, 20$ .

Figure 7: Directions of transition of the object at each interval.

Figure 8: Relationship between the adjacency distance $d$ and the localized range. For example, $d = 1$ represents that a value at the $i$ -th grid at time $t$ is only affected by neighbor $3 \times 3$ values at time $t - 1$ .

3.2.2 Experimental Results and Considerations

Figure 9 shows the time transition of the MSE of the observations and estimated results used by KF and SLOCK. From this figure, since the MSE of SLOCK is lower than that of KF and the observations, this method is a powerful tool for estimating the spatially uniform model.

In addition, we calculated the MSE for the transition matrices and executed sensitivity analysis for $d$ , $η$ , and $c$ (as shown in Appendix B.2).

Figure 9: Time transition of the MSE of the observations and estimated results used by KF and SLOCK.

3.3 Global Flow

3.3.1 Experimental Conditions

We generated “global flow” data, where various objects exist in images and move in each direction in each interval, as shown in Figure 10. A more detailed description of the generating process is stated in Appendix C.1. We added Gaussian noise $N (0, 20^{2})$ at each grid and obtained pseudo-observation data. We set the parameters of LLOCK to $τ = 50$ , $η = 0.8$ , $c = 1.0$ , and $d = 1$ and the state space model given by $F_{0} = H = V_{0} = I$ and $Q = R = {0.2}^{2} \times I$ .

Figure 10: (Left) Initial image of the global flow data. (Other panels) Flows of objects at each interval.

3.3.2 Experimental Results and Considerations

Figure 11 shows the time transition of the MSE of the observations and estimated results used by KF and LLOCK. As shown in this figure, LLOCK provides better estimated results, excluding near the changing points. In this data setting, since the changes are rapid, the MSE of the proposed method is worse than the observation in the vicinity of the changes.

In addition, we conducted short-term prediction by KF and LLOCK. The results are shown in Figure 12. We simulate each method until $t = 200$ and perform prediction thereafter, where the predicted MSE of the observations represent the MSE between the true state at a time and the observation at $t = 200$ . From this, the predictive ability of LLOCK is shown to be superior to that of KF and ad hoc observation for this data.

Figures regarding the MSE for the transition matrices and sensitivity analysis for $τ$ , $η$ , and $c$ are shown in Appendix C.2.

Moreover, we measured the calculation time of updating the transition matrices for the case in which images are $10 \times 10$ , $20 \times 20$ , and $30 \times 30$ . We applied SLOCK and LLOCK to data in this setting, as shown in Figure 13. From this figure, the calculation times of both methods are low enough to be executed in a real-time situation.

In addition, we calculated ideal memory cost and ad hoc memory cost for the proposed methods and the EM algorithm. Figure 14 shows these results when we assume a float-type array. Here, “ad hoc” means that we execute no memory-saving code, i.e., that we memorize the $l \times l$ matrices as is. However, thanks to localization and the spatially uniform assumption, we need only manage $l \times N_{l}$ arrays, referred to as the “ideal” cost, where $N_{l}$ represents the local dimension. According to this figure, the ideal memory cost is much less than that of the EM algorithm and the methods are easier to apply.

Figure 11: Time transition of the MSE of the observations and estimated results used by KF and LLOCK.

Figure 12: Time transition of the MSE of the observations and predicted results used by KF and LLOCK when we obtain data until $t = 200$ and perform prediction thereafter.

Figure 13: Calculation time when the observation dimensions are varied for (Left) SLOCK and (Right) LLOCK.

Figure 14: Memory usage when observation dimensions are varied assuming a float-type array.

3.4 Local Stationary Flow

3.4.1 Experimental Condition

We generated “local stationary flow” data for which objects spring up in the boundary and move in four directions corresponding to each field, as shown in Figure 15. The left-hand panel represents the flow of the data. For example, the flow direction of the upper-left part of the image is upward. The other panels show the pseudo-observation data at $t = 0$ , 5, and 10. A more detailed explanation of the generating process is shown in Appendix D.1. Generating the true data, we added Gaussian noise $N (0, 20^{2})$ and gained the observations.

We then applied LLOCK to these data to confirm that the method can capture local information. We set the parameters of LLOCK as $τ = 50$ , $η = 0.6$ , $c = 1.0$ , and $d = 1$ and the state space model given by $F_{0} = H = V_{0} = I$ and $Q = R = {0.2}^{2} \times I$ .

Figure 15: Generation of local stationary flow data. (Left) Flow of the data. For example, the flow direction of the upper-left part of the image is upward. (Other) Pseudo-observation data at $t = 0$ , $5$ , and $10$ .

3.4.2 Experimental Results and Considerations

Figure 16 shows the time transition of the MSE of the observations and the estimated results used by KF and LLOCK. From this figure, LLOCK has a lower MSE than that of KF and the observations.

In addition, we performed prediction for the short period used by KF and LLOCK as shown in Figure 17. Similar to the results for global flow data, the proposed method has better performance than the other methods for short-term prediction.

The MSE regarding the transition matrix and sensitivity analysis for $τ$ , $η$ , and $c$ are presented in Appendix D.2.

Figure 16: Time transition of the MSE of the observations and estimated results used by KF and LLOCK.

Figure 17: Time transition of the MSE between the true states and the observations and the predicted results used by KF and LLOCK when we obtain data until $t = 500$ and perform prediction thereafter.

4 Conclusion

In the present paper, we propose three real-time methods to estimate the states and state transition matrices in a linear Gaussian state space model. The first proposed method, namely, linear operator construction with the Kalman filter (LOCK), can approach the true transition matrices and the true states via application of the damped oscillation model. The advance methods, namely, SLOCK and LLOCK, achieve better performance in terms of noise reduction and short-term prediction through the three synthetic data: object moving, global flow, and local stationary flow. These methods are also superior to the EM algorithm in terms of computational and memory cost. In fact, the calculation time and memory usage of these methods are much less than those of the existing method. Therefore, these methods have the potential to estimate the transition of data, such as weather forecasting and object tracking, in real time.

Nevertheless, the proposed methods have three main drawbacks: dependence on a linear Gaussian framework, tuning of hyper-parameters, and tight assumption of SLOCK. First, as the proposed methods use a linear Gaussian formulation, we cannot directly apply these methods to nonlinear or non-Gaussian data. Then, the proposed methods need to provide interpretability for such data. Second, the proposed methods include the hyper-parameters $τ$ , $η$ , and $c$ . An automatic tuning method for these parameters is necessary for application to real data. Finally, although SLOCK has better performance for spatially uniform data, this assumption is rather limiting for real data. In order to overcome this issue, we are eager to develop a combined method of LLOCK and SLOCK because these assumptions and precisions are in trade-off relationships.

References

[1] Ballas, N., Yao, L., Pal, C., and Courville, A. (2015), “Delving Deeper into Convolutional Networks for Learning Video Representations,” arXiv preprint arXiv:1511.06432.
[2] Bertalmío, M., (2019), Denoising of Photographic Images and Video -Fundamentals, Open Challenges and New Trends-, Springer.
[3] Batson, J. and Royer, L. (2019), “Noise2self: Blind Denoising by Self-supervision,” arXiv preprint arXiv:1901.11365.
[4] Bishop, C. M. (2006), Pattern Recognition and Machine Learning, Springer.
[5] Chen, S. Y. (2012), “Kalman Filter for Robot Vision: A Survey,” IEEE Transactions on Industrial Electronics, 59, 11, 4409-4420, DOI: 10.1109/TIE.2011.2162714.
[6] Comaniciu, D., Ramesh, V., and Meer, P. (2003), “Kernel-based object tracking,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 25, 5, 564-577, DOI: 10.1109/TPAMI.2003.1195991.
[7] Elad, M. (2010), Sparse and Redundant Representations -From Theory to Applications in Signal and Image Processing-, Springer.
[8] Evensen, G. (2003), “The ensemble Kalman filter: Theoretical formulation and practical implementation,” Ocean dynamics, 53, 343-367, DOI: 10.1007/s10236-003-0036-9.
[9] Ghahramani, Z. and Hinton, G.E. (1996), “Parameter Estimation for Linear Dynamics Systems,” Technical Report CRG-TR-96-2, Department of Computer Science, University of Toronto.
[10] Gordon, N. J., D. J. Salmond, and A. F. M. Smith (1993), “A novel approach to non-linear and non-Gaussian Bayesian state estimation,” IEE Proceedings F - Radar and Signal Processing, 140, 107-113, DOI: 10.1049/ip-f-2.1993.0015.
[11] Hunt, B. R., Koestelich, E.J., and Szunyogh, I (2007), “Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter,” Phydica D: Nonlinear Phenomena, 230, 112-126, DOI: 10.1016/j.physd.2006.11.008.
[12] Kalman, R.E. (1960), “A New Approach to Linear Filtering and Prediction Problem”, Transactions of the ASME-Journal of Basic Engineering, 82(D), 35-45, DOI: 10.1115/1.3662552.
[13] Kitagawa, G. (1998), “A self-organizing state-space model,” Journal of the American Statistical Association, 93, 443, 1203-1215, DOI: 10.2307/2669862.
[14] Krull, A., Buchholz, T., and Jug, F. (2019), “Noise2Void - Learning Denoising from Single Noisy Images,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2129-2137, DOI:10.1109/CVPR.2019.00223.
[15] Lefkimmiatis, S. (2018) “Universal denoising networks: A Novel CNN Architecture for Image Denoising,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 3204-3213, DOI: 10.1109/CVPR.2018.00338.
[16] Lehtinen, J., Munkberg, J., Hasselgren, J., Laine, S., Karras, T., Aittala, M., and Aila, T. (2018), “Noise2Noise: Learning Image Restoration without Clean Data,” arXiv preprint arXiv:1803.04189.
[17] Lotter, W., Kreiman, G., and Cox, D. (2016), “Deep Predictive Coding Networks for Video Prediction and Unsupervised Learning,” arXiv preprint arXiv:1605.08104.
[18] Shi, X., Chen, Z., Wang, H., and Yeung, D. (2015), “Convolutional LSTM Network: A Machine Learning Approach for Precipitation Nowcasting,” Advances in neural information processing systems, 802-810.
[19] Shi, X., Gao, Z., Lausen, L., Wang, H., Yeung, D., Wong, W., and Woo, W. (2017), “Deep Learning for Precipitation Nowcasting: A Benchmark and A New Model,” Advances in neural information processing systems, 5617-5627.
[20] Shumway, R. H., Stoffer, D. S. (2017), Time Series Analysis and Its Applications With R Examples. Springer.
[21] Tian, C., Xu, Y., Fei, L., and Yan, K. (2018), “Deep Learning for Image Denoising: A Survey,”
International Conference on Genetic and Evolutionary Computing
, 563-572, DOI: 10.1007/978-981-13-5841-8_59.
[22] Yang, H., Shao, L., Zheng, F., Wang, L., and Song, Z. (2011), “Recent advances and trends in visual tracking: A review,” Neurocomputing, 74, 18, 3823-3831, DOI: 10.1016/j.neucom.2011.07.024.
[23] Wang, Y., Zhang, J., Zhu, H., Long, M., Wang, J., and Yu, P., H.. (2019), “Memory In Memory: A Predictive Neural Network for Learning Higher-Order Non-Stationarity from Spatiotemporal Dynamics,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, DOI:10.1109/CVPR.2019.00937.
[24] Zhang, K., Zuo, W., Chen, Y., Meng, D., and Zhang, L. (2017), “ Beyond A Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising,” IEEE Transactions on Image Processing, 26, 7, 3142-3155, DOI: 10.1109/TIP.2017.2662206.

Appendix A Damped Oscillation Model

Here, we present the residual results of the damped oscillation model. Figures 18 and 19 show the estimated transition matrices for experiment 2 by LOCK and EMKF, respectively. From these figures, the estimated elements of LOCK approach the true elements, whereas the correspondence of the EMKF method does not.

In addition, Figures 20 and 21 correspond to experiment 5. Figure 20 indicates that the estimated matrices by LOCK can match the true matrix, excluding the last few updates, the noise ratio of which are higher than in other matrices. Figure 21 indicates that the 100 simulated results are similar in terms of the time transition of the matrices.

Figure 18: Time transition of true transition matrices (red) and estimated transition matrices by LOCK (blue) in experiment 2. The four figures correspond to four elements of the matrices. For instance, the upper-left panel corresponds to $F_{11}$ .

Figure 19: Time transition of the true transition matrices (red) and the estimated transition matrices by EMKF (blue) in experiment 2. The four panels correspond to the four elements of the matrices. For instance, the upper-left panel corresponds to $F_{11}$ .

Figure 20: Time transition of the true transition matrices (red) and the estimated transition matrices by LOCK (blue) in experiment 5. The four panels correspond to the four elements of the matrices. For instance, the upper-left panel corresponds to $F_{11}$ .

Figure 21: Time transition of the true transition matrices (red) and the estimated transition matrices by EMKF (blue) in experiment 5. The four panels correspond to the four elements of the matrices. For instance, the upper-left panel corresponds to $F_{11}$ .

Appendix B Object Moving Data

b.1 How to Generate Data

First, we constructed a $25 \times 25$ array, the elements of which are 20. Then, we created random $N = 15$ core points, the $x$ and $y$ coordinates of which exist from 7 to 16. Second, we selected two random core points from the points and generated a link, which links the two points, the width of which is $w = 2$ . Third, we assigned the coordinates in the link linear $z$ values from $(V_{m a x} - V_{m i n}) \times m / M$

, added a Gaussian noise, mean and standard deviation of which are both 10, where

V_{m i n} = 100

V_{m a x} = 150

m

, and

M

are the minimum

z

value, the maximum

z

value, the iteration number, and the number of iterations, respectively. We iterated the second and third processes for

M = 10

times and obtained the base true image as shown in Figure 22.

Fourth, we set the direction list and direction change list as shown in Figure 7. In each direction, we constructed the translation matrix regarding the direction and then produced the image vector and the matrix. As a result, we obtained the true data. Finally, we added a Gaussian noise having a mean and standard deviation of 0 and 20, respectively, to the true data. We applied the absolute operator to the noise because the image data are greater than 0. Figure 22 shows images of the true data and the pseudo-observation data. This generating process is also shown on our GitHub page.

Figure 22: (Left) True image of the object moving data. (Other) Observation data at the initial time.

b.2 Residual Results

Figure 23 shows the Substantial MSE (SMSE) between the estimated transition matrices and the true transition matrices formulated by

SMSE(F,^F)=1nz∑(i,j)∈{(i,j)|Lij>0∨Pij>0}(Fij−^Fij)2

(41)

where $n_{z}$ represents the number of non-zero elements of the localization matrix, i.e.,

n_{z} = l \sum i, j = 1 L_{i, j} = l \sum i, j = 1 δ (P_{i, j} > 0) .

(42)

From this figure, the estimated matrices are close to the true matrix, especially the elements having true values of zero.

In addition, we conducted sensitivity analysis as shown in Figure 24. This figure indicates that the hyper-parameters of SLOCK are robust for this data.

Figure 23: Time transition of the SMSE between the true transition matrices and the estimated correspondence. Here, “true $= 0$ ” indicates the calculation of the SMSE for elements having true values that are zero.

Figure 24: Sensitivity analysis regarding adjacency distance $d$ , learning rate $η$ , and cutoff distance $c$ . During analysis of the other parameters, we set $d = 1$ , $η = 1.0$ , and $c = 1.0$ .

Appendix C Global Flow

c.1 How to Generate Data

First, we made a $30 \times 30$ matrix, the elements of which are 20. Then, we created random objects, the sizes of which are randomly chosen from two to four and the values of which are i.i.d. $N (150, 20^{2})$ . Second, we set the directions at each interval and move the objects corresponding to the direction at each time interval. Finally, we added a zero-mean Gaussian noise, the standard deviation of which is 20. This generating process is also shown on our GitHub page.

c.2 Residual Results

We calculated the SMSE between the true transition matrices and the estimated matrices as shown in Figure 25. The results indicate the estimated matrices are close to the true matrices, excluding rapid change points.

In addition, Figure 26 represents the results of sensitivity analysis for $τ$ , $η$ , and $c$ . From this figure, LLOCK has robustness for the synthetic data.

Figure 25: Time transition of the SMSE between the true transition matrices and the estimated correspondence. Here, “true $= 0$ ” indicates the calculation of the SMSE for elements having true values of zero.

Figure 26: Sensitivity analysis regarding update interval $τ$ , learning rate $η$ , and cutoff distance $c$ . During analysis of the other parameters, we set $τ = 50$ , $η = 0.8$ , and $c = 1.0$ .

Appendix D Local Stationary Flow

d.1 How to Generate Data

First, we constructed a $(T + L - 1) \times L$ matrix, the elements of which are 20, where $T = 1000$ and $L = 15$ indicate the number of time steps and the block length, respectively. Then, we created random objects, the sizes of which were randomly chosen from two to four and the values of which follow Gaussian distribution $N (150, 20^{2})$ , that were substituted into the object values for the random elements of the matrix. Second, we set the directions of four square blocks, the size of which are $15 \times 15$ , divided from a $30 \times 30$ field, as shown in Figure 15. Then, we set the initial value from the first 15 columns of the source matrix. Third, we moved the objects following the block flow and set new source values from the $t + L$ column of the matrix. Finally, we added zero-mean Gaussian noise having a standard deviation of 20. This generating process is also shown on our GitHub page.

d.2 Residual Results

Figure 27 shows time transition of the SMSE between the true transition matrices and estimated correspondence. This figure shows that the estimated results by LLOCK are close to the values of the true matrix.

In addition, we conducted sensitivity analysis for $τ$ , $η$ , and $c$ , as shown in Figure 28. This figure shows that LLOCK is a robust method for this local stationary flow data.

Figure 27: Time transition of the SMSE between the true transition matrices and the estimated correspondence. Here, “true $= 0$ ” indicates calculation of the SMSE for elements having true values of zero.

Figure 28: Sensitivity analysis regarding update interval $τ$ , learning rate $η$ , and cutoff distance $c$ . During the analysis for other parameters, we set $τ = 50$ , $η = 0.8$ , and $c = 1.0$ .

Appendix E Computer Environment

We used a GPU resource and Python. The following is the detailed environment of the present study.

Calculation machine
- CPU: Intel Xeon E5-2670 2.6GHz (8core) x 2
- Memory: 64 GB
- HDD: SAS 300 GB $\times$ 2 (RAID 1)
- OS: SuSE 12.0 Enterprise LINUX
- GPU: NVIDIA Tesla K20, NVIDIA Quadro K5000
- Other: CUDA8.0
Python
- Python 3.7.3
- Anaconda 4.6.14
- cupy-cuda80 5.4.0
- numpy 1.16.2

$x_{t \| t - 1}$	$= F_{t} x_{t - 1 \| t - 1},$	(7)
$V_{t \| t - 1}$	$= F_{t} V_{t - 1 \| t - 1} F_{t}^{T} + Q_{t},$	(8)
$K_{t}$	$= V_{t \| t - 1} H_{t}^{T} (H_{t} V_{t \| t - 1} H_{t}^{T} + R_{t})^{- 1},$	(9)
$x_{t \| t}$	$= x_{t \| t - 1} + K_{t} (y_{t} - H_{t} x_{t \| t - 1}),$	(10)
$V_{t \| t}$	$= V_{t \| t - 1} - K_{t} H_{t} V_{t \| t - 1},$	(11)

Real-time Linear Operator Construction and State Estimation with The Kalman Filter

Authors

Related Research

Real-time Linear Operator Construction and State Estimation with Kalman Filter

Layered Based Augmented Complex Kalman Filter for Fast Forecasting-Aided State Estimation of Distribution Networks

Real-time Ionospheric Imaging of S4 Scintillation from Limited Data with Parallel Kalman Filters and Smoothness

A Framework for Recognizing and Estimating Human Concentration Levels

Automated Offside Detection by Spatio-Temporal Analysis of Football Videos

Combining Ensemble Kalman Filter and Reservoir Computing to predict spatio-temporal chaotic systems from imperfect observations and models

A New Adaptive Noise Covariance Matrices Estimation and Filtering Method: Application to Multi-Object Tracking

Code Repositories

lock

This week in AI

1 Introduction

2 Methods

2.1 Notation

2.2 Kalman Filter

2.3 Linear Operator Construction with Kalman Filter (LOCK)

2.4 Local LOCK (LLOCK)

2.5 Spatially Uniform LOCK (SLOCK)

3 Results

3.1 Damped Oscillation Model

3.1.1 Model Equation

3.1.2 Experimental Conditions

3.1.3 Experimental Results and Considerations

3.2 Object Moving

3.2.1 Experimental condition

3.2.2 Experimental Results and Considerations

3.3 Global Flow

3.3.1 Experimental Conditions

3.3.2 Experimental Results and Considerations

3.4 Local Stationary Flow

3.4.1 Experimental Condition

3.4.2 Experimental Results and Considerations

4 Conclusion

References

Appendix A Damped Oscillation Model

Appendix B Object Moving Data

b.1 How to Generate Data

b.2 Residual Results

Appendix C Global Flow

c.1 How to Generate Data

c.2 Residual Results

Appendix D Local Stationary Flow

d.1 How to Generate Data

d.2 Residual Results

Appendix E Computer Environment