# PAPR Reduction for SVD/MIMO-OFDM via CS-PTS with Novel Low-Complexity Swarm
Intelligence

Luis F. Abanto-Leon

###### Abstract

Despite of the superior performance of cyclically shifted partial transmit
sequences (CS-PTS), this technique demands excessive complexity owing to the
exhaustive search. Customarily used in SISO systems solely, in this paper we
adapt CS-PTS to operate in MIMO-OFDM systems with singular value decomposition
(SVD) precoding. Moreover, to lessen the computational complexity due to the
exhaustive search, we introduce a novel search method based on swarm
intelligence, which can attain a fair balance between optimality and
complexity compared to other competing techniques. Our new approach—called
_best triad mating cuckoo search algorithm_ (BTM-CSA)—takes over the task of
searching sub-optimal solutions while maintaining the complexity affordable.
Computer simulations show that the proposed approach can attain near-optimal
performance at much lower complexity.

###### Index Terms:

MIMO-OFDM, SVD, PAPR reduction, cuckoo search algorithm, swarm intelligence.

##  I Introduction

Orthogonal frequency division multiplexing (OFDM) has been adopted by several
communication systems as pivotal technology due to its capability to provide
high data rate transmission, high spectral efficiency, robustness to multipath
fading [1] and narrow-band interference [2]. The inherent downside of OFDM is,
however, the large peak-to-average power ratio (PAPR) [3] provoked by
constructive combination of modulated sub-carriers which may occasionally
align in phase, thus generating large amplitudes. OFDM signals with large PAPR
experience inefficient power transmission [4] and distortion upon
amplification—caused by the nonlinearity of RF amplifiers. To accomplish
distortionless amplification, the input signal power may be backed-off so as
to force the amplifier to operate in the linear region. Nevertheless, this
imposition will cause increased power consumption and compromise energy
efficiency of the amplifier. In light of this fact, it is desirable to reckon
with signals that have low PAPR. Moreover, a low PAPR entails additional
merits such as ($i$) improvement of signal-to-quantization noise ratio (SQNR)
and ($ii$) simplification of analog-to-digital and digital-to-analog
converters [6].

In pursuance of reducing the PAPR of OFDM signals, several approaches have
been proposed. Clipping [4] consists on saturating the amplitude of a signal
wherein it exceeds a preset threshold. Nevertheless, it causes in-band
distortion and out-of-band radiation, which deteriorate the bit error rate
(BER) [5]. Companding [6] consists on signal compression at the transmitter—to
avoid nonlinear amplification and therefore distortion—and signal expanding at
the receiver, which aims at restoring the original amplitude of the signal.
However, the expanding process at the receiver also magnifies small-valued
noise added in the transmission channel causing BER degradation. Other
techniques such as tone reservation (TR) [7] and tone injection (TI) [7] are
capable of reducing the PAPR without incurring distortion. TR preselects a
subset of subcarriers that will not carry data information and exclusively
serve for canceling large peaks. In Tone Injection, the conventional
constellation is expanded such that one PSK/QAM symbol can be mapped into one
of several possible representations. Though effective, TI increases the
average transmit power and requires an exhaustive search process to find the
best mapping choice that minimizes the PAPR.

Another subgroup of distortionless techniques prevents the large peaks by
modifying the statistics of the OFDM signal through the application of
suitable phase rotations to each subcarrier. The most representative
techniques of this subgroup are Selective Mapping (SLM) [8] and Partial
Transmit Sequences (PTS) [9]. In SLM a set of several phase rotation patterns
that effect all the subcarriers are tested, and the one that achieves the
lowest PAPR is selected for signal transformation. The disadvantage of this
method is, nonetheless, the high complexity for computing several full-length
IFFT operations that derive from evaluating each pattern. On the other hand,
PTS divides the modulated subcarriers into disjoint compositions—called
partial transmit sequences—which are phase-rotated and linearly combined in
order to produce a signal with a PAPR lower than the original one. The
drawback, as well, is the excessive complexity caused by the exhaustive search
that is executed over all the possible combinations of phase rotations. Both
SLM and PTS need to transmit side information (SI), which represents
additional data needed to revert the phase rotations applied at the
transmitter.

Cyclically Shifted Partial Transmit Sequences (CS-PTS) [10] is a technique
that although allowing superior PAPR reduction than PTS and SLM has not
attracted much attention. In principle, CS-PTS is based on PTS but its
distinctive feature is the incorporation of time-domain shifting, which
provides an additional degree of freedom that enables linear phase rotation
per subcarrier and makes possible further PAPR reduction. On the downside, if
compared to PTS, CS-PTS exhibits substantial increase of search complexity
caused by the expanded dimensionality as a consequence of the additional
degree of freedom.

In order to address the high complexity search of CS-PTS, we resort to
employing a metaheuristic approach so as to control the number of evaluations
and therefore maintain the search complexity affordable without affecting the
performance severely. A metaheuristic approach called Cuckoo Search Algorithm
(CSA) was proposed in [11], wherein extensive simulations showed that CSA is
capable of outperforming other previous methods such as Particle Swarm
Optimization (PSO) [12] and Genetic Algorithms (GA) [13], in terms of
convergence rate. Metaheuristics are characterized by two main properties:
exploration and exploitation. The former is concerned with the capability of
effectively sampling/scanning the search space. As opposed to exhaustive
search where every possible solution is tested, the goal of metaheuristics is
to find an acceptable sub-optimal solution at lower computational expense
without inspecting every solution. The second property, i.e., exploitation,
denotes the capability of capitalizing on information in previous generations
to produce more fit solutions. CSA has remarkable exploration capability,
which may be attributed to Levy flights. Nevertheless, the exploitation
capability of CSA is limited. We have realized that if the exploitation
property is boosted, the convergence rate can be improved. In light of this
fact, we have coalesced additional mechanisms with CSA, purposing to ($i$)
boost the exploitation capability and ($ii$) achieve fairer balance between
exploration and exploitation, which will impact on the convergence rate. The
outcome of merging CSA with the proposed mechanisms results in Best Triad
Mating Cuckoo Search Algorithm (BTM-CSA), which will take over the task of
minimizing the PAPR.

Our contributions on this work can be summarized as follows:

  * •

We have adapted CS-PTS—standardly used in SISO systems only—to operate in
MIMO-OFDM systems with SVD precoding.

  * •

We propose a new metaheuristic approach called BTM-CSA, which will take over
the task of finding near-optimal solutions at enhanced convergence rate.

The remainder of this paper is organized as follows. In section II, we first
revisit the CS-PTS technique in SISO systems. We also describe the
modifications applied to CS-PTS in order to operate in SVD/MIMO-OFDM systems.
In section III, we present the BTM-CSA algorithm and describe in detail the
new modifications. In section IV, we explain how to use BTM-CSA to reduce the
PAPR in our system. In section V, the computational complexity of the proposed
approach as well as that of other techniques are shown. We devote section VI
to show and discuss the results of the proposed approach. Finally, in section
VII, we summarize our conclusions.

Figure 1: Transmitter of a SISO/CS-PTS system

##  II Cyclically Shifted Partial Transmit Sequences

###  II-A CS-PTS in SISO systems

The primal structure of the transmitter side of a SISO/CS-PTS system using
$N_{c}$ subcarriers is depicted in Fig. 1. The mapping of binary data into
PSK/QAM symbols produces a stream
$\mathbf{\bar{X}}=\big{[}\bar{X}[0],\bar{X}[1],\dots,\bar{X}[N_{c}-1]\big{]}^{\top}$
whose elements have a duration of $\Delta T_{s}$. Upon serial to parallel
conversion, the modulation symbols
$\mathbf{\breve{X}}=\big{[}\breve{X}[0],\breve{X}[1],\dotsc,\breve{X}[N_{c}-1]\big{]}^{\top}$
with enlarged duration $N_{c}\Delta T_{s}$ are obtained. The modulation
symbols $\breve{X}[k],k=0,1,\dotsc,N_{c}-1$ are divided into a set of $M$
disjoint partitions $\big{\\{}\mathbf{\breve{X}}_{m}\big{\\}}_{m=0}^{M-1}$
which are zero-padded by $L(N_{c}-1)$ zeros as shown in (1)

$\mathbf{\breve{X}}_{m}=\big{[}\underbrace{0\dotsc 0}_{\textrm{$k_{1}$
zeros}},\underbrace{\breve{X}[k_{1}],\dots,\breve{X}[k_{1}+\nicefrac{{N_{c}}}{{M}}-1]}_{\textrm{$\frac{N_{c}}{M}$
symbols}},\underbrace{0\dotsc 0}_{\textrm{$k_{2}$ zeros}}\big{]}^{\top}$ (1)

where $k_{1}=m\frac{N_{c}}{M}$ and $k_{2}=N_{c}L-(m+1)\frac{N_{c}}{M}$. The
parameter $L\geq 4$ represents the oversampling factor, whose purpose is
improving the time-bin resolution of the IFFT output so as to obtain a more
accurate representation of the continuous-time domain signal. This enables to
work with a good estimate of the true peak amplitude and therefore permits a
fairer calculation of the PAPR. Each partition $\mathbf{\breve{X}}_{m}$ is
inputted to the IFFT, which performs the modulation of the subcarriers by the
symbols assigned to that partition. At the output of each IFFT, a time-domain
partial transmit sequence $\mathbf{\breve{x}}_{m}$ is obtained via
$\mathbf{\breve{x}}_{m}=\text{IFFT}\big{\\{}\mathbf{\breve{X}}_{m}\big{\\}}$

$\Big{[}\big{\\{}\tilde{\gamma}_{m}\big{\\}}_{m=0}^{M-1},\big{\\{}\tilde{\delta}_{m}\big{\\}}_{m=0}^{M-1}\Big{]}=\underset{\big{[}\gamma_{0},\dots,\gamma_{M-1},\delta_{0},\dots,\delta_{M-1}\big{]}}{\arg\min}\left(\frac{\underset{0\leq
k\leq
N_{c}L}{\max}\Bigg{|}\displaystyle\sum_{m=0}^{M-1}\gamma_{m}\breve{x}_{m}\Bigg{[}k+\delta_{m}-N_{c}L\bigg{\lfloor}\frac{k+\delta_{m}}{N_{c}L}\bigg{\rfloor}\Bigg{]}\Bigg{|}^{2}}{\displaystyle\frac{1}{N_{c}L}\displaystyle\sum_{k=0}^{N_{c}L-1}\Bigg{|}\displaystyle\sum_{m=0}^{M-1}\gamma_{m}\breve{x}_{m}\Bigg{[}k+\delta_{m}-N_{c}L\bigg{\lfloor}\frac{k+\delta_{m}}{N_{c}L}\bigg{\rfloor}\Bigg{]}\Bigg{|}^{2}}\right)$
(6)

Each partial transmit sequence $\mathbf{\breve{x}}_{m}$ is simultaneously
weighted by a phase rotation $\gamma_{m}$ and cyclically delayed by a time
shift $\delta_{m}$, in order to reduce the amplitude correlation among
sequences and thus reduce the PAPR. As depicted in Fig. 1, each phase rotation
$\gamma_{m}$ is applied to every subcarrier of its corresponding $m$-th
partition. Similarly, each time shift $\delta_{m}$ effects a cyclic delay to
every subcarrier in the $m$-th partition. The selection of phase rotations
$\\{\gamma_{m}\\}_{m=0}^{M-1}$ is constrained to a set of elements so as to
avoid exponential growth of the search complexity. This set is defined by
$\Phi_{U}=\Set{e^{j\frac{2\pi}{U}u}\mid u=0,1,\dotsc,U-1}$, where $U$
represents the number of possible phase rotations. In a similar fashion, the
selection of time shifts is also restricted to a number of values in order to
maintain the search complexity within reasonable limits. The set of time
shifts is defined by (LABEL:e4) $\Phi_{D}=\Set{\frac{N_{c}}{D}d\mid
d=0,1,\dotsc,D-1}$, where $D$ represents the number of allowed time shifts.
From (4) we may realize that the time shifts are integer values equally
distributed over the range $\big{[}0,\frac{D-1}{D}N_{c}\big{]}$. Recall that
cyclic time-domain shifting produces linear variation of the phase response of
a signal. Therefore, as a result of cyclically delaying the partial transmit
sequences $\mathbf{\breve{x}}_{m}$, additional phase rotation per subcarrier
can be enabled. That is to say, every subcarrier $k$ belonging to the same
$m$-th partition will be rotated by an additional linearly increasing or
decreasing phase rotation $\theta_{m}^{(k)}$, defined by (LABEL:e5)
$\theta_{m}^{(k)}=\frac{2\pi k}{N_{c}}\delta_{m}$.

By coalescing the effect of both—phase rotations and time shifts—the overall
rotation in each subcarrier is equivalent to $\gamma_{m}$ \+
$\theta_{m}^{(k)}$. Time shifting adds a degree of freedom that is not
available in PTS, thus providing CS-PTS with more flexibility to accomplish
further PAPR reduction. Note that in order to reduce the PAPR of the transmit
signal, CS-PTS optimizes 2$M$ parameters: $M$ phase rotations and $M$ time
shifts. The process to find the optimal combination of $\gamma_{m}$ and
$\delta_{m}$ for the $M$ partitions is accomplished by means of an exhaustive
search over all the possible combinations. The PAPR minimization problem is
described by (6). Upon evaluation of every pattern, the best combination of
parameters, $\big{\\{}\tilde{\gamma}_{m}\big{\\}}_{m=0}^{M-1}$ and
$\big{\\{}\tilde{\delta}_{m}\big{\\}}_{m=0}^{M-1}$, capable of achieving the
lowest PAPR is selected for signal transformation and transmission. The
transformed OFDM signal
$\mathbf{x}=\big{[}x[0],x[1],\dots,x[N_{c}-1]\big{]}^{\top}$ ready for
transmission is expressed by (7)

$x[k]=\displaystyle\sum_{m=0}^{M-1}\tilde{\gamma}_{m}\breve{x}_{m}\Bigg{[}k+\tilde{\delta}_{m}-N_{c}L\bigg{\lfloor}\frac{k+\tilde{\delta}_{m}}{N_{c}L}\bigg{\rfloor}\Bigg{]}$
(7)

where $\lfloor.\rfloor$ is the floor function. Additionally, bear in mind that
the optimal parameters $\big{\\{}\tilde{\gamma}_{m}\big{\\}}_{m=0}^{M-1}$,
$\big{\\{}\tilde{\delta}_{m}\big{\\}}_{m=0}^{M-1}$ are to be sent as side
information to the receiver for enabling data recovery. At the receiver, phase
rotations are reversed by multiplying by their conjugate.

###  II-B CS-PTS in SVD/MIMO-OFDM systems

Fig. 2 depicts the transmitter side of the CS-PTS approach when Singular Value
Decomposition (SVD) precoding is employed with $N_{T}$ antennas at the
transmitter. Let
$\mathbf{\bar{X}}^{(i)}=\big{[}\bar{X}^{(i)}[0],\bar{X}^{(i)}[1],\dots,\bar{X}^{(i)}[N_{c}-1]\big{]}^{\top}$
(for $i=1,2,\dots,N_{T}$) denote the data symbols at the $i$-th antenna, where
each symbol has a duration of $\Delta T$ seconds. Upon serial to parallel
conversion, the symbols
$\mathbf{\hat{X}}^{(i)}=\big{[}\hat{X}^{(i)}[0],\hat{X}^{(i)}[1],\dots,\hat{X}^{(i)}[N_{c}-1]\big{]}^{\top}$
are obtained. SVD precoding requires channel state information, which is
assumed to be known at the transmitter and receiver.

Let $\mathbf{H}\in\mathbb{C}^{N_{R}\times N_{T}\times N_{c}}$ be the channel
matrix representing the state of all the links between $N_{R}$ antennas at the
receiver and $N_{T}$ antennas at the transmitter on all the $N_{c}$
subcarriers. Also let $\mathbf{H}_{k}\in\mathbb{C}^{N_{R}\times N_{T}}$
represent the channel state at the $k$-th subcarrier. By means of SVD, the
channel matrix can be decomposed into the product of three matrices, i.e.,
$\mathbf{H}_{k}=\mathbf{U}_{k}\mathbf{\Sigma}_{k}\mathbf{V}^{\dagger}_{k}$,
where $\dagger$ represents the Hermitian transpose. The complex unitary
matrices $\mathbf{U}_{k}$ and $\mathbf{V}_{k}$ have dimensions $N_{R}\times
N_{R}$ and $N_{T}\times N_{T}$, respectively, whereas $\mathbf{\Sigma}_{k}$ is
a rectangular diagonal matrix of dimensions $N_{R}\times N_{T}$ with non-
negative real entries. The precoding process consists in pre-multiplying the
input data in all the antennas by a matrix $\mathbf{V}_{k}$, in a subcarrier-
wise fashion, i.e.,
$\mathbf{\breve{X}}_{k}=\mathbf{V}_{k}\mathbf{\hat{X}}_{k}$, where
$\mathbf{\hat{X}}_{k}=\big{[}\hat{X}^{(1)}[k],\hat{X}^{(2)}[k],\dots,\hat{X}^{(N_{T})}[k]\big{]}^{\top}$
and
$\mathbf{\breve{X}}_{k}=\big{[}\breve{X}^{(1)}[k],\breve{X}^{(2)}[k],\dots,\breve{X}^{(N_{T})}[k]\big{]}^{\top}$.
After precoding has been accomplished in every subcarrier, the precoded
vectors
$\mathbf{\breve{X}}^{(i)}=\big{[}\breve{X}^{(i)}[0],\breve{X}^{(i)}[1],\dotsc,\breve{X}^{(i)}[N_{c}-1]\big{]}^{\top}$
at each antenna $i$ are split into $M$ non-overlapping sub-blocks, in a
similar fashion as has been detailed in (1), i.e.,
$\mathbf{\breve{X}}^{(i)}=\sum_{m=1}^{M-1}\mathbf{\breve{X}}^{(i)}_{m}$. After
IFFT, $M$ partial transmit sequences $\mathbf{\breve{x}}^{(i)}_{m}$ (for
$m=1,2,\dots,M-1$) per each transmit antenna are obtained. In order to
preserve the orthogonality of SVD, the same phase rotations and time shifts
must be applied to every $m$-th partition. Thus, the optimization problem is
reduced to finding the optimal phase rotations and time shifts that minimize
the maximum PAPR of all the antennas, i.e., the MIN-MAX value as shown in (8).

$\Big{[}\big{\\{}\tilde{\gamma}_{m}\big{\\}}_{m=0}^{M-1},\big{\\{}\tilde{\delta}_{m}\big{\\}}_{m=0}^{M-1}\Big{]}=\underset{\big{[}\gamma_{0},\dots,\gamma_{M-1},\delta_{0},\dots,\delta_{M-1}\big{]}}{\arg\min}\left\\{\underset{1\leq
i\leq N_{T}}{\max}\left(\frac{\underset{0\leq k\leq
N_{c}L}{\max}\Bigg{|}\displaystyle\sum_{m=0}^{M-1}\gamma_{m}\breve{x}^{(i)}_{m}\Bigg{[}k+\delta_{m}-N_{c}L\bigg{\lfloor}\frac{k+\delta_{m}}{N_{c}L}\bigg{\rfloor}\Bigg{]}\Bigg{|}^{2}}{\displaystyle\frac{1}{N_{c}L}\displaystyle\sum_{k=0}^{N_{c}L-1}\Bigg{|}\displaystyle\sum_{m=0}^{M-1}\gamma_{m}\breve{x}^{(i)}_{m}\Bigg{[}k+\delta_{m}-N_{c}L\bigg{\lfloor}\frac{k+\delta_{m}}{N_{c}L}\bigg{\rfloor}\Bigg{]}\Bigg{|}^{2}}\right)\right\\}$
(8)
Figure 2: Transmitter of a SVD/CS-PTS system

##  III Best Triad Mating Cuckoo Search Algorithm

Cuckoo Search Algorithm is an optimization technique inspired on the
reproduction strategy of some species of cuckoo birds that engage brood
parasitism in order to ensure the species survival. The strategy adopted by
parasite cuckoo birds is a type of kleptoparasitism, i.e., an individual gains
a benefit effortlessly, such as a prey or object that it could not have
obtained itself [14]. Cuckoo birds deceive other (host) birds in such a
fashion that the latter ones look after the cuckoo offspring unknowingly. This
tactic relieves cuckoo birds from laborious and lengthy tasks such as nest
construction and offspring feeding. Thus, cuckoo birds can devote more time
for food foraging and reproduction, which in overall will improve their
chances of survival [15]. These birds have evolved to such a sophisticated
level that they are able to breed eggs with similar patterns and colors to the
eggs of other bird species—intending to decrease the chances of the cuckoo
eggs being identified by the host birds. However, host birds are sometimes
able to identify the alien eggs and may decide either abandoning the nest or
ejecting the parasite eggs from the nest.

Cuckoo Search Algorithm captures the core of the reproduction strategy of
cuckoo birds. It is manifest to think that the reproductive behavior of cuckoo
birds has reached such a level of sophistication as a consequence of natural
evolution, and therefore it is efficient and effective in maximizing its goal.
For this reason, CSA is endowed with good search features that achieve
outstanding performance. For instance, CSA has affordable computational
complexity and remarkable exploration capability. However, we have identified
that the exploitation capability of CSA can be further developed. Exploration
and exploitation are two inherent properties of metaheuristics which must be
fairly balanced in order to attain fast convergence. Exploration is directly
related to the scanning or sampling efficiency of the solution search space.
Exploitation is concerned with the capability of availing information in
previous solutions to generate more fit solutions that can satisfy the problem
to a better extent. The advantage of metaheuristics is that in lieu of
obtaining the global optimum at the expense of high complexity—as occurs in
exhaustive search—it is possible to procure a near-optimal solution at lower
complexity provided that exploitation and exploration are well balanced. CSA
heavily depends on exploration since the candidate solutions are generated by
means of either Levy flights or using and approach similar to Differential
Evolution (DE) [16]. In order to overcome the strong dependence on exploration
and improve the balance between the two properties, four mechanisms have been
coalesced with CSA resulting in BTM-CSA. Although BTM-CSA reckons with
additional features, the average complexity of generating a candidate solution
is not increased because the proposed improvements are not applied to every
agent but only to a limited number with high potential to improve the
solutions. BTM-CSA exhibits fairer balance of properties than CSA, and this is
verifiable because of its faster convergence rate.

The pseudo-code of BTM-CSA is shown in Algorithm III. Preliminary
considerations to take into account for complete understanding of BTM-CSA are
summarized in the following statements.

  * •

There is an initial population of $N$ cuckoo birds and equal number of host
nests.

  * •

Each nest represents a potential solution to the objective function and the
suitability of each is defined by a fitness function $f$. Oftentimes, the
terms egg and nest are used interchangeably.

  * •

The nests that have the highest quality, i.e., highest fitness, will carry
over the next generation of birds.

  * •

Host birds may discover that their nest has been parasitized with a
probability $p_{a}$. In such an event, the nest will be abandoned and a new
one will be built in a different position on the space.

Algorithm 1 Best Triad Mating Cuckoo Search Algorithm

1:$N$ : number of host nests

2:$T$ : maximum number of evaluations

3:$p_{a}$ : fraction of the total nests that represent the worst solutions

4:$n$ : number of dimensions

5:$\tau$ : counter of evaluations

6:$N_{\text{G}}$ : number of nests in the neighborhood of the best nest

7:$N_{\text{H}}$ : number of nests to be mutated

8:begin

9: Define the objective function $f:\mathbb{R}^{n}\to\mathbb{R}$

10: Generate N host nests in $\Omega=\\{\mathbf{z}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N\\}$

11: Evaluate the fitness $F_{i}=\frac{1}{1+f(\mathbf{z}_{i})}$ of each nest
$\mathbf{z}_{i}$

12: while $(\sim stopcriterion)$ or $(\tau<T)$ do

13: Find the current fittest nest $\mathbf{z}_{i}$ in $\Omega$

14: Choose randomly another nest $\mathbf{z}_{i}$ from $\Omega$ avoiding
$\mathbf{z}_{\text{best}}$

15: Generate a new nest $\mathbf{z}_{j}$ via a Levy flight from
$\mathbf{z}_{i}$

16: Evaluate the fitness $F_{j}=\frac{1}{1+f(\mathbf{z}_{j})}$

17: if ($F_{j}>F_{i}$) then

18: Replace $\mathbf{z}_{i}$ by the new solution $\mathbf{z}_{j}$

19: end if

20: Generate a new nest $\mathbf{z}_{k}$ via a Levy flight from
$\mathbf{z}_{\text{best}}$

21: Evaluate the fitness $F_{k}=\frac{1}{1+f(\mathbf{z}_{k})}$

22: if ($F_{k}>F_{\text{best}}$) then

23: Replace $\mathbf{z}_{\text{best}}$ by the new solution $\mathbf{z}_{k}$

24: end if

25: Update the counter $\tau$

26: Let $\Omega_{\text{Q}}=\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$
be defined by the three nests with the highest fitness in $\Omega$, such that
$\displaystyle\hskip
28.45274pt\begin{cases}\Omega_{\text{Q}}\subset\Omega,\Omega_{\text{P}}\subset\Omega\\\
\Omega_{\text{Q}}\cup\Omega_{\text{P}}=\Omega\\\
\Omega_{\text{Q}}\cap\Omega_{\text{P}}=\emptyset\\\
\Omega_{\text{P}}=\\{\mathbf{p}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N-3\\}\end{cases}$

27: Let the hyperplane $\mathcal{H}$ be defined by
$\mathbf{q}_{1},\mathbf{q}_{2}$ and $\mathbf{q}_{3}$

28: Compute the parameters $\mathbf{\Gamma}$ and $\varepsilon$ of
$\mathcal{H}$

29: Compute the parameter $\mathbf{\Gamma}_{*}$

30: Define the set
$\Omega_{\text{R}}=\\{\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}\\}$ that
contains three potentially fitter nests obtained via Gaussian random walks
using $\mathbf{\Gamma}_{*}$ as a reference

31: Define
$\Omega_{\text{S}}=\Omega_{\text{Q}}\cup\Omega_{\text{R}}=\\{\mathbf{s}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,6\\}$ and sort the elements such that $F_{\mathbf{s}_{i}}\geq
F_{\mathbf{s}_{i+1}}$

32: Replace the elements of $\Omega_{\text{Q}}$ by the first 3 elements of
$\Omega_{\text{S}}$, i.e.,
$\Omega_{\text{Q}}=\\{\mathbf{s}_{i}\in\Omega_{\text{S}}\mid i=1,2,3\\}$

33: Update the set $\Omega$, such that
$\Omega=\Omega_{\text{P}}\cup\Omega_{\text{Q}}$

34: Update the counter $\tau$

35: Build the subset $\Omega_{\text{A}}=\\{\mathbf{a}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{A}}\\}$ consisting of potentially $N_{\text{A}}$ worse
nests in $\Omega$, such that  $\displaystyle\hskip
28.45274pt\begin{cases}\Omega_{\text{A}}\subset\Omega,\Omega_{\text{B}}\subset\Omega\\\
\Omega_{\text{A}}\cup\Omega_{\text{B}}=\Omega\\\
\Omega_{\text{A}}\cap\Omega_{\text{B}}=\emptyset\\\
\Omega_{\text{B}}=\\{\mathbf{b}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{B}}\\}\\\ N_{\text{A}}=p_{a}N\\\
N_{\text{B}}=N-N_{\text{A}}\end{cases}$

36: Build a subset $\Omega_{\text{C}}=\\{\mathbf{c}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{A}}\\}$ consisting of $N_{A}$ potentially fitter nests
than those contained in $\Omega_{\text{A}}$

37: Build the set $\Omega_{\text{D}}=\Omega_{\text{A}}\cup\Omega_{\text{C}}$
consisting of $2N_{\text{A}}$ elements

38: Sort the elements of
$\Omega_{\text{D}}=\\{\mathbf{d}_{i}\in\mathbb{R}^{n}\mid i=\hskip
24.18501pt1,2,\dots,2N_{\text{A}}\\}$ such that $F_{\mathbf{d}_{i}}\geq
F_{\mathbf{d}_{i+1}}$

39: Replace the elements of $\Omega_{\text{A}}$ by the first $N_{\text{A}}$
elements of $\Omega_{\text{D}}$, such that
$\Omega_{\text{A}}=\\{\mathbf{d}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{A}}\\}$

40: Update the set $\Omega$, such that
$\Omega_{\text{A}}\cup\Omega_{\text{B}}=\Omega$

41: Update the counter $\tau$

42: Find the current best nest $\mathbf{z}_{\text{best}}$

43: Build a set $\Omega_{\text{G}}=\\{\mathbf{g}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{G}}\\}$ composed of $N_{\text{G}}$ new nests generated
around the neighborhood of $\mathbf{z}_{\text{best}}$ using uniform random
walks

44: Select randomly $N_{\text{G}}$ nests from $\Omega$ avoiding
$\mathbf{z}_{\text{best}}$ and replace them by the newly generated eggs in
$\Omega_{\text{G}}$

45: Update the counter $\tau$

46: Build the set $\Omega_{\text{H}}=\\{\mathbf{h}_{i}\in\mathbb{R}^{n}\mid
i=1,2,\dots,N_{\text{H}}\\}$ containing the $\Omega_{\text{H}}$ least fit
nests and apply mutation

47: Keep the current best solutions for the next generation

48: Update the counter $\tau$

49: end while

50:end

In the following lines, we describe each of the modifications incorporated to
CSA. The first mechanism is explained in lines $20\sim 24$ of Algorithm III
and it deals with appropriate placement of the best nest in the search space.
As opposed to CSA, in BTM-CSA the best nest can be replaced by a more suitable
one if higher fitness is met. The second mechanism—detailed in the lines
$26\sim 34$—aims at fostering data sharing among the best solutions in order
to enhance the convergence rate by means of a concept that we call best triad
mating. Best triad mating exploits information on the top three solutions and
its objective is to intensify the search on a localized search space with high
probability of containing more fit solutions. This idea stems from the fact
that in CSA, cuckoo birds show limited social interaction. Thus, by promoting
collaborative sharing of relevant information, the convergence can be
improved. In the third mechanism—described in lines $42\sim 44$—we assume that
the best nest can accommodate more than one egg simultaneously, which will
intensify exploration in the vicinity of it in the following iterations.
Conversely, in CSA each nest accommodates a single egg only. Finally, the
purpose of the fourth mechanism—lines $46\sim 48$—is to create new nests in
different locations in order to replace a subset of nests that exhibit low
fitness by others with possibly higher expectancy of achieving superior
fitness. To wit, appropriate balance between exploitation and exploration is
highly desirable as this will regulate diversification of new solutions and
reusage of old ones.

In this work we use interchangeably the terms iterations and generations,
which must be understood as a complete run of the algorithm. CSA starts with
an initial set of $N$ nests or solutions, however, within a single generation
$2N$ solutions are generated from which only $N$ are preserved to the next
generation. BTM-CSA is as well initialized with $N$ eggs at the beginning and
within a generation $2N+N_{\nu}$ eggs are generated, although only $N$ are
preserved to the next generation as observed in Fig. 3. For the time being it
suffices to know that the initial $N$ solutions combine through different
procedures to temporarily obtain $2N+N_{\nu}$ eggs. We explain later in this
section why $2N+N_{\nu}$ solutions are generated. Note that within a
generation many solutions can be generated and each solution requires its
fitness to be evaluated. Therefore the number of generated solutions is equal
to the number of fitness evaluations, or simply evaluations. The counter
$\tau$ will increase by $2N$ after each generation—if CSA is used—or by
$2N+N_{\nu}$ in the case of BTM-CSA, until a maximum of $T$. Note that the
number of evaluations per generation is completely unrelated to the
computational complexity of an algorithm. It is more consistent to consider
the computational complexity per generated solution together with the total
number of evaluations in order to obtain a fair estimate of the overall
complexity.

………… $\Longleftarrow N$ initial solutions $\Longleftarrow$ $N+N_{\nu}$
discarded solutions $\Longleftarrow$ $N$ solutions preserved to the next
generation Input layer Processing layer Output layer Selection layer
Generation $t$ Figure 3: Generalized structure of BTM-CSA

###  III-A Mechanism 1: Update of the best solution

In CSA, the nest with the highest fitness in each generation is used as a
reference for generating Levy flights for the rest of nests. However, the best
nest is not allowed to update until the next iteration. And unless a solution
with higher fitness is generated, the best nest may remain unchanged for
several iterations rendering slow convergence. We therefore, consider herein
that the best solution also performs a Levy flight so that it can be shifted
to a new position in each iteration if higher fitness can be accomplished. By
means of this modification, all the available nests are given the opportunity
to explore the search space and possibly find more fit solutions.

As observed in (9), the generation of solutions is accomplished by means of
$n$-dimensional random walks. The best solution
$\mathbf{z}^{(t)}_{\text{best}}=\big{[}z^{(t)}_{\text{best}}(1),z^{(t)}_{\text{best}}(2),\dots,z^{(t)}_{\text{best}}(n)\big{]}^{\top}$
at iteration $t$, is used as reference for deriving new candidate solutions
$\mathbf{z}^{(t+1)}_{i}=\big{[}z^{(t+1)}_{i}(1),z^{(t+1)}_{i}(2),\dots,z^{(t+1)}_{i}(n)\big{]}^{\top}$
(for $i=1,2,\dots,N$), which may replace the older solutions
$\mathbf{z}^{(t)}_{i}$ if higher fitness is achieved. The random walks are
based on steps
$\mathbf{w}_{Levy}=\big{[}w_{Levy}(1),w_{Levy}(2),\dots,w_{Levy}(n)\big{]}^{\top}$,
which are drawn from a Levy distribution function. Additionally, a scaling
factor $\alpha$ must be properly chosen so as to be consistent with the scale
of the problem. Equation (9) contemplates two cases where the operator $\odot$
represents element-wise multiplication.

$\small\mathbf{z}^{(t+1)}_{i}=\begin{cases}\mathbf{z}^{(t)}_{i}+\alpha\Big{(}\mathbf{z}^{(t)}_{\text{best}}-\mathbf{z}^{(t)}_{i}\Big{)}\odot\mathbf{w}_{Levy}&\mbox{if
}\mathbf{z}^{(t)}_{i}\neq\mathbf{z}^{(t)}_{\text{best}}\\\
\mathbf{z}^{(t)}_{i}+\alpha\mathbf{w}_{Levy}&\mbox{if
}\mathbf{z}^{(t)}_{i}=\mathbf{z}^{(t)}_{\text{best}}\end{cases}\normalsize$
(9)

The objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ takes as input a
$n$-dimensional candidate solution and maps it to a real value. The fitness is
defined by $F_{\mathbf{z}_{i}}=\frac{1}{1+f(\mathbf{z}_{i})}$

###  III-B Mechanism 2: Best triad mating

Figure 4: Best triad mating process

As already explained, this mechanism aims at combining the top three
solutions—at the ongoing iteration $t$—fine a smaller search space wherein a
better solution is very likely to lie. The essence of this approach is to
concentrate efforts in a potentially richer but smaller subspace. The process
consists of the following five steps.

• Step 1: Let
$\Omega_{\text{Q}}=\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$ denote
a subset that consists of the fittest three nests, which have been sorted
based on their fitness. Further, let $\Omega_{\text{P}}$ be a subset of
$\Omega$ representing the complement of $\Omega_{\text{Q}}$. Using the
elements of $\Omega_{\text{Q}}$, we define a hyperplane $\mathcal{H}$ with
point-to-point distances
$\|\mathbf{q}_{1}-\mathbf{q}_{2}\|,\|\mathbf{q}_{2}-\mathbf{q}_{3}\|,\|\mathbf{q}_{3}-\mathbf{q}_{1}\|$
as shown in Fig.4. These three solutions
$\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$ achieve the highest
fitness in the ongoing generation $t$. However, potentially fitter solutions
might lie in a neighboring area to them.

• Step 2: Given the hyperplane $\mathcal{H}$, we calculate the parameters
$\mathbf{\Gamma}$ and $\varepsilon$ of the circumference inscribed in it. We
can readily calculate these values employing (10) and (11), where $\|.\|$
represents the Euclidean distance.

$\small\mathbf{\Gamma}=\frac{\displaystyle\sum^{2}_{i=0}\bigg{\|}\mathbf{q}_{i-3\big{\lfloor}\frac{i+1}{3}\big{\rfloor}+2}-\mathbf{q}_{i-3\big{\lfloor}\frac{i+2}{3}\big{\rfloor}+3}\bigg{\|}\mathbf{q}_{i-3\big{\lfloor}\frac{i}{3}\big{\rfloor}+1}}{\displaystyle\sum^{2}_{i=0}\bigg{\|}\mathbf{q}_{i-3\big{\lfloor}\frac{i}{3}\big{\rfloor}+1}-\mathbf{q}_{i-3\big{\lfloor}\frac{i+1}{3}\big{\rfloor}+2}\bigg{\|}}\normalsize$
(10)

$\varepsilon=\sqrt{8\displaystyle\prod_{i=0}^{2}\left\\{\frac{1}{2}-\frac{\bigg{\|}\mathbf{q}_{i-3\big{\lfloor}\frac{i+1}{3}\big{\rfloor}+2}-\mathbf{q}_{i-3\big{\lfloor}\frac{i+2}{3}\big{\rfloor}+3}\bigg{\|}}{\displaystyle\sum^{2}_{i=0}\bigg{\|}\mathbf{q}_{i-3\big{\lfloor}\frac{i}{3}\big{\rfloor}+1}-\mathbf{q}_{i-3\big{\lfloor}\frac{i+1}{3}\big{\rfloor}+2}\bigg{\|}}\right\\}}$
(11)

• Step 3: Let $F_{1}$, $F_{2}$ and $F_{3}$ be the fitness values that result
upon evaluating $\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$ in
(LABEL:e10), respectively. Since $\mathbf{q}_{i}$ are sorted according to
their fitness, then $F_{1}\geq F_{2}\geq F_{3}$. Also, let
$\mathbf{\Gamma}_{*}$ be computed by the weighted average shown in (12).

$\mathbf{\Gamma}_{*}=\frac{F_{2}+F_{3}}{F_{1}+F_{2}+F_{3}}\mathbf{\Gamma}+\frac{F_{1}}{F_{1}+F_{2}+F_{3}}\mathbf{q}_{1}$
(12)

The purpose of calculating $\mathbf{\Gamma}_{*}$ is to generate a suitable
reference that will allow to find more fit solutions in the vicinity. Notice
that when $F_{1}\geq F_{2}+F_{3}$, the solution $\mathbf{q}_{1}$ is weighted
with higher importance than $\mathbf{\Gamma}$. This means that because of the
higher fitness of $\mathbf{q}_{1}$ compared with the other two solutions,
$\mathbf{\Gamma}_{*}$ will be dragged nearer towards $\mathbf{q}_{1}$. An
extreme situation occurs when $F_{2}\approx 0$, $F_{3}\approx 0$ in which case
$\mathbf{\Gamma}\approx\mathbf{q}_{1}$. Because $F_{1}\geq F_{2}\geq F_{3}$
the other extreme situation occurs when $F_{1}=F_{2}=F_{3}$ in which case the
weight assignation becomes $2/3$ and $1/3$ for $\mathbf{\Gamma}$ and
$\mathbf{q}_{1}$, respectively.

• Step 4: Three new solutions
$\Omega_{\text{R}}=\\{\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}\\}$ are
generated using $n$-dimensional Gaussian random walks
$\mathbf{w}_{\text{Gaussian}}=[w_{\text{Gaussian}}(1),w_{\text{Gaussian}}(2),\dots,w_{\text{Gaussian}}(n)]^{\top}$
as shown in (13).

$\begin{cases}\mathbf{r}^{(0)}_{\ell}=\mathbf{\Gamma}_{*}\\\
\mathbf{r}^{(t+1)}_{\ell}=\mathbf{r}^{(t)}_{\ell}+\kappa_{\ell}\mathbf{w}_{\text{Gaussian}}&\mbox{if
}t>0\end{cases}$ (13)

Each of the newly generated solutions will be scaled by a weighting factor
$\kappa_{\ell}$ which is defined in (14)

$\kappa_{\ell}=\ell\frac{\varepsilon}{2}\left(v+\phi\left(\max\bigg{\\{}\frac{F_{1}-(F_{2}+F_{3})}{F_{1}+F_{2}+F_{3}},0\bigg{\\}}\right)\right)$
(14)

where $\ell=1,2,3$ and $\varepsilon$ was already defined in (11). In addition,
suitable values for $\phi$ and $\nu$ are comprehended in the ranges
$0.10\leq\phi\leq 0.45$, $0.55\leq\nu\leq 0.90$111These ranges of values have
been obtained via Monte Carlo simulation on the following benchmark functions:
Hyperdimensional sphere [17], Ackley [18], Michalewicz [19], Griewank [20],
and Easom[21]. . By generating new solutions
$\\{\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}\\}$ in the environs of
$\mathbf{\Gamma}_{*}$ the search is confined to a smaller but potentially
richer space. The process described herein aims to increase the exploitation
capability of the algorithm and thus accelerate the convergence rate.

• Step 5: Let the set $\Omega_{\text{S}}$ be defined by
$\Omega_{\text{S}}=\Omega_{\text{Q}}\cup\Omega_{\text{R}}$. Notice that
$\Omega_{\text{S}}$ is composed of both, the solutions
$\\{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}\\}$ and the newly generated
$\\{\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}\\}$. Let the elements of
$\Omega_{\text{S}}$ be sorted in ascending order of fitness value, such that
$F_{\mathbf{s}_{i}}\geq F_{\mathbf{s}_{i+1}}$ for all
$\mathbf{s}_{i}\in\Omega_{\text{S}}$. We redefine the subset
$\Omega_{\text{Q}}$ such that it contains the best three solutions of
$\Omega_{\text{S}}$, i.e.
$\Omega_{\text{Q}}=\\{\mathbf{s}_{1},\mathbf{s}_{2},\mathbf{s}_{3}\\}$.
Finally, we let $\Omega$ be the union of $\Omega_{\text{Q}}$ and
$\Omega_{\text{P}}$, i.e., $\Omega=\Omega_{\text{Q}}\cup\Omega_{\text{P}}$.
Notice that the cardinality of $\Omega$ has not changed, i.e., within each
subprocess several candidate solutions are generated and evaluated but only
$N$ are retained to the next generation.

###  III-C Mechanism 3: Exploitation of the best nest

The aim of the mechanism proposed here is to obtain a benefit from exploiting
the best known solution in each generation in order to approach the global
optimum as much as possible. However, it is impossible to determine whether
the best known solution—at any iteration $t$—is more proximate to a local or
global optimum. If it is closer to a local optimum, the convergence will be
local and the process will probably lead to a biased solution. Nonetheless, we
expect the fourth mechanism—described in the following subsection—to break the
local tendency and help the solution escape towards the global optimum or a
more fit local optimum. On the other hand, if the best known solution is
located in the proximity of the global optimum, the third mechanism will
improve the solution towards the global optimum. We assume that a number of
$N_{\text{G}}$ eggs are dumped in the best nest so that the future broods—in
the following iteration—will search new nests in the neighborhood. Conversely
to the method described in the second mechanism, we consider here the
information of the best nest only. Another difference is that uniform random
walks are used instead of Gaussian random walks. In this process, the worst
performance nests are not necessarily replaced by better ones as this will
cause degradation of the diversity of solutions and therefore produce biasing.
Instead, we randomly choose a number of $N_{\text{G}}$ solutions from
$\Omega$—avoiding the best solution $\mathbf{z}_{\text{best}}$—and replace
them by the newly generated nests
$\Omega_{\text{G}}=\\{\mathbf{g}_{1},\mathbf{g}_{2},\dots,\mathbf{g}_{N_{\text{G}}}\\}$.
This process is divided into two steps which are described below.

• Step 1: A random integer $N_{\text{G}}$ is drawn from the probability mass
function defined in (15). $N_{\text{G}}$ is the number of solutions—randomly
selected from $\Omega$—which are to be replaced by new solutions
$\Omega_{\text{G}}=\\{\mathbf{g}_{1},\mathbf{g}_{2},\dots,\mathbf{g}_{N_{\text{G}}}\\}$.
The solutions $\mathbf{g}_{i}\in\Omega_{\text{G}}$ are generated by (16).

$p(N_{\text{G}}=\zeta)=\begin{cases}0.1&\mbox{, }\zeta=0\\\ 0.2&\mbox{,
}\zeta=1\\\ 0.6&\mbox{, }\zeta=2\\\ 0.1&\mbox{, }\zeta=3\end{cases}$ (15)

• Step 2: By means of $n$-dimensional random walks
$\mathbf{w}_{\text{Uniform}}=[w_{\text{Uniform}}(1),w_{\text{Uniform}}(2),\dots,w_{\text{Uniform}}(n)]^{\top}$,
new $N_{\text{G}}$ solutions are generated as shown in (16), for
$\ell=1,2,\dots,N_{\text{G}}$. Also, consider that
$\mathbf{\Psi}_{1}=[\psi_{1}(1),\psi_{1}(2),\dots,\psi_{1}(n)]^{\top}$ is a
vector whose elements are obtained from a Gaussian probability density
function with mean zero and unit variance.

$\begin{cases}\mathbf{g}^{(0)}_{\ell}=\mathbf{z}_{\text{best}}\\\
\mathbf{g}^{(t+1)}_{\ell}=\mathbf{g}^{(t)}_{\ell}+\big{(}\mathbf{g}^{(0)}+0.25\mathbf{\Psi}_{1}\big{)}\odot\mathbf{w}_{\text{Uniform}}&\mbox{if
}t>0\end{cases}$ (16)

###  III-D Mechanism 4: Mutation of the worst nests

This is a simple method of randomizing and diversifying the candidate
solutions. We form the set
$\Omega_{\text{H}}=\\{\mathbf{h}_{1},\mathbf{h}_{2},\dots,\mathbf{h}_{N_{\text{H}}}\\}$
with the worst $N_{\text{H}}$ solutions from the set $\Omega$. Each element of
$\Omega_{\text{H}}$ is to be updated by means of (17) only if the fitness of
the new solution $\mathbf{h}^{(t+1)}_{\ell}$ is higher than that of the
previous $\mathbf{h}^{(t)}_{\ell}$, where $\ell=1,2,\dots,N_{\text{G}}$.

$\mathbf{h}^{(t+1)}_{\ell}=\frac{1}{2}\Big{(}\mathbf{h}^{(t)}_{\ell}+\mathbf{h}^{(t)}_{\text{worst}}\Big{)}\odot\mathbf{\Psi}_{2}\odot\mathbf{\Psi}_{3}$
(17)

Consider that
$\mathbf{h}_{\text{worst}}=[h_{\text{worst}}(1),h_{\text{worst}}(2),\dots,h_{\text{worst}}(n)]^{\top}$
represents the solution with the lowest fitness at iteration $t$. The elements
constituting
$\mathbf{\Psi}_{2}=[\psi_{2}(1),\psi_{2}(2),\dots,\psi_{2}(n)]^{\top}$ are
obtained from a random variable uniformly distributed in the range [1,2]. On
the other hand, the elements constituting
$\mathbf{\Psi}_{3}=[\psi_{3}(1),\psi_{3}(2),\dots,\psi_{3}(n)]^{\top}$ are
discrete and can be either -1 or 1 with equal probability.

At the beginning of this section, we mentioned that BTM-CSA temporarily
generates $2N+N_{\nu}$ solutions at the processing layer after which only $N$
are preserved to the next generation. Bear in mind that $2N$ nests are due to
its compatibility with CSA. The difference
$N_{\nu}=3+N_{\text{G}}+N_{\text{H}}$ is due to the mechanisms proposed
herein.222Mechanism 1 does not introduce additional solutions.

  * •

$N$ nests due to Levy flights (borrowed from CSA)

  * •

3 nests due to mechanism 2

  * •

$N_{\text{G}}$ nests due to mechanism 3333Note that
$\mathop{\mathbb{E}}\\{N_{\text{G}}\\}=1.7$ but, for practicality
$N_{\text{G}}=2$. See (15).

  * •

$N_{\text{H}}$ nests due to mechanism 4

  * •

$N$ due to selective randomization (borrowed from CSA)

##  IV Adapting BTM-CSA to SVD/CS-PTS

Most of the time, metaheuristics cannot be applied straightforwardly. It is
necessary to make adjustments so that the approach matches with the system
under evaluation. In SVD/CS-PTS the number of dimensions $n$ is strictly equal
to $2M$, $M$ phase rotations and $M$ time shifts. However, without loss of
generality one of the phase rotations and one of the time shifts can be fixed.
The reason is that reduction of PAPR is directly related to phase difference
amongst subcarriers rather than absolute phase. We therefore assume that
$\tilde{\gamma}_{0}=1$ and $\tilde{\delta}_{0}=0$. Thus, the unknown
parameters to be estimated are $\\{\tilde{\gamma}_{m}\\}_{m=1}^{M-1}$ and
$\\{\tilde{\delta}_{m}\\}_{m=1}^{M-1}$, i.e., $n=2(M-1)$. Then any candidate
solution at iteration $t$ has the following structure,
$\mathbf{z}^{(t)}_{i}=[\gamma^{(t)}_{0},\gamma^{(t)}_{1},\dots,\gamma^{(t)}_{M-1},\delta^{(t)}_{0},\delta^{(t)}_{1},\dots,\delta^{(t)}_{M-1}]^{\top}$
with $\tilde{\gamma}^{(t)}_{0}=1$ and $\tilde{\delta}^{(t)}_{0}=0$.

The objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ takes as input an
$n$-dimensional candidate solution and maps it to a real value, which in our
case is the maximum PAPR across all the $N_{T}$ transmit antennas, after phase
rotation and time shifting (See argument of (8)) . The fitness of each
candidate solution is evaluated via (LABEL:e10).

As observed in (8), the evaluation of each candidate solution requires $M$
complex operations due to weighting by $\gamma_{m}$. In strict sense, $M-1$
complex multiplication because we have assumed that $\gamma_{0}$ is fixed and
equal to one. Note that time shifting by $\delta_{m}$ can be accomplished by
altering the summation index only. Therefore, much of the computational
complexity is due to complex multiplications by $\gamma_{m}$. So, in order to
partly avoid these operations we have divided BTM-CSA into two stages as shown
in Fig. 5. In order to decrease the number of complex multiplications, we have
created $\rho_{1}$ intermediate nodes which bear the effect of the phase
rotations $\gamma_{m}$ only. From each intermediate node, $\rho_{2}$ solutions
bearing the added effect of time shifting are generated, amounting a total of
$\rho=\rho_{1}\rho_{2}$. Had we chosen to generate the candidate solutions in
a straightforward way as shown in Fig. 3, we would have had to evaluate
complex multiplications for each of the $\rho$ candidate solutions. Thus, by
adopting the described configuration, complex multiplications are associated
to only $\rho_{1}$ candidate solutions. For the sole purpose of further
reducing complexity—we have assumed that the intermediate nodes are generated
via Algorithm 1 yielding $\rho_{1}=2N+N_{\nu}$ intermediate nodes whereas the
values for $\delta_{m}$ are assigned randomly from the set $\Phi_{\text{D}}$.
Notice that at the intermediate stage, $\\{{\gamma}_{m}\\}_{m=0}^{M-1}$ are
the only parameters involved while $\\{{\delta}_{m}\\}_{m=0}^{M-1}=0$
temporarily, i.e., shifting is not considered. The $\rho_{2}$ solutions that
stem from each intermediate node share the same $\gamma_{m}$ and therefore
only $\delta_{m}$ are varied. We have assumed in our simulations that
$\rho_{2}=3$. We have also considered that $N=11$, $N_{\text{H}}=3$ which
makes $N_{\nu}=8$, $\rho_{1}=2N+N_{\nu}=30$ and $\rho=\rho_{1}\rho_{2}=90$.
That is to say, we evaluate 90 different solutions per each generation. Note
that the complex multiplications in this modified approach are only associated
to the $\rho_{1}=30$ intermediate nodes. Had we used BTM-CSA straightforwardly
or CSA to obtain 90 solutions in the processing layer, we would have required
to perform complex multiplications for each of the 90 candidate solutions.
Thus for the scenario described above, the complex multiplication are reduced
to one third of the total had we used BTM-CSA or CSA directly.

………………………………$N$ nests$\rho_{2}$ nests$\rho_{2}$ nests$\rho_{2}$
nests$\rho_{2}$ nests$\rho=\rho_{1}\rho_{2}$ nests$\rho_{1}$ nests Input layer
Processing layer Intermediate node Figure 5: Adaptation of BTM-CSA

##  V Computational Complexity

The computational complexity of the proposed approach as well as that of other
conventional techniques—such as SLM, PTS and CS-PTS—are described in this
section. Table I, Table II and Table III show the complexity of SLM, PTS and
CS-PTS, respectively. Table IV shows the computational complexity of BTM-CSA.
Table V shows the complexity for other metaheuristic approaches. In these
tables, the number of multiplications and additions of every subprocess
involved is given in detail. Additionally, the nature of
multiplications/additions, i.e., real or complex, has also been specified.

In SLM, $V$ different phase rotation patterns are applied to the input data so
as to generate $V$ different representations of it but with different PAPR
values. In this approach, the number of bits necessary to carry the side
information is equal to $\log_{2}V$. Table I summarizes the processes involved
in SLM and the associated complexity of each.

TABLE I: Computational Complexity of SLM for $V$ Different Patterns Process | Complexity
---|---
Multiplications | Additions
Zero-padded IFFTs | $V\Big{[}\frac{N_{c}L}{2}\log_{2}N_{c}+\frac{N_{c}L}{2}\Big{]}$ | $V\Big{[}N_{c}L\log_{2}N_{c}\Big{]}$
Type | Complex | Complex
Phase patterns | $V[N_{c}L]$ | $0$
Type | Complex | Complex
PAPR computation | $V[2N_{c}L]$ | $V[N_{c}L]$
Type | Real | Real

In PTS, each of the $M$ partitions is weighted by one of the $U$ possible
phase rotations of the set $\Phi_{\text{U}}$. It can be noticed that the total
number of phase patterns that yield from such a configuration is $U^{M}$.
However, without loss of performance, only one value can be fixed, yielding
$U^{M-1}$ different patterns. Each pattern will produce a different OFDM
signal whose PAPR is to be computed. The amount of side information to be
transmitted per OFDM symbol is equivalent to $(M-1)\log_{2}U$ bits. Table II
shows the processes required to perform PTS.

TABLE II: Computational Complexity of PTS for $U^{M-1}$ different patterns Process | Complexity
---|---
Multiplications | Additions
Zero-padded IFFTs | $M\Big{[}\frac{N_{c}L}{2}\log_{2}\Big{(}\frac{N_{c}}{M}\Big{)}+\frac{N_{c}L}{2}\Big{]}$ | $M\Big{[}N_{c}L\log_{2}\Big{(}\frac{N_{c}}{M}\Big{)}\Big{]}$
Type | Complex | Complex
Phase patterns | $U^{M-1}\big{[}(M-1)N_{c}L\big{]}$ | $U^{M-1}\big{[}(M-1)N_{c}L\big{]}$
Type | Complex | Complex
PAPR computation | $U^{M-1}\big{[}2N_{c}L\big{]}$ | $U^{M-1}\big{[}N_{c}L\big{]}$
Type | Real | Real

The last conventional technique considered is CS-PTS, which requires the same
number of IFFTs as PTS. The difference, however, is the number of patterns
that can be generated. Besides using phase rotations, CS-PTS admits cyclic
time-domain delay by values that are drawn from the set $\Phi_{\text{D}}$. The
total number of combinations is $(DU)^{M-1}$. The side information requires
$(M-1)[\log_{2}U+\log_{2}D]$ bits. Table III summarizes the processes that CS-
PTS entails.

TABLE III: Computational Complexity of CS-PTS for $(DU)^{M-1}$ different patterns Process | Complexity
---|---
Multiplications | Additions
Zero-padded IFFTs | $M\Big{[}\frac{N_{c}L}{2}\log_{2}\Big{(}\frac{N_{c}}{M}\Big{)}+\frac{N_{c}L}{2}\Big{]}$ | $M\Big{[}N_{c}L\log_{2}\Big{(}\frac{N_{c}}{M}\Big{)}\Big{]}$
Type | Complex | Complex
Phase patterns | $U^{M-1}\big{[}(M-1)N_{c}L\big{]}$ | $(DU)^{M-1}\big{[}(M-1)N_{c}L\big{]}$
Type | Complex | Complex
PAPR computation | $(DU)^{M-1}\big{[}2N_{c}L\big{]}$ | $(DU)^{M-1}\big{[}N_{c}L\big{]}$
Type | Real | Real

The complexity of our proposed approach which merges CS-PTS and BTM-CSA is
shown in Table IV. The number of solution patterns to be generated are
controlled by $\rho_{1}$ and $\rho_{2}$, as explained in Section IV. The
average computational cost for generating a single solution via BTM-CSA is
approximately $3(M-1)$ real multiplications5 and $2(M-1)$ real additions. Note
that there is no cost for generating a solution when using exhaustive search
because the patterns are assembled in a structured manner without additional
processing. When using metaheuristics there is always an extra cost for
generation of solutions. On the other hand, this extra cost in metaheuristics
is compensated by its merit of accomplishing near-optimal solutions without
exploring all solutions. Also, as can be observed in Table IV, the use of BTM-
CSA does not affect the complexity of the IFFTs but only those processes
related to the search process itself.

TABLE IV: Computational Complexity of CS-PTS Using BTM-CSA for $\rho=\rho_{1}\rho_{2}$ Solution Patterns Process | Complexity
---|---
Multiplications | Additions
Zero-padded IFFTs | $M\Big{[}\frac{N_{c}L}{2}\log_{2}(\frac{N_{c}}{M})+\frac{N_{c}L}{2}\Big{]}$ | $M\Big{[}N_{c}L\log_{2}(\frac{N_{c}}{M})\Big{]}$
Type | Complex | Complex
Phase patterns | $\rho_{1}[(M-1)N_{c}L]$ | $\rho_{1}\rho_{2}[(M-1)N_{c}L]$
Type | Complex | Complex
Generation of solutions | $\rho[3(M-1)]$ | $\rho[2(M-1)]$
Type | Real | Real
PAPR computation | $\rho[2N_{c}L]$ | $\rho[N_{c}L]$
Type | Real | Real

Table V shows the complexity for other competing techniques such as GA, PSO
and the original CSA. Notice that the complexity of BTM-CSA differs from that
of the other approaches. Recall that in order to reduce the complexity, we
have divided BTM-CSA to operate in a decoupled manner, i.e. the phase rotation
stage is independent from the time-domain shifting stage. Hence, in BTM-CSA
the total number of multiplications due to generation of patterns is
proportional to $\rho_{1}$ whereas it is proportional to
$\rho=\rho_{1}\rho_{2}$ in the case of GA, PSO and CSA. By inspecting Table V,
we realize that the generation of solutions is more costly in CSA than it is
for GA and PSO, which is a tentative reason to explain why CSA exhibits faster
convergence rate, i.e. the candidate solutions generated by CSA are more
refined but computationally more expensive. Regarding the process of
generation of solutions, GA has the lowest computational cost among the three
competing approaches. PSO has an average complexity between CSA and GA. On the
other hand, Table IV shows the average complexity of BTM-CSA which is even
lower than that of GA. Complexity is of course a critical factor when deciding
among several approaches but this decision should also be supported by the
convergence capabilities of the chosen method. For instance, although GA seems
to be a fair election due to its lower complexity, we will corroborate in the
next section that GA is indeed the approach with the worst convergence
performance.

TABLE V: Comparison of Computational Complexity of Metaheuristic Techniques
for $\rho$ Solution Patterns

Process | Complexity
---|---
Multiplications | Additions
Zero-padded IFFTs | $M\Big{[}\frac{N_{c}L}{2}\log_{2}(\frac{N_{c}}{M})+\frac{N_{c}L}{2}\Big{]}$ | $M\Big{[}N_{c}L\log_{2}(\frac{N_{c}}{M})\Big{]}$
Type | Complex | Complex
Phase patterns | $\rho[(M-1)N_{c}L]$ | $\rho[(M-1)N_{c}L]$
Type | Complex | Complex
PAPR computation | $\rho[(M-1)N_{c}L]$ | $\rho[(M-1)N_{c}L]$
Type | Complex | Complex
GA |  Generation of solutions | $\rho[3(M-1)]$ | $\rho[3(M-1)]$
2-4 |  Type | Complex | Complex
PSO |  Generation of solutions | $\rho[5(M-1)]$ | $\rho[5(M-1)]$
2-4 |  Type | Complex | Complex
CSA |  Generation of solutions | $\rho[8(M-1)]$ | $\rho[4(M-1)]$
2-4 |  Type | Complex | Complex

##  VI Numerical Experiments

The effect of varying the number of antennas on the PAPR distribution of
SVD/MIMO-OFDM systems can be observed in Fig. 6. For any $\text{PAPR}_{0}$ in
the abscissa, its corresponding value in the ordinate represents the
probability that the PAPR of an OFDM signal exceeds the threshold
$\text{PAPR}_{0}$—this holds true for SISO systems. In the case of SVD/MIMO-
OFDM, since there are $N_{T}$ transmit antennas, $\text{PAPR}_{0}$ represents
the maximum PAPR of all the transmit antennas. Note that as the number of
antennas grows, the maximum PAPR across all the antennas increases as well.
Bear in mind that no PAPR reduction technique was employed in this scenario.

Since the PAPR defined as
$\max\\{\text{PAPR}^{(1)},\text{PAPR}^{(2)},\dots,\text{PAPR}^{(N_{T})}\\}$
does not vary much with increasing number of transmit antennas, we consider in
the rest of our simulations the 4x4 configuration only.

Fig. 7 displays the PAPR reduction performance of SLM employing a wide number
of patterns $V=4,16,64,256,1024,4096$. The SLM patterns have been generated
with phase rotations pseudo-randomly selected from the set
$\\{e^{0},e^{\pi}\\}$. Fig. 8 shows the performance of PTS for several
configurations of $M=2,4,8$ and $U=2,4,8$. On the other hand, Fig. 9 shows
other scenarios of CS-PTS with $M=2,4,8$, $U=2,4,8$ and $D=2,4,8$.

$4$$5$$6$$7$$8$$9$$10$$11$$12$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$SISO2x2 SVD-MIMO3x3 SVD-MIMO4x4 SVD-
MIMO6x6 SVD-
MIMO$4$$5$$6$$7$$8$$9$$10$$11$$12$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$SISO2x2 SVD-MIMO3x3 SVD-MIMO4x4 SVD-
MIMO6x6 SVD-MIMO Figure 6: Effect of the number of transmit antennas on the
PAPR
$4$$5$$6$$7$$8$$9$$10$$11$$12$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$4x4 SVD w/o optimization4x4 SVD / SLM /
V = 44x4 SVD / SLM / V = 164x4 SVD / SLM / V = 644x4 SVD / SLM / V = 2564x4
SVD / SLM / V = 10244x4 SVD / SLM / V = 4096 Figure 7: PAPR reduction
performance via SLM
$4$$5$$6$$7$$8$$9$$10$$11$$12$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$4x4 SVD w/o optimization4x4 SVD / PTS /
M = 2 / U = 24x4 SVD / PTS / M = 2 / U = 44x4 SVD / PTS / M = 2 / U = 84x4 SVD
/ PTS / M = 4 / U = 24x4 SVD / PTS / M = 4 / U = 44x4 SVD / PTS / M = 8 / U =
24x4 SVD / PTS / M = 8 / U = 4 Figure 8: PAPR reduction performance via PTS
$4$$5$$6$$7$$8$$9$$10$$11$$12$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$4x4 SVD w/o optimization4x4 SVD / CS-
PTS / M = 2 / U = 2 / D = 24x4 SVD / CS-PTS / M = 2 / U = 2 / D = 44x4 SVD /
CS-PTS / M = 2 / U = 4 / D = 44x4 SVD / CS-PTS / M = 4 / U = 2 / D = 24x4 SVD
/ CS-PTS / M = 4 / U = 2 / D = 44x4 SVD / CS-PTS / M = 4 / U = 4 / D = 44x4
SVD / CS-PTS / M = 4 / U = 4 / D = 8 Figure 9: PAPR reduction performance via
CS-PTS

Considering the trade-off between complexity and performance, CS-PTS can
achieve superior results if compared with PTS and SLM. The reasons is that in
CS-PTS several of the patterns are generated by solely shifting the partial
transmit sequences, which does not incur additional complex multiplications.

In Fig. 10 we have considered the results of CS-PTS using $M=4$, $U=4$, $D=8$
(optimal) as a reference for comparison. In addition, the performance curves
of CSA and BTM-CSA are shown in order to analyze the difference between the
achievable performance by metaheuristics. We have considered CSA as the
competing approach because although it is computationally more complex than
its counterparts PSO and GA, its convergence is faster in terms of iterations.
We have considered a variety of iterations $\rho=270,540,1080,2160,4320$. We
can observe that for the same number of iterations BTM-CSA outperforms CSA.
Although the difference of performance among these two approaches may seem
negligible we shall recall that the complexity of BTM-CSA is less than half of
CSA’s. Also, to attain the optimal PAPR curve $(4\times 8)^{3}$ patterns have
been tested using exhaustive search whereas only $4320$ patterns have been
evaluated in BTM-CSA approach. At the probability of $10^{-4}$ the performance
difference between the optimal and BTM-CSA is roughly 0.1dB.

We may not expect PSO or GA to outperform CSA since they have lower converge
rate. However, it is interesting to compare these approaches to observe to
what extent their curves fall away from the optimal and BTM-CSA. Fig. 11
displays the performance of BTM-CSA against CSA, PSO and GA. We can observe
that GA exhibits the worst performance among the techniques whereas PSO shows
a significant improvement over GA but still behind CSA. Although BTM-CSA and
GA have both comparable complexities, BTM-CSA still outperforms GA by 0.15 dB
at a probability of $10^{-4}$ and this is steadily maintained for
$\rho=270,1080,4320$.

$6$$6.2$$6.4$$6.6$$6.8$$7$$7.2$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$Exhaustive Search (Optimal)CSA /
$\rho=270$BTM-CSA / $\rho=270$CSA / $\rho=540$BTM-CSA / $\rho=540$CSA /
$\rho=1080$BTM-CSA / $\rho=1080$CSA / $\rho=2160$BTM-CSA / $\rho=2160$CSA /
$\rho=4320$BTM-CSA / $\rho=4320$ Figure 10: Comparison of metaheuristic
approaches
$6$$6.2$$6.4$$6.6$$6.8$$7$$7.2$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$
$\rho=4320$ $\rho=1080$ $\rho=270$ $\text{PAPR}_{0}$
(dB)Prob$(\text{PAPR}>\text{PAPR}_{0})$Exhaustive Search (Optimal)BTM-
CSACSAPSOGA Figure 11: Comparison of metaheuristic approaches

##  VII Conclusions

In this paper we proposed a modified optimization technique called BTM-CSA
which exhibits faster convergence rate than the conventional CSA. BTM-CSA was
also tested against GA and PSO and it was verified that it outperforms both
approaches in terms of convergence rate and at lower complexity. We can
conclude that the proposed approach achieves a performance which is close to
the optimal and without exploring all the solution space as occurs when using
exhaustive search. If we count with a system with limited computational
capability, it is more suitable to explore the solution space in a smart and
randomized fashion rather than straightforwardly as occurs with exhaustive
search. In this regard, metaheuristics can provide a good balance between the
optimality of the solutions and computational complexity.

## References

  * [1] A. Batra and J. R. Zeidler, “Narrowband interference mitigation in OFDM systems,” in Proceedings of the IEEE Military Communications Conference (MILCOM), San Diego, CA, Nov 2008.
  * [2] S. C. Thompson, A. U. Ahmed, J. G. Proakis, J. R. Zeidler, and M. J. Geile, “Constant envelope OFDM,” IEEE Transactions on Communications, Vol. 56, No. 8, August 2008.
  * [3] M. Ergen, Mobile Broadband Including WiMAX and LTE, pp. 158-160, Springer, 2009.
  * [4] L. Cimini, “Analysis and Simulation of a digital mobile channel using OFDM,” IEEE Transactions on Communications, Vol. 33, pp. 665-675, July 1985.
  * [5] E. Costa, M. Midro, and S. Pupolin, “Impact of amplifier nonlinearities onOFDM transmission system performance,” IEEE Communications Letters, Vol. 3, pp. 37–39, February 1999.
  * [6] X. B. Wang, T. T. Tjhung, and C. S. Ng, “Reduction of peak-to-average power ratio of OFDM system using A companding technique,” IEEE Transactions on Broadcasting, Vol. 45, No. 3, pp. 303–307, September 1999.
  * [7] J. Tellado, “Peak to Average Power Ratio Reduction for Multicarrier Modulation,” PhD thesis, University of Stanford, Stanford, 1999.
  * [8] R. W. Baml, R. F. H. Fisher, and J. B. Huber, “Reducing the peak-to average power ratio of multi-carrier modulation by selective mapping,” Electronics Letters, Vol. 32, No. 22, pp. 2056-2057, October 1996.
  * [9] S. H. Muller and J. B. Huber, “OFDM with reduced peak-to average power ratio by optimum combining of partial transmit sequences,” Electronic Letters, Vol. 33, No. 5, pp. 368–369, February 1997.
  * [10] G. R. Hill, M. Faulkner, and J. Singh, “Reducing the peak-to-average power ratio in OFDM by cyclically shifting partial transmit sequences,” Electronics Letters, Vol. 36, No. 6, pp. 560–561, March 2000.
  * [11] X.-S. Yang and S. Deb, “Cuckoo Search via Levy Fights,” World Congress on Natural and Biologically Inspired Computing, pp. 210–214, 2009.
  * [12] J. Kennedy and R. Eberhart, “Particle Swarm Optimization,” Proceedings of IEEE International Conference on Neural Networks IV, pp. 1942–1948, 1995.
  * [13] J. H. Holland. Adaptation in natural and artificial systems. The University of Michigan Press, Ann Arbor, 1975.
  * [14] W. D Hamilton, Geometry for the Selfish Herd, Journal of Theoretical Biology, Vol. 31, pp. 295–311, 1970.
  * [15] X.-S. Yang and S. Deb, “Engineering optimization by cuckoo search,” International Journal of Mathematical Modeling and Numerical Optimization, Vol. 1, No. 4, 330–343, 2010.
  * [16] R. Storn and K. Price, ”Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces,” Journal of Global Optimization, pp. 341–359, 1997.
  * [17] L. A. Rastrigin, ”The convergence of the random search method in the extremal control of a many parameter system”, Automation and Remote Control, Vol. 10, pp. 1337–1342, 1963.
  * [18] D. H. Ackley, “A Connectionist Machine for Genetic Hill-Climbing,” Kluwer, 1987.
  * [19] Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs. Berlin, Heidelberg, New York: Springer-Verlag, 1992.
  * [20] A. O. Griewank, ”Generalized Descent for Global Optimization” Journal of Optimization Theory and Applications, Vol. 34, No. 1, pp. 11-39, 1981.
  * [21] E. E. Easom, A survey of global optimization techniques, M. Eng. thesis, University of Louisville, Louisville, KY, 1990.