Variational Temporal Deep Generative Model for Radar HRRP Target Recognition

For a high resolution radar (HRR), the wavelength of the radar signal is far smaller than the target size. Intuitively, HRRP is a representation of the time domain response of the target to an HRR pulse as a one-dimensional signature, which is the expression of the distribution of radar scattering centers along with the radar LOS [16][20][22], as shown in Fig. 1. Suppose the transmitted signal of an HRR is $s(t)e^{j2\pi f_{c}t}$, where $s(t)$ denotes the complex envelop and $f_{c}$ is the radar carrier frequency. Following the literatures [13, 19], the $n$th complex echo from the $d$th range cell ($n=1,2,...,N,d=1,2,...,D$) in the baseband is defined as follows

where $P_{d}$ represents the number of target scatterers in the $d$th range cell, $\lambda$ is the wavelength of HRR, $\sigma_{di}$ and $\theta_{di0}$ represent the intensity and initial phase of the $i$th scatterer in the $d$th range cell, respectively. $R_{di}(n)$ can be regarded as the radial distance between the $i$th scatterer in the $d$th range cell of the $n$th returned echo and the radar.

We set $R(n)$ as the radial distance between the target reference center in the $n$th returned echo and the radar, and set $L_{x}$ as the radial length of the target, usually $R(n)\gg L_{x}$. Due to the target rotation, there are radial displacements for the scatterers. Then $R_{di}(n)\approx R(n)+\Delta r_{di}(n)$, where $\Delta r_{di}(n)$ represents the radial displacement of the $i$th scatterer in the $d$th range cell in the $n$th returned echo. When $s(\cdot)$ is a rectangular pulse signal with unit intensity, it is usually omitted. Equation (1) can be approximated as $x_{d}(t,n)\approx x_{d}(n)=\sum_{i}^{P_{d}}\sigma_{di}e^{-j\frac{4\pi}{\lambda}R(n)}e^{j\phi_{di}(n)}$, where $\phi_{di}(n)=\theta_{di0}-\frac{4\pi}{\lambda}\Delta r_{di}(n)$ denotes the remained echo phase of the $i$th scatterer in the $d$th range cell of the $n$th returned echo, and $e^{-j\frac{4\pi}{\lambda}R(n)}$ represents the initial phase of the $n$th returned echo related to the target distance and radar wavelength. Since $e^{-j\frac{4\pi}{\lambda}R(n)}$ does not contain the target discriminative information, we can eliminate it and define the $n$th real HRRP sample $\boldsymbol{x}(n)$ as

where $|\cdot|$ means taking absolute value, and $D$ is the dimensionality of $\boldsymbol{x}(n)$.

Denoting a sample as $T$ sequentially observed $V$-dimensional vectors, we present a short review of the existing dynamical models used to model the temporal dependence across the range cells in HRRP [16, 17, 25, 26, 31].

Hidden Markov Model: The joint likelihood of the observation and underlying state sequence can be expressed as

where $\boldsymbol{s}_{n}=\{\boldsymbol{s}_{1,n},..,\boldsymbol{s}_{t,n},..,\boldsymbol{s}_{T,n}\}$ denotes the latent states of the $n$th sample, $\boldsymbol{\omega}_{0}$ is the initial state probability, $\boldsymbol{\Pi}$ is the state transition distribution and $\boldsymbol{\phi}$ are a set of parameters governing the emission probability.

Linear Dynamical System: It consists of the following observation and state equations

where $\boldsymbol{\Pi}$ is the transition matrix, $\Sigma$ and $\Delta$ are covariance matrices, and $\boldsymbol{s}_{t,n}$ denotes the latent state that is projected to the observed space via the factor loading matrix $\boldsymbol{\Phi}$.

Recurrent Neural Network: At timestep $t$, RNN reads the symbol $\boldsymbol{x}_{t,n}$ and updates its hidden state $\boldsymbol{h}_{t,n}^{(l)}$ at layer $l$

where $f$ is a deterministic non-linear function, ${\bf W}_{hh}^{(l)}$ is the transition matrix, ${\bf W}_{xh}^{(l)}$ is the weight matrix at layer $l$, and the bias vectors are omitted for conciseness.

Although HMM and LDS have been widely studied, their representation power may be limited when modeling the HRRP sequential samples. Compared with LDS and HMM, RNN typically owns extra expressive power due to the existence of deep hidden states and flexible non-linear transition function. However, the internal structure of RNN is in general entirely deterministic, with the only source of variability provided by its conditional output probability model, which may be inappropriate to model the kind of variability observed in the HRRP data [33].

According to Xu et al. [31], to discover the temporal dependence between the range cells within the single HRRP, we divide the $n$th $D$-dimensional real HRRP $\boldsymbol{x}(n)$ in (2) into $\boldsymbol{x}_{n}\in\mathbb{R}_{+}^{V\times T}$, which consists of $T$ sequentially observed $V$-dimensional vector. Shown in Fig. 2, we denote the HRRP sequence as $\boldsymbol{x}_{n}=[\boldsymbol{x}_{1,n},..,\boldsymbol{x}_{t,n},..\boldsymbol{x}_{T,n}]$, where $\boldsymbol{x}_{t,n}\in\mathbb{R}_{+}^{V\times 1}$ is the time sequential feature at timestep $t$ of sample $n$ and can be defined as

where $V$ is the size of window function, intercepting $V$ range cells from the real HRRP sample, $a_{t}=(V\!-\!O)\ast(t\!-\!1)$ is the index at timestep $t$, and $O$ denotes the overlap length across the windows, determining the degree of correlation between adjacent timesteps. Thus $T=\left\lfloor{(D-O)/(V-O)}\right\rfloor$. We denote ${\bf X}=\{\boldsymbol{x}_{n}\}_{n=1}^{N}$ as the HRRP dataset, which consists of $N$ independent and identically distributed (IID) HRRP sequences, and ${\bf Y}=\{y_{n}\}_{n=1}^{N}$ as corresponding labels, where $y_{n}$ is the label of $\boldsymbol{x}_{n}$. Note sequential inputs $\boldsymbol{x}_{1:T,n}$ from one HRRP sample $\boldsymbol{x}_{n}$ share the same label. The dataset $\{{\bf X},{\bf Y}\}$ can be fed into a dynamic model as $N$ IID samples, each of which contains $T$ sequentially observed $V$-dimensional vectors (with identical label). For simplicity, we only exhibit the modeling process for the $n$th HRRP sequence.

To characterize the sequential feature within a single HRRP sample, we generalize the deep Poisson gamma dynamical system (DPGDS) [40] to rGBN, whose generative model is sketched in Fig. 3 (a). Specifically, we consider the deep architecture with $L$ layers, and denote $\boldsymbol{s}_{t,n}^{(l)}\in\mathbb{R}_{+}^{K_{l}}$ as the latent state of $\boldsymbol{x}_{t,n}$ in (6) at layer $l$, time step $t$, where $K_{l}$ is the number of states at layer $l$. Different from HMM and LDS whose single-layer latent state at time step $t$ only relies on the state at previous time step $t-1$, the latent state $\boldsymbol{s}_{t,n}^{(l)}$ of rGBN, from the top layer to bottom, is formulated as

where $x\sim\mbox{Gam}(a,1/b)$ denotes the gamma distribution with shape $a$, scale $1/b$, and mean $a/b$; $\boldsymbol{\Pi}^{(l)}\in\mathbb{R}_{+}^{K_{l}\times K_{l}}$ the transition matrix of layer $l$, $\boldsymbol{\Phi}^{(l)}\in\mathbb{R}_{+}^{K_{l-1}\times K_{l}}$ the weight matrix connecting layers $l-1$ and $l$, $K_{0}=V$, and $1/b^{(l)}$ the gamma scale parameter at layer $l$. When $t=1$, $\boldsymbol{s}_{1,n}^{(l)}\sim\mbox{Gam}\left(b^{(l)}\boldsymbol{\Phi}^{(l+1)}\boldsymbol{s}_{1,n}^{(l+1)},1/b^{(l)}\right)$ for $1\leq l<L$ and $\boldsymbol{s}_{1,n}^{(L)}\sim\mbox{Gam}\left(b^{(L)}\textrm{1}_{K_{L}},1/b^{(L)}\right)$. When $t>1$, $\boldsymbol{s}_{t,n}^{(L)}\sim\mbox{Gam}\left(b^{(L)}\boldsymbol{\Pi}^{(L)}\boldsymbol{s}_{t-1,n}^{(L)},1/b^{(L)}\right)$.

The gamma shape parameter of $\boldsymbol{s}_{t,n}^{(l)}$ can be divided into two parts. One is the product of the transition matrix $\boldsymbol{\Pi}^{(l)}$ and latent state $\boldsymbol{s}_{t-1,n}^{(l)}$, capturing the temporal dependence at the current layer, while the other is the product of the connection weight matrix $\boldsymbol{\Phi}^{(l+1)}$ and latent state $\boldsymbol{s}_{t,n}^{(l+1)}$, capturing the hierarchical dependence at the current time. Moving beyond RNN using deterministic non-linear transition functions, we construct a dynamic probabilistic model using gamma distributed non-negative hidden units. Therefore, the proposed model is characterized by its expressive structure, which not only captures the correlations between the latent states across all layers and time steps, but also models the variability of the latent states, improving its ability to model sequential HRRP.

As a generative model with the deep temporal structure, the observed HRRP sequence at time step $t$ can be drawn from $p(\boldsymbol{x}_{t,n}\,|\,\boldsymbol{\Phi}^{(1)},\boldsymbol{s}_{t,n}^{(1)})$. To consider the non-negative nature of the HRRP data and facilitate inference, we introduce the Poisson randomized gamma (PRG) distribution defined as $\mbox{PRG}\left(\boldsymbol{x}_{t,n}\,|\,\boldsymbol{\Phi}^{(1)}\boldsymbol{s}_{t,n}^{(1)},c\right)$ in Zhou et al. [46], which has a point mass at $\boldsymbol{x}_{t,n}=0$ and is continuous for $\boldsymbol{x}_{t,n}>0$. Since the PRG distribution is generated as a Poisson mixed gamma distribution, the data likelihood can be expressed as

where $c>0$, $x\sim\mbox{Pois}(\lambda)$ represents the Poisson distribution with mean $\lambda$ and variance $\lambda$, and $\boldsymbol{\Phi}^{(1)}\in\mathbb{R}_{+}^{V\times K_{1}}$ is the weight matrix of layer $1$.

For scale identifiability and ease of inference and interpretation, we place the Dirichlet distribution prior on each column of $\boldsymbol{\Phi}^{(l)}$ and $\boldsymbol{\Pi}^{(l)}$ , $i.e.$, $\boldsymbol{\phi}_{k}^{(l)}$ and $\boldsymbol{\pi}_{k}^{(l)}$, by letting

for $l\in{1,...,L}$, which makes the elements of each column be non-negative and sum to one. Note $\boldsymbol{\pi}_{k}^{(l)}=(\boldsymbol{\pi}_{1,k}^{(l)},..,\boldsymbol{\pi}_{k_{1},k}^{(l)},..,\boldsymbol{\pi}_{K_{l},k}^{(l)})$ and $\boldsymbol{\pi}_{k_{1},k}^{(l)}$ describes how much the weight of state $k$ of the previous time at layer $l$ is transited to influence state $k_{1}$ of the current time at the same layer. Under the hierarchical dynamical model defined by (7) and (8), the joint likelihood of the observation HRRP and the temporal latent states can be constructed as

The parameters of the generative model comprise the transition and weight matrices, which we write as $\{\boldsymbol{\Pi}^{(l)},\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$. Fig. 3 (a) shows the graphical representation of the proposed generative model and Fig. 3(b) is the unfolded representation of the model structure.

Structure Analysis: As discussed above, if $x\sim\mbox{Gam}(a,1/b)$, the mean of $x$ is $a/b$; while if $x\sim\mbox{Pois}(\lambda)$, the mean of $x$ is $\lambda$. Therefore, the expected value of $\boldsymbol{x}_{t,n}$ in (8) and $\boldsymbol{s}_{t,n}^{(l)}$ in (7) can be expressed as

Based on (12) and (13), for a three-hidden layer rGBN shown in Fig. 3(a), we have

where the expected value of $\boldsymbol{x}_{t,n}$ depends on the latent states at the previous time step in each layer, indicating our proposed model captures and transmits long-range temporal information through its higher hidden layers. In addition, the proposed model can be viewed as the generalization of LDS and HMM with deep gamma distributed latent representations, and also can be considered as a probabilistic construction of the traditionally deterministic RNN by adding uncertainty into the latent space via a deep generative model.

Ignoring the temporal structure of the equation, we notice that the information of the whole HRRP data set can be compressed into the inferred network $\{\boldsymbol{\Phi}^{(1)},\boldsymbol{\Phi}^{(2)},\boldsymbol{\Phi}^{(3)}\}$, which depicts the global structure of the target in a single HRRP. To be more specific, we can visualize the $\boldsymbol{\phi}_{k}^{(l)}$ at layer $l$ as $\left[\prod_{p=1}^{l-1}\boldsymbol{\Phi}^{(p)}\right]\boldsymbol{\phi}_{k}^{(l)}$, which are often quite specific at the bottom layer and become increasingly more general when moving upwards. We will examine the weights of different layers to understand the general and specific aspects of the HRRP data, and illustrate how the weights of different layers are related to each other.

Quantizing the HRRP data: While the non-negative real HRRP dataset can be modeled by PRG distribution, the latent count $\boldsymbol{r}_{t,n}$ in (8) need to be sampled in each iteration of the training stage. Specifically, the conditional posterior of $\boldsymbol{r}_{t,n}$ given $\boldsymbol{x}_{t,n}$, $\boldsymbol{\Phi}^{(1)}\boldsymbol{s}_{t,n}^{(1)}$, and $c$ can be expressed as

More details about the PRG distribution can be found in Zhou et al. [46] and are omitted here for brevity. Despite the desired ability to depict non-negative continuous data, the PRG distribution may be time-consuming for the iterative procedure to infer latent counts in real-world applications. Considering the limited computation power in the training stage, HRRP data can be modeled with Poisson distribution, by directly discretizing the HRRP data to counts ($i.e.$, non-negative integers) before training, where the sampling step of the counts in each iteration can be avoided. Therefore, instead of modeling the observed data with (8), we assume that

where the scaling factor $\mu$ controls the fineness of this discretization. There is a trade-off between the benefit of the PRG distribution and restriction of computation resources. Note that modeling the observed HRRP with the Poisson distribution can still obtain similar temporal dependencies in (III-B) and hierarchical structure $\prod_{p=1}^{l}\boldsymbol{\Phi}^{(p)}$, omitted here for brevity.

In this section, we first infer the generative model global parameters $\{\boldsymbol{\Pi}^{(l)}\}_{l=1}^{L}$ and $\{\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$ with MCMC, then introduce a recurrent inference model to infer the latent states $\{\boldsymbol{s}_{t,n}^{(l)}\}_{t=1,l=1,n=1}^{T,L,N}$. Finally, we provide a hybrid SG-MCMC and recurrent variational inference, which is scalable at the training stage and fast at the testing stage.

MCMC inference for the generative network: Given the HRRP sequences, the inference task here is to find the weight matrices $\{\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$, transition matrices $\{\boldsymbol{\Pi}^{(l)}\}_{l=1}^{L}$, and latent states $\{\boldsymbol{s}_{t,n}^{(l)}\}_{t=1,l=1,n=1}^{T,L,N}$. While it is difficult to infer the introduced model for the coupling of $\{\boldsymbol{s}_{t,n}^{(l)}\}_{t=1,l=1,n=1}^{T,L,N}$ with $\{\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$ and $\{\boldsymbol{\Pi}^{(l)}\}_{l=1}^{L}$, the latent variables of the proposed model can be trained with a backward-upward–forward-downward (BUFD) Gibbs sampler described in Guo et al. [40], based on a variety of variable augmentation techniques. However, the Gibbs sampler needs to process all HRRP data samples in each iteration and hence has limited scalability. For scalable inference, we adopt the topic-layer-adaptive stochastic gradient Riemannian (TLASGR) MCMC algorithm [47, 43], which is proposed to update simplex-constrained global parameters [48] in a mini-batch learning setting. Relying on the Fisher information matrix (FIM) [49] to automatically adjust relative learning rates for different parameters across all layers, TLASGR-MCMC has proven its improved sampling efficiency. Given the latent states, we first sample augmented latent counts, then use TLASGR-MCMC to sample generative model parameters $\{\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$ and $\{\boldsymbol{\Pi}^{(l)}\}_{l=1}^{L}$. To be more specific, we can sample $\boldsymbol{\pi}_{k}^{(l)}$, the $k$th column of the transition matrix $\boldsymbol{\Pi}^{(l)}$, as

where ${\left[.\right]_{\angle}}$ denotes the simplex constraint that $\boldsymbol{\pi}_{k_{1},k}^{(l)}\geqslant 0$ and $\sum_{k_{1}=1}^{K_{l}}\boldsymbol{\pi}_{k_{1},k}^{(l)}=1$, $M_{k}^{(l)}$ is calculated using the estimated FIM, $\varepsilon_{i}$ denotes the learning rate at the $i$th iteration, both ${\tilde{}\boldsymbol{z}_{:k\boldsymbol{\cdot}}^{\left(l\right)}}$ and ${\tilde{z}_{\boldsymbol{\cdot}k\boldsymbol{\cdot}}^{\left(l\right)}}$ come from the augmented latent counts ${\bf Z}^{(l)}$, and ${\boldsymbol{\nu}^{\left(l\right)}}$ denotes the prior of ${\boldsymbol{\pi}_{k}^{\left(l\right)}}$. More details of TLASGR-MCMC can be found in Cong et al. [47] and Guo et al. [40].

Despite the attractive properties, both the proposed Gibbs sampler and TLASGR-MCMC usually rely on an iterative procedure to learn the temporal latent states of a new HRRP sample at the testing stage, which hinders real-time processing of the HRRP based RATR. To allow fast out-of-sample prediction, we further build an inference network to learn the latent states, as described below.

Recurrent variational inference model: Our inference model is motivated by variational auto-encoders (VAEs) [42]. As an example of a directed graphical model, the joint distribution over the observed variables $\boldsymbol{x}$ and latent variables $\boldsymbol{z}$ can be defined as $p(\boldsymbol{x},\boldsymbol{z})=p(\boldsymbol{x}\,|\,\boldsymbol{z})p(\boldsymbol{z})$, where $p(\boldsymbol{z})$ is the prior placed on the latent variables. To admit efficient inference, VAEs approximate $p(\boldsymbol{z}\,|\,\boldsymbol{x})$ with a variational family of distributions $q(\boldsymbol{z}\,|\,\boldsymbol{x})$, which adopts an inference network to map the observations directly to their latent space by optimizing the evidence lower bound (ELBO) [33][50][51], expressed as

The approximate posterior $q(\boldsymbol{z}\,|\,\boldsymbol{x})$ is often taken to be a Gaussian distribution as $N(\boldsymbol{u},diag(\boldsymbol{\sigma}))$, where the mean $\boldsymbol{u}$ and standard deviation $\boldsymbol{\sigma}$ are the output of highly non-linear functions of the input $\boldsymbol{x}$. Given the Gaussian variational posterior, the second term of the ELBO in (18) is analytic. Using the reparameterization trick [42], VAEs sample $\boldsymbol{z}$ with $\boldsymbol{z}=\boldsymbol{u}+\boldsymbol{\sigma}\odot\boldsymbol{\epsilon}$, where $\boldsymbol{\epsilon}$ is a vector of standard Gaussian variables. Therefore, the gradient of the first term of the ELBO with respect to the parameters of the inference network can be constructed as

whose estimation via Monte Carlo integration has low variance. The parameters of the inference model can be typically optimized via SGD.

VAEs provide an effective modeling paradigm for complex data distributions. However, their success so far is mostly restricted to non-sequential data with Gaussian distributed latent variables and does not generalize well to model non-negative and sequential HRRP representations. In this section, we propose a recurrent variational inference method to efficiently produce the multilayer temporal HRRP representations with (7) and (8) or (16) as the generative model. Given the global parameters $\{\boldsymbol{\Pi}^{(l)},\boldsymbol{\Phi}^{(l)}\}_{l=1}^{L}$, the task here is to infer the latent states $\{\boldsymbol{s}_{t,n}^{(l)}\}_{t=1,l=1,n=1}^{T,L,N}$ via an inference network. We first introduce a fully factorized distribution as

With (III-B) and (20), the objective function becomes

where $\boldsymbol{a}_{t,n}^{\left(l\right)}:=b^{(l)}\left({{\boldsymbol{\Phi}^{(l+1)}}\boldsymbol{s}_{t,n}^{(l+1)}+{\boldsymbol{\Pi}^{(l)}}\boldsymbol{s}_{t{\rm{-}}1,n}^{(l)}}\right)$ denotes the shape parameter of $\boldsymbol{s}_{t,n}^{(l)}$. Note under the BUFD Gibbs sampler [40], the conditional posterior of $\boldsymbol{s}_{t,n}^{(l)}$ given augmented latent counts is gamma distributed, and hence it might be more appropriate to use the gamma rather than Gaussian based distributions to construct the variational distribution $q\left({\boldsymbol{s}_{t,n}^{(l)}}\right)$. However, the gamma random variable is not reparameterizable with respect to its shape parameter. This motivates us to choose a surrogate distribution, which can be easily reparameterized, to approximate the gamma distribution. The Weibull distribution is a desirable choice for this purpose, as its probability density function resembles that of the gamma distribution, and the second term in the ELBO shown in (III-C) becomes analytic if it is used to construct $q\left({\boldsymbol{s}_{t,n}^{(l)}}\right)$ [43].

We may directly follow Zhang et al. [43] to construct a Weibull distribution based inference network as

whose parameters $\boldsymbol{k}_{t,n}^{(l)}$ and $\boldsymbol{\lambda}_{t,n}^{(l)}$ are both deterministically transformed from the hidden unit $\boldsymbol{h}_{t,n}$ and specified as

where $\boldsymbol{h}_{t,n}^{(l)}$ denotes the output of highly non-linear function of the observed $\boldsymbol{x}_{t,n}$. However, this construction does not take into consideration the temporal information transmitted from the previous time step. To exploit the temporal information, we propose a recurrent inference network which induces temporal dependencies across time steps, as illustrated in Fig. 3 (c). Therefore, similar to the RNN in (5), we define $\boldsymbol{h}_{t,n}^{(l)}$ as

where at layer $l$, ${\bf W}_{xh}^{(l)}\in\mathbb{R}^{K_{l}\times K_{l-1}}$ denotes the upward weight matrix, ${\bf W}_{hh}^{(l)}\in\mathbb{R}^{K_{l}\times K_{l}}$ the forward weight matrix connecting the hidden states, and $\boldsymbol{b}_{3}^{(l)}\in\mathbb{R}^{K_{l}}$ the bias vector. Therefore, the variational parameters of $\boldsymbol{s}_{t,n}^{(l)}$ are both non-linearly transformed with neural networks from the hidden state $\boldsymbol{h}_{t,n}^{(l-1)}$ of layer $l-1$ at time $t$ and the hidden state $\boldsymbol{h}_{t-1,n}^{(l)}$ of layer $l$ at time $t-1$, which are helpful for the inference network to take into account the temporal structure of the observed data.

While the sequential latent states of the inference model are similar to that of an RNN, the internal transition structure of an RNN is in general entirely deterministic. By contrast, the proposed model introduces uncertainty into the latent space to help better model the variability observed in highly structured HRRP data. One can sample the Weibull distributed latent state in (22) using the reparameterization trick as

For standard VAE, the generative model parameters and the corresponding inference network parameters can be typically jointly optimized via SGD, seeking to maximize the ELBO in (18) with standard backpropagation technique. Instead of finding a point estimate of the global parameters of generative model like in VAEs, we adopt a hybrid MCMC/VAE inference algorithm by combining TLASGR-MCMC and the proposed recurrent variational inference network. In specific, the generative model parameters $\{\boldsymbol{\Pi}^{(l)},\boldsymbol{\Phi}^{(l)}\}_{1,L}$ can be sampled with TLASGR-MCMC in (III-C) and the neural network parameters $\boldsymbol{\Omega}=\{{\bf W}_{xh}^{(l)},{\bf W}_{hh}^{(l)},{\bf W}_{hk}^{(l)},{\bf W}_{h\lambda}^{(l)},\boldsymbol{b}_{1}^{(l)},\boldsymbol{b}_{2}^{(l)},\boldsymbol{b}_{3}^{(l)}\}_{1,L}$ can be updated via SGD by maximizing the ELBO in (III-C). Applying the reparameterization trick of the Weibull distribution in (26), the gradient of the ELBO with respect to $\boldsymbol{\Omega}$ can be evaluated with low variance. In practice, a single Monte Carlo sample from $q\left(\boldsymbol{s}_{t,n}^{(l)}\right)$ is enough to obtain satisfactory performance. The proposed inference network is depicted in Fig. 3 (c) and the training strategy is outlined in Algorithm 1. For a new HRRP sample at the testing stage, given the generative model parameters and inferred network parameters $\boldsymbol{\Omega}$, we can directly obtain the conditional posteriors of latent states using the inference network, without performing iterations.

To illustrate the proposed model and compare it to existing methods, here we consider several three-dimensional synthetic datasets. Each toy data discussed below can be denoted as $\boldsymbol{x}\in\mathbb{R}^{V\times T}$ with $V=3,T=100$.

$\bullet$ Toy data 1: To verify whether our proposed model can learn the transition matrix accurately, we generate the data with rGBN-PRG, where we assume the transition matrix as $\boldsymbol{\Pi}=\begin{bmatrix}0.65&0.20\\ 0.35&0.80\\ \end{bmatrix}$, initial latent state as $\boldsymbol{s}_{1}\!=\!\begin{bmatrix}100&0\\ \end{bmatrix}^{\mathrm{T}}$, and the Dirichlet distributed weight matrix as $\boldsymbol{\Phi}\in\mathbb{R}^{3\times 2}$. The latent states can be generated with $\boldsymbol{s}_{t}\sim\mbox{Gam}\left(\boldsymbol{\Pi}\boldsymbol{s}_{t-1},1\right)$, where $\boldsymbol{s}_{t}\in\mathbb{R}^{2\times 1},T=100$. The non-negative observed data can be sampled from $\boldsymbol{r}_{t}\sim\mbox{Pois}\left(\boldsymbol{\Phi}\boldsymbol{s}_{t}\right),\boldsymbol{x}_{t}\sim\mbox{Gam}\left(\boldsymbol{r}_{t},1\right)$, where $\boldsymbol{x}_{t}\in\mathbb{R}^{3\times 1}$.

$\bullet$ Toy data 2: $x_{1,t}=t$, $x_{2,t}=2\exp(-t/15)+\exp(-((t-25)/10)^{2})$, and $x_{3,t}=5\sin(t^{2})+6$ for $t=1,\ldots,100$.

$\bullet$ Toy data 3: $x_{1,t}=t$, $x_{2,t}=2\mbox{mod}(t,3)$, $x_{3,t}=20\exp(-t/3)$ for $t=1,\ldots,50$, and $x_{1,t}=2t+30$, $x_{2,t}=3\mbox{mod}(t,2)+5$, and $x_{3,t}=30t\exp(-t)+10$ for $t=51,\ldots,100$, where $\mbox{mod}(t,k)$ denotes the modulo operation that returns the remainder after division of $t$ by $k$.