Achievable Rate Maximization for Underlay Spectrum Sharing MIMO System with Intelligent Reflecting Surface

In this letter, the achievable rate maximization problem is considered for intelligent reflecting surface (IRS) assisted multiple-input multiple-output (MIMO) systems in an underlay spectrum sharing scenario, subject to interference power constraints at the primary users. The formulated non-convex optimization problem is challenging to solve due to its non-convexity as well as coupling design variables in the constraints. Different from existing works that are mostly based on alternating optimization (AO), we propose a penalty dual decomposition based gradient projection (PDDGP) algorithm to solve this problem. We also provide a convergence proof and a complexity analysis for the proposed algorithm. We benchmark the proposed algorithm against two known solutions, namely a minimum mean-square error based AO algorithm and an inner approximation method with block coordinate descent. Specifically, the complexity of the proposed algorithm grows linearly with respect to the number of reflecting elements at the IRS, while that of the two benchmark methods grows with the third power of the number of IRS elements. Moreover, numerical results show that the proposed PDDGP algorithm yields considerably higher achievable rate than the benchmark solutions.

IRS-assisted single-input single-output (SISO) and multipleinput single-output (MISO) spectrum sharing systems were considered in [4]- [8]. However, studies on SU rate maximization for the more challenging IRS-assisted multiple-input multiple-output (MIMO) underlay spectrum sharing systems is still limited. In order to maximize the weighted sum rate of the secondary network in an IRS-assisted MIMO system in an underlay spectrum sharing scenario with multiple secondary receivers (SRs) and a single primary receiver (PR), a block coordinate descent algorithm along with an inner approximation method (BCD-IAM) was proposed in [9]. Similarly, in [10], the authors considered the problem of weighted sum rate maximization of the secondary network for an IRS-assisted MIMO spectrum sharing system consisting of multiple SRs as well as multiple PRs, where they adopted a weighted minimum mean-square error (WMMSE) method and an alternating optimization (AO)-based algorithm to find a suboptimal solution. On the other hand, a suboptimal solution to the problem of achievable secrecy rate maximization for an IRS-assisted MIMO system using AO, barrier method, SDR and minorization-maximization (MM), was presented in [11].
It is noteworthy that the methods proposed to maximize the SU rate for MISO systems in [4]- [8] are not applicable to MIMO systems, unless each receive antenna in the MIMO system is treated as a separate user. However, this reduces the beamforming/multiplexing gains which are major benefits of MIMO systems. The method presented in [9] for MIMO systems is proposed only for the case of a single PR. Similarly, the system considered in [11], when omitting the eavesdropper, reduces to the conventional IRS-assisted MIMO spectrum sharing system. However, the solution proposed in [11] is limited to the scenario where the direct links between the (secondary) transmitter and (primary and secondary) receivers are absent. On the other hand, the AO-based solution proposed in [10] can be inefficient due to the possibly high coupling of design variables, especially when the SU rate performance is dominated by the interference constraints at the PRs. Motivated by the drawbacks of existing methods, in this letter we propose a low-complexity but efficient algorithm to maximize the achievable rate of the SU in an IRS-assisted MIMO underlay spectrum sharing system consisting of multiple PRs, where direct links between the (secondary) transmitter and (primary and secondary) receivers are also present. Our main contributions in this letter are as follows: (i) we propose a penalty dual decomposition based gradient projection (PDDGP) algorithm to maximize the achievable rate; (ii) we provide a convergence proof and complexity analysis of the proposed algorithm; and (iii) we provide extensive numerical results to demonstrate the superiority of the proposed solution over the benchmark Algorithm 1: Gradient projection algorithm to solve (3) for fixed υ and ρ. Input solutions provided in [9] and [10], and to show the effect of different system parameters on the achievable rate.
Notations: Bold lowercase and uppercase letters denote vectors and matrices, respectively. X * , X T , X † , |X|, X and Tr(X) denote the complex conjugate, (ordinary) transpose, Hermitian transpose, determinant, Frobenius norm and trace of a matrix X, respectively. The vector space of complex matrices of size M ×N is denoted by C M×N . diag(x) denotes a square diagonal matrix where the elements of x comprise the main diagonal, vec d (X) denotes the column vector containing the elements of the main diagonal of the matrix X, and ∇ X f (·) denotes the complex-valued gradient of f (·) with respect to (w.r.t.) X * . E{·} denotes the expectation operator. X Y represents that X − Y is positive semidefinite, and Π X (x) denotes the Euclidean projection of the point x onto the set II. SYSTEM MODEL AND PROBLEM FORMULATION Let us first consider an IRS-assisted MIMO system in an underlay spectrum sharing scenario shown in Fig. 1, consisting of one secondary transmitter (ST), one SR, and K PRs. The number of antennas at the ST and SR are denoted by N T and N R , respectively, while the K PRs are each equipped with N P antennas. The communication between the ST and SR is assisted by an IRS consisting of N I low-cost passive reflecting elements. The ST-SR, ST-IRS, and IRS-SR channels are denoted by H TR ∈ C NR×NT , H TI ∈ C NI×NT , and H IR ∈ C NR×NI , respectively. On the other hand, the channel between the ST and the k-th PR (k ∈ K {1, 2, . . . , K}) is denoted by H Tk ∈ C NP×NT , while the channel between the IRS and the k-th PR is denoted by H Ik ∈ C NP×NI . Let θ [θ 1 θ 2 . . . θ NI ] T ∈ C NI×1 denote the phase-shift vector at the IRS, where θ l = e jφ l , l ∈ N I {1, 2, . . . , N I }, and φ l ∈ [0, 2π). Here the ST, SR and IRS constitute the secondary network, and in order to communicate with the SR, the ST uses the spectrum band whose license is owned by a primary network. However, with the help of properly designed beamformers at the ST and IRS, it is ensured that the interference inflicted by the secondary network on the PRs is below a predefined threshold. 1 The signal received at the SR is given by is the transmitted signal vector from the ST (which is intended for the SR), n R ∼ CN (0, σ 2 I) denotes the complex additive white Gaussian noise (AWGN) vector at the SR, and σ 2 is the noise power. Without loss of generality, we assume that the noise at all of the PRs is also distributed as CN (0, σ 2 I). With a slight abuse of notation, in the rest of letter, we make the replacements H TI ← H TI /σ, H TR ← H TR /σ and H Tk ← H Tk /σ. This normalization step will also mitigate potential numerical issues since we avoid dealing with extremely small quantities.
Problem Formulation: Denote the transmit covariance matrix at the ST by X E{xx † }, the maximum transmit power budget at the ST by P max , and the maximum tolerable interference power (normalized to the noise power) at the kth PR by P k . Then, the achievable rate maximization problem for the secondary network can be formulated as follows: subject to Tr(X) ≤ P max , (2b) Note that the constraints in (2b)-(2d) represent the transmit power constraint at the ST, the interference power constraints (IPCs) at the PRs and the unit-modulus constraints (UMCs) at the IRS, respectively. Furthermore, the feasible sets for X and θ are defined as X {X ∈ C NT×NT : Tr(X) ≤ P max ; X 0} and ϑ {θ ∈ C NI×1 : |θ l | = 1, l ∈ N I }, respectively.

A. PDDGP Algorithm
It is easy to see that problem (2) is non-convex due to the non-convexity of the objective as well as constraints (2c) and (2d). If the IPCs were not present, (2) would reduce to the conventional MIMO capacity problem for which the gradient projection algorithm in [13] can be applied to obtain a stationary solution. In fact, it is the coupling of X and θ due to the IPCs in (2c) that makes the problem much more difficult to solve and the gradient projection algorithm proposed in [13] is no longer applicable. The main idea of the existing optimization methods in [9], [10] for solving (2) is to alternately optimize X and θ. However, due to the coupling of X and θ, these methods are not guaranteed to return a stationary solution. Indeed we numerically demonstrate in Section IV that these methods are inferior to the PDDGP algorithm that we propose in the following.
Specifically, in order to cope with the IPCs, we apply the penalty dual decomposition (PDD) method recently presented in [14]. First, let g k (X, θ, P k , s k ) Then it is easy to see that (2c) is equivalent to g k (X, θ, P k , s k ) = 0, s k ≥ 0, ∀k ∈ K. Next, the augmented Lagrangian function can be written as follows: υ k is the Lagrangian multiplier corresponding to the constraint g k (X, θ, P k , s k ) = 0 and ρ is the penalty parameter. Hence, for a given (υ, ρ) an equivalent optimization problem can be formulated as maximize subject to s ≥ 0, (2b), (2d). (3b) Since the IPCs (and thus the coupling of X and θ) is now included in the augmented objective, the optimization variables in (3) are decoupled. As a result, a simple but efficient method can be applied to find a stationary solution to (3) which is based on the alternating projected gradient method (APGM). The APGM is motivated by the fact that the feasible set with respect to individual variables is simple in the sense that the projection onto this set can be expressed in closed form.
First, a projected gradient step for θ is performed while other variables are held fixed. In this regard, the gradient of R υ,ρ (X, θ, s) w.r.t. θ is given by Next using (4) (5) shown on the next page. Suppose that, from the current point θ, after moving along ∇ θRυ,ρ (X, θ, s) with some step size one arrives at a pointθ = [θ 1 ,θ 2 , . . . ,θ NI ] T ∈ C NI×1 . To obtain the next iterate, the projection ofθ onto the set ϑ is needed. It is easy to see that this projection is given by Π ϑ (θ) = [θ 1 ,θ 2 , . . . ,θ NI ] T , where for each l ∈ N Ī Next, a projected gradient step for X is carried out. Following a similar line of argument and using [15, eqns. (6.207), (6.195) and Table 4.3], a closed-form expression for the gradient ofR υ,ρ (X, θ, s) w.r.t. X can be obtained in (7), shown on the next page. Moreover, the projection of a given point Algorithm 2: The PDDGP algorithm. Input: θ 0 , X 0 , s 0 , υ 0 , µ 0 , α 0 , ρ, κ < 1 Output: θ ⋆ , X ⋆ 1 repeat 2 Solve problem (3) using Algorithm 1; ρ ← κρ; 6 until convergence; W ∈ C NT×NT onto the feasible set X can be shown to admit a water-filling solution. Due to space constraints, the details regarding the projection of W onto X are omitted, while the interested reader may refer to [13, Sec. III-C] for details.
We now turn to the optimization over s. It is easy to see that the optimal solution to (3), when θ and X are fixed, is simply given by s n,k = max{0, P k −Tr(Z k X n Z † k )}, ∀k ∈ K and thus a projection step for s is not necessary. The proposed APGM to find a stationary solution to (3) for fixed υ and ρ is summarized in Algorithm 1. Note that in Algorithm 1, µ n and α n denote the step size corresponding to the gradient step for the θ and X update in iteration n, respectively. Appropriate values of µ n and α n can be obtained by a backtracking line search. Specifically, we can set µ n = µ n−1 γ in where i n is the smallest positive integer such that R υ,ρ (X n−1 , θ n , s n−1 ) ≥R υ,ρ (X n−1 , θ n−1 , s n−1 ) where x, y 2ℜ(x † y) and γ < 1. A similar routine can be followed to obtain α n .
After solving (3) for given (υ, ρ), the penalty parameter, ρ, is decreased and the Lagrange multipliers are updated. The overall description of the proposed PDDGP algorithm to find a stationary solution to (3) is outlined in Algorithm 2.

B. Extension to Multiple SRs
We have presented the system model where there is a single SR for the sake of simplicity. However, we remark that the extension to deal with multiple SRs is straightforward. In this case, we can consider the maximization of the weighted sum rate (WSR) of all SRs as in [9], [10] and the transmit power constraint (2b) becomes the sum power constraint. It is easy to check that the projected gradient step for θ requires minimal modifications. Also, the projection onto the set defined by the sum power constraint still admits a water-filling solution. Thus, it is trivial to see that Algorithm 2 can be easily modified to solve the resulting WSR maximization problem. The detailed derivations are skipped here due to space limitation and only the numerical results are presented for the multiple SR case in the next section.

C. Proof of Convergence
The convergence proof of Algorithm 2 follows the same arguments as given in [14], and thus we only summa- rize the main points here. First, for given (υ, ρ), Algorithm 1 generates a strictly increasing objective sequence, R υ,ρ (X, θ, s). Since the feasible set is bounded, the iterates returned by Algorithm 1 converge to a limit point of (3a) to some accuracy ǫ k . Next, it can be shown that the sequence υ k+1 − υ k is bounded. Thus, from the dual update it holds that k∈K g 2 k (X, θ, P k , s k ) = ρ υ k+1 − υ k → 0 as ρ → 0.

D. Complexity Analysis
It is obvious that the complexity of the proposed PDDGP algorithm is dominated by that of Algorithm 1. Thus the complexity of Algorithm 1 is obtained by counting the required number of complex multiplications using the big-O notation. First, the computational complexity associated with ∇ θRυ,ρ (X, θ, s) in (5) is as follows. The complexity for computing the terms Therefore, the total complexity for calculating ∇ θRυ,ρ (X, θ, s) is . We remark that, since the projection onto ϑ has negligible complexity, the additional complexity to find an appropriate value of µ n using the backtracking line search in (8) can be ignored.
Next, the complexity to compute ∇ XRυ,ρ (X, θ, s) is O(N 2 T (N R +KN P )). Additionally, the complexity for projecting X n−1 +α n−1 ∇ XRυ,ρ (X n−1 , θ n , s n−1 ) onto the feasible set X is O(N 3 T ). Therefore, the computational complexity for the projected gradient step for X is O(N R N 2 T + KN P N 2 T + N 3 T ). Moreover, it is straightforward to observe that computing s n requires O(K(N P N I N T + N P N 2 T )) complex multiplications. Therefore, the per-iteration complexity of Algorithm 1 is It is noteworthy that in a practical IRS-assisted MIMO underlay spectrum sharing system, N I ≫ max{N T , N R , N P }, and therefore, the per-iteration complexity of Algorithm 1 can be approximated by O(N I (KN 2 T + 2N R N T + 2KN P )), which is linear w.r.t. N I . Note that the complexity of the algorithms proposed in both [9] and [10] grows as O(N 3 I + KN 2 I max{N P , N T }) which is significantly higher than that of our proposed method. The small-scale fading is assumed to be Rayleigh distributed, and that the path loss between two nodes is modeled as PL = (−30−10ξ log 10 (d/d 0 )) dB, where ξ is the path loss exponent, d is the distance between the nodes, and d 0 (= 1 m) is the reference distance. We assume that the path loss exponent for the ST-SR and ST-PR links is 3.75, whereas that for the ST-IRS, IRS-SR and IRS-PR links is 2.2 [9]. Furthermore, it is assumed that N R = 4, N P = 4, P k = 10 −13 W (∀k ∈ K), the noise power spectral density N 0 = −174 dBm/Hz and the system bandwidth is B = 10 MHz. In Figs. 3-5, the average achievable rate is computed over 100 channel realizations. To run Algorithm 2, we start with ρ = 10 and decrease it as ρ ← κρ where κ = 0.1 when the relative objective progress ofR υ,ρ (X, θ, s) [14,Eqn. (48)] in Algorithm 1 is less than 10 −5 . For Algorithm 1 we set γ = 0.5 in the backtracking line search. Algorithm 2 terminates when the difference between R(X, θ) andR υ,ρ (X, θ, s) is less than 10 −5 .
In Fig. 2, the convergence of Algorithms 1 and 2 are shown for N I = 64 and different values of N T . In the figure, each iteration is comprised of one update of X, θ and s. Each color in Fig. 2 in fact represents the convergence of Algorithm 1 for a fixed ρ whose value is clearly indicated in the figure. We recall that Algorithm 1 aims to maximize the augmented objectiveR υ,ρ (X, θ, s), and thus, for a given ρ, monotonic increase is only guaranteed forR υ,ρ (X, θ, s), but not for the original objective R(X, θ). This point is clearly seen in Fig. 2. We also observe that as the number of iterations increases, the difference between R(X, θ) and the augmented objectiveR υ,ρ (X, θ, s) decreases and approaches zero for sufficiently small ρ, in accordance with the principle of a penalty method. In the figure, we also show the run time for Algorithm 2, where the algorithm converges before 150 ms for N T = 8. Note that we have implemented the algorithm using Python v3.9.7 on a Linux PC with 7.5 GiB memory and Intel Core i5-7200U CPU.
In Fig. 3, we demonstrate the superiority of the proposed PDDGP algorithm over the BCD-IAM and WMMSE approaches proposed in [9] and [10], respectively, for N I = 64. Note that here only the first PR located at (0 m, 0 m) is considered since the algorithm in [9] was proposed for a single PR. It is evident from the figure that all the schemes result in a similar performance for small values of P max ; however, the proposed PDDGP algorithm offers a significantly higher average achievable rate for large values of P max compared to that achieved via the solutions proposed in [9] and [10]. Also, it is expected that the SU rate performance achieved via  the BCD-IAM and WMMSE schemes is the same because both of these solutions are derived from similar principles. The main reason for the inferior performance of the BCD-IAM and WMMSE for large P max is that as P max increases, the SU rate performance is mostly dominated by the IPCs at the PRs. In other words, the IPCs become binding for large P max . Consequently, the coupling between X and θ becomes stronger, and it is known that AO-based optimization cannot guarantee a stationary solution under such a scenario.
Next, in Fig. 4, the average achievable sum rate versus the maximum transmit power budget at the ST (P max ) is shown for N I = 64. In particular, we compare the proposed PDDGP algorithm to the WMMSE-based AO algorithm introduced in [10] for the case of two SRs located at (600 m, 0 m) and (600 m, 5 m). In Fig. 4, it can be clearly seen that the proposed PDDGP always outperforms the WMMSE-based algorithm. Note that the BCD-IAM and the WMMSE-based AO are similar in principle and thus share the same drawbacks. This explains why the WMMSE-based AO in [10] is inferior to our proposed algorithm. However, note that different from the case of a single PR in Fig. 3, the PDDGP algorithm outperforms the WMMSE algorithm even in the lower P max regime for the multiple PRs case, and the gains are more significant for larger P max . The benefit of an IRS-assisted MIMO system over that of the system without IRS can also be observed from Fig. 4, as the former results in a considerably higher achievable rate.
To investigate the effect of the number of reflecting elements at the IRS, in Fig. 5 we plot the average achievable rate achieved by the proposed PDDGP algorithm for one SR as a function of N I for N T = 16, N R = 4, N P = 2 and K = 4. Interestingly, for an exponential increase in N I , the achievable rate of the MIMO system also increases near-exponentially. This occurs because for large values of N I , the IRS can create a highly-focused beam to increase the signal-to-noise ratio (SNR) at the SR while satisfying the interference constraints at the PRs, resulting in a higher achievable rate at the SR.

V. CONCLUSION
In this letter, we have presented the optimization problem for achievable rate maximization in an IRS-assisted SU MIMO system in underlay spectrum sharing with multiple PRs, subject to a transmit power constraint at the ST, unit modulus constraints at the IRS and multiple interference power constraints at the PRs. A PDDGP algorithm has been proposed to jointly optimize the transmit beamforming and the IRS phase shifts. The numerical results established the superiority of the proposed algorithm over two existing approaches in terms of the achievable rate. We also showed that the complexity of the proposed PDDGP algorithm is considerably lower than that of the existing approaches.