On-the-Fly Communication-and-Computing for Distributed Tensor Decomposition Over MIMO Channels

TL;DR

This paper introduces FlyCom$^2$, an on-the-fly framework for distributed tensor decomposition over MIMO channels, enabling streaming, low-complexity computation with progressive aggregation and improved accuracy.

Contribution

The paper proposes a novel on-the-fly framework for distributed tensor decomposition that integrates streaming computation, progressive aggregation, and a reliability-based sketch selection scheme.

Findings

FlyCom$^2$ achieves comparable or better accuracy than one-shot methods.

The framework reduces device complexity and communication costs.

Error bounds show inverse relation between decomposition error and received sketch observations.

Abstract

Distributed tensor decomposition (DTD) is a fundamental data-analytics technique that extracts latent important properties from high-dimensional multi-attribute datasets distributed over edge devices. Conventionally its wireless implementation follows a one-shot approach that first computes local results at devices using local data and then aggregates them to a server with communication-efficient techniques such as over-the-air computation (AirComp) for global computation. Such implementation is confronted with the issues of limited storage-and-computation capacities and link interruption, which motivates us to propose a framework of on-the-fly communication-and-computing (FlyCom$^2$) in this work. The proposed framework enables streaming computation with low complexity by leveraging a random sketching technique and achieves progressive global aggregation through the integration of progressive uploading and multiple-input-multiple-output (MIMO) AirComp. To develop FlyCom$^2$, an on-the-fly sub-space estimator is designed to take real-time sketches accumulated at the server to generate online estimates for the decomposition. Its performance is evaluated by deriving both deterministic and probabilistic error bounds using the perturbation theory and concentration of measure. Both results reveal that the decomposition error is inversely proportional to the population of sketching observations received by the server. To further rein in the noise effect on the error, we propose a threshold-based scheme to select a subset of sufficiently reliable received sketches for DTD at the server. Experimental results validate the performance gain of the proposed selection algorithm and show that compared to its one-shot counterparts, the proposed FlyCom$^2$ achieves comparable (even better in the case of large eigen-gaps) decomposition accuracy besides dramatically reducing devices' complexity costs.

Equations

$$\mathcal{X}\approx \mathcal{G}{\times }_1{\mathbf{U}}_1{\times }_2{\mathbf{U}}_2\cdots {\times }_N{\mathbf{U}}_N\stackrel{△}{=}\tilde{\mathcal{X}},$$

$${\mathbf{Y}}_t={\mathbf{A}}_t\sum _{k=1}^K{\mathbf{H}}_{t,k}{\mathbf{B}}_{t,k}{\mathbf{S}}_{t,k}^{\top }+{\mathbf{A}}_t{\mathbf{Z}}_t,$$

$${\mathbf{B}}_{t,k}={\left({\mathbf{A}}_t{\mathbf{H}}_{t,k}\right)}^H{\left({\mathbf{A}}_t{\mathbf{H}}_{t,k}{\mathbf{H}}_{t,k}^H{\mathbf{A}}_t^H\right)}^{-1}.$$

$${\mathbf{Y}}_t=\sum _{k=1}^K{\mathbf{S}}_{t,k}^{\top }+{\mathbf{A}}_t{\mathbf{Z}}_t.$$

$$\mathsf{E}\left[\Vert {\mathbf{B}}_{t,k}{\mathbf{S}}_{t,k}^{\top }{\Vert }_F^2\right]$$

$$\le IP.$$

$$d\left({\tilde{\mathbf{U}}}_n,{\mathbf{X}}^{(n)}\right)=\Vert ({\mathbf{I}}_{I_n}-{\tilde{\mathbf{U}}}_n{\tilde{\mathbf{U}}}_n^{\top }){\mathbf{X}}^{(n)}{\Vert }_F^2.$$

$${\mathbf{Y}}_t^{\top }=\sum _k{\mathbf{X}}_k{\mathbf{\Omega }}_{t,k}+{\mathbf{Z}}_t^{\top }{\mathbf{A}}_t^{\top }.$$

$$\mathbf{X}=[{\mathbf{X}}_1,{\mathbf{X}}_2,\cdots ,{\mathbf{X}}_K].$$

$${\mathbf{Y}}_t^{\top }=[{\mathbf{X}}_1,\cdots ,{\mathbf{X}}_K][{\mathbf{\Omega }}_{t,1}^{\top },\cdots ,{\mathbf{\Omega }}_{t,K}^{\top }{]}^{\top }\stackrel{△}{=}\mathbf{X}{\mathbf{F}}_t,$$

$${\mathbf{F}}_t=[{\mathbf{\Omega }}_{t,1}^{\top },\cdots ,{\mathbf{\Omega }}_{t,K}^{\top }{]}^{\top }.$$

$$\tilde{\mathbf{U}}=f(\{{\mathbf{Y}}_{\ell }{\}}_{\ell \le t},\{{\mathbf{A}}_{\ell }{\}}_{\ell \le t}),$$

$${\tilde{\mathbf{Y}}}_t=\Re \{{\mathbf{Y}}_t^{\top }\}=\mathbf{X}{\mathbf{F}}_t+{\tilde{\mathbf{Z}}}_t^{\top }.$$

$$\mathbf{X}={\mathbf{U}}_{\mathbf{X}}{\mathbf{\Sigma }}_{\mathbf{X}}{\mathbf{V}}_{\mathbf{X}}^{\top },$$

$${\hat{\mathbf{Y}}}_t={\mathbf{C}}^{\frac{1}{2}}\mathbf{W}{\mathbf{D}}^{\frac{1}{2}},$$

$$\mathbf{C}\stackrel{△}{=}{\mathbf{U}}_{\mathbf{X}}\mathbf{\Lambda }{\mathbf{U}}_{\mathbf{X}}^{\top },$$

$${\mathbf{C}}^{\frac{1}{2}}={\mathbf{U}}_{\mathbf{X}}{\mathbf{\Lambda }}^{\frac{1}{2}}.$$

$$\boxed{{\mathbf{\Phi }}_t={\hat{\mathbf{Y}}}_t{\mathbf{D}}^{-\frac{1}{2}}={\mathbf{U}}_{\mathbf{X}}{\mathbf{\Lambda }}^{\frac{1}{2}}\mathbf{W}.}$$

$$\mathsf{E}[{\mathbf{S}}_{t,k}^{\top }{\mathbf{S}}_{t,k}]=\mathsf{E}[{\mathbf{\Omega }}_{t,k}^{\top }{\mathbf{X}}_k^{\top }{\mathbf{X}}_k{\mathbf{\Omega }}_{t,k}]=\mathsf{Tr}({\mathbf{X}}_k^{\top }{\mathbf{X}}_k){\mathbf{I}}_M.$$

$$\mathsf{Pr}\left({\mathbf{\Phi }}_t|{\mathbf{U}}_{\mathbf{X}},\mathbf{\Lambda }\right)=\frac{\exp \left(-\frac{tM}{2}\mathsf{Tr}\left({\mathbf{\Phi }}_t^{\top }{\mathbf{U}}_{\mathbf{X}}{\mathbf{\Lambda }}^{-1}{\mathbf{U}}_{\mathbf{X}}^{\top }{\mathbf{\Phi }}_t\right)\right)}{(2\pi {)}^{ItM/2}\mathsf{det}(\mathbf{\Lambda }{)}^{tM/2}}.$$

$$\mathcal{L}\left({\mathbf{U}}_{\mathbf{X}};{\mathbf{\Phi }}_t,\mathbf{\Lambda }\right)=$$

$$=$$

$$-\frac{ItM}{2}\ln (2\pi )-\frac{tM}{2}\ln (\mathsf{det}(\mathbf{\Lambda })).$$

$${\min }_{\tilde{\mathbf{U}}}$$

$$s.t.$$

$${\tilde{\mathbf{U}}}_{\mathbf{X}}=[\tilde{\mathbf{U}},{\tilde{\mathbf{U}}}^{\perp }].$$

$$\mathsf{Tr}\left({\mathbf{\Lambda }}^{-1}{\mathbf{U}}_{\mathbf{X}}^{\top }{\mathbf{\Phi }}_t{\mathbf{\Phi }}_t^{\top }{\mathbf{U}}_{\mathbf{X}}\right)=\sum _{i=1}^I\sum _{j=1}^I{\lambda }_j^{-1}{\gamma }_i({\mathbf{q}}_i^{\top }{\mathbf{u}}_j{)}^2.$$

$${\min }_{\{x_{ij}\}}$$

$$s.t.$$

$$\sum _{i=1}^I\sum _{j=1}^I{\lambda }_j^{-1}{\gamma }_i{x}_{ij}^2\ge \sum _{i=1}^I{\lambda }_i^{-1}{\gamma }_i.$$

$${\tilde{\mathbf{U}}}^{\star }=[{\mathbf{q}}_1,\cdots ,{\mathbf{q}}_r]={\mathcal{S}}_r\left({\mathbf{\Phi }}_t{\mathbf{\Phi }}_t^{\top }\right),$$

$$d(\tilde{\mathbf{U}},\mathbf{X})=\sum _{i=1}^r\sum _{j\ge r+1}({\sigma }_i^2-{\sigma }_j^2)⟨{\tilde{\mathbf{u}}}_i,{\mathbf{u}}_j{⟩}^2+\sum _{i\ge r+1}{\sigma }_i^2.$$

$$({\sigma }_i^2-{\sigma }_j^2)⟨{\tilde{\mathbf{u}}}_i,{\mathbf{u}}_j{⟩}^2\le \max \left\{4,{\delta }_{ij}^2\right\}\frac{\Vert \mathbf{\Delta }{\mathbf{u}}_j{\Vert }_2^2}{{\sigma }_i^2-{\sigma }_j^2},i\le r<j,$$

Full text

I Introduction

In mobile networks, enormous amounts of data are continuously being generated by billions of edge devices. Data analytics can be performed on distributed data to support a broad range of mobile applications, ranging from e-commerce to autonomous driving to IoT sensing [1, 2]. One basic class of techniques is tensor decomposition, which extracts a low-dimensional structure from large-scale multi-attribute data with a tensor representation (a high-dimensional counterpart of a matrix) [3, 4]. A popular technique in this class, the Tucker decomposition, is a higher-dimensional extension of the singular-value decomposition (SVD) that has supported diverse applications such as Google’s image recognition and Cynefin’s spotting of anomalies. In mobile networks, tensor decomposition can be implemented in a centralized manner, which requires uploading of high-dimensional data from many devices to a central server. However, such implementation is stymied not only by a communication bottleneck but also by data privacy issues [5, 6].

In view of these issues, we focus on distributed tensor decomposition (DTD) that avoids direct data uploading and reduces communication overhead by distributing the computation of data tensors to the devices. A direct distributed implementation would call for parallel iterative methods such as alternating least squares [7] and stochastic gradient descent [8, 9] over edge devices, which, however, results in high communication overhead due to slow convergence. On the other hand, DTD can be realized via one-shot distributed matrix analysis techniques [10, 11, 12], since the desired orthogonal factor matrices can be estimated as the principal eigenspaces of unfolding matrices of the tensor along different modes [13]. These one-shot methods improve communication efficiency at a slight cost of decomposition accuracy by following two steps: 1) computing local estimates of the desired factor matrix at each of the devices using local data; 2) uploading and aggregating the local estimates at the server to compute a global estimate. Though alleviated, the communication bottleneck still exists due to the required aggregation of high-dimensional local tensors over potentially many devices. This multi-access problem can be addressed by using a technique called over-the-air computation (AirComp), which exploits the waveform superposition property of a multi-access channel to realize over-the-air data aggregation in one shot [14, 5]. AirComp finds applications in communication-efficient distributed computing and learning, and has been especially popular for federated learning (see, e.g., [15, 16, 17]).

Considering a DTD system with AirComp, this work aims to solve two open problems. The first is the prohibitive cost and latency of computation at resource-constrained devices. A traditional one-shot DTD algorithm requires each device to perform eigenvalue decomposition of a potentially high-dimensional local dataset. The resulting computation complexity increases super-linearly with the data dimensions [18], and the consequent latency makes it difficult for DTD to support emerging mission-critical applications [19]. The second problem is that the one-shot transmissions by devices are susceptible to link disruption. Specifically, a loss of connection during the transmission of high-dimensional local principal components can render the already received partial data useless. In other words, the existing designs lack the feature of graceful performance degradation due to fading.

To solve these problems, we propose the novel framework on-the-fly communication-and-computing (FlyCom2). Underpinning this framework is the use of a technique from randomized linear algebra, randomized sketching, that generates low-dimensional random representations, called sketches, of a high-dimensional data sample by projecting it onto randomly generated low-dimensional sub-spaces [20, 21]. This technique has been successfully used in diverse applications ranging from online data tracking [22] to matrix approximation [21]. In FlyCom2, in place of the traditional high-dimensional local eigenspaces, each device generates a stream of low-dimensional sketches for uploading to the server. Considering a multiple-input-multiple-output (MIMO) channel, the simultaneous transmission of local sketches is enabled by spatially multiplexed AirComp [23]. Upon its arrival at the server, each aggregated sketch is immediately used to improve the global tensor decomposition, giving the name of FlyCom2. Since random sketches serve as independent observations of the tensor, the server can produce an estimate for tensor decomposition in every time slot based on the sketches already received. The FlyCom2 framework addresses the above-mentioned open problems in several aspects. First, random sketching that involves matrix multiplication has much lower complexity than eigen-decomposition and helps reduce the complexity of on-device computation. Second, the DTD accuracy depends on the number of successfully received sketches and hence is robust to loss of sketches in the transmission. This gives FlyCom2 a graceful degradation property in the event of link disruptions or packet losses. Third, as the principal components of a high-dimensional tensor are usually low-dimensional, the progressive DTD at the server is shown to approach its optimal performance quickly as the number of aggregated sketches increases, thereby reining in the communication overhead. Last, the parallel streaming communication and computation in FlyCom2 are more efficient than the sequential operations of the traditional one-shot algorithms due to the communication-computation separation.

In designing the FlyCom2 framework, this work makes the following key contributions.

• •

On-the-Fly Sub-space Detection: One key component of the framework is an optimal on-the-fly detector at the server to estimate the tensor’s principal eigenspace from the received, noisy (aggregated) sketches. To design a maximum likelihood (ML) detector, a whitening technique is used to pre-process the sketches so as to yield an effective global observation, which is shown to have the covariance matrix sharing the same eigenspace as the tensor. Using the result, the ML estimation problem is formulated as a sub-space alignment problem and is solved in closed form. It is observed from the solution that the optimal estimate of the desired principal eigenspace approaches its ground truth as the said observation’s dimensionality grows (or equivalently more sketches are received).

• •

DTD Error Analysis: The end-to-end performance of the FlyCom2-based DTD system is measured by the squared error of the estimated principal eigenspace with respect to (w.r.t.) its ground truth. Using perturbation theory and concentration of measure, bounds are derived on both the error and its expectation. These results reveal that the error consists of one residual component contributed by non-principal components and the other component caused by random sketching. Moreover, the error is observed to be linearly proportional to the number of received sketches, validating the earlier claims on the progressive nature of the designed DTD as well as its feature of graceful degradation. This also suggests a controllable trade-off between the decomposition accuracy and communication overhead, which is useful for practical implementation.

• •

Threshold-Based Sketch Selection: Removing severely channel distorted sketches from use in the sub-space detection can lead to performance improvements. This motivates the design of a sketch-selection scheme that applies a threshold on a scaling factor in MIMO AirComp that reflects the received signal-to-noise ratio (SNR) of an aggregated sketch. We show that such a threshold can be efficiently optimized by an enumeration method whose complexity is polynomial in the population of received sketches.

The remainder of the paper is organized as follows. Section II introduces system models and metrics. Section III gives an overview of the proposed FlyCom2 framework. Then, Section IV presents the design of the on-the-fly sub-space estimator and its error analysis. The sketch-selection scheme is proposed in Section V. Numerical results are provided in Section VI, followed by concluding remarks in Section VII.

II Models, Operations, and Metrics

We consider the support of DTD in a MIMO system, as illustrated in Fig. 1. The relevant models, operations and metrics are described in what follows.

II-A Distributed Tensor Decomposition

We consider the distributed implementation of the popular Tucker method for tensor decomposition [13]. For ease of notation, the tensor is assumed to have $N$ modes; these modes generalize the concepts of columns and rows in matrices with the first $(N-1)$ modes corresponding to data features and mode $N$ indexing data samples. For instance, in a surveillance system, images captured by multiple cameras are expressed as local tensors with three modes indicating pixels, colors, and data sample indices, respectively. Let the samples collected by device $k$ be represented by a local tensor $\mathcal{X}_{k}\in\mathbb{R}^{I_{1}^{(k)}\times I_{2}^{(k)}\cdots\times I_{N}^{(k)}}$, where $I_{n}^{(k)}$ denotes the dimensionality of mode $n$ of local tensor $k$. To simplify notation, we assume that local tensors have the same dimensions for their feature modes: $I_{n}^{(k)}=I_{n}$, $\forall k,1\leq n\leq N-1$. Next, these local tensors are aggregated from $K$ devices to form a global tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times I_{N}}$ with $I_{N}=\sum_{k}I_{N}^{(k)}$. The Tucker decomposition of $\mathcal{X}$ can be written as [13]

[TABLE]

where $\mathcal{G}\in\mathbb{R}^{r_{1}\times r_{2}\cdots\times r_{N}}$ represents a core tensor [that generalizes the singular values in SVD], $\mathbf{U}_{n}\in\mathbb{R}^{I_{n}\times r_{n}}$ is an orthogonal factor matrix corresponding to the $n$-th mode, satisfying $\mathbf{U}_{n}^{\top}\mathbf{U}_{n}=\mathbf{I}_{r_{n}}$ with $r_{n}$ ($r_{n}\leq I_{n}$) representing the number of principal dimensions, and $\times_{n}$ denotes the mode-$n$ matrix product [13]. In the sequel, we pursue these factor matrices $\{\mathbf{U}_{n}\}$ as they reveal the characteristics of the data tensor at different modes. Given $\{\mathbf{U}_{n}\}$, the computation of $\mathcal{G}$ is straightforward [3]. In centralized computation with full data aggregation, $\{\mathbf{U}_{n}\}$ can be computed by using the higher-order SVD approach [3, 20]. In this approach, the tensor is first flattened along a chosen mode $n$ to yield a matrix $\mathbf{X}^{(n)}\in\mathbb{R}^{I_{n}\times J_{n}}$ with $J_{n}=\prod_{j=1,j\neq n}^{N}I_{n}$, termed mode-$n$ unfolding; then the desired factor matrix is computed as $\mathbf{U}_{n}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{r_{n}}]$ where $\mathbf{u}_{i}$ is given by the $i$-th principal eigenvector of the mode-$n$ unfolding. Let this operation be represented by $\mathcal{S}_{r_{n}}(\cdot)$ and hence $\mathbf{U}_{n}=\mathcal{S}_{r_{n}}(\mathbf{X}^{(n)}(\mathbf{X}^{(n)})^{\top})$.

In contrast with its centralized counterpart, DTD computes the eigenspaces of different unfolding matrices distributively. This avoids the aggregation of raw data and preserves the data ownership [11, 10]. Considering the computation of $\mathbf{U}_{n}$, DTD goes through the following procedure: 1) local tensors are flattened along chosen mode $n$ to generate local unfoldings, denoted by $\{\mathbf{X}_{k}^{(n)}\}$; 2) devices compute low-dimensional component $\{\mathbf{S}_{k}\}$ from local unfoldings $\{\mathbf{X}_{k}^{(n)}\}$ through dimensionality reduction techniques; 3) the server gathers these local components from devices and aggregates them into a global component, denoted by $\mathbf{S}$, to yield a global estimate of the ground truth, $\mathbf{U}_{n}$. It is worth mentioning that the computation results $\{\mathbf{S}_{k}\}$ depend on a particular dimensionality reduction technique. For example, when using principal component analysis (PCA) [11, 10], $\{\mathbf{S}_{k}\}$ are computed as the principal eigenspaces of $\{\mathbf{X}_{k}^{(n)}\}$ at the devices, and then the server averages them to estimate $\mathbf{U}_{n}$. In this work, a random sketching approach is adopted as elaborated in Section III.

II-B MIMO Over-the-Air Computation

FlyCom2 builds on MIMO AirComp to aggregate local results over the air, which is described as follows. First, let $N_{\text{r}}$ and $N_{\text{t}}$, with $N_{\text{r}}\geq N_{\text{t}}$, denote the numbers of antennas at the edge server and each device, respectively. We assume perfect transmit channel state information (CSI) as well as symbol-level and phase synchronization between the devices [23]. Time is slotted and then grouped to form coherence blocks with $t$ denoting the block index. In each time slot, an $N_{\text{t}}\times 1$ vector of complex scalar symbols is transmitted over $N_{\text{t}}$ antennas. Then a coherence block spans at least $I$ symbol slots to support the transmission of an $N_{\text{t}}\times I$ matrix. In an arbitrary block, say $t$, all edge devices transmit simultaneously their $I\times M$ real matrices, denoted as $\{\mathbf{S}_{t,k}\}$, each of which is termed a matrix symbol. As a result, the server receives an over-the-air aggregated matrix symbol, $\mathbf{Y}_{t}$, as

[TABLE]

where $\mathbf{H}_{t,k}\in\mathbb{C}^{N_{\text{r}}\times N_{\text{t}}}$ denotes the channel matrix corresponding to device $k$, $\mathbf{Z}_{t}$ models additive Gaussian noise with independent and identically distributed (i.i.d.) elements of $\mathcal{CN}(0,\sigma^{2})$, $\mathbf{A}_{t}\in\mathbb{C}^{M\times N_{\text{r}}}$ and $\mathbf{B}_{t,k}\in\mathbb{C}^{N_{\text{t}}\times M}$ ($M\leq N_{\text{t}}$) denote receive and transmit beamforming matrices, respectively. To realize AirComp, we consider zero forcing (ZF) transmit beamforming that inverts individual MIMO channels [23]. Mathematically, conditioned on a fixed receive beamformer, transmit beamforming matrices are given as

[TABLE]

The received matrix $\mathbf{Y}_{t}$ is then rewritten as

[TABLE]

In the absence of noise, the AirComp in (3) provides the one-shot realization of the desired aggregation operation for DTD. The average transmission power of each device is constrained to not exceed a power budget of $P$ per slot, i.e. $\forall t,k$

[TABLE]

The transmit SNR is then given by $\gamma=\frac{P}{\sigma^{2}}$.

II-C Error Metric

With the above-described MIMO AirComp, a noisy version of $\mathbf{U}_{n}$, denoted by $\tilde{\mathbf{U}}_{n}$, will be computed progressively from a set of received matrices $\{\mathbf{Y}_{t}\}$ (see Section III). The $\tilde{\mathbf{U}}_{n}$ deviates from the ground truth due to both distributed computation and channel noise. The resulting error can be a performance metric of DTD in the wireless system. Mathematically, given $\tilde{\mathcal{X}}$ as the tensor derived from $\{\tilde{\mathbf{U}}_{n}\}$, the error is measured as $\|\mathcal{X}-\tilde{\mathcal{X}}\|_{F}^{2}$ that can be bounded as $\|\mathcal{X}-\tilde{\mathcal{X}}\|_{F}^{2}\leq\sum_{n=1}^{N}\|(\mathbf{I}_{I_{n}}-\tilde{\mathbf{U}}_{n}\tilde{\mathbf{U}}_{n}^{\top})\mathbf{X}^{(n)}\|_{F}^{2}$ [20]. This suggests that the overall error is determined by the error of independently decomposing each unfolding matrix. Hence, we define the DTD error as

[TABLE]

III Overview of On-the-Fly Communication-and-Computing

To support DTD over edge devices with limited computation power, we propose a FlyCom2 framework as shown in Fig. 1. Next, we first briefly introduce the random approach exploited in FlyCom2 and then explain how to use FlyCom2 to support DTD.

III-A Data Dimensionality Reduction via Random Sketching

Recall that the DTD requires data dimensionality reduction on devices prior to transmission. For high-dimensional tensors, the traditional PCA technique becomes too complex for resource-constrained devices. To address this issue, we adopt a technique for random dimensionality reduction, known as random sketching, which is simpler than PCA as it only relies on matrix multiplication and also requires a smaller number of passes over the datasets [21]. Specifically, given an $I\times J$ data matrix $\mathbf{X}$, random sketching uses a $J\times M$ random matrix, termed dimension reduction mapping (DRM) and denoted by $\mathbf{\Omega}$, to map $\mathbf{X}$ to an $I\times M$ sketch matrix $\mathbf{S}$ with $J\gg M$: $\mathbf{S}=\mathbf{X}\mathbf{\Omega}$. The mapping $\mathbf{\Omega}$ can be composed of i.i.d. Gaussian elements and projects the high-dimensional $\mathbf{X}$ to random directions in a space of low dimensionality. Despite the random projection, the mutual vector distances between the rows of $\mathbf{X}$ can be approximately preserved such that the principal (column) eigenspace of the sketch, $\mathbf{S}$, constitutes a good approximation of $\mathbf{X}$. The approximation accuracy grows as $M$ increases and becomes perfect when $M$ is equal to $J$ [21]. Importantly, to estimate an $r$-dimensional principal eigenspace, random sketching has the complexity of $\mathcal{O}(IJM)$ and requires a single data pass of memory, as opposed to the complexity of $\mathcal{O}(\min\{I,J\}^{2}\times\max\{I,J\})$ and $\mathcal{O}(r)$ memory passes in PCA [18].

III-B FlyCom2-Based DTD

Based on the preceding random-sketching technique, we propose the FlyCom2 framework that decomposes the high-dimensional DTD into on-the-fly processing and transmission of streams of low-dimensional random sketches. Thereby, we not only overcome devices’ resource constraints but also achieve a graceful reduction of DTD error as the communication time increases. Without loss of generality, we focus on the computation of the principal eigenspace $\mathbf{U}_{n}$ for an arbitrary data-feature mode $n$ with $n\in[1,\cdots,N-1]$. To simplify notation, the superscript $(n)$ and subscript $n$ are omitted. The detailed operations of FlyCom2 are described as follows.

III-B1 On-the-Fly Computation at Devices

Each device streams a sequence of low-dimensional local sketches to the server by generating and transmitting them one by one in data packets. First, the progressive computation of local sketches at devices is introduced as follows. Let each local tensor, say $\mathcal{X}_{k}$ at device $k$, be flattened along the desired mode to generate the unfolding matrix $\mathbf{X}_{k}$. Then, in the (matrix-symbol) slot $t$, each device $k$ draws i.i.d. $\mathcal{N}(0,1)$ entries to form a $J\times M$ DRM, denoted by $\mathbf{\Omega}_{t,k}$, or retrieves it efficiently from memory [24]. Then an $M$-dimensional local sketch for $\mathbf{X}_{k}$ can be computed as $\mathbf{S}_{t,k}=\mathbf{X}_{k}\mathbf{\Omega}_{t,k}$, which is then uploaded to the server immediately before computing the next sketch $\mathbf{S}_{t+1,k}$. This allows efficient communication-and-computation parallelization as shown in Fig. 2.

III-B2 On-the-Fly Global Random Sketching

MIMO AirComp is used for low-latency aggregation of the local sketches simultaneously streamed by devices. Local temporal sketches are progressively aggregated at the server by linearly modulating them as MIMO AirComp symbols. Consider the uploading of the $t$-th local sketches. It follows from (3) that the matrix symbol received at the server can be written as

[TABLE]

To explain how to use $\mathbf{Y}_{t}$ in estimating the principal eigenspace of the global unfolding matrix $\mathbf{X}$, we first consider the case without channel noise, in which $\mathbf{Y}_{t}^{\top}=\sum_{k}\mathbf{X}_{k}\mathbf{\Omega}_{t,k}$. Since the global tensor $\mathcal{X}$ is given by assembling local tensors along mode $N$, the corresponding global unfolding matrix, denoted by $\mathbf{X}$, is related to the local unfoldings $\{\mathbf{X}_{k}\}$ as

[TABLE]

It follows that

[TABLE]

where we define

[TABLE]

As $\{\mathbf{\Omega}_{t,k}\}$ are mutually independent, $\mathbf{F}_{t}$ has i.i.d. $\mathcal{N}(0,1)$ elements and can be used as an $M$-dimensional DRM for randomly sketching $\mathbf{X}$. Therefore, in the absence of channel noise, $\mathbf{Y}_{t}$ gives an $M$-dimensional global sketch for $\mathbf{X}$. The dimension of the global sketch grows, thereby improving the DTD accuracy, as more aggregated local sketches are received (or equivalently as $t$ progresses), giving the name of on-the-fly global sketching.

III-B3 On-the-Fly Sub-space Detection at the Server

In the case with channel noise, the server can produce an estimate of the desired principal eigenspace, $\mathbf{U}$, based on the noisy observations accumulated up to the current symbol slot. Specifically, in slot $t$, given the current and past received matrix symbols, $\{\mathbf{Y}_{\ell}\}_{\ell\leq t}$, and the receive beamformers $\{\mathbf{A}_{\ell}\}_{\ell\leq t}$ (discussed in the sequel), the server estimates $\mathbf{U}$ as

[TABLE]

where the estimator $f(\cdot)$ is optimized in the sequel to minimize the DTD error in (5).

Following the above discussion, the procedure for FlyCom2-based DTD is summarized as follows.

The key component of the FlyCom2 framework, the on-the-fly sub-space estimator $f(\cdot)$, is designed in Section IV. The performance of FlyCom2-based DTD is enhanced using a sketch selection algorithm designed in Section V.

IV Optimal Sub-space Detection for FlyCom2

In this section, we design the sub-space detection function of the FlyCom2 framework, namely $f(\cdot)$ mentioned in the preceding section. It consists of two stages – pre-processing of received symbols and the subsequent sub-space estimation, which are summarized in Algorithm 1 and designed in the following sub-sections. Furthermore, the resultant DTD error is analyzed.

IV-A Pre-Processing of Received Matrix Symbols

The pre-processing function is to accumulate received matrix symbols from slot $1$ to the current slot, $t$, and generate from them an effective matrix for the ensuing sub-space detection. The operation is instrumental for on-the-fly detection to obtain a progressive performance improvement. The design of the pre-processing takes several steps. First, since the transmitted symbol $\mathbf{X}\mathbf{F}_{t}$ is real but the channel noise is complex, the real part of the received symbols, namely $\mathbf{Y}_{t}$ in (6), gives an effective observation of the transmitted symbol111It is possible to transmit the coefficients of $\mathbf{X}\mathbf{F}_{t}$ over both the in-phase and quadrature channels, which halves air latency. The extension is straightforward (see, e.g. [10, Section II]) but complicates the notation without providing new insights. Hence, only the in-phase channel is used in this work.. Let $\tilde{\mathbf{Y}}_{t}$ denote the effective observation in slot $t$ and $\tilde{\mathbf{Z}}_{t}$ the real part of $\mathbf{A}_{t}\mathbf{Z}_{t}$. It follows that

[TABLE]

Second, the relation between the eigenspace of $\mathbf{X}$ and the accumulated observations up to the current slot is derived as follows. To this end, let the SVD of $\mathbf{X}$ be expressed as

[TABLE]

where $\mathbf{\Sigma}_{\mathbf{X}}$ comprises descending singular values along its diagonal. Then, the accumulation of the current and past observations, denoted by $\hat{\mathbf{Y}}_{t}=[\tilde{\mathbf{Y}}_{1},\tilde{\mathbf{Y}}_{2},\cdots,\tilde{\mathbf{Y}}_{t}]$, is a random Gaussian matrix as shown below.

Lemma 1**.**

The accumulated aggregations, $\hat{\mathbf{Y}}_{t}$, can be decomposed as

[TABLE]

where the left covariance matrix $\mathbf{C}=\mathbf{X}\mathbf{X}^{\top}+\frac{1}{2tM}\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})\mathbf{I}_{I}$, the right covariance matrix $\mathbf{D}=\frac{\mathsf{Tr}(\mathbf{X}^{\top}\mathbf{X})\mathbf{I}_{tM}+\frac{1}{2}I\sigma^{2}\mathsf{diag}(\mathbf{A}_{1}\mathbf{A}_{1}^{H},\cdots,\mathbf{A}_{t}\mathbf{A}_{t}^{H})}{\mathsf{Tr}(\mathbf{X}^{\top}\mathbf{X})+\frac{1}{2tM}I\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})}$, and $\mathbf{W}$ is a random Gaussian matrix with i.i.d. $\mathcal{N}(0,1)$ entries.**

Proof.

See Appendix -A ∎

Third, based on (10), the covariance matrix, $\mathbf{C}$, in Lemma 1 can be rewritten as

[TABLE]

where we define $\mathbf{\Lambda}=\mathbf{\Sigma}_{\mathbf{X}}^{2}+\frac{1}{2tM}\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})\mathbf{I}_{I}$. Hence, the square root, $\mathbf{C}^{\frac{1}{2}}$, is given as

[TABLE]

Remark 1** (Effective Sketching with Channel Noise).**

According to Lemma 1 and (12), the accumulated observations, $\hat{\mathbf{Y}}_{t}$, gives a sketch of the matrix $\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}^{\frac{1}{2}}$ using a Gaussian DRM with the covariance of $\mathbf{D}$. The matrix $\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}^{\frac{1}{2}}$ and the unfolding matrix $\mathbf{X}$ share the eigenspace, $\mathbf{U}_{\mathbf{X}}$. Furthermore, as $\mathbf{\Lambda}$ retains the singular values in descending order, the top-$r$ principal eigenspace of $\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}^{\frac{1}{2}}$ is identical to that of $\mathbf{X}$ for any $1\leq r\leq M$.**

Finally, according to the preceding discussion, the desired principal eigenspace of $\mathbf{X}$ can be estimated from the sketch $\hat{\mathbf{Y}}_{t}$. It is known that randomized sketching prefers DRMs with i.i.d. entries [21]. To improve the performance, $\hat{\mathbf{Y}}_{t}$ can be further “whitened” to equalize the right covariance $\mathbf{D}$. Specifically, let $\hat{\mathbf{Y}}_{t}$ be right-multiplied by $\mathbf{D}^{-\frac{1}{2}}$ to yield the final effective observation in time slot $t$ as

[TABLE]

To compute the covariance matrix $\mathbf{D}$, the server needs to acquire the value of $\mathsf{Tr}(\mathbf{X}\mathbf{X}^{\top})=\sum_{k}\mathsf{Tr}(\mathbf{X}_{k}\mathbf{X}_{k}^{\top})$. Note that each term in the summation, say $\mathsf{Tr}(\mathbf{X}_{k}\mathbf{X}_{k}^{\top})$, relates to the covariance of transmitted symbols, $\mathsf{E}[\mathbf{S}_{t,k}^{\top}\mathbf{S}_{t,k}]$, as

[TABLE]

Then, $\mathsf{Tr}(\mathbf{X}_{k}^{\top}\mathbf{X}_{k})$ can be acquired at the server by one-time feedback.

IV-B Optimal Sub-space Estimation

In this sub-section, the principal eigenspace of the unfolding matrix $\mathbf{X}$ with dimensions fixed as $r$, is estimated from the effective observation given in (13) under the ML criterion. First, using (13), the distribution of the observation $\mathbf{\Phi}_{t}$ conditioned on $\mathbf{U}$ and $\mathbf{\Lambda}$ is given as

[TABLE]

This yields the logarithm of the likelihood function, required for ML estimation, as

[TABLE]

Let $\mathbf{U}$ denote the desired $r$-dimensional principal components of $\mathbf{X}$, as obtained from splitting $\mathbf{U}_{\mathbf{X}}=[\mathbf{U},\mathbf{U}^{\bot}]$. It is observed from (IV-B) that only the first term depends on the variable $\mathbf{U}$. Then letting $(\tilde{\mathbf{U}},\tilde{\mathbf{U}}_{\mathbf{X}})$ denote an estimate of $(\mathbf{U},\mathbf{U}_{\mathbf{X}})$, the ML-estimation problem can be formulated as

[TABLE]

Despite the non-convex orthogonality constraints, the problem in (15) can be solved optimally in closed form, as follows. First, define the eigenvalue decomposition $\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}=\mathbf{Q}\mathbf{\Gamma}\mathbf{Q}^{T}$ with $\mathbf{Q}=[\mathbf{q}_{1},\cdots,\mathbf{q}_{I}]$ and $\mathbf{\Gamma}=\mathsf{diag}(\gamma_{1},\cdots,\gamma_{I})$ with eigenvalues arranged in a descending order. Then, given $\mathbf{\Lambda}=\mathsf{diag}(\lambda_{1},\cdots,\lambda_{I})$ and $\mathbf{U}_{\mathbf{X}}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{I}]$, the objective function of (15) can be rewritten as

[TABLE]

Next, define $x_{ij}=\mathbf{q}_{i}^{\top}\mathbf{u}_{j}$ and rewrite the constraints in (15) as $\sum_{i=1}^{I}x_{ij}^{2}=\mathbf{u}_{j}^{\top}\mathbf{Q}\mathbf{Q}^{\top}\mathbf{u}_{j}=1$ and $\sum_{j=1}^{I}x_{ij}^{2}=\mathbf{q}_{i}^{\top}\mathbf{U}\mathbf{U}^{\top}\mathbf{u}_{i}=1$. Without loss of optimality, such constraints can be further relaxed as $\sum_{i=1}^{I}x_{ij}^{2}\geq 1$ and $\sum_{j=1}^{I}x_{ij}^{2}\geq 1$. This allows the problem in (15) to be reformulated as a convex problem:

[TABLE]

Since $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{I}$, the objective of (16) subject to the constraints is lower bounded as

[TABLE]

The lower bound can be achieved by letting $x_{ii}=1$, $\forall i$ and $x_{ij}=0$, $\forall i\neq j$. The optimal solution for (16) follows as shown below.

Proposition 1**.**

Based on the ML criterion, in slot $t$, the optimal on-the-fly estimate of the $r$-dimensional principal components of the unfolding matrix, $\mathbf{X}$, is denoted as $\tilde{\mathbf{U}}^{\star}$ and given as

[TABLE]

where $\mathbf{\Phi}_{t}$ is the effective observation in slot $t$ as given in (13) and we recall $\mathcal{S}_{r}(\cdot)$ to yield the $r$-dimensional principal eigenspace of its argument.**

Remark 2** (Minimum Number of FlyCom2 Operations).**

For the result in (18) to hold, the dimensions of the current effective observations $\mathbf{\Phi}_{t}$ should be larger than those of $\mathbf{U}$, i.e. $tM\geq r$. This implies that the FlyCom2 should run at least $t\geq r/M$ rounds to enable the estimation of an $r$-dimensional principal eigenspace of the tensor.**

IV-C DTD Error Analysis

Based on the optimal sub-space detection designed in the preceding sub-section, we mathematically quantify the key feature of FlyCom2 that the DTD error gracefully decreases with communication rounds. The existing error analysis for random sketching does not target distributed implementation and hence requires no communication links [21, 20]. The new challenge for the current analysis arises from the need to account for the distortion increased by the MIMO AirComp transmission. In what follows, we derive deterministic and probabilistic bounds on the DTD error defined in (5).

IV-C1 Deterministic Error Bound

As the unfolding matrix comprises $r$ principal components, its singular values can be represented as $\mathbf{\Sigma}_{\mathbf{X}}=\mathsf{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{I})$ with $\sigma_{1}=\cdots=\sigma_{r}\gg\sigma_{r+1}\geq\cdots\geq\sigma_{I}$, where we assume the same principal singular values following the literature (see, e.g. [20]).

Lemma 2**.**

Consider the DTD of the unfolding matrix $\mathbf{X}$ in tensor decomposition that has an $r$-dimensional principal eigenspace $\mathbf{U}=[\mathbf{u}_{1},\cdots,\mathbf{u}_{r}]$ and the singular values $\mathbf{\Sigma}_{\mathbf{X}}$. The estimation of $\mathbf{U}$ as in (18) yields the DTD error given as

[TABLE]

Proof.

See Appendix -B. ∎

On the right hand side, the first term, $\sum_{i=1}^{r}\sum_{j\geq r+1}(\sigma_{i}^{2}-\sigma_{j}^{2})\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle^{2}$, represents the error due to random sketching; the second term $\sum_{i\geq r+1}\sigma_{i}^{2}$ represents the residual error due to non-zero non-principal components of $\mathbf{X}$.

Next, we make an attempt to characterize the behavior of each error term, $(\sigma_{i}^{2}-\sigma_{j}^{2})\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle^{2}$. Let $\tilde{\mathbf{u}}_{i}$ and $\mathbf{u}_{j}$ denote the $i$-th and $j$-th ($i\leq r<j$) eigenvectors of the sample covariance matrix $\frac{1}{tM}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}$ and the covariance matrix $\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}\mathbf{U}_{\mathbf{X}}^{\top}$, respectively. The error in (19) is caused by the perturbation $\mathbf{\Delta}=\frac{1}{tM}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}-\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}\mathbf{U}_{\mathbf{X}}^{\top}$. Using this fact allows us to obtain the following desired result.

Lemma 3**.**

Consider a fixed realization $\mathbf{W}$ in the DRM, $\mathbf{\Phi}_{t}$, in (13) and the error term, $(\sigma_{i}^{2}-\sigma_{j}^{2})\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle^{2}$, in Lemma 2 is upper bounded as

[TABLE]

*where $\delta_{ij}\overset{\triangle}{=}\frac{\min\{2|\tilde{\lambda}_{i}-\lambda_{i}|,(\sigma_{i}^{2}-\sigma_{j}^{2})\}}{|\tilde{\lambda}_{i}-\lambda_{j}|}$ with $\lambda_{i}$ and $\tilde{\lambda}_{i}$ being the $i$-th eigenvalues of $\mathbf{U}_{\mathbf{X}}\mathbf{\Lambda}\mathbf{U}_{\mathbf{X}}^{\top}$ and $\frac{1}{tM}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}$, respectively. ***

Proof.

See Appendix -C. ∎

The upper bound in Lemma 3 suggests two scaling regions of the DTD error, namely $\delta_{ij}\geq 2$ and $\delta_{ij}<2$. Invoking the well-known Weyl’s theorem (see, e.g. [25]), the norm of the perturbation $\mathbf{\Delta}$ and hence the value of $\delta_{ij}$ reduce as FlyCom2 progresses in time. This is aligned with the result in Fig. 3(a) where the average value of $\{\delta_{ij}\}$ is observed to decrease with increasing communication time $t$. To simplify the analysis, we focus on the case of $\delta_{ij}\leq 2$, $\forall i\leq r<j$, by assuming sufficiently large $t$. In this case, the upper bound in Lemma 3 simplifies to

[TABLE]

Next, based on Lemma 2 and (20), the desired deterministic error bound is derived as follows.

Theorem 1** (Expected Error Bound).**

Given the receive beamformers $\{\mathbf{A}_{\ell}\}_{\ell\leq t}$ of FlyCom2-based DTD and $\delta_{ij}\leq 2$, $\forall i\leq r<j$, the expected error can be bounded as

[TABLE]

where $\delta_{ij}$ and $\lambda_{j}$ follow those in Lemma 3.**

Proof.

See Appendix -D. ∎

The error bound in the Theorem 1 is compared numerically with the exact error and that of centralized tensor decomposition in Fig. 3(b). One can observe the bound to capture the trend of decreasing DTD error as $t$ progresses. In particular, it shows that under a small perturbation,

[TABLE]

IV-C2 Probabilistic Error Bound

We derive in the sequel a probabilistic bound on the DTD error using the method of concentration of measure. A relevant useful result is given below.

Lemma 4** (McDiarmid’s Inequality [26]).**

Let $g$ be a positive function on independent variables $\{W_{m}\}$ satisfying the bounded difference property:

[TABLE]

with constants $\{c_{m}\}$ and $W_{M}$ being i.i.d. as $W_{m}$. Then, for any $\epsilon>0$,

[TABLE]

Using Lemma 4, the desired result is obtained as shown below.

Theorem 2** (Probabilistic Error Bound).**

Given receive beamformers $\{\mathbf{A}_{\ell}\}_{\ell\leq t}$ and $\delta_{ij}\leq 2$, $\forall i\leq r<j$, for any $\epsilon\geq 0$, the error of the FlyCom2-baded DTD can be upper bounded as

[TABLE]

with the probability of at least $\left[1-e^{-\frac{\epsilon^{2}}{2\kappa^{8}}}\right]\mathsf{erf}\left(\frac{\kappa}{\sqrt{2}}\right)^{tM(I-r)}$, where $\mathsf{erf}(\cdot)$ denotes the error function defined as $\mathsf{erf}(y)=\frac{2}{\sqrt{\pi}}\int_{0}^{y}e^{-x^{2}}\mathrm{d}x$ and $\kappa\geq 1$.**

Proof.

See Appendix -E. ∎

The upper bound on the DTD error in Theorem 2 holds almost surely if the constant $\epsilon$ is sufficiently large. Comparing Theorem 1 and Theorem 2, one can make the important observation that as the communication time ($t$) progresses, both the DTD error and its expectation vanish at the same rate of

[TABLE]

Another observation is that non-principal components of the tensor contribute to the DTD error but the effect is negligible when the eigen-gap is large.

V Optimal Sketch Selection for FlyCom2

Discarding aggregated sketches that have been transmitted under unfavourable channel conditions can improve the FlyCom2 performance. This motivates us to design a sketch-selection scheme in this section.

V-A Threshold Based Sketch Selection

First, we follow the approach in [23] to design the receive beamforming, $\{\mathbf{A}_{t}\}$, for MIMO AirComp. To this end, we decompose $\mathbf{A}_{t}$ as $\mathbf{A}_{t}=\eta_{t}\mathbf{U}_{\mathbf{A}_{t}}$, where the positive scalar $\eta_{t}$ is called a denoising factor and $\mathbf{U}_{\mathbf{A}_{t}}$ is an $M\times N_{\text{r}}$ unitary matrix. Following similar steps as in [23], we can show that to minimize the DTD error bounds in Theorem 1 and 2, the beamformer component should be aligned with the channels of devices as

[TABLE]

where $\lambda_{\mathbf{H}_{t,k}}$ and $\mathbf{U}_{\mathbf{H}_{t,k}}$ denote the $N_{\text{t}}$-th eigenvalue and the first $N_{\text{t}}$ eigenvectors of $\mathbf{H}_{t,k}\mathbf{H}_{t,k}^{\top}$, respectively. Furthermore, the denoising factor $\eta_{t}$ should cope with the weakest channel by being

[TABLE]

It follows from (22) and (23) that $\mathsf{Tr}(\mathbf{A}_{t}^{H}\mathbf{A}_{t})=\eta_{t}M$ and $\lambda_{j}$ in DTD error bounds in Theorem 1 and 2 can be expressed as

[TABLE]

which shows that the error relies on only the denoising factor up to the current time slot. The result also suggests that it is preferable to select from the received sketches $\{\tilde{\mathbf{Y}}_{n}\}_{\ell\leq t}$ those associated with small $\eta_{\ell}$ that reflects a favorable channel condition. Naturally, we can derive a threshold based selection scheme as follows:

[TABLE]

where the threshold $\eta_{\mathrm{th}}$ is optimized in the sequel.

V-B Threshold Optimization

The threshold, $\eta_{\mathrm{th}}$, in (25), needs to be optimized to minimize the error in (21). Solving the problem is hindered by that the singular values $\{\sigma_{j}\}$ are not available at the server in advance. We tackle this problem by designing a practical optimization scheme. To this end, we resort to using an upper bound on the DTD error as shown below.

Lemma 5**.**

Let $\tilde{M}$ denote the number of aggregated sketches selected from $\{\tilde{\mathbf{Y}}_{\ell}\}_{\ell\leq t}$ based on (25) with the threshold $\eta_{\mathrm{th}}$, the DTD error in (21) satisfies

[TABLE]

*where $c$ is a constant.

---

Proof.

See Appendix -F. ∎

Lemma 5 suggests that a sub-optimal threshold can be obtained by minimizing the error upper bound. Let $S$ denote a set and $|S|$ its cardinality. Then, the threshold-optimization problem can be formulated as

[TABLE]

where $\eta_{\ell}$ follows the definition in (23). One can observe that the objective function of (26) is a monotonically increasing function w.r.t. the variable $\eta_{\mathrm{th}}$ if $\tilde{M}$ is fixed. Such piecewise monotonicity of the objective function (26) renders a linear-search solution method (e.g. bisection search) infeasible, but allows the optimal solution, $\eta_{\mathrm{th}}^{\star}$, to be restricted into a finite set, say $\eta_{\mathsf{th}}^{\star}\in\left\{\eta_{1},\eta_{2},\cdots,\eta_{t}\right\}$. Then, finding $\eta_{\mathsf{th}}^{\star}$ is simple by exhausted enumeration as follows. Let $\tilde{M}_{\ell}$ denote the number of selected sketches corresponding to the threshold fixed as $\eta_{\mathsf{th}}=\eta_{\ell}$. Define

[TABLE]

The optimal threshold solving the problem in (26) is $\eta_{\mathsf{th}}=\eta_{\ell^{\star}}$. Two remarks are offered as follows. First, the above research for the optimal threshold has complexity linearly proportional to $t$, the population of accumulated sketches at the server. Second, the implementation of the optimization at the server requires feedback of a scalar from each device, namely $\mathsf{Tr}\left(\mathbf{X}_{k}^{\top}\mathbf{X}_{k}\right)$ from device $k$.

VI Experimental Results

VI-A Experimental Settings

First, the MIMO AirComp system is configured to have the following settings. There are $K=20$ edge devices connected to the server. The array sizes at each device and the server are set as $N_{\text{t}}=4$ and $N_{\text{r}}=16$, respectively. The Rayleigh channel with shadow fading is adopted, in which each MIMO channel is given as $\mathbf{H}_{t,k}=\sqrt{\beta_{t,k}}\hat{\mathbf{H}}_{t,k}$ with $\beta_{t,k}$ following a Gamma distribution $\Gamma(1.2,0.83)$ (see, e.g. [27]) and $\hat{\mathbf{H}}_{t,k}$ comprising i.i.d. $\mathcal{CN}(0,1)$ entries. Different channels are independent. Second, we use a synthetic data model following the DTD literature (see, e.g. [20]). Considering the computation of the $n$-th factor matrix, the unfolding matrix of the data tensor has the size of $I_{n}\times(\prod_{j=1,j\neq n}^{N})=100\times 1500$, and its columns are uniformly distributed over devices. Under such settings, each local sketch has the length of $I_{n}=100$ that is smaller than a single channel coherence block (see, e.g. [28], for justification). To demonstrate the performance of the proposed FlyCom2 for data with a range of parameterized spectral distribution, the singular values of the unfolding matrix are set to decay with a polynomial rate:

[TABLE]

where the first $r=12$ principal singular values are fixed as $1$ and $\xi>0$ controls the decay rate of residual values. Furthermore, the left and right eigenspaces of the unfolding matrix are generated as those of random matrices with i.i.d. $\mathcal{N}(0,1)$ entries [20].

Third, we consider two benchmarking schemes that are variants of SVD-based DTD.

• •

Centroid SVD-DTD: Devices compute local principal eigenspaces $\{\hat{\mathbf{U}}_{k}\}$ of their on-device data samples by using SVD and the server then aggregates these local results as $\mathbf{P}=\frac{1}{K}\sum_{k}\hat{\mathbf{U}}_{k}\hat{\mathbf{U}}_{k}^{\top}$. The principal eigenspace of $\mathbf{P}$ represents the centroid of all local estimates $\{\hat{\mathbf{U}}_{k}\}$ on the Grassmannian manifold and is extracted to form a global estimate of the ground truth [11, 10].

• •

Alignment SVD-DTD: The scheme follows a similar procedure as above except for aggregating local results, $\{\hat{\mathbf{U}}_{k}\}$, as $\mathbf{P}=\frac{1}{K}\sum_{k}\hat{\mathbf{U}}_{k}\mathbf{J}_{k}$, where the orthogonal matrices $\{\mathbf{J}_{k}\}$ are alignment matrices that are to be optimized by using past global estimates to improve the system performance [12].

The aggregation operations in both benchmark schemes are implemented using MIMO AirComp [14, 23] as FlyCom2 for fair comparison.

VI-B Performance Gain of Sketch Selection

In Fig. 4, we compare the error performance of FlyCom2 between the cases with and without sketch selection. The communication time is measured by the total number of symbol slots used in uploading local sketches, namely $tI_{n}$, where $I_{n}$ is the number of rows of local sketches. We vary the dimension of the receive beamformer, $M$, from $1$ to $4$ to achieve different tradeoffs between channel diversity and multiplexing. It is observed from Fig. 4 that the proposed selection scheme helps reduce the expected DTD error for different $M$. The gain emerges when the communication time exceeds $M$-dependent threshold (e.g. $6000$ time slots for $M=1$) as the total number of available sketches becomes sufficiently large. Another observation from Fig. 4 is that the DTD error without selection decreases at an approximately linear rate w.r.t. the communication time, which is aligned with our conclusion in (21). In the sequel, FlyCom2 is assumed to have sketch selection with $M$ fixed as $2$.

VI-C Error Performance of FlyCom2

While FlyCom2 requires much simpler on-device computation than benchmarking schemes (see Section VI-D), we demonstrate in Fig. 5 that it can achieve comparable or even better error performance than the latter. Fig. 5 shows the expected DTD error versus the communication time. The performance of the benchmark schemes with one-shot computation and communication appears as single points in the figure. The results in Fig. 5 show that FlyCom2-based DTD achieves comparable decomposition accuracy as the benchmark schemes with progressing time. Furthermore, its performance is improved by increasing the decay rate ($\xi$) of the singular values, which validates our conclusion that large eigen-gaps help distinguish principal from non-principal eigenvectors during random sketching. For instance, for $\xi=2$, the proposed scheme approaches the centroid and alignment based SVD-DTD in performance for communication time larger than $2600$ and $4000$ symbol slots, respectively. As $\xi$ increases to $3$, the former achieves the same error performance as the alignment-based method while outperforming the centroid based method. Furthermore, one can observe from Fig. 5 that the proposed on-the-fly framework realizes a flexible trade-off between the decomposition accuracy and communication time, which is the distinctive feature of the design.

VI-D Device Computation Costs of FlyCom2

In Fig. 6, we compare two kinds of computational costs at devices, namely complexity and memory passes, between FlyCom2 and benchmark schemes. The complexity refers to the flop count of computation, and the memory passes are equal to the number of memory visits for reading data entries. The computational advantage of FlyCom2 is demonstrated by comparing the cost of matrix-vector multiplication in random sketching with that of deterministic SVD used in the one-shot benchmarking schemes. Specifically, given $I\times J$ local unfolding matrices, deterministic SVD has the complexity proportional to $\min\{I,J\}^{2}\times\max\{I,J\}$ [18]; based on matrix multiplication, the complexity of FlyCom2 to yield an $M$-dimensional sketch at each time slot is $IJM$. For the schemes in comparison, their curves of computation complexity versus sample size are plotted in Fig. 6(a). One can observe that the proposed FlyCom2 dramatically reduces devices’ complexity by more than an order of magnitude. On the other hand, Fig. 6(b) displays the curves of the number of memory passes versus the principal dimensionality, $r$. The proposed design keeps a constant memory pass for matrix multiplication, as opposed to that of SVD which increases linearly with the principal dimensionality. For example, the number of memory passes is reduced using FlyCom2 by $30$ times for $r=30$.

VII Conclusion

We have proposed the FlyCom2 framework, that supports the progressive computation of DTD in mobile networks. Through the use of random sketching techniques at devices, the traditional one-shot high-dimensional mobile communication and computation is reduced to low-dimensional operations spread over multiple time slots. Thereby, the resource constraints of devices are overcome. Furthermore, FlyCom2 obtains its distinctive feature of progressive improvement of DTD accuracy with increasing communication time, providing robustness against link disruptions. To develop the FlyCom2 based DTD framework, we have designed an on-the-fly sub-space estimator and a sketch-selection scheme to ensure close-to-optimal system performance.

Beyond DTD, high-dimensional communication and computation pose a general challenge for machine learning and data analytics in wireless networks. We expect that FlyCom2 can be further developed into a broad approach for efficient deployment of relevant algorithms such as federated learning and distributed optimization. For the current FlyCom2 targeting DTD, its extension to accommodate other wireless techniques such as broadband transmission and radio resource management is also a direction worth pursuing.

-A Proof of Lemma 1

In (9), both the $\mathbf{\mathbf{F}}_{t}$ and $\tilde{\mathbf{Z}}_{t}$ have i.i.d. zero-mean Gaussian entries, thereby enforcing the observation $\tilde{\mathbf{Y}}_{t}$ to be a Gaussian matrix. This conclusion holds for all observations and thus the aggregation $\hat{\mathbf{Y}}_{t}$ is a Gaussian matrix and can be decomposed into $\mathbf{C}^{\frac{1}{2}}\mathbf{W}\mathbf{D}^{\frac{1}{2}}$. Therein, $\mathbf{W}$ has i.i.d. $\mathcal{N}(0,1)$ entries. The covariance matrices $\mathbf{C}$ and $\mathbf{D}$ are computed as follows. First, there is $\mathsf{E}[\hat{\mathbf{Y}}_{t}\hat{\mathbf{Y}}_{t}^{\top}]=\mathsf{E}[\mathbf{C}^{\frac{1}{2}}\mathbf{W}\mathbf{D}\mathbf{W}^{\top}\mathbf{C}^{\frac{1}{2}}]$, where the right hand side of the equation equals to $\mathbf{C}\mathsf{Tr}(\mathbf{D})$ while the left hand side is given as

[TABLE]

On the other hand, using $\mathsf{E}[\hat{\mathbf{Y}}_{t}^{\top}\hat{\mathbf{Y}}_{t}]=\mathsf{E}[\mathbf{D}^{\frac{1}{2}}\mathbf{W}^{\top}\mathbf{C}\mathbf{W}\mathbf{D}^{\frac{1}{2}}]=\mathbf{D}\mathsf{Tr}(\mathbf{C})$, we have

[TABLE]

where an arbitrary diagonal block, say $\mathsf{E}[\tilde{\mathbf{Y}}_{\ell}^{\top}\tilde{\mathbf{Y}}_{\ell}]$, $\forall\ell\leq t$, can be expressed as

[TABLE]

Concluding the above results yields the covariance matrices as $\mathbf{C}=\mathbf{X}\mathbf{X}^{\top}+\frac{1}{2tM}\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})\mathbf{I}_{I}$ and $\mathbf{D}=\frac{tM\mathsf{Tr}(\mathbf{X}^{\top}\mathbf{X})\mathbf{I}_{tM}+\frac{1}{2}ItM\sigma^{2}\mathsf{diag}(\mathbf{A}_{1}\mathbf{A}_{1}^{H},\cdots,\mathbf{A}_{t}\mathbf{A}_{t}^{H})}{tM\mathsf{Tr}(\mathbf{X}^{\top}\mathbf{X})+\frac{1}{2}I\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})}$, respectively. This completes the proof.

-B Proof of Lemma 2

First, rewrite the error as

[TABLE]

where the last step is due to $\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{i}\rangle^{2}=1-\sum_{j\neq i}\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle^{2}$.

Then, under the assumption of $\sigma_{i}=\sigma_{j}$, $\forall i,j\leq r$, the second term on the right side of the above equation can be rewritten as

[TABLE]

Putting the results above together, the conclusion in Lemma 2 follows.

-C Proof of Lemma 3

According to Remark 1, we have $\lambda_{i}>\lambda_{j}$, $\forall i\leq r<j$. Then, using the perturbation theory [29, Theorem 3.1 & Theorem 3.2], the $\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle$ can be upper bounded as

[TABLE]

where $\lambda_{i}$ and $\tilde{\lambda}_{i}$ denote the $i$-th eigenvalues of $\mathbf{\Lambda}$ and $\frac{1}{tM}\mathbf{\Phi}_{t}\mathbf{\Phi}_{t}^{\top}$, respectively. Based on $\min\{\max\{a,b\},c\}=\max\{\min\{a,c\},\min\{b,c\}\}\leq\max\{a,\min\{b,c\}\}$, the upper bound can be further written as

[TABLE]

Recall that $\lambda_{i}=\sigma_{i}^{2}+\sigma^{2}\sum_{\ell\leq t}\mathsf{Tr}(\mathbf{A}_{\ell}^{H}\mathbf{A}_{\ell})/2tM$, we have $\lambda_{i}-\lambda_{j}=\sigma_{i}^{2}-\sigma_{j}^{2}$, which yields

[TABLE]

which completes the proof.

-D Proof of Theorem 1

First, the square vector norm, $\|\mathbf{\Delta}\mathbf{u}_{j}\|_{2}^{2}$, can be rewritten as

[TABLE]

where we define $\mathbf{e}_{j}=[0,\cdots,0,1,0,\cdots,0]$ with the $j$-th element being $1$ and other elements being [math]. Since $\mathbf{W}$ has i.i.d. $\mathcal{N}(0,1)$ entries, the expectation of the upper bound of $(\sigma_{i}^{2}-\sigma_{j}^{2})\langle\tilde{\mathbf{u}}_{i},\mathbf{u}_{j}\rangle^{2}$ can be expressed as

[TABLE]

The result of $\mathbf{e}_{j}\mathsf{E}[\mathbf{W}\mathbf{W}^{\top}\mathbf{\Lambda}\mathbf{W}\mathbf{W}^{\top}]\mathbf{e}_{j}^{\top}$ can be derived by representing $\mathbf{W}=[\mathbf{w}_{1}^{\top},\mathbf{w}_{2}^{\top},\cdots,\mathbf{w}_{I}^{\top}]^{\top}$, where $\mathbf{w}_{i}$ has i.i.d. $\mathcal{N}(0,1)$ entries and is independent with $\mathbf{w}_{j}$ with $i\neq j$. In specific, there is

[TABLE]

Putting the above results together yields

[TABLE]

which completes the proof.

-E Proof of Theorem 2

Using Lemma 3, we rewrite the error bound as

[TABLE]

where the involved random part is defined as $g(\mathbf{W})\overset{\triangle}{=}\mathsf{Tr}\left(\mathbf{W}\mathbf{W}^{\top}\mathbf{\Lambda}\mathbf{W}\mathbf{W}^{\top}\mathbf{\Lambda}_{1}\right)-\mathsf{Tr}\left(\mathbf{W}\mathbf{W}^{\top}\mathbf{\Lambda}_{2}\right)$ with $\mathbf{\Lambda}_{1}=\sum_{i=1}^{r}\sum_{j\geq r+1}\frac{4}{\sigma_{i}^{2}-\sigma_{j}^{2}}\frac{\lambda_{j}}{(tM)^{2}}\mathbf{e}_{j}^{\top}\mathbf{e}_{j}$ and $\mathbf{\Lambda}_{2}=\sum_{i=1}^{r}\sum_{j\geq r+1}\frac{8\lambda_{j}^{2}}{(\sigma_{i}^{2}-\sigma_{j}^{2})tM}\mathbf{e}_{j}^{\top}\mathbf{e}_{j}$. It then follows that the probabilistic error bound is upper bounded as

[TABLE]

Next, rewrite $\mathbf{W}=[\tilde{\mathbf{w}}_{1},\tilde{\mathbf{w}}_{2},\cdots,\tilde{\mathbf{w}}_{tM}]$ and the random variable $g(\mathbf{W})$ can be rewritten as

[TABLE]

Let $\hat{\mathbf{w}}_{m}$ be an independent copy of $\tilde{\mathbf{w}}_{m}$ and then there is

[TABLE]

Note that the above equation does not have an upper bound since the Gaussian random variables go from minus infinity to infinity. To endow on $g(\mathbf{W})$ the bounded difference property, the Gaussian concentration is exploited as follows. Specifically, a $\mathcal{N}(0,1)$ variable $w$ can be smaller than a threshold, say $\kappa>1$, with the probability of $p(|w|\leq\kappa)=\mathsf{erf}\left(\frac{\kappa}{\sqrt{2}}\right)\overset{\triangle}{=}p_{\kappa}$, where $\mathsf{erf}\left(\cdot\right)$ denotes the error function. Hence, let the event that the abstract value of the last $I-r$ elements in the vectors $\tilde{\mathbf{w}}_{m}$ and $\hat{\mathbf{w}}_{m}$ are bounded by $\kappa$ be denoted by $\mathsf{BD}$ and its complement by $\mathsf{UBD}$. Then, there are $\mathsf{Pr}(\mathsf{BD})=p_{\kappa}^{tM(I-r)}$ and $\mathsf{Pr}(\mathsf{UBD})=1-p_{\kappa}^{tM(I-r)}$. Hence, with the probability of $\mathsf{Pr}(\mathsf{BD})$, $|g^{\prime}\left(\{\tilde{\mathbf{w}}_{M}\}_{M\neq m},\tilde{\mathbf{w}}_{m}\right)-g^{\prime}\left(\{\tilde{\mathbf{w}}_{M}\}_{M\neq m},\hat{\mathbf{w}}_{m}\right)|$ can be upper bounded as

[TABLE]

which allows us to leverage the concentration theorem shown in Lemma 4 to give

[TABLE]

Concluding both cases of $\mathsf{BD}$ and $\mathsf{UBD}$, the probabilistic error bound can be expressed as

[TABLE]

where the last inequality is due to the fact that a conditional probability is always smaller than $1$. Finally, with proper algebraic substitution, we have

[TABLE]

which completes the proof.

-F Proof of Lemma 5

First of all, for any $i\leq r$ and $j\geq r+1$, we have $\sigma_{i}^{2}-\sigma_{j}^{2}\geq\sigma_{r}^{2}-\sigma_{r+1}^{2}$, which gives

[TABLE]

Then, define $C=\sum_{k}\mathsf{Tr}\left(\mathbf{X}_{k}^{\top}\mathbf{X}_{k}\right)=\sum_{i=1}^{I}\sigma_{i}^{2}$ and $\eta_{\mathrm{th}}=\zeta_{\mathsf{th}}C$. It follows that $\sigma_{j}^{2}\leq\sigma_{i}^{2}\leq C/r$, $\forall j\geq r+1>i$, which further gives $r\lambda_{j}\leq r\sigma_{j}^{2}+\frac{r\sigma^{2}}{2}\zeta_{\mathsf{th}}C\leq(1+\frac{r\sigma^{2}}{2}\zeta_{\mathsf{th}})C$. As a result, one can obtain

[TABLE]

where the summation term can be upper bounded as

[TABLE]

Putting the above results together yields

[TABLE]

where $\frac{(I-r)C^{2}}{(\sigma_{r}^{2}-\sigma_{r+1}^{2})(I+1)r}$ can be treated as a constant independent from the variable $\zeta_{\mathsf{th}}$. This completes the proof.