# Model Order Reduction in Neuroscience

Bülent Karasözen Institute of Applied Mathematics & Department of Mathematics,
Middle East Technical University, Ankara-Turkey<EMAIL_ADDRESS>
###### Abstract

Human brain contains approximately $10^{9}$ neurons, each with approximately
$10^{3}$ connections, synapses, with other neurons. Most sensory, cognitive
and motor functions of our brains depend on the interaction of a large
population of neurons. In recent years, many technologies are developed for
recording large numbers of neurons either sequentially or simultaneously.
Increase in computational power and algorithmic developments have enabled
advanced analyses of neuronal population parallel to the rapid growth of
quantity and complexity of the recorded neuronal activity. Recent studies made
use of dimensionality and model order reduction techniques to extract coherent
features which are not apparent at the level of individual neurons. It has
been observed that the neuronal activity evolves on low-dimensional subspaces.
The aim of model reduction of large-scale neuronal networks is accurate and
fast prediction of patterns and their propagation in different areas of brain.
Spatiotemporal features of the brain activity are identified on low
dimensional subspaces with methods such as dynamic mode decomposition (DMD),
proper orthogonal decomposition (POD), discrete empirical interpolation (DEIM)
and combined parameter and state reduction. In this paper, we give an overview
about the currently used dimensionality reduction and model order reduction
techniques in neuroscience.

Keywords:neuroscience, dimensionality reduction, proper orthogonal
decomposition, discrete empirical interpolation, dynamic mode decomposition,
state and parameter estimation.

Classification[MSC 2010]: 93A15,92C55, 37M10,37M99,37N40,65R32.

##  1 Introduction

Due to the advances in recording and imaging technologies, the number of
recorded signals from brain cells increased significantly in the last few
years. The recorded spatio-temporal neural activity give rise to networks with
complex dynamics. In neuroscience, molecular and cellular level details are
incorporated in large-scale models of the brain in order to reproduce
phenomena such as learning and behavior. The rapid growth of simultaneous
neuronal recordings in scale and resolution brings challenges with the
analysis of the neuronal population activity. New computational approaches
have to be developed to analyze, visualize, and understand large-scale
recordings of neural activity. While algorithmic developments and the
availability of significantly more computing power have enabled analysis of
larger neuronal networks, these advances cannot keep pace with increasing size
and complexity of recorded activity. The activity of complex networks of
neurons can often be described by relatively few distinct patterns. Model
order reduction techniques enable us to identify the coherent spatial–temporal
patterns.

The presence or absence of a neural mechanism can be analyzed for neuronal
populations. Dimensionality reduction methods [6] are data-driven statistical
techniques for forming and evaluating hypotheses about population activity
structure, which are summarized in Section 2. One of the goals of
neurosciences is fast and accurate predictions of the potential propagation in
neurons. The differential equations describing the propagation of potential in
neurons were developed and are described by Hodgkin and Huxley equations [12].
They consists of a coupled system of ordinary and partial differential
equations (ODEs and PDEs). The dimension of the associated discretized systems
is very large for accurately simulating neurons with realistic morphological
structure and synaptic inputs. In Section 3 we present two model order
reduction approaches based on POD and DEIM [5] which can predict accurately
the potential propagation in large scale neuronal networks leading to
important speedups [17, 16, 2]. Using the functional neuroimagining data from
electroencephalography (EEG) or functional magnetic resonance imaging (fMRI),
different regions of the brain can be inferred by dynamic causal modeling
(DCM) [7]. Effective connectivity is parameterised in terms of coupling among
unobserved brain states and neuronal activity in different regions. In Section
4 we describe a combined state and parameter reduction for parameter
estimation and identification [10] to extract effective connectivity in
neuronal networks from measured data, such as data from electroencephalography
(EEG) or functional magnetic resonance imaging (fMRI). In Section 5 the data-
driven, equation free, model order reduction method dynamic mode decomposition
(DMD) is described for identifying sleep spindle networks [3]. Reduced order
models with POD and DEIM and four variants of them are presented for neuronal
synaptic plasticity and neuronal spiking networks in Section 6.

##  2 Dimensionality reduction methods

Coordination of responses across neurons exists only at the level of the
population and not at the level of single neurons. The presence or absence of
a neural mechanism can be analyzed for neuronal populations. Dimensionality
reduction methods are data-driven statistical techniques for forming and
evaluating hypotheses about population activity structure. They produce low-
dimensional representations of high-dimensional data with the aim to extract
coherent patterns which preserve or highlight some feature of interest in the
data [6]. The recorded neurons of dimension $D$ are likely not independent of
each other, because they belong to a common network of neuronal populations.
From the high-dimensional data of neuronal recordings, a smaller number of
explanatory variables $K$ ( $K<D$) are extracted with the help of
dimensionality reduction methods. The explanatory variables are not directly
observed, therefore they are referred to as latent variables. The latent
variables define a $K$-dimensional space representing coherent patterns of the
noisy neural activity of $D$ neurons.

There exists several dimensionality reduction methods which differ in the
statistical interpretation of the preserved and discarded features of the
neuronal populations. We summarize the commonly used statistical methods of
dimensionality reduction methods in [6], where further references about the
methods can be found.

Principal component and factor analysis; The most widely known technique to
extract coherent patterns from high dimensional data is the modal
decomposition. A particularly popular modal decomposition technique is
principal component analysis (PCA), which derives modes ordered by their
ability to account for energy or variance in the data. In particular, PCA is a
static technique and does not model temporal dynamics of time-series data
explicitly, so it often performs poorly in reproducing dynamic data, such as
recordings of neural activity. The low-dimensional space identified by PCA
captures variance of all types, including firing rate variability and spiking
variability, whereas factor analysis (FA) discards the independent variance
for each neuron. and preserves variance that is shared across neurons.

Time series methods: The temporal dynamics of the population activity can be
identified if the data comes from a time series. The most commonly used time
series methods for dimensionality reduction neural recordings are: hidden
Markov models (HMM) [18], kernel smoothing followed by a static dimensionality
reduction method, Gaussian process factor analysis (GPFA) [35], latent linear
dynamical systems (LDS) [4] and latent nonlinear dynamical systems (NLDS)
[26]. They produce latent neural trajectories that capture the shared
variability across neurons. The HMM is applied when a jump between discrete
states of neurons exists, other methods identify smooth changes in firing
rates over time.

Methods with dependent variables: On many neuronal recordings the high-
dimensional firing rate space is associated with labels of one or more
dependent variables, like stimulus identity, decision identity or a time
index. The dimensionality reduction aims in this case to project the data such
that differences in these dependent variables are preserved. The linear
discriminant analysis (LDA) can be used to find a low-dimensional projection
in which the $G$ groups to which the data points belong are well separated.

Nonlinear dimensionality reduction methods: All the previous methods assume a
linear relationship between the latent and observed variables. When the data
lies on a low-dimensional, nonlinear manifold in the high-dimensional space, a
linear method may require more latent variables than the number of true
dimensions of the data. The most frequently used methods to identify nonlinear
manifolds are Isomap79 [31] and locally linear embedding (LLE) [28]. Because
the nonlinear methods use local neighborhoods to estimate the structure of the
manifold, population responses may not evenly explore the high-dimensional
space. Therefore theses methods should be used with care.

##  3 Proper orthogonal decomposition (POD) and discrete empirical
interpolation (DEIM) for Hodgin-Huxley model

One of the goals of neurosciences is fast and accurate predictions of the
potential propagation in neurons. The differential equations describing
propagation of potential in neurons are described by Hodgkin and Huxley (HH)
cable equations [12]. They consists of a coupled system of ordinary (ODEs) and
partial differential equations PDEs. Accurate simulation of morphology,
kinetics and synaptic inputs of neurons requires solution of large systems of
nonlinear ODEs. The complexity of the models are determined by the synapse
density of the dentritic length ($1\mu$ one micron). In simulations, for one
synapse per micron on a cell $5$ mm dendrite at $5,000$ compartments each with
$10$ variables are needed, which results in $50,000$ coupled nonlinear system
of ODEs [17, 16]. To recover complex dynamics, efficient reduced order
neuronal methods are developed using the POD and DEIM from the snapshots of
the in space and time discretized coupled PDEs and ODEs [17, 16, 2]. In this
section we describe two of them. They differ in the formulation of the HH
cable equation and of the equations for the gating variables.

###  3.1 Morphologically accurate reduced order modeling

The neuronal full order models (FOMs) in [17, 16] consists of $D$ branched
dendritic neurons $B=\sum_{d=1}^{D}B_{d}$ meeting at the soma, where the
$d^{th}$ has $B_{d}$ branches. It is assumed that the branch $b$ carries $C$
distinct ionic currents with associated densities and $G_{bc}(x)$ and reversal
potentials $E_{c},c=1,\ldots,C$. The kinetics of current $c$ on branch $b$ is
governed by the $F_{c}$ gating variables, $w_{bcf},f=1,\ldots,F_{c}$. When
subjected to input at $S_{b}$ synapses, the nonlinear HH cable equation for
the transmembrane potential $v_{b}(x,t)$ with the equation for the gating
variables $w_{bcf}$ is given by (see [2] Fig.1, model network with three
cables)

$\displaystyle a_{b}C_{m}\frac{\partial v_{b}}{\partial t}=$
$\displaystyle\frac{1}{2R_{i}}\frac{\partial}{x}\left(a_{b}^{2}\frac{\partial
v_{b}}{\partial x}\right)$ (1) $\displaystyle-
a_{b}\sum_{c=1}^{C}G_{bc}(x)(v_{b}-E_{c})\prod_{f=1}^{F_{c}}w_{bcf}^{q_{cf}}$
$\displaystyle\frac{1}{2\pi}\sum_{s=1}^{S_{b}}g_{bs}(t)\delta(x-x_{bs})(v_{b}-E_{bs})$
$\displaystyle\frac{\partial w_{bcf}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{w_{cf,\infty}(v_{b})-w_{bcf}}{\tau_{cf}(v_{b})},\quad
0<x<l_{b},\;t>0,$ (2)

where $g_{bs}(nS)$ is the time course, $x_{bs}$ is the spatial location and
$E_{bs}$ is the reversal potential of the $s$th synapse on branch $b$. The
variables and parameters in (1) are described in [17, 16].

These branch potentials interact at $J$ junction points, where junction $J$
denotes the soma. The $D$ dendrites join at soma. Continuity of the potential
at the soma leads to a common value at current balance denoted by
$v_{\sigma}(t)$. Then the networked form of (1) becomes

$\displaystyle a_{b}C_{m}\frac{\partial v_{\sigma}}{\partial t}=$
$\displaystyle\frac{\pi}{A_{\sigma}R_{i}}\sum_{d=1}^{D}\frac{\partial}{\partial
x}\left(a_{b_{J}^{d}}^{2}(l_{b_{J}^{d}})\frac{\partial
v_{b_{J^{d}}}(l_{b_{J^{d}}},t)}{\partial x}\right)$ (3) $\displaystyle-
a_{b}\sum_{c=1}^{C}G_{\sigma
c}(x)(v_{\sigma}-E_{c})\prod_{f=1}^{F_{c}}w_{\sigma cf}^{q_{cf}}(t)$
$\displaystyle\frac{1}{A_{\sigma}}\sum_{s=1}^{S_{b}}g_{\sigma
s}(t)(v_{\sigma}(t)-E_{\sigma s})$ $\displaystyle\frac{\partial w_{\sigma
cf}(t)}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{w_{cf,\infty}(v_{\sigma}(t))-w_{\sigma
cf}(t)}{\tau_{cf}(v_{\sigma})(t)},\quad 0<x<l_{b},\;t>0.$ (4)

When the cell is partitioned into $N$ compartments, with $C$ distinct ionic
currents per compartment and with $F$ gating variables per current, the
following nonlinear ODEs are obtained

$\displaystyle v^{\prime}(t)=$ $\displaystyle
Hv(t)-(\Phi(w(t)e).v(t)+\Phi(w(t))E_{i}$ (5)
$\displaystyle+G(t).(v(t)-E_{s}),\quad v(t)\in\mathbb{R}^{N}$ $\displaystyle
w^{\prime}(t)=$ $\displaystyle(A(v(t))-w(t))./B(v(t)),\quad
w(t)\in\mathbb{R}^{N\times C\times F}$ (6)

where $H\in\mathbb{R}^{N\times N}$ is the Hines matrix [11], $e=[1\;1\cdots
1]^{T}\in\mathbb{R}^{C}$ and the ‘dot’ operator, $a.b$, denotes element-wise
multiplication. $E_{i}$ and $E_{s}$ are respectively the vector of channel
reversal potentials and the vector of synaptic reversal potentials,
respectively Eq. (5) is discretized in time by the second order discretized
Euler scheme [11].

In [16] using the snapshots of $v(t)$ and of the nonlinear term
$N(v(t),w(t))\equiv(\Phi)w(t))e).v(t)-\Phi(w(t)))E_{i}$ at times
$t_{1},t_{2},\ldots,t_{n}$ the POD and DEIM modes are constructed.

The reduced membrane potential $v_{r}$ are constructed using the POD basis,
the reduced gating variables $w_{r}$ are obtained after applying the DEIM to
the nonlinear terms. The reduced order model in [16] preserves the spatial
precision of synaptic input, captures accurately the subthreshold and spiking
behaviors.

In [17] a linearized quasi active reduced neuronal model is constructed using
balanced truncation and ${\mathcal{H}}_{2}$ approximation of transfer
functions in time. ROMs preserve the input-output relationship and reproduce
only subthreshold dynamics.

###  3.2 Energy stable neuronal reduced order modeling

In [1, 2] a different form of the HH cable equation and ODEs for gating
variables is considered. The intracellular potential $v(x,t)$ and three gating
variables $m(x,t),\;h(x,t)$, and $n(x,t)$ describe the activation and
decativation of the ion channels, of the sodium channels and of the potassium
channels, respectively. A single cable in the computational domain
$(x,t)\in[0,L]\times(0,T]$ describing the distribution of the potential
$u(x,t)$ is given by [1, 2]

$\frac{\partial u}{\partial
t}=\frac{\mu}{a(x)}\left(a(x)^{2}u_{x}\right)_{x}-\frac{1}{C_{m}}g(m,h,n)u+\frac{1}{C_{m}}f(m,h,n,x,t),$
(7)

where $a(x)$ the radius of the neurons and $C_{m}$ is specific membrane
capacitance, $\mu=\frac{1}{2C_{m}R_{i}}>0$ the ratio with $R_{i}$ the axial
resistivity. The conductance $g(x,t)$ is a polynomial of the gating variables

$g(x,t)=g_{1}m^{3}h+g_{2}n^{4}+g_{3}>0,$ (8)

with the source term

$f(m,h,n,x,t)=g_{1}E_{1}m^{2}h+g_{2}E_{2}n^{4}+g_{3}E_{3}-i(x,t),$ (9)

where $E_{l},\;l=1,2,3$ are equilibrium potentials and $i(x,t)$ input current
at $x$

$i(x,t)=\sum_{s=1}^{N_{s}}i_{s}(x,t),\quad x\in[0,L].$ (10)

The nonlinear ODEs for the gating variables are given by

$\displaystyle\frac{\partial m}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{m}(v(x,t)(1-m(x,t))-\beta_{m}v(x,t))m(x,t),$
$\displaystyle\frac{\partial h}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{h}(v(x,t))(1-h(x,t))-\beta_{h}v(x,t))h(x,t),$ (11)
$\displaystyle\frac{\partial n}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{n}(v(x,t))(1-n(x,t))-\beta_{n}v(x,t))n(x,t),$

Expression for
$\alpha_{m},\;\alpha_{h},\;\alpha_{n},\;\beta_{m},\;\beta_{h},\;\beta_{n}$ and
boundary conditions can be found in [2].

In [1, 2], a model network with three cables connected to a soma is used. The
equations governing the potential propagation in the network $N_{c}$ neuron
cables-dentrites and /or axons with the superscript ${}^{(c)},\;c=1,\ldots
N_{c}$, are given as

$\displaystyle\frac{\partial v^{(c)}}{\partial t}=$
$\displaystyle\frac{\mu}{a^{(c)}(x^{(c)})}\left(\left(a^{(c)}\left(x^{(c)}\right)^{2}\right)v^{(c)}_{x}\right)_{x}-\frac{1}{C_{m}}g\left(m^{(c)},h^{(c)},n^{(c)}\right)u^{(c)}$
$\displaystyle+$
$\displaystyle\frac{1}{C_{m}}f\left(m^{(c)},h^{(c)},n^{(c)},x^{(c)},t\right)$
(12)

$\displaystyle\frac{\partial m^{(c)}}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{m}(v^{(c)}(1-m^{(c)})-\beta_{m}v^{(c)})m^{(c)},$
$\displaystyle\frac{\partial h^{(c)}}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{h}(v^{(c)})(1-h^{(c)})-\beta_{h}v{(c)})h^{(c)},$ (13)
$\displaystyle\frac{\partial n^{(c)}}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha_{n}(v^{(c)}))(1-n^{(c)})-\beta_{n}v^{(c)})n^{(c)},$

for $x^{(c)}\in\Omega^{(c)}=[0,L^{(c)}]$ together with the boundary
conditions.

The semi-discrete form of these equations are approximated using energy stable
summation by parts [1, 2] for the model network. The reduced order bases (ROB)
for multiple cables of identical lengths are assembled into a network in block
form. The block structure of the ROB allows a flexible structure-preserving
model reduction approach with an independent approximation in each cable and
energy stability and accuracy properties follow from this block structure.
Computation of the time varying reduced variables in the gating equations at
every time $t$ is costly because they scale with dimension of FOM. A
nonnegative variant of the discrete empirical interpolation method (NNDEIM) is
developed in [2] to preserve the structure and energy stability properties of
the equations.

The capability of the greedy-based approach to generate accurate predictions
in large-scale neuronal networks is demonstrated for system with more than
$15,000$ degree of freedoms (dofs). The state variable ROB of dimension $l=15$
with POD modes together with the nonnegative ROBs of dimension $p=60$ with
NNDEIM modes are constructed using a greedy approach to predict the potential
variation at the soma. The speedup of simulations is about $20$ larger than
Galerkin projection only is $1.3$ without using the NDEIM.

##  4 Combined state and parameter reduction for dynamic causal modelling

In neuroscience different regions of the brain are inferred using
neuroimagining data from EEG or fMRI recordings using the method od dynamic
causal modeling (DCM) [7]. Effective connectivity is parameterised in terms of
coupling among unobserved brain states and neuronal activity in different
regions. In DCM the neuronal activity is of the observed brain region is
represented as a SISO (single input single output) linear state-space system

$\dot{x}=A_{\mathrm{d}yn}(\mu)x+B_{\mathrm{d}yn}u$ (14)

with the parameterized connectivity $A_{\mathrm{d}yn}(\mu)$ and external input
matrices $B_{\mathrm{d}yn}$.

Linearization of the nonlinear DCM hemodynamic forward sub-model (balloon
model) [7] transforms the neuronal activity to the measured BOLD (blood oxygen
level dependent) response. Linearization around the equilibrium results in the
following single input, single output (SISO) system:

$\displaystyle B_{obs}$ $\displaystyle:=$ $\displaystyle(1\;0\;0\;0)^{T},\quad
C_{obs}=(0\;0\;V_{0}k_{1}\;V_{0}k_{2}),$ (15) $\displaystyle\dot{z}_{i}$
$\displaystyle=$ $\displaystyle A_{obs}z_{i}+B_{obs}x_{i},$ (16)
$\displaystyle y_{i}$ $\displaystyle=$ $\displaystyle C_{obs}z_{i},$ (17)
$\displaystyle z_{0}$ $\displaystyle=$ $\displaystyle(0\;0\;0\;0)^{T},$ (18)

$A_{\mathrm{o}bs}:=\left(\begin{array}[]{cccc}\frac{1}{\tau_{s}}&\frac{1}{\tau_{f}}&0&0\\\
1&0&0&0\\\
0&\frac{1}{\tau_{0}E_{0}}(1-(1-E_{0})(1-\ln(1-E_{0})))&\frac{1}{\tau_{0}}&\frac{1-\alpha}{\tau_{0}\alpha}\\\
0&\frac{1}{\tau_{0}}&0&\frac{1}{\tau_{0}\alpha}\end{array}\right).$ (19)

The fMRI measurements at the $i^{th}$ brain region are reflected by the output
variables $y_{i}$. For the meaning of the variables and parameters in (15) and
(19) we refer to [10, 9]. The linearized forward sub-models are embedded into
the fMRI connectivity model

$\left(\begin{array}[]{c}\dot{x}\\\ \dot{z}_{1}\\\ \dot{z}_{2}\\\ \vdots\\\
z_{N_{dyn}}\end{array}\right)=\left(\begin{array}[]{ccccc}A_{dyn}(\mu)&0&0&\cdots&0\\\
\delta_{1,1}&A_{obs}&0&&\\\ \delta_{2,1}&0&A_{obs}&&\\\ \vdots&&\ddots&\\\
\delta_{1,N_{dyn}}&&&A_{obs}\end{array}\right)\left(\begin{array}[]{c}x\\\
z_{1}\\\ z_{2}\\\ \vdots\\\
z_{N_{dyn}}\end{array}\right)+\left(\begin{array}[]{c}B_{dyn}\\\ 0\\\ 0\\\
\vdots\\\ 0\end{array}\right)v,$ (20)

$y=\left(0\left(\begin{array}[]{ccc}C_{obs}&&\\\ &\ddots&\\\
&&C_{obs}\end{array}\right)\right)\left(\begin{array}[]{c}x\\\ z_{1}\\\
z_{2}\\\ \vdots\\\ z_{N_{dyn}}\end{array}\right),$ (21)

where $\delta_{ij}\in\mathbb{R}^{4\times N_{\mathrm{d}yn}}$ denotes the
Kronecker matrix with non-zero elements located at the $(i,j)^{th}$ component.

The linearized state-space forward model (20) and (21) corresponds to a
multiple input, multiple output (MIMO) system

$\dot{x}(t)=A(\mu)x(t)+Bu(t),\qquad y(t)=Cx(t),$ (22)

where $x\in\mathbb{R}^{N}$ is the internal state, $u\in\mathbb{R}^{J}$ the
external input, $y\in\mathbb{R}^{O}$ the observed output and $\mu$ are the
parameters describing different conditions.

For large number of $M:=N^{2}$ parameters, the computational cost for
inferring the parameters and states is very high. In [10, 8] a combined state
and parameter model order reduction is developed for parameter estimation and
identification to extract effective connectivity. The inversion procedure
consists of two phases, the offline and online phases. In the offline phase,
the underlying parameterized model is reduced jointly in states and
parameters. In online phase, the reduced order model’s parameters are
estimated to fit the observed experimental data. Using state and parameter
reduction, the computational cost is reduced in the offline phase. The
simultaneous reduction of state and parameter space is based on Galerkin
projections with the orthogonal matrices for the state
$V\in\mathbb{R}^{N\times n}$ and for the parameters $P\in\mathbb{R}^{M\times
m}$. The reduced model is of lower order $n<<N,\;m<<M$ than the original full
order model. The reduced states $x_{r}(t)\in\mathbb{R}^{n}$ and the reduced
parameters $\mu\in\mathbb{R}^{m}$ are computed as

$\dot{x}_{r}(t)=A_{r}(\mu_{r})x_{r}(t)+B_{r}u(t),\qquad y_{r}(t)=C_{r}x(t)$
(23)

with a reduced initial condition $x_{r,0}=V^{T}x_{0}$ and the reduced
components

$\displaystyle\mu_{r}$ $\displaystyle=$ $\displaystyle
P^{T}\mu\in\mathbb{R}^{m},$ $\displaystyle A_{r}(\mu_{r})$ $\displaystyle=$
$\displaystyle V^{T}A(P\mu_{r})V\in\mathbb{R}^{n\times n},$ $\displaystyle
B_{r}$ $\displaystyle=$ $\displaystyle V^{T}B\in\mathbb{R}^{n\times J},$
$\displaystyle C_{r}$ $\displaystyle=$ $\displaystyle CV\in\mathbb{R}^{O\times
m}.$

In the online phase, the optimization based inverse problem is combined with
the reduction of state and parameter space. The inversion is based on
generalized data-driven optimization approach to construct the ROMs in [23]
and enhanced with the Monte-Carlo method to speed up the computations. The
state projection $V\in\mathbb{R}^{N\times n}$ and parameter projection
$P\in\mathbb{R}^{m\times m}$ are determined iteratively based on a greedy
algorithm by maximizing the error between the high-fidelity original and the
low-dimensional reduced model in the Bayesian setting.

Numerical experiments with the DCM model in [23] show highly dimensional
neuronal network system is efficiently inverted due to the short offline
durations. In the offline phase, Monte-Carlo enhanced methods require more
than an order of magnitude less offline time compared to the original and
data-misfit enhanced methods. In the online phase the reduced order model has
a speedup factor about an order of magnitude more compared to the full-order
inversion. The output error of the data-misfit enhanced method is close to
full order method. The output-errors of Monte-Carlo decrease with increasing
number of simulation but does not reach the output error of the full order
system. The source code is available in MATLAB [8].

##  5 Dynamic mode decomposition

Dynamic mode decomposition (DMD) is a data-driven, equation free ROM technique
[20]. It was initially developed to reduce the high dimensional dynamic data
obtained from experiments and simulations in the fluid mechanics into a small
number of coupled spatial–temporal modes [29, 30]. DMD was applied to explore
spatial–temporal patterns in large-scale neuronal recordings in [3]. DMD can
be interpreted as combination of discrete Fourier transform (DFT) in time and
principal component analysis (PCA) [15] in space. Both PCA and independent
component analyses (ICA) [13] are static techniques, which perform poorly in
reproducing dynamic data, such as recordings of neural activity.

The data is taken from electrocorticography (ECoG) recordings. Voltages from
$n$ channels of an electrode array sampled every $\Delta t$. These
measurements are arranged at snapshot $k$ to the column vectors
${\mathbf{x}}_{k}$. The $m$ snapshots in time construct to data matrices
shifted in time with $\Delta t$

${\mathbf{X}}=\left[\begin{array}[]{cccc}|&|&&|\\\
{\mathbf{x}}_{1}&{\mathbf{x}}_{2}&\cdots&{\mathbf{x}}_{m-1}\\\
|&|&&|\end{array}\right],\quad{\mathbf{X}}^{\prime}=\left[\begin{array}[]{cccc}|&|&&|\\\
{\mathbf{x}}_{2}&{\mathbf{x}}_{3}&\cdots&{\mathbf{x}}_{m}\\\
|&|&&|\end{array}\right]$ (24)

These matrices are assumed to be related linearly in time

${\mathbf{X}}^{\prime}={\mathbf{A}}{\mathbf{X}}.$ (25)

The DMD of the data matrix pair ${\mathbf{X}}$ and ${\mathbf{X}}^{\prime}$ is
given by the eigendecomposition of ${\mathbf{A}}$ using the singular value
decomposition (SVD) of the data matrix ${\mathbf{X}}=U\Sigma V^{*}$ by
computing the pseudoinverse
${\mathbf{A}}\approx{\mathbf{X}}^{\prime}{\mathbf{X}}^{\dagger}\equiv{\mathbf{X}}^{\prime}{\mathbf{V}}{\mathbf{\Sigma}}^{-1}{\mathbf{U}}^{*}.$
The spatio-temporal modes are computed by the exact DMD algorithm [32].

Because DMD does not contain explicit spatial relationship between neighboring
measurements, traveling waves occurring in neuronal networks can not be
captured well with a few coherent modes. DMD was also used as a windowed
technique with a temporal window size constrained by lower bound as for the
discrete Fourier transformation (DFT). In contrast to the fluid dynamics where
$n>>m$, in neuroscience the electrode arrays that have tens of channels $n$ in
the recordings with $m$ number of snapshots in the windows data per second, so
that $n<m$. The number of singular values $v$ of ${\mathbf{X}}$ are less than
$n$ and $m-1$, which restricts the maximum number of DMD modes and eigenvalues
to $n$. Because of this the dynamics can be captured over $m$ snapshots. The
rank mismatch is resolved by appending to the snapshot measurements with $h-1$
time-shifted versions of the data matrices. The augmented data matrix
${\mathbf{X}}_{\mathrm{a}ug}$ is given as

${\mathbf{X}}_{\mathrm{a}ug}=\left[\begin{array}[]{cccc}|&|&&|\\\
{\mathbf{x}}_{1}&{\mathbf{x}}_{2}&\cdots&{\mathbf{x}}_{m-h}\\\ |&|&&|\\\
|&|&&|\\\ {\mathbf{x}}_{2}&{\mathbf{x}}_{3}&\cdots&{\mathbf{x}}_{m-h-1}\\\
|&|&&|\\\ &&\cdots&\\\ |&|&&|\\\
{\mathbf{x}}_{h}&{\mathbf{x}}_{h+1}&\cdots&{\mathbf{x}}_{m-1}\\\ |&|&&|\\\
\end{array}\right].$ (26)

The augmented matrices ${\mathbf{X}}_{{\mathrm{a}ug}}$ and
${\mathbf{X}}^{\prime}_{{\mathrm{a}ug}}$ are Hankel matrices, which are
constant along the skew diagonal, as in the Eigenvalue Realization Algorithm
(ERA) [14]. The number of the stacks $h$ is chosen such that $hn>2m$. A
measure to determined the optimal number of stacks $h$ is the approximation
error

$E=\frac{||{\mathbf{X}}-\hat{\mathbf{X}}||_{F}}{||{\mathbf{X}}||_{F}}$

where $||\cdot||_{F}$ is the Frobenius norm. The approximation error $E$ is
decreasing with increasing number of stacks $h$ and reaches a plateau, so that
the DMD accuracy does not significantly increases.
DMD is applied in [3] as an automated approach to detect and analyze reliably
spatial localization and frequencies of sleep spindle networks from human ECoG
recordings. A MATLAB implementation is available at github.com/bwbrunton/dmd-
neuro/.

##  6 Reduced order modeling of biophysical neuronal networks

Recently reduced order models for ODEs

$\dot{x}(t)=A(t)x(t)+f(x(t))+Bu(t)$ (27)

are constructed using POD and DEIM to investigate input-output behavior of the
neuronal networks in the brain [22, 21], where $x(t)$ are state, and $u(t)$
are input variables, respectively.

The model in [22] is based on the chemical reactions of molecules in synapses,
that are the intercellular information transfer points of neurons. The
signaling pathways in striatal synaptic plasticity are modeled in [19]. This
model describes how certain molecules, which are a prerequisite for learning
in the brain, act in synapses. The stoichiometric equations obey the law of
mass action, which leads to a deterministic system of $44$ ODEs of the form
(27) . The state $x(t)$ of the control system (27) is a collection of ions,
molecules, and proteins that act in neuronal synapses. The linear part of (27)
is sparse, the nonlinearities are quadratic. The time dependent stimulus
$u(t)$ consists of molecules that are important for neuronal excitability and
plasticity, calcium and glutamate.

In [21], a nonlinear biophysical network model is considered, describing
synchronized population bursting behavior of heterogeneous pyramidal neurons
in the brain [27]. Neurons communicate by changing their membrane voltage to
create action potentials (spikes), propagating from cell to cell. Spiking is
the fundamental method of sensory information processing in the brain, and
synchronized spiking is an emergent property of biological neuronal networks.
The ODE system (27) in [21] consists of the states $x(t)$ as a collection of
$50$ neurons, each modeled with $10$ ODEs, leading totally to a system of ODEs
of dimension $500$. Each cell is modeled with Hodgkin-Huxley equations, where
each cell has only two compartments (dendrites and soma) and these
compartments have different ion channels. The state variables $x(t)$ include
the voltages of somatic and dendritic compartments, dendritic calcium
concentration, synaptic and ion channel gating variables and the input $u(t)$
is an injected current. Additionally, the soma compartment voltages are
coupled to dentritic compartments of randomly chosen cells. This networking of
the output of cells as input to other cells is key for producing synchronized
population behavior. The nonlinearities are Hodgkin-Huxley type, i.e.
exponential functions as well as cubic and quartic polynomials.

In [22], POD+DEIM was applied to a data-driven biological model of plasticity
in the brain (27). The ROMs with POD-DEIM reduce significantly the simulation
time and error between the original model and reduced order solutions can be
tuned by adjusting the number POD and DEIM bases independently. When the ROMs
are trained in a matching time interval of $10000$ seconds, accurate results
are obtained. However, generalizing the reduced model to longer time intervals
is challenging, which is characteristic for all nonlinear models.

In [21], the network model (27) is reduced with localized DEIM (LDEIM) [24],
discrete adaptive POD (DAPOD) [33, 34], and adaptive DEIM [25]. DEIM and the
variations are used here in combination with POD. ROMs require large number
POD and DEIM bases, to accurately simulate the input-output behavior in the
ROMs. In this model, every cell is heterogeneous in parameters and there are
also jump/reset conditions, which are factors that pose additional challenges
to the reduced order methods. However, the ROMs in [21] were able to replicate
the emergent synchronized population activity in the original network model.
DAPOD and ADEIM perform best in preserving the spiking activity of the
original network model. ADEIM is too slow and does not allow low enough
dimensions to offset the computational costs of online adaptivity. DAPOD is
able to find a lower dimensional POD basis online than the other methods find
offline, but the runtime is close to the original model.

## References

  * [1] D. Amsallem and J. Nordström.  High-order accurate difference schemes for the Hodgkin–Huxley equations.  Journal of Computational Physics, 252:573 – 590, 2013.
  * [2] D. Amsallem and J. Nordström.  Energy stable model reduction of neurons by nonnegative discrete empirical interpolation.  SIAM Journal on Scientific Computing, 38(2):B297–B326, 2016.
  * [3] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz.  Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition.  Journal of Neuroscience Methods, 258:1 – 15, 2016.
  * [4] L. Buesing, J. H. Macke, and M. Sahani.  Spectral learning of linear dynamics from generalised-linear observations with application to neural population data.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 25, pages 1682–1690. Curran Associates, Inc., 2012.
  * [5] S. Chaturantabut and D. Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
  * [6] John P Cunningham and Byron M Yu.  Dimensionality reduction for large-scale neural recordings.  Nature Neuroscience, 17(11):1500–1509, 2014.
  * [7] K.J. Friston, L. Harrison, and W. Penny.  Dynamic causal modelling.  NeuroImage, 19(4):1273 – 1302, 2003.
  * [8] C. Himpe.  optmor - optimization-based model order reduction (version 1.2), 2015\.
  * [9] C. Himpe.  Combined State and Parameter Reduction: For Nonlinear Systems with an Application in Neuroscience.  Internationaler Fachverlag für Wissenschaft & Praxis, 2016.
  * [10] C. Himpe and M. Ohlberger.  Data-driven combined state and parameter reduction for inverse problems.  Advances in Computational Mathematics, 41(5):1343–1364, 2015.
  * [11] M. Hines.  Efficient computation of branched nerve equations.  Int J Biomed Comput., 15(1):69–76, 1984.
  * [12] A. L. Hodgkin and A. F. Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  Bulletin of Mathematical Biology, 52(1):25–71, Jan 1990.
  * [13] A. Hyvärinen and E. Oja.  Independent component analysis: algorithms and applications.  Neural Networks, 13(4):411 – 430, 2000.
  * [14] Juang J.N. and Pappa R.S.  An eigensystem realization algorithm for modal parameter identification and model reduction.  Journal of Guidance, Control, and Dynamics, 8(5):620–627, 1985\.
  * [15] I. T. Jolliffe.  Principal component analysis.  Springer Series in Statistics. Springer-Verlag, New York, 2005.
  * [16] A. R. Kellems, S. Chaturantabut, D. C. Sorensen, and S. J. Cox.  Morphologically accurate reduced order modeling of spiking neurons.  Journal of Computational Neuroscience, 28(3):477–494, 2010.
  * [17] A. R. Kellems, D. Roos, N. Xiao, and S. J. Cox.  Low-dimensional, morphologically accurate models of subthreshold membrane potential.  Journal of Computational Neuroscience, 27(2):161, 2009.
  * [18] C. Kemere, G. Santhanam, B. M. Yu, A. Afshar, S. I. Ryu, T. H. Meng, and K. V. Shenoy.  Detecting neural-state transitions using hidden Markov models for motor cortical prostheses.  Journal of Neurophysiology, 100(4):2441–2452, 2008.
  * [19] B. Kim, SL. Hawes, F. Gillani, LJ. Wallace, and KT. Blackwell.  Signaling pathways involved in striatal synaptic plasticity are sensitive to temporal pattern and exhibit spatial specificity.  PLoS Comput Biol, 9(3):e1002953, 2013.
  * [20] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor.  Dynamic mode decomposition: Data-driven modeling of complex systems.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016.
  * [21] M. Lehtimäki, L. Paunonen, and M.-L. Linne.  Projection-based order reduction of a nonlinear biophysical neuronal network model.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 1–6, 2019.
  * [22] M. Lehtimäki, L. Paunonen, S. Pohjolainen, and M.-L. Linne.  Order reduction for a signaling pathway model of neuronal synaptic plasticity.  IFAC-PapersOnLine, 50(1):7687 – 7692, 2017.  20th IFAC World Congress.
  * [23] C. Lieberman, K. Willcox, and O. Ghattas.  Parameter and state model reduction for large-scale statistical inverse problems.  SIAM Journal on Scientific Computing, 32(5):2523–2542, 2010.
  * [24] B. Peherstorfer, D. Butnaru, K. Willcox, and H.-J. Bungartz.  Localized discrete empirical interpolation method.  SIAM J. Sci. Comput., 36(1):A168–A192, 2014.
  * [25] B. Peherstorfer and K. Willcox.  Online adaptive model reduction for nonlinear systems via low-rank updates.  SIAM J. Sci. Comput., 37(4):A2123–A2150, 2015.
  * [26] B. Petreska, B. M Yu, J. P Cunningham, S. Gopal, Stephen I. Y., K. V. Shenoy, and M. Sahani.  Dynamical segmentation of single trials from population neural data.  In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 756–764. Curran Associates, Inc., 2011.
  * [27] P. F. Pinsky and J. Rinzel.  Intrinsic and network rhythmogenesis in a reduced traub model for CA3 neurons.  Journal of Computational Neuroscience, 1(1):39–60, 1994.
  * [28] S. T. Roweis and L. K. Saul.  Nonlinear dimensionality reduction by locally linear embedding.  Science, 290(5500):2323–2326, 2000.
  * [29] C. W. Rowley, I. Mezić, S. Bahheri, P. Schlatter, and D. S. Hennigson.  Spectral analysis of nonlinear flows.  Journal of Fluid Mechanics, 641:115–127, 2009.
  * [30] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, 2010.
  * [31] J. B. Tenenbaum, V. de Silva, and J. C. Langford.  A global geometric framework for nonlinear dimensionality reduction.  Science, 290(5500):2319–2323, 2000.
  * [32] J. H. Tu, C. W. Rowley., D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz.  On dynamic mode decomposition: Theory and applications.  Journal of Computational Dynamics, 1(2158-2491_2014_2_391):391, 2014.
  * [33] M. Yang and A. Armaou.  Dissipative distributed parameter systems on-line reduction and control using DEIM/APOD combination.  In 2018 Annual American Control Conference (ACC), pages 2557–2562, 2018.
  * [34] M. Yang and A. Armaou.  Revisiting APOD accuracy for nonlinear control of transport reaction processes: A spatially discrete approach.  Chemical Engineering Science, 181:146 – 158, 2018.
  * [35] B. M. Yu, J. P. Cunningham, G. Santhanam, S. I. Ryu, K. V. Shenoy, and M. Sahani.  Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity.  Journal of Neurophysiology, 102(1):614–635, 2009.


# RIS-Aided D2D Communications Relying on Statistical CSI with Imperfect
Hardware

Zhangjie Peng, Tianshu Li, Cunhua Pan, Hong Ren,
and Jiangzhou Wang (_Corresponding author: Cunhua Pan_.)Z. Peng is with the
College of Information, Mechanical, and Electrical Engineering, Shanghai
Normal University, Shanghai 200234, China, also with the National Mobile
Communications Research Laboratory, Southeast University, Nanjing 210096,
China, and also with the Shanghai Engineering Research Center of Intelligent
Education and Bigdata, Shanghai Normal University, Shanghai 200234, China
$(\text{e-mail: pengzhangjie@shnu.edu.cn})$. T. Li is with the College of
Information, Mechanical and Electrical Engineering, Shanghai Normal
University, Shanghai 200234, China $(\text{e-mail:
1000479056@smail.shnu.edu.cn})$. C. Pan is with the School of Electronic
Engineering and Computer Science at Queen Mary University of London, London E1
4NS, U.K. $(\text{e-mail: {c.pan}@qmul.ac.uk})$.H. Ren is with the National
Mobile Communications Research Laboratory, Southeast University, Nanjing
210096, China $(\text{e-mail: {hren}@seu.edu.cn})$.J. Wang is with the School
of Engineering and Digital Arts, University of Kent, CT2 7NT Canterbury, U.K.
(e-mail: j.z.wang@kent.ac.uk).

###### Abstract

In this letter, we investigate a reconfigurable intelligent surfaces
(RIS)-aided device to device (D2D) communication system over Rician fading
channels with imperfect hardware including both hardware impairment at the
transceivers and phase noise at the RISs. This paper has optimized the phase
shift by a genetic algorithm (GA) method to maximize the achievable rate for
the continuous phase shifts (CPSs) and discrete phase shifts (DPSs). We also
consider the two special cases of no RIS hardware impairments (N-RIS-HWIs) and
no transceiver hardware impairments (N-T-HWIs). We present closed-form
expressions for the achievable rate of different cases and study the impact of
hardware impairments on the communication quality. Finally, simulation results
validate the analytic work.

###### Index Terms:

Reconfigurable intelligent surface (RIS), hardware impairment, D2D
communication, intelligent reflecting surface (IRS).

$\displaystyle{\rm{\gamma}}_{i}=\\!\frac{p_{i}\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{ai}\right|^{2}}{\sum\limits_{j=1,j\not=i}^{K}{\\!\left(p_{j}\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{aj}\right|^{2}\\!\right)}\\!+\\!|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\\!+\\!\varUpsilon_{ri}\\!+\\!\sigma_{i}^{2}}.$
(10) $\displaystyle
R_{i}\approx\\!\log_{2}\\!\left(\\!1\\!+\frac{p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\tilde{\Gamma}_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}}{\left(1+\kappa_{r}\right)\left(1+\kappa_{t}\right)\sum\limits_{j=1}^{K}{\left(p_{j}\alpha_{bi}\alpha_{aj}\frac{\varepsilon_{i}\beta_{j}\tilde{\Gamma}_{i,j}+L\left(\varepsilon_{i}+\beta_{j}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{j}+1\right)}\right)}-p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\tilde{\Gamma}_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}+\sigma_{i}^{2}}\right)$
(13)
$\vspace{-0.2cm}R_{i}^{\text{N\\!-RIS\\!-HWIs}}\approx\\!\log_{2}\\!\left(\\!1\\!+\frac{p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\Gamma_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}}{\left(1+\kappa_{r}\right)\left(1+\kappa_{t}\right)\sum\limits_{j=1}^{K}{\left(p_{j}\alpha_{bi}\alpha_{aj}\frac{\varepsilon_{i}\beta_{j}\Gamma_{i,j}+L\left(\varepsilon_{i}+\beta_{j}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{j}+1\right)}\right)}-p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\Gamma_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}+\sigma_{i}^{2}}\right)$
(15)

##  I Introduction

Reconfigurable intelligent surface (RIS) is a new transmission technology that
can configure the radio channel in a desired manner by optimizing the phase
shift of each reflecting element [1]. Due to their appealing properties of low
power consumption and low cost, RIS-aided wireless communications have
attracted much research attentions [2, 3, 4, 5]. Some initial attempts to
study RIS-aided communication systems include RIS-aided full-duplex systems
[2], RIS-aided physical layer security [3], RIS-aided wireless power transfer
[4], RIS-aided mobile edge computing [5], and RIS-aided multiuser transmission
[6]. The cascaded channel estimation was studied in [7].

D2D communication has received great attention to meet the rapidly increasing
demand for data traffic [8]. In D2D communication in cellular networks,
cellular and D2D users transmit signals simultaneously using the same spectrum
of cellular users. However, there is a paucity of contributions on RIS-aided
D2D communications. In [9], the authors proposed to deploy an RIS that can
assist pairs of devices in their communication when the direct links are
blocked by high buildings, plants and walls.

However, most of the existing contributions on D2D communication are based on
the assumption of perfect hardware in the radio frequency chains [10] and the
RISs [9]. Existing studies have shown that hardware impairments adversely
affect the system performance of the multiple-input multiple-output systems
[11, 12]. In addition, an RIS has low hardware cost [2, 3, 4, 5, 9], which is
prone to hardware imperfections.

Figure 1: System model for RIS-aided D2D communications.

Based on above, a natural question is whether we can use non-ideal low-cost
hardware to assist the D2D communications. Specifically, our contributions are
summarized as follows: 1) We consider an RIS-aided D2D communication system
over Rician fading channels, considering hardware impairments both at the
terminals and at the RISs. We assume that the hardware impairment is coupling
with the transmission signal at the transceiver. Moreover, the random phase
error caused by the phase noise at the RIS follows the Von Mises distribution;
2) we present the closed-form expressions for the achievable rate in the
general case of both hardware impairments both at the terminals and at the
RISs, which is then maximized by a genetic algorithm (GA) method. To obtain
more design insights, we also study the two special cases of no RIS hardware
impairments (N-RIS-HWIs) and no transceiver hardware impairments (N-T-HWIs);
3) we provide simulation results to verify our derived expressions. In
addition, the exhaustive search method is used to show that our results can
achieve the globally optimal solution.

The rest of the letter is organized as follows. The model for RIS-aided D2D
communication system with hardware impairments is depicted in Section II. We
derive the sum achievable rate’s expression and maximize it in Section III.
Numerical results are presented in Section IV. And conclusions are drawn in
Section V.

##  II System Model

We consider an RIS-aided D2D communication system, where an RIS assists $K$
pairs of devices to exchange information, as shown in Fig. 1. It is assumed
that there are a total of $K$ D2D links in the network, where the $i$th
single-antenna transmitter (receiver) is denoted by $\text{DT}_{i}$
($\text{DR}_{i}$) for $i=1,\cdots,K$. In the overlay mode, the radio resources
occupied by D2D links are orthogonal to that of the cellular links, which
ensures no interference between D2D links and cellular links. The RIS includes
$L$ scattered-reflection elements. And the phase shift matrix is denoted by
$\bm{\Theta}=diag(e^{j\theta_{1}},\cdots,e^{j\theta_{\ell}},\cdots,e^{j\theta_{L}})$,
where $\theta_{\ell}$ is the phase shift of the $\ell$th scattered-reflection
element.

The channel between $\text{DT}_{i}$ and the RIS (the RIS and $\text{DR}_{i}$)
can be written as

$\vspace{-0.1cm}\textbf{g}_{ai}=\sqrt{\alpha_{ai}}\textbf{h}_{ai},$ (1)
$\vspace{-0.2cm}\textbf{g}_{bi}=\sqrt{\alpha_{bi}}\textbf{h}_{bi},$ (2)

where $\alpha_{ai}$ and $\alpha_{bi}$ denote the large-scale fading
coefficients, and $\textbf{h}_{ai}\in{{\mathbb{C}}^{L\times 1}}$ and
$\textbf{h}_{bi}\in{{\mathbb{C}}^{L\times 1}}$ denote Rician fading channels,
which can be expressed as

$\vspace{-0.1cm}\textbf{h}_{ai}=\sqrt{\frac{\varepsilon_{i}}{\varepsilon_{i}+1}}\overline{\textbf{h}}_{ai}+\sqrt{\frac{1}{\varepsilon_{i}+1}}\tilde{\textbf{h}}_{ai},$
(3)
$\vspace{-0.2cm}\textbf{h}_{bi}=\sqrt{\frac{\beta_{i}}{\beta_{i}+1}}\overline{\textbf{h}}_{bi}+\sqrt{\frac{1}{\beta_{i}+1}}\tilde{\textbf{h}}_{bi},$
(4)

where $\varepsilon_{i}$ and $\beta_{i}$ denote the Rician factors,
$\tilde{\textbf{h}}_{ai}\in{{\mathbb{C}}^{L\times 1}}$ and
$\tilde{\textbf{h}}_{bi}\in{{\mathbb{C}}^{L\times 1}}$ are non-line-of-sight
components, whose entries are standard Gaussian random variables with
independent and identical distribution, i.e., $\mathcal{CN}(0,1)$, and
$\overline{\textbf{h}}_{ai}\in{{\mathbb{C}}^{L\times 1}}$ and
$\overline{\textbf{h}}_{bi}\in{{\mathbb{C}}^{L\times 1}}$ are line-of-sight
(LoS) components. Particularly, $\overline{\textbf{h}}_{ai}$ and
$\overline{\textbf{h}}_{bi}$ can be written as

$\vspace{-0.3cm}\\!\\!\\!\overline{\mathbf{h}}_{ai}\left(\varphi_{i}^{a},\varphi_{i}^{e}\right)\\!\\!=\\!\\!\left[\\!\\!\begin{array}[]{c}\\!\\!\\!\\!\\!1,\cdots,e^{j2\pi\frac{d}{\lambda}\left(x\sin\varphi_{i}^{a}\sin\varphi_{i}^{e}+y\cos\varphi_{i}^{e}\right)},\cdots,\\\
e^{j2\pi\frac{d}{\lambda}\left(\left(\sqrt{L}-1\right)\sin\varphi_{i}^{a}\sin\varphi_{i}^{e}+\left(\sqrt{L}-1\right)\cos\varphi_{i}^{e}\right)}\\\
\end{array}\\!\\!\\!\\!\\!\right]^{T}\\!\\!\\!\\!\\!,$ (5)
$\vspace{-0.2cm}\overline{\mathbf{h}}_{bi}\left(\varsigma_{i}^{a},\varsigma_{i}^{e}\right)\\!\\!=\\!\\!\left[\begin{array}[]{c}\\!\\!\\!\\!\\!1,\cdots,e^{j2\pi\frac{d}{\lambda}\left(x\sin\varsigma_{i}^{a}\sin\varsigma_{i}^{e}+y\cos\varsigma_{i}^{e}\right)},\cdots,\\\
e^{j2\pi\frac{d}{\lambda}\left(\left(\sqrt{L}-1\right)\sin\varsigma_{i}^{a}\sin\varsigma_{i}^{e}+\left(\sqrt{L}-1\right)\cos\varsigma_{i}^{e}\right)}\\\
\end{array}\\!\\!\\!\\!\\!\right]^{T}\\!\\!\\!\\!\\!,$ (6)

where $0\leqslant x,y\leqslant\sqrt{L}-1$,
$\varphi_{i}^{a},\varphi_{i}^{e}(\varsigma_{i}^{a},\varsigma_{i}^{e})$
respectively denote the $i$th pair of devices’ azimuth and elevation angles of
arrival (angle of departure). We assume $d=\frac{\lambda}{2}$ in our letter
[1].

The signal transmitted from $\text{DT}_{i}$ is given by [11]

$\vspace{-0.2cm}x_{i}=\hat{x}_{i}+\eta_{ti},$ (7)

where $\hat{x}_{i}$ $\sim\mathcal{CN}(0,1)$ represents the data symbol,
$\eta_{ti}\sim\mathcal{C}\mathcal{N}\left(0,\kappa_{t}\mathbb{E}\left\\{|\hat{x}_{i}|^{2}\right\\}\right)$
denotes the impact of impairments in the transmitter hardware, and
$\kappa_{t}\in$ (0, 1) is the coefficient which characterizes the impairments
at the transmitter.

The RIS’s hardware impairments can take the form of phase noise because the
RIS is a passive device and highly accurate reflection adjustment is
impractical. The phase noise matrix is denoted by
$\mathbf{\Phi}=diag(e^{j\varDelta\theta_{1}},\cdots,e^{j\varDelta\theta_{\ell}},\cdots,e^{j\varDelta\theta_{L}})$,
for $\ell=1,2,\cdots,L$. $\varDelta\theta_{\ell}$ is random phase error that
follows the Von Mises distribution with zero mean and concentration parameter
$\kappa_{\varDelta\theta}$. Additionally, $\varDelta\theta_{\ell}$ has a
characteristic function
$\chi=\frac{I_{1}(\kappa_{\varDelta\theta})}{I_{0}(\kappa_{\varDelta\theta})}$
with $I_{p}(\kappa_{\varDelta\theta})$ being the modified Bessel function of
the first kind and order $p$ [12].

We exploit statistical channel state information (CSI), which varies much
slowly than the instantaneous CSI and can be readily obtained. The signal
received at $\text{DR}_{i}$ is given by

$\displaystyle y_{i}$
$\displaystyle=\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi}\sum_{j=1}^{K}{\sqrt{p_{j}}}\mathbf{g}_{aj}x_{j}+\eta_{ri}+n_{i}$
$\displaystyle=\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\mathbf{x}+\eta_{ri}+n_{i},$ (8)

where
$\mathbf{G}_{a}\triangleq\left[\mathbf{g}_{a1},\cdots,\mathbf{g}_{aK}\right]$,
$\mathbf{x}\triangleq\left[x_{1},\cdots,x_{K}\right]^{\text{T}}$, $p_{j}$
denotes the transmit power of $\text{U}_{Aj}$ with
$\mathbf{\Lambda}\triangleq\text{diag}\left(p_{1},\cdots,p_{K}\right)$,
$\eta_{ri}\sim\mathcal{C}\mathcal{N}\left(0,\varUpsilon_{ri}\right)$ denotes
the impact of impairments in the receiver hardware with
$\varUpsilon_{ri}=\kappa_{r}\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\mathbf{x}|^{2}\right\\}$, $\kappa_{r}\in$ (0, 1)
is the coefficient which characterizes the impairments at the receiver, and
$n_{i}$ $\sim\mathcal{CN}(0,\sigma_{i}^{2})$ is the additive white Gaussian
noise at $\text{DR}_{i}$, for $i=1,\cdots,K$.

To analyze the achievable rates, we consider the expansion of (II) as follows

$\displaystyle y_{i}$
$\displaystyle=\sqrt{p_{i}}\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{ai}\hat{x}_{i}+\sum_{j=1,j\not=i}^{K}{\sqrt{p_{j}}}\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{aj}\hat{x}_{j}$ $\displaystyle+\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}+\eta_{ri}+n_{i},$ (9)

where $\bm{\eta}_{t}\triangleq\left[\eta_{t1},\cdots,\eta_{tK}\right]^{T}$.

From (II), we obtain the signal-to-interference plus noise ratio at
$\text{DR}_{i}$ in (10) shown at the bottom of this page.

Therefore, the achievable rate for $\text{DR}_{i}$ can be expressed as

$\vspace{-0.2cm}R_{i}=\mathbb{E}\\{\text{log}_{2}(1+{\rm{\gamma}}_{i})\\},$
(11)

and the sum achievable rate can be written as

$\vspace{-0.1cm}C=\sum_{i=1}^{K}R_{i}.$ (12)

##  III Achievable Rate Analysis

In this section, we first derive an approximate expression for the sum
achievable rate in the following theorem. Then, we consider two special cases
of N-RIS-HWIs and N-T-HWIs. Finally, we maximize the achievable rate with GA
(genetic algorithm) method.

######  Theorem 1.

In a D2D communication system aided by RIS, the ergodic sum achievable rate of
$\text{DR}_{i}$ can be approximated as (13), where $\tilde{\Gamma}_{i,i}$ and
$\tilde{\Gamma}_{i,j}$ can be given by

$\displaystyle\tilde{\Gamma}_{i,h}\\!=$
$\displaystyle\mathbb{E}\left\\{\\!\left|\sum_{\ell=1}^{L}{e^{j\left[\left(\theta_{\ell}+\varDelta\theta_{\ell}\right)+\pi
T_{i,h}^{n,m}\right]}}\right|^{2}\\!\right\\}\\!\\!$ $\displaystyle=$
$\displaystyle L+2\chi^{2}\\!\\!\\!\\!\\!\\!\\!\\!\sum_{1\leq m<n\leq
L}{\\!}\\!\\!\\!\\!\\!\\!\\!\cos\left(\theta_{n}\\!-\\!\theta_{m}+\pi
T^{n,m}_{i,h}\right),h\in\left\\{i,j\right\\}$ (14)

where
$T^{n,m}_{i,h}=\left(x_{n}-x_{m}\right)p_{i,h}+\left(y_{n}-y_{m}\right)q_{i,h}$,
with $x_{z}=\lfloor\frac{\left(z-1\right)}{\sqrt{L}}\rfloor$,
$y_{z}=\left(z-1\right)mod\sqrt{L}$, $z\in\left\\{m,n\right\\}$,
$p_{i,h}=\sin\varphi_{i}^{a}\sin\varphi_{i}^{e}+\sin\varsigma_{h}^{a}\sin\varsigma_{h}^{e}$,
and $q_{i,h}=\cos\varphi_{i}^{e}+\cos\varsigma_{h}^{e}$.

Proof: See Appendix A.  $\Box$

######  Corollary 1.

Assuming that the RIS hardware is ideal and thus there is no phase noise,
i.e., $\varDelta\theta_{\ell}=0$, for $i=1,2,\cdots,N$. It follows
$\mathbf{\Phi}=\mathbf{I}_{L}$. The rate of $\text{DR}_{i}$ on N-RIS-HWIs can
be approximated as (15), where $\Gamma_{i,i}$ and $\Gamma_{i,j}$ can be
defined as

$\vspace{-0.1cm}\Gamma_{i,h}\\!=\\!L+2\sum_{\mathclap{1\leq m<n\leq
L}}\cos\left(\theta_{n}\\!-\theta_{m}+\pi
T^{n,m}_{i,h}\right),h\in\left\\{i,j\right\\}.$ (16)

Proof: When $\mathbf{\Phi}=\mathbf{I}_{L}$, we can use [9, Eq. (20) Eq. (21)],
and write
$\mathbb{E}\left\\{\left|\textbf{h}_{bi}^{T}\bm{\Theta}\textbf{h}_{aj}\right|^{2}\right\\}$
and
$\mathbb{E}\left\\{\left|\textbf{h}_{bi}^{T}\bm{\Theta}\textbf{h}_{ai}\right|^{2}\right\\}$
as

$\vspace{-.3cm}\mathbb{E}\left\\{\left|\textbf{h}_{bi}^{T}\bm{\Theta}\textbf{h}_{ah}\right|^{2}\right\\}=\frac{\varepsilon_{i}\beta_{h}\Gamma_{i,h}+L(\varepsilon_{i}+\beta_{h})+L}{(\varepsilon_{i}+1)(\beta_{h}+1)},h\in\left\\{i,j\right\\},$
(17)

where $\Gamma_{i,h}$ is defined in (16).

By substituting (17) into (30), we are able to obtain the final result in
(15). This completes the proof. $\Box$

We aim to solve the achievable rate maximization problem by optimizing the
phase shifts matrix $\bm{\Theta}$ in the special case with only one pair of
users, i.e., $K=1$. Without loss of generality, the user pair is referred to
as user pair $i$. We can rewrite the achievable rate as

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!R_{i}^{\text{N-RIS-
HWIs}}\approx$
$\vspace{-0.2cm}\\!\log_{2}\\!\\!\left(\\!\\!\\!1\\!+\\!\frac{\frac{p_{i}\alpha_{bi}\alpha_{ai}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\left(\varepsilon_{i}\beta_{i}\Gamma_{i,i}\\!+\\!L\left(\\!\varepsilon_{i}+\beta_{i}\right)\\!+\\!L\right)}{\\!\\!\\!\\!\left(\kappa_{t}\kappa_{r}\\!+\\!\kappa_{t}\\!+\\!\kappa_{r}\\!\right)\\!\\!\frac{p_{i}\alpha_{bi}\alpha_{ai}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\\!\\!\left(\\!\varepsilon_{i}\beta_{i}\Gamma_{i,i}\\!+\\!L\\!\left(\varepsilon_{i}\\!+\\!\beta_{i}\right)\\!+\\!L\\!\right)\\!+\\!\sigma_{i}^{2}}\\!\\!\right)\\!\\!.$
(18)

The rate depends on the phase shifts matrix $\bm{\Theta}$ only through the
intermediate variable $\Gamma_{i,i}$. Since
$\Gamma_{i,i}\\!=\left|\sum_{\ell=1}^{L}{e^{j\left[\theta_{\ell}+\pi
T_{i,h}^{n,m}\right]}}\right|^{2}\leqslant\left(\sum_{\ell=1}^{L}{|e^{j\left[\theta_{\ell}+\pi
T_{i,i}^{n,m}\right]}|}\right)^{2}=L^{2}$ and $\Gamma_{i,i}\gg 0$, the
optimization problem can be formulated as follows

$\displaystyle\max_{\Gamma_{i,i}}\quad$ $\displaystyle R_{i}^{\text{N-RIS-
HWIs}}$ (19a) s.t. $\displaystyle 0\leqslant\Gamma_{i,i}\leqslant L^{2}.$
(19b)

The expression for the first-order derivative of
$\text{SINR}_{i}\left(\Gamma_{i,i}\right)$ with respect to $\Gamma_{i,i}$ is

$\vspace{-.8cm}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{\partial\text{SINR}_{i}\left(\Gamma_{i,i}\right)}{\partial\Gamma_{i,i}}=$
$\displaystyle\frac{\frac{p_{i}\alpha_{bi}\alpha_{ai}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\varepsilon_{i}\beta_{i}\sigma_{i}^{2}}{\left[\\!\left(\\!\kappa_{t}\kappa_{r}\\!+\\!\kappa_{t}\\!+\\!\kappa_{r}\\!\right)\\!\frac{p_{i}\alpha_{bi}\alpha_{ai}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\left(\varepsilon_{i}\beta_{i}\Gamma_{i,i}\\!+\\!L\left(\varepsilon_{i}+\beta_{i}\right)\\!+\\!L\right)\\!+\\!\sigma_{i}^{2}\right]^{2}}$
$\displaystyle\geqslant 0.$ (20)

Thus, the optimal design for $\bm{\Theta}$ corresponds to setting
$\Gamma_{i,i}=L^{2}$, where the optimal phase shifts of all the RIS elements
are given by $\theta_{\ell}=-\pi T_{i,h}^{n,m}+C_{0},\forall\ell$, and $C_{0}$
is an arbitrary constant.

Considering this single-user system with optimal phase shift, we assume the
transmit power is scaled as $p_{i}=\frac{E_{u}}{L^{2}}$ and
$p_{i}=\frac{E_{u}}{L}$ with $L\rightarrow\infty$. If the transmit power is
scaled as $p_{i}=\frac{E_{u}}{L^{2}}$, the rate becomes

$\displaystyle R_{i}^{\text{N}-\text{RIS}-\text{HWIs}}\rightarrow$ (21)
$\displaystyle\log_{2}\\!\\!\left(1\\!+\\!\frac{\frac{E_{u}\alpha_{bi}\alpha_{ai}\varepsilon_{i}\beta_{i}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}}{\\!\\!\\!\\!\left(\kappa_{t}\kappa_{r}\\!+\\!\kappa_{t}\\!+\\!\kappa_{r}\\!\right)\frac{E_{u}\alpha_{bi}\alpha_{ai}\varepsilon_{i}\beta_{i}}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\,\,\\!\\!\\!+\\!\sigma_{i}^{2}}\\!\\!\right),\
\text{as}\ L\rightarrow\infty.$

If the transmit power is scaled as $p_{i}=\frac{E_{u}}{L}$, the rate becomes

$\vspace{-.2cm}R_{i}^{\text{N}-\text{RIS}-\text{HWIs}}\rightarrow\\!\log_{2}\\!\\!\left(1\\!+\\!\frac{1}{\kappa_{t}\kappa_{r}\\!+\\!\kappa_{t}\\!+\\!\kappa_{r}\,\,\\!\\!\,\,\\!\\!\\!}\right),\;\text{as}\;L\rightarrow\infty.$
(22)

In this case, the achievable rate only depends on the impairment coefficients
at the transceiver when $L$ grows into infinity.

######  Corollary 2.

Assuming that the transceiver hardware is ideal, i.e.,
$\kappa_{t}=\kappa_{r}=0$. The rate of $\text{DR}_{i}$ with N-T-HWIs can be
approximated as

$\displaystyle R_{i}^{\text{N-T-HWIs}}\approx$
$\displaystyle\log_{2}\\!\left(\\!1\\!+\\!\frac{p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\tilde{\Gamma}_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}}{\\!\\!\\!\\!\\!\\!\\!\sum\limits_{j=1,j\not=i}^{K}{\left(p_{j}\alpha_{bi}\alpha_{aj}\frac{\varepsilon_{i}\beta_{j}\tilde{\Gamma}_{i,j}+L\left(\varepsilon_{i}+\beta_{j}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{j}+1\right)}\right)}\\!+\\!\sigma_{i}^{2}}\\!\right)\\!.$
(23)

We design the phase shift with one pair of users, the achievable rate is

$\vspace{-.25cm}R_{i}^{\text{N-T-
HWIs}}\approx\log_{2}\\!\left(\\!1\\!+\\!\frac{p_{i}\alpha_{bi}\alpha_{ai}\frac{\varepsilon_{i}\beta_{i}\tilde{\Gamma}_{i,i}+L\left(\varepsilon_{i}+\beta_{i}\right)+L}{\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}}{\\!\\!\sigma_{i}^{2}}\\!\right),$
(24)

which is a monotonically increasing function of $\tilde{\Gamma}_{i,i}$. The
optimal design for $\bm{\Theta}$ corresponds to setting
$\tilde{\Gamma}_{i,i}=L^{2}$, where the optimal phase shifts of all the RIS
elements are given by $\theta_{\ell}=-\pi T_{i,h}^{n,m}+C_{1},\forall\ell$,
and $C_{1}$ is an arbitrary constant.

Considering this single-user system with optimal phase shift, we assume the
transmit power is scaled as $p_{i}=\frac{E_{u}}{L^{2}}$ with
$L\rightarrow\infty$. If the transmit power is scaled as
$p_{i}=\frac{E_{u}}{L^{2}}$, the rate becomes

$\vspace{-0.2cm}R_{i}^{\text{N}-\text{T}-\text{HWIs}}\\!\\!\\!\\!\rightarrow\\!\log_{2}\\!\\!\left(\\!\\!1\\!+\\!\frac{E_{u}\alpha_{bi}\alpha_{ai}\varepsilon_{i}\beta_{i}}{\\!\\!\sigma_{i}^{2}\left(\varepsilon_{i}+1\right)\left(\beta_{i}+1\right)}\\!\right)\\!\\!,\;\text{as}\;L\rightarrow\infty.$
(25)

Due to the first-order derivative of the sum rate with respect to the phase
shift $\bm{\Theta}$ is quite hard to obtain, we maximize the rate with the GA
method adopted in [9] by optimizing the phase shifts. The complexity of the
proposed GA algorithm is $N_{t}*n$, where $N_{t}$ is the population size, and
$n$ is the number of generations evaluated. We take into account both
continuous phase shifts (CPSs) and discrete phase shifts (DPSs). In practice,
only a limited number of phase shifts can be used. We assume that the phase
shift of the reflecting elements is quantized with $B$ bits when considering
DPS, and thus $2^{B}$ phase shifts can be chosen for each reflecting element.

##  IV  Numerical Results

We evaluate the impact of various parameters on the data rate performance. We
set SNR = $p_{i}$, $\varepsilon_{i}$ = 10 dB, $\kappa=\kappa_{t}=\kappa_{r}$,
$\kappa_{\varDelta\theta}$ = 4, $\sigma_{i}^{2}=1$, and
$\varphi_{i}^{a}=\varphi_{i}^{e}$ and $\varsigma_{i}^{a}=\varsigma_{i}^{e}$
are respectively set as $\varphi_{i}$ and $\varsigma_{i}$ in [9], for
$i=1,\cdots,K$. Other parameters are set the same as [9].

Figure 2: Sum achievable rate versus SNR with $L$ = 16, $K$ = 6,
$\kappa=0.05$, $B$ = 2.

In Fig. 2, we depict the rate in (12) versus the SNR obtained from the
approximate expression in (30) and Monte-Carlo simulations from (11), when $L$
= 16, $K$ = 6, $\kappa$ = 0.05, $B$ = 2. The Monte-Carlo simulation results
match the analytical expressions well, which verifies the data rate
performance. Furthermore, the rate with N-RIS-HWIs performs better than the
case with RIS-HWIs. As $\kappa_{\varDelta\theta}$ decreases, the rate
decreases. Additionally, when $\kappa_{\varDelta\theta}=0$, the Von Mises
distribution degenerates into the uniform distribution, i.e.,
$\varDelta\theta_{\ell}\sim\mathcal{U}[-\pi,\pi)$, for $\ell=1,2,\cdots,L$. In
this case, the random scheme curve has the same performance as that of our
derived results, which demonstrate that there is no need to optimize the phase
shift. The phenomenon is consistent with the result in [12].

Figure 3: Sum achievable rate versus SNR with $L$ = 9, $K$ = 2 and $B$ = 2.

Fig. 3 shows the sum rate versus the SNR when $L$ = 8, $K$ = 2, and $B$ = 2.
We can reduce the harmful impact of the T-HWIs by tuning the phase shifts of
the RIS. Compared with the random scheme, the optimal phase shift can achieve
higher rate. It is worth noting that the proposed GA method achieves similar
performance to the exhaustive search, which implies that our proposed
algorithm can achieve almost the globally optimal solution. Moveover, we
observe that different levels of hardware impairment obtain different rate:
the higher the level of the hardware impairment, the worse the performance of
the rate. The rate with the ideal transceiver hardware (N-T-HWIs, $\kappa$ =
0) performs the best, as expected.

Figure 4: Sum achievable rate versus Rician factor with $B$ = 2 and SNR = 20
dB.

Fig. 4 illustrates the rate versus the Rician factor when SNR = 20 dB and $B$
= 2. In all cases, with different values of $L$ and $K$, the rates approach
the fixed value as Rician factor $\varepsilon_{i}\rightarrow\infty$. This is
because the channels are mainly influenced by LoS component when Rician factor
is large.

Figure 5: Sum achievable rate versus $B$ with $\varepsilon_{i}$ = 10 dB and
SNR = 20 dB.

Fig. 5 shows the sum rate versus $B$. The figure shows that, under the
condition of DPS, three quantization bits are sufficient, which provides
useful insight into the engineering design of systems aided by RISs. In
addition, to resolve the issue caused by imperfect hardware, we can increase
the number of low-cost reflection elements.

##  V Conclusion

In this letter, we investigated an RIS-aided D2D communication system over
Rician fading channels, considering hardware impairments both at the terminals
and the RISs. We considered two special cases of N-RIS-HWIs and N-T-HWIs. We
derived closed-form expressions for the achievable rate of different cases.
The impact of the hardware impairment on the system performance have been
observed. To resolve the issue caused by imperfect hardware, we can increase
the number of low-cost reflection elements. Additionally, three quantization
bits are sufficient for the DPSs setup, which provides useful insight into the
engineering design of systems aided by RISs. Moreover, the extension to
jointly optimizing power allocation and the phase shifts of the RIS will be
left for our future work.

##  Appendix A Proof of Theorem 1

$\vspace{-0.4cm}R_{i}\approx\log_{2}\left(1+\frac{p_{i}\mathbb{E}\left\\{\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{ai}\right|^{2}\right\\}}{\sum\limits_{j=1,j\not=i}^{K}{\left(p_{j}\mathbb{E}\left\\{\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{aj}\right|^{2}\right\\}\right)}+\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}+\varUpsilon_{ri}+\sigma_{i}^{2}}\right)$
(23)
$\displaystyle\varUpsilon_{rj}=\kappa_{r}\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\mathbf{x}|^{2}\right\\}=\kappa_{r}\mathbb{E}\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\left(\mathbf{\hat{x}}+\bm{\eta}_{t}\right)\left(\mathbf{\hat{x}}^{H}+\bm{\eta}_{t}^{H}\right)\sqrt{\mathbf{\Lambda}}^{H}\mathbf{G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\right\\}$
$\displaystyle\\!\\!\\!\\!\\!\\!\overset{\text{(a)}}{=}\\!\kappa_{r}\mathbb{E}\\!\left\\{\\!\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\mathbb{E}\\!\left\\{\\!\mathbf{\hat{x}\hat{x}}^{H}\\!+\\!\mathbf{\hat{x}}\bm{\eta}_{t}^{H}\\!+\\!\bm{\eta}_{t}\mathbf{\hat{x}}^{H}\\!+\\!\bm{\eta}_{t}\bm{\eta}_{t}^{H}\\!\right\\}\\!\sqrt{\mathbf{\Lambda}}^{H}\mathbf{G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\\!\right\\}\\!\overset{\text{(b)}}{=}\\!\kappa_{r}\left(1\\!+\\!\kappa_{t}\right)\mathbb{E}\\!\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\mathbf{\Lambda
G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\\!\right\\}\\!$
(26)

$\vspace{-0.2cm}R_{i}\approx\\!\log_{2}\\!\left(\\!1\\!+\frac{p_{i}\alpha_{bi}\alpha_{ai}\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{ai}\right|^{2}\right\\}}{\sum\limits_{j=1,j\not=i}^{K}{\left(p_{j}\alpha_{bi}\alpha_{aj}\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}\right|^{2}\right\\}\right)}+\left(\kappa_{t}\kappa_{r}+\kappa_{t}+\kappa_{r}\right)\alpha_{bi}\sum\limits_{j=1}^{K}{\left(\alpha_{aj}p_{j}\mathbb{E}\left\\{|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}|^{2}\right\\}\right)}+\sigma_{i}^{2}}\right)$

$\vspace{-0.25cm}=\log_{2}\\!\left(\\!1\\!+\frac{p_{i}\alpha_{bi}\alpha_{ai}\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{ai}\right|^{2}\right\\}}{\left(1+\kappa_{r}\right)\left(1+\kappa_{t}\right)\sum\limits_{j=1}^{K}{\left(p_{j}\alpha_{bi}\alpha_{aj}\mathbb{E}\left\\{|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}|^{2}\right\\}\right)}-p_{i}\alpha_{bi}\alpha_{ai}\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{ai}\right|^{2}\right\\}+\sigma_{i}^{2}}\right)$ (30)

Using Lemma 1 in [13], $R_{i}$ in (11) is approximated as (23).

Next, we derive $\mathbb{E}\left\\{\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{ai}\right|^{2}\right\\}$,
$\mathbb{E}\left\\{\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{aj}\right|^{2}\right\\}$,
$\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}$, and
$\varUpsilon_{rj}$. To begin with, we have

$\vspace{-.2cm}\mathbb{E}\left\\{\left|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
g}_{ah}\right|^{2}\right\\}=\alpha_{bi}\alpha_{ah}\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{ah}\right|^{2}\right\\},h\in\left\\{i,j\right\\}.$ (24)

Both $\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}$ and
$\varUpsilon_{rj}$ contain the terms of
$\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{ai}\right|^{2}\right\\}$ and
$\mathbb{E}\left\\{\left|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}\right|^{2}\right\\}$. We derive these items later.

Then, $\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}$ can be written as
follows

$\vspace{-.2cm}\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left\\{\\!|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}\\!\\!=\\!\\!\kappa_{t}\mathbb{E}\\!\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\mathbf{\Lambda}\mathbf{G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\right\\}.\\!\\!\\!$
(25)

$\varUpsilon_{rj}$ also contains the terms of
$\mathbb{E}\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\mathbf{\Lambda}\mathbf{G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\right\\}$,
which we will derive later.

Moreover, $\varUpsilon_{rj}$ is calculated in (A) at the top of this page.
Equality in (a) is obtained because $\mathbf{\hat{x}}$ is independent of both
$\mathbf{g}_{bi}$ and $\mathbf{G}_{a}$; $\bm{\eta}_{t}$ is also independent of
both $\mathbf{g}_{bi}$ and $\mathbf{G}_{a}$. Additionally, equality in (b)
uses the following results

$\vspace{-0.2cm}\mathbb{E}\left\\{\mathbf{\hat{x}\hat{x}}^{H}\right\\}=\mathbf{I}_{K},\
\mathbb{E}\left\\{\bm{\eta}_{t}\bm{\eta}_{t}^{H}\right\\}=\kappa_{t}\mathbf{I}_{K},$
(27)
$\vspace{-0.15cm}\mathbb{E}\left\\{\mathbf{\hat{x}}\bm{\eta}_{t}^{H}\right\\}=\mathbf{0},\
\mathbb{E}\left\\{\bm{\eta}_{t}\mathbf{\hat{x}}^{H}\right\\}=\mathbf{0}.$ (28)

We can find $\mathbb{E}\left\\{|\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\sqrt{\mathbf{\Lambda}}\bm{\eta}_{t}|^{2}\right\\}$ and
$\varUpsilon_{rj}$ contain the item
$\mathbb{E}\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi G}_{a}\mathbf{\Lambda
G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\right\\}$,
given by

$\displaystyle\mathbb{E}\left\\{\mathbf{g}_{bi}^{T}\mathbf{\Theta\Phi
G}_{a}\mathbf{\Lambda
G}^{H}_{a}\mathbf{\Phi}^{H}\mathbf{\Theta}^{H}\mathbf{g}_{bi}^{*}\right\\}$
$\displaystyle\overset{\text{(c)}}{=}\alpha_{bi}\mathbb{E}\left\\{\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
H}_{a}\sqrt{\mathbf{D}_{a}}\mathbf{\Lambda}\sqrt{\mathbf{D}_{a}}^{H}\mathbf{H}_{a}^{H}\mathbf{\Theta\Phi}^{H}\mathbf{h}_{bi}^{*}\right\\}$
$\displaystyle=\alpha_{bi}\mathbb{E}\left\\{\sum_{j=1}^{K}{\alpha_{aj}p_{j}|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}|^{2}}\right\\}$
$\vspace{-.2cm}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\sum_{j=1}^{K}{\alpha_{bi}\alpha_{aj}p_{j}\mathbb{E}\left\\{|\mathbf{h}_{bi}^{T}\mathbf{\Theta\Phi
h}_{aj}|^{2}\right\\}}.$ (29)

We obtain equality in (c) by defining
$\mathbf{D}_{a}\triangleq\text{diag}\left(\alpha_{a1},\cdots,\alpha_{aK}\right)$.

Substituting (24), (25), (A) and (29) into (23), we obtain (30). Employing the
results of [9, Eq. (20) Eq. (21)] one obtains

$\vspace{-0.25cm}\mathbb{E}\left\\{\left|\textbf{h}_{bi}^{T}\bm{\Theta\Phi}\textbf{h}_{ah}\right|^{2}\right\\}=\frac{\varepsilon_{i}\beta_{h}\tilde{\Gamma}_{i,h}+L(\varepsilon_{i}+\beta_{h})+L}{(\varepsilon_{i}+1)(\beta_{h}+1)},h\in\left\\{i,j\right\\},$
(31)

where $\tilde{\Gamma}_{i,h}$ is defined by

$\vspace{-0.2cm}\tilde{\Gamma}_{i,h}\\!\\!=\\!\\!L\\!+\\!2\\!\\!\\!\\!\\!\\!\\!\\!\sum_{1\leq
m<n\leq
L}\\!\\!\\!\\!\\!\\!\\!{\mathbb{E}\left\\{\cos\left[\left(\theta_{n}+\varDelta\theta_{n}\right)\\!-\left(\\!\theta_{m}+\varDelta\theta_{m}\right)+\pi
T^{n,m}_{i,h}\right]\right\\}}.$ (32)

With the help of

$\mathbb{E}\left\\{\cos\left[\left(\theta_{n}+\varDelta\theta_{n}\right)\\!-\left(\\!\theta_{m}+\varDelta\theta_{m}\right)+\pi
T^{n,m}_{i,h}\right]\right\\}$
$\vspace{-.2cm}=\chi^{2}\cos\left[\left(\theta_{n}\\!-\theta_{m}+\pi
T^{n,m}_{i,h}\right)\right],$ (33)

$\tilde{\Gamma}_{i,h}$ can be further given in (1).

By substituting (31) into (30), we are able to obtain the final result in
(13). The proof of theorem 1 is completed.

## References

  * [1] C. Pan _et al._ , “Reconfigurable intelligent surfaces for 6G systems: Principles, applications, and research directions,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 14–20, Jun. 2021.
  * [2] Z. Peng, Z. Zhang, C. Pan, L. Li, and A. L. Swindlehurst, “Multiuser full-duplex two-way communications via intelligent reflecting surface,” _IEEE Trans. Signal Process._ , vol. 69, pp. 837–851, Jan. 2021.
  * [3] H. Shen, W. Xu, S. Gong, Z. He, and C. Zhao, “Secrecy rate maximization for intelligent reflecting surface assisted multi-antenna communications,” _IEEE Commun. Lett._ , vol. 23, no. 9, pp. 1488–1492, Sep. 2019.
  * [4] C. Pan _et al._ , “Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1719–1734, Aug. 2020.
  * [5] 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, no. 11, pp. 2666–2682, Nov. 2020.
  * [6] Z. Yang, W. Xu, C. Huang, J. Shi, and M. Shikh-Bahaei, “Beamforming design for multiuser transmission through reconfigurable intelligent surface,” _IEEE Trans. Commun._ , vol. 69, no. 1, pp. 589–601, Jan. 2021.
  * [7] B. Zheng and R. Zhang, “Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 4, pp. 518–522, Apr. 2020.
  * [8] J. Lee and J. H. Lee, “Performance analysis and resource allocation for cooperative D2D communication in cellular networks with multiple D2D pairs,” _IEEE Commun. Lett._ , vol. 23, no. 5, pp. 909–912, May 2019.
  * [9] Z. Peng, T. Li, C. Pan, H. Ren, W. Xu, and M. D. Renzo, “Analysis and optimization for RIS-aided multi-pair communications relying on statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 4, pp. 3897–3901, Apr. 2021.
  * [10] N. Serafimovski _et al._ , “Practical implementation of spatial modulation,” _IEEE Trans. Veh. Technol._ , vol. 62, no. 9, pp. 4511–4523, Nov. 2013.
  * [11] A. Afana and S. Ikki, “Analytical framework for space shift keying mimo systems with hardware impairments and co-channel interference,” _IEEE Commun. Lett._ , vol. 21, no. 3, pp. 488–491, Mar. 2017.
  * [12] A. Papazafeiropoulos, C. Pan, P. Kourtessis, S. Chatzinotas, and J. M. Senior, “Intelligent Reflecting Surface-assisted MU-MISO Systems with Imperfect Hardware: Channel Estimation, Beamforming Design,” 2021. [Online]. Available: https://arxiv.org/abs/2102.05333v1
  * [13] Q. Zhang, S. Jin, K.-K. Wong, H. Zhu, and M. Matthaiou, “Power scaling of uplink massive MIMO systems with arbitrary-rank channel means,” _IEEE J. Sel. Topics Signal Process._ , vol. 8, no. 5, pp. 966–981, Oct. 2014.