Robust Transceiver Design for IRS-Assisted Cascaded MIMO Systems

{Robust transceiver design against unresolvable system uncertainties is of crucial importance for reliable communication. We consider a MIMO multi-hop system, where the source, the relay, and the destination are equipped with multiple antennas. Further, an intelligent reconfigurable surface (IRS) is established to cancel the RSI as much as possible. The considered decode-and-forward (DF) hybrid relay can operate in either half-duplex or full-duplex mode, and the mode changes adaptively depending on the RSI strength. We investigate a robust transceiver design problem, which maximizes the throughput rate corresponding to the worst-case RSI under a self-interference channel uncertainty bound constraint. To the best of our knowledge, this is the first work that uses the IRS for RSI cancellation in MIMO full-duplex DF relay systems. The yielded problem turns out to be a non-convex optimization problem, where the non-convex objective is optimized over the cone of semidefinite matrices. We propose a closed-from lower bound for the IRS worst case RSI cancellation. Eventually, we show an important result that, for the worst case scenario, IRS can be helpful only if the number of IRS elements are at least as large as the size of the interference channel. Moreover, a novel method based on majorization theory is proposed to find the best response of the transmitters and relay against worst case RSI. Furthermore, we propose a multi-level water-filling algorithm to obtain a locally optimal solution iteratively. Finally, we obtain insights on the optimal antenna allocation at the relay input-frontend and output-frontend, for relay reception and transmission, respectively.


I. INTRODUCTION
Reliability and throughput are two of the most crucial requirements for the next generation of wireless networks. Optimally relaying the signal from a source to a destination can help enhance reliability and capacity of networks and is currently an active research area [1].
Another emerging candidate for relaying signals are reconfigurable intelligent surfaces (IRSs) [2]. IRS is a device equipped with multiple passive reconfigurable reflectors that can reflect the colliding waves with an adjustable phase. One of the biggest advantages of IRSs is that they work in a real time manner without consuming noticeable amount of power [3]. However, the characteristics of the IRS (e. g. the lack of signal amplification and decode-and-forward processes) can potentially limit its functionality. As a result, in cases where reliability and throughput are of greater importance than the power consumption, conventional relays might still be a better option than IRSs. For instance, authors in [4] show that a simple full-duplex relay can outperform an IRS in terms of throughput under certain conditions. In this paper, we exploit a full-duplex (FD) relay as a primary means of communication, and utilize an IRS to help it through self-interference (SI) cancelation.
IRSs can be utilized in various ways to help the direct links enhance the performance of the system.
In [5], authors have utilized an IRS to maximize the weighted sum rate MISO system. Authors in [6] proposed a method to minimize the power consumption in a MISO system equipped with an IRS and authors in [7] investigate the problem of energy efficiency in an IRS-assisted MISO downlink system. The problem of rate maximization in a MIMO system has been presented in [8], where the authors propose an iterative algorithm to find the best IRS pattern assuming the perfect channel state information (CSI) is given. While the performance of IRS communication systems has been extensively studied, there are not much researches that are considering robust design when the perfect CSI is not available [9].
One of the earliest studies of the robust transmission designs of IRS-assisted systems was done in [10], where a bounded CSI error model is applied to a problem of power miminization in a MISO transmission system. There, by the virtue of semidefinite progtramming (SDP), the authors turn the original problem into a sequence of convex sub-problems. The robust power minimization subject to the outage probability constraints under statistical cascaded channel error model is considered in [11], where the aims is to optimize the system under worst-case rate constraint. Authors in [12] have proposed a robust algorithm for mean squared error (MSE) minimization for a single user MISO system equipped with an IRS. Their method provides a closed form solution for each iteration. However, it can be used only for the case of single user system and cannot be extended to more general cases where there are multiple users. Recently, a rubust algorithm based on penalty dual decomposition (PDD) technique is proposed in [9] for sum-rate maximization where they assumed that the channel estimation error follows a complex normal distribution.
Exploiting a relay to improve communication throughput rate is a classic alternative for IRS in communication systems. However, utilizing a relay in a network raises some important questions to be answered.
For instance, how should the relay process the received signal before dispatching it to the destination? Now, the relay can receive a signal from the source, process it and transmit it towards the destination in a successive manner. This type of relaying technique is known as half-duplex relaying. Alternatively, while receiving a signal at a certain time instant, the relay can simultaneously transmit the previously received signals. This technique is known as full-duplex relaying [13].
As a consequence of transmitting and receiving at a common resource unit, the relay is confronted with SI. Note that, full-duplex relaying potentially increases the total throughput rate of the communication compared to the half-duplex counterpart, only if the SI is handled properly at the relay input. By physically isolating the transmitter and receiver frontends of the relay, a significant portion of SI can be reduced [14].
Moreover, analog and/or digital signal processing at the relay input can be utilized to cancel a portion of SI [15]- [18]. This can be realized if the estimate of the SI can be obtained at the relay. These SI cancellation procedures can effectively mitigate the destructive impact of SI up to a certain level. Hence, the remaining portion, the so-called residual self-interference (RSI), is still present at the relay input. The distribution of the RSI is investigated in [19], [20]. This RSI is mainly due to the channel estimation uncertainties and transmitter noise. Therefore, the quality of channel estimation plays an important role for limiting RSI if the conventional modulation techniques are utilized. Interestingly, the authors in [21] employ a superimposed signaling procedure (asymmetric modulation constellation) in the basic point-topoint FD communication for cancelling the SI and further retrieving the desired information contents without requiring channel estimates. They show that for the same average energy per transmission block, the bit error rate of their proposed method is better than that of conventional ones. The RSI degrades the performance of the communication quality evidently. To this end, the authors in [22] study the degreesof-freedom (DoF), i.e., the slope of the rate curve at asymptotically high SNR, and its relation to the performance of FD cellular network in the presence of RSI. Moreover, the authors in [23] and [24] investigate the joint rate-energy and delivery time optimization of FD communication, respectively, when RSI is still present. Furthermore, the authors in [25] study the sum rate capacity of the FD channel with and without such degradation. In the presence of RSI, the authors in [26] study the capacity of Gaussian two-hop FD relay.
Robust transceiver design against the worst-case RSI channel helps find the threshold for switching between HD and FD operating modes. This set up is commonly known as hybrid relay systems [27]. The authors in [28] investigate a robust design for multi-user full-duplex relaying with multi-antenna DF relay.
In that work, the sources and destinations are equipped with single antennas. Moreover, the authors in [29] investigate a robust transceiver design for FD multi-user MIMO systems for maximizing the weighted sum-rate of the network.
Contribution: Motivated by above, in this work, we consider a DF multi-hop system with multiple antennas at the source, relay and destination along with an IRS to provide additional links. Then we try to maximize the sum rate throughput for the worst case RSI scenario. To the best of our knowledge, this is the first time that the throughput rate maximization against the worst case RSI is evaluated for IRSassisted DF full-duplex relay in MIMO systems. First, we simplify the problem by finding an analytical lower bound for the performance of the IRS. Then, the optimization of maximum achievable rate of the DF full-duplex relaying is cast as a non-convex optimization problem. Thereafter, we propose a low complexity method to find the solution using majorization theory. We propose an efficient algorithm to solve this problem in polynomial time. Finally, the transmit signal covariances at the source and the relay are designed efficiently to improve robustness against worst-case RSI channel in a given uncertainty bound.
Notice that once the covariances are known, one can easily find the precoders using conventional methods such as SVD decomposition etc.

II. SYSTEM MODEL
We consider the communication from a source equipped with N t antennas to a destination with N r antennas. The reliable communication from the transmitter to the destination is assumed to be only feasible by means of a relay with K t transmitter and K r receiver antennas at the output and input frontends, respectively. This means that the direct link from the transmitter to the destination and the link from the transmitter to the IRS and to the destination has a negligible impact on the throughput. This assumption is realistic for the scenarios where the path loss is high due to the high frequency ranges such as mmWave and Terahertz or due to far distances [30]. An IRS consisting of M elements is established to either cancel the RSI or help enhance one of the transmitter-relay/relay-destination links. The system model is given in Fig. 1. As depicted, there are three different scenarios at which the IRS can perform. First, IRS can be used to eliminate the RSI. In this case, the received signals at the relay and destination are given by respectively, where κ ∈ {0, 1}. Notice that, κ = 0 coincides with HD relaying and κ = 1 denotes FD relaying. The transmit signal of the source is denoted by x s ∈ C Nt with the covariance matrix Q s = E[x s x H s ], and the transmit signal of the relay is represented by x r ∈ C Kt , with the covariance matrix The additive noise vectors at the relay and destination are denoted by n r ∈ C Kr and n d ∈ C Nr , respectively, which are assumed to follow zero-mean Gaussian distributions with covariance matrices σ 2 t I and σ 2 d I respectively, where I is identity matrix. The source-relay channel is represented by H 1 ∈ C Kr×Nt and the relay-destination channel is denoted by H 2 ∈ C Nr×Kt . These channels are assumed to be partially known. Similarly, the partially known channels from the relay's transmitter to the IRS and from the IRS to the relay's receiver are denoted by H rI ∈ C M ×Kt and H Ir ∈ C Kr×M respectively, see Fig.   1. Θ ∈ C M ×M is a diagonal matrix representing the IRS phase profile. Throughout this paper, we assume that the IRS is passive and therefore for each element we have |θ m | ≤ 1 [31]. Furthermore, the channels from the source to the IRS and from the IRS to the destination are also assumed to be partially known and denoted as H SI and H ID respectively. Finally, the SI channel at the relay is represented by H r , which is assumed to be known only imperfectly. In what follows, we present the achievable throughput rates for the HD and FD relaying. We start with the HD relay, in which κ = 0. In the second case, IRS can be exploited to enhance the quality of the channel between the source and the relay. In such a case, the received signals at the relay and destination can be expressed as where H SI ∈ C NtM is the channel from the source to the IRS. Finally, IRS can be used to help the channel from the relay to the destination. In this case, the received signals are going to be where H ID ∈ C M ×Nr is the channel from the source to the IRS.
In this paper, we consider the channel uncertainty model used for FD communication in [32]- [34] in which some information about the channel state is provided to the receiver. Therefore, the receiver can use the estimation of the channel, denoted asĤ x where x ∈ {1, 2, SI, IR, RI, ID}. The channel estimation error is then expressed asH x = H x −Ĥ x , whereĤ x andH x are uncorrelated. Further we assume that Fig. 1: System model of a IRS assisted full-duplex relay. In our model both source and destination are equipped with N t and N r antennas respectively. Also relay is equipped with K t transmitting and K t receiving antennas and IRS has M passive elements.
Notice that this assumption is more genereal than zero mean circularly symmetric complex Gaussian entries or non-i.i.d channel estimation errors applied in [35] and [36], [37] respectively.
In what follows, we find the achievable rate for three aforementioned cases and compare them to see under what conditions each of them should be applied.

A. Overview
Suppose that the relay employs a DF strategy. In the full-duplex scenario both source-relay and relaydestination links are active at the same time. As a result, the signals from the relay transmitter interfere with the receiving signal at the relay receiver. We assume that an estimate of the SI channel H r is available at the relay denoted byĤ r . Hence, the RSI is represented byH r is given as

B. Mathematical Preliminaries
Considering a FD DF relay, the following rates are achievable [38], in which depending on how the IRS is applied to the system, the three following sets of rates are possible. First, where H tot 1 = H r + H rI ΘH Ir when the IRS is used to cancel the self interference. Second, whereĤ tot 2 = Ĥ 1 + H rI ΘĤ Ir andH tot 2 = H 1 +H rI ΘH Ir if the IRS is established to help the source-relay channel and finally whereĤ tot 3 = Ĥ 2 +Ĥ rI ΘĤ Ir andH tot 3 = H 2 +H rI ΘH Ir if the IRS is utilized to enhance the rate of relay-destination channel. Notice that, assuming that the RSI remains uncanceled, a robust transceiver against the worst-case RSI channel is required which is formulated as an optimization problem as follows max Qs,Qr,Θ min Hr min R FD sr , R FD rd (15) subject to x ∈ {1, 2, r, RI, IR, ID, SI} (15c) in which the throughput rate with respect to the worst-case RSI channel is maximized. Two constraints P s and P r represent the transmit power budgets at the source and the relay respectively. In constraint (16a), T x represents the RSI or the channel estimation error bound corresponding to H x . Notice that, Tr(H xH H x ) represents the sum of the squared singular values of H x . It should be noted that, using a bounded matrix norm is the most common way for modeling the uncertainty of a matrix [39], [40]. In practice, T x can be found using stochastic methods when the distribution of the channel error is known. Otherwise, one may find it using sample average approximation method. Finally, constraints (15d) are due to the unit modulus limitation of the IRS elements.
The problem (15) is non-convex and hard to solve. As a result, for each of the above-mentioned scenarios, we propose a simplified version of the optimization problem and try to solve it instead. Note that as we are interested in finding the throughput corresponding to the worst case RSI, any simplification in the optimization problem should be in favor of the RSI and interference. In the next, we first analyse the performance of the system when the IRS is helping the relay to cancel the RSI. Consequently, we show that the problem (15) can be simplified to the following optimization problem max Qs,Qr .
and where V ec(·) denotes the vector of all non-zero elements of its input matrix. We can equivalently write T as where * denotes column-wise Khatri-Rao product defined as below and where A i is the i'th column of A and ⊗ denotes the Kronecker product. See Appendix I for proof.
One can show that T ≤ ( √ T r −σ min (H Ir * H T rI )) 2 . As mentioned before, problem (16) is a simplification of the problem (15). This means every achievable rate which is inside the feasible set of (16) is also inside the feasible set of (15) 1 . The intuition behind that is, in problem (15) the minimization over RSI happens only one time and the worst case RSI simultaneously tries to cancel the effect of best configuration of IRS and the best covariance matrices. Whereas in (16), first, the RSI does its worst damage on the performance of the best IRS configuration, and after that, performs another optimization to bring the worst power allocation against the best covariance matrices 2 . In what follows, we provide our proposed ways to deal with optimization problems (17) and (16) respectively.
Proof. We begin the proof with an intuitive example and then extend it to the more general case. Assume that K t = 1, K r = 2 and M = 3. Then we have Also, consider the following optimization problem Here, notice that H Ir diag(H T rI ) is a linear map from a three-dimensional into a two-dimensional space. One simple example of such a mapping can be found in Fig 2. Here, an example of mapping from threedimensional to two-dimensional space is shown. The left shape shows the feasible set for the IRS with three elements in a real valued space. The cube belongs to the case of T , i. e. constraints −1 ≤ θ m ≤ 1, ∀m, while the sphere shows the constraint θ 2 1 + θ 2 2 + θ 2 3 ≤ 1 which belongs to T . In the right, the feasible sets belonging to two aforementioned regions after performing mapping f are presented as an example. It can be seen that the mapping of the first set of constraints (the hexagon) covers the whole area of that of the second one (the ellipse). One important key is, as mapping is a linear function, we have where A and B are two arbitrary sets and f is the mapping.
With the above considerations in mind, now we are ready to deal with the optimization problem. and (13c) imposed by the IRS constraints respectively. Notice that, the latter is always a subset of the former, i.e, S ellipse ⊂ S parallelogram , where S represents the area of the set. Now, considering the max-min problem (20), we are able to have a geometrical interpretation of the solution. Note that (20) indicates the distance betweenh r and −H Ir diag(H T rI )Θ. Whileh r wants to maximize this distance, the job of Θ is to minimize it as much as possible. It can be seen in the figure that the best choice of RSI for h r is the point that is farthest from the parallelogram. It is evident that this point is somewhere on the circumference of the outer circle. On the other hand, the best choice for the Θ is to be chosen in a way to be as close as possible toh r . These points are indicated respectively as G and F on the figure. In general, as the number of IRS elements or the dimensions ofh r increase, the mapping of the hypercube becomes more and more complicated and finding the optimal distance becomes more difficult. However, there is an upper bound for this distance. As shown in the figure, if instead of the cube, we limit the feasible set of IRS elements to the sphere inside the cube, i.e. replacing (12c) with (13c), the solution to the problem becomes GE ≥ GF . It turns out that finding GE is very simple as by the definition we have σ min (H Ir diag(H T rI )) = OE, and also we know that √ T r = GO. Therefore, we can conclude ). Finally, we use one last upper bound to make the original problem even easier to solve. Note that, if instead of the ellipse, we consider the circle inscribed in it, we will have It is worth mentioning that the geometrical representation for the optimization problem (15) is different, because there, considering that the RSI wants to bring the worst representation against the IRS configuration and covariance matrices simultaneously, the RSI cannot freely span the whole circle. This is due to the fact that some regions in the cicle might not be a good choice when it comes to RSI design against covariance matrices. However, if the best representation of RSI against the covariance matrices also provides the best RSI against the IRS configuration, the solution to (16) and (15) will be the same. Eventually, instead of optimizaton problem (15) one can solve optimization probem (16).
The solution to the new problem is guaranteed to be achievable by the original problem as well.
Notice that one can readily extend this interpretation into the complex domain, as the constraint (21c) will still be a subset of constraints (20c). It remains to show one can generalize the geometrical proof for arbitrary large dimensions. This means that the channel dimensions and the number of IRS elements can be any natural numbers. Interestingly, it is enough to show that the geometrical proof based on 2 norms and Euclidean distance exists for higher dimensions. This proof is given in Appendix I where it is shown The next is to solve problem (16). Solving this problem is hard in general as it is non-convex. Hence, we first use the following lemma and theorem to solve it. There, it is shown that for every possible choice of H 1 and H 2 , there exists at least one set of simultaneously diagonalizable matrices H tot , Q s and Q r that are the solutions to the problem (16).

Lemma 1.
For two positive semi-definite and positive definite matrices A and B with eigenvalues Proof. Consider Fiedler's inequality given by [41], Furthermore, given B as a positive definite matrix, the followings hold, Now, dividing the sides of (24) by B , one can readily obtain (23).
Note that, in (23) the inequalities hold with equality if and only if A and B are diagonalizable over a common basis. Using the result of Lemma 1, R FD sr can be lowerbounded as Also it holds that λ i σ 2 Note that, the inequality holds with equality whenever share a common basis. Next, we use the following inequality The above inequality holds true since λ i (H 1 Q sH H 1 ) ≤ T 1 P s . Now, instead of doing the minimization over the left-hand side (LHS) of equation (28), we can first minimize the right-hand side (RHS) of equation (29) to find an achievable rate. Similarly, for R FD rd we have Remark 1. Having the equality C =H r Q rH H r , one can generally conclude the rule of multiplication is determinant i. e. det (C) = det H H rH r det (Q r ). Further, using the properties of determinants we can is a random permutation of i and indicates that there is no need for λ ρ(i) H H rH r to be in decreasing order. However, one cannot generally conclude λ i (C) = λ ρ(i) H H rH r λ i (Q r ) , ∀i, unlessH H rH r and Q r share common basis.
As a result of Remark 1, in general, we cannot rewrite (28) in terms of λ i H H rH r , λ i (Q r ), λ i H H 1 H 1 and λ i (Q s ). However, if we show that for every choice of Q s , there exists a matrix Q s with properties: 1) 3) Tr(Q s ) ≤ Tr(Q s ); then we can use Q s instead and rewrite (28) in terms of λ i H H 1 H 1 and λ i (Q s ) to simplify the problem. The first property implies that both Q s and Q s have the exact same impact on the capacity. Hence, if we find a Q s which is the solution to the problem (16), its corresponding Q s will also be a solution. The second property means, unlike Q s , Q s actually shares the common basis with H H 1 H 1 . The last property implies that Q s is at least as good as Q s in terms of power consumption. Observe that if we show for every feasible Q s there exists at least one such Q s , then we can solve the problem (16) in a much easier way. The reason is, in such a case, instead of searching for optimal Q s over the whole feasible set, we can search for the optimal Q s . Unlike Q s , finding Q s does not need a complete search over the whole feasible set since Q s shares a common basis with H H 1 H 1 . Therefore, we can limit our search only to the portion of the feasible set in which the matrices have eigendirections identical to those of H H 1 H 1 . Similarly, if we show for every choice ofH r , there exist at least oneH r for which we have three conditions where ρ (i) is a random permutation of i and indicates that there is no need for λ ρ(i) (Q s ) to be in decreasing order.
Proof. The proof is given in Appendix II.
For the sake of simplicity, we use following notions for the rest of the paper, Now, using Theorem 2 alongside Lemma 1, we infer that with no loss of generality, instead of optimising over matrices, one can do the optimization over eigenvalues to find the optimal value for RSH of (29).
Then we have subject to γ s 1 ≤ P s , (39a) Note that, the two additional constraints (39d) and (39e) need to be satisfied due to the conditions of Lemma 1 (i.e. eigenvalues have to be in decreasing order). Interestingly, these two additional constraints are affine. The above optimization problem can further be simplified using the following lemma, Lemma 2. The objective function of the optimization problem (39) is optimized when the constraints (39a) and (39c) are satisfied with equality.
Proof. Intuitively, as the objective function is an increasing and decreasing function of each element of γ s and σ 2 r respectively, at convergence, the constraints are met with equality. See Appendix III for the proof.

C. Algorithm Description
Now, we need to solve the optimization problem (39). It can be readily shown that R FD rd is a monotonically increasing function of P r . Furthermore, one can show that R FD sr is an increasing function w.r.t. P s and decreasing function w.r.t. T and P r (See Appendix IV). Consequently, the worst-case RSI chooses a strategy to reduce the spectral efficiency, while the relay and the source cope with such strategy for improving the system robustness. That means, on one hand the RSI hurts the stronger eigendirections of the received signal space more than the weaker ones. However, on the other hand the source tries to cope with this strategy adaptively by smart eigen selection. This process clearly makes the optimization problem complicated at the source-relay hop. Unlike the source-relay hop, the resource allocation problem at the relay-receiver hop is rather easy. Since at the relay-receiver hop there is only one maximization and we can find the sum capacity simply by using the well-known water-filling algorithm.
Observe that although finding each of R FD sr and R FD rd separately is a convex problem, the problem (39) as a whole is not convex. Therefore in this paper, we find the optimal R FD sr by keeping R FD rd fixed. Then we use the resulted R FD sr to find optimal R FD rd and again, using the new resulted R FD rd to find optimal R FD sr . This iterative process repeats until the convergence. Our simulation showed that the algorithm has a very fast convergence and only in rare cases it takes more than 20 iterations for the algorithm to converge. This is mainly due to the fact that inequalities (39d) and (39e) restrict the eigenvalues to vary up to a certain limit, which in turn, makes the whole outputs more stable. As it can be seen, only less than 3% of cases did not converge until 50 iterations.
Notice that the optimum values for the transmission power on relay hop may not sum to P r . The reason is that R FD sr is a monotonically decreasing function of P r and as we are interested in the min(R FD sr , R FD rd ), with R FD sr < R FD rd we will have min(R FD sr , R FD rd ) = R FD sr . Therefore, it is in our interest to keep P r as low as possible to increase R FD sr as much as possible. Analogously, in the case of R FD sr > R FD rd we have min(R FD sr , R FD rd ) = R FD rd which can be increased by increasing the total power usage of relay's transmitter. As a result the well-known bisection method can be used to find the optimal rate where we have R FD sr = R FD rd ,

Algorithm 2
The optimal γ s 1: Find power allocation P 0 using water-filling algorithm 2: while |P (q) − P (q−1) | is large do unless the case R FD sr ≥ R FD rd happens even if the maximum allowed power is used at the relay transmitter. In such a case, the relay-destination link becomes the bottleneck. Now we focus on how to find R FD sr . In order to find the sum rate for the source-relay hop, we assume that we are already given γ r which is the vector of relay input powers that maximizes the sum rate at relay-destination hop. The next step is to do the minimization over σ r and the maximization over γ s . One approach to solve this problem is to solve it iteratively. With this method, first one finds the optimal γ s by solving the maximization part of (39) under the assumption that the optimal σ r is given, and then, having the optimal γ s the minimization part of (39) can be solved efficiently. This process goes on until the convergence of γ s and/or σ r . The maximization part is performed using water-filling method. However, the additional conditions ∀ i ≤ min(M, K r ), σ 2 1 i γ s ρ(i) ≥ σ 2 1 i+1 γ s ρ(i+1) should be taken into account. For instance, if the optimal value for γ s i turns out to be equal to zero, then we should have γ s j = 0 for all j > i irrespective of their SNR. Fig. 5 depicts two different examples of multi-level water-filling algorithm. As it can be seen, first, regular water-filling algorithm is considered where for each subchannel we have as its channel gain. After finding the water-level in this way, we need to impose . These additional restrictions act like caps on top of the water and create multilevel water-filling which can be interpreted as a cave. Fig. 5 (a) shows the case where these caps do not make any subchannel to have zero power. However, Fig. 5 (b) shows the case where subchannel i = 13 has to be zero as a result of the cap imposed by the additional constraints (39d). In this case we have γ s ρ(13) = 0 and as a result min 1≤i ≤13 {σ 2 1 i γ s ρ(i ) } = 0. Thus, this condition forces all other subchannels (i.e. i > 13) to get no power. Algorithm 2 provides the detail of multilevel water-filling. For the minimization part, Lagrangian multiplier is used. We have Calculating ∂L ∂σ 2 r i = 0 we arrive at where λ is the water level. Similarly to the maximization case, there are additional constraints γ r i σ 2 s ρ(i) ≥ γ r i+1 σ 2 r i+1 that must be considered during the minimization process. However, it can be shown that if the constraints γ r i ≥ γ r i+1 and σ 2 1 i γ s ρ(i) ≥ σ 2 1 i+1 γ s i+1 are met, then the constraint γ r i σ 2 s ρ(i) ≥ γ r i+1 σ 2 r i+1 becomes redundant. Please refer to Appendix V for proof.   Next we deal with the optimization for the cases where IRS is utilized to help either the source-relay or relay-destination channels. In such cases, the optimization part over the covariance matrices remains the same as the above-mentioned case. Also, the optimization of the IRS elements can be done using eigenvalue decomposition and the algorithm introduced in [8]. Notice that for the case in which IRS is assisting the source relay link, the term T 1 P s in (29) should be replaced with (T 1 + T SI T IR )P s , and for the case where IRS helps the relay-destination link, the term T 2 P s in (47) should be replaced with (T 2 + T RI T ID )P r .

IV. ACHIEVABLE RATE (HALF-DUPLEX RELAY)
We consider a simple HD relay, where the source and the relay transmit in two subsequent time instances.
Notice that for the case of HD, IRS can be used to assist both the source-relay and relay-destination channels as the singal is being sent over each of these channels in a different time slot. therefore, the received signals at the relay and the destination can respectively be expressd as Consequently, the achievable rates for the transmitter-relay and relay-destination links can be expressed as bellow where H 1 =Ĥ 1 +Ĥ IR ΘĤ SI and H 2 =H 2 +H ID ΘH RI . Also, R HD sr and R HD rd are the achievable rates on the source-relay and relay-destination links, respectively Using time sharing, the achievable rate between the source and destination nodes is given by where α is the time-sharing parameter. Note that, in half-duplex relaying the source and relay transmissions are conducted in separate channel uses. Hence, the transmit covariance matrices Q s ∈ H Nt×Nt and Q r ∈ H Kt×Kt are optimized by maximizing the achievable rate from the source to the destination. Here, the convex cone of Hermitian positive semidefinite matrices of dimensions N t ×N t and K t ×K t are represented by H Nt×Nt and H Kt×Kt , respectively. Importantly, for maximizing this achievable rate, the time-sharing parameter, i.e., α needs to be optimized alongside the system parameters, e.g., power allocation. Readily, optimal α occurs at αR HD sr = (1 − α)R HD rd . Therefore, the achievable rate becomes as follows Notice that as the objective function of the above optimization problem is a monotonically increasing The above optimization problem follows the same approach applied for the optimization of relay-destination link in FD scenario. As a result, the same method could be applied to find it. In other words, the wellknown water-filling algorithm can be used to find the optimal covariance matrices along with the algorithm introduced in [8] to find the best IRS pattern. This process continues iteratively until it finally converges.
The solution to R HD rd is the same as well and the same procedure can be applied to find Q r .

V. NUMERICAL RESULTS
We assume the transmit power budgets at the source and at the relay are P s = 5 and P r = 1 respectively.
Moreover, the AWGN spectral density is assumed to be −175dBm and the bandwidth is BW = 180 MHz.
In this section, we investigate the performance of IRS-assisted full-duplex relaying with RSI channel uncertainty bound T r , i.e., Tr(H rH H r ) ≤ T r . We consider all the channels to follow the Rician distribution with the factor = 0.1 and the specificaiton given in Table I. We also assume T x = 0.001, x ∈ {1, 2, SI, IR, RI, ID}. We perform Monte-Carlo simulations with L = 10 3 realizations from random channels and noise vectors. Hence, the average worst-case throughput rate is defined as the average of worst-case rates for L randomizations, i.e., R av = 1 L L l=1 R l . Notice that, for each set of realizations, we solve the robust transceiver design as is elaborated in Algorithm 1. We run different sets of simulations as described in the following subsections.

A. Antenna Array Increment with no IRS
In this part, first we assume that there is no IRS installed. Then we evaluate the performance of the system usign different strategies. Thereafter, we examine how installing an IRS can help increase the throughput. We consider two cases, where the source, relay and destination are equipped with (a)-small antenna array and (b)-large antenna arrays. In order to see the impact of IRS, we first assume that there is no IRS installed. For these cases, we have (a)-N t = 4, K r + K t = 10, N r = 4, (b)-N t = 10, K r + K t = 24, N r = 10.
These cases are considered to highlight the performance of full-duplex DF relaying as a function of the number of antennas with the worst-case RSI. Interestingly, as the number of antennas at the source, relay and destination increase, full-duplex relaying achieves a higher throughput rate even with strong RSI. This can be seen by comparing rates from Fig. 6(a) to those from Fig. 6(b). Furthermore, notice that the worst-case RSI casts strong interference on the strong streams from the source to the destination. With very low RSI power T r → 0, full-duplex almost doubles the throughput rate compared to the half-duplex counterpart. This can be seen in Fig. 6, where the curves have their intercept point with the vertical axis.
However, as T r increases, the efficiency of full-duplex operation drops. It is worth noting that at low RSI power the DoF plays the most important role to have higher sum rate. For instance, consider Fig. 6(a), in which the cases FD={4,5,5,4}, FD={4,4,6,4} and FD={4,6,4,4} have DoF total = 4 -DoF total is the minimum of the DoF of source-relay and relay-destination channels, i.e., DoF total = min (DoF sr , DoF rd )while the cases FD={4,7,3,4} and FD={4,3,7,4} have DoF total = 3. At T r = 0, there is a noticable gap between the first three cases and the last two ones, while the difference of the first three cases from each other is small. The big gap is due to the difference in DoF total and the small one is due to the difference in SNR. Similarly, in Fig. 6(b) the three cases FD={10,12,12,10}, FD={10,10,14,10} and FD={10,10,14,10} with DoF total = 10 have higher rates than the two cases FD={10,6,18,10} and FD={10,18,6,10} with DoF total = 6. Finally, it can be seen in both Fig. 6(a) and Fig. 6(b) that at T r = 0 there is no difference between cases that have the same DoF total but different DoF sr and DoF rd . As it can be seen in both Fig. 6(a) and (b), for cases with K t > K r the sum rate drops quickly as RSI increases. In fact, the more relative antennas at the relay transmitter compared to its receiver, the faster the sum rate drops with the rise in RSI.
To understand this behaviour of the system better, again, consider case {10,18,6,10} and also, suppose T r → ∞. As discussed before, we have DoF sr = 10 and DoF rd = 6. Moreover, we have DoF I = 6 for the interference channel (H r ). Unlike the case with no interference, in this case the bottleneck is no longer the relay-destination link. This is due to the fact that interference can act to the detriment of some six of the source-relay subchannels. As we have DoF I = 6, interference can choose at most six independent subchannels and as we assumed T r → ∞, for those subchannels we get SINR → 0. Therefore no information can be conveyed from those links and the bottleneck becomes the source-relay link with 4 usable subchannels. It can be seen in Fig. 6(a) and (b) that as T r increases, the cases with the same sumrate at T r = 0 start to diverge because of the different characteristics of the interference they experience.
We explain the effect of interference in the following subsection in more detail.

B. The Impact of IRS
In this part, we evaluate the impact of IRS on the throughput rate when it is used to performe different tasks. Fig. 7 shows the throughput for three different scenarios, namely, when IRS is used to help the transmitter-relay link, when it is applied to cancel RSI and when the IRS job is to help the realy-destination link and then the results are compared with two cases where the system is working in HD wiht IRS and the case where the system is working in FD with no IRS. It is also assumed that Tr (HrH H r ) = 75% i. e. the system works at high RSI range. As it can be seen, the highst performance is achieved when the IRS is utilized to deal with the RSI. As a result, for the rest of the paper we use the IRS for this purpose. Also, it can be seen that when the number of     all. This is mainly due to the fact that, unlike the average case, for the case of worst case scenario, the number of IRS elements should be at least as large as the dimension ofH r , otherwise, the the IRS feasible set cannot span into all dimensions ofH r . Therefore, there is always least one representation forH r in which IRS cannot do any RSI cancellation. Also, comparing two figures Fig. 8(a) and Fig. 8(b), one can conclude that, when the dimension ofH r increases, the effort that IRS has to make in order to cancel RSI remarkably increases which is consistent with the previous statement.

C. Relay Tx/Rx Antenna allocation
Suppose that the relay has K t + K r = 8 in total. Furthermore, following cases in which the number of antenna at the source and destination are {N t , N r } = {4, 4}. The question is, from 8 antennas at the relay, how many should be used for reception for the robust design? Fig. 9 shows the sum rate as a function of K r for different values of T where there is no IRS and there is an IRS with M = 60 elements respectively.
As it can be seen, by using more antennas for reception than for transmission, i.e., K r > K t , at the relay, the throughput rate is maximized. This is due to the fact that, increasing the signal-to-noise ratio (SNR) of the source-relay streams enhances the overall throughput rate more than increasing the number of antennas for transmission in order to enhance the DoF of the relay-destination link. Furthermore, notice that in this setup the overall DoF from the source to destination is limited by the DoF of the source-relay link, i.e., the bottleneck is in the first hop.
By comparing two scenarios we see that having an IRS not only improves the rates in all cases, but also it mat change the best antenna allocation. For instance, for the case of T = 15%, it is best to have 6 antennas at the relay receiver and 4 antennas at the relay transmitter. However, after establishing the IRS, the best antenna allocation changes to 5 antennas at each end. For instance, the results show that although the DoF total for both {N t , K t , K r , N r } = {4, 3, 5, 4} and {N t , K t , K r , N r } = {4, 5, 3, 4} is 3, the sum rate capacity of the latter is much better than that of the former at high interference. This is because of the fact when DoF sr > DoF rd , the source relay link enjoys DoF sr − DoF rd subchannels with no interference. Therefore, the source can manage to obtain higher sum rate by choosing its power allocation wisely. However, in the case of DoF sr ≤ DoF rd , no matter how well the power allocation is done, all sub channels suffer from interference at the source-relay end.

D. Full-Duplex vs Half-Duplex
In this subsection, we determine the thresholds where the HD relaying outperforms the FD relaying.
This thresholds provide a mode-switching threshold in hybrid HD/FD relay systems. As it can be seen in Fig. 10, for each case of K r , there is a maximum value of T P above which the HD mode outperforms the FD mode in terms of sum rate maximization. Furthermore, Fig. 10 shows the threshold for different different IRS configurations. For this part we continued with the case of N t = 4, K r = 5, K t = 5, N r = 4.
As it can be seen, by increasing the number of antennas, the threshold occurs at higher RSI. This is in fact a direct result of getting better performance by having more antennas at the relay's receiver. It is worth noting that the IRS has a great impact on the performance of FD relaying. For instance, by having an IRS consists of only 60 elements, FD mode outperforms HD mode in almost all cases.

VI. CONCLUSION
In this paper, we investigated a multi-antenna source communicating with a multi-antenna destination through a multi-antenna relay. The relay is assumed to exploit a decode-and-forward (DF) strategy. An IRS is installed to hep the relay cope with the RSI. The transceivers are designed in order to be robust against  In contrast, points below the curve belong to cases where FD performs better than HD. the worst-case residual self-interference (RSI). To this end, the worst-case achievable throughput rate is maximized. This optimization problem turns out to be a non-convex problem. Assuming that the degreesof-freedom (DoF) of the source-relay link is less than the DoF of the relay-destination link, we determined the left and right matrices of the singular vectors of the worst-case RSI channel. Then, the problem is simplified to the optimal power allocation at the transmitters, which guarantees robustness against the worst-case RSI singular values. This simplified problem is still non-convex. Based on the intuitions for optimal power allocation at the source and relay, we proposed an efficient algorithm to capture a stationary point. Hence, in a DF relay with multi-stream beamforming, we determine the threshold where the halfduplex relaying outperforms the full-duplex relaying. This threshold provides a mode-switching threshold in hybrid half-duplex full-duplex relay systems.
Definition 1. For a vector a, we denote vector a ↓ which has the same components as a except that they are sorted in a decreasing order.
The vector a is said to be majorized by vector b and denoted by a ≺ b if: where a ↓ i is the i'th component of a ↓ , N r is the number of vector components and K ≤ N . If the last equality does not hold, a is said to be weakly majorized by b and denoted by a ≺ w b. Definition 3. The vector a is said to be multiplicatively majorized by vector b and denoted by a ≺ × b if: ). Adding all above inequalities, we get . (77) The above equation indicates R FD sr (σ 2 r ) > R FD sr (σ 2 r ) which contradicts the first assumption R FD sr (σ 2 r ) ≤ R FD sr (σ 2 r ). This completes the proof of the minimization part. For the maximization part, the general idea is the same. Again, the proof is by contradiction. We assume the optimal vector γ s , for which we have R FD sr (γ s ) ≥ R FD sr (γ s ), does not sum to P s . Therefore, we have γ s 1 < P s . Then there exists ε > 0 for which we have γ s 1 + ε = P s . Now we define where, η = i 1+σ 2 r i γr i σ 2 1 i + γ s i . Now we define the new source power allocation as below One can check that i γ s ρ(i) = P s and σ 2 1 i γ s ρ(i) ≥ σ 2 1 i+1 γ s ρ(i+1) . Thus, the new source power allocation is in the feasible set. Now the remaining is to make sure the new allocation does not change the corresponding σ 2 r . Using Lagrangian multiplier we have L = i log 2 1 + σ 2 1 i γ s ρ(i) 1 + γ r i σ 2 = i log 2 1 + σ 2 1 i (γ s ρ(i) + ε i ) 1 + γ r i σ 2 Now notice that as i log 2 (1 + ε η ) is a constant we have ∂ i log 2 (1+ ε η ) ∂σ 2 r i = 0 and ∂ i log 2 (1+ ε η ) ∂λ = 0. As a result, the optimum interference allocation for γ r is the same as that of γ r . Similarly to the case of minimization, here we have i ε i = ε. Also we have ε i ≥ 0 and there exist at least one i for which we have ε i > 0. Finally as f i (x) = log 1 + σ 2 1 i x 1+σ 2 r ρ(i) γr i is a monotonically increasing function of x, we conclude R FD sr (γ s ) < R FD sr (γ s ) which contradicts the first assumption of γ s being the optimal source power allocation, and the proof is complete. where and φ 2 . = min i −σ 2 1 i γ s ρ(i) γ r i 1 + σ 2 r ρ(i) γ r i 1 + σ 2 1 i γ s ρ(i) + σ 2 r ρ(i) γ r i respectively.
Next, we show g(P r ) = R FD sr (P r ) − R FD rd (P r ) is a monotonically decreasing function of P r . It is sufficient to show dR FD sr dPr ≤ 0 and dR FD rd dPr > 0. We have Now we define It is obvious that ψ 1 > 0 and ψ 2 ≥ 0. Now we have Finally, one can conclude dg dP r = dR FD sr dP r − dR FD rd dP r ≤ −ψ 2 − ψ 1 < 0.
One can check that for positive values x, y and a, we always have 2x+ay 2 √ x 2 +axy ≥ 1. As a result f is an increasing function of both x and y. The rest of the proof is as follows in which (a) holds because γ r i ≥ γ r i+1 and (b) holds because σ 2 1 i γ s ρ(i) ≥ σ 2 1 i+1 γ s ρ(i+1) .