Joint Optimization of Secure Over-the-Air Computation and Reliable Multicasting Assisted by a MIMO Untrusted Two-Way Relay

In this article, we investigate the joint optimization of secure Over-the-Air Computation (AirComp) and reliable multicasting assisted by a multiple-input–multiple-output-untrusted two-way relay, where artificial noise is employed at the access point (AP) to interfere the relay for ensuring secure AirComp. We aim to minimize the computation distortion at the AP by jointly designing the transmit variables of all the nodes and the aggregation variables at the AP and relay, under the secure AirComp constraint, the reliable multicasting constraint, and the transmit power constraints of all the nodes. We consider two scenarios that perfect and imperfect channel state information (CSI) are available. For the former, the formulated optimization problem is highly nonconvex, and we propose an effective block coordinate descent (BCD)-penalty successive convex approximation (penaltySCA) method to solve the nonconvex problem. For the latter, we model the imperfect CSI by using the worst case criterion and the formulated robust optimization problem is much more challenging than its counterpart with perfect CSI. To solve the robust problem effectively, we first transform it into a deterministic optimization problem by employing some powerful mathematical lemmas and then apply the proposed BCD-penaltySCA method to solve the reformulated deterministic problem. The proposed methods are shown by simulations to significantly reduce the computation distortion compared with other benchmarks under considering secure AirComp and reliable multicasting.

Joint Optimization of Secure Over-the-Air Computation and Reliable Multicasting Assisted by a MIMO Untrusted Two-Way Relay Quanzhong Li , Hualiang Luo , and Liang Yang Abstract-In this article, we investigate the joint optimization of secure Over-the-Air Computation (AirComp) and reliable multicasting assisted by a multiple-input-multiple-output-untrusted two-way relay, where artificial noise is employed at the access point (AP) to interfere the relay for ensuring secure AirComp.We aim to minimize the computation distortion at the AP by jointly designing the transmit variables of all the nodes and the aggregation variables at the AP and relay, under the secure AirComp constraint, the reliable multicasting constraint, and the transmit power constraints of all the nodes.We consider two scenarios that perfect and imperfect channel state information (CSI) are available.For the former, the formulated optimization problem is highly nonconvex, and we propose an effective block coordinate descent (BCD)-penalty successive convex approximation (penal-tySCA) method to solve the nonconvex problem.For the latter, we model the imperfect CSI by using the worst case criterion and the formulated robust optimization problem is much more challenging than its counterpart with perfect CSI.To solve the robust problem effectively, we first transform it into a deterministic optimization problem by employing some powerful mathematical lemmas and then apply the proposed BCD-penaltySCA method to solve the reformulated deterministic problem.The proposed methods are shown by simulations to significantly reduce the computation distortion compared with other benchmarks under considering secure AirComp and reliable multicasting.Index Terms-Beamforming, channel state information (CSI), multiple-input multiple-output (MIMO), over-the-air computation (AirComp), security optimization, untrusted two-way relay.

I. INTRODUCTION
O VER-THE-AIR computation (AirComp) is a promising technique to aggregate the sensing data from the distributed sensors in the Internet of Things Quanzhong Li is with the School of Computer Science and Engineering, Sun Yat-sen University, Guangzhou 510006, China (e-mail: liquanzh@ mail.sysu.edu.cn).
Liang Yang is with the College of Computer Science and Electronic Engineering, Hunan University, Changsha 410082, China (e-mail: liangy@ hnu.edu.cn).
Motivated by the great success of relaying in wireless communications, relays have been introduced in the evolving AirComp networks [26], [27], [28], [29], [30], [31], [32], [33], [34].However, to the best of our knowledge, an untrusted two-way relay-assisted AirComp is not in the literature and its security issue is worth studying.In this article, we investigate the joint optimization of secure AirComp and reliable multicasting in a MIMO untrusted two-way relay-assisted communication and computation network, where the AP and the relay have multiple antennas and the sensors have a single antenna.The signal transmission between the AP and the sensors is finished in two time slots.In the first time slot, the sensors send their preprocessed sensing data to the relay and meanwhile the AP transmits the multicasting signal plus the artificial noise (AN) with beamformers to the relay, and in the second time slot, the relay broadcasts the received signal multiplied by a beamforming matrix to the AP and sensors.To ensure secure AirComp and reliable multicasting, we jointly optimize the transmit variables of all the nodes and the aggregation variables at the AP and relay, with the objective to minimize the computation distortion at the AP under the constraints that the minimum computation distortion at the relay is larger than some given threshold, the signal-to-interferenceplus-noise ratio (SINR) at each sensor is higher than some predetermined value, and the transmit power at all the nodes cannot exceed their budget.

A. Novelties and Contributions
The novelty of our work comes from the following three main aspects.First, compared with the existing relay-assisted AirComp [26], [27], [28], [29], [30], [31], [32], [33], [34], the proposed untrusted relay-assisted AirComp wants to prevent the relay from wiretapping the aggregation of the sensors' data and thus has an additional security constraint based on the computation distortion, which is highly nonconvex and make the optimization problem much more challenging to deal with.Second, conventional untrusted relay-assisted communication needs to decode each data of the sources (sensors), aiming to maximize the secrecy rate or improve the communication reliability [36], [37], [38], [39], [40], [41], while our proposed untrusted relay-assisted AirComp focuses on the secure aggregation of the sensors' data instead of decoding individual data and its performance metric is the computation distortion in the received value of the desired function, which is a totally different objective from convectional untrusted relay-assisted communication.Third, the existing relay-assisted AirComp [26], [27], [28], [29], [30], [31], [32], [33], [34] only consider the case of perfect channel state information (CSI) which is only applicable to a specific scenario, while our work considers both cases of perfect and imperfect CSI and has a greater flexibility to adapt to different application scenarios.
The main contributions of our work are summarized as follows.
1) We propose a MIMO untrusted two-way relay-assisted secure AirComp and reliable multicasting network, where we employ the multicasting signal and the AN to interfere the relay to ensure secure AirComp.We formulate the MSE minimization problem with the nonconvex secure AirComp constraint and reliable multicasting constraints, considering both perfect and imperfect CSI. 2) With perfect CSI available, we propose an effective block coordinate descent (BCD)-penalty successive convex approximation (penaltySCA) method to solve the nonconvex MSE minimization problem, where the aggregation variables are found in closed form and the transmit variables are optimized by the SCA or penal-tySCA method in which a second-order cone programming (SOCP) or a semidefinite programming (SDP) is solved.
3) When the CSI is imperfect, we employ the worst case criterion to model the CSI, and the corresponding robust MSE minimization problem is much more challenging than its counterpart with perfect CSI.To make the robust problem tractable, we transform it to a deterministic optimization problem by employing some powerful mathematical lemmas.Since he reformulated deterministic problem is still nonconvex, we apply the proposed BCD-penaltySCA method to deal with it.4) We provide numerical results to show that compared with other benchmarks, the proposed methods have a significant improvement in terms of the computation distortion under the constraints of secure AirComp and reliable multicasting.

B. Organization and Notations
The remainder of this article is organized as follows.In Section III, the secure AirComp and reliable multicasting network model is described.In Section IV, we formulate the optimization problem for the case of perfect CSI and propose an effective BCD-penaltySCA algorithm.In Section V, we formulate the robust optimization problem for the case of imperfect CSI and solve the robust problem by the proposed BCD-penaltySCA algorithm.Simulation results are presented in Section VI, and we conclude this article in Section VII.
Notations: U † , U * , U T , tr(U), and U denote the conjugate transpose, conjugate, transpose, trace, and Frobenius norm of the matrix U, respectively.vec(V) means stacking the columns of the matrix V into a single vector.denotes the Hadamard product.W 0 mean that W is positive semidefinite.{x} denotes the real part of the complex number x.

II. RELATED WORKS
Since the relay can substantially overcome the disadvantage of short coverage of single-hop AirComp networks, relay-assisted (i.e., two-hop) AirComp networks have been studied recently [26], [27], [28], [29], [30], [31], [32], [33].Wang et al. [26] studied a hierarchical Aircomp network assisted by multiple single-antenna relays and proposed a centralized algorithm and a decentralized algorithm with global and partial CSI, respectively.In [27], a part of sensors transmit their sensing signals to the AP with the assistance of a singleantenna cooperative relay, while other sensors directly send their sensing signals to the AP, and the impact caused by the relay position has been studied.A cooperative multiantenna amplify-and-forward (AF) relay-aided AirComp network has been considered in [30], while the corresponding case of cooperating with multiple multiantenna AF relays has been investigated in [31].The works [32], [33] mainly focused on the performance analysis of relay-assisted AirComp networks.Jiang et al. [32] analyzed the MSE outage probability and diversity order in the context of AirComp and developed a relay selection scheme to achieve the full diversity order, and the work [33] studied wireless edge federated learning system based on relay-assisted AirComp and characterized the decoding failure probability at the AP.Different from the above two-hop AirComp networks, the work [34] proposed a multihop digital AirComp network where the authors derived the computation rate and proposed a time allocation and power control scheme to improve the network performance.
In the aforementioned scenarios, the relay is considered to be trusted, in the sense that, the relay will not infer the aggregation result of the sensing data from the sensors.However, when the relay is interested in the aggregation result and tends to decipher the result for the unallowed operations, despite operating with the desired relaying protocols, the relay becomes untrusted [35].To maximize the secrecy rate, the security issue of an untrusted relay-assisted communication has been widely investigated in [36], [37], [38], [39], [40], and [41].In [36], the secrecy outage probability and ergodic secrecy rate in a cooperative network in presence of multiple untrusted single-antenna relays and multiple passive eavesdroppers has been investigated.A nonorthogonal multiple access scheme has been considered in [37], where the single-antenna relay that assists the communications between the source and the far user is assumed to be untrusted.In [38], both scenarios of single-antenna and multiple-antenna relay are investigated, and the analytical expression for a lower bound on the ergodic secrecy sum rate has been derived.In [39], both cooperative and noncooperative secure beamforming schemes with an untrusted MIMO AF relay have been investigated, and the result in [39] has been extended to MIMO two-way untrusted relay systems [40].Later, the MIMO twoway untrusted relay in [40] has been considered to perform full-duplex operation in [41].
With regard to the security design for AirComp, in [42], a friendly jammer is used to against passive eavesdropping, where the jammer's signal is reconstructed and fully canceled by the legitimate receiver but deteriorates the eavesdropper's signal-noise ratio (SNR), and thus inhibits the illegitimate receiver's ability to estimate the value of the objective function.Different from the secrecy rate of the physical layer security defined in traditional wireless communications, a δ-semantically secure based on the total variation norm on signed measures and a V-MSE-secure are defined in [42], which means that the estimation MSE at the eavesdropper is at least V under a uniformly distributed objective, regardless of which estimator the eavesdropper uses.In [43], a multiantenna full-duplex AP utilizes the AN to degrade the eavesdropper's links while receiving sensors' preprocessed sensing signals.

III. SYSTEM MODEL
Consider a MIMO untrusted two-way relay-assisted AirComp and multicasting networks as shown in Fig. 1, where K sensors want to send their sensing data (e.g., temperature and humidity) to the AP for aggregation and meanwhile the AP transmits a common signal to all the sensors for some multicasting applications (e.g., updating the status of sensors), both with the help of an untrusted two-way relay.We assume that the sensors are equipped with a single antenna, and the relay and AP has N r and N p antennas, respectively [30].
The transmission of data is finished in two consecutive time slots.In the first time slot, all the sensors send their data to the relay simultaneously, in which the transmitted data by the kth sensor is denoted as where s k is the sensing data with zero mean and unit variable, and b k is the corresponding transmit coefficient.Meanwhile, the AP transmits a common signal s c to the relay, and in order to protect the aggregation of sensors' data from being wiretapped by the untrusted relay, the AP also transmits an AN signal z to interfere the relay.Thus, the transmitted signal by the AP is expressed as where w and V is the beamforming vector/matrix for s c and z, respectively, and s c ∈ CN (0, 1), z ∈ CN (0, I).
From ( 1) and ( 2), the received signal at the relay in the first time slot is given by where h k and H p denote the channel from the kth sensor and AP to the relay, respectively, and n r ∈ CN (0, σ 2 r I) is the additive noise at the relay.
In the second time slot, the signal received at the relay is multiplied by a beamforming matrix F and then forwarded to the sensors and AP.From (3), the signal received at the AP is given by where g k and G p denotes the channel from the relay to the kth sensor and AP, respectively, and n k ∈ CN (0, σ 2 k ) and n p ∈ CN (0, σ 2 p I) are the additive noise at the kth sensor and AP, respectively.
To achieve secure AirComp, the minimum computation distortion measured by MSE at the untrusted relay should be larger than some given threshold, such that the relay cannot compute the aggregation of the sensors' data accurately enough.Considering the relay is interested to compute the arithmetic sum of the sensors' data, i.e., s = K k=1 s k , and applies the aggregation beamforming vector a to estimate the value of s as ŝ = a † y a . ( Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. Based on (3) and ( 5), the computation distortion at the relay is given by [26], [27], [28] MSE from which we can see that the interference from the AP can effectively increase the computation distortion at the relay and thus provide computation security against the untrusted relay.At the AP, the sensors' data are aggregated by a received beamforming vector d.We assume that the CSI about H p and G p are accurate, which is reasonable when the AP and relay are fixed and the channels are estimated through pilot signals with high transmit power [26], [27], [28].Then, the AP can completely remove the self-interference signal ws c +Vz in (4), and the remaining signal at the AP is given by From ( 7), the MSE for the AP to estimate the arithmetic sum of the sensors' data is expressed as From ( 6) and ( 8), we can see that the signal transmitted by the AP will only affect the computation distortion at the relay not that at the AP, which possibly provides an effective method for computation security against the relay while not interfering the data aggregation at the AP.
At the kth sensor, the received signal given by Note that y k contains the signal s k sent by the kth sensor in the first time slot, and s k is known at the kth sensor.Thus, the kth sensor can eliminate the self-interference caused by s k before decoding the common signal s c .Let r k denote the achievable information rate of the kth sensor, and we use the achievable information rate of each sensor to measure the quality of the multicasting.Obviously, the effect of self-interference elimination depends on the corresponding CSI, i.e., h k and g k , and the cases of perfect and imperfect CSI will be discussed in the rest of this article.Based on the above, we formulate the optimization problem as min In problem (10), we aim to minimize the computation distortion at the AP [i.e., (10a)], based on the conditions that the minimum computation distortion at the untrusted relay should be larger than a given threshold ξ r [i.e., (10b)] to protect the computation security, the achievable information rate of kth sensor is greater than a target transmission rate R k [i.e., (10c)] to ensure the reliable multicasting, and the transmit power constraints at the sensors, AP and relay [i.e., (10d)-(10f)] with the corresponding power budget {P k }, P ap and P r .
Remark 1 (Encryption Technique): Besides the PLS, another security technology, called the encryption technique, is also widely used [44], [45], [46], [47].The homomorphic encryption, which has the property that the result of decryption after adding ciphertext is equivalent to the addition of plaintext, could be applied in this AirComp network.Combining with the PLS technology, the ciphertext received by the relay may be scrambled.Due to the avalanche effect [44], it is more difficult for the untrusted relay to recover the plaintext from the received ciphertext, and the protection of the aggregation result has been strengthened.As for some classic homomorphic encryption methods, such as the elliptic curve cryptography [45] and Paillier [46], we need to change the computing function of AirComp into a polynomial to adapt to their operation on the ciphertext.In addition, since the encryption methods, especially some nonlightweight methods, may bring a large computational burden to the sensors [47], and whether the encryption method can be applied to the network depends on the sensors' computing resources and power budget.

IV. JOINT OPTIMIZATION WITH PERFECT CSI
In this section, we assume that the CSI about {h k } and {g k } is perfect, which is reasonable when the sensors are static and the channels vary slow [29], [30], [31].With perfect CSI, the kth sensor can eliminate its own data s k in (9), and the remaining signal at the kth sensor is given by From (11), the SINR to decode the multicasting signal s c at the kth sensor can be expressed as which should be designed to ensure the reliable multicasting.Combining (12) and let ρ k = 2 R k −1 denote the target SINR corresponding to the given target transmission rate of the kth sensor, the optimization problem (10) The challenges of solving problem (13) lie in the following two main aspects.First, the optimization variables in problem (13) are highly coupled in the objective and constraints.Second, the computation security constraint (13b) and the reliable multicasting constraint (13c) are highly nonconvex even with decoupled optimization variables.In general, there is no standard method for solving such nonconvex optimization problems optimally.Thus, in the following, we propose an effective BCD-penaltySCA algorithm to solve problem (13), which decouples the coupling optimization variables and makes problem (13) more traceable by splitting the problem into several subproblems, and then solves them alternatively.

A. Optimizing the Aggregation Vectors {d, a}
With given the variables {F, V, w, b k }, the optimization problem with respective to the aggregation vectors {d, a} are two unconstrained convex quadratic problems whose optimal solution can be obtained just by setting the first-order derivations of the objectives to being zero as where

B. Optimizing the Relay Beamforming Matrix F
With given the variables {d, a, V, w, b k }, the optimization problem with respective to the relay beamforming matrix F can be expressed as min where Problem ( 20) is a nonconvex optimization problem since the left hand side of each constraint in (20b) is the difference of two convex quadratic functions.To solve problem (20) effectively, we employ the SCA method to (20).
According to the property of the convex function, we have the following inequality: where F (n) is the optimal solution at the nth iteration.By replacing g † k FR (1/2) 2 in (20b) by the lower bound [i.e., the right hand side of ( 24)], we solve the following convex problem at the (n + 1)th iteration: where k is a constant.Note that constraint (25b) holds is a sufficient condition for constraint (20b) to hold, and the two constraints are equivalent when the optimization variable F converges.
The optimization problem ( 25) is a convex quadratically constrained quadratic problem (QCQP), and we can recast it as an SOCP as min which is solved effectively by using the interior-point algorithm [48] or the CVX software [49].
C. Optimizing the Transmit Variables {V,w, b k } Inserting ( 17) into (15) and given the other variables, the optimization problem with respective to the transmit variables Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
{V, w, b k } can be expressed as min where ξr = K−ξ r and Although the objective (27a) and the transmit power constraints (13d)-(13f) are convex, problem ( 27) is still nonconvex since the secure AirComp constraint (27b) and reliable multicasting constraint (27c) are nonconvex.To solve problem (27) effectively, in this section, we will propose an efficient penal-tySCA method that converts the nonconvex problem into a differential convex (DC) form and then linearizes the nonconvex terms.
Because the nonconvex constraint (27b) contains the inverse operation of the matrix, it is very difficult to solve it directly.In order to facilitate the expression and handling, we introduce some slack variables {X, Y, T, t} and combine ( 18) and ( 19), and then the constraint (27b) can be represented as where Note that the equality constraints (28a)-(28c) are linear and thus are convex, and the left-hand side of the inequality constraint (28e) is a matrix fractional function [48], which means that (28e) is also convex and can be represented as a linear matrix inequality (LMI) However, the equality constraint (28d) is quadratic over T and thus is a nonconvex constraint.To tackle this equality constraint, we need the following result.
Lemma 1: The equality constraint Y = TT † is equivalent to where the first constraint is an LMI and the second is a difference of convex constraint.
Proof: See the Appendix.
Substituting (27b) by (28), the optimization problem ( 27) can be expressed as min (31), (32). (34d) Now, we can employ the penalty method [48] to solve problem (34), where a penalty term is incorporated into the objective function in order to handle the nonconvex constraint min where ς > 0 is a sufficiently large penalty factor.Heretofore, we have transformed the nonconvex constraint (27b) into a series of convex constraints and a DC penalty term.Nevertheless, the optimization problem ( 35) is nonconvex due to the last penalty term in (35a) and the nonconvex constraint (35b) [i.e., (27c)].Since both (35a) and (35b) are in the DC form, we employ the SCA method to solve it.By the convexity, we first have where (V (n) , is the optimal solution at the nth iteration.

D. Proposed Algorithm
In the proposed BCD-penaltySCA algorithm, we solve problem (13) by alternately updating the blocks {d, a}, {F} and {V, w, b k } until convergence.The solution of updating {d, a} is given in closed form.The subproblems of updating {F} and {V, w, b k } are solved by applying the SCA and the proposed penaltySCA, respectively.The details of the proposed BCD-penaltySCA algorithm are summarized in Algorithm 1.

Remark 3 (Complexity of Algorithm 1):
The main computational burden of Algorithm 1 is from solving the SOCP (26) and SDP (40), whose complexity is about O(K 0.5 N 6 m log(1/ε)) and O(K 0.5 N 7 m log(1/ε)) where N m = max(N r , N p ) [50].Thus, the complexity of Algorithm 1 is about O(L 1 K 0.5 N 7 m log(1/ε)) where L 1 is the iterative number for the convergence of Algorithm 1.
Remark 4 (Limitation of the Proposed Algorithm): In general, the proposed BCD-penaltySCA algorithm can only obtain the local optimal solution to problem (13).Although the outer polyblock approximation algorithm in [51] is employed to obtain a global optimal solution for a DC program, it cannot be applied to obtain a globally optimal solution to the DC program (35) due to the nonconvex quadratic and LMI constraints in problem (35).Since problem ( 35) is a subproblem of the original problem (13), how to obtain a global optimal solution to problem ( 13) is still unknown.In addition, the penalty factor ς needs to be selected appropriately according to a specific problem to obtain a faster convergence speed.

V. ROBUST OPTIMIZATION WITH IMPERFECT CSI
In the above section, we consider the available CSI about {h k } and {g k } is perfect.Here, a more practical scenario is considered where the available CSI about {h k } and {g k } is imperfect, which may occur when the sensors are located in the mobile devices.As in [52], [53], and [54], we employ the worst case criterion to model the imperfect CSI, i.e., h k and g k are given by where hk and ḡk is the estimated channel, and h,k and g,k represents the channel error, which is norm bounded with a corresponding given radius h,k and g,k .
In the second time slot, because of the imperfect CSI, the self-interference term at the received signal y k given in (9) cannot be removed completely.Only part of the self-interference term is removed, thus the remaining received signal at the kth sensor is From (43), the SINR to decode the multicasting signal s c at the kth sensor is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where 13), ( 41), (42), and ( 44), the optimization problem (10) can be further expressed as min Problem ( 46) is a robust optimization problem and is more challenging than its counterpart with perfect CSI, i.e., problem (13), mainly due to the following two reasons.One is that the existing channel errors { h,k } and { g,k } make the problem highly nonconvex even when the BCD algorithm is applied.The other is that the constraints become semiinfinite and we cannot deal with them like the deterministic constraints.

A. Deterministic Problem Transformation
To make the robust problem (46) tractable, a key step is to eliminate the channel errors { h,k } and { g,k }.In what follows, we successively handle the terms with channel uncertainty and convert the robust problem into a deterministic problem.
1) Transformation of the Objective Function (46): By using the epigraph reformulation [48], the objective (46a) can be rewritten as min Obviously, the constraint (47b) is still hard to deal with because it is a semi-infinite constraint.In the following, we make our efforts to convert this semi-infinite constraint to the finite one by eliminating the channel errors { h,k }.
For proceeding, we reexpress the MSE p in the constraint (47b) as MSE p = ψ p 2 (48) where and By employing the Schur complement lemma [55], the constraint (47b) is represented as Constraint ( 52) still contains the channel errors.In order to eliminate these errors, we need the following sign-definiteness lemma.
Lemma 2 [54]: Assume C is a Hermitian matrix and given arbitrary matrices {U i , V i } N i=1 , the following semi-infinite LMI: holds if and only if there are Choosing the parameters appropriately in Lemma 2 as we can rewrite the constraint (52), i.e., (47b), as an LMI where

2) Transformation of the Constraint (46b):
To make the constraint (46b) tractable, we introduce slack variables and transform the constraints involving channel uncertainty into a more straightforward form min a max Due to the infinite number of constraints parameterized by h k , we need to borrow the following S-Lemma.
Lemma 3 [56]: Define the functions Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. where provided that there exists a point y 0 such that h 1 (y 0 ) < 0. Applying Lemma 3 again to the constraint (56), we have where τk = τ k + 2 {b k a † hk } − 1, and 3) Transformation of (46c): Since there are fractions in constraint (46c) and both the numerator and denominator contain channel uncertainty, we cannot rewrite the constraints as the form that the norm of a vector is less than or equal to some value and apply the Lemma 2 to convert them into LMIs.By introducing the slack variables {ω k }, we rewrite (46c) as For (59a), neglecting higher order error terms, we can rewrite it as where and Then, (60), i.e., (59a), is rewritten as an LMI ⎡ where T .(66) For (59b), its form is similar to the constraint ( 56), and we apply the Lemma 3 again.For ease of exposition, we insert g k = ḡk + g,k into (59b) and recast it as where Applying Lemma 3 to (67), we covert the constraint (59b) to an LMI 4) Transformation of (46f): Similar to (47b), we can use Lemma 2 again and rewrite the robust transmit power constraint (46f) as an LMI where and 5) Deterministic Problem: Until now, we have eliminated all the channel errors, and we can transform the robust problem ( 46) into a deterministic one ), (46e), ( 54), ( 55), ( 57), ( 58), ( 64), ( 68), (69), (70).
Although problem (74) have finite constraints, it is still hard to be solved because the variables are highly coupled.In the following, we propose an effective BCD-penaltySCA method to solve problem (74).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Proposed BCD-penaltySCA Method
Similar to the case of perfect CSI, we decouple the deterministic problem (74) into several subproblems and alternately solve these subproblems by convex technique or penaltySCA method.
1) Optimizing d and {τ k }: With given other variables, the optimization problem with respective to d and {τ k } are two separated subproblems, shown as follows: min d,t,{μ p,k} t s.t.(54) (75) and max The two optimization problems are convex SDPs and can be effectively solved by using the interior-point algorithm [48] or the CVX software [49].
2) Optimizing a and ϒ: With given other variables, the optimization problem with respective to a and ϒ is min Although the objective and ( 57) are convex, problem (77) is still nonconvex due to the equality constraint ϒ = aa † .Using Lemma 1, the equality constraint ϒ = aa † is equivalent to Replacing the nonconvex constraint (77b) with (78) and (79), we can recast problem (77) as Now, we can apply the proposed penaltySCA method to solve (80), and in the (n+1)th iteration, we solve the following convex SDP: where a Optimizing w, V, and {b k }: With given other variables, the optimization problem with respective to w, V and 54), ( 57), ( 64), (69), (70) where ξr = ξ r − σ 2 r a 2 − K k=1 τ k .There are two nonconvex constraints (85c) and (85d) in the optimization problem (85).Following the same technique above, we apply the proposed penaltySCA method to solve problem (85).In the (n + 1)th iteration, the convex SDP that needs to be solved is given by min where a (n) w = FH p w (n) .5) Proposed Algorithm: The proposed BCD-penaltySCA algorithm for solving problem (46) is summarized in Algorithm 2.
Remark 5 (Convergence Analysis of Algorithm 2): Similar to the proposed algorithm (Algorithm 1), updating the variables in steps 3-6 of the proposed algorithm (Algorithm 2), we will obtain a monotonically decreasing sequence of the objective Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Remark 6 (Complexity of Algorithm 2): The computational complexity of Algorithm 2 is mainly from solving the SDPs, whose complexity is about O(K 4 N 7 m log(1/ε)) [50].Thus, the computational complexity of Algorithm 2 is about , where L 2 is the iterative number for the convergence of Algorithm 2.
Remark 7 (The Choice Between Robust and Approximate CSI): Another way to handle the channel uncertainty is to use the approximated CSI.As in [57], [58], and [59], the AP can act as a central processor to perform the optimization process.Thus, the AP has the ability to estimate the error caused by using the approximate CSI and decides whether to use Algorithm 1 to obtain a approximate solution with a lower computational complexity or Algorithm 2 to obtain a robust solution.

VI. SIMULATION RESULTS
In this section, we provide the numerical results to verify the effectiveness of the proposed schemes for an untrusted relayassisted AirComp network.As in [26], we assume that all the entries of the channel matrices/vectors (H p , G p ) and (h k , g k ) are independent and identically distributed complex Gaussian random variables with zero mean and variances d −n 0 and d −n k , respectively, where d 0 and d k denote the distance from AP to relay and from relay to the kth sensor, respectively, where n = 3.5 is the path-loss exponent.For simplicity, we also assume that d 0 = 1, d k is randomly generated from [0.5, 1.5], and the powers of sensors and noises, the SINR thresholds, and the channel error radius are the same, i.e., P k = P s , σ 2 r = σ 2 p = σ 2 k = σ 2 , ρ k = ρ, and h,k = g,k = .We set the antenna number of the AP and relay is N p = N r = 20, and the number of sensors is K = 10 [30], [31].

A. Performances of Proposed Scheme for Perfect CSI
With perfect CSI available, we compare our proposed scheme with three benchmarks, including the zero-forcing reception and zero-forcing transmission scheme proposed in [60] (denoted as "ZFR-ZFT") and the maximal-ratio reception and maximal-ratio transmission scheme proposed in [61] (denoted as "MRR-MRT") and our proposed scheme with randomly generated AN (denoted as "RandAN").If not specified, we set the transmit power budget of the sensors, relay and AP is P s /σ 2 = 10 dB, P r /σ 2 = 30 dB, and P ap /σ 2 = 40 dB,  the minimum MSE constraint at the relay is ξ r = 3, and the SINR threshold at the sensors is ρ = 2 dB.
In Fig. 2 , we present the average computation MSE at the AP under different minimum MSE constraint at the relay.From Fig. 2 , we see that the proposed scheme has a significant reduction of the computation MSE at the AP over the "RandAN" scheme for different minimum MSE constraint at the relay, which indicates that the AN plays an important role to satisfy the secure AirComp constraint at the relay.Furthermore, the RandAN scheme outperforms the "ZFR-ZFT" and "MRR-MRT" schemes, which means that using specified relay beamforming schemes cannot reduce the computation MSE at the AP effectively.In Fig. 2 , we also see that when the minimum MSE constraint at the relay is not too large, the average computation MSE at the AP by all the schemes grows not too fast, however, when the minimum MSE constraint at the relay is larger, the average computation MSE at the AP has a relatively rapid increase, which is maybe because that when the minimum MSE constraint at the relay is too large, most of the power is allocated to satisfy the secure AirComp constraint and thus significantly increases the average computation MSE at the AP.In Fig. 3, we present the average computation MSE at the AP under different SINR threshold at the sensors.From Fig. 3, we can see that the proposed scheme is better than all the other benchmarks for different reliable multicasting requirements.We can also see from Fig. 3 that as the SINR threshold at the sensors increases, the average computation MSE at the AP by all the schemes grows, which is because that as the SINR threshold at the sensors increases, more power has to be allocated to satisfy the reliable multicasting requirement and then less power is allocated to reduce the computation MSE at the AP.
In Fig. 4, we present the average computation MSE at the AP under different transmit power budget of the sensors.From Fig. 4, it is seen that for different transmit power budget of the sensors, the proposed scheme outperforms all the other benchmarks in terms of the average computation MSE at the AP.We can also see from Fig. 4 that as the transmit power budget of the sensors increases, the average computation MSE at the AP by all the schemes decreases, which is because as the transmit power of the sensors increases, the AP is easier to align the signals of all the sensors to compress the calculation error.
In Fig. 5, we investigate the average computation MSE at the AP under different network configuration.If we only increase the antenna number of AP and relay, N p and N r , the network obtains more antenna resources and has greater space freedom to optimize the beamforming design at the AP and relay, thus achieving better network performance [62].But, if we only increase the number of sensors, K, the network needs to support more sensors to achieve reliable multicast, and at this time, fewer transmit antennas at the AP and relay are difficult to support more sensors [63].Therefore, as the number of sensors increases, the antenna number of AP and relay should also increase accordingly, and in Fig. 5, we set N p = N r = 1.1K.From Fig. 5, we can see that as the network size increases, the average computation MSE at the AP by our proposed algorithm slowly increases.Thus, in order to meet the accuracy requirements of AirComp in practical applications, it is necessary to choose an appropriate network configuration.In Table I, we show the time to compute the solution by our proposed algorithm (Algorithm 1) under different network configuration.The simulations were performed in MATLAB on a Windows desktop with four Intel i3 cores and 16 GB of RAM.From Table I, we can see that as the network size increases, the computing time of our proposed algorithm (Algorithm 1) increases, which is because Algorithm 1 needs to solve a larger scale SOCP or SDP.In order to accelerate the computing of the optimization solutions, a feasible approach is to design a parallel algorithm for Algorithm 1, such as using the consensus alternating direction method of multipliers (ADMMs) [64] to solve the SOCP and using parallel primal-dual interior-point method [65] to solve the SDP.This is reserved for an interesting future work.

B. Performances of Proposed Scheme for Imperfect CSI
With imperfect CSI available, we compare our proposed robust scheme (dented as "Proposed Robust") with the proposed robust scheme with randomly generated AN (denoted as "Robust RandAN"), and the proposed scheme with perfect CSI is given which serves as a performance upper (or MSE lower) bound (dented as "Perfect CSI").If not specified, the simulation parameter setting is the same as that in the case of perfect CSI, and a normalized channel error radius is used and set as 0.02.
In Fig. 6, we show the average worst case computation MSE at the AP under different minimum MSE constraint at the relay.From Fig. 6, we can see that the proposed robust scheme is close to the scheme with perfect CSI, which indicates that the proposed robust scheme can effectively reduce  the performance loss raised by the channel errors.We can also see from Fig. 6 that the proposed robust scheme is better than the Robust RandAN scheme, which again shows the importance of the AN to satisfy the secure AirComp constraint even with the channel errors.
In Fig. 7, we show the average worst case computation MSE at the AP under different normalized channel error radius.From Fig. 7, it is seen that the average worst case computation MSE at the AP by the proposed robust scheme and the Robust RandAN scheme increases with the channel error radius increasing.We can also see from Fig. 7 that the performance gap between the proposed robust scheme, the Robust RandAN scheme and the Perfect CSI scheme tends to be greater when the channel error radius goes to be larger.This is because as the channel error radius increases, the robust secure AirComp constraint and reliable multicasting requirement become more strict and thus increase the worst case computation MSE.
In Fig. 8, we investigate the average worst case computation MSE at the AP under different network configuration, where we set N p = N r = 1.1K.From Fig. 8, we can see that as the network size increases, the average worst case computation MSE at the AP by our proposed robust scheme increases slowly.But the average worst case computation MSE by the  Robust RandAN scheme grows more quickly, which is due to the fact that randomly generated AN becomes harder to satisfy the secure AirComp constraint at the relay and easier to raise the computation MSE at the AP.From Fig. 8, we can also see that the proposed robust scheme is close to the scheme with perfect CSI, which indicates again that the proposed robust scheme can effectively reduce the performance loss due to the channel errors.
In Table II, we show the time to compute the solution by our proposed robust algorithm (Algorithm 2) under different network configuration.The simulations were performed on the same desktop as Algorithm 1. From Table II, we can see that as the network size increases, the proposed robust algorithm (Algorithm 2) needs to solve larger scale SDPs, and thus the computing time increases.Similar to Algorithm 1, in order to accelerate the computing of the solutions by robust algorithm (Algorithm 2), a feasible approach is to design a parallel algorithm for Algorithm 2, e.g., using parallel primal-dual interior-point methods [65] to solve the SDPs.

VII. CONCLUSION
In this article, we have studied the joint optimization of secure AirComp and reliable multicasting in a MIMO untrusted two-way relay-assisted computation and communication networks, where AN is employed at the AP to interfere the relay for ensuring secure AirComp.We have formulated the optimization problems for two cases of perfect and imperfect CSI.We propose an effective BCD-penaltySCA algorithm Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
for the case of perfect CSI, while for the case of imperfect CSI, we first transform the robust optimization problem to a deterministic problem and then employ the proposed BCD-penaltySCA algorithm to solve the reformulated deterministic problem.Numerical results have demonstrated that the proposed schemes outperform other benchmarks in terms of the computation distortion under the constraints of secure AirComp and reliable multicasting.
(IoT) network by exploiting the supposition of wireless signals over the multiple Manuscript received 8 March 2023; revised 6 May 2023; accepted 9 May 2023.Date of publication 11 May 2023; date of current version 25 September 2023.This work was supported in part by the National Natural Science Foundation of China under Grant 62272493 and Grant 61802447; and in part by the Guangdong Basic and Applied Basic Research Foundation under Grant 2023A1515011201.(Quanzhong Li and Hualiang Luo contributed equally to this work.)(Corresponding author: Quanzhong Li.)

Fig. 2 .
Fig. 2. Average computation MSE at the AP under different minimum MSE constraints at the relay.

Fig. 3 .
Fig. 3. Average computation MSE at the AP under different SINR thresholds at the sensors.

Fig. 4 .
Fig. 4. Average computation MSE at the AP under different transmit power budgets of the sensors.

Fig. 5 .
Fig. 5.Average computation MSE at the AP under different network configurations where N p = N r = 1.1K.

Fig. 6 .
Fig. 6.Average worst case computation MSE at the AP under different minimum MSE constraints at the relay.

Fig. 7 .
Fig. 7. Average worst case computation MSE at the AP under different normalized channel error radii.

Fig. 8 .
Fig. 8. Average worst case computation MSE at the AP under different network configurations where N p = N r = 1.1K.

TABLE II COMPUTING
TIME OF ALGORITHM 2 UNDER DIFFERENT NETWORK CONFIGURATIONS