The purpose of this project is to forecast sales volume and sales revenue. In order to lay the ground work for an approach based on a bayesian reactive message passing approach, the Part 1 project will work through the following article to gain a thorough understanding of its contents:
Podusenko, A., Kouw, W. M., & de Vries, B. Message Passing-Based Inference for Time-Varying Autoregressive Models
Follow-up parts will make use of implementation code in Julia to address the client’s problem.
It is not usually possible to make effective analysis annotations and expansions inside a pdf document. This notebook attempts to provide a remedy. Its contents is in the form of:
My additions are usually included within {{my addition}} or [[my addition]]. If in doubt about any part of the contents of this notebook, please consult the original paper/article.
Time-varying autoregressive (TVAR) models are widely used for modeling of non-stationary signals. Unfortunately, online joint adaptation of both states and parameters in these models remains a challenge. In this paper, we represent the TVAR model by a factor graph and solve the inference problem by automated message passing-based inference for states and parameters. We derive structured variational update rules for a composite “AR node” with probabilistic observations that can be used as a plug-in module in hierarchical models, for example, to model the time-varying behavior of the hyper-parameters of a time-varying AR model. Our method includes tracking of variational free energy (FE) as a Bayesian measure of TVAR model performance. The proposed methods are verified on a synthetic data set and validated on real-world data from temperature modeling and speech enhancement tasks.
Keywords: Bayesian inference; free energy; factor graph; hybrid message passing; model selection; non-stationary systems; probabilistic graphical models
Autoregressive (AR) models are capable of describing a wide range of time series patterns [1,2]. The extension to
Time-Varying AR (TVAR) models, where the
AR coefficients are allowed to vary over time, supports
tracking of non-stationary signals.
TVAR models have been successfully applied to a wide range of applications, including speech signal processing [3–5], signature verification [6], cardiovascular response modeling [7], acoustic signature recognition of vehicles [8], radar signal processing [9], and EEG analysis [10,11].
The realization of TVAR models in practice often poses some computational issues. For many applications, such as speech processing in a hearing aid, both a low computational load and high modeling accuracy are essential.
The problem of parameter tracking in TVAR models has been extensively explored in a non-Bayesian setting. For example, ref. [12] employs over-determined modified Yule-Walker equations and [13] applies the covariance method to track the parameters in a TVAR model. In [14], expressions for the mean vector and covariance matrix of TVAR model coefficients are derived and [15] uses wavelets for TVAR model identification. Essentially, all these approaches provide maximum likelihood estimates of coefficients for TVAR models without measurement noise. In [16], autoregressive parameters were estimated from noisy observations by using a recursive least-squares adaptive filter.
Bayesian evidence can be used to track the optimal AR model order or more generally, to compare the performance of a TVAR model to an alternative model. To date,
Bayesian parameter tracking in AR models has mostly been achieved by Monte Carlo sampling methods [17–21].
In this paper, we develop a low-complexity variational message passing-based (VMP) realization for tracking of states, parameters and free energy (an upper bound on Bayesian evidence) in TVAR models. All update formulas are closed-form and the complete inference process can easily be realized.
VMP is a low-complexity distributed message passing-based realization of variational Bayesian inference on a factor graph representation of the model [22,23]. Previous work on message passing-based inference for AR models include [24], but their work describes maximum likelihood estimation and therefore does not track proper posteriors and free energy. In [25], variational inference is employed to estimate the parameters of a multivariate AR model, but their work does not take advantage of the factor graph representation.
The factor graph representation that we employ in this paper provides some distinct advantages from other works on inference in TVAR models. First, a factor graph formulation is by definition completely modular and supports re-using the derived TVAR inference equations as a plug-in module in other factor graph-based models. In particular, since we allow for measurement noise in the TVAR model specification, the proposed TVAR factor can easily be used as a latent module at any level in hierarchical dynamical models. Moreover, due to the modularity, VMP update rules can easily be mixed with different update schemes such as belief propagation and expectation [26,27] in other modules, leading to hybrid message passing schemes for efficient inference in complex models. We have implemented the TVAR model in the open source and freely available factor graph toolbox ForneyLab [28].
The rest of this paper is organized as follows. In Section 2, we specify the TVAR model as a probabilistic state space model and define the inference tasks that relate to tracking of states, parameters, and Bayesian evidence. After a short discussion on the merits of using Bayesian evidence as a model performance criterion (Section 3.1), we formulate Bayesian inference in the TVAR model as a set of sequential prediction-correction processes (Section 3.2). We will realize these processes as VMP update rules and proceed with a short review of Forney-style factor graphs and message passing in Section 4. Then, in Section 5, the VMP equations are worked out for the TVAR model and summarized in Table 1. Section 6 discusses a verification experiment on a synthetic data set and applications of the proposed TVAR model to temperature prediction and speech enhancement problems. Full derivations of the closed-form VMP update rules are presented in Appendix A.
In this section, we first specify TVAR as a state-space model. This is followed by an inference problem formulation.
{{
We use an ‘asymmetric’ definition (AR additions are positive) of the AR process as described in https://en.wikipedia.org/wiki/Autoregressive_moving-average_model:
\[ s_t = \sum_{m=1}^M \theta_{m,t}s_{t-m} + \epsilon_t \]
where
\(M\) is the order of the autoregression process
\(s_{t-m} \in \mathbb{R}\) is state component \(m\) at time \(t\)
\(\theta_{m,t} \in \mathbb{R}\) is parameter component \(m\) at time \(t\)
\(\epsilon_t \in \mathbb{R}\) is independent and identically distributed (i.i.d.) white noise at time \(t\)
}}
A TVAR model is specified as
\[ \begin{aligned} \theta_{m, t} &\sim \mathcal{N}\left(\theta_{m, t-1},\, \omega\right) \qquad\qquad\qquad (1a)\\ s_t &\sim \mathcal{N}\left(\sum_{m=1}^M \theta_{m, t} s_{t-m},\, \gamma^{-1}\right) \qquad(1b)\\ x_t &\sim \mathcal{N}\left(s_t,\, \tau^{-1}\right) \qquad\qquad\qquad (1c),\qquad \{\{ \text{invert tau} \}\} \end{aligned} \]
where
\(M\) is the order of the autoregression process
\(x_t \in \mathbb{R}\)
\(s_t \in \mathbb{R}\) is the state
\(\theta_{m,t} \in \mathbb{R}\) is the parameters at time \(t\), respectively.
\(M\) denotes the order of the AR model. As a notational convention, we use \(\mathcal{N}(\boldsymbol\mu, \mathbf\Sigma)\) to denote a Gaussian distribution with mean \(\boldsymbol\mu\) and co-variance matrix \(\mathbf\Sigma\). We can re-write (1) in state-space form as
\[\begin{aligned} & \boldsymbol{\theta}_t \sim \mathcal{N}\left(\boldsymbol{\theta}_{t-1},\, \omega \mathbf{I}_M\right)\qquad\qquad\qquad(\text{2a})\\ & \mathbf{s}_t \sim \mathcal{N}\left(\mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1},\, \mathbf{V}(\gamma)\right)\qquad(\text{2b}) \\ & x_t \sim \mathcal{N}\left(\mathbf{v}^{\top}_u \mathbf{s}_t,\, \tau^{-1}\right),\qquad\qquad\qquad(\text{2c}) \end{aligned}\]where
\(\boldsymbol{\theta}_t = (\theta_{1,t}, \theta_{2,t}, ..., \theta_{M,t})^\top\) \(\qquad\text{\{\{changed to this\}\}}\)
\(\mathbf{s}_{t-1} = (s_{t-1}, s_{t-2}, ..., s_{t-M})^\top\)
\(\mathbf{s}_t = (s_{t}, s_{t-1}, ..., s_{t-M+1})^\top\)
\(\mathbf{v}_u = (1,0,...,0)^\top\) is an \(M\)-dimensional unit vector,
\(\mathbf{V}(\gamma) = (1/ \gamma)\mathbf{v}_u\mathbf{v}_u^\top\)
\(\mathbf{B}(\boldsymbol{\theta}_t) = \begin{bmatrix} \boldsymbol{\qquad\theta}_t^\top \\ \mathbf{I}_{M-1} & \mathbf{0} \end{bmatrix} \qquad\qquad\qquad\text{(3)}\)
{{
For example, with \(M=5\) and
\[ \boldsymbol{\theta}_{t} = \begin{bmatrix} 0.30 \\ 0.31 \\ 0.32 \\ 0.33 \\ 0.34 \\ \end{bmatrix}, \qquad \mathbf{s}_{t-1} = \begin{bmatrix} 0.07 \\ 0.17 \\ 0.27 \\ 0.37 \\ 0.47 \\ \end{bmatrix} \]
\[ \begin{aligned} \mathbf{B}(\boldsymbol{\theta}_t) \cdot \mathbf{s}_{t-1} &= \begin{bmatrix} \theta_{1,t} & \theta_{2,t} & \theta_{3,t} & \theta_{4,t} & \theta_{5,t} \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{bmatrix} \cdot \begin{bmatrix} s_{t-1} \\ s_{t-2} \\ s_{t-3} \\ s_{t-4} \\ s_{t-5} \\ \end{bmatrix} \\ &= \begin{bmatrix} 0.30 & 0.31 & 0.32 & 0.33 & 0.34 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{bmatrix} \cdot \begin{bmatrix} 0.07 \\ 0.17 \\ 0.27 \\ 0.37 \\ 0.47 \\ \end{bmatrix} \\ &= \begin{bmatrix} 0.30\cdot0.07+0.31\cdot0.17+0.32\cdot0.27\cdot0.33\cdot0.37+0.34\cdot0.47 \\ 0.07 \\ 0.17 \\ 0.27 \\ 0.37 \\ \end{bmatrix} \\ &= \begin{bmatrix} 0.442 \\ 0.07 \\ 0.17 \\ 0.27 \\ 0.37 \\ \end{bmatrix} \end{aligned} \]
}}
{{
NOTE ABOUT Eq. 2(b):
\(\mathbf{s}_t \sim \mathcal{N}\left(\mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}, \mathbf{V}(\gamma)\right)\qquad(\text{2b})\)
This is a multivariate normal distribution. The mean is a vector, given by \(\mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}\). The covariance matrix is given by \(\mathbf{V}(\gamma)\). It has a single non-zero value in the top-left corner. This top-left value is the variance of the first element of the random vector \(\mathbf{s}_t\). Because the remaining diagonal values are zero, it means that the variance of all he remaining vector elements are zero, i.e. they are deterministic. The paper’s associated code does not implement this equation with a MvNormal distribution. I think the reason is due to a problem with positive definiteness. I tried to implement this equation with a MvNormal distribution but got an error in this regard (see below). In any case, you only use the FIRST element of the mean vector \(\mathbf{B}(\boldsymbol{\theta}_t) \mathbf{s}_{t-1}\) and the TOP-LEFT value of the covariance matrix \(\mathbf{V}(\gamma)\) to draw a random sample. The rest of the transition matrix simply selects the appropriate new state values.
}}
versioninfo() ## Julia versionJulia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)
Environment:
JULIA_NUM_THREADS = import Pkg
Pkg.add(Pkg.PackageSpec(;name="RxInfer", version="3.6.0"))
Pkg.add(Pkg.PackageSpec(;name="Distributions"))
Pkg.add(Pkg.PackageSpec(;name="Random"))
using Random, Distributions, LinearAlgebra
using RxInfer Resolving package versions...
No Changes to `~/.julia/environments/v1.10/Project.toml`
No Changes to `~/.julia/environments/v1.10/Manifest.toml`
Resolving package versions...
No Changes to `~/.julia/environments/v1.10/Project.toml`
No Changes to `~/.julia/environments/v1.10/Manifest.toml`
Resolving package versions...
No Changes to `~/.julia/environments/v1.10/Project.toml`
No Changes to `~/.julia/environments/v1.10/Manifest.toml`Pkg.status()Status `~/.julia/environments/v1.10/Project.toml`
[31c24e10] Distributions v0.25.111
[86711068] RxInfer v3.6.0
[9a3f8284] Randomrng = MersenneTwister(123)
M = 5
γ̃ = 1.7
θₜ = [0.30; 0.31; 0.32; 0.33; 0.34]
vᵤ = Vector(vcat(1.0, zeros(M-1)))
sₜ₋₁ = [0.07; 0.17; 0.27; 0.37; 0.47]5-element Vector{Float64}:
0.07
0.17
0.27
0.37
0.47Iₘ₋₁ = Diagonal(ones(M-1))4×4 Diagonal{Float64, Vector{Float64}}:
1.0 ⋅ ⋅ ⋅
⋅ 1.0 ⋅ ⋅
⋅ ⋅ 1.0 ⋅
⋅ ⋅ ⋅ 1.0V0 = Vector(zeros(M-1)) ##vector of zeros4-element Vector{Float64}:
0.0
0.0
0.0
0.0transpose(θₜ)1×5 transpose(::Vector{Float64}) with eltype Float64:
0.3 0.31 0.32 0.33 0.34hcat(Iₘ₋₁, V0) 4×5 Matrix{Float64}:
1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0BIθₜI = vcat(transpose(θₜ), hcat(Iₘ₋₁, V0))5×5 Matrix{Float64}:
0.3 0.31 0.32 0.33 0.34
1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0sₜ₋₁5-element Vector{Float64}:
0.07
0.17
0.27
0.37
0.47BIθₜI*sₜ₋₁5-element Vector{Float64}:
0.442
0.07
0.17
0.27
0.37## for sake of interest, the first element can be selected like this:
transpose(vᵤ)*BIθₜI*sₜ₋₁0.442## for sake of interest, a distribution can be created like this:
γ̃_std = sqrt(inv(γ̃))
Normal(transpose(vᵤ)*BIθₜI*sₜ₋₁, γ̃_std)Normal{Float64}(μ=0.442, σ=0.7669649888473704)V = (1/γ̃)*vᵤ*transpose(vᵤ) ## covariance matrix5×5 Matrix{Float64}:
0.588235 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0## for sake of interest, the top-left element can be selected like this:
transpose(vᵤ) * V * vᵤ0.5882352941176471## Try to create a MvNormal distribution. The error seems to be due
## to the lack of data in the matrix.
d = MvNormal(BIθₜI*sₜ₋₁, V)PosDefException: PosDefException: matrix is not positive definite; Cholesky factorization failed.## When I add a tiny value in the other diagonal elements, it works:
for i in 2:M
V[i, i] = tiny
end
V5×5 Matrix{Float64}:
0.588235 0.0 0.0 0.0 0.0
0.0 1.0e-12 0.0 0.0 0.0
0.0 0.0 1.0e-12 0.0 0.0
0.0 0.0 0.0 1.0e-12 0.0
0.0 0.0 0.0 0.0 1.0e-12d = MvNormal(BIθₜI*sₜ₋₁, V)FullNormal(
dim: 5
μ: [0.442, 0.07, 0.17, 0.27, 0.37]
Σ: [0.5882352941176471 0.0 … 0.0 0.0; 0.0 1.0e-12 … 0.0 0.0; … ; 0.0 0.0 … 1.0e-12 0.0; 0.0 0.0 … 0.0 1.0e-12]
)## now we can draw a vector sample from this Multivariate distribution:
rand(rng, d)5-element Vector{Float64}:
1.354893792066023
0.0700020481797078
0.17000114265090288
0.2700004594156204
0.3699996033209207Notice how only \(\mathbf{V}(\gamma)[1,1]\) has a distribution. The remaining components from the \(\mathbf{B}(\boldsymbol{\theta}_t) \mathbf{s}_{t-1}\) vector each have a practically \(0\) variance and keep their values.
Technically, a TVAR model usually assumes \(\tau = 0\) indicating that there is no measurement noise. Note that the presence of measurement noise in (2c) “hides” the states \(s_t\) in the generative model (2) from the observation sequence \(x_t\), yielding a latent TVAR. We add measurement noise explicitly, so the model is able to accept information from likelihood functions that are not constrained to be delta functions with hard observations. As a result, the AR model that we define here can be used at any level in deep hierarchical structures such as [29] as a plug-in module.
In a time-invariant AR model, \(\boldsymbol\theta\) are part of the parameters of the system. In a timevarying AR model, we consider \(\boldsymbol\theta_t\) and \(\mathbf{s}_t\) together the set of time-varying states. The parameters of the TVAR model are \(\{\boldsymbol{\theta}_0, \mathbf{s}_0, \omega, \gamma, \tau \}\).
At the heart of the TVAR model is the transition model (2b), where \(\mathbf{B}(\boldsymbol{\theta}_t)\) is a companion matrix with AR coefficients. The multiplication \(\mathbf{B}(\boldsymbol{\theta})\mathbf{s}_{t-1}\) performs two operations: a dot product \(\boldsymbol{\theta}^\top_t\mathbf{s}_{t-1}\) and a vector shift of \(\mathbf{s}_{t-1}\) by one time step. The latter operation can be interpreted as bookkeeping, as it shifts each entry of \(\mathbf{s}_{t-1}\) one position down and discards \(s_{t-M}\).
For a given time series \(x_: = (x_1, x_2, . . . , x_T)\), we are firstly interested in recursively updating posteriors for the states \(p(\mathbf{s}_t \mid x_{1:l})\) and \(p(\boldsymbol\theta_t \mid x_{1:l})\). In this context, prediction, filtering and smoothing are recovered for \(l < t\), \(l = t\) and \(l > t\), respectively.
Furthermore, we are interested in computing posteriors for the parameters \(p(\boldsymbol{\theta}_0 \mid x)\), \(p(\mathbf{s}_0 \mid x)\), \(p(\omega \mid x)\), \(p(\gamma \mid x)\) and \(p(\tau \mid x)\).
Finally, we are interested in scoring the performance of a proposed TVAR model \(m\) with specified priors for the parameters. In this paper, we take a full Bayesian approach and select Bayesian evidence \(p(x_: \mid m)\) as the performance criterion. Section 3.1 discusses the merits of Bayesian evidence as a model performance criterion.
In this section, we first discuss some of the - merits of using Bayesian evidence as a model performance criterion. This is followed by an exposition of - how to compute Bayesian evidence and the desired posteriors in the TVAR model.
Consider a model \(m\) with parameters \(\boldsymbol{\theta}\) and observations \(x\). Bayesian evidence \(p(x_: \mid m)\) is considered an excellent model performance criterion. Note the following decomposition [30]:
\[ \begin{aligned} \log p(x_: \mid m) & =\log \frac{p(x_: \mid \boldsymbol{\theta}, m) p(\boldsymbol{\theta} \mid m)}{p(\boldsymbol{\theta} \mid x_:, m)} \quad \text { (use Bayes rule) } \\ & =\int p(\boldsymbol{\theta} \mid x_:, m) \cdot \underbrace{\log \frac{p(x_: \mid \boldsymbol{\theta}, m) p(\boldsymbol{\theta} \mid m)}{p(\boldsymbol{\theta} \mid x_:, m)}}_{\log p(x_: \mid m) \text { is not a function of } \boldsymbol{\theta}} \mathrm{d} \boldsymbol{\theta} \\ & =\underbrace{\int p(\boldsymbol{\theta} \mid x_:, m) \log p(x_: \mid \boldsymbol{\theta}, m) \mathrm{d} \boldsymbol{\theta}}_{\text {data fit }}-\underbrace{\int p(\boldsymbol{\theta} \mid x_:, m) \log \frac{p(\boldsymbol{\theta} \mid x_:, m)}{p(\boldsymbol{\theta} \mid m)} \mathrm{d} \boldsymbol{\theta}}_{\text {complexity }}\qquad(4) \end{aligned} \]
The first term (data fit or sometimes called accuracy) measures how well the model predicts the data \(x_:\), after having learned from the data. We want this term to be large (although only focusing on this term could lead to over-fitting). The second term (complexity) quantifies the amount of information that the model absorbed through learning by moving parameter beliefs from \(p(\boldsymbol{\theta}\mid m)\) to \(p(\boldsymbol{\theta}\mid x_:,m)\). To see this, note that the mutual information between two variables \(\boldsymbol{\theta}\) and \(x_:\), which is defined as
\[ I[\boldsymbol{\theta} ; x_:]=\iint p(\boldsymbol{\theta}, x_:) \log \frac{p(\boldsymbol{\theta} \mid x_:)}{p(\boldsymbol{\theta})} \mathrm{d} \boldsymbol{\theta} \mathrm{d} x_: \]
can be interpreted as expected complexity. The complexity term regularizes the Bayesian learning process automatically. Preference for models with high Bayesian evidence implies a preference for models that get a good data fit without the need to learn much from the data set. These types of models are said to generalize well, since they can be applied to different data sets without specific adaptations for each data set. Therefore, Bayesian learning automatically leads to models that tend to generalize well.
Note that Bayesian evidence for a model \(m\), given a full time series \(x_: = (x_1, x_2, ..., x_T)\), can be computed by multiplication of the sample-based evidences:
\[ p(x_: \mid m)=\prod_{t=1}^T p\left(x_t \mid x_{1:t-1}, m\right) \tag{5} \]
To illustrate the type of calculations that are needed for computing Bayesian model evidence and the posteriors for states and parameters, we now proceed to write out the needed calculations for the TVAR model in a filtering context.
Assume that at the beginning of time step \(t\), we are given the state posteriors \(q(\mathbf{s}_{t-1}\mid x_{1:t-1})\), \(q(\boldsymbol{\theta}_{t-1}\mid x_{1:t-1})\). We will denote the inferred probabilities by \(q(\cdot)\), in contrast to factors from the generative model that are written as \(p(\cdot)\). We start the procedure by setting the state priors for the generative model at step \(t\) to the posteriors of the previous time step
\[ p\left(\mathbf{s}_{t-1} \mid x_{1:t-1}\right):=q\left(\mathbf{s}_{t-1} \mid x_{1:t-1}\right) \tag{6} \] \[ p\left(\boldsymbol{\theta}_{t-1} \mid x_{1:t-1}\right):=q\left(\boldsymbol{\theta}_{t-1} \mid x_{1:t-1}\right) \tag{7} \]
Given a new observation \(x_t\), we are now interested in inferring the evidence \(q(x_t \mid x_{t-1})\), and in inferring posteriors \(q(\mathbf{s}_t \mid x_{1:t})\) and \(q(\boldsymbol{\theta}_t \mid x_{1:t})\).
This involves a prediction (forward) pass {{‘generence’}} through the system that leads to the evidence update, followed by a correction (backward) pass {{inference}} that updates the states. We work this out in detail below. For clarity of exposition, in this section we call \(\mathbf{s}_t\) states and \(\boldsymbol\theta_t\) parameters. Starting with the forward pass (from latent variables toward observation) {{‘generence’}}, we first compute a parameter prior predictive:
\[ \underbrace{q\left(\boldsymbol{\theta}_t \mid x_{1:t-1}\right)}_{\substack{\text{parameter} \\ \text{prior predictive}}} =\int \underbrace{p\left(\boldsymbol{\theta}_t \mid \boldsymbol{\theta}_{t-1}\right)}_{\substack{\text{parameter}\\\text{transition}}} \underbrace{p\left(\boldsymbol{\theta}_{t-1} \mid x_{1:t-1}\right)}_{\substack{\text { parameter } \\ \text { prior }}} \mathrm{d} \boldsymbol{\theta}_{t-1}. \tag{8} \]
Then the prior predictive for the state transition becomes:
\[ \underbrace{q\left(\mathbf{s}_t \mid \mathbf{s}_{t-1,} x_{1:t-1}\right)}_{\substack{\text{state transition}\\\text { prior predictive }}} =\int \underbrace{p\left(\mathbf{s}_t \mid x_{t-1}, \boldsymbol{\theta}_t\right)}_{\text {state transition }} \underbrace{q\left(\boldsymbol{\theta}_t \mid x_{1:t-1}\right)}_{\substack{\text { parameter } \\\text { prior predictive }}} \mathrm{d} \boldsymbol{\theta}_t. \tag{9} \]
Note that the state transition prior predictive, due to its dependency on time-varying \(\boldsymbol\theta_t\), is a function of the observed data sequence. The state transition prior predictive can be used together with the state prior to inferring the state prior predictive:
\[ \underbrace{q\left(\mathbf{s}_t \mid x_{1:t-1}\right)}_{\substack{\text { state } \\\text { prior predictive }}} =\int \underbrace{q\left(\mathbf{s}_t \mid \mathbf{s}_{t-1}, x_{1:t-1}\right)}_{\substack{ \text { state transition } \\ \text { prior predictive }}} \underbrace{p\left(\mathbf{s}_{t-1} \mid x_{1:t-1}\right)}_{\text {state prior }} \mathrm{d} \mathbf{s}_{t-1}. \tag{10} \]
The evidence for model \(m\) that is provided by observation \(x_t\) is then given by
\[ \underbrace{q\left(x_t \mid x_{1:t-1}\right)}_{\text {evidence }}=\int \underbrace{p\left(x_t \mid \mathbf{s}_t\right)}_{\substack{\text { state } \\\text { likelihood }}} \underbrace{q\left(\mathbf{s}_t \mid x_{1:t-1}\right)}_{\substack{\text { state prior } \\\text { predictive }}} \mathrm{d} \mathbf{s}_t. \tag{11} \]
When \(x_t\) has not yet been observed, \(q(x_t \mid x_{1:t-1})\) is a prediction for \(x_t\). After plugging in the observed value for \(x_t\), the evidence is a scalar that scores how well the model performed in predicting \(x_t\). As discussed in (5), the results \(q(x_t \mid x_{1:t-1})\) for \(t = 1, 2, ..., T\) in (11) can be used to score the model performance for a given time series \(x = (x_1, x_2, ..., x_T)\). Note that to update the evidence, we need to integrate over all latent variables \(\boldsymbol{\theta}_{t-1},\boldsymbol{\theta}_{t},\mathbf{s}_{t-1}\) and \(\mathbf{s}_t\) (by (8)–(11)). In principle, this scheme needs to be extended with integration over the parameters \(\omega\), \(\gamma\) and \(\tau\).
Once we have inferred the evidence, we proceed by a backward corrective pass through the model to update the posterior over the latent variables given the new observation \(x_t\). The state posterior can be updated by Bayes rule:
\[ \underbrace{ q\left(\mathbf{s}_t \mid x_{1:t}\right)}_{\text {state posterior }} =\frac { \overbrace{p\left(x_t \mid \mathbf{s}_t\right)}^ {\substack{\text{state}\\ \text{likelihood}}} \overbrace{q\left(\mathbf{s}_t \mid x_{1:t-1}\right)}^ {\substack{\text{state prior}\\ \text{predictive}}} } { \underbrace{q\left(x_t \mid x_{1:t-1}\right)}_{\text{evidence}} } \tag{12} \]
Next, we need to compute a likelihood function for the parameters. Fortunately, we can re-use some intermediate results from the forward pass. The likelihood for the parameters is given by
\[ \underbrace{ q\left(x_t \mid \boldsymbol{\theta}_t, x_{1:t-1}\right)}_{ \substack{\text{parameter}\\\text{likelihood}}} =\int \underbrace{p\left(x_t \mid \mathbf{s}_t\right)}_{ \substack{\text { state } \\\text { likelihood }}} \underbrace{ q\left(\mathbf{s}_t \mid \boldsymbol{\theta}_t, x_{1:t-1}\right)}_{ \substack{\text { state prior } \\ \text { predictive }}} \mathrm{d} \mathbf{s}_t \tag{13} \]
The parameter posterior then follows from Bayes rule:
\[ \underbrace{ q\left(\boldsymbol{\theta}_t \mid x_{1:t}\right)}_{\text {parameter posterior }} = \frac {\overbrace{ q\left(x_t \mid \boldsymbol{\theta}_t, x_{1:t-1}\right)}^{\substack{\text{parameter}\\\text{likelihood}}} \overbrace{ q\left(\boldsymbol{\theta}_t \mid x_{1:t-1}\right)}^{\substack{\text{parameter}\\\text{prior predictive}}} } {\underbrace{ q\left(x_t \mid x_{1:t-1}\right)}_{\text{evidence}} } \tag{14} \]
Equations (11), (12) and (14) contain the solutions to our inference task. Note that the evidence \(q(x_t)|x_{1:t-1}\) is needed to normalize the latent variable posteriors in (12) and (14). Moreover, while we integrate over the states by (11) to compute the evidence, (14) reveals that the evidence can alternatively be computed by integrating over the parameters through
\[ \underbrace{ q\left(x_t \mid x_{1:t-1}\right)}_{ \text {evidence }} =\int \underbrace{ q\left(x_t \mid \boldsymbol{\theta}_t, x_{1:t-1}\right)}_{ \substack{\text { parameter } \\\text { likelihood }}} \underbrace{ q\left(\boldsymbol{\theta}_t \mid x_{1:t-1}\right)}_{\substack{\text { parameter } \\\text { prior predictive }}} \mathrm{d} \boldsymbol{\theta}_t. \tag{15} \]
This latter method of evidence computation may be useful if re-using (11) in (14) leads to numerical rounding issues.
Unfortunately, many of Equations (8) through (14) are not analytically tractable for the TVAR model. This happens due to (1) integration over large state spaces, (2) non-conjugate prior-posterior pairing, and (3) the absence of a closed-form solution for the evidence factor.
To overcome this challenge, we will perform inference by a hybrid message passing scheme in a factor graph. In the next section, we give a short review of Forney-Style Factor Graphs (FFG) and Message-Passing (MP) based inference techniques.
In this section, we make a brief introduction of Forney-Style Factor graph (FFG) and sum-product (SP) algorithm. After that we review the minimization of variational free energy and Variational Message Passing (VMP) algorithm.
A Forney-style Factor graph is a representation of a factorized function where the factors and variables are represented by nodes and edges, respectively. An edge is connected to a node if and only if the (edge) variable is an argument of the node function. In our work, we use FFGs to represent factorized probability distributions. FFGs provide both an attractive visualization of the model and a highly efficient and modular inference method based on message passing. An important component of the FFG representation is the equality node. If a variable \(x\) is shared between more than two nodes, then we introduce two auxiliary variables \(x'\) and \(x''\) and use an equality node
\[ f_=\left(x, x^{\prime}, x^{\prime \prime}\right)=\delta\left(x-x^{\prime}\right) \delta\left(x-x^{\prime \prime}\right) \tag{16} \]
to constrain the marginal beliefs over \(x\), \(x'\), \(x''\) to be equal. With this mechanism, any factorized function can be represented as an FFG.
An FFG visualization of the TVAR model is depicted in Figure 3, but for illustrative purposes, we first consider an example factorized distribution
\[ p\left(x_1, x_2, x_3, x_4\right)=p\left(x_1\right) p\left(x_2 \mid x_1\right) p\left(x_3 \mid x_2\right) p\left(x_4 \mid x_3\right) \tag{17} \]
This distribution can be visualized by an FFG shown in Figure 1. An FFG is in principle an undirected graph but we often draw arrows on the edges in the “generative” direction, which is the direction that describes how the observed data is generated. Assume that we are interested in computing the marginal for \(x_2\), given by
\[ p\left(x_2\right)=\iiint p\left(x_1, x_2, x_3, x_4\right) \mathrm{d} x_1 \mathrm{~d} x_3 \mathrm{~d} x_4 \tag{18} \]
We can reduce the complexity of computing this integral by rearranging the factors over the integration signs as
\[ \begin{aligned} p\left(x_2\right) & =\underbrace{\int \underbrace{p\left(x_1\right)}_{\vec{\mu}_1\left(x_1\right)} p\left(x_2 \mid x_1\right) \mathrm{d} x_1}_{\vec{\mu}_2\left(x_2\right)} \cdot \underbrace{(\int p\left(x_3 \mid x_2\right) \underbrace{\left.\left(\int p\left(x_4 \mid x_3\right) \mathrm{d} x_3\right)\right)}_{\overleftarrow{\mu}_3\left(x_3\right)} \mathrm{d} x_3)}_{\overleftarrow{\mu}_2\left(x_2\right)} \qquad{\text{(19a)}}\\ & =\overrightarrow{\mu}_2\left(x_2\right) \cdot \overleftarrow{\mu}_2\left(x_2\right). \qquad (\text{19b}) \end{aligned} \]
These re-arranged integrals can be interpreted as messages that are passed over the edges, see Figure 1. It is a notational convention to call a message \(\overrightarrow{\mu}(\cdot)\) that aligns with the direction of the edge arrow a forward message and similarly, a message \(\overleftarrow{\mu}(\cdot)\) that opposes the direction of the edge is called a backward message.
This message passed-based algorithm of computing the marginal is called belief propagation (BP) or the sum-product algorithm. As can be verified in (19), for a node with factor \(f(y, x_1, . . . , x_n)\), the outgoing BP message \(\overrightarrow{\mu}(y)\) to variable \(y\) can be expressed as
\[ \overrightarrow{\mu}_y(y)=\int \cdots \int f\left(y, x_1, \ldots, x_n\right) \prod_{i=1}^n \overrightarrow{\mu}_i\left(x_i\right) \mathrm{d} x_i. \tag{20} \]
where \(\overrightarrow{\mu}_i(x_i)\) is an incoming message over edge \(x_i\). If the factor graph is a tree, meaning that the graph contains no cycles, then BP leads to exact Bayesian inference. A more detailed explanation of belief propagation message passing in FFGs can be found in [26].
Technically, BP is a message passing algorithm that belongs to a family of message passing algorithms that minimize a constrained variational free energy functional [31]. Unfortunately, the sum-product rule (20) only has a closed-form solution for Gaussian incoming messages \(\overrightarrow{\mu}_i(x_i)\) and linear variable relations in \(f(x, s_1, . . . , s_n)\). Another important member of the free energy minimizing algorithms is the Variational Message Passing (VMP) algorithm [22]. VMP enjoys a wider range of analytically computable message update rules.
We shortly review variational Bayesian inference and VMP next. Consider a model \(p(x_:, \mathbf{s}_:)\) with observations \(x_:\) and unobserved (latent) variables \(\mathbf{s}_:\). We are interested in inferring the posterior distribution \(p(\mathbf{s}_:\mid x_:)\). In variational inference we introduce an approximate posterior \(q(\mathbf{s}_:)\) and define a variational free energy functional as
\[ F[q] \triangleq \int q(\mathbf{s}_:) \log \frac{q(\mathbf{s}_:)}{p(x_:, \mathbf{s}_:)} \mathrm{d}\mathbf{s}_:=\underbrace{\int q(\mathbf{s}_:) \log \frac{q(\mathbf{s}_:)}{p(\mathbf{s}_: \mid x_:)} \mathrm{d} \mathbf{s}_:}_{\text {KL divergence } D_{\mathrm{KL}}(q, p)}-\underbrace{\log p(x_:)}_{\text {log-evidence }} \tag{21} \]
The second term in (21) (log-evidence) is not a function of the argument of \(F\). The first term is a KL-divergence, which is by definition non-negative and only equals zero for \(q(\mathbf{s}_:) = p(\mathbf{s}_: \mid x_:)\). As a result, variational inference by minimization of \(F[q]\) provides
\[ q^*(\mathbf{s}_:)=\arg \min _q F[q] \tag{22} \]
which is an approximation to the Bayesian posterior \(p(\mathbf{s}_: \cdot x_:)\). Moreover, the minimized free energy \(F[q^*]\) is an upper bound for minus log-evidence and in practice is used a model performance criterion. Similarly to (4), the free energy can be decomposed as
\[ F[q]= \underbrace{\int q(\mathbf{s}_:) \log p(x_: \mid \mathbf{s}_:, m) \mathrm{d} \mathbf{s}_:}_{\text {accuracy }} - \underbrace{\int q(\mathbf{s}_:) \log \underbrace{ \frac {q(\mathbf{s}_:)} {p(\mathbf{s}_: \mid m)} }_{\text{prior}} \mathrm{d} \mathbf{s}_:}_{\text {complexity }} \tag{23} \]
which underwrites its usage as a performance criterion for model \(m\), given observations \(x_:\).
In an FFG context, the model \(p(x_:, \mathbf{s}_:)\) is represented by a set of connected nodes. Consider a generic node of the FFG, given by \(f(x, s_1, . . ., s_n)\) where in the case of VMP, the incoming messages are approximations to the marginals \(q_i(s_i), i = 1, . . . , n\), see Figure 2.
It can be shown that the outgoing VMP message of \(f\) towards edge \(x\) is given by [32]
\[ \overrightarrow{v}(x) \propto \exp \left(\int \cdots \int \log f\left(x, s_1, \ldots, s_n\right) \prod_{i=1} q\left(s_i\right) \mathrm{d} s_i\right). \tag{24} \]
In this paper, we adopt the notational convention to denote belief propagation messages (computed by (20)) by \(\mu\) and VMP messages (computed by (24)) by \(\nu\). The approximate marginal \(q(x)\) can be obtained by multiplying incoming and outgoing messages on the edge for \(x\)
\[ q(x) \propto \overrightarrow{\nu}(x) \overleftarrow{\nu}(x). \tag{25} \]
This process (compute forward and backward messages for an edge and update the marginal) is executed sequentially and repeatedly for all edges in the graph until convergence. In contrast to BP-based inference, the VMP and marginal update rules (24) and (25) lead to closed-form expressions for a large set of conjugate node pairs from the exponential family of distributions. For instance, updating the variance parameter of a Gaussian node with a connected inverse-gamma distribution node results in closed-form VMP updates.
In short, both BP- and VMP-based message passing can be interpreted as minimizing variational free energy, albeit under a different set of local constraints [31]. Typical constraints include factorization and form constraints on the posterior such as \(q(\mathbf{s}_:) = \prod_i q_i(s_i)\) and \(q(\mathbf{s}_:) = \mathcal{N}(\mathbf{s}_:|\boldsymbol{\mu}, \mathbf\Sigma)\), respectively. Since the constraints are local, BP and VMP can be combined in a factor graph to create hybrid message passing-based variational inference algorithms. For a more detailed explanation of VMP in FFGs, we refer to [32]. Note that hybrid message passing does in general not guarantee to minimize variational free energy [33]. However, in our experiments in Section 6 we will show that iterating our stationary solutions by message passing does lead to free energy minimization.
In this section, we focus on deriving message passing-based inference in the TVAR model. We develop a TVAR composite factor for the FFG framework and specify the intractable BP messages around the TVAR node. Then we present a message passing-based inference solution.
The TVAR model at time step \(t\) can be represented by an FFG as shown in Figure 3. We are interested in providing a message passing solution to the inference tasks as specified by Equations (8)–(14). At the left-hand side of Figure 3, the incoming messages are the priors \(p(\boldsymbol{\theta}_{t-1}|x_{1:t-1})\) and \(p(\mathbf{s}_{t-1}|x_{1:t-1})\). At the bottom of the graph, there is a new observation \(x_t\). The goal is to pass messages in the graph to compute posteriors \(q(\boldsymbol\theta_t|x_{1:t})\) (message 16 ) and \(q(\mathbf{s}_t|x_{1:t})\) (message 11 ). In order to support smoothing algorithms, we also want to be able to pass incoming prior messages from the right-hand side to outgoing messages 13 and 18 at the left-hand side. Forward and backward messages are drawn as open and closed circles respectively.
Technically, the generative model (2) at time step \(t\) for the TVAR model can shortly be written as \(p(x_t|\mathbf{z}_t)p(\mathbf{z}_t)|\mathbf{z}_{t-1}\)), where \(\mathbf{z}_t = \{\mathbf{s}_t, \boldsymbol{\theta}_t, \omega, \gamma, \tau \}\) are the latent variables. On this view, we can write the free energy functional for the TVAR model at time step \(t\) as
\[ F\left[q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:t}\right)\right]=\iint q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:t}\right) \text{log} \frac {\overbrace{ q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:t}\right) }^{\text{posterior}} } { \underbrace{ p\left(x_t \mid \mathbf{z}_t\right) p\left(\mathbf{z}_t \mid \mathbf{z}_{t-1}\right)}_{\text{generative model}} \underbrace{ p\left(\mathbf{z}_{t-1} \mid x_{1:t-1}\right) }_{\text{prior from past}} } \text{d}\mathbf{z}_{t-1}\text{d}\mathbf{z}_t. \tag{26} \]
and minimize \(F[q]\) by message passing. In a smoothing context, we would include a prior from the future \(p(\mathbf{z}_t|x_{t+1:t+T}) := q(\mathbf{z}_t|x_{t+1:t+T})\), yielding a free energy functional
\[ F\left[q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:T}\right)\right]=\iint q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:T}\right) \text{log} \frac {\overbrace{ q\left(\mathbf{z}_{t-1}, \mathbf{z}_t \mid x_{1:T}\right) }^{\text{posterior}} } { \underbrace{ p\left(x_t \mid \mathbf{z}_t\right) p\left(\mathbf{z}_t \mid \mathbf{z}_{t-1}\right)}_{\text{generative model}} \underbrace{ p\left(\mathbf{z}_{t} \mid x_{t+1:t+T}\right) }_{\text{prior from future}} } \text{d}\mathbf{z}_{t-1}\text{d}\mathbf{z}_t. \tag{27} \]
In a filtering context, \(q(\mathbf{z}_t \mid x_{t+1:t+T}) \propto 1\) and the functional (27) simplifies to (26).
The modularity of message passing in FFGs allows us to focus on only the intractable message and marginal updates. For instance, while there is no problem with the analytical computation of the backward message 12 , the corresponding forward message 4 ,
\[ \overrightarrow{\mu}\left(\mathbf{s}_t\right) = \int \mathcal{N}\left(\mathbf{s}_t \mid \mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}, \mathbf{V}(\gamma)\right) \underbrace{\overrightarrow{\mu}\left(\mathbf{s}_{t-1}\right) \overrightarrow{\mu}\left(\boldsymbol{\theta}_t\right) \overrightarrow{\mu}(\gamma)}_{\text {Gaussian messages }} \mathrm{d} \gamma \mathrm{d} \boldsymbol{\theta}_t \mathbf{s}_{t-1} \tag{28} \]
cannot be analytically solved [34]. Similarly, some other messages 13 , 14 and 15 do not have a closed-form solution in the constrained free energy minimization framework. For purpose of identification, in Figure 3 intractable messages are marked in red color.
In an FFG framework, we can isolate the problematic part of the TVAR model (Figure 3) by introducing a “composite” AR node. Composite nodes conceal their internal operations from the rest of the graph. As a result, inference can proceed as long as each composite node follows proper message passing communication rules at its interfaces to the rest of the graph. The composite AR node
\[ f_{\mathrm{AR}}\left(\mathbf{s}_t, \mathbf{s}_{t-1}, \boldsymbol{\theta}_t, \gamma\right)=\mathcal{N}\left(\mathbf{s}_t \mid \mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}, \mathbf{V}(\gamma)\right) \tag{29} \]
is indicated in Figure 3 by a dashed box. Note that the internal shuffling of the parameters \(\boldsymbol{\theta}_t\) and \(\gamma\), respectively by means of \(B(\boldsymbol{\theta}_t)\) and \(\mathbf{V}(\gamma)\), is hidden from the network outside the composite AR node.
We isolate the composite AR node by the specification
\[ f_{\mathrm{AR}}(x_:, \mathbf{s}, \boldsymbol{\theta}, \gamma)=\mathcal{N}(x_: \mid \mathbf{B}(\boldsymbol{\theta}) \mathbf{s}, \mathbf{V}(\gamma)), \tag{30} \]
where, relative to (29), we used substitutions \(x_: = \mathbf{s}_t, \mathbf{s} = \mathbf{s}_{t-1}, \boldsymbol{\theta} = \boldsymbol{\theta}_t\).
Under the structural factorization constraint (See Appendix A.1 for more on structural VMP).
\[ q(x_:, \mathbf{s}, \boldsymbol{\theta}, \gamma)=q(x_:, \mathbf{s}) q(\boldsymbol{\theta}) q(\gamma), \tag{31} \]
and consistency constraints
\[ q(x_:)=\int q(x_:, \mathbf{s}) \mathrm{d} \mathbf{s}, \quad q(\mathbf{s})=\int q(x_:, \mathbf{s}) \mathrm{d} x_: \tag{32} \]
the marginals \(q(\boldsymbol{\theta}), q(\mathbf{s}), q(x_:)\) and \(q(\gamma)\) can be obtained from the minimisation of the composite-AR free energy functional
\[ F_{\mathrm{AR}}[q]=\int q(x_:, \mathbf{s}) q(\boldsymbol{\theta}) q(\gamma) \log \frac {\overbrace{ q(x_:, \mathbf{s}) q(\boldsymbol{\theta}) q(\gamma) }^{\text{posterior}}} {\underbrace{ f_{\mathrm{AR}}(x_:, \mathbf{s}, \boldsymbol{\theta}, \gamma) }_{\text{AR node}}} \mathrm{d} x_: \mathrm{d} \mathbf{s} \mathrm{d} \boldsymbol{\theta} \mathrm{d} \gamma. \tag{33} \]
Recalling (25), we can write the minimizer of FE functional (33) with respect to \(\boldsymbol{\theta}\) as
\[ q(\boldsymbol{\theta}) \propto \overrightarrow{\nu}(\boldsymbol{\theta}) \overleftarrow{\nu}(\boldsymbol{\theta}) \tag{34} \]
where \(q(\boldsymbol{\theta})\) is associated with the incoming message to AR node and \(\overrightarrow{\nu}(\boldsymbol{\theta})\) is a variational outgoing message. Hence, the outgoing message from the AR node toward \(\boldsymbol{\theta}\) can be written as
\[ \overrightarrow{\nu}(\boldsymbol{\theta}) \propto \exp \left(\mathbb{E}_{q(x_:, \mathbf{s}) q\left(\boldsymbol{\theta}_t\right) q(\gamma)} \log [\mathcal{N}(x_: \mid \mathbf{B} \mathbf{s}, \mathbf{V})]\right) \tag{35} \]
In Appendix A we work out a closed-form solution for this and all other update rules plus an evaluation of free energy for the composite AR node. The results are reported in Table 1. With these rules in hand, the composite AR node can be plugged into any factor graph and take part in a message passing-based free energy minimization process.
In this section, we first verify the proposed methodology by a simulation of the proposed TVAR model on synthetic data, followed by validation experiments on two real world problems. We implemented all derived message passing rules in the open source Julia package ForneyLab.jl [28]. The code for the experiments and for the AR node can be found in public Github repositories. (https://github.com/biaslab/TVAR_FFG, accessed on 27 May 2021, https://github.com/biaslab/LAR, accessed on 27 May 2021). We used the following computer configuration to run the experiments. Operation system: macOS Big Sur, Processor: 2,7 GHz Quad-Core Intel Core i7, RAM: 16 GB.
To verify the proposed TVAR inference methods, we synthesized data from two generative models \(m_1\) and \(m_2\), as follows:
{{
}}
\[ \begin{aligned} & \boldsymbol{\theta}_t \sim \begin{cases}\delta\left(\boldsymbol{\theta}_t-\boldsymbol{\theta}_{t-1}\right) & \text { if } m=m_1 \qquad\text{\{\{AR\}\}} \qquad(\text{36a})\\ \mathcal{N}\left(\boldsymbol{\theta}_{t-1},\, \omega \mathbf{I}_M\right) & \text { if } m=m_2 \qquad\text{\{\{TVAR\}\}} \end{cases} \\ & \boldsymbol{s}_t \sim \mathcal{N}\left(B\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}, \mathbf{V}(\gamma)\right) \qquad(\text{36b})\\ & x_t \sim \mathcal{N}\left(\mathbf{v}^T_u \mathbf{s}_t, \tau\right) \qquad(\text{36c}) \end{aligned} \]
with priors
\[ \begin{aligned} p(M=k) & =\prod_{k=1}^{10} 0.1^{[M=k]} \qquad(\text{37a})\\ \boldsymbol{\theta}_0 & \sim \begin{cases}\mathcal{N}(\mathbf{0}, \mathbf{I}) & \text { if } m=m_1 \\ \mathcal{N}(\mathbf{0}, 1\mathrm{e}12 \,\mathbf{I}) & \text { if } m=m_2\end{cases} \qquad(\text{37b})\\ \mathbf{s}_0 & \sim \mathcal{N}(\mathbf{0}, 1\mathrm{e}12 \,\mathbf{I}) \qquad(\text{37c})\\ \gamma & \sim \mathrm\Gamma(1.0, 1\mathrm{e}-5) \qquad(\text{37d}) \\ \tau & =1.0 \qquad(\text{37e})\\ \omega & =0.01 \qquad(\text{37f}) \end{aligned} \]
where \(M\) is the number of AR coefficients. Although these models differ only with respect to the properties of the AR coefficients \(\boldsymbol{\theta}\), this variation has an important influence on the data generative process.
One-time segment FFGs corresponding to the Equation (36) are depicted in Figure 4.
For each model, we generated a data set of 100 different time series, each of length 100 (so we have \(2\times100\times100\) data points). For each time series, as indicated by (37a), the AR order \(M\) of the generative process was randomly drawn from the set \(\{1, 2, . . . , 10\}\). We used rather non-informative/broad priors for states and parameters for both models, see (37). This was done to ensure that the effect of the prior distributions is negligible relative to the information in the data set.
These time series were used in the following experiments. We selected two recognition models \(m_1\) and \(m_2\) with the same specifications as were used for generating the data set. The recognition models were trained on time series that were generated by models with the same AR order.
We proceeded by computing the quantities \(q(\mathbf{s}_{1:T}|x_{1:T})\), \(q(\boldsymbol{\theta}_{1:T}|x_{1:T})\), \(q(\gamma|x_{1:T})\) and \(F[q(\mathbf{z}_{t-1}|x_{1:T})]\) (where \(\mathbf{z}\) comprises all latent states and parameters) for both models, following the proposed rules from Table 1.
As a verification check, we first want to ensure that inference recovers the hidden states \(s_t\) for each \(t \in (1, 2, . . ., 100)\). Secondly, we want to verify the convergence of FE. As we have not used any approximations along the derivations of variational messages, we expect a smoothly decreasing curve for FE until convergence. The results of the verification stage are highlighted for a typical case in Figure 5. The figure confirms that states \(s_t\) are accurately tracked and that a sliding average of the AR coefficients \(\boldsymbol{\theta}_t\) is also nicely tracked. Figure 5 also indicates that the FE uniformly decreases towards lower values as we spend more computational power.
We note that the FE score by itself does not explain whether the model is good or not, but it serves as a good measure for model comparison. In the following subsection, we demonstrate how FE scores can be used for model selection.
AR models are well-known for predicting different weather conditions such as wind, temperature, precipitation, etc. Here, we will revisit the problem of modeling daily temperature. We used a data set of daily minimum temperatures (in \(\degree\) C) in Melbourne, Australia, 1981–1990 (3287 days) (https://www.kaggle.com/paulbrabban/daily-minimumtemperatures-in-melbourne, accessed on 27 May 2021). We then corrupted the data set by adding random noise sampled from \(\mathcal{N}(0, 10.0)\) to the actual temperatures. A fragment of the time-series is depicted in Figure 6.
To estimate the actual temperature based on past noisy observations by computing \(q(\mathbf{s}_t|x_{1:t})\), we use a TVAR model with measurement noise to simulate uncertainty about corrupted observations. The model is specified by the following equation set
\[ \begin{aligned} \boldsymbol{\theta}_t & \sim \mathcal{N}\left(\boldsymbol{\theta}_{t-1}, \mathbf{I}_M\right) \qquad\qquad\qquad\text{(38a)}\\ \mathbf{s}_t & \sim \mathcal{N}\left(\mathbf{B}\left(\boldsymbol{\theta}_t\right) \mathbf{s}_{t-1}+\mathbf{v}_u \eta, \mathbf{V}(\gamma)\right) \qquad\text{(38b)}\\ x_t & \sim \mathcal{N}\left(\mathbf{v}_u^{\top} \mathbf{s}_t, \tau\right)\qquad\qquad\qquad\text{(38c)} \end{aligned} \]
with priors
\[ \begin{gathered} \boldsymbol{\theta}_0 \sim \mathcal{N}(\mathbf{0}, \mathbf{I}) \quad \mathbf{s}_0 \sim \mathcal{N}(\mathbf{0}, \mathbf{I}) \quad \eta \sim \mathcal{N}(0.0,10.0) \qquad \text{(39a)}\\ \gamma \sim \Gamma(1.0,1.0) \quad \tau \sim \Gamma(0.1,1.0) \qquad \text{(39b)} \end{gathered} \]
Since the temperature data is not centered around \(0 \degree\text{C}\), we added a bias term \(\eta\) to the state \(\mathbf{s}_t\). The corresponding FFG is depicted in Figure 7.
Note that we put a Gamma prior on the measurement noise precision \(\tau\), meaning that we are uncertain about the size of the error of the thermometer reading. The inference task for the model is computing \(q(\mathbf{s}_{t}|x_{1:t})\), in other words, we track the states based only on past data. Of course, after training, we could use the model for temperature prediction by tracking \(q(\mathbf{s}_{t+k}|x_{1:t})\) for \(k\geq 1\). We compare the performance of four TVAR models with AR orders \(M = \{1, 2, 3, 4\}\). To choose the best model, we computed the average FE score for each TVAR(M) model.
Figure 8 shows that on average TVAR(3) outperforms its competitors. The complexity vs accuracy decomposition (23) of FE explains why a lower order model may outperform higher order models. TVAR(4) maybe as accurate or more accurate than TVAR(3) but the increase in accuracy is more than offset by the increase in complexity. For the lower order models, it is the other way around: they are less complex and involve fewer computations than TVAR(3), but the loss in model complexity leads to too much loss in data modeling accuracy. Overall, TVAR(3) is the best model for this data set. Practically, we always favor the model that features the lowest FE score. In the next subsection we will use this technique (scoring FE) for online model selection.
Single-channel speech enhancement (SCSE) is a well-known challenging task that aims to enhance noisy speech signals that were recorded by a single microphone. In single microphone recordings, we cannot use any spatial information that is commonly used in beamforming applications. Much work has been done to solve the SCSE task, ranging fromWiener filter-inspired signal processing techniques [35,36] to deep learning neural networks [37]. In this paper, we use data from the speech corpus (NOIZEUS) (https://ecs.utdallas.edu/loizou/speech/noizeus/, accessed on 27 May 2021) [38] and corrupted clean speech signals with white Gaussian noise, leading to a signal-to-noise ratio (SNR)
\[ \operatorname{SNR}\left(s_{1:T}, y_{1:T}\right)=10 \log_{10}\left[\frac{\sum_t^T s_t^2}{\sum_t^T\left(s_t-y_t\right)^2}\right] \approx 13.36 \mathrm{~dB} \tag{40} \]
where \(s_{1:T} = (s_1, . . ., s_T)\) and \(y_{1:T} = (y_1, . . ., y_T)\) are clean and corrupted speech signals. \(s_t\) is a speech signal at time \(t\) and \(T\) is the length of the signal.
Historically, AR models have shown to perform well for modeling speech signals in the time (waveform) domain [39,40]. Despite the fact that speech is a highly nonstationary signal, we may assume it to be stationary within short time intervals (frames) of about 10 [ms] each [41]. Since voiced, unvoiced and silence frames have very different characteristics, we used 5 different models (a random walk model (RW), AR(1), AR(2), TVAR(1) and TVAR(2)) for each frame of 10 [ms] with 2.5 [ms] overlap. Given a sampling frequency of 8 [kHZ], each frame results into 80 samples with 20 samples overlap. The AR and TVAR models were specified by Equation (36). For each frame, we evaluated the model performance by minimized FE and selected the model with minimal FE score. We used identical prior parameters for all models where possible. To recover the speech signal we computed the mean values of \(q(\mathbf{s}_t|x_{1:T})\) of the selected model for each frame. The SNR gain of this SCSE system was
\[ \operatorname{SNR}\left(s_{1:T}, \mathbf{s}_{1:T}\right)-\operatorname{SNR}\left(s_{1: T}, y_{1:T}\right) \approx 3.7 \mathrm{~dB}. \tag{41} \]
Figure 9 show the spectrograms of the clean, noisy and filtered signal respectively.
Next, we analyze the inference results in a bit more detail. Table 2 shows the percentage of winning models for each frame based on the free energy score.
As we can see, for more than 30% of all frames, the random walk model performs best. This happens mostly because for a silent frame the AR model gets penalized by its complexity term. We recognize that in about 90% of the frames the best models are AR(1) and RW. On the other hand, for the frames where the speech signal transitions from silent or unvoiced to voiced, these fixed models start to fail and the time-varying AR models perform better. This effect can be seen in Figure 10.
Figure 11 shows the performance of the AR(2) and RW models on a frame with a voiced speech signal. For this case, the AR(2) model performs better.
Finally, Figure 12 shows how the TVAR(2) model compares to the RWmodel on one of the unvoiced/silence frames. While the estimates of TVAR(2) appear to be more accurate, it pays a bigger “price” for the model complexity term in the FE score and the RW model wins the frame.
We have introduced a TVAR model that includes efficient joint variational Bayesian tracking of states, parameters and free energy. The system can be used as a plug-in module in factor graph-based representations of other models. At several points in this paper, we have made some design decisions that we shortly review here.
While FE computation for the AR node provides a convenient performance criterion for model selection, we noticed in the speech enhancement simulation that separate FE tracking for each candidate model leads to a large computational overhead. There are ways to improve the model selection process that we used in the speech enhancement simulation. One way is to consider a mixture model of candidate models and track the posterior over the mixture coefficients [42]. Alternatively, a very cheap method for online Bayesian model selection may be the recently developed Bayesian Model Reduction (BMR) method [43]. The BMR method is based on a generalization of the Savage-Dickey Density Ratio and supports tracking of free energy of multiple nested models with almost no computational overhead. Both methods seem to integrate well with a factor graph representation and we plan to study this issue in future work.
In this paper, the posterior factorization (31) supports the modeling of temporal dependencies between input and output of the AR node in the posterior. Technically, (31) corresponds to a structural VMP assumption, in contrast to the more constrained meanfield VMP algorithm that would be based on \(q(\mathbf{z}) = \prod_i q_i(z_i)\), where \(\mathbf{z}\) is the set of all latent variables [44]. We could have also worked out alternative update rules for the assumption of a joint factorization of precision \(\gamma\) and AR coefficients \(\boldsymbol{\theta}\). In that case, the prior (incoming message \(\overrightarrow{\nu}(\boldsymbol{\theta},\gamma)\) to AR node) would be in the form of a Normal-Gamma distribution. While any of these these assumptions are technically valid, each choice accepts a different trade-off in the accuracy vs. complexity space. We review structural VMP in Appendix A.1.
In the temperature modelling task, we added some additional random variables (bias, measurement noise precision). To avoid identifiability issues, in (38a) we fixed the covariance matrix of the time-varying AR coefficient to the identity matrix. In principle, this constraint can be relaxed. For example, an Inverse-Wishart prior distribution can be added to the covariance matrix.
In our speech enhancement experiments in Section 6.3, we assume that the measurement noise variance is known. In a real-world scenario, this information is usually not accessible. However, online tracking of measurement noise or other (hyper-)parameters is usually not a difficult extension when the process is simulated in a factor graph toolbox such as ForneyLab [28]. If so desired, we could add a prior on the measurement noise variance and track the posterior. The online free energy criterion (23) can be used to determine if the additional computational load (complexity) of Bayesian tracking of the variance parameter has been compensated by the increase in modeling accuracy.
The realization of the TVAR model in ForneyLab comes with some limitations. For large smoothing problems (say, >1000 data points), the computational load of message passing in ForneyLab becomes too heavy for a standard laptop (as was used in the paper). Consequently, in the current implementation it is difficult to employ the AR node for processing large time series on a standard laptop. To circumvent this issue, when using ForneyLab, one can combine filtering and smoothing solutions into a batch learning procedure. In future work we plan to remedy this issue by some ForneyLab refactoring work. Additionally, the implemented AR node does not provide a closed-form update rule for the marginal distribution when the probability distribution types of the incoming messages (priors) are different from the ones used in our work. Fortunately, ForneyLab supports resorting to (slower) sampling-based update rules when closed-form update rules are not available.
We presented a variational message passing approach to tracking states and parameters in latent TVAR models. The required update rules have been summarized and implemented in the factor graph package ForneyLab.jl, thus making transparent usage of TVAR factors available in freely definable stochastic dynamical systems. Aside from VMP update rules, we derived a closed-form expression for the variational free energy (FE) of an AR factor. Free Energy can be used as a proxy for Bayesian model evidence and as such allows for model performance comparisons between the TVAR models and alternative structures. Owing to the locality and modularity of the FFG framework, we demonstrated how AR nodes can be applied as plug-in modules in various dynamic models. We verified the correctness of the rules on a synthetic data set and applied the proposed TVAR model to a few relatively simple but different real-world problems. In future work, we plan to extend the current factor graph-based framework to efficient and transparent tracking of AR model order and to online model comparison and selection with alternative models.