A Probabilistic Reformulation Technique for Discrete RIS Optimization in Wireless Systems

TL;DR

This paper introduces a probabilistic reformulation method for discrete RIS optimization in wireless systems, transforming the problem into a continuous domain to improve scalability and performance.

Contribution

It develops a novel probabilistic technique that converts discrete RIS optimization into a continuous problem, enabling analytical and stochastic solutions with better performance.

Findings

The method outperforms traditional approaches like CPP and SDR.

It provides both stochastic sampling and analytical gradient descent solutions.

The approach is applicable to various discrete RIS optimization problems.

Abstract

The use of reconfigurable intelligent surfaces (RIS) can improve wireless communication by modifying the wireless link to create virtual line-of-sight links, bypass blockages, suppress interference, and enhance localization. However, enabling the RIS to modify the wireless channel requires careful optimization of the RIS phase-shifts. Although discrete RIS is more practical given hardware limitations, continuous RIS phase-shift optimization has attracted significantly more attention than discrete RIS optimization, which suffers from issues like quantization error and scalability. We develop a comprehensive probabilistic technique to transform discrete optimization problems into optimization problems of continuous domain probability parameters by interpreting the discrete optimization variable as a categorical random vector and computing expectations with respect to those parameters. We rigorously establish that for the unconstrained case, the optimal points of the reformulation and the original problem coincide. For the constrained case, we prove that the transformed problem is a relaxation of the original problem. We apply the proposed technique to two canonical discrete RIS applications: SINR maximization and overhead-aware rate and energy efficiency (EE) maximization. The reformulation enables both stochastic and analytical interpretations of the original problems. The former interpretation yields a stochastic sampling technique, whereas the latter yields an analytical gradient descent (GD) approach that employs closed-form approximations for the expectation. The numerical results demonstrate that the proposed technique is applicable to a variety of discrete RIS optimization problems and outperforms other general approaches, such as closest point projection (CPP) and semidefinite relaxation (SDR) methods.

Tables

Table 1. TABLE I : Average number of RIS elements to maximize EE in presence of interference with varying p 𝑝 p .

	0 dBm	10 dBm	20 dBm	30 dBm	40 dBm	50 dBm
UA ($T_0=0.2$ms)	$70.8$	$59.7$	$51.7$	$46.4$	$49.2$	$94.5$
UA ($T_0=1$ms)	$21.1$	$17.3$	$14.64$	$13$	$13.4$	$23.6$
SSA-T ($T_0=0.2$ms)	$52.7$	$47.4$	$46$	$46.1$	$47.9$	$40.7$
SSA-T ($T_0=1$ms)	$17.4$	$14.6$	$13.5$	$13.4$	$13.2$	$14.7$

Equations

$$\mathbb{P}(\mathbf{x}|\mathbf{P})=\prod _{i=1}^n\sum _{j=1}^b\delta (x_i-c_j)p_{i,j},p_{i,j}\in [0,1],\sum _{j=1}^b{p}_{i,j}=1,$$

$${\Omega }_{\mathbf{x}}\subseteq {\Omega }_{\mathbf{P}}.$$

$$\mathrm{Degen}(\mathbf{x})=\mathbf{P}=\left[\begin{array}{cc}0& 1\\ 1& 0\\ 0& 1\end{array}\right].$$

$$\begin{array}{r}\underset{\mathbf{x}}{\min }f(\mathbf{x})\le \xi (\mathbf{P})=\sum _{k=1}^{b^n}f(\mathbf{x}\{k\})\mathbb{P}(\mathbf{x}=\mathbf{x}\{k\}|\mathbf{P})\le \underset{\mathbf{x}}{\max }f(\mathbf{x}),\end{array}$$

$$\xi ({\mathbf{P}}_0)=\sum _{k=1}^{b^n}f(\mathbf{x}\{k\})\mathbb{P}(\mathbf{x}=\mathbf{x}\{k\}|{\mathbf{P}}_0)=f(\mathbf{x}\{k_0\})=p^{\ast }$$

$$⟹\sum _{k=1,k\ne k_0}^{b^n}\left(f(\mathbf{x}\{k\})-f(\mathbf{x}\{k_0\})\right)\mathbb{P}(\mathbf{x}=\mathbf{x}\{k\}|{\mathbf{P}}_0)=0.$$

$$d^{\ast }\le p_e^{\ast }\le p_c^{\ast }.$$

$$g_c(\mathbf{\lambda },\mathbf{v})=$$

$$=$$

$$g_c(\mathbf{\lambda },\mathbf{v})=$$

$$\sum _{j=1}^r{v}_j{\mathrm{E}}_{\mathbf{x}\sim \mathbb{P}(\mathbf{x}|\mathbf{P})}[h_j(\mathbf{x})],$$

$${\mathrm{E}}_{\mathbf{x}\sim \mathbb{P}(\mathbf{x}|{\mathbf{p}}_x)}\left[\sum _{j=1}^n\sum _{i=1}^n{x}_i{x}_j{W}_{ij}\right]$$

$$={\sum _{j=1}^n\sum _{i=1}^n}_{i\ne j}(2p_{x,i}-1)(2p_{x,j}-1)W_{ij}+\sum _{j=1}^n{W}_{ii},$$

$$\mathrm{E}[\mathbf{x}{\mathbf{x}}^T]=(\mathbf{y}{\mathbf{y}}^T)\odot {\mathbf{E}}_m+{\mathbf{I}}_N,$$

$${\mu }_{qf}(\mathbf{G},\mathbf{z},\mathbf{y})=$$

$${\vartheta }_{qf}(\mathbf{G},\mathbf{z},\mathbf{y})=$$

$${\mu }_{ql}(\mathbf{G},\mathbf{z},\mathbf{y})=\mathrm{E}[{\mathbf{x}}^T\mathbf{G}\mathbf{x}{\mathbf{z}}^T\mathbf{x}]=2{\mathbf{y}}^T{\mathbf{G}}_{wd}\mathbf{z}+{\mathbf{z}}^T\mathbf{y}\mathrm{Tr}(\mathbf{G})+$$

$$1^T\{({\mathbf{G}}_{wd}{\mathbf{Y}}_{wd})\odot {\mathbf{Y}}_{wd}\}(\mathbf{y}\odot \mathbf{z}),$$

$${\vartheta }_{ql}(\mathbf{G},\mathbf{z},\mathbf{y})=2{\mathbf{G}}_{wd}\mathbf{z}+\mathbf{z}\mathrm{Tr}(\mathbf{G})+(({\mathbf{G}}_T^T\odot {\mathbf{E}}_m)\mathbf{y})\odot \mathbf{z}+$$

$$\mathrm{diag}({\mathbf{G}}_T\mathrm{diag}(\mathbf{y}\odot \mathbf{z}){\mathbf{E}}_m)+({\mathbf{G}}_T\odot {\mathbf{E}}_m)(\mathbf{y}\odot \mathbf{z}),$$

$${\mu }_{qs}(\mathbf{G},\mathbf{y})=\mathrm{E}[({\mathbf{x}}^T\mathbf{G}\mathbf{x}{)}^2]={\mathbf{y}}^T\left({\mathbf{G}}_s-\mathbf{F}(\mathbf{y})\right)\mathbf{y}+\mathrm{Tr}(\mathbf{G}{)}^2+$$

$$2\mathrm{Tr}(\mathbf{Z})+({\mathbf{y}}^T\mathbf{G}\mathbf{y}{)}^2-{\mathbf{d}}^T{\mathbf{G}}_g\mathbf{d},$$

$$\mathbf{B}=\left[\begin{array}{c}{\mathbf{b}}_{1,1},\dots ,{\mathbf{b}}_{1,N}\\ \cdots ,\cdots ,\cdots ,\\ {\mathbf{b}}_{N,1},\dots ,{\mathbf{b}}_{N,N}.\end{array}\right],$$

$${\vartheta }_{qs}(\mathbf{G},\mathbf{y})=({\mathbf{G}}_s+{\mathbf{G}}_s^T)\mathbf{y}+2{\mathbf{y}}^T\mathbf{G}\mathbf{y}(\mathbf{G}+{\mathbf{G}}^T)\mathbf{y}-$$

$$2{\mathbf{y}}^T(\mathbf{G}+{\mathbf{G}}_{wd})\mathbf{y}(\mathrm{diag}(\mathbf{G})\odot \mathbf{y})-{\mathbf{d}}^T\mathrm{diag}(\mathbf{G})(\mathbf{G}+{\mathbf{G}}_{wd})\mathbf{y}$$

$$-\mathrm{diag}(\mathbf{G}{)}^T\mathbf{d}(\mathbf{G}+{\mathbf{G}}_{wd}{)}^T\mathbf{y}-2(({\mathbf{G}}_g+{\mathbf{G}}_g^T)\mathbf{d})\odot \mathbf{y}-$$

$$8\mathbf{y}\odot {\mathbf{b}}_s-4({\mathbf{U}}_{wd}+{\mathbf{U}}_{wd}^T)\mathbf{y},$$

$$y_r=(h_{d_0}+{\mathbf{h}}_0^H\mathrm{diag}(\mathbf{\theta }){\mathbf{f}}_0)x_{s,0}+$$

$$\sum _{i=1}^{N_I}(h_{d_i}+{\mathbf{h}}_i^H\mathrm{diag}(\mathbf{\theta }){\mathbf{f}}_i)x_{s,i}+w,$$

$$\gamma =\frac{{\beta }_0|h_{d_0}+{\mathbf{h}}_{c_0}^H\mathbf{\theta }{|}^2}{\sum _{i=1}^{N_I}{\beta }_i|h_{d_i}+{\mathbf{h}}_{c_i}^H\mathbf{\theta }{|}^2+{\sigma }_w^2}=\frac{f_s(\mathbf{\theta })}{f_I(\mathbf{\theta })}=\frac{{\mathbf{\theta }}^T{\mathbf{R}}_0\mathbf{\theta }+{\mathbf{c}}_0^T\mathbf{\theta }}{{\mathbf{\theta }}^T\mathbf{K}\mathbf{\theta }+{\mathbf{s}}^T\mathbf{\theta }},$$

$$\underset{\mathbf{\theta }\in \{-1,+1{\}}^N}{\min }$$

$$\underset{{\mathbf{y}}_s\in [-1,+1{]}^N}{\min }$$

$${\mathbb{P}}_B(\mathbf{\theta }|{\mathbf{p}}_s)=\prod _{n=1}^N\left(\delta ({\theta }_n-1)p_{s,n}+\delta ({\theta }_n+1)(1-p_{s,n})\right),$$

$${\mathcal{J}}_1({\mathbf{y}}_s)=\frac{\mathrm{E}[f_s(\mathbf{\theta })]}{\mathrm{E}[f_I(\mathbf{\theta })]}=\frac{{\mu }_{qf}({\mathbf{R}}_0,{\mathbf{c}}_0,{\mathbf{y}}_s)}{{\mu }_{qf}(\mathbf{K},\mathbf{s},{\mathbf{y}}_s)},$$

$${\mathcal{J}}_2({\mathbf{y}}_s)={\mathcal{J}}_1({\mathbf{y}}_s)-\frac{\mathrm{E}[f_s(\mathbf{\theta })f_I(\mathbf{\theta })]}{{\mathrm{E}}^2[f_I(\mathbf{\theta })]}+\frac{\mathrm{E}[f_I^2(\mathbf{\theta })]\mathrm{E}[f_s(\mathbf{\theta })]}{{\mathrm{E}}^3[f_I(\mathbf{\theta })]}.$$

$$c_v({\mathbf{y}}_s)=\mathrm{E}[f_s(\mathbf{\theta })f_I(\mathbf{\theta })]=\frac{{\mu }_{qs}({\mathbf{R}}_0+\mathbf{K},{\mathbf{y}}_s)-{\mu }_{qs}({\mathbf{R}}_0+\mathbf{K},{\mathbf{y}}_s)}{4}+{\mu }_{ql}({\mathbf{R}}_0,\mathbf{s},{\mathbf{y}}_s)+{\mu }_{ql}(\mathbf{K},{\mathbf{c}}_0,{\mathbf{y}}_s)+{\mathbf{c}}_0^T\left(({\mathbf{y}}_s{\mathbf{y}}_s^T)\odot {\mathbf{E}}_m+{\mathbf{I}}_N\right)\mathbf{s},$$

$${\vartheta }_{cv}={\nabla }_{{\mathbf{y}}_s}{c}_v({\mathbf{y}}_s)=\frac{{\vartheta }_{qs}({\mathbf{R}}_0+\mathbf{K},{\mathbf{y}}_s)-{\vartheta }_{qs}({\mathbf{R}}_0+\mathbf{K},{\mathbf{y}}_s)}{4}+{\vartheta }_{ql}({\mathbf{R}}_0,\mathbf{s},{\mathbf{y}}_s)+{\vartheta }_{ql}(\mathbf{K},{\mathbf{c}}_0,{\mathbf{y}}_s)+\mathbf{s}\odot ({\mathbf{E}}_m({\mathbf{c}}_0\odot {\mathbf{y}}_s))+$$

Full text

I Introduction

An RIS is a large array composed of low-cost reflecting elements, each of which can impart controllable phase-shifts to the incident signal, thereby modifying the propagation channel. However, because of hardware constraints, the phase-shift induced by each reflective element is normally limited to a set of discrete values. When configured appropriately, RISs can create multiple virtual LoS links [2], improve channel rank [3], transform a fast-fading channel to a slow-fading one [4], suppress co-channel interference [5], enhance localization performance [6, 7, 8], etc. However, optimizing the RIS phase-shifts is the first step to reaping these benefits. Despite the fact that RIS optimization has been extensively investigated in the literature, the majority of its attention has been directed toward the scenario of continuous phase-shifts as this scenario allows easier insights and upper-bounds on the performance of a wireless network. As a consequence, the discrete RIS optimization techniques often appear as an afterthought and the existing techniques that deal with discrete RIS optimizations suffer from various issues, such as scalability and arbitrarily bad performance due to quantization error.

Motivated by the scarcity of scalable and reliable discrete RIS optimization techniques [9], stochastic interpretation of semidefinite relaxation technique [10], and recent efforts to approach binary optimization problems [11] with a lens of probability, we develop a comprehensive technique to transform optimization problems of discrete variables into optimization problems of continuous domain probability parameters. We also rigorously prove that in terms of optimal points, the transformed problem is mathematically equivalent in the unconstrained case and a relaxation with respect to the original problem in the constrained case. Moreover, we gain further insights into our reformulation by investigating the simple two-way partitioning problem and report several moment and gradient results for quadratic forms in binary optimization problems. Ultimately, we apply this reformulation in two different canonical discrete RIS optimization problems demonstrating both the stochastic and analytical approaches. The numerical results confirm that the expectation-based algorithms outperform the conventional approaches. Note that, even though the proposed reformulation is inspired by discrete RISs, the scope of the reformulation is more general and could potentially find applications in other domains as well.

I-A Related Work and Motivation

Although the study and design of discrete RIS phase-shifts have sparked some interest recently, a large portion of the literature relaxes the discrete constraint to a continuous one, solves the approximate problem, and then quantizes the solution to the closest discrete point. This two-fold approximation is shown to provide arbitrarily bad solutions in the worst-case scenario [9]. Yet, continuous RIS optimization remains a big part of the discrete RIS optimization literature. In light of this, some of the most used optimization strategies for continuous RIS phase-shits are combinations of a) SDR, b) minorize-maximization (MM) algorithm, c) penalty methods, d) manifold optimization, e) alternating direction method of multipliers (ADMM), and f) treating phase-shifts as optimization variables instead of the complex gains they provide.

In [12, 2], the authors jointly optimized the active beamforming vector and RIS-based passive beamforming vector in multiple-antenna systems employing SDR. The authors of [13] utilized RISs to enhance the physical layer security of a multiuser multiple-input-single-output (MU-MISO) wireless system. In particular, they used a combination of a penalty-based approach, SDR, and successive convex approximation (SCA) to address the unit modulus constraint of RIS phase-shifts. Energy-efficient RIS designs for a MU-MISO wireless system are developed in [14]. In this paper, the authors developed two algorithms to maximize energy efficiency of the network. One of them uses gradient descent method whereas the other uses MM algorithm while considering unit modulus constraint and a realistic power consumption model. Similarly, in [15], the authors leveraged the MM algorithm and complex circle manifold (CCM) method to propose two algorithms that maximizes weighted sum-rate of a multicell MIMO network. In an RIS-assisted backscatter system, the RIS is optimized with a combination of SDR and ADMM technique in [16]. In [17], weighted sum-rate is maximized in an RIS-aided cell-free network through ADMM. The authors of [18] tackled the resource allocation problem in an RIS-assisted wireless network. They developed an algorithm that uses the SCA and penalty method to jointly optimize phase-shifts and on-off status of RISs to maximize energy efficiency under a total power constraint. The authors of [19] enhanced the physical layer security by optimizing RIS phase-shifts with fractional programming and manifold optimization techniques. Using a similar technique, the authors of [20] developed an algorithm combining both SDR and manifold optimization techniques to optimize an RIS-aided edge caching system. Recently, the authors of this paper treated the vector of phase-shifts itself as an optimization variable instead of the vector of the complex gains they provide and optimized the RIS phase-shifts with the GD method to maximize the SINR [21] similar to the methodology used in [22].

In the realm of discrete RIS optimization literature, the optimization tools often used are exhaustive search, CPP from a continuous relaxation, and branch-and-bound (BB) methods. In [23], the authors investigated a practical discrete RIS-aided wideband orthogonal frequency division multiplexing (OFDM) system and optimized the discrete RIS element-wise exhaustive search in an alternating optimization framework. The authors of [24, 25] optimized the RIS phase-shifts in a MISO wireless network using BB methods that scale exponentially with the number of RIS elements. An RIS-aided MIMO system with low-resolution digital-to-analog converters (DACs) is jointly optimized with particle swarm optimization (PSO) algorithm in [26]. This algorithm is shown to work with both continuous and discrete RIS phase-shifts. For maximizing the achievable rate in a MIMO system, the authors in [27] relaxed the discrete RIS constraints to continuous ones and used a projected gradient method (PGM) to solve the problem. Similar continuous relaxation and CPP method was conducted in [28]. Specifically, the authors in [28] investigated the spectral efficiency and energy efficiency trade-off in a MU-MIMO setup. In that paper, the discrete RIS phase-shift constraint was relaxed to a continuous one and then solved by MM and accelerated gradient method. The mentioned optimization problems are either not scalable or suffer from arbitrarily bad performance due to the mentioned two-fold approximation. However, there are some recent efforts [29, 9, 30] that provide scalable optimal discrete RIS beamforming optimization for single-input-single-output systems as the resulting objective function has a low-rank matrix with $d\leq 2$ based on the fixed-rank result of [31]. These strategies suffer from being too specialized for single-antenna scenarios and do not work in multi-antenna scenarios that are more practical.

Going beyond the RIS literature, there has been some interesting approaches to binary quadratic optimization problems [10, 32, 33, 11, 34]. The authors of [10] show that SDR formulation is actually a stochastic version of the original non-convex quadratic program. Even when the non-convexity originates from binary variables, the stochastic interpretation still works. The authors of [32, 33] approach binary quadratic programs with a probabilistic data association (PDA) algorithm that treats the optimization variable as a binary random variable and iteratively updates the probabilities with Gaussian noise approximation. This is shown to achieve near-optimal results. Recently, the authors of [11] provide a stochastic gradient descent framework for binary optimization problems by similarly treating the optimization variable as a random variable and then taking expectation on it. Inspired by these probabilistic optimization techniques along with the lack of novel generalized discrete RIS optimization techniques, we develop a comprehensive probabilistic technique to transform discrete optimization problems that opens up new avenues to approach these problems in the continuous probability parameter domain. This technique is then used in discrete RIS cases to showcase its general nature and effectiveness.

I-B Contributions

We approach the general discrete optimization problems with a different perspective of probability. This results in a comprehensive probabilistic reformulation technique with a wide applicability, including to the discrete RIS problems, which was our original inspiration behind this work. Our key contributions in this paper are listed next.

I-B1 A comprehensive probabilistic technique for general discrete optimization problems

We develop a comprehensive probabilistic technique to reformulate general discrete optimization problems (that are not limited to binary programs) into continuous domain problems. In particular, we re-imagine the entries of the optimization variable as independently but not identically distributed (i.n.i.d) categorical random variables and replace the objective function and constraints, if any, with their expectations. We rigorously establish the equivalence between a general unconstrained problem with a unique optimal solution and the reformulated problem in terms of the optimal point. Additionally, when the original problem is constrained, we prove that the primal solution of the transformed problem is bounded between the dual and primal solution of the original problem. We also show that when strong duality holds, the transformed problem has the same optimal objective value as the original problem. Utilizing this technique, random sampling from a non-degenerate probability parameter solution can provide a better solution with the number of samples similar to Gaussian randomization in SDR.

I-B2 Derivation of various analytical moments and their gradients associated with the quadratic form and binary random vectors

As discrete RIS problems deal with binary phase-shifts often, using our reformulation technique naturally gives rise to expectations associated with the quadratic form and the binary random vectors. For example, both the denominator and numerator of the SINR or secrecy rate often contain quadratic forms [22]. The quadratic form is a canonical construct that appears in the wireless literature frequently. For a gradient-based optimization approach, the gradients of these expectations will also be required. For this reason, we also derive the first and second moments of the said quadratic forms along with their gradients. These key intermediate results are later used in one of our algorithms demonstrating their importance.

I-B3 GD algorithm for the SINR maximization

As the first canonical case study, we apply this technique to an SINR maximization problem and propose a stochastic GD and an analytical GD approach to solve the reformulated problem. We derive and use the first and second-order Taylor approximations of the expectation of the SINR in the analytical GD algorithm while an estimator of the gradient is used in the stochastic approach. The expectation-based algorithms are shown to perform better than the conventional practical approaches evaluated.

I-B4 Stochastic sampling approach for ternary random vectors for EE and rate maximization

We also apply this technique to our second case study, an overhead-aware rate and EE maximization problem which leads to expectations associated with a ternary random vector. As deriving the analytical expectation was challenging for this specific case study, we develop a stochastic sampling approach for such a ternary random vector where the gradient is estimated with Monte Carlo (MC) samples, thereby demonstrating the versatility of the proposed approach. We demonstrate that this framework is well-suited for non-smooth objective forms and performs well in both interference-free and interference-rich scenarios. Moreover, the developed stochastic approach is demonstrated to work with different objective functions like rate and EE without the need for changing the algorithm.

I-B5 Computational complexity discussion

We have also derived worst-case computational complexities with big-O notation for all the proposed algorithms.

Notations

The distribution of a standard complex normal random variable is denoted by $\mathcal{CN}(0,1)$. The matrix, scalar and vector entities are denoted by $\bf X$, $x$, and $\bf x$, respectively. All the vectors are column vectors unless defined explicitly. For a vector $\bf x$, ${\rm diag}\left({\bf x}\right)$ denotes a diagonal matrix with the entries of $\bf x$ as its diagonal elements. For a matrix $\bf X$, ${\bf X}^{H}$, ${\bf X}^{T}$, ${\rm Re}\left({\bf X}\right)$, ${\rm Tr}\left({\bf X}\right)$, ${\rm diag}({\bf X})$, and ${\bf X}\succeq 0$ denote its conjugate transpose, transpose, real part, trace, diagonal elements as a vector, and positive semidefiniteness, respectively. Additionally, ${\bf X}_{wd}={\bf X}-{\rm diag}({\bf X})$. The expectation operation is denoted by $\rm E[\cdot]$, ${\rm var}(\cdot)$ denotes a total variance operator which evaluates the trace of the variance-covariance matrix of the random vector argument, and the operator $\odot$ denotes element-wise multiplication between two matrices. The L0 and L2 norm are denoted by $\|\cdot\|_{0}$ and $\|\cdot\|_{2}$, respectively. The identity matrix and all-one column vector of dimension $N$ are denoted by ${\bf I}_{N}$ and ${\bf 1}_{N}$, respectively.

II Probabilistic Reformulation for Discrete Optimization

II-A The Case of Unconstrained Discrete Optimization Problem

We begin with a general unconstrained discrete optimization problem where we make no assumptions about the objective function’s convexity. The optimization variable is a vector of length $n$ and each of the entry can take a discrete value among the set $\mathcal{C}=\{c_{1},c_{2},\ldots,c_{b}\}$.

{mini}

—s— x ∈C^nf(x).

Our main goal is to reformulate the problem in a form that does not deal with the discrete domain and shares the optimal solution with the original problem. To that end, we propose to re-imagine entries of $\bf x$ as i.n.i.d categorical random variables with the following joint probability density function (PDF):

[TABLE]

where the $(i,j)$-th entry of the matrix $\bf P$ is denoted by $p_{i,j}$, the $i$-th entry of $\bf x$ is denoted by $x_{i}$, and $\delta(\cdot)$ is the Dirac delta function. We then reformulate the original problem into a stochastic optimization problem: {mini}—s— p_i,j ∈Fξ(P)=E_x∼P(x—P)[f(x)],

where $\mathcal{F}$ is the set of possible $p_{i,j}$’s defined by (1). The connection between (II-A) and (II-A) and their solution sets are summarised in the following lemma.

Lemma 1**.**

The solution sets of the problems (II-A) and (II-A) are denoted by $\Omega_{\bf x}$ and $\Omega_{\bf P}$ and,

[TABLE]

Moreover if the unique optimal solution of (II-A) is ${\bf x}_{\rm opt}$, then ${\bf P}_{\rm opt}={\rm Degen}({\bf x}_{\rm opt})$ is the unique optimal solution of (II-A), where the ${\bf P}={\rm Degen}({\bf x})$ operation implies that the $(i,j)$-th entry of $\bf P$ is defined as $p_{i,j}=1$ only when $x_{i}=c_{j}$ while all the other entries are zero.

Proof:

We observe that $\Omega_{\bf x}$ has $b^{n}$ elements and each of them corresponds to one of the possible $b^{n}$ combinations that ${\bf x}$ can take. In (II-A), the same objective values can be attained by the corresponding ${\bf P}={\rm Degen}({\bf x})$ which is the parameter matrix of $n$ degenerate categorical distributions. Let us illustrate this with an example. Suppose we have a vector ${\bf x}=[-1,1,-1]^{T}$. Each element of ${\bf x}$ can adopt either $c_{1}=1$ or $c_{2}=-1$. By referring to the definition of the ${\rm Degen}(\cdot)$ function given in Lemma 1, the resulting parameter matrix for this vector is:

[TABLE]

The interpretation of this matrix is as follows: the first element of ${\bf x}$ is $-1$ with a probability of $1$, the second element is $1$ with a probability of $1$, and so forth. The parameter matrix thereby represents the degenerate distributions of the elements of ${\bf x}$ if it was reimagined as a random vector. Importantly, this parameter matrix can be translated back into the original ${\bf x}$ vector. Consequently, for every possible permutation of ${\bf x}$ (as given in equation (II-A)), there exists a corresponding parameter matrix in equation (II-A) that provides the same objective value. Upon reflecting on this, it becomes evident that the set $\Omega_{\bf x}$ must be included in $\Omega_{\bf P}$. In mathematical terms, this relationship can be represented as $\Omega_{\bf x}\subseteq\Omega_{\bf P}$.

For any feasible $\bf P$, it can be shown that,

[TABLE]

where ${\bf x}\{k\}$ denotes the $k$-th combination out of possible $b^{n}$ combinations of $\bf x$. This stems from the observations that the expectation is nothing but a convex combination of all the possible values of $f({\bf x})$. This is possible because the probability terms, which are always nonnegative, sum up to one, enabling the expression of the expectation as this sum. Such a convex combination of scalar values essentially represents a probability-weighted average. Each scalar is weighed by its corresponding probability or chance of occurrence. These probabilities fundamentally dictate the placement of the weighted average on the line between the minimum and maximum scalar values. Due to the constraint of the probabilities adding up to one, this average cannot exist outside this range. For instance, when the probabilities tend to favor larger scalar values, the resulting combination leans closer toward the maximum and vice versa. It is crucial to clarify that equating the expectation to a convex combination does not imply that the expectation is a convex function. We are not discussing the convexity of the expectation itself. Yet, due to the inherent properties of convex combinations, inequality (2) consistently holds, as indicated in [35]. Additionally, Carathéodory’s theorem offers further proof of this fact [11, 35].

Now assume that ${\bf x}_{\rm opt}$ is the unique optimal solution of $\eqref{eq:ori}$. It follows that, ${\bf P}_{\rm opt}=\rm{Degen}({\bf x}_{\rm opt})$ is an optimal solution of (II-A). Consider that $\exists{\bf P}_{0}\neq{\bf P}_{\rm opt}$, such that, $\xi({\bf P}_{0})=\xi({\bf P}_{\rm opt})=f({\bf x}_{\rm opt})$. The parameter matrix $\bf P$ cannot denote $n$ degenerate categorical distributions as the corresponding ${\bf x}_{0}={\rm Degen}^{-1}({\bf P}_{0})$ would violate the uniqueness assumption on ${\bf x}_{\rm opt}$. We then consider the non-degenerate distribution case. As the optimal value $p^{*}$ is shown to be the same for both of these problems, we can assume that $p^{*}=f({\bf x}\{k_{0}\})$ without any loss of generality. Then,

[TABLE]

As for some $k$, the value $f({\bf x}\{k\})$ needs to be equal to $f({\bf x}\{k_{0}\})$ for (4) to be true, this would also violate the uniqueness assumption on $k_{0}$. ∎

II-B The Case of Constrained Discrete Optimization Problem

Next, we explore whether such a coincident optimal solution through such a reformulation is valid for constrained problems as well. To that end, we write a general optimization problem with constraints without assuming convexity below: {mini}—s— x ∈C^nf_0(x), \addConstraintf_i(x)≤0  ∀i=1, 2, …, m \addConstrainth_j(x)= 0  ∀j=1, 2, …, r,

where the optimal value and solution set of this problem are denoted by $p_{c}^{*}$, and $\Psi_{\bf x}$, respectively. The transformed formulation is expressed as: {mini}—s— p_i,j ∈FE_x∼P(x—P)[f_0(x)], \addConstraintE_x∼P(x—P)[f_i(x)]≤0  ∀i=1, 2, …, m \addConstraintE_x∼P(x—P)[h_j(x)]= 0  ∀j=1, 2, …, r,

with optimal value $p_{e}^{*}$, and solution set $\Psi_{\bf P}$.

Lemma 2**.**

The original problem (II-B) and the transformed problem $\eqref{eq:eqns}$ share the same dual problem with the dual solution $d^{*}$. Moreover,

[TABLE]

Proof:

Similar to the proof of Lemma 1, it can be readily seen that for every feasible $\bf x$ in (II-B), there is a corresponding feasible parameter matrix ${\bf P}={\rm Degen}({\bf x})$ in (II-B). It directly follows from this observation with similar reasoning in the previous proof that the solution set $\Psi_{\bf x}$ is a subset of $\Psi_{\bf P}$. Consequently, the transformed optimization problem can be seen as a relaxation of the original constrained problem (II-B) and provides a better optimal value. Thus it can be established that $p_{e}^{*}\leq p_{c}^{*}$.

Next, we investigate the dual function of the original problem by expressing it as the infimum of the Lagrangian [36, 37]:

[TABLE]

where $\lambda_{i}$ is the $i$-th entry of $\boldsymbol{\lambda}$ and $v_{j}$ is the $j$-th entry of $\bf v$. The vectors discussed here are the dual variables related to our problem. We can think of the dual function as a softened form of equation (II-B), which consists of more stringent or hard constraints [36]. Crucially, for all non-negative vectors $\boldsymbol{\lambda}$, the dual function serves as a consistent lower bound for the optimal value of the primal problem, denoted $p_{c}^{*}$. For a more in-depth treatment of this well-known result, readers can refer to [36]. Next, using Lemma 1, we can reformulate the above dual function into the following expression:

[TABLE]

where the optimal dual solution after maximizing the concave dual function is denoted by $d^{*}$. We note that the dual function of (II-B) is equivalent to (7). Given the lower bound characteristic of the dual function, it is logical to conclude $d^{*}\leq p_{e}^{*}$. Moreover, strong duality ensures equality. Compiling these inequalities, we deduce $d^{*}\leq p_{e}^{*}\leq p_{c}^{*},$ which consequently proves the Lemma. It further implies that the relaxation (II-B) is non-trivial and it is bounded by $d^{*}$, given that the dual solution is bounded. ∎

II-C Discussion on the Two-way Partitioning Example

In this subsection, we will focus on the simple two-way partitioning problem to demonstrate the technique. We focus on this foundational example to facilitate a more comprehensive understanding and to draw parallels with other probabilistic methods more effortlessly. While the selected problem covers a wide array of applications, such as binary phase beamforming in an RIS-aided network [9], we delve into two more complex applications in Section III. We begin with the description of the two-way partitioning problem below: {maxi}—s— x ∈{-1,1}^nx^TWx,

where ${\bf W}\in\mathbb{R}^{n\times n}$ is a symmetric matrix. Next, we derive our reformulation of (II-C) based on Lemma 1 below starting with the following result:

[TABLE]

where $\bf p$ is the vector of parameters with $p_{x,i}=\mathbb{P}[x_{i}=1]$ denoting the $i$-th entry, $x_{i}$ denotes the $i$-th entry of $\bf x$, and $W_{ij}$ denotes the $i,j$-th entry of $\bf W$. This result directly follows from the facts that ${\rm E}[x_{i}]=2p_{x,i}-1$, ${\rm E}[x_{i}^{2}]=1$, and the entries are i.n.i.d.

Remark 1**.**

Note that, if all the entries of $\bf y$ are either $+1$ or $-1$, $\bf y=\rm{Degen}^{-1}({\bf p}_{x})$. In other words, in that case, $\bf y$ is a feasible $\bf x$ and vice versa.

With this result, the transformed problem is as follows: {maxi}—s— y ∈[-1,1]^ny^TW_wdy,

where ${\bf y}=2{\bf p}_{x}-{\bf 1}$, and ${\bf W}_{wd}$ is the matrix ${\bf W}$ with its diagonal elements set to zero. Note that, we effectively converted a binary quadratic program (BQP) to a non-convex box-constrained quadratic program (BoxQP) emphasizing the ability of our reformulation to change the structure of a problem while being equivalent in terms of optimal point. However, this is a known result in the optimization community [38, 39] and our reformulation provides a probabilistic proof. Other than this structural change, we can obtain more insights about our reformulation by focusing on the SDR of the original problem which can be derived by considering the following stochastic program by taking expectation on the objective value and the domain of (II-C) [10]: {maxi}—s— X ⪰0E_ζ∼N(0,X)[ζ^TWζ], \addConstraintE_ζ∼N(0,X)[ζ⊙ζ]=1_n,

where ${\bf X}$ is an arbitrary symmetric positive semidefinite matrix, and $\boldsymbol{\zeta}$ is a random vector drawn from a normal distribution with zero mean and covariance $\bf X$. Through the simple observation ${\rm E}_{{\boldsymbol{\zeta}}\sim\mathcal{N}({\bf 0},{\bf X})}\left[{\boldsymbol{\zeta}}{\boldsymbol{\zeta}}^{T}\right]={\bf X}$, this is equivalent to the classic SDR problem described below: {maxi}—s— XTr(WX) \addConstraintX⪰0 \addConstraintX_ii=1,  i=1,2,…,n.

In the above formulation, the addition of the rank-one constraint ${\rm rank}({\bf X})=1$ would make the problem equivalent to (II-C). However, this relaxed formulation is solvable in polynomial time unlike (II-C). SDR can also be seen as a relaxation of the original problem when we allow $x_{i}$ to be a multidimensional vector with a unit norm. These vectors can be found from The Cholesky decomposition of the solution of SDR. If the angle between two such vectors is really small, that implies which those two entries of $\bf x$ are more likely to be in the same group [40]. In other words, the SDR provides us with pairwise correlation information. In contrast, a relaxed version of our reformulation (1) will provide us the probabilities with which each entry of the original vector $\bf x$ will be $+1$ or $-1$.

The findings of this subsection reveal that our proposed reformulation has the potential to alter the structure of optimization problems. Even if the change is trivial in the case of simple objective forms, it is expected that more complex objective forms will yield non-trivial changes, which can have a significant impact on the efficiency and effectiveness of the optimization process. We will demonstrate these non-trivial structural changes in canonical case studies related to discrete RIS optimization in Section III. Furthermore, the results indicate that although our reformulation differs from SDR in terms of the information it provides, they share similarities in the formulation from a stochastic standpoint. These results encourage further exploration of our approach in more complex objective forms.

II-D Some Useful Results for Quadratic Expressions for Binary Random Vectors

As most discrete RIS applications deal with binary phase-shift RIS $\{-1,+1\}$ due to its simplicity in operation and implementation, it is only appropriate to derive the analytical moments and their gradients associated with the binary random vectors defined in (1). They can be used in different optimization contexts with such expectation-based formulations. The higher moment results are motivated by the previous subsection and will be heavily used in the next section. We begin with the covariance matrix next.

Remark 2**.**

For a random vector ${\bf x}\in\{-1,+1\}^{n}$ with i.n.i.d entries and expectation $\mathrm{E}[{\bf x}]={\bf y}$, the covariance matrix is

[TABLE]

where ${\bf E}_{m}$ is the all-one matrix with a hollow diagonal and ${\bf p}$ is defined similarly to (21).

Now, we state the first moment and its gradient in Lemma 3 without proof due to its trivial nature and partial proof in (8).

Lemma 3**.**

For a random vector ${\bf x}\in\{-1,+1\}^{n}$ with i.n.i.d entries and expectation $\mathrm{E}[{\bf x}]={\bf y}$, the expectation and the gradient of a sum between a quadratic form and a linear form are

[TABLE]

where $\bf G$ is a real symmetric matrix.

Next, we derive an expectation that is very important for covariance calculations between a quadratic form and a linear form in the next theorem.

Theorem 1**.**

For a random vector ${\bf x}\in\{-1,+1\}^{n}$ with i.n.i.d entries and expectation $\mathrm{E}[{\bf x}]={\bf y}$, the expectation of a product between a quadratic form and a linear form is

[TABLE]

where $\bf G$ is a real symmetric matrix and ${\bf Y}={\bf y}{\bf 1}^{T}$.

Proof:

See Appendix -A. ∎

We just state the gradient of the above expectation without proof in Corollary 1.

Corollary 1**.**

The gradient of the derived expectation in Theorem 1 can be calculated as:

[TABLE]

where ${\bf G}_{T}={\bf G}_{wd}{\bf T}_{0}$, and ${\bf T}_{0}={\rm diag}({\bf y}){\bf E}_{m}$.

Next, we focus on the second moment of a quadratic form in Theorem 2.

Theorem 2**.**

For a random vector ${\bf x}\in\{-1,+1\}^{n}$ with i.n.i.d entries and expectation $\mathrm{E}[{\bf x}]={\bf y}$, the second moment of a quadratic form is

[TABLE]

where $\bf G$ is a real symmetric matrix, ${\bf d}={\bf y}\odot{\bf y}$, ${\bf G}_{s}=2{\rm Tr}({\bf G}){\bf G}_{wd}+4{\bf Z}_{wd}$, ${\bf Z}={\bf G}_{wd}{\bf G}_{wd}^{T}$, ${\bf F}({\bf y})=({\bf y}\odot{\bf y})^{T}{\rm diag}({\bf G})({\bf G}+{\bf G}_{wd})+4{\bf U}_{wd}$, ${\bf U}=[{\bf I}_{N}\otimes({\bf y}\odot{\bf y})^{T}]{\bf B}$, and ${\bf G}_{g}=2{\bf G}_{wd}\odot{\bf G}_{wd}$. The matrix ${\bf B}$ is defined through blocks as

[TABLE]

where the $i$-th element of ${\bf b}_{k,j}$ is ${\bf b}_{k,j}^{i}=G_{{wd}_{ij}}G_{{wd}_{ki}}$.

Proof:

See Appendix -B. ∎

Now, we derive the gradient of the second moment in the Corollary 2.

Corollary 2**.**

The gradient of the derived expectation in Theorem 2 can be calculated as:

[TABLE]

where ${\bf d}={\bf y}\odot{\bf y}$, and $i$-th entry of ${\bf b}_{s}$ is ${\bf y}^{T}{\bf B}_{t}[i]{\bf y}-{\rm Tr}({\bf B}_{t}[i])$. The matrix ${\bf B}_{t}[i]$ can be derived by multiplicating the $i$-th column of ${\bf G}_{wd}$ with the $i$-th row of ${\bf G}_{wd}$.

Proof:

See Appendix -C. ∎

III Applications of The Proposed Reformulation

In a MIMO communication scenario, optimizing RIS phase shifts can be a challenging task, particularly when dealing with discrete RISs. Discrete RIS problems are generally more difficult to solve, making it necessary to split the original problem into smaller sub-problems that can be handled separately. Therefore, we focus on the canonical forms of discrete RIS sub-problems that frequently appear in the literature. For a unified treatment, we have chosen a signal model capable of representing a range of RIS-aided scenarios and sub-problems, including a device-to-device communication link, a cellular network where each antenna serves a different user through antenna selection, and a wireless communication link with interferers while the receive beamformer vector remains fixed [21]. This signal model can be expressed in the following point-to-point representation:

[TABLE]

where $y_{r}$ is the received signal from the Tx of interest (denoted by $i=0$), $h_{d_{i}}$ denotes the direct channel between the $i$-th Tx and Rx, ${\bf h}_{i}$ is the Tx-RIS channel, ${\bf f}_{i}$ denotes the RIS-Rx channel, $x_{s,i}$ is the data for the $i$-th Tx, ${\rm E}[x_{s,i}^{2}]=\beta_{i}$, $\boldsymbol{\theta}$ is the $N$-element discrete RIS phase configuration vector, $N_{I}$ is the number of interferers, and $w$ is the additive noise. Note that the users and interferers always transmit at their maximum power $p$. For a general MIMO communication scenario, these channels can be seen as the actual channels pre-multiplied and post-multiplied by precoding and receiver beamformer vectors, respectively. With such a versatile signal model, two use cases for RIS-assisted wireless communication systems are explored, namely, SINR maximization and overhead-aware RIS optimization. Note that, in both cases, RISs are assumed to be controlled by the receiver through an RIS controller [41].

III-A SINR Maximization with RIS Optimization

III-A1 System model

We consider a generic system model dictated by the signal model (17). We consider that the RIS phase vector $\boldsymbol{\theta}=[\theta_{1}\quad\theta_{2}\ldots\theta_{n}\ldots\theta_{N}]^{T}$ and $\theta_{n}\in\{-1,+1\}$. For ease of notation, we also define ${\bf h}_{c_{i}}=\left({\bf h}_{i}^{H}{\rm diag}({\bf f}_{i})\right)^{H}$ With this discrete RIS, the SINR can be expressed as,

[TABLE]

where ${\bf R}_{i}=\beta_{i}{\rm Re}\left({\bf h}_{c_{i}}{\bf h}_{c_{i}}^{H}+\frac{|h_{d_{i}}|^{2}}{N}{\bf I}_{N}\right)$, ${\bf K}=\sum\limits_{i=1}^{N_{I}}{\bf R}_{i}+\frac{\sigma_{w}^{2}}{N}\mathbf{I}_{N}$, $\sigma_{w}^{2}$ is the variance of the additive Gaussian noise, ${\bf c}_{i}=2\beta_{i}{\rm Re}({\rm conj}(h_{d_{i}}{\bf h}_{c_{i}}))$, and ${\bf s}=\sum\limits_{i=1}^{N_{I}}{\bf c}_{i}$.

III-A2 RIS optimization

In this subsection, our objective is to maximize the SINR given in (18) while the RIS elements are discrete in nature. The optimization problem is described below:

[TABLE]

As the domain of this problem is discrete and the problem is a fractional quadratic program, a common way to solve this problem is to relax the discrete domain and then project the solution to the closest discrete point. The relaxed version is solved through GD in [21]. Note that, we also consider SDR in the simulation results. We approach this problem with our reformulation (II-A) and transform this problem into a continuous domain problem. The reformulated problem is as follows:

[TABLE]

where ${{\bf y}_{s}}=2{\bf p}_{s}-1$ and ${\boldsymbol{\theta}}$ is assumed to be distributed with the joint PDF

[TABLE]

where $p_{s,n}\in[0,1]$ is the $n$-th entry of ${\bf p}_{s}$ and $\theta_{n}\in\{-1,+1\}$. We propose two approaches to solve (20): a) stochastic sampling approach (SSA), and b) analytical gradient descent approach. The former approach generally does not require an explicit expression of the gradient whereas the latter does. In SSA, a typical gradient estimator, which is based on the log-derivative trick and Monte Carlo sampling as described in [42], is often used in the GD algorithm. Following that, we have formulated an approach, SSA-B, which is a special case of SSA for binary variables using the same log-derivative trick and Monte Carlo sampling. We have omitted the details of SSA-B here, as this variant of SSA for binary variables has already appeared in a different context - the Bayesian optimal design of experiments, in [11]. For a detailed understanding of SSA-B, readers can refer to [11, Algorithm 3.1]. Next, we detail the analytical optimization approach for this case study. We explicitly develop the SSA for a non-binary random vector in the next case study, which is a direct result of the general probabilistic reformulation technique that we rigorously devised in Section II.

In the analytical gradient descent approach, calculating the direct expectation of a ratio of correlated random variables is difficult. So, we consider the Taylor series approximations of such an expectation [43]. Both the first-order approximation $\mathcal{J}_{1}({\bf y}_{s})$ and second-order approximation $\mathcal{J}_{2}({\bf y}_{s})$ are stated below:

[TABLE]

The second-order approximation requires two additional expectations that are derived along with their gradients in (23).

Using the definitions in (23), we can express the gradients of the Taylor series approximations as follows:

[TABLE]

Note that, they are stated without proof as they can be derived trivially with the basic chain rule. Armed with these gradients, we can develop simple update rules of a projected GD algorithm next:

[TABLE]

where $\mathbf{y}_{s}^{(t)}=2\mathbf{p}_{s}^{(t)}-{\bf 1}$ is the transformed probability vector at the $t$-th iteration, $\beta^{(t)}$ is the step-size and $\nabla_{\bf y}\mathcal{J}_{l}({\bf y}_{s}^{(t)})$ is the gradient of the $l$-th order Taylor approximation of the true expectation where $l\in\{1,2\}$. The steps (26) and the (27) are considered gradient step and projection step, respectively. For our box constraints, the projection turns out to be clipping the vector $\mathbf{y}_{s}^{(t+\frac{1}{2})}$ to $-1$ and $+1$. We also use Armijo-Goldstein (AG) line search [44] to find a good step-size while avoiding saddle points due to its diminishing nature [45]. Complete details of the GD approach are shown in Algorithm 1. Note that, according to the Remark 1, a feasible discrete $\boldsymbol{\theta}$ is also a feasible $\bf x$ and corresponds to the degenerate PDF itself that generates $\boldsymbol{\theta}$. So, we find the vector that aligns the phases of the reflected signals with the phase of the direct signal:

[TABLE]

where $({\bf h}_{c_{0}})_{n}$ denotes the $n$-th element of ${\bf h}_{c_{0}}$ and project it to $\{-1,+1\}$ for a feasible ${\bf y}_{\rm init}$. After the projected gradient descent, we sample $G$ feasible solutions and choose the best one. The complete procedure is described in Algorithm 1. Note that the numerical results associated with the case studies will be discussed in the next section.

III-B Overhead-aware Rate and EE maximization in an RIS-aided system

In order to tackle another canonical setting, we now focus on the RIS sub-problems where configuring the RIS to optimize the chosen performance metric requires finding the optimal number of reflecting elements $N_{\rm opt}$ first. However, this is only possible for simpler objective forms [46]. To address this limitation, we introduce a comprehensive stochastic sampling approach that optimizes more complex objective forms, including rate and EE while taking interference into account, circumventing the explicit calculation of $N_{\rm opt}$.

III-B1 System Model

Our system model is solely dictated by the signal model in (17). Along with that, we include the overhead and power consumption models developed in [46, 47]. We also assume that each RIS element has the ability to turn off or $\boldsymbol{\theta}\in\{-1,0,+1\}^{N}$. Note that, this allows us to avoid explicit derivation of $N_{\rm opt}$. We also assume that the estimated channels are reliable. Now, we define the rate of the system considering interference and channel estimation overhead below:

[TABLE]

where the bandwidth is denoted by $B$, the noise variance is $\sigma_{w}^{2}=BN_{0}$, $N_{0}$ is the noise power spectral density, total duration of the time slot is denoted by $T$, $T_{E}(\|\boldsymbol{\theta}\|_{0})$ denotes the time taken to estimate the channels, and $T_{F}(\|\boldsymbol{\theta}\|_{0})$ is the feedback duration of the RIS configuration. The channel estimation time and feedback duration time are dependent on the number of RIS elements $\|\boldsymbol{\theta}\|_{0}$ and are expressed next,

[TABLE]

where $T_{0}$ is the duration of each pilot tone, and $b_{F}=2$ is the number of bits used to represent the states of each RIS element. Additionally, $B_{F}$, $p_{F}$, and $h_{F}$ refer to the communication bandwidth, transmit power, and effective channel, respectively, in the feedback phase. Subsequently, the total power consumption can be expressed as

[TABLE]

where $P_{E}(\|\boldsymbol{\theta}\|_{0})=\frac{P_{0}T_{E}(\|\boldsymbol{\theta}\|_{0})}{T}$ is the power consumption in the channel estimation phase, $P_{0}$ denotes the power of each pilot tone, $p$ is the maximum transmit power in the data transmission phase, $P_{c,n}$ is the power required to operate each RIS element and $P_{c,0}$ is the static hardware power for the remaining system components. Additionally, $\frac{1}{\mu}$ and $\frac{1}{\mu_{F}}$ denote the transmit amplifier efficiency in the data transmission and feedback phase, respectively. Finally, the EE of the system is defined by,

[TABLE]

In this system model, the number of RIS elements is chosen to be $N=\min(N_{max},N_{0})$, where $N_{max}$ is a parameter and $N_{0}$ is the maximum integer for which $T_{E}(N_{0})+T_{F}(N_{0})<T$. This condition assures that the rate is realistic.

III-B2 RIS optimization

The general optimization problem can be expressed as, {mini}—s— θ ∈{-1,0,+1}^NJ(θ),

where $\mathcal{J}(\boldsymbol{\theta})$ can be $-R(\boldsymbol{\theta})$ or $-EE(\boldsymbol{\theta})$. As the objective function is non-smooth and non-convex due to the presence of the interference term and L0 norm, we resort to the Lemma 1 and reformulate the problem below: {mini}—s— p,q ∈(0,1)^NE_θ∼P_E( θ—p,q)[J(θ)], \addConstraintp+q≤1_N.

We assume that the $n$-th element of $\boldsymbol{\theta}$ or $\theta_{n}$ is an independent categorical random variable and can take the value $+1$ with probability $q_{n}$ that denotes the $n$-th entry of $\bf q$ and $-1$ with probability $p_{n}$ that denotes the $n$-th entry of $\bf p$. The joint PDF for the ternary random vector can be expressed as:

[TABLE]

The presence of coupled optimization variables in (III-B2) complicates the problem. To circumvent this, we implement the block coordinate descent (BCD) framework. This method decouples the problem with two coupled variable sets into two tractable sub-problems, each addressing a single set of variables while considering the other fixed. Not only does this approach simplify the problems, but it is also inspired by the ties between Dykstra’s algorithm for projections onto intersections of convex sets and BCD [48]. Ultimately, the SSA is used to resolve the sub-problems emerging from the BCD structure, as illustrated in the previous case study.

We start by taking the gradient of the objective function in (III-B2) assuming $\bf q$ is fixed. Note that, all the following derivations can easily be derived when $\bf p$ is fixed by substituting $\nabla_{\bf p}$ with $\nabla_{\bf q}$ and are not derived explicitly. However, we will provide those results in the appropriate lemmas.

[TABLE]

where $(a)$ comes from the definition of expectation, $\boldsymbol{\theta}\{k\}$ denotes the $k$-th possible combination out of the possible $3^{N}$ in an arbitrary indexing order, $(b)$ comes from the identity $\nabla_{\bf p}\log\mathbb{P}_{E}(\boldsymbol{\theta}|{\bf p},{\bf q})=\frac{1}{\mathbb{P}_{E}(\boldsymbol{\theta}|{\bf p},{\bf q})}\nabla_{\bf p}\mathbb{P}_{E}(\boldsymbol{\theta}|{\bf p},{\bf q})$, and $(c)$ converts the summation into expectation. The MC approximation of this gradient for a stochastic optimization approach is:

[TABLE]

where $N_{e}$ is the number of samples used. For completeness and to reduce the variance of this estimator without increasing $N_{e}$ drastically, we introduce a baseline $b_{\bf p}$ in the objective function. Such an estimator has the following form,

[TABLE]

where ${\bf d}_{\bf p}=\frac{1}{N_{e}}\sum\limits_{j=1}^{N_{e}}\nabla_{\bf p}\log\mathbb{P}_{E}(\boldsymbol{\theta}\{j\}|{\bf p},{\bf q})$. Note that ${\rm E}[{\bf d}_{\bf p}]=0$ as the following results stands:

[TABLE]

where $(a)$ comes from writing out the expectation in a summation. Using this result, we can also show that both estimators are also unbiased and ${\rm E}[\tilde{\bf g}_{\bf p}]={\rm E}[\hat{\bf g}_{\bf p}]={\bf g}_{\bf p}$. In the next lemma, we include the key gradient results for the ternary random vector that is instrumental in the stochastic sampling approach.

Lemma 4**.**

The gradient of the log of joint PDF with respect to $\bf p$ and $\bf q$ are:

[TABLE]

where ${\bf e}_{n}$ is the $n$-th unit vector of length $N$.

Proof:

We start by substituting the joint PDF:

[TABLE]

Similarly, the gradient $\nabla_{\bf q}\log\mathbb{P}_{E}(\boldsymbol{\theta}\{j\}|{\bf p},{\bf q})$ can be derived with ease. ∎

With these important gradients available, we find the optimal baseline for the estimator defined in (36) in the next lemma.

Lemma 5**.**

The optimal baselines when with respect to $\bf p$ and $\bf q$ are

[TABLE]

where $\hat{\bf g}_{\bf q}$ and ${\bf d}_{\bf q}$ can be found by replacing the $\nabla_{\bf q}$ in place of $\nabla_{\bf p}$ in (35) and (36), respectively.

Proof:

See Appendix -D. ∎

Remark 3**.**

We can also approximate ${\rm E}[\hat{\bf g}_{\bf p}^{T}{\bf d}_{\bf p}]$ by taking $b_{m}$ batches of $N_{e}$ data points and average them to get $\hat{b}_{\bf p}^{*}$ to use in the algorithm.

Now we have all the information to develop the stochastic sampling approach for ternary random variables. The algorithm to solve the sub-problems is demonstrated in Algorithm 2 and the BCD architecture is illustrated in Algorithm 3. In the Algorithm 2, the entries of ${\bf r}_{\rm init}$ are independent and identically distributed (i.i.d) with uniform distribution $\mathcal{U}(0,r_{max})$, where $0<r_{max}\leq 1$. By choosing a small $r_{max}$, we control the initial sparsity of the solution.

III-C Worst-case computational complexity discussion

In this subsection, we derive the worst-case computational complexities for the proposed algorithms in terms of big-O notation. However, we would like to note that the complexity of gradient descent-based algorithms cannot be trivially expressed in the big-O notation, as the number of iterations for convergence heavily depends on the initial point and cannot be precisely determined [49]. In the literature, the number of iterations is regarded as a parameter, and subsequently, the complexity is represented using big-O notation [50, 51]. In this subsection, we follow the same approach, while also preserving more terms in the big-O expression for a better comparison among the proposed algorithms. Building upon the previous discussion, The algorithms are based on five fundamental operations: gradient calculation, descent-projection, inner looping, outer looping, and sampling. The descent-projection operation has a complexity of $O(N)$ across all algorithms. Regarding the inner looping operation, we need $I_{E_{1}},I_{E_{2}},I_{1},$ and $I_{2}$ iterations for the gradient descent algorithms to converge for E-GD using first and second-order Taylor approximations, and for the stochastic sampling approaches with binary and ternary variables respectively. There is typically no need for outer looping iterations except for the ternary variable stochastic sampling method due to the BCD framework. In this case, we assume that $I_{BCD}$ iterations are needed to achieve convergence.

III-C1 E-GD with first-order Taylor series approximation

For this algorithm, the gradient calculation step is primarily dictated by the matrix multiplications inherent in the quadratic forms of (24), with a complexity of $O(N^{2})$. The sampling step adds an $O(GN^{2})$ complexity due to $G$ evaluations of the objective function, making the total complexity $O(I_{E_{1}}(N^{2}+N)+GN^{2})$.

III-C2 E-GD with second-order Taylor series approximation

For this variant of the algorithm, the gradient calculation step is primarily affected by the matrix multiplications required for computing the matrix ${\bf U}$ as per Theorem 2, and has a complexity of $O(N^{4})$. The other steps share the same complexities as the first-order version, yielding a total complexity of $O(I_{E_{2}}(N^{4}+N)+GN^{2})$.

III-C3 Stochastic sampling for binary variables

The first step of this algorithm, the gradient estimator calculation, is dominated by the $N_{ens}$ objective evaluations resulting in a complexity of $O(N_{ens}N^{2})$. The sampling step carries a complexity of $O(G_{s}N^{2})$ due to $G_{s}$ objective function evaluations, which makes the overall complexity $O(I_{1}(N_{ens}N^{2}+N)+G_{s}N^{2})$.

III-C4 Stochastic sampling for ternary variables

The main differences between this algorithm and the binary variant lie in the objective function evaluation, which has a complexity of $O(N^{2}+N)$ due to the additional L0 norm calculation. This yields a total complexity of $O(I_{BCD}I_{2}(N_{ens}(N^{2}+N)+N)+G_{s}N^{2})$.

IV Simulation Results

For our simulation results, we focus on a canonical (and perhaps most practically relevant) RIS scenario where RIS can significantly enhance performance: the creation of virtual line-of-sight (LoS) links when direct paths are obstructed, as highlighted in [14, 21, 52]. We maintain this assumption throughout our simulation. We also consider that we operate in a high interference regime where one interferer exists with average power similar to our user. This also highlights the ability of our developed algorithms to cope with high interference. The common simulation parameters used in both the cases are $\beta_{i}=p\delta_{PL}$, $p$ is the transmit power, $\delta_{PL}=-110$ dB, $B=5$ MHz, and $N_{0}=-174$ dBm/Hz [46]. Additionally, all the channels are Rician distributed with the Rician factor of $4$ while all the results in this section are averaged over $1000$ independent channel realizations.

IV-A SINR maximization with RIS optimization

In this application, we compare the achievable capacity $C_{cap}=\log_{2}(1+\gamma)$ of our developed algorithms with the popular SDR method and the CPP methods. The transmit power $p$ is considered to be [math] dBm. The algorithm parameters are $\varrho=0.5,$, $\varepsilon=0.0005$, $\epsilon_{th}=10^{-2}$, $\beta_{\rm init}=0.01$, and $G=100$. Our proposed first-order and second-order analytical GD algorithms are denoted by E-GD-1 and E-GD-2, respectively while our proposed stochastic sampling approach is denoted by SSA-B. The solution of the GD algorithm developed in [21] for continuous phase shifts projected to the discrete phase-shifts also acts like a baseline and is denoted by CPP-1. The CPP of the solution of (19) when the constraint is relaxed to be continuous is denoted by CPP-2. Note that, the only difference between E-GD-1 and CPP-2 is the final sampling step as the former treats the solution as a probability vector, and the latter projects it to $\{-1,+1\}$ for a solution. The CPP of the simple signal alignment scheme in (28) is denoted by SA. CPP methods are considered comparison baselines as they are more practical in terms of speed and are often used in the literature over the traditional branch-and-bound methods that do not scale well with the number of elements.

In Fig. 1(a), we can observe that all the expectation-based algorithms perform better than the CPP algorithms, for all $N$, and the SDR for $N>20$. Along with that, SSA-B outshines the expectation-based EGD algorithms that utilize approximations for expectation computation. The edge of SSA-B lies in its robust gradient estimates, derived without reliance on Taylor series approximations, hence providing more precise results. Moreover, the accuracy of SSA-B’s gradient estimates can be enhanced by increasing the sample size, although this incurs a higher computational cost. Moreover, the scheme CPP-1 performs worse compared to CPP-2 due to its design for continuous RIS phase-shifts with unit-modulus constraints, which means its RIS optimization variable domain spans all angles from [math] to $2\pi$ corresponding to the set of all unit-modulus complex gains. In contrast, the domain of CPP-2 ranges from $-1$ to $1$, making it closer to the original domain of $\{-1,+1\}$. This difference gets more prominent as the number of RIS elements grows and CPP-2 provides a sharper increase in achievable capacity than CPP-1.

in Fig. 1(b), we plot the run-time for a single iteration of all the algorithms with varying numbers of RIS elements. These results are taken from the simulations needed to create Fig. 1(a) on a 3.6GHz Intel Core i7-4790 8-CPU system with 16GB RAM. From this plot, we note that the runtime of SSA-B is between the E-GD-1 and E-GD-2 methods while SDR is prohibitively slow. The runtime of our proposed E-GD-2 method is better than SDR but still slower than its first-order counterpart due to the complex gradient calculation. The overall performance of our analytical GD algorithms is dependent on the trade-off between the complexity of the gradient and the accuracy of the approximation for the expectation. These simulation results demonstrate the superiority of the expectation-based algorithms in discrete optimization problems providing important insights into such analytical expectation derivation.

IV-B Overhead-aware rate and EE maximization in an RIS-aided system

In this application, we maximize the rate and the EE of the system with our stochastic sampling approach. As a baseline, we compare it with the solution in [46] without interference projected to the discrete RIS phase-shifts. It should be noted that when interference exists, this baseline is no longer relevant because the unimodality required to compute the optimal number of RIS elements is dependent on the simple objective structure without interference. The simulation parameters are set according to [46]: $B_{F}=1$ MHz, $P_{c,0}=45$ dBm, $P_{c,n}=10$ dBm, $\mu=\mu_{F}=1$, $T_{0}=1$ ms, $p_{F}=30$ dBm, $P_{0}=10$ dBm, $T=100$ ms, and $N_{max}=300$. The optimization algorithm parameters are, $\epsilon=10^{-6}$, $N_{e}=200$, $b_{m}=10$, $r_{max}=0.1$, $t_{max}=300$, $G_{s}=10000$, $\epsilon_{t}=10^{-8}$, $\beta_{s}=0.5$ for EE and $\beta_{s}=0.01$ for rate optimization. In the simulation figures, the upper bound is calculated with the optimum continuous phase shifts without interference through the unimodal approach (UA) devised in [46]. The CPP of this approach also acts as a baseline and is denoted by UA while our algorithm is denoted by SSA-T or stochastic sampling approach for the ternary variable.

In Fig. 2(a), we plot the average EE achieved with the transmit power when interference is not present. For $T_{0}=1$ ms, our algorithm performs very similarly to the unimodal approach. However, for $T_{0}=0.2$ ms, our algorithm achieves an EE that is $0.18$ Mbit/J less at $p=30$ dBm than the UA. While the UA method offers optimal results in the continuous RIS case where no interference is present - a scenario that can be viewed as a special instance of the general formulation with zero interference - it naturally extends well to the discrete case as well. In contrast, our proposed algorithm has a broader scope, demonstrating its capability to handle any form of objective function. Despite this versatility, the trade-off is a guarantee of optimality, hence the observed performance is completely expected. Our algorithm uniquely excels in managing interference and can adapt to any general objective form, an area where the UA method notably underperforms. Therefore, the superior performance of the UA approach in this specific case is anticipated, as it was designed precisely for such interference-free scenarios. This distinction underscores the unique use case of the stochastic sampling approach: when a reliable solution for the continuous problem exists, discrete projection may be sufficient. However, when the objective function becomes complex, even in its continuous form, our proposed algorithm shines, providing high-quality solutions where other methods might fall short. We can also observe that at the high transmit power region, the EE drops as power consumption dominates and our algorithm approaches the unimodal approach and the upper bound.

In Fig. 2(b), we plot the average EE achieved with the transmit power when interference is present. In situations where interference is present, as noted in Section III, the transmit power displayed on the x-axis represents the maximum transmit powers of both the user and the interferer. This mimics a non-cooperative scenario where each entity maximizes its own performance by transmitting at peak power simultaneously. We observe that our algorithms continue to perform close to the upper bound whereas the unimodal approach fails as expected. We can also observe the general trend of better performance with decreasing $T_{0}$. The reasoning is two-fold: a) we have more time to transmit data due to lower channel estimation time, and b) more overhead for RIS elements can be supported, resulting in the utilization of more RIS elements. This can be verified in Table I where we report the average number of RIS elements to maximize EE in presence of interference varying with the transmit power. We can also observe the effect of $G_{s}$ in this figure. With a lower $G_{s}=100$, the EE achieved is around $0.62$ Mbit/J less than the default parameter $G_{s}=10000$ case at $p=30$ dBm. Finally, in Fig. 2(c), we observe that in the rate maximization problem, the proposed algorithm performs closer to the upper bound than the UA approach, irrespective of interference. Without interference, the performance of SSA-T marginally exceeds that of UA. However, in the presence of interference, the advantage of SSA-T over UA becomes substantial. Furthermore, as we increase the maximum transmit power, it becomes clear that the achievable rate via our proposed approach reaches a saturation limit. This limit is imposed by the proportional increase in interference power which is not completely suppressed by the discrete RIS along with user transmit power.

V Conclusion

In this paper, we developed a novel probabilistic reformulation technique for general discrete optimization problems. In particular, we interpret the discrete optimization variable as a categorical random vector and take expectations on the objective function along with any constraints present. We provide rigorous mathematical justification that the corresponding degenerate PDF of the unique optimal solution of an unconstrained problem is the unique optimal solution of the transformed problem and for a constrained problem, the primal solution of the transformed problem is bounded between the dual and primal solutions of the original problem implying that it is a relaxation of the original problem. However, if strong duality holds, the transformed problem provides the same objective value as the original constrained problem. We also explored a simple two-way partitioning problem to gain more insights into our reformulation such as its similarity to SDR, and capability to change the problem structure. We ultimately used this technique to tackle two canonical discrete RIS applications: a) SINR maximization, and b) overhead-aware rate and EE maximization. As demonstrated in our RIS applications, the reformulation allows for both stochastic and analytical interpretations of the original problems. For the SINR maximization problem, an analytical GD technique based on closed-form approximations for the expectation is proposed, while a stochastic sampling approach is proposed for both applications. The numerical results reveal that there is a fundamental trade-off between the complexity of the gradient and the accuracy of the approximation in our proposed analytical GD methods, and the expectation-based algorithms outperform the other algorithms evaluated. The simulation results also demonstrate that our proposed framework is very general and performs well for both rate and EE maximization problems without much change in the algorithm. In particular, we show that it performs at par with the algorithm specifically developed for the interference-free case when interference is not present and keeps performing well even when interference exists. We also explicitly calculate the worst-case computational complexities for our proposed algorithms. As the scope of this technique is very general, utilizing this technique to develop a more sophisticated projected gradient descent framework and a general methodology to deal with constrained problems are left as future work.

-A Proof of Theorem 1

The expectation can be calculated by converting the matrix expressions into series sums as shown below,

[TABLE]

In step $(a)$, we use the fact that $\bf G$ is real symmetric, $x_{i}^{2}=1$, and ${\rm E}[x_{i}]=y_{i}$. The third term can be expressed in the following form:

[TABLE]

where ${\bf G}_{i}$ denotes the matrix ${\bf G}$ with the $i$-th row and column set to zeros. In matrix form, this can be expressed as ${\bf 1}^{T}\{({\bf G}_{wd}{\bf Y}_{wd})\odot{\bf Y}_{wd}\}({\bf y}\odot{\bf z})$, where ${\bf Y}={\bf y}{\bf 1}^{T}$. This completes the proof.

-B Proof of Theorem 2

We start this proof by expanding a generic quadratic term in (44).

In (44), the matrix ${\bf U}=[{\bf I}_{N}\otimes({\bf y}\odot{\bf y})^{T}]{\bf B}$, and the matrix ${\bf B}$ is defined through blocks as

[TABLE]

where the $i$-th element of ${\bf b}_{k,j}$ is ${\bf b}_{k,j}^{i}=G_{{wd}_{ij}}G_{{wd}_{ki}}$.

Considering that the $d$-th term in the final expression without the numeric coefficient is denoted by $S_{d}({\bf x})$, the above expression can be expressed as

[TABLE]

Considering $x_{i}^{2}=1$, the second moment of a quadratic form can be expressed as,

[TABLE]

where ${\bf Z}={\bf G}_{wd}{\bf G}_{wd}^{T}$. This is readily found by taking expectation on (44). Note that the final term or $S_{5}({\bf y})$ can be found from the following observation:

[TABLE]

The theorem is proved by combining the final two expressions.

-C Proof of Corollary 2

As the gradient of most of the terms can be trivially calculated [53, 54], we focus on the nontrivial gradient calculations here. In particular, the gradient of ${\bf y}^{T}{\bf U}_{wd}{\bf y}$ with respect to y is derived next. We can write the following expression due to the chain rule:

[TABLE]

where $({\bf U}_{wd})_{jk}$ is $(j,k)$-th element of the matrix ${\bf U}_{wd}$ and the matrix ${\bf U}$ can be expressed as

[TABLE]

With this formulation, the inner derivative is $\frac{\partial}{\partial y_{i}}\left(({\bf U}_{wd})_{jk}\right)=2y_{i}({\bf b}_{j,k})_{i}\quad\forall j\neq k$, where $({\bf b}_{j,k})_{i}$ is the $i$-th element of the vector ${\bf b}_{j,k}$. Finally, the gradient can be written as,

[TABLE]

where $(a)$ is obtained by some matrix manipulations and the definition of ${\bf b}_{j,k}$ vectors.

-D Proof of Lemma 5

We begin by calculating the total variance of the estimator $\tilde{\bf g}_{\bf p}$ below.

[TABLE]

where ${\rm var}({\bf d}_{\bf p})=\frac{1}{N_{e}^{2}}\sum\limits_{j=1}^{N_{e}}{\rm var}\left(\nabla_{\bf p}\log\mathbb{P}_{E}(\boldsymbol{\theta}\{j\}|{\bf p},{\bf q})\right)$ and as the variance does not depend on the $j$-th index, we can calculate the variance of the inner quantity next ignoring the index.

[TABLE]

Using (52) in (51), we can write that,

[TABLE]

where $(a)$ is a result of the following observations: ${\rm E}[{\theta_{n}^{2k_{0}}}]=q_{n}+p_{n}$ and ${\rm E}[{\theta_{n}^{2k_{0}+1}}]=q_{n}-p_{n}$, where $k_{0}$ is a non-negative integer. Note that this is a convex quadratic expression in $b_{\bf p}$ and we can find the minimum by equating the derivative of this variance equal to zero. The optimal baseline is $b_{\bf p}^{*}=\frac{N_{e}}{\sum\limits_{n=1}^{N}\frac{1}{p_{n}}+\frac{1}{1-p_{n}-q_{n}}}{\rm E}[\hat{\bf g}_{\bf p}^{T}{\bf d}_{\bf p}]$ and the lemma is proved.