STAR-RIS-Aided Mobile Edge Computing: Computation Rate Maximization with Binary Amplitude Coefficients

We consider a STAR-RIS-aided MEC system as shown in Fig. 1, in which there are an AP with $N$ antennas, $K$ single-antenna users, and an $M$-element STAR-RIS. The AP is attached to a MEC edge server. Since the STAR-RIS can provide full-space coverage by allowing simultaneous transmission and reflection of the incident signal, it serves both $T$ users in the transmission space and $R$ users in the reflection space.

Each user has a limited energy budget but has intensive computation tasks to deal with, and the available energy for computing the task in user $k$ is denoted as $E_{k}$ in Joules (J). In practice, if the computation task is too complicated to be completed by a user (possibly due to excessive energy or time involved), part of the task should be delegated to the edge server. Hence, a partial offloading mode is adopted [3]. Partial offloading is designed to cope with computation tasks that can be arbitrarily divided to facilitate parallel operations at users for local computing and the AP for edge computing. Being grid-powered, the MEC server has strong computation and storage capabilities for helping users to compute their offloaded tasks and report the computation results in a negligible time.²²2In practice, as long as the overall computation workload of users doesn’t exceed the MEC server’s computing capacity, the computation time can be disregarded as the computation at MEC servers is conducted simultaneously with the uplink data transmission[5, 19, 20].

For uplink transmission, let $\mathcal{T}=\{1,\ldots,T\}$, $\mathcal{R}=$ $\{T+1,\ldots,K\}$ and $\mathcal{K}=\{\mathcal{T}\cup\mathcal{R}\}$ denote the index sets of the $T$ users in transmission space, the $R=K-T$ users in reflection space, and all the users, respectively. Let $s_{t}$ and $s_{r}\in\mathbb{C}$ with zero mean and unit variance denote the information symbols of user $t\in\mathcal{T}$ and user $r\in\mathcal{R}$ for the offloading task, and $p_{t}$ and $p_{r}$ denote the transmit power of user $t$ and user $r$, respectively. Note that all the users with offloading requirements transmit their signals simultaneously, and thus we can express the corresponding received signal $\boldsymbol{y}\in\mathbb{C}^{N\times 1}$ at the AP as

where $\boldsymbol{g}_{t}$, $\boldsymbol{g}_{r}$ $\in\mathbb{C}^{N\times 1}$ are respectively the equivalent baseband channels from the user $t$ and $r$ to the AP, and $\boldsymbol{z}\sim\mathcal{CN}\left(\mathbf{0},\sigma^{2}\boldsymbol{I}_{N}\right)$ is the receiver noise at the AP with $\sigma^{2}$ being the noise power.

With the deployment of a STAR-RIS, the equivalent baseband channels from the user $t$ or $r$ to the AP consists of both the direct link and transmitted or reflected link from the STAR-RIS. Therefore, $\boldsymbol{g}_{t}$ can be modeled as

where $\boldsymbol{h}_{\mathrm{d},t}\in\mathbb{C}^{N\times 1},\boldsymbol{h}_{\mathrm{s},t}\in\mathbb{C}^{M\times 1}$, and $\boldsymbol{G}\in\mathbb{C}^{M\times N}$ are the narrow-band quasi-static fading channels from user $t$ to the AP, from user $t$ to the STAR-RIS, and from the STAR-RIS to the AP, respectively. The STAR-RIS phase shift matrix $\boldsymbol{\Theta}=\operatorname{diag}\left(e^{j\theta_{1}},\ldots,e^{j\theta_{M}}\right)\in\mathbb{C}^{M\times M}$ is a diagonal matrix with $\theta_{m}\in[0,2\pi)$ being the phase shift of the $m^{th}$ element and $m\in\mathcal{M}=\{1,2,\ldots,M\}$. The STAR-RIS transmission amplitude coefficient matrix $\boldsymbol{\Omega}_{t}=\operatorname{diag}\left(\rho^{\mathrm{t}}_{1},\ldots,\rho^{\mathrm{t}}_{M}\right)\in\mathbb{R}^{M\times M}$ is a diagonal matrix with $\rho^{\mathrm{t}}_{m}$ being the transmission amplitude coefficient of the $m^{th}$ element.

Similar to (2), $\boldsymbol{g}_{r}$ is expressed as

where $\boldsymbol{h}_{\mathrm{d},r}\in\mathbb{C}^{N\times 1}$, $\boldsymbol{h}_{\mathrm{s},r}\in\mathbb{C}^{M\times 1}$ denote the narrow-band quasi-static fading channels from the user $r$ to the AP, and from the user $r$ to the STAR-RIS, respectively. The STAR-RIS reflection amplitude coefficient matrix $\boldsymbol{\Omega}_{r}=\operatorname{diag}\left(\rho^{\mathrm{r}}_{1},\ldots,\rho^{\mathrm{r}}_{M}\right)\in\mathbb{R}^{M\times M}$ is a diagonal matrix with $\rho^{\mathrm{r}}_{m}$ being the reflection amplitude coefficient of the $m^{th}$ element. In this paper, we consider the mode switching (MS) protocol, where $\rho_{m}^{\mathrm{t}},\rho_{m}^{\mathrm{r}}\in\{0,1\}$ and $\left(\rho_{m}^{\mathrm{t}}\right)^{2}+\left(\rho_{m}^{\mathrm{r}}\right)^{2}=1$. Such an “on-off” type operating protocol is much easier to implement compared to the energy splitting (ES) protocol and therefore more practical in the application scenarios[12]. Moreover, from the users’ perspective, the MS protocol has lower latency than time switching (TS) protocols, which separates two orthogonal time slots for reflection and transmission users and inevitably incurs higher latency for users assigned to the second time slot.

For the offloading task, we introduce an energy partition parameter $a_{k}\in[0,1]$ for user $k\in\mathcal{K}$, and $a_{k}E_{k}$ represents the energy used for computation offloading. Correspondingly, the transmit power of user $k$ for computation offloading is given as

where $L$ is the length of the time slot.

We consider the linear beamforming strategy and denote $\boldsymbol{v}_{k}\in$ $\mathbb{C}^{N\times 1}$ as the receive beamforming vector of the AP for decoding $s_{k}$. Based on (1), the received signal at the AP for user $k$, denoted by ${\hat{s}}_{k}\in\mathbb{C}$, is then given by

The uplink signal-to-interference-plus-noise ratio (SINR) observed at the AP for user $k$ is thus given by

where we denote an energy partition vector $\boldsymbol{a}=\left[a_{1},\ldots,a_{K}\right]^{\mathrm{T}}$ and a reflection and transmission amplitude coefficient vector $\boldsymbol{\rho}=[\rho_{1}^{\mathrm{t}},\cdots,\rho_{M}^{\mathrm{t}},\rho_{1}^{\mathrm{r}},\cdots,\rho_{M}^{\mathrm{r}}]^{\mathrm{T}}$. Then the computation rate of user $k$ due to offloading is

where $B$ is the system bandwidth.

Regarding local computing, the computation rate is $f_{k}/C_{k}$, where $C_{k}$ is the number of required CPU cycles for users to compute 1-bit of input data, and $f_{k}$ is the user’s CPU frequency. We adopt the dynamic voltage and frequency scaling technique for increasing the computation energy efficiency for users through adaptively controlling the CPU frequency for local computing. Specifically, the energy consumption for the user can be calculated as $L\kappa_{k}f_{k}^{2}$ with $\kappa_{k}$ being the effective capacitance coefficient of user $k$. Also, as $(1-a_{k})E_{k}$ represents the energy for local computing, we have $(1-a_{k})E_{k}=L\kappa_{k}f_{k}^{2}$. By expressing $f_{k}$ in terms of other parameters in this equation, we finally obtain the local computation rate of user $k$ as

We aim to maximize the computation rate (both due to offloading and local computation) of all the users in a given time slot of duration $L$ through jointly optimizing the phase shifts in $\boldsymbol{\Theta}$, the reflection and transmission amplitude coefficients in $\boldsymbol{\rho}=[\rho_{1}^{\mathrm{t}},\cdots,\rho_{M}^{\mathrm{t}},\rho_{1}^{\mathrm{r}},\cdots,\rho_{M}^{\mathrm{r}}]^{\mathrm{T}}$, the receive beamforming vectors in $\boldsymbol{V}=\left[\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{K}\right]$, and the energy partition parameters in $\boldsymbol{a}$. The corresponding optimization problem is formally written as

However, problem (9) is a non-convex optimization problem without a closed-form solution as $\boldsymbol{V}$, $\boldsymbol{a}$, $\boldsymbol{\Theta}$, and $\boldsymbol{\rho}$ are highly coupled in (9a). Therefore, we apply the BCD framework to $\mathscr{P}_{0}$ so that each block is handled iteratively. Furthermore, due to the binary nature of reflection and transmission amplitude coefficients $\rho_{m}^{\mathrm{t}}$ and $\rho_{m}^{\mathrm{r}}$, problem (9) is a nonlinear mixed-integer problem, and it usually requires an exponential time complexity to find the optimal solution. To address the above issue, we apply the logarithmic smoothing-based method to solve the binary variables of $\mathscr{P}_{0}$ in the next section.

When other variables are fixed, the subproblem for optimizing $\boldsymbol{\rho}$ is given as

where we replaced (9e) with (10c) since they are equivalent when $\rho_{m}^{\mathrm{t}},\rho_{m}^{\mathrm{r}}\in\{0,1\}$, and $R_{k}(\boldsymbol{\rho})=B\log_{2}\Big{(}1+\gamma_{k}\left(\boldsymbol{\rho}\right)\Big{)},k\in\{t,r\}$ is obtained from (7) by fixing all variables except $\boldsymbol{\rho}$. The binary constraint (10b) can be equivalently transformed into the following constraints[24]

Based on (12), we always have $\rho_{m}^{\mathrm{x}}-\left(\rho_{m}^{\mathrm{x}}\right)^{2}\geq 0$, where the equality holds if and only if $\rho_{m}^{\mathrm{x}}$ is $0$ or $1$, i.e., a binary variable, and the term $\rho_{m}^{\mathrm{x}}-\left(\rho_{m}^{\mathrm{x}}\right)^{2}$ attains its maximum at $\rho_{m}^{\mathrm{x}}=\frac{1}{2}$. Rather than directly handling the nonlinear constraints (11), we add a penalty term $\rho_{m}^{\mathrm{x}}-\left(\rho_{m}^{\mathrm{x}}\right)^{2}$ with a penalty parameter $\gamma>0$. The problem (10) then becomes

where the penalty function introduced this way is “exact” in the sense that the problem (13) and (10) have the same global optimum for a sufficiently large penalty parameter $\gamma$. In terms of the penalty parameter, we employ an increasing penalty strategy to enforce the exact constraints successively.

The above idea of applying a penalty term to enforce binary constraint generally performs well, as verified by the simulations in [25, 26]. However, the penalty term might introduce many undesired stationary points and local optima. To avoid getting stuck in undesired stationary points or local optima during the optimization process, we further consider a smoothing technique.

A popular method to eliminate poor local optima is using smoothing methods to suppress the unsmoothness of the objective function [27]. It aims to transform the original problem with many local optima into one with fewer local optima and thus has a higher chance of obtaining the global optimal solution.

The basic idea of smoothing is to add a strictly convex function (or concave, depending on maximizing or minimizing objective function) to the original objective, i.e.,

where $f(\boldsymbol{x})$ is the original objective function, $\Phi(\boldsymbol{x})$ is a strictly convex function, and $\mu$ is the barrier parameter. If $\Phi(\boldsymbol{x})$ is chosen to have a Hessian that is sufficiently positive definite for all $\boldsymbol{x}$, then $\mathcal{F}(\boldsymbol{x},\mu)$ is strictly convex when $\mu$ is large enough. This is important in smoothing-based methods because the basic idea is to solve the problem (14) iteratively for a decreasing sequence of $\mu$ starting with a large value. With a sufficiently large parameter $\mu$ at the beginning, any local optima of $\mathcal{F}(\boldsymbol{x},\mu)$ is also the unique global optimum and thus is trivial to optimize.

For a binary vector, a strongly convex function $\Phi(\boldsymbol{x})$ is [28]

This function is well-defined when $\boldsymbol{0}<\boldsymbol{x}<\boldsymbol{1}$ and attains its maximum at $x_{c}=\frac{1}{2}$. Based on (14) and (15), we can derive a smoothing-based algorithm for handling (13). Applying (15) to the optimization problem (13), it can be formulated as

Since (16) has a continuously differentiable objective function and a convex feasible set, the projected-gradient (PG) method can be exploited to tackle this problem. It alternatively performs an unconstrained gradient descent step and computes the projection of the update onto the feasible set of the optimization problem. To be specific, the update of the $i^{th}$ iteration is given by

where $\mathcal{Q}\left(\boldsymbol{\rho}\right)$ is the objective function in (16), $\tau(i)$ is the Armijo step size to guarantee convergence, and $\nabla_{\boldsymbol{\rho}}\mathcal{Q}\left(\boldsymbol{\rho}\right)$ is the gradient of $\mathcal{Q}\left(\boldsymbol{\rho}\right)$, with its explicit expression shown in (37) of Appendix B.

On the other hand, to project $\boldsymbol{\rho}\left(i+\frac{1}{2}\right)$ onto the feasible set determined by constraints (16b) and (16c), an update point $\boldsymbol{\rho}(i+1)$ can be obtained by solving the following optimization problem

where $\left[\begin{array}[]{cc}\!\!I_{M}&\!\!\!\!I_{M}\!\!\!\\ \end{array}\right]\in\mathbb{R}^{M\times 2M}$ and we transformed the constraints (16b) into its equivalent compact matrix form (18c). Since the feasible set of $\boldsymbol{\rho}$ is the intersection of a hyperplane and a box, a closed-form solution can be derived based on the Karush-Kuhn-Tucker (KKT) condition and is given by the following Proposition 1, which is proved in Appendix C. Based on (17), (37) and (19), we can iteratively update $\boldsymbol{\rho}$, where the convergent point is guaranteed to be a stationary point of (16).

In order to enforce the smoothness induced by the last term of (16a), the problem (16) is solved for a sequence of decreasing values of $\mu$ since the solution to (16), $\boldsymbol{\rho}^{\star}(\mu)$, is a continuously differentiable function [29]. Initially, it starts with a suitably large $\mu=\mu_{0}$ since it is vital that the iterates move away from an undesired local optimum. After obtaining the solution of (16) with $\mu=\mu_{0}$, it is used as the starting point of solving (16) but with $\mu=\mu_{1}=\eta_{\mu}\mu_{0}$ where $\eta_{\mu}<1$. The iteration goes on until $\mu_{k}$ reaches the tolerance for barrier value $\epsilon_{\mu}$[30]. In terms of $\gamma$, it can be changed synchronously with $\mu$ by a multiplicative factor $\eta_{\gamma}>1$ since both smoothing-based and penalty-based methods use solution $\boldsymbol{\rho}$ from the previous iteration as the starting point of the next iteration. The overall algorithm for solving $\boldsymbol{\rho}$ under the proposed logarithmic smoothing framework is summarized in Algorithm 1. Since the barrier parameter ends with a value that is close to zero[23], for a sufficiently large value of the penalty parameter $\gamma$, Algorithm 1 will obtain a stationary point of the original problem (10)[25]. The performance improvement of Algorithm 1 compared with the conventional penalty-based method will be illustrated later in Section V-A.

When other variables are fixed, the subproblem for updating $\boldsymbol{a}$ is

Conventionally, since $R_{k}\left(\boldsymbol{a}\right)$ in (21a) can be re-expressed as the difference of two concave functions, the DC programming method can be applied to convexify the problem (21) and further solved numerically by the existing convex solvers such as CVX (the second term in (21a) is already concave). However, the above method is a multi-stage iterative optimization algorithm with the outer loop involving DC programming, and the inner loop still requires iterative numerical methods, thus incurring high computational complexity³³3A common practice is to use the interior-point method, and its complexity order is at least $\mathcal{O}(K^{3.5})$.. To overcome this difficulty, we propose a new approach based on a Lagrangian dual reformulation of the energy partition maximization problem and subsequently apply the quadratic transform method. This leads to an algorithm in which each iteration is performed in closed-forms and bisection search rather than being solved with numerical convex solvers. Thus the proposed method is more desirable and computationally efficient.

First, to establish the legitimacy of applying Lagrangian dual reformulation of (21), we provide the following property, which is proved in Appendix D.

Based on Proposition 2, variables $\boldsymbol{a}$ and $\{\eta_{k}\}_{k=1}^{K}$ are updated iteratively. When $\boldsymbol{a}$ is fixed, the function $\mathcal{D}$ is concave with respect to $\eta_{k}$. Therefore, the global optimal $\eta_{k}$ is obtained by setting $\partial\mathcal{D}/\partial\eta_{k}$ to zero, i.e.,

To optimize $\boldsymbol{a}$ when $\{\eta_{k}\}_{k=1}^{K}$ is fixed, we apply the quadratic transform on the fractional term in (23). In particular, when $\eta_{k}$ is held fixed, only the last two terms of $\mathcal{D}$, which are a sum-of-ratio form and a sum of square roots, are related to the optimization of $a_{k}$. In general, if the objective function is a sum of fractional forms with the fraction’s numerator being concave and the denominator being convex (known as the concave-convex form), one can employ the quadratic transform to convert the concave-convex fractional programming into a sequence of convex subproblems[31, 32]. Since the sum of square roots can be considered as a special case of sum-of-ratio form with a denominator equal to 1, by the quadratic transform, we further recast $\mathcal{D}$ to

where $\{y_{k},z_{k}\}_{k=1}^{K}$ are the auxiliary variables. Then according to the iterative update rules of the quadratic transform, with fixed $\boldsymbol{a}$, the optimal $\{y_{k},z_{k}\}_{k=1}^{K}$ can be directly given as

On the other hand, since $\mathcal{D}_{q}$ is concave with respect to $\boldsymbol{a}$ with fixed $\{y_{k},z_{k}\}_{k=1}^{K}$, the iterative update for $a_{k}$ is obtained by setting $\partial\mathcal{D}_{q}/\partial a_{k}$ to zero, i.e.,

where

Note that equation (28) is a strictly decreasing function with respect to $a_{k}\in[0,1]$ since its derivative is derived as $\frac{3}{4}s_{k}(1-a_{k})^{-\frac{7}{4}}-\frac{o_{k}}{2}a^{-\frac{3}{2}}<0$. Furthermore, if $a_{k}$ is 0 or 1, the left-hand side of the equation (28) will go to infinity and minus infinity, respectively. The above characteristics imply that equation (28) must have a unique solution, and a bisection search algorithm can be used to find the optimal $a_{k}^{\star}$. The details are omitted here for brevity.

These updating steps amount to an iterative optimization as stated in Algorithm 2. The objective function (23) is bounded above and monotonically nondecreasing after each iteration. The solution of Algorithm 2 arrives at a stationary point of the reformulated problem (22). Since the problem (22) is equivalent to (21), it is a stationary point of the original problem (21)[32].

When other variables are fixed, the subproblem of $\mathscr{P}_{0}$ for updating $\boldsymbol{V}$ is

Since the uplink SINR of user $k$ only contains its own receive beamforming vector $\boldsymbol{v}_{k}$, we can optimize $\boldsymbol{v}_{k}$ in a parallel manner, and the objective function is simply the SINR in (6). This gives the $k^{th}$ subproblem as $\max_{\boldsymbol{v}_{k}}\gamma_{k}\left(\boldsymbol{v}_{k}\right)=\frac{\boldsymbol{v}_{k}^{\mathrm{H}}\boldsymbol{A}_{k}\boldsymbol{v}_{k}}{\boldsymbol{v}_{k}^{\mathrm{H}}\boldsymbol{B}_{k}\boldsymbol{v}_{k}},$ where $\boldsymbol{A}_{k}=p_{k}\boldsymbol{g}_{k}\left(\boldsymbol{g}_{k}\right)^{\mathrm{H}}$ and $\boldsymbol{B}_{k}=\sum_{i=1,i\neq k}^{K}p_{i}\boldsymbol{g}_{i}\left(\boldsymbol{g}_{i}\right)^{\mathrm{H}}+$ $\sigma^{2}\mathbf{I}_{N}$. It is easy to recognize that this is a generalized eigenvector problem, and its optimal solution $\boldsymbol{v}_{k}^{\star}$ should be the eigenvector corresponding to the largest eigenvalue of the matrix $\left(\boldsymbol{B}_{k}\right)^{-1}\boldsymbol{A}_{k}$.

When other variables are fixed, the subproblem for updating STAR-RIS phase shifts is given by

To handle the modulus constraints in (31b), the semi-definite relaxation (SDR) method and manifold method are widely adopted. However, SDR-based methods only yield an approximate solution without an optimality guarantee and incur a heavy computational burden as the variable’s dimension increases. Furthermore, manifold optimization methods slow down the convergence due to their nested loop architecture[33]. Therefore, we propose a gradient descent (GD) method algorithm as follows. In contrast to the SDR and manifold optimization-based methods, this method avoids convex relaxation, does not need to solve high-dimensional SDP problems, and gets rid of the nested loops. Thus, it significantly improves computational efficiency.

Observe that the unknown variable $\boldsymbol{\Theta}$ in the feasible set (31b) is in fact $\left\{\theta_{m}\right\}_{m=1}^{M}\in[0,2\pi)$. Therefore, problem (31) can be recast into an unconstrained optimization problem as

where we drop the constraints $\theta_{m}\in[0,2\pi),\forall m\in\mathcal{M}$ since the objective function is periodic and taking modulus of $2\pi$ can recover a $\theta_{m}\in[0,2\pi)$. Since (32) is an unconstrained optimization problem with a differentiable objective function, it can be readily solved by the GD method to obtain a stationary point [34]. To be specific, the update of $\theta_{m}$ at the $i^{th}$ iteration is given by

where $\tau(i)$ is the Armijo step size to guarantee convergence and $\nabla_{\theta_{m}}\mathcal{H}\left(\{\theta_{m}\}_{m=1}^{M}\right)$ is the gradient of $\mathcal{H}\left(\{\theta_{m}\}_{m=1}^{M}\right)$ given by

The proposed four-block BCD optimization algorithm for solving the computation rate maximization problem $\mathscr{P}_{0}$ in (9) is summarized in Algorithm 3, where the converged point is guaranteed to a stationary point [35].

Notice that the iteration with respect to $\{\theta_{m}\}_{m=1}^{M}$ is based on the GD method, which only involves first-order differentiation. Therefore, it has $\mathcal{O}(M)$ complexity order. On the other hand, the computational complexity of the iteration with respect to $\boldsymbol{\rho}$ is dominated by the gradient step. Hence, the complexity order for updating $\boldsymbol{\rho}$ is also $\mathcal{O}(M)$. Furthermore, the complexity to update $\boldsymbol{V}$ and $\boldsymbol{a}$ in Algorithm 2 is $\mathcal{O}(KN^{3})$ and $\mathcal{O}(K)$, respectively [36]. Based on the above discussions, the complexity order of solving (9) is linear in $M$ and $K$. This makes the proposed algorithm suitable for massive elements and user networks. The computation complexity of the proposed algorithm is summarised in Table I.

First, we demonstrate the convergence behavior of Algorithm 1 compared with the penalty-based method in Section III for solving (10). Note that both the smoothing-based and penalty-based methods are solved by the PG method for a fair comparison.

The superiority of the proposed smoothing-based method in terms of the objective values is shown in Fig. 3(a) under a specific channel realization. It is observed that the proposed algorithm achieves a higher computation rate compared with the penalty-based method. This shows the effectiveness of the proposed smoothing-based method in eliminating the poor local optima and therefore achieving a higher computation rate.

Next, in Fig. 3(b), we compare the convergence behavior of Algorithm 2 and the DC programming with CVX for solving (21). It shows that Algorithm 2 and the DC programming reach almost the same objective value at the end. Although DC programming converges in fewer iterations than Algorithm 2, Algorithm 2 is much more efficient than DC programming on a per-iteration basis since Algorithm 2 updates variables in closed-forms or via bisection search, while DC programming requires solving a convex optimization numerically in each iteration. This is reflected by the computation time of Algorithm 2 to reach a computation rate of $4.61\times 10^{7}$ being only 0.057 s while that of the DC programming is 3.116 s.

In addition, we demonstrate the overall convergence of the proposed algorithm (Algorithm 3) and compare with the following two benchmarks: SDR-DC method for solving ${\{\theta_{m}\}_{m=1}^{M}}$ via SDR and solving $\boldsymbol{a}$ with DC programming, and the ‘Penalty-based Method’ for solving (10) with the penalty-based method as in Section III with exact binary constraints.

The superiority of the proposed algorithm in terms of the computation rate is shown in Fig. 4 under a specific channel realization. It is observed that the proposed algorithm achieves a higher computation rate than the penalty-based method. On the other hand, as the proposed algorithm include an extra smoothing parameters update step, it takes more steps in outer iterations in the BCD algorithm to converge. It is also observed that both the proposed and the penalty-based methods outperform the SDR-DC method in convergence speed and computation rate performance.

In Fig. 5, we show the average performance gap between the proposed and penalty-based methods versus the number of elements in STAR-RIS. We can see that the proposed method outperforms the penalty-based method in terms of the objective values obtained. Moreover, as the number of elements increases, the performance gap widens in general. It illustrates the potential of the proposed smoothing-based method in solving a massive STAR-RIS-aided wireless system.

To further show the low computational complexity of the proposed algorithm, we compare the computation time with SDR-DC and penalty-based methods. As shown in Fig. 6(a) and Fig. 6(b), with the number of elements or users increases, the proposed method and penalty-based method save the computation time to a large extent compared with the SDR-DC method (e.g., more than 10 times differences with 60 elements and 8 times differences with 14 users, in Fig. 6(a) and Fig. 6(b), respectively), and the advantage becomes more prominent as $K$ or $M$ increases. On the other hand, it shows that the proposed algorithm achieves almost the same computation time compared with the penalty-based method for small numbers of $M$ and $K$. As $M$ and $K$ increase, the computation time of the proposed algorithm is only marginally higher than that of the penalty-based method. Both the proposed and penalty-based methods are suitable for large-scale STAR-RIS optimization.

To verify the effectiveness of the STAR-RIS in the proposed MEC system and the performance of the proposed algorithm, we compare the performance with the following schemes.