# Toward Energy Efficient Multiuser IRS-Assisted URLLC Systems: A Novel Rank
Relaxation Method

Jalal Jalali†, Filip Lemic♯, Hina Tabassum⋆, Rafael Berkvens†, and Jeroen
Famaey† †IDLab, University of Antwerp - imec, Sint-Pietersvliet 7, 2000
Antwerp, Belgium
♯AI-Driven Systems Lab, i2Cat Foundation, Barcelona, Spain
⋆Department of Electrical Engineering and Computer Science, York University,
Toronto, ON, Canada
Email<EMAIL_ADDRESS>
###### Abstract

This paper proposes an energy efficient resource allocation design algorithm
for an intelligent reflecting surface (IRS)-assisted downlink ultra-reliable
low-latency communication (URLLC) network. This setup features a multi-antenna
base station (BS) transmitting data traffic to a group of URLLC users with
short packet lengths. We maximize the total network’s energy efficiency (EE)
through the optimization of active beamformers at the BS and passive
beamformers (a.k.a phase shifts) at the IRS. The main non-convex problem is
divided into two sub-problems. An alternating optimization (AO) approach is
then used to solve the problem. Through the use of the successive convex
approximation (SCA) with a novel iterative rank relaxation method, we
construct a concave-convex objective function for each sub-problem. The first
sub-problem is a fractional program that is solved using the Dinkelbach method
and a penalty-based approach. The second sub-problem is then solved based on
semi-definite programming (SDP) and the penalty-based approach. The iterative
solution gradually approaches the rank-one for both the active beamforming and
unit modulus IRS phase-shift sub-problems. Our results demonstrate the
efficacy of the proposed solution compared to existing benchmarks.

###### Index Terms:

Alternating optimization (AO), intelligent reflecting surface (IRS), semi-
definite programming (SDP), ultra-reliable low-latency communication (URLLC).

##  I Introduction

Intelligent reflecting surfaces (IRSs) have recently emerged as a cost-
effective and energy efficient technology to reconfigure wireless environments
through passive reflections. Unlike conventional relays, IRSs do not require
signal decoding or amplification [1]. The integration of IRSs into the
communication network offers enhanced reliability, reduced packet re-
transmissions, and latency. Consequently, IRSs can efficiently realize ultra-
reliable low-latency communications (URLLC) with finite block length regime
[2]. Unlike infinite blocklength regime, the achievable rate in the finite
blocklength regime is a complex function of the signal-to-interference-plus-
noise ratio (SINR), blocklength, and the probability of decoding errors [3].

Recently, a handful of research studies investigated resource allocation in
IRS-assisted URLLC networks. For instance, Hashemi _et al._ presented the
performance analysis of the average achievable data-rate and error probability
over an IRS-aided URLLC transmission with and without phase noise [4].
Considering non-linear energy harvesting, the end-to-end performance of the
IRS-assisted wireless system was analyzed in [5] for industrial URLLC
applications, and the approximate closed-form expression of the block error
rate was derived. Xie _et al._ studied an IRS-assisted downlink multiuser
URLLC system and jointly optimized the user grouping and the blocklength
allocation at the base station (BS), as well as the reflective beamforming at
the IRS for latency minimization [6]. Li _et al._ investigated maximizing
uplink data-rate in a URLLC system by jointly optimizing user transmit power,
BS receive beamformer, and IRS reflection phase shifts under given error
probability and blocklength conditions [7]. Very recently, an IRS-assisted
machine-type communication network was considered in [8] where multi-objective
optimization was applied to simultaneously maximize EE and the number of
admitted users with short packet length. The problem is solved by using
fractional programming approach, i.e., quadratic transformation. where a non-
practical signal-to-interference-plus-noise ratio (SINR) bound was considered
with no hard constraint on the delay [8] None of the aforementioned research
works considered an energy efficient design of IRS-assisted URLLC networks,
although other performance metrics were considered in [9, 10, 11, 12, 7, 13,
14].

In this paper, we design an energy efficient beamforming and phase-shift
optimization algorithm for downlink IRS-assisted URLLC services, wherein a
multi-antenna BS serves multiple single-antenna URLLC receivers with the
support of an IRS. To this end, we formulate an energy efficiency maximization
problem, subject to the minimum required data-rate (in the finite blocklength
regime) for each URLLC user and unit-modulus constraints at the IRS. We
decompose the problem into two sub-problems and apply an alternating
optimization (AO) approach to optimize both active and passive beamformers at
the BS and IRS, respectively. Using an innovative iterative rank relaxation
and successive convex approximation (SCA) method, we construct a concave-
convex objective function for each sub-problem. Then, the first sub-problem
(which is also a fractional program) is solved by applying the Dinkelbach
method and a penalty-based approach; whereas, the second sub-problem is solved
by applying the penalty-based approach and semi-definite programming (SDP).
The iterative solution gradually approaches the rank-one for both the active
beamforming and unit modulus IRS phase-shift sub-problems. Our simulation
results indicate the benefits of the IRS in enhancing the energy efficiency of
URLLC systems and provide a comparison with benchmarks.

Notations: In this work, matrices and vectors are denoted by boldface capital
letters, like $\mathbf{A}$, and lower-case letters, such as $\mathbf{a}$,
respectively. If we consider a square matrix $\mathbf{A}$, $\mathbf{A}^{T}$,
$\mathbf{A}^{H}$, ${\rm{rank}}(\mathbf{A})$, $\text{Tr}(\mathbf{A})$, and
$||\mathbf{A}||_{*}$ denote its transpose, Hermitian conjugate transpose,
rank, trace, and norm, respectively. $\mathbf{I}_{N}$ represents the $N\times
N$ identity matrix. The diagonalization operation is represented by
$\rm{diag}(\cdot)$. $\Re\\{\cdot\\}$ is used to denote the real part of a
complex function. We express the absolute value of a complex scalar and the
Euclidean norm of a complex vector by $|\cdot|$. The distribution of a
circularly symmetric complex Gaussian (CSCG) random vector with mean
$\boldsymbol{\mu}$ and covariance matrix $\mathbf{C}$ is denoted as
$\mathcal{C}\mathcal{N}(\boldsymbol{\mu},\mathbf{C})$. $Q^{-1}(\cdot)$
represents the inverse of the Gaussian Q-function. $\mathbb{C}^{M\times N}$
denotes an $M\times N$ dimensional complex matrix, and $\nabla_{\mathbf{x}}$
expresses the gradient vector with respect to $\mathbf{x}$. Finally, we
express ${N\times N}$ positive semi-definite matrices as
$\mathbf{A}\in\mathbb{S}^{N}_{+}$ and read $\mathbf{A}\succeq\mathbf{0}$.

This paper is structured as follows: In Section II, we introduce the system
model and outline the proposed EE optimization problem. The resource
allocation algorithm to solve the EE problem is detailed in Section III. In
Section IV, we assess the performance of our novel algorithm. Finally, Section
V draws conclusions.

Figure 1: Illustration of a multi-user IRS-assisted downlink network
comprising one BS and $K$ URLLC users employing finite block length
transmissions. ${\mathbf{H}}$ is the channel matrix between the BS and the
IRS, ${\mathbf{h}}_{\rm{IRS},1}$ denotes the channel response vector from the
IRS to user 1, and $\mathbf{h}_{\rm{BS},1}$ indicates the channel between the
BS and user 1.

##  II System Model and Problem Formulation

We consider a downlink multiuser URLLC-enabled IRS-assisted system,
represented in Fig. 1, which comprises an IRS with $N$ elements, a BS equipped
with $M$ antennas, and $K$ users, each with a single antenna. The group of IRS
elements, BS antennas, and users are denoted by the sets
$\mathcal{N}=\\{1,...,N\\}$, $\mathcal{M}=\\{1,...,M\\}$, and
$\mathcal{K}=\\{1,...,K\\}$, respectively. Additionally, we suppose that
$B_{k}$ bits of information are assigned to user $k$. In this setup, the BS
transforms these information bits into a block code in time slot $l$,
characterized by a length of $m_{d}$ symbols. We decode the $k$-th user block
code as $x_{k,l}$, with $l$ being part of the set
$\mathcal{L}=\\{1,2,...,m_{d}\\}$.

The transmission signal emanating from the BS can be expressed as
$\mathbf{s}{[l]}=\sum_{k\in\mathcal{K}}{{\boldsymbol{\omega}_{k,l}}{x_{k,l}}}$,
where ${\boldsymbol{\omega}}_{k,l}\in\mathbb{C}^{M\times 1}$ represents the
beamforming vector for user $k$. The channel links exhibit time invariance
(i.e., slow fading). Additionally, we assume that the BS possesses full access
to the channel state information (CSI) and URLLC users’ delay requirements
[10, 8, 15, 16]. The baseband equivalent channel responses for BS-to-IRS, IRS-
to-user $k$, and BS-to-user $k$ are denoted as
${\mathbf{H}}\in\mathbb{C}^{N\times M}$,
${\mathbf{h}}_{\rm{IRS},k}\in\mathbb{C}^{N\times 1}$, and
$\mathbf{h}_{\rm{BS},k}\in\mathbb{C}^{M\times 1}$, respectively. Also, we
define the matrix of reflection coefficients at the IRS as
$\boldsymbol{\Psi}_{l}=\text{diag}(\alpha_{1,l}e^{j\phi_{1,l}},\alpha_{2,l}e^{j\phi_{2,l}},...,\alpha_{N,l}e^{j\phi_{N,l}})$,
where $\alpha_{n,l}\in[0,1]$ and $\phi_{n,l}\in(0,2\pi]$, $\forall
n\in\mathcal{N},\forall l\in\mathcal{L}$ are the reflection amplitude and
phase shift of the $n$-th reflection coefficient at the IRS during time slot
$l$, respectively.

By defining the combined channel link seen by the $k$-th URLLC user as
$\mathbf{h}_{k}^{H}\triangleq{\mathbf{h}}_{\rm{IRS},k}^{H}\boldsymbol{\Psi}_{l}{\mathbf{{\mathbf{H}}}}+\mathbf{h}_{\rm{BS},k}^{H}$,
$\forall k\in\mathcal{K},$ $\forall l\in\mathcal{L}$, the SINR of that user
can be expressed as:

$\Gamma_{k,l}=\frac{{{\left|\mathbf{h}_{k}^{H}{\boldsymbol{\omega}}_{k,l}\right|}^{2}}}{\sum\limits_{i\neq
k,i\in\mathcal{K}}{{\left|\mathbf{h}_{k}^{H}{\boldsymbol{\omega}}_{i,l}\right|}^{2}}+\sigma_{k}^{2}},\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (1)

where $\sigma_{k}^{2}$ represents the noise variance received at user $k$. In
the context of URLLC systems, an accurate estimate of the achievable data-rate
for each user can be delineated as follows [17]:

$\displaystyle\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})=\mathcal{U}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})-\mathcal{V}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}),~{}\forall
k\in\mathcal{K},$ (2)

where

$\displaystyle\mathcal{U}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
$\displaystyle=\sum\limits_{l\in\mathcal{L}}\log_{2}(1+{\Gamma_{k,l}}),\ \ \ \
\ \ \ \ \forall k\in\mathcal{K},$ (3)
$\displaystyle\mathcal{V}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
$\displaystyle=\sum\limits_{l\in\mathcal{L}}\beta_{k}\sqrt{\Delta_{k,l}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},$ (4) $\displaystyle\beta_{k}$
$\displaystyle=\frac{Q^{-1}(\epsilon_{k})}{\sqrt{m_{d}}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K}.$ (5)

Here, $\epsilon_{k}$ refers to the decoding error probability, $m_{d}$, as
defined earlier, is representative of the blocklength, and $\Delta_{k,l}$ the
channel dispersion, computed as:

$\Delta_{k,l}=(\log_{2}e)^{2}\left(1-\frac{1}{\left(1+\Gamma_{k,l}\right)^{2}}\right),\forall
k\in\mathcal{K},\forall l\in\mathcal{L}.$ (6)

To ensure users’ quality of service (QoS) regarding the received number of
bits, the reliability, and the latency, a minimum data-rate denoted by
$R_{\min,k}$ should be satisfied for each user as follows:

$\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\\!\boldsymbol{\Psi}_{l})\geq
R_{\min,k},\>\forall k\in\mathcal{K}.$ (7)

Next, we describe the network’s energy efficiency (EE) as the ratio of the
total system data-rate over the associated network power consumption in
$\rm{[bits/Joule/Hz]}$:

$\mathscr{\eta}_{\mathrm{eff}}(\boldsymbol{\omega}_{k,l},\\!\boldsymbol{\Psi}_{l})=\frac{\sum\limits_{k\in\mathcal{K}}\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\\!\boldsymbol{\Psi}_{l})}{\sum\limits_{k\in\mathcal{K}}{{\sum\limits_{l\in\mathcal{L}}{\left\|{\boldsymbol{\omega}}_{k,l}\right\|}^{2}}}+P_{\rm{IRS}}+NP_{d}+P_{\rm{c}}^{\rm{BS}}},$
(8)

where $P_{\rm{IRS}}$ represents the static power consumption necessary for
sustaining the basic circuit operations of the IRS, while $P_{d}$ is the
dynamic power expended per reflecting component, and $P_{\rm{c}}^{\rm{BS}}$
denotes the circuit power at the BS.

Our objective is to maximize the total EE by optimizing the active beamformers
at the BS and the phase shifts at the IRS. Consequently, the problem can be
articulated as follows:

$\displaystyle\text{P}_{1}$
$\displaystyle:\mathop{\max}\limits_{\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}}\>\>\>\mathscr{\eta}_{\mathrm{eff}}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
$\displaystyle s.t.$
$\displaystyle:\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\\!\boldsymbol{\Psi}_{l})\geq
R_{\min,k},\forall k\in\mathcal{K},$ (9a)
$\displaystyle~{}~{}|\boldsymbol{\Psi}_{{nn},l}|=1,~{}~{}\ \ \ \ \ \ \ \
~{}~{}~{}~{}\forall n\in\mathcal{N},\forall l\in\mathcal{L},$ (9b)
$\displaystyle~{}~{}\sum\limits_{l\in\mathcal{L}}\sum\limits_{k\in\mathcal{K}}{\left\|{\boldsymbol{\omega}}_{k,l}\right\|}^{2}\leq
p_{\max},$ (9c) $\displaystyle~{}~{}\boldsymbol{\omega}_{k,l}=0,~{}~{}\ \ \ \
\ \ \ \ ~{}~{}~{}~{}~{}~{}~{}\>\forall k\in\mathcal{K},\forall
l\geq{\rm{T}}_{k},$ (9d)

where (9a) is the minimum data-rate requirement ($R_{\min,k}$) of the $k$-th
URLLC user. Constraint (9b) seeks that the diagonal phase shift matrix
contains $N$ unit-modulus elements along its main diagonal. The constraint
(9c) defines the limitation of the BS transmission power budget, with
$p_{\max}$ signifying the maximum permissible BS transmission power.
Constraint (9d) is imposed to safeguard the real-time URLLC service
functionality, ensuring that user $k$ receives service within the first
${\rm{T}}_{k}$ time slots to satisfy its delay requirements.

Given the non-convex nature of the objective function and constraints in
$\text{P}_{1}$, in the following, we propose an AO-based low-complexity
solution.

##  III Energy Efficient Design for IRS-Assisted URLLC Network

In this section, we decompose $\text{P}_{1}$ into two sub-problems. The first
sub-problem focuses on determining the optimal active beamformers at the BS,
while the second sub-problem concentrates on optimizing the phase shifts at
the IRS. We adopt the SCA technique to initially derive a lower bound for
$\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$, which will
subsequently be utilized when addressing each sub-problem.

###### Lemma 1

Let $(j)$ denote the subscript associated with the feasible solution procured
in the $j$-th iteration of the SCA algorithm. A lower bound approximation of
$\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$ is given
by:

$\displaystyle\mathcal{R}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})\geq\tilde{\mathcal{R}}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})=\sum\limits_{l\in\mathcal{L}}\tilde{\mathcal{R}}_{k,l}^{(j)}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}),$
(10)
$\displaystyle\\!\\!\\!\tilde{\mathcal{R}}_{k,l}^{(j)}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})\triangleq\\!-\gamma_{k,l}\mathfrak{b}_{k,l}\\!+\\!2\rho_{k,l}\Re\left\\{\\!\\!\sum\limits_{i\neq
k,i\in\mathcal{K}}\mathfrak{a}_{i,k,l}^{(j)}\mathfrak{a}_{i,k,l}^{H}\\!+\\!\sigma^{2}_{k}\\!\right\\}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\>+2\Re\left\\{\frac{\mathfrak{a}_{k,k,l}^{(j)}\mathfrak{a}_{k,k,l}^{H}}{{\mathfrak{b}_{k,l}^{(j)}(\mathfrak{b}_{k,l}^{(j)}-|\mathfrak{a}_{k,k,l}^{(j)}|^{2})}}\right\\}+\xi_{k,l},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (11)

where

$\displaystyle\mathfrak{a}_{k,i,l}$
$\displaystyle=\mathbf{h}_{k}^{H}{\boldsymbol{\omega}}_{i,l},$ (12)
$\displaystyle\mathfrak{b}_{k,l}$
$\displaystyle={\sum\limits_{i\in\mathcal{K}}{{\left|\mathbf{h}_{k}^{H}{\boldsymbol{\omega}}_{i,l}\right|}^{2}}+\sigma_{k}^{2}},$
(13) $\displaystyle\rho_{k,l}$
$\displaystyle=\frac{\beta_{k}(\log_{2}e)^{2}}{\mathfrak{b}_{k,l}^{(j)}\sqrt{(\Gamma^{(j)}_{k,l})^{2}+2\Gamma^{(j)}_{k,l}}}$
(14) $\displaystyle\gamma_{k,l}$
$\displaystyle=\frac{|\mathfrak{a}_{k,k,l}^{(j)}|^{2}}{\mathfrak{b}_{k,l}^{(j)}(\mathfrak{b}_{k,l}^{(j)}-|\mathfrak{a}_{k,k,l}^{(j)}|^{2})}+\frac{\rho_{k,l}(\mathfrak{b}_{k,l}^{(j)}-|\mathfrak{a}_{k,k,l}^{(j)}|^{2})}{\mathfrak{b}_{k,l}^{(j)}},$
(15)

$\displaystyle\xi_{k,l}$
$\displaystyle=\log_{2}(1+\Gamma_{k,l}^{(j)})-\frac{|\mathfrak{a}_{k,k,l}^{(j)}|^{2}}{\mathfrak{b}_{k,l}^{(j)}(\mathfrak{b}_{k,l}^{(j)}-|\mathfrak{a}_{k,k,l}^{(j)}|^{2})}$
$\displaystyle-\frac{\beta_{k}(2\Delta_{k,l}^{(j)}+(\log_{2}e)^{2})}{2\sqrt{\Delta_{k,l}^{(j)}}}-\frac{\rho_{k,l}\mathfrak{b}_{k,l}^{(j)}}{2(1+\Gamma^{(j)}_{k,l})}.$
(16)

In the $(j+1)$-th iteration, the data-rate function
$\tilde{\mathcal{R}}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
exhibits concavity with respect to each variable, reaching its boundary at the
point $(\boldsymbol{\omega}_{k,l}^{(j)},\boldsymbol{\Psi}_{l}^{(j)})$.

###### Proof 1

See Appendix A.

The lower bound presented is more manageable compared to the original data-
rate function
$\tilde{\mathcal{R}}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$ in
(2). However, this bound still involves coupled optimization variables. To
address this, we implement the AO approach. Specifically, we optimize
$\boldsymbol{\omega}_{k,l}$ and $\boldsymbol{\Psi}_{l}$ by alternatively
refining each variable while keeping the others constant.

###  III-A Sub-problem (1): Optimizing $\boldsymbol{\omega}_{k,l}$ with given
$\boldsymbol{\Psi}_{l}$

In this sub-problem, we presume that the passive reflecting elements at the
IRS, i.e., $\boldsymbol{\Psi}_{l}$ are fixed, and we proceed with designing
the active beamformers, $\boldsymbol{\omega}_{k,l}$, at the BS. By employing
the principles of semi-definite programming (SDP), we derive the following:
$\boldsymbol{\Omega}_{k,l}=\boldsymbol{\omega}_{k,l}\boldsymbol{\omega}_{k,l}^{H}\in\mathbb{S}^{M}_{+}$
and $\mathbf{H}_{k}=\mathbf{h}_{k}\mathbf{h}_{k}^{H}\in\mathbb{S}^{M}_{+}$,
$\forall k\in\mathcal{K}$. Leveraging on the aforementioned SDP relaxation, we
can rewrite the SINR in (1) as follows:

$\Gamma_{k,l}=\frac{{{\text{Tr}({\mathbf{{\rm{{\mathbf{H}}}}}}_{k}^{H}\boldsymbol{\Omega}_{k,l})}}}{{\sum\limits_{i\in\mathcal{K},i\neq
k}{{\text{Tr}({\mathbf{{\rm{{\mathbf{H}}}}}}_{k}^{H}\boldsymbol{\Omega}_{i,l})}+\sigma_{k}^{2}}}},\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (17)

Consequently, the data-rate function defined in (2) can be reshaped as
$\mathcal{R}_{k}(\boldsymbol{\Omega}_{k,l})=\mathcal{U}_{k}(\boldsymbol{\Omega}_{k,l})-\mathcal{V}_{k}(\boldsymbol{\Omega}_{k,l})$.
With the aid of the SDP transformations, the EE optimization problem
$\text{P}_{1}$ can be reformulated as follows:

$\displaystyle\text{P}_{2}$
$\displaystyle:\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l}}\>\>\>\mathscr{\eta}_{\mathrm{eff}}(\boldsymbol{\Omega}_{k,l})$
$\displaystyle s.t.$
$\displaystyle:~{}\mathcal{R}_{k}(\boldsymbol{\Omega}_{k,l})\geq
R_{\min,k},\forall k\in\mathcal{K},$ (18a)
$\displaystyle~{}~{}~{}\boldsymbol{\Omega}_{k,l}\succeq\mathbf{0},\ \ \ \ \ \
\ \ \ \ \ \ \>\>\forall k\in\mathcal{K},\forall l\in\mathcal{L},$ (18b)
$\displaystyle~{}~{}~{}{\rm{rank}}(\boldsymbol{\Omega}_{k,l})\leq 1,\ \ \ \>\
\ \forall k\in\mathcal{K},\forall l\in\mathcal{L},$ (18c)
$\displaystyle~{}~{}~{}\sum\limits_{l\in\mathcal{L}}\sum\limits_{k\in\mathcal{K}}\text{Tr}(\boldsymbol{\Omega}_{k,l})\leq
p_{{\max}},$ (18d)
$\displaystyle~{}~{}~{}\text{Tr}(\boldsymbol{\Omega}_{k,l})=0,~{}\ \ \ \ \ \ \
~{}\forall k\in\mathcal{K},\forall l\geq{\rm{T}}_{k}.$ (18e)

It is important to note that the constraint (18a) in $\text{P}_{2}$ does not
exhibit concavity. To circumvent this non-concavity, we employ the outcome
from Lemma 1. A proposed surrogate lower bound of
$\mathcal{R}_{k}(\boldsymbol{\Omega}_{k,l})$ is introduced, where it is
guaranteed that
$\mathcal{R}_{k}(\boldsymbol{\Omega}_{k,l})\geq\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})$.
Hence, we can re-express $\text{P}_{2}$ as:

$\displaystyle\text{P}_{3}$
$\displaystyle:\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l}}\>\>\>\tilde{\mathscr{\eta}}_{\mathrm{eff}}(\boldsymbol{\Omega}_{k,l})=\frac{\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})}{\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})}$
$\displaystyle s.t.$
$\displaystyle:~{}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})\geq
R_{\min,k},\ \ \ \ \ \forall k\in\mathcal{K},$ (19a)
$\displaystyle~{}~{}~{}\eqref{p2-2}-\eqref{p2_5}.$

where
$\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})=\sum_{l\in\mathcal{L}}\sum_{k\in\mathcal{K}}\text{Tr}(\boldsymbol{\Omega}_{k,l})+P_{\rm{IRS}}+NP_{d}+P_{\text{c}}^{\text{BS}}$.
The numerator of the objective function and constraint (19a) exhibit the form
of convex-concave functions. Nevertheless, the overall problem remains non-
convex due to the presence of the non-convex rank-one constraint (18c) and the
fractional objective function. Our forthcoming focus firstly centers on
resolving the non-convex rank-one constraint (18c), which emerges within the
context of our optimization problem. Fortunately, an innovative method has
been proposed in [18], specifically devised to address these types of rank
constraints. We hereby direct our attention to the ensuing proposition.

###### Proposition 1

Consider a nonzero positive semi-definite beamforming matrix, denoted as
$\boldsymbol{\Omega}_{k,l}\in\mathbb{S}^{M}_{+}$, where $M$ represents the
dimensionality of the matrix. We assert that $\boldsymbol{\Omega}_{k,l}$ is a
rank one matrix if and only if the inequality
$\varpi_{k,l}\mathbf{I}_{M-1}-\mho^{T}\boldsymbol{\Omega}_{k,l}\mho\succeq\mathbf{0}$
holds true, where $\varpi_{k,l}=0$, $\mathbf{I}_{M-1}$ denotes an identity
matrix with dimension $M-1$, and $\mho\in\mathbb{R}^{M\times(M-1)}$
corresponds to the eigenvectors corresponding to the $M-1$ smallest
eigenvalues of $\boldsymbol{\Omega}_{k,l}$.

###### Proof 2

See Appendix B.

According to proposition 1, we can replace the rank constraint (19a) with a
positive semi-definite constraint:

$\varpi_{k,l}\mathbf{I}_{M-1}-\mho^{{(q)}^{T}}\boldsymbol{\Omega}_{k,l}\mho^{(q)}\succeq\mathbf{0},\forall
k\in\mathcal{K},\forall l\in\mathcal{L}.$ (20)

Given the unavailability of $\mho^{{(q)}^{T}}$, we resort to the SCA method
and employ the smallest eigenvectors of $\boldsymbol{\Omega}_{k,l}$. To
achieve the goal of ultimately having $\varpi_{k,l}=0$ while facilitating the
attainment of an initial feasible point, we introduce a penalty term for
$\varpi_{k,l}$ in the objective function. Thus, at SCA iteration $(q)$, the
following convex problem is addressed:

$\displaystyle\text{P}_{4}$
$\displaystyle:\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l},\varpi_{k,l}}\>\>\>\tilde{\mathscr{\eta}}_{\mathrm{eff}}(\boldsymbol{\Omega}_{k,l})-\varsigma^{(q)}\sum_{l\in\mathcal{L}}\sum_{k\in\mathcal{K}}\varpi_{k,l}$
$\displaystyle s.t.$
$\displaystyle:~{}\eqref{p2-2},\eqref{p2_4},\eqref{p2_5},\eqref{p3-c1},\eqref{rank_one}.$

Note that the penalty factor’s updates follow an interior iterative sub-
algorithm with the formula
$\varsigma^{(q)}=\min(\nu\varsigma^{(q-1)},\vartheta_{\max})$, where $\nu$ is
a positive coefficient factor, and $\vartheta_{\max}$ is the maximum allowable
penalty factor.

The objective function in $\text{P}_{4}$ is a fractional function, rendering
the optimization problem to be still non-convex. To tackle this issue, we
adopt Dinkelbach method [19], with $(d)$ as its iteration, to address the
fractional form of the optimization problem. Denoting
$\tilde{\mathscr{\eta}}_{\mathrm{eff}}^{*}(\boldsymbol{\Omega}_{k,l})$ as the
optimal energy efficiency point within the feasible solution set defined by
its constraints, we then transform the problem into a non-fractional
optimization problem for further resolution, i.e.,

$\displaystyle\text{P}_{5}$
$\displaystyle:\\!\\!\\!\\!\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l},\varpi_{k,l}}\\!\\!\sum\limits_{k\in\mathcal{K}}\\!\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{(d)})\\!-\\!\varrho^{(d)}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{(d)})\\!-\\!\varsigma^{(d)}\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}$
$\displaystyle s.t.$
$\displaystyle:~{}\eqref{p2-2},\eqref{p2_4},\eqref{p2_5},\eqref{p3-c1},\eqref{rank_one},$

where
$\varrho^{(d)}=\underset{\boldsymbol{\Omega}_{k,l}^{(d)}}{{\max}}~{}\tilde{\mathscr{\eta}}_{\mathrm{eff}}^{(d)}(\boldsymbol{\Omega}_{k,l}^{(d)})$,
and $\varsigma^{(d)}$ is the new penalty factor update rules
$\varsigma^{(d)}=\min(\nu\varsigma^{(d-1)},\vartheta_{\max})\times\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{(d)})$.
Finally, the convexity of the objective function in $\text{P}_{5}$ with
respect to the active beamforming variables can be demonstrated through a
formal proof based on the following proposition.

Algorithm 1 Iterative SCA-based Resource Allocation Algorithm

0: Set ${d}=0$ for the Dinkelbach procedure, set the maximum number of
iteration ${D}_{\max}$, initialize the beamformer matrix
$\boldsymbol{\Omega}_{k,l}=\boldsymbol{\Omega}_{k,l}^{0}$, initialize
$\varsigma^{0}$, set $\vartheta_{\max}\gg 1$, $\nu>1$, and set the tolerance
$\varepsilon=10^{-3}$.

1: repeat

2: Calculate
$\tilde{\mathcal{R}}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$ for
a given $\boldsymbol{\Psi}_{l}$.

3: Calculate the rank-one relaxation constraint in (20).

4: Solve $\text{P}_{5}$ for $\varrho^{(d-1)}$.

5: if
$\big{|}\sum\limits_{k\in\mathcal{K}}\\!\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{(d)})-\varrho^{(d-1)}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{(d)})$
$-\varsigma^{(d)}\sum\limits_{l\in\mathcal{L}}\sum\limits_{k\in\mathcal{K}}\varpi_{k,l}\big{|}\\!\leq\\!\varepsilon$

6: return $\boldsymbol{\Omega}_{k,l}=\boldsymbol{\Omega}_{k,l}^{(d)}$,
$\varrho^{*}=\varrho^{(d-1)}$.

7: else
$\varrho^{(d)}=\tilde{\mathscr{\eta}}_{eff}^{(d)}(\boldsymbol{\Omega}_{k,l}^{(d)})$,
end if.

8: Update $\varsigma^{(d)}=\min(\nu\varsigma^{(d-1)},\vartheta_{\max})$.

9: ${d}\leftarrow{d}+1$.

10: until ${d}={D}_{{\max}}$.

11: return $\boldsymbol{\Omega}_{k,l}=\boldsymbol{\Omega}_{k,l}^{*}$.

###### Proposition 2

The optimal energy efficiency, denoted as
$\tilde{\mathscr{\eta}}_{eff}^{*}(\boldsymbol{\Omega}_{k,l}^{*})$, serves as a
means to derive the resource allocation policy if and only if

$\displaystyle\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l}}\>$
$\displaystyle\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{(d)})-\varrho^{(d)}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{(d)}-\\!\varsigma^{(q)}\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}=$
$\displaystyle\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{*})-\varrho^{*}~{}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{*})-\varsigma^{(q)}\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}=0,$
(23)

for
$\sum_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{*})\geq
0$ and $\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{*})\geq 0$, where
$\boldsymbol{\Omega}_{k,l}^{*}$ provides the optimal solution to
$\text{P}_{5}$.

###### Proof 3

See Appendix C.

The optimization problem $\text{P}_{5}$ is now a convex problem and can be
efficiently solved by standard convex optimization solvers such as CVX [20,
21, 22]. Ultimately, we detail our proposed algorithm in Algorithm 1.

###  III-B Sub-problem (2): Optimizing $\boldsymbol{\Psi}_{l}$ with given
$\boldsymbol{\omega}_{k,l}$

In the second sub-problem, equipped with the optimal active beamforming
matrices $\boldsymbol{\omega}_{k,l}$ from the preceding sub-problem, we
progress with the optimization of the passive reflecting elements at the IRS,
denoted by $\boldsymbol{\Psi}_{l}$. Given that the semi-definite matrix
$\boldsymbol{\Omega}_{k,l}$ is provided, the optimization problem in
$\text{P}_{1}$ shifts towards the maximization of the data-rate. The major
hurdle in optimizing the phase shifts at the IRS arises due to the constraint
(9b). Particularly, the constraint (9b) imposes a unit-modulus constraint,
which presents a significant challenge in the quest to solve the problem. As
such, we start by defining
${\boldsymbol{e}}=(e^{j\phi_{1}},...,e^{j\phi_{N}})^{H}\in\mathbb{C}^{N\times
1}$ and
$\tilde{{\boldsymbol{e}}}=[{\boldsymbol{e}}^{T}\>\psi]^{T}\in\mathbb{C}^{(N+1)\times
1}$, where $\psi\in\mathbb{C}$ is a dummy variable with $|\psi|=1$. To aid in
crafting the solution, we further define
${\mathbf{E}}=\tilde{{\boldsymbol{e}}}\tilde{{\boldsymbol{e}}}^{H}\in\mathbb{C}^{(N+1)\times(N+1)}$.
Consequently, we derive the following:

$\displaystyle\left|({\mathbf{h}}_{\rm{IRS},k}^{H}\boldsymbol{\Psi}_{l}{\mathbf{\mathbf{H}}}+\mathbf{h}_{\rm{BS},k}^{H})\boldsymbol{\Omega}_{k,l}\right|^{2}$
$\displaystyle\triangleq\text{Tr}(\mathbf{E}\mathbf{Z}_{k}\boldsymbol{\Omega}_{k,l}\mathbf{Z}_{k}^{H})$
$\displaystyle=\text{Tr}(\boldsymbol{\Omega}_{k,l}\mathbf{Y}_{k}),\forall
k\in\mathcal{K},$ (24)

where
$\mathbf{Z}_{k}=[(\text{diag}({\mathbf{h}}_{\rm{IRS},k}^{H})\mathbf{H})^{T}\>\>\>\mathbf{h}_{\rm{BS},k}^{*}]^{T},\mathbf{Y}_{k}=\mathbf{Z}_{k}^{H}\mathbf{E}\mathbf{Z}_{k}.$
In a manner akin to $\text{P}_{3}$, we address the non-convex constraint (9a)
and the objective function. To achieve this, we employ Lemma 1 and
subsequently rewrite the data-rate function as
$\tilde{\mathcal{R}}_{k}(\boldsymbol{\Psi}_{l})=\mathcal{U}_{k}(\boldsymbol{\Psi}_{l})-\mathcal{V}_{k}(\boldsymbol{\Psi}_{l})$.
Now, we restate the optimization problem $\text{P}_{1}$ as follows:

$\displaystyle\text{P}_{6}$
$\displaystyle:\mathop{\max}\limits_{\mathbf{E},\boldsymbol{\Psi}_{l}}\>\>\>\sum\limits_{k\in\mathcal{K}}\tilde{R}_{k}(\boldsymbol{\Psi}_{l})$
$\displaystyle s.t.$
$\displaystyle:\tilde{\mathcal{R}}_{k}(\boldsymbol{\Psi}_{l})\geq
R_{\min,k},\>\forall k\in\mathcal{K},$ (25a)
$\displaystyle~{}~{}{\rm{diag}}(\boldsymbol{\Psi}_{l})=\mathbf{1}_{N+1},~{}\forall
l\in\mathcal{L},$ (25b) $\displaystyle~{}~{}{\rm{rank}}(\mathbf{E})\leq 1,$
(25c) $\displaystyle~{}~{}\boldsymbol{\Psi}_{l}\succeq 0,\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall l\in\mathcal{L},$ (25d)
$\displaystyle~{}~{}\mathbf{E}\succeq\mathbf{0}.$ (25e)

Following a similar approach as in $\text{P}_{6}$, and drawing on the insights
from proposition 1, we can substitute the rank constraint (25c) with a
positive semi-definite constraint:

$\zeta_{l}\mathbf{I}_{N}-\Upsilon^{{(g)}^{T}}\boldsymbol{\Psi}_{l}\Upsilon^{(g)}\succeq\mathbf{0},\forall
l\in\mathcal{L},$ (26)

where $\Upsilon^{(g)}$ is an $N\times(N+1)$ matrix at the $(g)$-th SCA
iteration with its columns corresponding to the smallest $N$ eigenvectors of
$\boldsymbol{\Psi}_{l}$. Furthermore, to sustain a rank of one for
$\boldsymbol{\Psi}_{l}$, equation (26) must hold true, assuming $\zeta_{l}=0$.
Hence, by incorporating $\zeta_{l}$ as a penalty term into the objective
function of $\text{P}_{6}$, we write the subsequent convex problem:

$\displaystyle\text{P}_{7}$
$\displaystyle:\underset{{\mathbf{E}},\boldsymbol{\Psi}_{l}}{\text{max}}\sum\limits_{k\in\mathcal{K}}{\tilde{\mathcal{R}}_{k}(\boldsymbol{\Psi}_{l})}-\varkappa^{(g)}\sum\limits_{l\in\mathcal{L}}\zeta_{l}$
$\displaystyle s.t.$
$\displaystyle:~{}\eqref{p6-1},\eqref{p6-2},\eqref{p6-4},\eqref{p6-5},\eqref{rank_one_2},$

where $\varkappa^{(g)}$ is a sequences of increasing weights. The optimization
problem $\text{P}_{7}$ now can be efficiently solved just as $\text{P}_{5}$
[20, 23]. The solution of these two sub-problems yields the suboptimal
solution of $\text{P}_{1}$.

###### Proposition 3

The objective function of $\text{P}_{1}$ is ensured to be monotonically non-
decreasing throughout the iterations of the proposed algorithm.

###### Proof 4

See Appendix D. $\blacksquare$

Ultimately, the final iterative AO algorithm that solves sub-problems is
presented in Algorithm 2.

Algorithm 2 Proposed Iterative AO algorithm

0: Set $s=0$, maximum number of iteration ${S}_{\text{max}}$, and initialize
the beamformer matrix and phase shifts as
$\boldsymbol{\Omega}_{k,l}=\boldsymbol{\Omega}_{k,l}^{0}$ and
$\boldsymbol{\Psi}_{l}=\boldsymbol{\Psi}_{l}^{0}$, respectively.

1: repeat

2: Solve problem $\text{P}_{5}$ for given $\boldsymbol{\Psi}_{l}^{s}$, and
obtain the optimal solution $\boldsymbol{\Omega}_{k,l}^{s}$ based on Algorithm
1.

3: Solve problem $\text{P}_{7}$ for given $\boldsymbol{\Omega}_{k,l}^{s}$, and
obtain the optimal solution $\boldsymbol{\Psi}_{l}^{s}$.

4: $s\leftarrow s+1$.

5: until ${s}={S}_{\text{max}}$

6: return $\\{\boldsymbol{\Omega}_{k,l}^{{*}},\boldsymbol{\Psi}_{l}^{*}\\}=$
$\\{\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{s}\\}.$

##  IV Numerical Results and Discussions

This section provides simulation results to investigate the efficiency of the
proposed algorithm in IRS-supported downlink URLLC systems, utilizing finite
blocklength codes. A rectangular region with dimensions of $(100,100)$ meters
is considered. The BS is positioned at $(0,0)$ meters, the IRS is located at
$(50,0)$ meters, and all users are scattered randomly within the boundaries of
this rectangle region. The path loss model is expressed as
$35.3+37.6\log_{10}(d_{k})~{}\rm{[dBm]}$, where $d_{k}$ refers to the distance
in kilometers between the BS and the $k$-th URLLC user [8]. The convergence
tolerance for the proposed AO algorithm, which is based on rank relaxation,
SCA, and the Dinkelback method, is set at $10^{-2}$. It is assumed that the
thermal noise density stands at $-174~{}\rm{[dBm/Hz]}$. In addition, the
maximum probability of decoding error for URLLC user $k$, represented as
$\epsilon_{k}$, is set to be $10^{-7}$. Furthermore, the simulation parameters
are configured with $K=4$, $M=5$, and $R_{\min,k}=1.6~{}\rm{[bits/Sec/Hz]}$ ,
in accordance with the methodology outlined in references [10] and [8].

Figure 2: Impact of tolerable decoding error, $\epsilon_{k}$, on the EE.

Fig. 2 depicts the average EE in correlation with diverse maximum decoding
error probability values, denoted as $\epsilon_{k,{\max}}$, keeping the block
code length constant at $250$ symbols, represented as $m_{d}=250$. The plotted
data brings to light an augmenting pattern in EE as the permissible decoding
error broadens. This increment is attributable to the inverse proportionality
of $Q^{-1}(\epsilon_{k})$ and $\epsilon_{k}$. The dominating contributor to
this rise is a diminished rate loss
$\mathcal{V}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$ of (2)
correlating with an enlargement in the maximum acceptable decoding error
probability. This correlation ensures that the minimal data-rate obligation is
met at reduced transmission power, thereby amplifying EE. The figure also
signifies an upsurge in EE when the number of reflecting elements is
increased.

Moreover, the figure also includes a comparative study with two basic models.
Baseline scheme 1 operates under a fixed beamforming policy at the IRS [8],
whereas Baseline scheme 2 considers a system without IRS. Our proposed
methodology’s performance, as demonstrated by the figure, surpasses that of
baseline scheme 2 by $78.9\%$. This superior performance results from the
incorporation of IRS and the harmonized optimization of active and passive
beamforming matrices at the BS and IRS. This strategy displays an improved
performance when compared to baseline scheme 1, $51.7\%$ and $48.8\%$ for
$N=30$ and $N=20$, respectively, thereby underscoring the substantial benefits
of integrating IRS into our system.

Figure 3: Impact of the blocklength, $m_{d}$, on the EE.

Fig. 3 illustrates the impact of the blocklength size, denoted as $m_{d}$, on
the achievable EE at a constant $\epsilon_{k}=10^{-7}$. With an increase in
$m_{d}$, the EE initially experiences a moderate boost, followed by a slow
elevation until it plateaus. In comparison to the scenario with an IRS, the
scenario without it consistently exhibits lower EE, independent of the
quantity of reflecting elements. Smaller $N$ values necessitate higher
transmit power to meet QoS standards, leading to a reduction in EE. This
reduction arises from the amplified interference in the data-rate function due
to multi-user interference in the SINR function, resulting in an increase in
transmit power, a decrease in data-rate, and consequently negatively pruned
EE.

##  V Conclusion

This paper investigated the resource allocation for a downlink multiuser
intelligent reflecting surface (IRS)-aided ultra-reliable low-latency
communication (URLLC) system. In particular, the resource allocation design
with active/passive beamforming was formulated to maximize the energy
efficiency (EE) while considering URLLC Quality of Service (QoS) requirements
for each user with adopting the short packet transmission. Subsequent to this,
we brought forth an innovative iterative rank relaxation and adopted the
successive convex approximation (SCA) approach, as well as a penalty-based
method to solve each sub-problem. Simulation results demonstrated the efficacy
of the proposed method compared to a system without or with unoptimized IRSs.
Building upon these findings, our research could serve as a foundation for
further studies on rank-constrained resource allocation strategies in
comparable IRS-aided URLLC systems.

##  VI Acknowledgement

This work was supported by the CHIST-ERA grant SAMBAS (CHIST-ERA20-SICT-003),
with funding from FWO, ANR, NKFIH, and UKRI, and, in part, by the framework of
the Recovery Plan, Transformation and Resilience (UNICO I+D 5G 2021, nr.
TSI-063000-2021-6-Open6G Joint Open 6G Communications and Sensing), funded by
the Spanish Ministry of Economic Affairs and Digital Transformation and
European Union - NextGeneration EU.

##  Appendix A Proof of the Lemma 1

According to the data-rate equation in (2), we can drive a lower bound on
$\mathcal{U}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$ and an
upper bound on
$\mathcal{V}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$. In order
to approximate the non-convex logistic function, we utilize the Majorization
Minimization (MM) algorithm [24, 25, 26]. This is achieved by constructing a
surrogate function using a first-order Taylor approximation, which is
expressed as follows:

$\displaystyle\mathcal{U}_{k}(\textbf{a})$
$\displaystyle\simeq\mathcal{U}_{k}(\textbf{a}^{(j)})+\nabla_{\textbf{a}}\mathcal{U}_{k}(\textbf{a}^{(j)}).(\textbf{a}-\textbf{a}^{(j)})\triangleq\tilde{\mathcal{U}}_{k}(\textbf{a}),$
(28) $\displaystyle\mathcal{V}_{k}(\textbf{b})$
$\displaystyle\simeq\mathcal{V}_{k}(\textbf{b}^{(j)})+\nabla_{\textbf{b}}\mathcal{V}_{k}(\textbf{b}^{(j)}).(\textbf{b}-\textbf{b}^{(j)})\triangleq\tilde{\mathcal{V}}_{k}(\textbf{b}),$
(29)

where
$\textbf{a},\textbf{b}=\\{\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}\\}$
and $j$ is the iteration number. Moreover, $\textbf{a}^{(j)}$ and
$\textbf{b}^{(j)}$ denote the solutions of the problem at $(j)^{{th}}$
iteration. Based on (12)-(16), we first rewrite the SINR formula as follows:

$\Gamma_{k,l}=\frac{\left|\mathfrak{a}_{k,k,l}\right|^{2}}{\mathfrak{b}_{k,l}-\left|\mathfrak{a}_{k,k,l}\right|^{2}},\forall
k\in\mathcal{K},\forall l\in\mathcal{L}.\vspace*{2mm}$ (30)

Equivalently, we can rewrite the channel dispersion as:

$\displaystyle\Delta_{k,l}=(\log_{2}e)^{2}\left(1-\left(\frac{\mathfrak{b}_{k,l}-\left|\mathfrak{a}_{k,k,l}\right|^{2}}{\mathfrak{b}_{k,l}}\\!\right)^{2}\right),$
$\displaystyle\forall k\in\mathcal{K},\forall l\in\mathcal{L}.$ (31)

Subsequently, we have:

$\displaystyle\vspace*{1mm}\mathcal{U}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
$\displaystyle=-\sum\limits_{l\in\mathcal{L}}\log_{2}\left(1-\frac{\left|\mathfrak{a}_{k,k,l}\right|^{2}}{\mathfrak{b}_{k,l}}\right)$
$\displaystyle=\sum\limits_{l\in\mathcal{L}}\mathcal{U}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}),~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},$ (32)
$\displaystyle\mathcal{V}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})$
$\displaystyle=\sum\limits_{l\in\mathcal{L}}\beta_{k}\sqrt{\Delta_{k,l}}$
$\displaystyle=\sum\limits_{l\in\mathcal{L}}\mathcal{V}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l}),~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K}.$ (33)

Thus, employing the MM technique allows us to derive the following lower limit
for $\mathcal{U}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})$:

$\displaystyle\mathcal{U}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})\geq$
$\displaystyle-\frac{\left|\mathfrak{a}_{k,k,l}^{(j)}\right|^{2}}{\mathfrak{b}_{k,l}^{(j)}\left(\mathfrak{b}_{k,l}^{(j)}-\left|\mathfrak{a}_{k,k,l}^{(j)}\right|^{2}\right)}\mathfrak{b}_{k,l}$
$\displaystyle+2\Re\left\\{\frac{\mathfrak{a}_{k,k,l}^{(j)}\mathfrak{a}_{k,k,l}^{H}}{\mathfrak{b}_{k,l}^{(j)}-\left|\mathfrak{a}_{k,k,l}^{(j)}\right|^{2}}\right\\}+\xi_{k,l}^{\mathcal{U}},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (34)

where

$\xi_{k,l}^{\mathcal{U}}=\log_{2}(1+\Gamma_{k,l}^{(j)})-\frac{|\mathfrak{a}_{k,k,l}^{(j)}|^{2}}{\mathfrak{b}_{k,l}^{(j)}(\mathfrak{b}_{k,l}^{(j)}-|\mathfrak{a}_{k,k,l}^{(j)}|^{2})}.$
(35)

Similarly, we use the first-order Taylor approximation to derive the following
upper limit for
$\mathcal{V}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})$:

$\displaystyle\mathcal{V}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})\leq\beta_{k}\sqrt{\Delta_{k,l}^{j}}+\frac{(\log_{2}e)^{2}\beta_{k}}{2\sqrt{\Delta_{k,l}^{j}}}$
$\displaystyle\left(\\!1-\frac{1}{\left(1+\Gamma_{k,l}\right)^{2}}\\!\right),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\forall k\in\mathcal{K},\forall
l\in\mathcal{L}.$ (36)

Now, by defining $\Xi_{k,l}\triangleq 1/(1+\Gamma_{k,l})$ and according to
[8], we can rewrite:

$\displaystyle\frac{1}{\left(1+\Gamma_{k,l}\right)^{2}}$
$\displaystyle=\Xi_{k,l}^{2}\geq
2\Xi_{k,l}^{(j)}\Xi_{k,l}-(\Xi_{k,l}^{(j)})^{2}$
$\displaystyle=\frac{2}{\left(1+\Gamma_{k,l}^{(j)}\right)\left(1\\!+\\!\Gamma_{k,l}\right)}-\frac{1}{\left(1\\!+\\!\Gamma_{k,l}^{(j)}\right)^{2}},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (37)

where the inequality holds due to convexity of $(\Xi_{k,l})^{2}$. The term
$1/(1+\Gamma_{k,l})$ in (37) is non-convex due to its fractional form. We
rewrite this term as follows and then apply MM:

$\displaystyle\frac{1}{1+\Gamma_{k,l}}=\frac{\mathfrak{b}_{k,l}-\left|\mathfrak{a}_{k,k,l}\right|^{2}}{\mathfrak{b}_{k,l}}\geq-\frac{\left(\mathfrak{b}_{k,l}^{(j)}-\big{|}\mathfrak{a}_{k,k,l}^{(j)}\big{|}^{2}\right)\mathfrak{b}_{k,l}}{\big{|}\mathfrak{b}_{k,l}^{(j)}\big{|}^{2}}$
$\displaystyle~{}~{}~{}~{}~{}+2\Re\left\\{\frac{\sum\limits_{i\neq
k,i\in\mathcal{K}}\mathfrak{a}_{i,k,l}^{(j)}\mathfrak{a}_{i,k,l}^{H}+\sigma^{2}_{k}}{\mathfrak{b}_{k,l}^{(j)}}\right\\},\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ (38)

where, again, the inequality holds due to the convexity with respect to the
$\mathfrak{a}_{i,k,l}$ and $\mathfrak{b}_{k,l}$. Thus, we derive an upper
bound as:

$\displaystyle\begin{split}\mathcal{V}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})&\leq\frac{(\log_{2}e)^{2}\beta_{k}}{\sqrt{(\Gamma^{(j)}_{k,l})^{2}+2\Gamma^{(j)}_{k,l}}}\left(\frac{\left(\mathfrak{b}_{k,l}^{(j)}-\big{|}\mathfrak{a}_{k,k,l}^{(j)}\big{|}^{2}\right)\mathfrak{b}_{k,l}}{\big{|}\mathfrak{b}_{k,l}^{(j)}\big{|}^{2}}\right.\\\
&\left.~{}-2\Re\left\\{\frac{\sum\limits_{i\neq
k,i\in\mathcal{K}}\mathfrak{a}_{i,k,l}^{(j)}\mathfrak{a}_{i,k,l}^{H}+\sigma^{2}_{k}}{\mathfrak{b}_{k,l}^{(j)}}\right\\}\right)+\xi_{k,l}^{\mathcal{V}},\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{K},\forall l\in\mathcal{L},\end{split}$ (39)

where

$\displaystyle\xi_{k,l}^{\mathcal{V}}$
$\displaystyle=\frac{(2\Delta_{k,l}^{(j)}+(\log_{2}e)^{2})\beta_{k}}{2\sqrt{\Delta_{k,l}^{(j)}}}$
$\displaystyle+\frac{(\log_{2}e)^{2}\beta_{k}}{2(1+\Gamma^{(j)}_{k,l})\sqrt{(\Gamma^{(j)}_{k,l})^{2}+2\Gamma^{(j)}_{k,l}}},\forall
k\in\mathcal{K},\forall l\in\mathcal{L}.$ (40)

Ultimately, the lower bound of the data-rate function in (2) can be
approximated as:

$\tilde{\mathcal{R}}_{k}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi}_{l})\triangleq\sum\limits_{l\in\mathcal{L}}\big{(}\mathcal{U}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})-\mathcal{V}_{k,l}(\boldsymbol{\omega}_{k,l},\boldsymbol{\Psi})\big{)},\forall
k\in\mathcal{K}.$ (41)

This completes the proof.

##  Appendix B Proof of the Proposition 1

Assuming the nonnegative eigenvalues of $\boldsymbol{\Omega}_{k,l}$ are
arranged in descending order as $[\kappa_{M};\kappa_{M-1};...;\kappa_{1}]$, we
can exploit the relationship between an eigenvector’s Rayleigh quotient and
its associated eigenvalue. Consequently, the matrix
$\varpi_{k,l}\mathbf{I}_{M-1}-\mho^{T}\boldsymbol{\Omega}_{k,l}\mho,$ $\forall
k\in\mathcal{K},\forall l\in\mathcal{L},$ takes on the form of a diagonal
matrix, with its diagonal elements set as
$[\varpi_{k,l}-\kappa_{M-1};\varpi_{k,l}-\kappa_{M-2};...;\varpi_{k,l}-\kappa_{1}]$.
In light of this, we observe that $\boldsymbol{\Omega}_{k,l}$ possesses $M-1$
smallest eigenvalues all being zero if and only if the conditions
$\varpi_{k,l}\mathbf{I}_{M-1}-\mho^{T}\boldsymbol{\Omega}_{k,l}\mho\succeq\mathbf{0}$
and $\varpi_{k,l}=0,$ $\forall k\in\mathcal{K},\forall l\in\mathcal{L},$ are
satisfied simultaneously. Consequently, $\boldsymbol{\Omega}_{k,l}$ qualifies
as a rank one matrix under these conditions [18].

##  Appendix C Proof of the Proposition 2

Let’s denote $\varrho^{*}$ be the optimal solution, corresponding to the
optimal resource allocation policy $\boldsymbol{\Omega}^{*}_{k,l}$ of the
objective function in $\text{P}_{5}$, that is:

$\displaystyle\varrho^{*}=\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l}}\>\frac{\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})}{\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})}-\varsigma^{(q)}\sum_{l\in\mathcal{L}}\sum_{k\in\mathcal{K}}\varpi_{k,l}.$
(42)

Therefore, the optimal EE must satisfy the following inequality:

$\displaystyle\varrho^{*}=$
$\displaystyle\frac{\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}^{*}_{k}(\boldsymbol{\Omega}^{*}_{k,l})}{\mathcal{E}^{*}_{\rm{tot}}(\boldsymbol{\Omega}^{*}_{k,l})}-\varsigma^{(q)}\sum_{l\in\mathcal{L}}\sum_{k\in\mathcal{K}}\varpi_{k,l}$
$\displaystyle\geq$
$\displaystyle\frac{\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})}{\mathcal{E}^{(d)}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})}-\varsigma^{(q)}\sum_{l\in\mathcal{L}}\sum_{k\in\mathcal{K}}\varpi_{k,l}.$
(43)

It becomes readily apparent that:

$\displaystyle\sum\limits_{k\in\mathcal{K}}\\!\\!\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})\\!-\\!\\!\varrho^{*}\\!\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})\\!-\\!\varsigma^{(q)}\\!\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})\\!\leq\\!0,$
(44)

Thus, it can be inferred that:

$\displaystyle\sum\limits_{k\in\mathcal{K}}\\!\\!\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l}^{*})\\!-\\!\varrho^{*}\\!\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{*})\\!-\\!\\!\varsigma^{(q)}\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l}^{*})\\!=\\!0,\\!$
(45)

Therefore, we have:

$\displaystyle\mathop{\max}\limits_{\boldsymbol{\Omega}_{k,l}}\>$
$\displaystyle\sum\limits_{k\in\mathcal{K}}\tilde{\mathcal{R}}_{k}(\boldsymbol{\Omega}_{k,l})-\varrho^{*}~{}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\>-\varsigma^{(q)}\\!\sum_{l\in\mathcal{L}}\\!\sum_{k\in\mathcal{K}}\\!\varpi_{k,l}\mathcal{E}_{\rm{tot}}(\boldsymbol{\Omega}_{k,l})=0,$
(46)

and this can be attained through the resource allocation policy. This
concludes the proof. $\blacksquare$

##  Appendix D Proof of the Proposition 3

Let’s denote the objective functions of $\text{P}_{1}$, $\text{P}_{5}$, and
$\text{P}_{7}$ respectively as $\mathcal{\daleth}_{\text{P}_{1}}$,
$\mathcal{\daleth}_{\text{P}_{5}}$, and $\mathcal{\daleth}_{\text{P}_{7}}$.
Also, in consideration of
$\\{\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{s}\\}$ and
$\\{\boldsymbol{\Omega}_{k,l}^{s-1},\boldsymbol{\Psi}_{l}^{s-1}\\}$ as the
feasible solutions of $\text{P}_{1}$ in the $s$-th and $(s-1)$-th iterations
respectively, we can establish the following inequalities:

$\displaystyle\hskip
17.07164pt\mathcal{\daleth}_{\text{P}_{1}}(\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{s})$
$\displaystyle=\mathcal{\daleth}_{\text{P}_{7}}(\boldsymbol{\Omega}_{k,l}^{*},\boldsymbol{\Psi}_{l}^{s})$
$\displaystyle\geq\mathcal{\daleth}_{\text{P}_{7}}(\boldsymbol{\Psi}_{l}^{s-1})$
$\displaystyle=\mathcal{\daleth}_{\text{P}_{1}}(\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{s-1}),$
(47)
$\displaystyle\mathcal{\daleth}_{\text{P}_{1}}(\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{s-1})$
$\displaystyle=\mathcal{\daleth}_{\text{P}_{5}}(\boldsymbol{\Omega}_{k,l}^{s},\boldsymbol{\Psi}_{l}^{*})$
$\displaystyle\geq\mathcal{\daleth}_{\text{P}_{5}}(\boldsymbol{\Omega}_{k,l}^{s-1})$
$\displaystyle=\mathcal{\daleth}_{\text{P}_{1}}(\boldsymbol{\Omega}_{k,l}^{s-1},\boldsymbol{\Psi}_{l}^{s}).$
(48)

By utilizing inequalities (47) and (48), we can guarantee the improvement in
the value of the objective function of $\text{P}_{1}$ after every iteration.
$\blacksquare$

## References

  * [1] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” IEEE Commun. Mag., vol. 58, pp. 106–112, Jan. 2020.
  * [2] M. Bennis, M. Debbah, and H. V. Poor, “Ultrareliable and low-latency wireless communication: Tail, risk, and scale,” Proc. IEEE, vol. 106, pp. 1834–1853, Oct. 2018.
  * [3] Y. Polyanskiy, H. V. Poor, and S. Verdu, “Channel coding rate in the finite blocklength regime,” IEEE Trans. Inf. Theory, vol. 56, pp. 2307–2359, May 2010.
  * [4] R. Hashemi, S. Ali, N. H. Mahmood, and M. Latva-aho, “Average rate and error probability analysis in short packet communications over RIS-aided URLLC systems,” IEEE Trans. Veh. Tech., vol. 70, pp. 10320–10334, Aug. 2021.
  * [5] S. Dhok, P. Raut, P. K. Sharma, K. Singh, and C.-P. Li, “Non-linear energy harvesting in RIS-Assisted URLLC networks for industry automation,” IEEE Trans. Commun., vol. 69, pp. 7761–7774, Nov. 2021.
  * [6] H. Xie, J. Xu, Y.-F. Liu, L. Liu, and D. W. K. Ng, “User grouping and reflective beamforming for IRS-Aided URLLC,” IEEE Wirel. Commun. Lett., vol. 10, pp. 2533–2537, Nov. 2021.
  * [7] Z. Li, H. Shen, W. Xu, P. Zhu, and C. Zhao, “Resource allocation for irs-assisted uplink URLLC systems,” IEEE Commun. Lett., vol. 27, pp. 1540–1544, Jun. 2023.
  * [8] J. Jalali, A. Khalili, A. Rezaei, R. Berkvens, M. Weyn, and J. Famaey, “IRS-based energy efficiency and admission control maximization for IoT users with short packet lengths,” IEEE Trans. Veh. Technol., vol. 72, pp. 12379–12384, Sept. 2023.
  * [9] V. D. P. Souto, S. Montejo-Sánchez, J. L. Rebelatto, R. D. Souza, and B. F. Uchôa-Filho, “IRS-Aided physical layer network slicing for URLLC and eMBB,” IEEE Access, vol. 9, pp. 163086–163098, Dec. 2021.
  * [10] J. Jalali, A. Rezaei, A. Khalili, and J. Famaey, “Power-efficient joint resource allocation and decoding error probability for multiuser downlink MISO with finite block length codes,” in 2022 Proc. IEEE 25th Int. Symp. Wirel. Personal Multimed. Commun. (WPMC), pp. 232–237, Oct. 2022.
  * [11] H. Ren, K. Wang, and C. Pan, “Intelligent reflecting surface-aided URLLC in a factory automation scenario,” IEEE Trans. Commun., vol. 70, pp. 707–723, Jan. 2022.
  * [12] T. Bai, C. Pan, Y. Deng, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Latency minimization for intelligent reflecting surface aided mobile edge computing,” IEEE J. Sel. Areas Commun., vol. 38, pp. 2666–2682, Nov. 2020.
  * [13] M. Darabi, W. R. Ghanem, V. Jamali, L. Lampe, and R. Schober, “Active IRS design for RSMA-based downlink URLLC transmission,” in 2023 Proc. IEEE Wirel. Commun. Netw. Conf. (WCNC), pp. 1–6, Mar. 2023.
  * [14] J. Jalali, A. Khalili, R. Berkvens, and J. Famaey, “Joint offloading policy and resource allocation in IRS-aided MEC for IoT users with short packet transmission,” in 2023 Proc. IEEE Veh. Technol. Conf. (VTC2023-Fall), pp. 1–7, Oct. 2023.
  * [15] W. R. Ghanem, V. Jamali, Y. Sun, and R. Schober, “Resource allocation for multi-user downlink MISO OFDMA-URLLC systems,” IEEE Trans. Commun., vol. 68, pp. 7184–7200, Nov. 2020.
  * [16] J. Jalali, A. Khalili, A. Rezaei, and J. Famaey, “Is active IRS useful for mmwave wireless networks or not?,” in 2023 Proc. IEEE Int. Conf. Computing, Netw. Commun. (ICNC), pp. 377–382, Feb. 2023.
  * [17] Y. Polyanskiy, H. V. Poor, and S. Verdu, “Channel coding rate in the finite blocklength regime,” IEEE Trans. Inf. Theory, vol. 56, pp. 2307–2359, May 2010.
  * [18] C. Sun and R. Dai, “An iterative rank penalty method for nonconvex quadratically constrained quadratic programs,” SIAM J. Control Optim., vol. 57, pp. 3749–3766, Jan. 2019.
  * [19] W. Dinkelbach, “On nonlinear fractional programming,” J. Manag. Sci., vol. 13, no. 3, p. 492–498, 1967.
  * [20] S. Boyd and L. Vandenberghe, Convex Optimization.  Cambridge University Press, Mar. 2004.
  * [21] J. Jalali, A. Khalili, H. Tabassum, R. Berkvens, J. Famaey, and W. Saad, “On the energy efficiency of THz-NOMA enhanced UAV cooperative network with SWIPT,” Sept. 2023, [Online]. Available: https://arxiv.org/pdf/2309.13836.pdf.
  * [22] J. Jalali and A. Khalili, “Optimal resource allocation for MC-NOMA in SWIPT-enabled networks,” IEEE Commun. Lett., vol. 24, pp. 2250–2254, Oct. 2020.
  * [23] J. Jalali, A. Khalili, A. Rezaei, J. Famaey, and W. Saad, “Power-efficient antenna switching and beamforming design for multi-user swipt with non-linear energy harvesting,” in 2023 Proc. IEEE 20th Consumer Commun. Netw. Conf. (CCNC), pp. 746–751, Jan. 2023.
  * [24] J. Lee and S. Leyffer, Mixed Integer Nonlinear Programming.  Springer Science & Business Media, Dec. 2011.
  * [25] J. Jalali, Resource allocation for SWIPT in multi-service wireless networks.  M.S. thesis, Dept. Telecommun. Inf. Process., TELIN/IMEC, Ghent Univ., Ghent, Belgium, Jun. 2020.
  * [26] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” IEEE Trans. Signal Process., vol. 65, pp. 794–816, Feb. 2017.


# Visualising and Explaining Deep Learning Models for Speech Quality
Prediction

Henrik Tilkorn, Gabriel Mittag1, Sebastian Möller1,2
1 _Quality and Usability Lab, TU Berlin_
2 _Language Technology, DFKI Berlin_

## Abstract

Estimating quality of transmitted speech is known to be a non-trivial task.
While traditionally, test participants are asked to rate the quality of
samples; nowadays, automated methods are available. These methods can be
divided into: 1) intrusive models, which use both, the original and the
degraded signals, and 2) non-intrusive models, which only require the degraded
signal. Recently, non-intrusive models based on neural networks showed to
outperform signal processing based models. However, the advantages of deep
learning based models come with the cost of being more challenging to
interpret. To get more insight into the prediction models the non-intrusive
speech quality prediction model NISQA is analyzed in this paper. NISQA is
composed of a convolutional neural network (CNN) and a recurrent neural
network (RNN). The task of the CNN is to compute relevant features for the
speech quality prediction on a frame level, while the RNN models time-
dependencies between the individual speech frames. Different explanation
algorithms are used to understand the automatically learned features of the
CNN. In this way, several interpretable features could be identified, such as
the sensitivity to noise or strong interruptions. On the other hand, it was
found that multiple features carry redundant information.

## Introduction

A common method to rate the quality of a speech sample is the evaluation by
test participants using a 5-point scale, where one corresponds to the lowest
and five to the highest quality. These ratings are averaged over all
participants to acquire the mean opinion score (MOS). This procedure is
inconvenient since it is very time and cost-intensive.
To overcome these disadvantages, multiple automated methods have been
developed and applied successfully in practice. With the bandwidth extension
to super-wideband, these methods additionally had to adapt to this more
challenging setting. The speech quality assessment model NISQA proposed by
Mittag and Möller [1] is able to estimate the MOS value of super-wideband
speech transmission in a non-intrusive fashion.

## Methodology

NISQA receives as inputs the degraded signals in form of a spectrogram. This
spectrogram is divided into multiple shorter segments before being processed.
To understand which features of the input have a high influence on the
decision making several different techniques were used.

### Explanation approaches

Occlusion Sensitivity is a straightforward approach, where small regions of
the input are masked. If the classification score, namely the estimated speech
quality, for the occluded input is lower than the score of the original, the
masked part has a positive influence on the overall quality. On the other hand
is an increasing classification score a sign of a hidden negative influence.
To get a complete overview of all the influences of different regions the mask
is gradually moved over the input. To use Occlusion Sensitivity a few design
choices have to be made. The size and the shape of the mask have to be
determined and secondly, the stride, which is how far the mask gets shifted
against its previous location, has to be set. To choose a proper size and
shape it is beneficial to include prior knowledge about the inputs. In this
case, spectrograms were scrutinized and therefore it is known that the width-
axis corresponds to time and the height-axis displays the presence of certain
frequencies. In general, it is possible that disruptions that lead to a
decrease of the speech quality occur only at a particular time and therefore
would be detected by a vertically oriented mask, or occur only on a certain
frequency and then could be detected by a horizontal mask. The size of the
stride mainly controls how fine-grained the overall result is. If the stride
is chosen smaller than the mask size the masks of consecutive inputs are
overlapping and an averaging step is needed to get the interpretable results.
Despite being an intuitive approach, Occlusion Sensitivity has a shortcoming
that is hard to avoid. To get a detailed enough explanation, the size and the
stride have to be chosen relatively small, which leads, in combination with
rather large input images, to too many masked versions to be evaluated. Since
disturbances can occur in different shapes it is also necessary to use
variously shaped and oriented masks which again leads to many evaluations.
While for Occlusion Sensitivity no knowledge about the actual model was
needed, the following methods are explicitly designed to explain the decisions
of neural networks.
DeepLIFT (DL) [2] is a backpropagation-based approach to analyze the decisions
of neural networks. To estimate the relevance scores the model receives the
sample of interest and a neutral reference input. The reference input has to
be chosen in such a way that the model can’t detect any features. In the
context of speech quality assessment, an empty spectrogram will suffice, since
it has no features that have a harmful or beneficial influence on the speech
quality. The model will produce a prediction for the sample and the neutral
reference each. The difference between these two predictions is called
difference-from-reference. To obtain the actual explanation, the difference-
from-reference is propagated back through the model using several propagation
rules.
Another explanation algorithm is Integrated Gradients (IG) proposed by
Sundararajan et al. in 2017 [3]. Many gradient-based explanations suffer from
several shortcomings, some of which are discussed in the section ’Challenges
of explanations’. To overcome these shortcomings IG was designed following the
axioms of sensitivity and implementation invariance, which are also described
in the following section. Similar to DL, IG uses a neutral reference input
called baseline. The model can be understood as a function mapping from the
input to an output space. The explanation is then obtained by calculating the
path-integral of the gradient along the straight path between the input sample
and the baseline. Since calculating the gradient directly is often difficult
or not possible a Riemann approximation is used. Further details can be found
in the the original work.
While DL and IG aim to create explanations on the input level, the Conductance
approach takes a closer look at the hidden units of the model. Conductance
[4], proposed by Dhamdhere et al., is an extension to IG and can be understood
as the flow of IG values through a hidden unit. The goal is to identify which
patterns of the input cause an activation of particular units.

### Challenges of explanations

It is often difficult to distinguish between not understandable reasoning of a
complex model and a bad explanation method. Therefore, it is helpful to assure
some properties by design and anticipate well known problems from the
beginning. Some of these problems and properties will be elaborated hereafter.
The thresholding problem can occur if the gradient is used as an indicator for
relevance. The gradient of some functions like ReLU is not continuous around
their thresholding point. Consequently, the assigned relevance for two
arbitrary close points can vary significantly.
The saturation problem can arise if a function depends on multiple inputs. If
one of these inputs is sufficiently large (or small) it can saturate the
function, i.e. the function reaches its maximal output value by only having
this one input. A backpropagation-based explanation method may now allocate
all the relevance to this dominant input and disregard the other inputs,
despite their potential influences.
Sensitivity is one of the axioms used to design Integrated Gradients. A model
is sensitive if (a) every difference between input and baseline that leads to
a change in classification is assigned a non-zero relevance and (b) every
variable that doesn’t influence the classification has a relevance score of
zero.
The second axiom used for IG is implementation invariance. It states that
functional equivalent models should yield the same explanation. Two models are
functionally equivalent if they produce the same output given the same input,
independent of the underlying implementation.
The completeness criterion assures that relevance scores attributed to the
input sum up to exactly the difference from reference in the output layer. In
other words, the amount of relevance is calculated in the forward pass and on
the backward pass no relevance is disregarded or artificially added.
The advanced explanation methods introduced before account for these problems
and fulfill these properties, except for DL not being implementation
invariant.

## Experiments

The experiments were conducted on the “SwissQual P.OLQA SWB 503” dataset [5].
It consists of 216 samples with 54 conditions. The samples from each condition
share a specific kind of disturbance. All samples belong to a test set and
therefore weren’t used during training.
Occlusion Sensitivity can be applied straightforwardly. The spectrograms from
the dataset are images with a height of 48px and a width of 1300px. The masks
were applied before the spectrograms were divided into smaller segments. To
discover as many reasons for decreased speech quality as possible masks with
different shapes and orientations were used. To detect very local events a
square mask of size 24x24px was used. To find temporal or frequency-related
disturbances horizontal masks with different heights (2,4,6,8px) and vertical
filters with a width of 13px were used respectively. For all these experiments
a stride identical to the mask size was used. The model then estimates a MOS
value for every masked version. To create the actual explanation, it is
necessary to keep track of the positions of the mask. Each region can then be
colored according to the influence it had while being masked.
The remaining methods, namely DL, IG and Conductance, are all implemented in
the Captum library [6]. Captum is a collection of many different explanation
algorithms designed for PyTorch models. These implementations however aim to
evaluate CNN. Since NISQA is consisting of a CNN and a LSTM, only the CNN
part, which is responsible for the feature extraction, can be examined.
NISQA’s CNN receives one segment of the spectrogram at a time and has 20
output channels. To keep the number of visual explanations in a processable
range, two different approaches were used. First, only 30 segments were
examined, but with respect to each of the 20 output channels. Secondly, every
1300 segments of a sample were evaluated, but only with respect to one output
channel at the time. Since these segments are highly overlapping an average
was calculated to get the final explanation.

## Results

Occlusion Sensitivity validates two expected behaviors of NISQA. First, only
non-empty regions affect the estimation of the speech quality. If this was not
true either NISQA or the implementation of Occlusion Sensitivity would be
flawed. The second immediate result is that the prior speech quality of a
sample affects how big the influence of masking some parts is. A sample with
very high quality, for example, a MOS value over 4, is likely to be negatively
influenced by the mask. On the other hand, the quality of a sample with low
quality most likely will be improved by masking. It could also be shown that
hiding noise can improve the estimated quality in some cases. This is not a
contradiction considering how noise influences the perceived quality. It may
be less disturbing if there is a constant background noise instead of a noisy
part sharply interrupted by a short moment of silence. Nevertheless, a
complete lack of noise will lead to a higher quality in general. It is also
worth mentioning that masking strong disturbances, which are clearly visible
on the spectrogram, does not necessarily lead to a quality improvement. Again
it is likely that the introduced interruption has a stronger negative effect
than the previously present perturbations. It remains to emphasize that
Occlusion Sensitivity is a rather simple approach and therefore the results
shouldn’t be over-interpreted. While it is useful to validate some general
assumptions a more fine-grained approach should be used for further insights.

Figure 1: Explanations per segment (top row) w.r.t. different features (each
one row). Only first 5 of 20 features depicted. Explanations were created with
DeepLIFT.

Integrated Gradients and DeepLIFT were used to check which parts of the input
sample contribute to the implicitly learned features of NISQAs CNN. To
identify those features the shape of the explanation is more important
compared to the kind of explanation, namely if it’s a positive or negative
contribution. As exemplary depicted in Figure 1, an explanation for every
segment with respect to every feature was created. It stands out that
explanations for some features are empty while others have red-blue colored
relevance maps. An empty explanation states that the respective feature is not
present in the examined segment. For the explanations of both algorithms, it
holds that some features are nearly always present, while some occur only on a
few samples and some features are never detected in the entire dataset. The
absence of some features can have two reasons. Either they occurred only on
the much larger training dataset and are simply not present in this test set
or during the training process these output nodes were discarded since only a
smaller number of features was necessary to sufficiently learn the problem.
To get a better understanding of what the present features stand for, the
explanations w.r.t. one feature for all segments of a sample were checked. An
example for these averaged explanations can be seen in Figures 2 and 3. While
DL and IG use different approaches to produce their explanations the results
came out qualitatively very similar. Before interpreting what the features
represent, it is to mention that the algorithms in some cases produce
explanations for empty segments. This behavior most likely results from a non-
optimal choice of the baseline. This artifact leads to some horizontal stripes
in the overall explanation due to the average-step.
The clearest interpretation has feature 17, for which an exemplary explanation
is shown in Figure 2. Overall it is rarely active and seems to detect sharp
disruptions in the samples. Since these kinds of disturbances lead to a strong
decrease in the perceived quality, this feature strongly correlates with a low
MOS value. Another feature that is plain to see belongs to the 16th output
node. As shown in Figure 3 it corresponds to the absence of voice,
nevertheless it is insensitive to the presence of noise. While features 16 and
17 are well distinguishable, the remaining features are either similar to
others or not visually clearly interpretable. The features 0, 13 and 19 also
show sensitivity for the regions without any speech but do not leave out the
areas as clearly as feature 16 does. Feature 7 on the other hand reacts to the
presence of speech and thus the opposite of the previously mentioned features,
as also shown in Figure 3.

Figure 2: Spectrogram of an exemplary voice sample (top). Averaged
explanations over all segments w.r.t. feature 17 (bottom). Explanation created
with DeepLIFT.

Besides the capability of detecting certain patterns of the input, some of the
features seem to have a focus on a specific vertical region. Figure 4 shows an
overview of these different orientations. Despite occurring in explanations of
most samples this vertical focus does not generalize to every sample.
Therefore it remains to clarify whether the features are indeed specialized to
certain frequencies or if this behavior is due to some hidden causes.

Figure 3: Spectrogram of an exemplary voice sample (top). Averaged explanation
over all segments w.r.t. feature 16 (middle), feature 7 (bottom). Both
explanations created with DeepLIFT. Figure 4: Cutouts of different
spectrograms and different explanations, showing vertical orientation of
features. From left to right [feature, vertical orientation, algorithm]: (13,
top, DL), (0, top, DL), (11, top; bottom, DL), (19, top; bottom, DL), (5,
center, IG), (8, center; bottom, IG), (3, bottom, DL).

While IG and DL create explanations on the input layer Conductance examines
the reasoning on the hidden convolutional layers of the CNN. In image
classification it is well-known, that the first few layers are responsible to
detect simple features, such as edges, enabling layers deeper in the network
to find more complex shapes. This behavior could not be found in this
experiment. However it is suspected that a similar process happens, but since
the features do not have such a prominent shape as in the object detection
case, they are just not visually detectable by a human observer. Nevertheless,
some patterns in the convolutional layers are recognizable. The explanations
of the first convolutional layer look quite similar to those of the input
layer. The most notable behavior is, as shown in Figure 5, that different
filters of the layer are responsible for different regions of the segment.
This suggests that the model decomposes the problem into several subproblems.

Figure 5: Segment of spectrogram (left) and explanations for different feature
maps of the first convolutional layer w.r.t. feature 0. Explanations created
with Conductance.

In general, the features with a similar explanation from IG and DL tend to
have a similar explanation from Conductance. On the second convolutional
layer, another interesting behavior appears. As depicted in Figure 6 sometimes
different features react to the same input pattern, but with an inverted
explanation. For the following convolutional layers deeper in the network, no
interpretable patterns could be identified. This is most likely due to the
fact that their feature map size is rather small while the representational
power lies within the high number of features maps.

Figure 6: Explanations for one segment of different feature maps of the second
convolutional layer w.r.t feature 8 (top row) and feature 11 (bottom row).
Both features focus on the same input structure, but weight them in opposite
ways.

## Conclusion

The goal of this investigation was to better understand the reasoning of the
speech quality assessment model NISQA. Using the simple Occlusion Sensitivity
approach it was shown that NISQAs basic behavior works as intuitively
expected. The more fine-grained techniques Integrated Gradients and DeepLIFT
were able to identify groups of implicitly learned features that are used to
estimate the speech quality. However, only the CNN component of NISQA was
evaluated. To get a more complete idea of the reasoning process it would be
beneficial to also include the LSTM part of NISQA to acquire an end-to-end
explanation of the estimation. The used techniques also have a set of
hyperparameters. Further fine-tuning may reveal additional details.

## References

  * [1] G. Mittag and S. Möller, “Non-intrusive speech quality assessment for super-wideband speech communication networks,” in ICASSP 2019, 2019.
  * [2] A. Shrikumar, P. Greenside, and A. Kundaje, “Learning important features through propagating activation differences,” in PMLR 2017, 2017.
  * [3] M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic attribution for deep networks,” in PMLR 2017, 2017.
  * [4] K. Dhamdhere, M. Sundararajan, and Q. Yan, “How important is a neuron?,” arXiv preprint arXiv:1805.12233, 2018.
  * [5] ITU-T Rec. P.863, Perceptual objective listening quality assessment.  Geneva, Switzerland: ITU-T, 2018.
  * [6] N. Kokhlikyan, V. Miglani, M. Martin, E. Wang, B. Alsallakh, J. Reynolds, A. Melnikov, N. Kliushkina, C. Araya, S. Yan, and O. Reblitz-Richardson, “Captum: A unified and generic model interpretability library for pytorch,” 2020\.