I Introduction
Attaining the vision of InternetofThings (IoT) will require the ubiquitous deployment of an enormous number of sensors (e.g., tens of billions) in our society [1, 2]. The bruteforce approach of “transmitthencompute” is obviously impractical for this largescale sensor network as the massive radio access would result in excessive network latency and low efficiency in spectrum utilization. The situation is exacerbated at high mobility where ultrafast data aggregation from many sensors is desired. This is the case when sensors and/or the accesspoint (AP) are mounted on ground vehicles or unmanned aerial vehicles (UAV) for ubiquitous city surveillance in the smartcity application [see Fig. 1(a)], or for wildarea monitoring to avoid natural disasters [see Fig. 1(b)]. Motivated by the need of ultrafast data aggregation, an intelligent solution, known as overtheair computation (AirComp), is proposed recently that exploits the signalsuperposition property of a multiaccess channel (MAC) to compute a class of so called nomographic functions of distributed sensing data via concurrent sensor transmissions (see Fig. 1), thereby integrating computation and communication [3, 4]. Examples of such functions are show in Table I. Unlike ratecentric wireless systems where simultaneous transmissions result in interference, the computation accuracy for a sensor network with AirComp may grow with the number of simultaneous sensors due to the sensingnoise averaging. In this paper, we aim to advance the area of AirComp by developing multipleinputmultiple output (MIMO) AirComp for nextgeneration multiantenna multimodal sensor networks. The technology supports highmobility multimodal (HMM) sensing by enabling multifunction computation via spatial multiplexing and accurate reception of the results by exploiting spatial diversity.
Ia OvertheAir Computation for Sensor Networks
The idea of AirComp can be tracked back to the pioneer work studying functional computation in sensor networks [4].^{1}^{1}1AirComp also appeared the same time as a key operation in the scheme named physicallayer network coding in [5]. In [4], structured codes (e.g., lattice codes) are designed for reliably computing at an AP a function of distributed sensing values transmitted over a MAC. The significance of the work lies in its counterintuitive finding that “interference” can be harnessed to help computing. Subsequently, it was proved that the simple analog transmission without coding, where transmitted signals are scaled versions of sensing values, can achieve the minimum functional distortion achievable by any scheme [6]. On the other hand, coding can be still useful for other settings such as sensing correlated Gaussian sources [7]. The satisfactory performance (with optimality in certain cases) of simple analog AirComp have led to an active area focusing on designing and implementing techniques for receiving a desired function of concurrent signals, namely a targeted coherent combination of the signal waveforms [8, 9, 10, 11, 12, 13, 14, 3, 15, 16]. In particular, considering analog AirComp, power control at sensors was optimized in [8, 9], the computation rate (defined as the number of functional values computed per time slot) analyzed in [10]
, and the effect of channel estimation error characterized in
[11].The implementation of AirComp faces several practical issues. One is the synchronization of all active sensors required for coherent combining at the AP. To cope with synchronization errors, a scheme was proposed in [12, 13] where the sensing value is modulated at each sensor as the power of transmitted signal and furthermore a random phase rotation is applied to the signal. The design transforms functional computation at the receiver to power detection while synchronization error appears as random noise. An alternative solution, called AirShare, was developed in [14] for synchronizing sensors by broadcasting a referenceclock signal and its effectiveness was demonstrated using a prototype. In addition, by applying appropriate data pre/postprocessing, later work extended and implemented AirComp to compute a variety of functions besides the linear ones (as summarized in Table I) [15, 16].
It is worth mentioning that the coding techniques designed for AirComp in computingcentric sensor networks [4] inspired researchers to adopt relevant principles and ideas in designing new schemes for ratecentric communication networks [17, 18, 19, 20, 21]. For relay assisted networks, the application of AirComp at relay nodes for decoding and forwarding linear functions of the transmitted messages led to the invention of the well known computeandforward relaying schemes [17, 18, 19]. Building on lattice coding, a novel so called integerforcing linear receiver was proposed for spatial multiplexing in a multipleinputmultipleoutput (MIMO) system that attempts to create an effective channel matrix with integer coefficients to facilitate lattice decoding. The key operation, integer forcing, is similar to AirComp and computes a desired set of linear functions with integer coefficients [20, 21]. In parallel with the above research, extensive progresses were made in the area of physical layer network coding where the celebrated network coding schemes invented for wired networks were extended to wireless networks with AirComp relays (see survey in [22]).
Sensing devices targeting emerging applications such as smart cities are sophisticated. Each typically contains multiple multimodal sensors monitoring different environmental parameters (e.g., pressure, light, humidity, and temperature) [23]. In particular, a smartphone recruited for crowdsensing typically contains seven or more sensors such as inertial, GPS, and light sensors [24]. In view of prior work, the existing solutions focus only singlefunction AirComp assuming singleantenna sensors having unimodal sensing capabilities. However, nextgeneration wireless networks equipped with largescale arrays will make it possible to simultaneously compute multiple functions of multimodal sensing data overtheair. This inspires the current work on developing the technology of AirComp for a MIMO MAC, which can simultaneously spatial multiplex multifunction computation and suppress computation errors by exploiting spatialdiversity gain. Thereby, the datafusion latency in sensor networks can be substantially reduced, meeting the ultralow latency requirement in nextgeneration networks especially when highmobility support is needed [1].
IB MIMO Beamforming: MultiAccess versus AirComp
Beamforming design for multiuser MIMO systems is a classic topic that has been extensively studied and there exists a rich relevant literature [25, 26, 27]. In terms of network topology, the multiantenna multimodal sensor network we consider is equivalent to a MIMO multiaccess communication network where a single AP supports simultaneous uplink transmissions of multiple users. For the communication network, the designs of multiuser MIMO beamforming at the AP can be largely grouped into capacityachieving nonlinear designs based on successive interference cancellation (SIC) [25] and lowcomplexity linear designs (minimum meansquared error (MMSE) or zeroforcing) [26, 27, 28]. All of the designs share the same objective of decoupling multiuser signals by interference suppression and spatial multiplexing of data streams for each user. In contrast, AirComp receive beamforming in the sensor network has a different objective of minimizing the total distortion in the received values of multiple functions combining multimodal data simultaneously sent by a set of sensors. Due to the difference in objective between communication and sensing, the known designs for the former are inapplicable for the latter. On the other hand, existing AirComp literature considers only unifunction computation targeting unimodal sensing as discussed earlier. This makes receive beamforming for multifunction AirComp for multimodal sensing an uncharted problem to be tackled in this work.
It is worth mentioning that the discovery of uplinkdownlink duality is a breakthrough in multiuser MIMO communication. The duality reveals similar structures in optimal beamformers for the MACs and broadcast channels that exist under various performance criteria ranging from capacity maximization [29, 30] to MMSE [31]. This allows beamforming designs derived for the MACs to be applied to the broadcast channels where the optimal beamforming design was largely an open problem prior to the finding of the duality. Inspired by this finding, we address a similar question in the current work: What is the downlink dual of the (uplink) AirComp beamforming for sensor networks?
Name  Expression 

Arithmetic Mean  
Weighted Sum  
Geometric Mean  
Polynomial  
Euclidean Norm 
IC Contributions and Organization
We consider a multiantenna multimodal sensor network where a multiantenna AP performing fusion of data transmitted by a cluster of multiantenna multimodal sensors. By measuring multiple timevarying parameters of the environment, each sensor generates multiple data streams. In each time slot, a sensor transmits a set of multimodal data values in the analog way, namely by amplitude modulation [12, 13], over multiple antennas. The transmissions of all sensors are simultaneous. The AP attempts to jointly receive multiple nomographic functions (such as those in Table I) of distributed sensor data by AirComp and spatial multiplexing. The AirComp of a nomographic function is implemented by three cascaded operations: 1) preprocessing at sensors, 2) weighted summation of preprocessed outputs by simultaneous transmissions, and 3) postprocessing at the AP [12, 13, 3]. In the current scenario, transmit and receive beamforming are applied to spatially multiplex multifunction AirComp as well as exploit spatial diversity to minimize the distortion of function values caused by noise, which is measured by sum meansquared error (MSE) over functions.
While the traditional unifunction AirComp is a simple technique, the proposed multifunction version is challenging with the optimization of receive beamforming proved to be NPhard. Specifically, for unifunction AirComp, channel inversion at each sensor yields a desired weighted sum of preprocessed data at the AP, giving the desired function value after postprocessing [3]. For multifunction AirComp, channel inversion remains optimal as shown in this work and is implemented by zeroforcing beamforming. Nevertheless, receive beamforming for multifunction AirComp, referred to as MIMOAirComp equalization, is a new design problem that finds no relevant result in the AirComp literature. The equalizer optimization is nonconvex but can be relaxed as a semidefinite programming (SDP) problem and thus solved using an iterative interior point algorithm. The standard approach does not yield any insight into the optimal equalizer structure and more importantly does not lead to an efficient channelfeedback design for acquiring the equalizer at the AP. To address these issues, we impose an orthogonality constraint on the AirComp equalizer, which is a technique for approximate beamformer design and limited channel feedback as widely applied in the literature [32, 33, 34, 35]. Concretely, this allows a closetooptimal equalizer to be derived in closedform using tools from differentialgeometry, revealing an interesting geometry structure in the design. Moreover, the closedform solution leads to an efficient channel feedback design that exploits the AirComp architecture for direct equalizer acquisition at the AP from simultaneous feedback by all sensors.
The main contributions of this work are summarized as follows.

MultiFunction AirComp Beamforming: As mentioned, while zeroforcing transmit beamforming is found to be optimal, the receivebeamformer optimization for sumMSE minimization under transmissionpower constraints can be proved to be NPhard. By tightening the constraints, the resultant approximate problem is found to involve optimization on a Grassmann manifold, which can be interpreted as the space of subspaces. This allows the application of differential geometry to solve the approximate problem. The derived solution shows the normalized receive beamformer to be the weighted centroid of a cluster of points on the manifold
, where each point represents the eigensubspace of an individual MIMO channel and the corresponding weight its smallest eigenvalue. In addition, the optimal beamformer norm is also derived in closedform. Such a beamformer design allowing efficient computation is verified by simulation to be closetooptimal.

AirCompMulticasting Duality: As a byproduct of our investigation, for the special case of AirComp with singleantenna (unimodal) sensors, the problem of receivebeamforming optimization is discovered to have the identical form as the classic problem of multicast transmitbeamforming, thereby establishing a novel AirCompmulticasting duality. The latter problem is known to be NPhard and typically solved using the semidefinite relaxation (SDR) method. The significance of our finding lies in allowing the solution method for AirComp beamforming to be transferable to multicastbeamforming. This yields a new solution method for the latter with complexity much lower than the existing SDR approach as the network scales up.

AirComp Channel Feedback: Last, we solve the AirComp feedback problem: How to efficiently acquire the derived AirComp beamformer at the AP, which depends on global channelstate information (CSI), by sensor distributed transmissions based on local CSI? Given channel reciprocity, it is discovered that the AirComp system architecture can be also used for efficient feedback. The resultant number of feedback rounds is independent of the sensor population, overcoming the drawback of traditional channel training. Novel feedback techniques are designed for sequential feedback of the normalized AirComp beamformer and beamformer norm based on their derived expressions, where each feedback round involves concurrent transmissions by all sensors. Essentially, the two techniques implement AirComp of two functions, namely the weighted centroid of a set of matrices and the maximum of a set of scalar values. They are hence exclusively for multifunction AirComp and may not be applicable for traditional multiuser MIMO communication systems where multiuser feedback requires orthogonal channels and focuses on precoder quantization [36, 37].
The remainder of the paper is organized as follows. Section II introduces the AirComp system model. Section III presents the problem formulation for the enabling beamforming design and channel feedback. The proposed beamforming design is presented in Section IV, and the duality between the uplink AirComp and downlink multicasting is discussed in Section V. Section VI proposes an efficient channelfeedback scheme that can be implemented by AirComp. Simulation results are provided in Section VII, followed by concluding remarks in Section VIII.
Ii System Model
We consider a wireless sensor network consisting of multimodal sensors and a single AP as illustrated in Fig.2. All nodes are equipped with antenna arrays. Specifically, antennas are deployed at each sensor and at the AP. Each multimodal sensor records the values of heterogeneous timevarying parameters of the environment, e.g., temperature, pollution, and humidity. The data from the record of the th parameter is referred to as type
data. For an arbitrary time slot, the measurement vector of the
th sensor, grouping sample values, is denoted by , where is the measurement of the parameter at sensor . Instead of collecting the whole data set, the AP aims at computing functions of corresponding types of data, denoted by , to support ultrafast HMM sensing. The class of functions that are computable by AirComp are called nomographic functions as defined below.Definition 1 (Nomographic Function [13]).
The function is said to be nomographic, if there exist preprocessing functions along with a postprocessing function such that it can be represented in the form:
(1) 
Some common nomographic functions are listed in Table I. Based on the nested form of (1), the AirComp of a nomographic function can be implemented in the sensor network by three operations as illustrated in Fig. 2: 1) preprocessing at each sensor specified by where operates on the type data at sensor , 2) summation of preprocessed data realized by multiaccess, and 3) postprocessing at the AP. Considering the computation of the geometric mean of type data as an example, the preprocessing , and the postprocessing . Let denote the multimodal symbol vector after preprocessing and
the AirComputed function values. For ease of transmissionpower control and without loss of generality, the symbols are assumed to be normalized to have unit variance, i.e.,
, where the normalization factor for each data type is uniform for all sensors and can be inverted at the AP. Given the onetoone mapping between and according to (1), we refer to as the targetfunction vector for ease of exposition.Iia AirComp Phase
Assuming symbollevel synchronization^{2}^{2}2This can be achieved by broadcasting a common reference clock from the AP to sensors using the AirShare solution developed in [14]., all users transmit their symbol vectors simultaneously using their arrays. The distortion of the received vector with respect to the targetfunction vector due to MIMO channels and noise is suppressed using transmit and receive beamforming. In other words, the joint beamforming attempts to attain coherent combining of symbol vectors at the AP. Let denote the receive beamforming matrix and the transmit beamforming matrix at sensor . Then the symbol vector received by the AP after receive beamforming is given as
(2) 
where represents the MIMO channel matrix for the link from the sensor to the AP, and is the additive white Gaussian noise (AWGN) vector with independent and identically distributed (i.i.d.) elements. The distortion of with respect to , which quantifies the AirComp performance, is measured by the MSE defined as follows:
(3) 
IiB Channel Feedback Phase
Consider the existence of channel reciprocity and assume that perfect local CSI is available at all sensors. One can infer from (4) that computing the MMSE receive beamformer requires global CSI, namely . As mentioned, the naive approach of estimating the global CSI would incur long latency and large overhead when the number of sensors is large, thus is impractical for the ultrafast HMM sensing applications. An intelligent channelfeedback design is proposed in the sequel to allow beamformer acquisition via concurrent transmissions by all sensors. Let denote the signal matrix transmitted by sensor where is the signal length in symbol. Given typical high transmission power for channel training and feedback, the feedback observation at the AP can be assumed to be noiseless and thus be represented by an matrix as follows:
(5) 
As is clear in the sequel, with proper design of , can serve as a sufficient surrogate of the global CSI in beamformer computation at the AP.
Iii Problem Formulation
Iiia AirComp Beamforming Problem
Consider the joint optimization of the transmit and receive beamformers under the MMSE criterion and the transmissionpower constraints. It is assumed that the average transmission power of each sensor cannot exceed a given positive value . Since the transmitted symbols have unit variance, the power constraints are given as
(6) 
Following a common approach in the MIMO beamforming literature (see e.g., [32, 33, 34, 35]), the receive beamformer is constrained to be orthonormal matrix. As mentioned, the constraint can lead to a closedform suboptimal solution with only marginal performance loss and furthermore facilitate efficient channel feedback design as presented in Section VI. As pointed out in [35], for most communication objectives, it is the subspace spanned by the beamformer but not the exact beamformer that has a crucial effect on the system performance, justifying the said constraint. Furthermore, under the MMSE criterion, a positive scaling factor , called denoising factor, is included in for regulating the tradeoff between noise reduction and transmissionpower control. To be specific, reducing suppresses noise but requires larger transmission power to maintain the MSE of computed function values. Mathematically, we can write with
being a tall unitary matrix and thus
. Then given the MSE in (4), the MMSE beamforming problem can be formulated as:The problem is solved in Section IV.
IiiB AirComp Channel Feedback Problem
We propose the use of the AirComp architecture in Fig. 2 to realize the mentioned receivebeamformer acquisition by concurrent transmissions by all sensors. Let denote the derived beamformer solution to problem P1 and and be the feedback counterparts of the AirComp operations and (see Fig. 2). The key design constraint is that the transmitted signal in (5) must be a function of local CSI only, denoted as . Then it follows that
and the problem of AirComp feedback design reduces to the design of the functions and . The solution is presented in Section VI.
Iv MultiFunction AirComp: Beamforming
In this section, the AirComp beamforming problem in Problem P1 is solved. While zeroforcing transmit beamforming can be proved to be optimal, the receive beamforming optimization is found to be NPhard. An approximate problem is obtained by tightening the power constraints. This problem allows a practical solution approach based on differential geometry. The solution reveals that the optimal receive beamformer can be approximated by the weighted centroids of a cluster of points on a Grassmann manifold, each corresponding to the subspace of an individual MIMO channel. To facilitate exposition, some basic definitions and principles of Grassmann manifolds are provided in Appendix A.
Problem P1 is difficult to solve due to its nonconvexity. The lack of convexity arises from the coupling between transmit and receive beamformers in the objective function, and the orthogonality constraint on the receive beamformer. To simplify the problem, zeroforcing transmit beamforming conditioned on a receive beamformer is first shown to be optimal as follows.
Lemma 1.
Given a receive beamformer , the MSE objective stated in (4) is minimized by the following zeroforcing transmit beamformers:
(7) 
Proof: See Appendix AA.
We note that the power constraint imposed on the precoder will be enforced in the sequel via regulating the norm of the equalizer , or equivalently, the denoising factor .
Remark 1 (Number of AirComputable Functions).
Note that to ensure matrix is invertible, it requires . This implies that, the number of functions that can be simultaneously computed by the proposed multifunction AirComp is limited by . The result is due to the underpinning limit of MIMO spatial multiplexing: the maximum number of spatial streams is .
By substituting (7) in Lemma 1, Problem P1 is transformed to the equivalent problem of minimizing the denoising factor of the receive beamformer:
where is defined earlier as the normalized receive beamformer. Though Problem P3 has a simpler form than P1, it remains nonconvex due to the nonconvex orthonormal constraint of the receive beamformer. In fact, Problem P3 is found in the next section to be NPhard via proving its equivalence to the NPhard multicast beamforming problem. To develop a tractable approximation of the problem, a reasonable modification of the power constraints is derived. To this end, a useful inequality is obtained as shown below.
Lemma 2.
Let denote the compact form of singular value decomposition (SVD) of . Then we have the following inequality:
(8) 
where the equality holds given a wellconditioned channel, i.e., for some constant . Proof: See Appendix AB.
Tightening the power constraints in Problem P3 using Lemma 2 gives the approximate problem:
Since the feasible set of Problem P4 is smaller than that of P3, the solution to P4 is a feasible solution though potentially a suboptimal one to P3. To solve Problem P4 using differential geometry, an equivalent form containing subspace distances between the receive beamformer and individual MIMO channels is obtained as follows.
Lemma 3.
Problem P5 is not yet in a ready form admitting the differentialgeometry solution approach and requires an additional approximation. For this purpose, the objective function is bounded below.
Lemma 4.
The objective function in Problem P5 can be bounded as follows:
(9) 
where we define which is a constant independent of the control variable .
The proof is straightforward and omitted for brevity. Approximating the objective function in P5 by either the lower or the upper bound in Lemma 4 both lead to the same approximate problem which is given by:
Remark 2 (Beamformer Geometric Interpretation).
Problem P6 allows a geometric interpretation of the desired receive beamformer . In fact, the problem is to find a weighted centroid of a set of points (each being a subspace) on a Grassmann manifold with the squared projection 2norm as the distance metric, where the weights are . This reveals that the receive beamformer makes the besteffort to be aligned with all the MIMO channel matrices with the alignment (subspace) distances adjusted by corresponding channel gains as specified by the smallest channel eigenvalues.
Problem P6 can be approximately solved by a closedform solution that can be efficiently computed without resorting to an iterative algorithm. Particularly, the closedform solution can be derived by replacing the projection 2norm distance with the projection Fnorm (see Appendix A). Note that for a small principal angle [38] and are exactly equivalent in the case of according to (27). Thus, the problem P6 can be approximated as
(10) 
which still seeks a weighted centroid of channel subspaces as before but based on a different subspace distance metric. According to the definition in (26), can be computed in a matrix form by
Then, substituting it to the objective function in Problem P7, the problem can be equivalently written as
(11) 
Problem P7 remains nonconvex due to 1) maximimizing a convex objective function and 2) the orthogonality constraints on the variable . Nevertheless, by intelligently constructing an equivalent unconstrained problem, we are able to derive a closedform solution for Problem P7 (see the following Lemma 5) via analyzing the stationary points of the unconstrained problem. For ease of exposition, define a matrix , called effective CSI, as follows:
(12) 
As shown shortly in Lemma 5, the normalized AirComp receive beamformer depends on the global CSI via the effective CSI. In other words, is sufficient for computing .
Lemma 5.
V AirCompMulticasting Duality
In this section, consider the case of singleantenna unimodal sensors. The AirComp receivebeamforming problem for uplink sensingdata collection is shown to be equivalent to the well known problem of transmit beamforming for downlink multicasting (see Fig. 3). This establishes the AirCompmulticasting duality, allowing the lowcomplexity beamforming design in the preceding section to be transferable to solve the NPhard multicast beamforming problem.
Va Review of the Multicast Beamforming Problem
Consider the scenario that multiple singleantenna users request the same data stream from a multiantenna AP equipped as shown in Fig. 3(a). Assuming global CSI is available at AP, the problem of multicast beamforming is to minimize the total transmission power subject to a set of signaltonoiseratio (SNR) constraints specifying the users’ qualityofservice requirements. Mathematically, the problem can be formulated as follows:
(15) 
where denotes the multicast beamforming vector, , and are the channel vector, noise variance and target SNR of user , respectively. The problem can be proved to be NPhard [39]. Nevertheless, it is known that a closetooptimal solution can be efficiently computed using the well known SDR technique. The key idea of SDR is to recast the problem as an equivalent rankone constrained SDP by denoting as follows:
(16) 
Then SDR drops the rankone constraint and solves the relaxed SDP. Finally, a rankone approximate solution of the original problem is retrieved by a Gaussian randomization strategy based on the solution of the relaxed SDP (or simply the principal eigenvector of it). For more details on the SDR algorithms, readers are referred to the key references in the area
[39, 40].VB Duality between AirComp and Multicast Beamforming
For the special case of a sensor network with singleantenna sensors () and a multiantenna AP. The receive beamforming reduces to a vector denoted by . Its design problem for MMSE AirComp can be directly simplified from problem P5 (which is equivalent to the original P3 when ) to the following form:
(17) 
where represents the orientation of the channel vector, i.e., .
Lemma 6.
Remark 3 (AirCompMulticasting Duality).
Comparing (15) and (18) reveals that the beamforming problems for the uplink AirComp and downlink multicasting share the same mathematical form. This establishes the AirCompmulticasting duality that is analogous to the famous uplinkdownlink duality for multiuser MIMO communication [29]. Intuitively, the AirCompmulticasting duality is a result of the fact that both the AirComp and multicast beamformers must make the same best effort to be aligned with multiple vector channels (see Fig. 3) though for different objectives: one is to minimize the distortion in the computed function value and the other maximize the minimum SNR.
The AirComp beamforming technique designed in the last section is based on computing the weighted centroid on the Grassmann manifold. The duality allows the technique to be applied to multicast beamforming. Compared with the classic SDR method discussed in the preceding subsection, the AirComp beamforming has the following two main advantages.

(Efficient CSI Feedback): The AirComp beamforming solution requires only the effective CSI in (12) and thus enables the efficient “oneshot” channel feedback scheme presented in the next section. In contrast, the SDR solution requires global CSI and thus requires all users to feed back their local CSI. This results in excessive channel training overhead when the number of users is large.

(Low Computation Complexity): As shown in Lemma 5, the weighted centroid solution requires only oneshot computation and has a relatively low complexity of arising from the principal eigenvector calculation. However, the SDR requires first solving a SDP of dimension , by an iterative interior point method, resulting in the complexity of where denotes the solution accuracy [40]. The complexity becomes overwhelming when the numbers of the receive antennas and/or users are large. This makes the AirComp solution preferable in practice. The lowcomplexity advantage of the proposed solution is also verified by simulation in Section VII.
Vi MultiFunction AirComp: Channel Feedback
In this section, the AirComp feedback problem stated in Problem P2 is solved by feedback technique design. Using the optimization results in Section IV, two novel techniques are designed in the following subsections for sequential feedback of two components of the AirComp receive beamformer, namely the normalized beamformer and the denoising factor. Essentially, the two techniques realize the AirComp of two functions: 1) the weighted centroid of a set of matrices and 2) the maximum of a set of scalars, which can be thus implemented using the AirComp architecture in Fig. 2.
Via Feedback of Normalized Beamformer
Assume that the feedback channel is noiseless due to high transmission power for channel feedback. The AirComp feedback scheme for acquisition of the normalized beamformer in (13) is derived as follows. As indicated by Lemma 5, the normalized beamformer can be directly computed from effective CSI matrix in (12). The key step for solving the AirComp feedback problem in Problem P2 is to enforce the equality . Then given and Let denote the compact SVD of , it is easy to verify using Problem P2 that designing the feedback signals as gives the desired equality . Then the normalized receive beamformer can be computed as the principal eigenvectors of the received signal . The AirComp feedback design for acquiring is summarized as follows.
Normalized Beamformer Feedback:
(19) 
where denote the dominant left eigenvectors of . It is important to note that scaling the feedback signal by a constant so as to meet a transmission power constraint has no effect on the received normalized beamformer .
The above AirComp feedback technique inherits the advantage of AirComp by turning interference from multiple access into useful signals for functional computing. In contrast with the traditional method of channel training, increasing the number of simultaneous sensors may even be beneficial to the computation accuracy via sensingnoise averaging.
ViB Feedback of Denoising Factor
Following (14), the denoising factor is the maximum of a set of scalars called feedback values and defined as:
(20) 
Since the maximum is not a nomographic function, it is not directly AirComputable. However, an intelligent feedback technique is presented shortly that shows the possibility of denoisingfactor acquisition by AirComp over a fixed number of feedback rounds. To begin with, it is assumed that the normalized beamformer is acquired at the AP using the technique in the preceding subsection and then broadcasts to all sensors. This allows each sensor to apply zeroforcing transmit beamforming for inverting the corresponding channel matrix. Specifically, the transmit beamformer at sensor is given as . Such beamforming creates a effective set of parallel MACs such that the signal vectors transmitted by sensors are summed at the AP. In other words, with denoting the signal vector for sensor , the receive signal vector at the AP is . Furthermore, we assume that the denoising factors lie in a fixed finite range . Given the parallel MACs and the assumption, the algorithm for denoisingfactor feedback with feedback rounds is described as follows.
Algorithm for DenoisingFactor Feedback:

(Initialization): Set the feedback counter and intialize the feedbackquantization range with , and .

(Feedback Quantization): A quantizer codebook with values, denoted by , with , is generated by uniformly partitioning the range . Thus, the maximum quantization error is bounded by half of each partition interval denoted as : . Quantizing the feedback value in (20) at sensor gives the codebook index .

(Concurrent Feedback): Each sensor transmits a signal vector comprising a single at the location specified by the corresponding codebook index and ’s at other locations. Specifically, the signal vector for sensor is
(21) Then all sensors transmit their signal vectors simultaneously over the said effective parallel MACs. The AP finds the largest index of a nonzero element in the received signal vector , denoted as the . Thereby, using the codebook , it can be inferred at the AP that the quantized value of the denoising factor is and the exact value lies in the range .

(Refining Quantization Range): To improve the quantization resolution in the next feedback round, the AP refines the quantization range as
Next, increase the counter by setting and go back to 2) if or otherwise stop the feedback process.
From the above algorithm, the key result on the feedback accuracy follows.
Proposition 1.
Given feedback rounds, the AP receives a quantized version of , denoted as , with the quantization error bounded as
(22) 
Proposition 1 implies that the feedback error reduces exponentially with the number of feedback rounds. Specifically, adding a feedback round improves the feedback resolution by bit. Given a target resolution with the maximum quantization error , the required number of feedback rounds by the proposed feedback scheme is given by
(23) 
As an example, for some practical settings of , , , , the required number of rounds (feedback slots) is according to (23). This is much smaller than the number of sensors in a dense network which determines the feedback rounds if the traditional method of channel training is adopted.
Remark 4 (Comparison with StateoftheArt).
The stateoftheart algorithm for AirComp of a maximum function was proposed in [3] based on a different principle from the current design. In [3], the maximum of a set of distributed feedback values is progressively computed at the AP by sequential detection of the bits in the binary representation of the desired value via scheduling transmitting sensors by broadcasting a threshold. As the result, each feedback round increases the feedback resolution by a single bit and the algorithm cannot be straighforwardly extended to exploit spatial multiplexing. In contrast, by exploiting spatial channels for implementing uniform quantization, the proposed feedback algorithm achieves multibit resolution improvement for each feedback round as mentioned earlier.
ViC Comparison with Conventional Channel Training
For conventional multiuser channel training, sensors take turns to transmit pilot signals to AP for uplink channel training to avoid collision (see e.g., [26, 27]). The pilot signal for each sensor should be a or larger matrix for estimating a channel matrix. Thus it takes at least symbol slots to complete the channel training process for a network comprising sensors. In contrast, the proposed AirComp feedback technique for normalized beamformer feedback involves simultaneous transmissions of all , each of size , which thus requires only symbol slots. This together with the slots (typically ) for the feedback of the denoising factor yields the total feedback slots of independent of the network size . Consider a typical dense sensor network with , and , it takes only slots for the proposed AirComp feedback scheme in contrast to slots required by the conventional channel training. Thus AirComp feedback achieves time of feedback overhead reduction in this example.
Vii Simulation Results
In this section, the performance of the proposed multifunction AirComp is evaluated by simulation. The simulation parameters are set as follows unless specified otherwise. The number of multimodal sensors is , the AP array size at AP , the sensor array size and the number of computed functions are equal to . Each MIMO channel are assumed to be i.i.d. Rician fading
, modelled as i.i.d. complex Gaussian random variables with nonzero mean
and variance . In addition, the average transmitSNR constraint, defined as , is set to be dB.ViiA Baseline Beamforming Schemes
Given that the optimization of AirComp beamforming is a NPhard problem, for the purpose of comparison, we consider two baseline AirComp beamforming schemes designed based on the classic approaches, namely antenna selection and eigenmode beamforming. Both schemes assume zeroforcing transmit beamforming in (7) and their difference lies in the receive beamformers. Define the sumchannel matrix . To enhance the receive SNRs, the antennaselection scheme selects the receive antennas observing the largest channel gains in the sum channel . Consequently, the effective channel matrix after beamforming consists of rows of with largest vector norms. On the other hand, to select the strongest eigenmodes of for AirComp, the normalized eigenmode receive beamformer consists of the dominant left eigenvectors of . The denoising factor of each type of beamforming design is computed following (14) with modified accordingly.
ViiB Performance of MultipleFunction AirComp
In Fig 4, the MSE performance of the proposed multifunction AirComp beamforming is compared with that of two baseline schemes introduced in the preceding subsections. A varying number of functions , size of receive array , and also number of sensors are considered in Fig. 4(a)  4(c), respectively. Several key observations can be made as follows. First, for all schemes, the MSE is a increasing function of and but an increasing function of . This coincides with our intuition that, higher computation throughput is at a cost of declining accuracy, and more connected sensors makes it harder to design one common receive beamformer to equalize all different users’ MIMO channels. Nevertheless, deploying more receive antennas compensate for the performance degradation by exploiting diversity gain. Second, under various parameter settings, the proposed scheme outperforms the other two baseline schemes, showing the effectiveness of the new design approach based on optimization on the Grassmann manifold. Furthermore, the performance gain of the proposed design is larger in the regime of large and , further confirming the effectiveness of the proposed design for multimodal sensing and dense networks. Last, one can observe that the performance between different schemes converge as grows. This suggests that the large diversity gain enhances the receive SNRs such that the optimzation of AirComp beamforming is less critical and simple designs suffice.
ViiC Comparison with the SDR Method
The discovered AirCompmulticasting duality leads to the availability of two methods, the proposed weighted centroid and the SDR methods, for designing beamforming in either type of systems. Their performance and complexity are compared by simulation as follows. Consider singleantenna unimodal sensors as in Section V. The comparison of the MSE performance and computation time between the proposed weighted centroid and the SDR solutions is provided in Fig. 5 for the varying receive array size and number of sensors
. The computation time is measured using MATLAB. It is observed that the weighted centroid solution can achieve comparable performance as the SDR solution, which is optimal with a high probability for the NPhard multicast beamforming problem as shown in
[40]. On the other hand, the former achieves dramatic computation time reduction with respect to the latter, ranging from 100x to 1000x in the considered ranges of and . The simulation results support our previous analysis in Section V that the complexity of the proposed solution is independent of the network size and furthermore insensitive to the variation of the array size while the complexity of the SDR solution is . Thereby, the proposed solution features low complexity and is preferred in the large scale sensor networks (or large scale multicast networks) or when the AP is equipped with a large scale array.Viii Concluding Remarks
In this paper, we have proposed the framework of multifunction AirComp for MIMO HMM sensor networks. In particular, we have developed an approach for designing receive beamforming using tools from differential geometry. This approach achieves dramatic complexity reduction than the stateoftheart SDR approach while maintaining comparable performance. Furthermore, building on the AirComp system architecture, intelligent channel feedback techniques have been designed for enabling AirComp beamforming. Unlike the traditional method of channel training, the techniques prevent feedback overhead from escalating with the number of sensors and thus are highly efficient for dense HMM sensor networks. Last, the discovery of AirCompmulticasting duality allows the lowcomplexity beamforming design to be transferable to multiantenna multicast systems, which traditionally relies on the computationintensive SDR method for beamforming optimzation. The work points to the promising new research area of MIMO AirComp where many interesting research issues warrant further investigation such as sensor scheduling, broadband AirComp, AirComp for multiAP cooperative sensor networks, and AirComp for supporting distributed learning and inference.
A Preliminaries on Grassmann Manifold
A1 Stiefel and Grassmann Manifolds
The Stiefel manifold is the set of all by tall orthonormal matrices for , denoted by . Mathematically, . On the other hand, the Grassmann manifold is a set of all dimensional subspaces in , denoted by . Thereby a Grassmann manifold can be seen as the quotient space of . To be specific, a point on the Grassmann manifold corresponds to a class of by orthonormal matrices on the Stiefel manifold that span the same column subspace defined by the point. Choosing an arbitrary matrix from this class and using it as a generator, the class, denoted as , can be mathematically written as where denotes the group of unitary matrices. This leads to a relation between the Grassmannian and the Stiefel : [38].
A2 Distance Metrics on Grassmann manifold
Algorithms on Grassmann manifold often involves the calculation of the distance between points on the manifold. There exist many different distance metrics and all of them are derived from a key notion called geodesic. Roughly speaking, a geodesic is the unique curve linking two points on a manifold that has the shortest length among all. The length of the geodesic is called the geodesic distance (or arc length). Mathematically, given , their geodesic distance, denoted as is calculated by where are called the principal angles, measuring the minimal angles among any two sets of orthonormal bases spanning the two subspaces[38]. An efficient way to compute the principal angles is to perform SVD on , i.e.,
(24) 
where the singular values in (24) yields the cosines of the principal angles. Based on these angles, a rich set of subspacedistance metrics can be defined. Two particular metrics of relevance in this paper are [38]:
(25)  
(26) 
where is the vector formed by , and matrices and follow those in (24). For , we have the following inequalities relating different distance metrics:
(27) 
Appendix A Supplementary Proof for the Derived Key Results
AA Proof of Lemma 1
Given the MSE objective provided in (4), it is easy to note that both the first and the second terms within, i.e., and are positive semidefinite matrix with nonnegative eigenvalues. As a result, for any given equalizer , we have the following inequality:
(28) 
It is easy to verify that setting to have the zeroforcing structure in (7) enforces
and thus achieves the equality in (28), which completes the proof.
AB Proof of Lemma 2
Utilizing the fact that, for any square matrix, the inequality holds, it is straightforward to show
(29) 
Note that is a tall unitary matrix, thus matrix and share the same eigenspectrum due to the wellknown unitary invariant property. Therefore, it is easy to see that the upper bound (29) becomes exact when the channel is wellconditioned.
Then given the compact eigenvalue decomposition , we have the following
(30)  
(31)  
(32) 
where the second equality uses the fact that, for arbitrary two matrices and of the same size, and have the same eigenspectrum, while the last inequality is due to [41, Corollary 11], and it is easy to verify that the equality holds when
is a scaled identity matrix, namely
which also implies the channel is wellconditioned.AC Proof of Lemma 3
Note that the set of power constraints in P4 can be rewritten as one single constraint by:
(33) 
It is easy to note that the minimum in P4 is achieved when the above constraint is active (the equality holds). Therefore, one can move the constraint in (33) to the objective function and have the following equivalent minmax problem:
(34) 
The problem in (34) can be further simplified by maxmin the inverse of the objective function and dropping the constant term which leads to the following form
(35) 
The objective function in (35) is related to the projection 2norm Grassmannian metric via
(36) 
The equation can be easily verified by using the definitions provided in (24) and (25).
AD Proof of Lemma 5
In this proof, we first construct an equivalent problem of P7 by modifying the objective function leveraging the orthogonality constraint, giving the problem . Then, we show that solving a relaxed version of without the orthogonality constraint, denoted by , yields a solution enforcing the constraint still. Therefore, one can conclude that problem P7 and share the same optimal solution, while the unconstrained problem can be solved easily by checking all stationary points of the objective. The detailed derivation is presented below.
After some simple algebra manipulation exploiting the linearity of the trace operation, the objective function in (P7) can be further simplified as
(37) 
where is the effective CSI that has been defined in (12). Then problem (P7) reduces to
(38) 
Starting with problem , let’s first define an alternative objective function given by
(39) 
It is straightforward to verify that the following problem is equivalent to problem .
(40) 
Now, let’s relax the orthogonality constraint, and consider the following unconstrained problem:
(41) 
Since the function is a smooth function with gradient defined everywhere, the solution to the problem should be a stationary point of , i.e., . It follows that
(42) 
From (42), one can note that is one of the stationary points.
To seek other stationary points that , we left multiply both sides of the equality in (42) with