# Parameter Estimation for Linear Dynamical Systems
This notebook details an algorithm for estimating the parameters of a linear dynamical system using the Expectation Maximization (EM) algorithm.

Given a set of observations (measurements), this offline algorithm iterates between the E-step that fixes the current parameters and computes posterior probabilities over the hidden states and an M-step that maximizes the expected log-likelihood of the parameters using the posterior distribution computed in the E-step.

For linear Gaussian systems, the E-step is the Kalman smoother and the M-step is a linear regression problem.

## Model
The state space is defined by the following:
$$ \begin{align*} x_{t+1} &= A x_t + w_t \\
y_t &= C x_t + v_t\end{align*}$$
where $w_t\sim \mathcal{N}(0,Q)$ is the process noise and $v_t\sim \mathcal{N}(0,R)$ is the measurement noise. It is assumed that $x_t \in \mathbb{R}^k$ and $y_t \in \mathbb{R}^p$

It is assumed that only the output $y_t$ is measured and the state variables and noise are hidden.

To simplify the analysis, we combine the state variable with the process noise and measurement with its noise into single Gaussian random variables.

$$ \begin{align*} P(x_{t+1}|x_t) &= \exp{\{-\frac{1}{2} \left[x_{t+1}-A x_t\right]^T Q^{-1} \left[x_{t+1}-A x_t\right]\}} (2 \pi)^{-\frac{k}{2}}|Q|^{-\frac{1}{2}} \\
P(y_{t}|x_t) &= \exp{\{-\frac{1}{2} \left[y_{t}-C x_t\right]^T R^{-1} \left[y_{t}-C x_t\right]\}} (2 \pi)^{-\frac{p}{2}}|R|^{-\frac{1}{2}} \end{align*}$$



## Kalman Smoother Estimator
The Kalman Smoother Estimator is an algorithm that seeks to estimate the trajectory of the latent state variables given a sequence of measurements. The state evolution and the output equations are modeled by the state space equation described above.

The optimal estimator $x_t^n = \mathbb{E}(x_t|y_1,y_2,...,y_n)$ by minimizing the mean squared error $P_t^n = \mathbb{E}((x_t-x_t^n)^T (x_t-x_t^n) | y_1,y_2,...,y_n)$.

### Forward Recursion
The forward recursion follows the Kalman Filter equations where the state estimate is updated based on current measurement.
#### Prediction
1. Predicted State Estimate: $x_{t|t-1} = A x_{t-1|t-1}$
2. Predicted Covariance Estimate: $V_{t|t-1} = A V_{t-1|t-1} A^T + Q $

#### Update
1. Measurement Residual: $z_t = y_{t} - C x_{t|t-1}$
2. Innovation Covariance: $S_t = C V_{t|t-1} C^T + R$
3. Optimal Kalman Gain: $K_t = P_{t|t-1} C^T S_t^{-1}$
4. Updated State Estimate: $x_{t|t} = x_{t|t-1} + K_t z_t$
5. Covariance Update: $V_{t|t} = (I-K_t C) V_{t|t-1}$

### Backward Recursion
The backward recursion performs another update to the state estimates based on future measurements.

Starting from time $t=n$, with $x_{n|n}$ & $V_{n|n}$ obtained from the Kalman filter:
1. Define $J_{t-1} = V_{t-1|t-1} A V_{t|t-1}^{-1}$
2. $x_{t-1|n} = x_{t-1|t-1} + J_{t-1}(x_{t|n}-A x_{t-1|t-1})$
3. $V_{t-1|n} = V_{t-1|t-1} + J_{t-1}(V_{t|n}-V_{t|t-1})J_{t-1}^T$

---

A useful quantity that is not used in the Kalman smoother, but is used in the EM algorithm is $V_{t,t-1}$. This can be computed during the backward recursion for convenience.
Starting from time $t=n$, with $V_{n,n-1|n} = (I-K_n C)A V_{n-1|n-1}$:
$$V_{i-1,i-2|n} = V_{i-1|i-1}J_{i-2}^T + J_{i-1}(V_{i,i-1|n}-A V_{i-1|i-1})J_{i-2}^T$$


## EM Algorithm
Assuming Markov property, and defining ${x} \equiv {x_1,x_2,...,x_n}$ & ${y} \equiv {y_1,y_2,...,y_n}$
$$ \begin{align} P(\{x\},\{y\}) = &P(x_0) \prod_{t=0}^{n-1} P(x_{t+1}|x_t) \prod_{t=0}^n P(y_t|x_t) \\
\Rightarrow \log{P(\{x\},\{y\})} = &-\sum_{t=1}^n{\frac{1}{2} \left[x_{t+1}-A x_t\right]^T Q^{-1} \left[x_{t+1}-A x_t\right]} \\
 &-\sum_{t=0}^n{\frac{1}{2} \left[y_{t}-C x_t\right]^T R^{-1} \left[y_{t}-C x_t\right]}  \\
 &-\frac{1}{2} \left[x_{0}- \pi_0\right]^T V_0^{-1} \left[x_{0}- \pi_0\right]  \\
 &-\frac{n-1}{2}\log{|Q|} - \frac{n}{2}\log{|R|} - \frac{n(p+k)}{2}\log{2 \pi}\end{align}$$
 where $P(x_0) \sim \mathcal{N}(\pi_0,V_0)$ with $\pi_0$ & $V_0$ to be estimated.

### E Step
The E-step of the algorithm requires computing the expected log likelihood $\mathcal{Q} = \mathbb{E}\left[\log P(\{x\},\{y\})|\{y\} \right]$ (Note: Its $\mathbb{E}\left(\log P(\cdot)|\{y\}\right)$).

To compute the above Expectation, we define:
$$\begin{align} \hat{x}_t &= \mathbb{E}(x_t|{y}) = x_{t|n}\\
P_t &= \mathbb{E}(x_t x_t^T|{y}) = V_{t|n} + x_{t|n} x_{t|n}^T \\
P_{t,t-1} &= \mathbb{E}(x_t x_{t-1}^T|{y}) = V_{t,t-1|n} + x_{t|n} x_{t-1|n}^T \end{align} $$
 and can be computed as part of the Kalman smoother estimator described above.

### M Step
The parameters of this system are $A$ , $C$ , $Q$ , $R$ , $\pi_0$ , $V_0$. Each of this can be computed by setting the corresponding partial derivative of the expected log likelihood to zero.

#### Output Matrix:
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial C} &= \sum_{t=0}^n R^{-1} \left[y_t - C x_t \right] (-x_t)^T \\
&= \sum_{t=0}^n \Bigg( R^{-1} \left[ y_t x_t^T - C P_t\right]\Bigg) \\
&= 0 \iff C_{\text{new}} = \Bigg(\sum_{t=0}^n y_t x_t^T \Bigg) \Bigg(\sum_{t=0}^n P_t \Bigg)^{-1}\end{align}$$

#### Output Noise Covariance:
Here we make use of the following property: $\frac{\partial{\log|R|}} {\partial{R^{-1}}} = R$
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial R^{-1}} &= \frac{n}{2} R -\sum_{t=0}^n{\left[ \frac{y_t y_t^T}{2} - C x_t y_t^T + \frac{Cx_t x_t^T C^T}{2}\right]} \\ &= \frac{n}{2} R -\sum_{t=0}^n{\left[ \frac{y_t y_t^T}{2} - C x_t y_t^T + \frac{C P_t C^T}{2}\right]}\\ &= \frac{n}{2} R -\sum_{t=0}^n{\left[ \frac{y_t y_t^T}{2} - \frac{C x_t y_t^T}{2} - \frac{C}{2} (x_t y_t^T - P_t C^T)\right]} \\ 
&= \frac{n}{2} R -\sum_{t=0}^n{\left[ \frac{y_t y_t^T}{2} - \frac{C x_t y_t^T}{2} \right] }\\
&= 0 \iff R_{\text{new}} =  \frac{1}{n} \left[\sum_{t=0}^n{ \left(y_t y_t^T - C_{\text{new}} x_t y_t^T\right)} \right] \end{align}$$

#### State Transition Matrix:
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial A} &= \sum_{t=0}^{n-1} Q^{-1} \left[x_{t+1} - A x_t \right] (-x_t)^T \\
&= \sum_{t=0}^{n-1} \Bigg( Q^{-1} \left[ x_{t+1} x_t^T - A P_t\right]\Bigg) \\
&= 0 \iff A_{\text{new}} = \Bigg(\sum_{t=0}^{n-1} P_{t+1,t} \Bigg) \Bigg(\sum_{t=0}^{n-1} P_t \Bigg)^{-1}\end{align}$$

#### Process Noise Covariance:
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial Q^{-1}} &= \frac{n-1}{2}Q -\sum_{t=0}^{n-1}{\left[ \frac{x_{t+1} x_{t+1}^T}{2} - \frac{A x_t x_{t+1}^T}{2}  - \frac{x_{t+1}x_t^T A^T}{2} + \frac{Ax_t x_t^T A^T}{2}\right]} \\ &= \frac{n-1}{2}Q - \frac{1}{2}\sum_{t=0}^{n-1}{\left[ P_{t+1} - A P_{t,t+1}  - P_{t+1,t} A^T + A P_t A^T\right]} \\
&= \frac{n-1}{2}Q - \frac{1}{2}\sum_{t=0}^{n-1}{\left[ P_{t+1} - A P_{t,t+1}  - (P_{t+1,t} - A P_t) A^T\right]} \\
&= \frac{n-1}{2}Q - \frac{1}{2}\sum_{t=0}^{n-1}{\left[ P_{t+1} - A P_{t,t+1} \right]} \\
&= 0 \iff Q_{\text{new}} = \frac{1}{n-1} \left(\sum_{t=0}^{n-1}{P_{t+1} }- A_{\text{new}}\sum_{t=0}^{n-1}{P_{t,t+1}}\right)\end{align}$$

#### Initial State Mean
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial \pi_0} &= (x_0-\pi_0)V_0^{-1} \\
&= 0 \iff \pi_{0,\text{new}} = x_0 \end{align}$$

#### Initial State Covariance
$$ \begin{align}\frac{\partial \mathcal{Q}}{\partial V_{-1}} &= \frac{V_1}{2} - \frac{1}{2}(P_0-x_0\pi_0^T-\pi_0x_p^T+\pi_0\pi_0^T) \\
&= 0 \iff V_{0,\text{new}} = P_0 - x_0 x_0^T \end{align}$$

## References
1. Ghahramani, Z., & Hinton, G. E. (1996a). Parameter estimation for linear dynamical systems (Tech. Rep. No. CRG-TR-96-2). Toronto: Department of Computer Science, University ofToronto. Available [online](http://www.gatsby.ucl.ac.uk/»zoubin/papers/tr-96-2.ps.gz).
2. Shumway, R. H., & Stoffer, D. S. (1982). An approach to time series smoothing and forecasting using the EM algorithm. J. Time Series Analysis, 3(4), 253–264.

59453a1de3beb3a9bfc4306e5981531265222957
ipynb
Jupyter Notebook
files/DynamicalSystems/LinearSystemIdentification_EM.ipynb
chandrusuresh/MyNotes
4e0f86195d6d9eb3168bfb04ca42120e9df17f0b
MIT
files/DynamicalSystems/LinearSystemIdentification_EM.ipynb
chandrusuresh/MyNotes
4e0f86195d6d9eb3168bfb04ca42120e9df17f0b
MIT
files/DynamicalSystems/LinearSystemIdentification_EM.ipynb
chandrusuresh/MyNotes
4e0f86195d6d9eb3168bfb04ca42120e9df17f0b
MIT
Qwen/Qwen-72B
1. YES
2. YES