Online Expectation-Maximization Based Frequency and Phase Consensus in Distributed Phased Arrays

We consider a directed array network of $N$ nodes described by a directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ as before, whereas the $(n,m)$-th edge in $\mathcal{E}$ is assigned herein a weight $w_{m,n}$. We assume that $\mathcal{G}$ represents a strongly connected directed graph, i.e., there exists at least one directed path from any node $n$ to any node $m$ in the graph with $n\neq m$. Furthermore, it is assumed that the $n$-th node in the network knows its in-neighbor set $\chi^{in}_{n}\triangleq\{v\in\mathcal{V}:w_{n,v}>0\}$ and its out-neighbor set $\chi^{out}_{n}\triangleq\{v\in\mathcal{V}:w_{v,n}>0\}$. The proposed algorithm is based on the push-sum consensus algorithm [24, 25], and is referred to herein as the push-sum frequency and phase consensus ($\text{P}_{\text{s}}$FPC) algorithm. In the $\text{P}_{\text{s}}$FPC algorithm, the $n$-th node maintains three temporary variables in the $k$-th iteration, i.e., $x^{f}_{n}(k)$, $x^{\theta}_{n}(k)$, and $s_{n}(k)$, which are updated as follows.

where $s_{n}(0)=1$ and the weight $w_{n,m}\triangleq\frac{1}{d^{out}_{m}}$ if $m\in\chi^{in}_{n}$ and is $0$ otherwise. The parameter $d^{out}_{m}$ is the out-degree of node $m$ which can be found by computing the cardinality of the set $\chi^{out}_{m}$. These weights make the $N\times N$ weighting matrix $\mathbf{W}=[w_{n,m}]$ a column stochastic matrix that supports an average consensus [32]. The parameters $\hat{f}_{m}(k)$ and $\hat{\theta}_{m}(k)$ represent the estimated frequency and phase of the $m$-th node in the $k$-th iteration, and they are multiplied by the weighting factor $s_{m}(k-1)$ from the previous iteration. This multiplication is required to recover the temporary variables $x^{f}_{n}(k-1)$ and $x^{\theta}_{n}(k-1)$ taking into account the fact that at the beginning of each iteration, the nodes observe only the estimated frequencies and phases after their updated values in the previous iteration are influenced by the oscillator drifts and the phase jitters. The updated frequency and phase of the $n$-th node in the $k$-th iteration are computed by

The above steps are repeated iteratively at each node until the algorithm converges. The pseudo-code of the $\text{P}_{\text{s}}$FPC algorithm is provided in Algorithm 1.

Herein, we theoretically analyze the residual phase errors of the nodes upon the convergence of the $\text{P}_{\text{s}}$FPC algorithm. To this end, let at the $I$-th iteration, the estimated frequencies and phases of the nodes be written in the short-hand vector form as $\hat{\mathbf{z}}(I)=\mathbf{z}(I-1)+\mathbf{e}_{I}$ in which $\mathbf{z}\in\{\mathbf{f},\bm{\theta}\}$ and $\mathbf{e}_{I}\sim\mathcal{N}\left(\mathbf{0},\sigma^{2}_{e}\mathbf{I}_{N}\right)$. Note that by using (II) and (II), it can be seen that when $\mathbf{z}=\mathbf{f}$ then the error variance is $\sigma^{2}_{e}=\sigma^{2}_{f}+\left(\sigma^{m}_{f}\right)^{2}$, and when $\mathbf{z}=\bm{\theta}$ then the variance is $\sigma^{2}_{e}=\pi^{2}T^{2}\sigma^{2}_{f}+\left(\sigma^{m}_{\theta}\right)^{2}+\sigma^{2}_{\theta}$. Thus, at the $I$-th iteration, the temporary variables in (II-A) can be alternatively written as

where $\odot$ is the Hadamard product between the two vectors, and the operation $[.]_{n}$ selects the $n$-th element of the resulting vector. Inserting $\hat{\mathbf{z}}(I)$ in (II-B) and using backward recursion, it can be reduced as follows

Similarly we have $s_{n}(I)=\left[\mathbf{W}^{I}\mathbf{s}(0)\right]_{n}$. For a column stochastic matrix $\mathbf{W}$, it is shown in [33] that as $I\rightarrow\infty$, $\mathbf{W}^{I}=\bm{\pi}\mathbf{1}^{T}$ where $\bm{\pi}>0$ is the stationary distribution of the random process with transition matrix $\mathbf{W}$ that satisfies $\mathbf{W}\bm{\pi}=\bm{\pi}$ with $\mathbf{1}^{T}\bm{\pi}=1$. Thus from (II-A), the updated frequency or phase of node $n$ at the $I$-th iteration can be written in the shorthand form as

where to get (II-B) we set $\mathbf{s}(0)=\mathbf{1}$ (as in Section II-A) that gives for a large $I$, $s_{n}(I)=\mathbf{W}^{I}\mathbf{s}(0)=\pi_{n}N$ in which $\pi_{n}$ is the $n$-th element of $\bm{\pi}$. Similarly, for larger $I$, the first summand in (II-B) reduces to $\frac{1}{N}\sum^{N}_{n=1}z_{n}(0)$, i.e., it gives the average of the initial frequencies and phases of the nodes, whereas the second summand in (II-B) represents the residual total phase error due to the frequency and phase offset errors induced at the nodes in every iteration. To quantify this residual phase error for node $n$ after $I$ iterations, we find the variance of the following term

where to get the second equality in the above equation, we used the following two facts. Firstly, a column stochastic $\mathbf{W}$ implies that $\mathbf{1}^{T}\mathbf{W}^{i}=\mathbf{1}$, and secondly, we have $\left(\mathbf{I}_{N}-\bm{\pi}\mathbf{1}^{T}\right)^{i}=\left(\mathbf{I}_{N}-\bm{\pi}\mathbf{1}^{T}\right)$ as it is a projection matrix. Now let $\lambda_{2}$ be the modulus of the second largest eigenvalue of the mixing matrix $\mathbf{W}$. Then it can be conveniently shown in a few steps that $\lambda^{i}_{2}\mathbf{I}_{N}-\left(\mathbf{W}-\bm{\pi}\mathbf{1}^{T}\right)^{i}$ is a positive semi-definite matrix. Thus, using this identity along with assuming that the frequency and phase offset errors induced at the nodes in every iteration are mutually independent across the array, the variance of the residual phase error for node $n$ can be solved as

Since for large $I$, the scalars $s_{n}(I)\rightarrow\pi_{n}N$, it implies that we can write the phase error in (9) as $\sigma^{2}_{n,\text{residual}}=\sigma^{2}_{e}\sum^{I^{{}^{\prime}}}_{i=1}\lambda^{2i}_{2}+C$ where the constant $C$ is close to zero for the larger arrays. Consequently, for sparsely connected arrays, $\lambda_{2}$ is close to $1$ and $\sum^{I}_{i=1}\lambda^{2i}_{2}\gg 1$ which results in a higher residual phase error, but as the connectivity increases, $\lambda_{2}\rightarrow 0$ and the residual phase error decreases. This variation in the total phase error as a function of the change in the array connectivity is illustrated through simulations in Section III-E.

For validation purposes, we examine the frequency and phase synchronization performances of the proposed $\text{P}_{\text{s}}$FPC algorithm through simulations. To this end, we consider an array network of $N$ nodes that is randomly generated with connectivity $c$ between the nodes. The parameter $c$ is defined as the total number of existing connection in the generated network divided by the total number of all possible connections $\frac{N}{2}(N-1)$. Thus, $c\in[0,1]$ where a smaller value of $c$ implies a sparsely connected array network with a few connections per node, and a larger value of $c$ defines a highly connected array network with many connections per node. The average number of connections per node in a network is given by $D=c(N-1)$. The carrier frequency of the nodes is chosen in this analysis to be $f_{c}=1$ GHz, and the sampling frequency is set to $f_{s}=10$ MHz throughout this paper. As the update interval for distributed phased arrays is usually on the order of ms, we set $T=0.1$ ms for all the simulated cases in this paper.

Figs. 1 and 2 show the frequencies and phases for all the $N$ nodes in the array relative to their average values vs. the number of iteration of the $\text{P}_{\text{s}}$FPC algorithm from a single trial, when a moderately connected array network with connectivity $c=0.2$ is generated for the different number of nodes $N$ in the array. The SNR of the observed signals is either set as $0$ dB or $30$ dB. These figures show that with the increase in the number of iterations, both frequencies and phases of the nodes converge to the average of their initial values for all the considered cases. The residual error upon convergence results from the oscillators’ frequency drifts and phase jitters, as well as the frequency and phase estimation errors induced at the nodes. Thus, when the SNR increases to $30$ dB for $N=100$ nodes, it is observed that the residual error also decreases significantly upon the convergence due to the decrease in the estimation errors.

Now let the total residual phase error of node $n$ be denoted by $\delta\phi_{n}=2\pi\delta f_{n}T+2\pi\varepsilon_{f}T+\delta\theta^{f}_{n}+\delta\theta_{n}+\varepsilon_{\theta}$. Thus, in Fig. 4, we plot the standard deviation of the total phase errors $\delta\phi_{n}$ for all the $N$ nodes in the array for the $\text{P}_{\text{s}}$FPC algorithm by varying the connectivity $c$ between the nodes and the number of nodes $N$ in the array, and when SNR$=30$ dB is assumed. For the comparison, we generate the bidirectional array networks for each $c$ and $N$ values and show the standard deviation of the total phase error of the DFPC algorithm proposed in [11] for such array networks. Note that the standard deviation of the DFPC-based algorithms in this work provides the lower bounds on the standard deviation of the total phase error of the $\text{P}_{\text{s}}$FPC-based algorithms, and the slight increase in the phase error of the $\text{P}_{\text{s}}$FPC-based algorithms as compared to the DFPC-based algorithms is due to the decrease in the number of connection per node in the directed networks which effects the accuracy of the local averages computed per node throughout the network as expected. For this figure, both $\text{P}_{\text{s}}$FPC and DFPC algorithms were run for $100$ iterations and the final standard deviations of the total phase errors were averaged over $10^{3}$ trials The center of each point in the error bar plot is the average value of the samples and the length of the bar defines the standard deviation around the average value. It is observed that for each $N$ value, as the connectivity $c$ between the nodes increases then the total phase error of the $\text{P}_{\text{s}}$FPC and DFPC algorithms decreases proportionately. Furthermore, an increase in the number of nodes $N$ in the array, for a given $c$ value, also decreases the total phase errors of both algorithms. As discussed in the previous subsection, the decrease in the residual phase error of $\text{P}_{\text{s}}$FPC with the increase in the $c$ or $N$ values is because the network algebraic connectivity $\lambda_{2}$ decreases which in turn reduces the residual phase of $\text{P}_{\text{s}}$FPC. The residual phase error of DFPC is also a function of $\lambda_{2}$ as shown in [11] and thus the same trend is observed for the DFPC algorithm as well.

Finally, in Fig. 4, we compare the convergence speeds of the $\text{P}_{\text{s}}$FPC and DFPC algorithms by plotting the standard deviation of the total phase error vs. the number of iteration of the algorithms when the two different $N$ and $c$ values are assumed and the SNR is set to $30$ dB. It is observed that either when $N$ increases from $20$ to $100$ for a given $c$ value, or when $c$ increases from $0.2$ to $0.5$ for a given $N$ value, the $\text{P}_{\text{s}}$FPC algorithm converges faster. Specifically, for $N=20$ and $c=0.2$, $\text{P}_{\text{s}}$FPC converges in $32$ iterations, but for this $N$ value, when $c$ increases to $0.5$, $\text{P}_{\text{s}}$FPC converges in about $11$ iterations. Moreover, for $c=0.5$ and $N=100$ nodes, $\text{P}_{\text{s}}$FPC converges in about $5$ iterations. This is because as either $c$ or $N$ increases then the average number of connections per node in the network, i.e., the $D$ value, increases which in turn quickly stabilizes the local averages computed at the nodes. Although the same trend is also observed for the DFPC algorithm in all the considered cases in this figure, particularly it is observed that the $\text{P}_{\text{s}}$FPC algorithm converges faster than the DFPC algorithm in all cases. This is because the network is bidirectional for the DFPC algorithm, and as such there are more connections per node in the network which in turn delays the convergence of the local averages computed at the nodes. Note that the $\text{P}_{\text{s}}$FPC algorithm can be deployed for a bidirectional network as well, in which case it results in the same synchronization and convergence performance as that of the DFPC algorithm; however, the DFPC algorithm does not converge in case of directed networks due to its dependence on the average consensus algorithm [15] as discussed before. In practice, since the failure of a link may not be easily detectable, this means that the $\text{P}_{\text{s}}$FPC algorithm is a preferred choice over the DFPC algorithm.

To write the state-space equation for the $n$-th node, let at the $k$-th time instant its state vector be given by $\mathbf{x}_{n}(k)=[f_{n}(k),\theta_{n}(k)]^{T}$ in which $f_{n}(k)$ and $\theta_{n}(k)$ represent its frequency and phase, respectively. Using the frequency drift and phase jitter models as described in Section II, the state-space equation for the $n$-th node is written as

in which the innovation noise vector is defined as $\mathbf{u}_{n}=[\delta f_{n},\delta\theta^{f}_{n}+\delta\theta_{n}]^{T}$, and we assume that $\mathbf{u}_{n}$ is normally distributed with zero mean and the correlation matrix $\mathbf{Q}_{n}$ which is given by

Next to write the observation equation, let the frequency and phase estimates of the signal for the $n$-th node be written in the vector form as $\mathbf{y}_{n}(k)=\left[\hat{f}_{n}(k),\hat{\theta}_{n}(k)\right]^{T}$, then in terms of the estimation errors, this observation vector is defined as

in which the measurement noise vector $\mathbf{v}_{n}$ is given by $\mathbf{v}_{n}=\left[\varepsilon_{f},\varepsilon_{\theta}\right]^{T}$ in which $\varepsilon_{f}$ and $\varepsilon_{\theta}$ denote the frequency and phase estimation errors, respectively. We assume that $\mathbf{v}_{n}$ is also normally distributed with zero mean and the correlation matrix $\bm{\Sigma}_{n}$ which is given by

in which the standard deviations $\sigma^{m}_{f}$ and $\sigma^{m}_{\theta}$ define the marginal distributions on the frequency and phase estimation errors, respectively.

In the above state-space model in (10) and (12), the measurement noise vector $\mathbf{v}_{n}$ is independent of the state vector $\mathbf{x}_{n}(k)$ and the innovation noise vector $\mathbf{u}_{n}$, thus the state vector estimate at each time instant can be easily obtained by using Kalman filtering as described below.

Kalman filtering is an online estimation algorithm that estimates the state vector $\mathbf{x}_{n}(k)$ of the $n$-th node at the current time instant $k$ by using all the observations up to the present time. This process involves sequentially computing the prediction-update step and the time-update step at each time instant. To this end, let the state vector estimate of the $n$-th node at time instant $k-1$ be given by the vector $\mathbf{m}_{n,k-1}(k-1)$ whereas its error covariance matrix is defined by the matrix $\mathbf{V}_{n,k-1}(k-1)$. As annotated in the subscripts, these parameters in KF are computed by using the observations up to time $k-1$. Thus, in the prediction-update step at time instant $k$, we use the linear transformation in (10) and predict the new state vector estimate and the corresponding error covariance matrix as follows

In the time-update step, an a priori normal distribution is assumed on the state vector $\mathbf{x}_{n}(k)$ where its mean and covariance are set equal to the predicted estimate vector and the error covariance matrix in (III-B), respectively. This a priori distribution is then used to compute the a posteriori mean and covariance of the state vector given the current observation $\mathbf{y}_{n}(k)$ by using

where the Kalman gain matrix $\mathbf{K}_{n}(k)$ is given by

The mean vector $\mathbf{m}_{n,k}(k)$ in (15) gives the MMSE estimate of the state vector at time instant $k$ and the matrix $\mathbf{V}_{n,k}(k)$ defines its error covariance matrix. These means and covariances of the state vectors are used in the EM algorithm derived below to compute the maximization step in every iteration.

The prediction update and time update steps of KF in (III-B) and (15) assume that the innovation noise covariance matrix $\mathbf{Q}_{n}$ and the measurement noise covariance matrix $\bm{\Sigma}_{n}$ are known a priori. As discussed earlier, these model parameters are unknown in practice and must be estimated for Kalman filtering. Let $\bm{\Theta}=\{\mathbf{Q}_{n},\bm{\Sigma}_{n}\}$ denote the set of these model parameters, then given all the instantaneous observations of the node $n$ up to the present time $K$ as denoted by $\mathbf{y}^{1:K}_{n}$, the marginal log-likelihood of $\bm{\Theta}$ is written as

where the conditional distributions in (III-C) can be easily found from (10) and (12) by shifting the distributions of the noise vectors to the given state vectors. Note that due to the normally distributed conditional distributions and due to the log of the summation in (III-C), the above objective function is a non-concave function of $\bm{\Theta}$, and directly maximizing it with respect to $\bm{\Theta}$ does not provide a closed-form solution for estimating $\mathbf{Q}_{n}$ and $\bm{\Sigma}_{n}$. An expectation-maximization (EM) algorithm [39, 28] is an iterative algorithm that often results in the closed-form update equations for the unknown parameters and maximizes the marginal log-likelihood function in every iteration until convergence to its local maximum or saddle point [40]. An EM algorithm starts with an initial estimate of the unknown parameters, thus if $\bm{\Theta}^{(l-1)}$ is the estimate at its $(l-1)$-st iteration, in the $l$-th iteration it computes an expectation step (E-step) and a maximization step (M-step). In the E-step, it computes the expectation of the complete data log-likelihood function using the estimate of $\bm{\Theta}$ from the previous iteration as follows.

where the expectation in (18) is a conditional expectation given the observations $\mathbf{y}^{1:K}_{n}$ and the estimate $\bm{\Theta}^{(l-1)}$. In the M-step, it maximizes this objective function with respect to $\bm{\Theta}$ to get its new estimate by solving

The above E-step and M-step are repeated iteratively until the convergence is achieved. The EM algorithm in (18) and (19) describes a batch-mode EM algorithm that is proposed in [39, 37] where first the entire data set up to the time instant $K$ is collected, and then the algorithm is run iteratively on this data set until convergence is achieved. This batch-mode EM algorithm is applicable in cases where the cost of collecting the entire data set at once is affordable, whereas in our synchronization task due to the oscillators instantaneous frequency drifts and phase jitters, we aim to instantaneously estimate $\bm{\Theta}$ and update the frequencies and phases of the nodes in every iteration to synchronize these parameters across the array.

To instantaneously estimate $\bm{\Theta}$ from the observed data, an online version of the EM algorithm is proposed in [28, 41] where the E-step in (18) is replaced by a stochastic approximation step in the $k$-th iteration as follows.

where $\gamma_{k}$ is a smoothing factor, and we define

Note that the iteration index in the E-step in (20) is the same as the time index due to the online nature of the algorithm. Assuming a constant smoothing factor with $\gamma_{k}=1-\alpha$, where $\alpha$ can be optimized on an initial dataset [41], the E-step in (20) can be written in compact form as

In the M-step of the online EM algorithm, we maximize the objective function in (III-C) with respect to $\bm{\Theta}$ as in (19).

Now we describe the computation of the auxiliary function in (III-C) as follows. By using the conditional distributions from (10), (12), and (III-C), this function can be computed as shown in (III-C), in which $m^{f}_{n,k}(k)$ and $m^{\theta}_{n,k}(k)$ are the first and second element of the vector $\mathbf{m}_{n,k}(k)$, respectively, and the scalars $\gamma^{n}_{k,1}$ and $\gamma^{n}_{k,2}$ define the $(1,1)$-th and $(2,2)$-th indexed elements of the matrix $\bm{\Gamma}^{n}_{k,k}$, respectively. The correlation matrices $\bm{\Gamma}^{n}_{k-1,k-1}$, $\bm{\Gamma}^{n}_{k,k}$, $\bm{\Gamma}^{n}_{k,k-1}$, and the mean vector $\mathbf{m}_{n,k}(k)$ in (III-C), in iteration $k$, are defined in the following. To begin, the matrix $\bm{\Gamma}^{n}_{k-1,k-1}$ is defined by

where by using the fixed-point Kalman smoothing equations from [34], we get the covariance matrix $\mathbf{V}_{n,k}(k-1)$ as

in which

and we get the mean $\mathbf{m}_{n,k}(k-1)$ in (III-C) as

Similarly, the correlation matrix $\bm{\Gamma}^{n}_{k,k-1}$ in (III-C) is defined as

and the correlation matrix $\bm{\Gamma}^{n}_{k,k}$ in (III-C) is given by using the Kalman filtering time-update equations as

where $\mathbf{V}_{n,k}(k)$ and $\mathbf{m}_{n,k}(k)$ are defined in (15). In fact, note that all the matrices in (III-C)-(III-C) are computed using the KF prediction-update and time-update equations defined in (III-B) and (15) in Section III-B.

Next we compute the M-step of the online EM algorithm in which we maximize the auxiliary function in (III-C) to recursively estimate $\bm{\Theta}=\{\mathbf{Q}_{n},\bm{\Sigma}_{n}\}$. To this end, to estimate $\mathbf{Q}_{n}$, we take the partial derivative of $L_{n,k}\left(\bm{\Theta};\bm{\Theta}^{(k-1)}\right)$ with respect to the matrix $\mathbf{Q}_{n}$ and set it equal to zero. We get the estimate of $\mathbf{Q}_{n}$ in the $k$-th iteration as

where to get the recursive form for updating $\mathbf{Q}_{n}$, we write the numerator in (30) as

thus the update equation for $\mathbf{Q}_{n}$ is written as

in which $\lambda_{k}\triangleq\frac{1-\alpha^{k}}{1-\alpha}$ and we set $\xi(0)=0$. Similarly the matrix $\bm{\Sigma}_{n}$ can also be estimated in the M-step through computing the partial derivatives. Since it is a diagonal matrix as shown in (13), it can be easily estimated by recursively estimating its diagonal elements $\sigma^{m}_{f}$ and $\sigma^{m}_{\theta}$ as follows. To estimate $\sigma^{m}_{f}$, we compute

where $\xi^{f}(0)=0$, and to estimate $\sigma^{m}_{\theta}$, we use