# Two competing populations with a common environmental resource ††thanks:

Keith Paarporn, James Nelson K. Paarporn and J. Nelson are with the Department
of Computer Science, University of Colorado, Colorado Springs. Contact:
<EMAIL_ADDRESS>This work is supported in part by the UCCS
Committee on Research and Creative Works.

###### Abstract

Feedback-evolving games is a framework that models the co-evolution between
payoff functions and an environmental state. It serves as a useful tool to
analyze many social dilemmas such as natural resource consumption, behaviors
in epidemics, and the evolution of biological populations. However, it has
primarily focused on the dynamics of a single population of agents. In this
paper, we consider the impact of two populations of agents that share a common
environmental resource. We focus on a scenario where individuals in one
population are governed by an environmentally “responsible” incentive policy,
and individuals in the other population are environmentally “irresponsible”.
An analysis on the asymptotic stability of the coupled system is provided, and
conditions for which the resource collapses are identified. We then derive
consumption rates for the irresponsible population that optimally exploit the
environmental resource, and analyze how incentives should be allocated to the
responsible population that most effectively promote the environment via a
sensitivity analysis.

##  I Introduction

Classical game theory describes strategic interactions between individuals
provided that the incentives for their choices are static, i.e. do not change
over time. However, their choices often have an impact on a shared
environment, which in turn may affect their incentives for future decisions.
The utilization of common resources such as water, fishing grounds, or traffic
networks best illustrates this – high individual utilization makes fewer
resources available to others in the future [1]. Clearly, a dynamic interplay
exists between users’ decisions and their impact on a shared environmental
state.

The recent framework of _feedback-evolving games_ addresses this interplay by
incorporating a changing environmental state coupled with existing
evolutionary game theoretic dynamics [2, 3, 4]. The core model considers
payoff functions to individuals that depend on whether the environment is in
an abundant or depleted state. It is primarily used to understand how the
payoff functions should be structured in order to avoid a “tragedy of the
commons”, an outcome in which the environment ultimately becomes depleted. In
other words, these payoff functions reflect environmental incentive policies
that can be designed to manage the common resource [5, 6, 7].

Much of the existing research focuses on a _single population_ of agents whose
decisions impact the environment in isolation [8, 9, 10, 7, 11, 12, 13]. These
models are not sufficient to describe multi-population interactions. For
example, individuals in neighboring countries follow different environmental
policies, yet utilize resources from the same common source (e.g. fish in the
ocean, water from rivers, clean air). Extensions of feedback-evolving games to
multiple populations, and possibly multiple local environments, would enable a
vastly richer set of scenarios for study. While some work has studied the
dynamics of multi-population models in recent years [14, 15, 16], an
interesting and understudied direction involves hierarchical decision-making,
i.e. the strategic selection of local incentive policies, under the context
that the activity from one population has externalities on resources utilized
by other populations.

In this paper, we extend the framework in this direction by considering two
populations that share the same local environmental resource. We focus on a
particular setting where one population, labelled “responsible”, implements a
pro-environmental incentive policy such that the common resource can be
sustained in the absence of other populations. The other population, labelled
“irresponsible”, does not restrain its consumption activity. Our study centers
on the following question, stated informally as

_To what extent can the irresponsible population exploit the environmental
resource?_

While too much consumption could cause the resource to collapse, too little
consumption leaves missed opportunities. Our highlighted contributions are:

  * •

A stability analysis for the dynamics of the two-population system.

  * •

Identify conditions for which the environment collapses, and when it can be
sustained.

  * •

Derivation of the optimal consumption rates for the irresponsible population.

  * •

A sensitivity analysis for the incentive policies chosen by the responsible
population.

The last item above illustrates how incentives should be allocated in order to
maximally promote the environmental resource. Interestingly, we find that
incentivizing mutual cooperation (i.e. coordinated) is far more effective in
promoting resource levels than incentivizing unilateral cooperation (anti-
coordinated).

The paper is organized as follows. Section II provides relevant background on
single population feedback-evolving games, before presenting the two-
population model. Section III states our main result regarding the dynamics of
the model. Section IV proposes and solves an optimization problem regarding
the irresponsible population’s choice of consumption levels. Simulations and
sensitivity analyses are given in Section V, followed by concluding remarks.

##  II Model

Before describing our two-population model, preliminary background on single-
population feedback-evolving games is provided.

###  II-A Preliminary on Feedback-evolving games

A single-population feedback-evolving game describes a population of agents
that have access to a degradable environmental resource. The relative
abundance of the resource is denoted by $n\in[0,1]$. At any given time, each
agent is choosing whether to use a low consumption action (strategy
$\mathcal{L}$), or a high consumption action (strategy $\mathcal{H}$). High
consumption degrades $n$, and low consumption improves $n$. The immediate
payoff experienced by an agent is described by the environment-dependent
$2\times 2$ payoff matrix,

$A_{n}=n\begin{bmatrix}R_{1}&S_{1}\\\
T_{1}&P_{1}\end{bmatrix}+(1-n)\begin{bmatrix}R_{0}&S_{0}\\\
T_{0}&P_{0}\end{bmatrix}$ (1)

Here, the first row and column corresponds to a low consumer, and the second
row and column corresponds to a high consumer. Entry $ij$
($i,j\in\\{\mathcal{L},\mathcal{H}\\}$) indicates the experienced payoff to an
agent using strategy $i$ when encountering an agent using strategy $j$. We
denote $x\in[0,1]$ as the fraction (or frequency) of agents in the population
using strategy $\mathcal{L}$. The payoff experienced by each type of agent is
then given by

$\pi_{\mathcal{L}}(x,n)=[A_{n}[x,1-x]^{\top}]_{1},\quad\pi_{\mathcal{H}}(x,n)=[A_{n}[x,1-x]^{\top}]_{2}$
(2)

The payoffs are determined by the parameters in the $A_{0}$ and $A_{1}$
matrices. The $A_{1}$ matrix is the payoff matrix when the environment is
abundant. Following the literature on feedback-evolving games, we make the
following assumption about the $A_{1}$ matrix.

Figure 1: Diagram of our two-population model. Activities from both
populations impact the shared environmental state, $n$. The rates
$\theta_{i}$, $\alpha_{i}$ denote the restoration and degradation rates from
the activities of population $i$, respectively.

######  Assumption 1.

High consumption is the dominant strategy in $A_{1}$, i.e.
$\delta_{TR1}:=T_{1}-R_{1}>0$ and $\delta_{PS1}:=P_{1}-S_{1}>0$

On the other hand, the $A_{0}$ matrix describes payoffs when the environment
is depleted. We interpret $A_{0}$ to be an “ _environmental policy_ ” that the
population has follows. For example, the low consumption strategy may be more
incentivized when the environment is bad (e.g. subsidies for using electric
vehicles). The payoff structure of the $A_{0}$ matrix is completely determined
from the parameters $\delta_{SP0}:=S_{0}-P_{0}$ and
$\delta_{RT0}:=R_{0}-T_{0}$.

We will use the replicator equation to describe how agents revise their
decisions over time. Moreover, the environment $n$ evolves over time, as it is
influenced by the decisions in the population. The overall system dynamics is
given by two coupled ODEs:

$\displaystyle\dot{x}$ $\displaystyle=x(1-x)g(x,n)$ (3) $\displaystyle\dot{n}$
$\displaystyle=\epsilon n(1-n)(\theta x-\alpha(1-x))$

where

$g(x,n):=\pi_{\mathcal{L}}(x,n)-\pi_{\mathcal{H}}(x,n)$ (4)

is the payoff difference between low and high consumers, $\theta>0$ is the
restoration rate from low consumption activity, $\alpha>0$ is the degradation
rate from high consumption activity, and $\epsilon>0$ is a time-scale
separation constant. The form of the $\dot{n}$ equation is often referred to
as the _tipping point_ dynamics, since $n$ is increasing only if there are
sufficiently high fraction of low consumers.

The state of the system $(x,n)$ evolves over the state space $[0,1]^{2}$.
Indeed, the set $(0,1)^{2}$ is forward-invariant under the dynamics (3). The
originating work [2] provided a full characterization of the asymptotic
outcomes of system (3) for all possible environmental policy matrices $A_{0}$.

The result below summarizes the variety of behaviors that system (3) can
exhibit.

######  Theorem 2.1 (adapted from [2]).

The environmental policy $(\delta_{SP0},\delta_{RT0})$ determines the
asymptotic properties of (3) as follows.

1) _Sustained resource:_ If $(\delta_{SP0},\delta_{RT0})\in\mathcal{V}$, where

$\mathcal{V}:=\\{(y_{1},y_{2})\in\mathbb{R}^{2}:y_{1}>0\text{ and
}-\frac{\theta}{\alpha}y_{1}<y_{2}<\frac{\delta_{TR1}}{\delta_{PS1}}y_{1}\\},$
(5)

then the fixed point
$(x^{*},n^{*})=(\frac{\alpha}{\alpha+\theta},\frac{g(x^{*},0)}{-\partial
g/\partial n(x^{*})})\in(0,1)^{2}$ is the only asymptotically stable fixed
point in the system.

2) _Oscillating Tragedy of the commons:_ If $\delta_{SP0}>0$ and
$\frac{\delta_{TR1}}{\delta_{PS1}}\delta_{SP0}<\delta_{RT0}$, then system (3)
exhibits a stable heteroclinic cycle between the four corner fixed points,
$(0,0)$, $(1,0)$, $(1,1)$, $(0,1)$.

3) _Tragedy of the commons:_ If $\delta_{SP0}\leq 0$, or $\delta_{SP0}>0$ and
$\delta_{RT0}<-\frac{\theta}{\alpha}\delta_{SP0}$, then the only
asymptotically stable fixed point has $n=0$.

The environmental policies belonging to $\mathcal{V}$ are ones that can
sustain a stable and non-zero resource level. In other words, they avert a
_tragedy of the commons_ [17], which refers to a fixed point where the
resource level is zero. The policies in $\mathcal{V}$ can be considered as
more desirable operating conditions than the policies described in items 2 and
3 above. The oscillating tragedy of the commons is considered an undesirable
outcome, in which the system cycles between states of nearly zero resource
levels, and states of fully abundant resource levels. Visually, the state
trajectory “spirals outward” from the interior, and its $\omega$-limit set is
the boundary of the state space.

###  II-B Model: Two-population feedback-evolving games

Now, we consider two populations that share the environmental resource. Figure
1 illustrates this scenario. Specifically, the consumption decisions of the
members of both populations have effects on the environmental resource. Now,
$x_{i}\in[0,1]$ denotes the fraction of low consumers in population
$i\in\\{1,2\\}$. The behavior of population $i$ is governed by the payoff
matrices $A_{0}^{(i)}$ (parameters $\delta_{SP0}^{i}$, $\delta_{RT0}^{i}$) and
$A_{1}^{(i)}$ (parameters $\delta_{TR1}^{i}$, $\delta_{PS1}^{i}$), restoration
rate $\theta_{i}>0$, degradation rate $\alpha_{i}>0$. In line with Assumption
1, we maintain the parameters in the abundant state satisfy
$\delta_{TR1}^{i},\delta_{PS1}^{i}>0$. The payoff differences for each
population are given by

$g_{i}(x_{i},n):=\pi_{\mathcal{L}}^{(i)}(x_{i},n)-\pi_{\mathcal{H}}^{(i)}(x_{i},n)$
(6)

where $\pi_{j}^{(i)}$ is the experienced payoff to a $j$-strategist in
population $i$. The system dynamics are now given by the set of three ODEs
with state variable $\boldsymbol{z}=(x_{1},x_{2},n)$:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{1}(1-x_{1})g_{1}(x_{1},n)$ (7)
$\displaystyle\dot{x}_{2}$ $\displaystyle=x_{2}(1-x_{2})g_{2}(x_{2},n)$
$\displaystyle\dot{n}$ $\displaystyle=\epsilon n(1-n)h(x_{1},x_{2})$

where $h(x_{1},x_{2}):=\sum_{i=1}^{2}(\theta_{i}x_{i}-\alpha_{i}(1-x_{i}))$,
and with initial condition $(x_{1}(0),x_{2}(0),n(0))\in(0,1)^{3}$. By
construction, $(0,1)^{3}$ is a forward-invariant set. We are interested in
studying how the dynamics of the two-population environmental coupled system
(7) behaves with respect to the population’s defining parameters – especially
how the environmental policies $A_{0}^{(1)}$, $A_{0}^{(2)}$ adopted by each
population affects environmental resources. We will make the following
assumptions on these environmental policies.

######  Assumption 2.

We assume that $\frac{\partial g_{i}}{\partial n}(x_{i})<0$ for all
$x_{i}\in[0,1]$, $i=1,2$. This is equivalent to
$\delta_{SP0}^{i}>-\delta_{PS1}^{i}$ and $\delta_{RT0}^{i}>-\delta_{TR1}^{i}$.

Assumption 2 asserts that the relative payoff to low consumers in both
populations monotonically decreases as the environmental state improves. The
next and final assumption focuses our study on a scenario where one population
is “environmentally responsible”, and the other population is “irresponsible”.

######  Assumption 3.

For population 1, we assume that
$(\delta_{SP0}^{1},\delta_{RT0}^{1})\in\mathcal{V}$ (with $\theta=\theta_{1}$
and $\alpha=\alpha_{1}$). For population 2, we assume that
$\delta_{SP0}^{2},\delta_{RT0}^{2}<0$.

The first part of Assumption 3 asserts that agents in population 1 are
cooperative enough to sustain the resource (Theorem 2.1, item 1) in the
absence of population 2. The second part asserts that agents in population 2
make no effort to conserve resources, always preferring the high consumption
strategy. We will refer to population 1 as _responsible_ , and population 2 as
_irresponsible_. The set of all feasible policies specified by Assumptions 2
and 3 is visually depicted in Figure 2.

Figure 2: The set of all sustainable policies for population 1 is shown as the
blue region. A sustainable policy maintains a stable, nonzero resource in the
absence of other populations (Theorem 2.1). This region is the focus of
Assumptions 2 and 3.

##  III Stability analysis of two-population game

In this section, we characterize the dynamical behavior of the two-population
coupled system by analyzing the local stability properties of all of its fixed
points.

######  Theorem 3.1.

The asymptotic dynamics of the two-population system (7) are summarized below.

1) Suppose $\alpha_{2}>\theta_{1}$. Then a tragedy of the commons is globally
asymptotically stable, i.e. $\lim_{t\rightarrow\infty}n(t)=0$.

2) Suppose $\alpha_{2}\leq\theta_{1}$.

  1. (a)

If
$\frac{\alpha_{2}-\theta_{1}}{\alpha_{1}+\alpha_{2}}\delta_{SP0}^{1}\leq\delta_{RT0}^{1}\leq\frac{\delta_{TR1}}{\delta_{PS1}}\delta_{SP0}^{1}$,
then the only asymptotically stable fixed point is of the form $(x_{1},0,n)$,
where

$x_{1}^{*}:=\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},\quad
n^{*}:=-\frac{g_{1}(x_{1}^{*},0)}{\frac{\partial g_{1}}{\partial
n}(x_{1}^{*})}.$ (8)

  2. (b)

If
$\delta_{RT0}^{1}<\frac{\alpha_{2}-\theta_{1}}{\alpha_{1}+\alpha_{2}}\delta_{SP0}^{1}$,
then the only asymptotically stable fixed point constitutes a tragedy of the
commons.

Several remarks are in order. In item 1 above, the irresponsible population
induces a tragedy of the commons if its consumption rate is higher than the
responsible population’s restoration rate. Item 2a provides a region of
sustainable policies for population 1. We note that this region gets smaller
as $\alpha_{2}$ increases while remaining less than $\theta_{1}$.

Moreover, we observe the sustained resource level $n^{*}$ in item 2a is a
decreasing function in $\alpha_{2}$. Taking the derivative with respect to
$\alpha_{2}$, we obtain

$\frac{\partial
n^{*}}{\partial\alpha_{2}}=\frac{g_{1}(x_{1}^{*},0)\frac{\partial^{2}g_{1}}{\partial
n\partial x}-\frac{\partial g_{1}}{\partial x}(x_{1}^{*},0)\frac{\partial
g_{1}}{\partial n}(x_{1}^{*})}{(\alpha_{1}+\theta_{1})(\frac{\partial
g_{1}}{\partial n}(x_{1}^{*}))^{2}}$

It is negative since the numerator being negative is equivalent to the
condition
$\delta_{RT0}^{1}\leq\frac{\delta_{TR1}}{\delta_{PS1}}\delta_{SP0}^{1}$.

Item 2b provides a region where the population 1 policy fails to sustain the
resource even though $\alpha_{2}<\theta_{1}$.

The proofs for each part of Theorem 3.1 are presented in Appendix -A.

Figure 3: (Left) This surface plot shows the resource level that results from
the optimal consumption rate $\alpha_{2}^{*}$ detailed in Theorem 4.1. The red
line indicates the bottom border of the set of feasible policies for
population 1. We observe that $R(\alpha_{2}^{*})$ is increasing in both
$\delta_{RT0}^{1}$ and $\delta_{SP0}^{1}$. (Center) The optimal consumption
rate detailed in Theorem 4.1. (Right) An example of the utility function
$U_{2}(\alpha_{2})$ under policy
$(\delta_{SP0}^{1},\delta_{RT0}^{1})=(3,-0.5)$. For all simulations, the
parameter values are: $\delta_{TR1}=10$, $\delta_{PS1}=6$, $\theta_{1}=0.75$,
$\alpha_{1}=1$, $\epsilon=0.1$.

##  IV Results: exploitation of resources

In this section, we study the following hierarchical decision problem posed
informally as: _how much consumption can the irresponsible population get away
with?_ To approach this question, we consider a local authority for population
2 that may set the consumption rate $\alpha_{2}\geq 0$. This decision is
representative of, for example, water usage or fishing regulations.

From Theorem 3.1, we succinctly summarize the resource level that results in
the asymptotically stable state for any choice of $\alpha_{2}$ by defining the
_resource function_

$R(\alpha_{2}):=\begin{cases}n^{*}(\alpha_{2}),&\text{if
}\alpha_{2}\leq\theta_{1}\text{ and}\\\
&\quad\frac{\alpha_{2}-\theta_{1}}{\alpha_{1}+\alpha_{2}}\delta_{SP0}^{1}\leq\delta_{RT0}^{1}\leq\frac{\delta_{TR1}}{\delta_{PS1}}\delta_{SP0}^{1}\\\
0,&\text{else }\end{cases}$ (9)

where the function
$n^{*}(\alpha_{2}):=-\frac{g_{1}(\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},0)}{\frac{\partial
g_{1}}{\partial n}(\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}})}$ is
taken directly from (8). The _utility_ that the authority seeks to maximize is
defined by

$U(\alpha_{2}):=\alpha_{2}\cdot R(\alpha_{2}).$ (10)

The choice to increase consumption comes at the cost of worsening or even
destroying the environmental resource. Thus, (10) captures the tension between
resource consumption and the stability of the resource. The optimal
consumption rate can be determined by solving the optimization problem

$\alpha_{2}^{*}:=\arg\max_{\alpha_{2}\geq 0}U(\alpha_{2}).$ (OC)

Note here that we are assuming a fixed environmental policy
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$ for population 1 that satisfies
Assumption 3. Our main result below provides a full characterization of the
optimal consumption rate and utility.

######  Theorem 4.1 (Optimal consumption rate).

The optimal solution and value of (OC) are given as follows.

  1. (a)

If
$C(\delta_{SP0}^{1})\leq\delta_{RT0}^{1}\leq\frac{\delta_{TR1}^{1}}{\delta_{PS1}^{1}}\delta_{SP0}^{1}$,
then

$\alpha_{2}^{*}=\theta_{1}$ (11)

and
$R(\alpha_{2}^{*})=\frac{\delta_{RT0}^{1}}{\delta_{RT0}^{1}+\delta_{TR1}^{1}}$.

  2. (b)

If
$\max\\{-\frac{\theta_{1}}{\alpha_{1}}\delta_{SP0}^{1},-\delta_{TR1}^{1}\\}\leq\delta_{RT0}^{1}\leq
C(\delta_{SP0}^{1})$, then

$\alpha_{2}^{*}=-\frac{(\alpha_{1}+\theta_{1})}{a_{1}}\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})\left(1-\sqrt{\frac{Y}{b_{1}\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})}}\right)$ (12)

and $R(\alpha_{2}^{*})=\frac{\sqrt{b_{1}\frac{\partial g_{1}}{\partial
n}(\bar{\alpha}_{1})}-\sqrt{Y}}{a_{1}\sqrt{b_{1}^{-1}\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})}}$, where we have defined

$\displaystyle
C(\delta_{SP0}^{1}):=\frac{1}{2}\left[-((1-\bar{\alpha}_{1})\delta_{PS1}^{1}+\delta_{TR1}^{1})+\right.$
(13)
$\displaystyle\quad\quad\left.\sqrt{((1-\bar{\alpha}_{1})\delta_{PS1}^{1}+\delta_{TR1}^{1})^{2}+4(1-\bar{\alpha}_{1})\delta_{TR1}^{1}\delta_{SP0}^{1}}\right]$
$\displaystyle a_{1}$
$\displaystyle:=\delta_{SP0}^{1}-\delta_{RT0}^{1}+\delta_{PS1}^{1}-\delta_{TR1}^{1}$
$\displaystyle b_{1}$ $\displaystyle:=\delta_{RT0}^{1}-\delta_{SP0}^{1}$

and

$Y:=\delta_{TR1}^{1}\delta_{SP0}^{1}-\delta_{RT0}^{1}\delta_{PS1}^{1}$ (14)

Item a) specifies the range of population 1 environmental policies where
population 2 benefits most from the maximal consumption rate
$\alpha_{2}^{*}=\theta$. Any higher consumption rate will cause the resource
to collapse. Interestingly, in this range, the resource (and consequently,
utility) is highly sensitive at the threshold, since $R(\theta_{1})>0$ and
discontinuously drops to zero for $\alpha_{2}>\theta_{1}$.

Item b) gives the range of environmental policies where population 2 benefits
most from a consumption rate that is not maximal, i.e.
$\alpha_{2}^{*}<\theta_{1}$ (12). Any higher consumption causes the resource
to degrade marginally faster, i.e. $U(\alpha_{2})$ becomes a decreasing
function for $\alpha_{2}>\alpha_{2}^{*}$. Unlike in the policies from item a),
the utility $U(\alpha_{2})$ maintains continuity on $\alpha_{2}\geq 0$, even
as it becomes zero for all $\alpha_{2}>\theta_{1}$.

The proof of Theorem 4.1 is given in Appendix -B. Numerical computations of
the optimal quantities are shown in Figure 3.

##  V Simulations: Sensitivity of cooperation incentives

From the perspective of population 2, the optimal consumption rate is chosen
to always result in a non-zero resource level, $R(\alpha_{2}^{*})$, given a
_fixed_ population 1 policy $(\delta_{SP0}^{1},\delta_{RT0}^{1})$. In this
section, we seek to understand how the choice of environmental policy
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$ impacts the resulting resource level,
under the optimal consumption rate for population 2. Intuitively, suppose a
local authority for population 1 has the option to administer incentives to
promote the low consumption strategy, given that the irresponsible population
will optimally exploit the updated policy.

The updated policy becomes $(\delta_{SP0}^{1}+u_{s},\delta_{RT0}^{1}+u_{r})$,
where $u_{s},u_{r}\geq 0$ are the added incentives. The payoff matrix
$A_{0}^{1}$ then reads as

$\begin{bmatrix}R_{0}+u_{r}&S_{0}+u_{s}\\\ T_{0}&P_{0}\end{bmatrix}$ (15)

The addition of the $u_{r}$ term represents added incentives for _mutual
cooperation_ – individuals that practice low consumption are rewarded when
many also practice the low consumption strategy. On the other hand, the
addition of the $u_{s}$ term represents added incentives for _unilateral
cooperation_ – individuals that practice low consumption are rewarded when
many practice the high consumption strategy. How should the authority select
$u_{s},u_{r}$? Which type of incentive is more effective in promoting the
resource level?

###  V-A Calculation of sensitivities

Our approach here is to evaluate the sensitivity of $R(\alpha_{2}^{*})$ with
respect to small changes in the policy. In other words, we calculate the
gradient of the (re-defined) function
$R^{*}(\delta_{SP0}^{1},\delta_{RT0}^{1}):=R(\alpha_{2}^{*})$ with respect to
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$.

The partial derivative of $R^{*}$ with respect to $\delta_{SP0}^{1}$, in the
region
$C(\delta_{SP0}^{1})\leq\delta_{RT0}^{1}\leq\frac{\delta_{TR1}^{1}}{\delta_{PS1}^{1}}\delta_{SP0}^{1}$
(part (a) of Theorem 4.1), is simply $\frac{\partial
R^{*}}{\partial\delta_{SP0}^{1}}=0$. In this case, the addition of any $u_{s}$
incentive has _no effect_ on the resource level. In the region
$\max\\{-\frac{\theta_{1}}{\alpha_{1}}\delta_{SP0}^{1},-\delta_{TR1}^{1}\\}\leq\delta_{RT0}^{1}\leq
C(\delta_{SP0}^{1})$ (part (b) of Theorem 4.1), the partial derivative with
respect to $\delta_{SP0}^{1}$ is calculated to be

$\displaystyle\frac{\partial R^{*}}{\partial\delta_{SP0}^{1}}$
$\displaystyle=\frac{\delta_{PS1}^{1}-\delta_{TR1}^{1}}{a_{1}^{2}}(1-\phi)+$
(16)
$\displaystyle\frac{-b_{1}}{2a_{1}}\phi\left(\frac{1}{-b_{1}}-\frac{\delta_{TR1}^{1}}{Y}+\frac{1-\bar{\alpha}_{1}}{-\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})}\right)$

where

$\phi:=\sqrt{\frac{Y}{b_{1}\frac{\partial g_{1}}{\partial
n}(\bar{\alpha}_{1})}},$ (17)

depends on the policy $(\delta_{SP0}^{1},\delta_{RT0}^{1})$, and $a_{1}$,
$b_{1}$, and $Y$ are variables that also depend on the policy
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$, and were defined in the statement of
Theorem 4.1. The sign of $b_{1}=\delta_{RT0}^{1}-\delta_{SP0}^{1}$ is
negative, since $C(\delta_{SP0}^{1})\leq\delta_{SP0}^{1}$ with equality if and
only if $\delta_{SP0}^{1}=0$.

The partial derivative of $R^{*}$ with respect to $\delta_{RT0}^{1}$, in the
region
$C(\delta_{SP0}^{1})\leq\delta_{RT0}^{1}\leq\frac{\delta_{TR1}^{1}}{\delta_{PS1}^{1}}\delta_{SP0}^{1}$
(part (a) of Theorem 4.1), is calculated to be

$\frac{\partial
R^{*}}{\partial\delta_{RT0}^{1}}=\frac{\delta_{TR1}^{1}}{(\delta_{RT0}^{1}+\delta_{TR1}^{1})^{2}}>0.$

In the region
$\max\\{-\frac{\theta_{1}}{\alpha_{1}}\delta_{SP0}^{1},-\delta_{TR1}^{1}\\}\leq\delta_{RT0}^{1}\leq
C(\delta_{SP0}^{1})$ (part (b) of Theorem 4.1), the partial derivative with
respect to $\delta_{RT0}^{1}$ is calculated to be

$\displaystyle\frac{\partial R^{*}}{\partial\delta_{RT0}^{1}}$
$\displaystyle=\frac{\delta_{TR1}^{1}-\delta_{PS1}^{1}}{a_{1}^{2}}(1-\phi)+$
$\displaystyle\frac{-b_{1}}{2a_{1}}\phi\left(\frac{1}{b_{1}}+\frac{\delta_{PS1}^{1}}{Y}+\frac{\bar{\alpha}_{1}}{-\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})}\right)$

The following basic property holds for both sensitivities.

######  Lemma 5.1.

For any set of fixed parameter values $\delta_{TR1}^{1}$, $\delta_{PS1}^{1}$,
$\bar{\alpha}_{1}$, and any feasible policy
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$, it holds that $\frac{\partial
R^{*}}{\partial\delta_{SP0}^{1}},\frac{\partial
R^{*}}{\partial\delta_{RT0}^{1}}\geq 0$.

Thus, increasing incentive for low consumption can never make the resource
level worse under the optimal consumption rate of population 2.

###  V-B Simulations: comparison of incentives

To address the question of which type of incentive, $u_{r}$ or $u_{s}$, is
more effective, we inspect the _sensitivity ratio_

$\rho(\delta_{SP0}^{1},\delta_{RT0}^{1}):=\frac{\partial
R^{*}/\partial\delta_{SP0}^{1}}{\partial R^{*}/\partial\delta_{RT0}^{1}}$ (18)

A ratio of $\rho<1$ indicates that the incentive $u_{r}$ is more effective
than $u_{s}$, and $\rho>1$ indicates the opposite. A numerical computation of
the sensitivity ratio is provided in Figure 4.

While a formal analysis is not yet provided in this paper, the numerical
computations strongly suggest that the $u_{r}$ incentive is generally more
effective than $u_{s}$. We observe that $\rho<1$ for a large space of policies
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$, and $\rho>1$ only when
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$ is close to the threshold curve
$(\delta_{SP0}^{1},C(\delta_{SP0}^{1}))$ and $\delta_{PS1}^{1}$ is
sufficiently small relative to $\delta_{TR1}^{1}$. For larger values of
$\delta_{PS1}^{1}$, we observe that $\rho>1$ for all policies
$(\delta_{SP0}^{1},\delta_{RT0}^{1})$. We also reiterate that $\rho=0$ for all
policies above the threshold curve $C(\delta_{SP0}^{1})$, i.e. only the
$u_{r}$ incentive can help improve the resource level.

##  VI Conclusion and Future Work

In this paper, we studied a feedback-evolving game with two populations that
share a common environmental resource. This marks initial steps in extending
the framework to multi-population interactions. We focused on a particular
scenario where one population is “responsible” about using environmental
resources, while the other population is “irresponsible”. We characterized the
asymptotic dynamic outcomes. We then evaluated to what extent the
irresponsible population can take advantage of the resources by deriving
optimal consumption rates. Lastly, a sensitivity analysis was provided that
suggests incentivizing mutual cooperation is more effective in promoting
resource levels under the two-population setting.

This work sets the stage for future studies involving hierarchical decision-
making with multiple populations. For example, how would multiple
irresponsible populations interact with one another? One can also consider
scenarios where populations have their own local environmental resources that
have externalities on others.

Figure 4: _Sensitivity ratios_. This series of plots shows the sensitivity
ratio $\rho$ as the population 1 policy $(\delta_{SP0}^{1},\delta_{RT0}^{1})$
varies. The red lines indicate the boundaries of the set of feasible policies
for population 1. In all plots, we observe the ratio is highest near the curve
$C(\delta_{SP0}^{1})$ (dashed red line), i.e. $u_{s}$ becomes effective, and
indeed can be more effective than $u_{r}$ for lower values of
$\delta_{PS1}^{1}$ (left and center plots). However, $u_{s}$ is generally not
as effective relative to $u_{r}$. For higher values of $\delta_{PS1}^{1}$,
$u_{s}$ will never be as effective as $u_{r}$, i.e. $\rho<1$ for all policies
(right plot). For all simulations, the parameter values are:
$\delta_{TR1}=10$, $\theta_{1}=0.75$, $\alpha_{1}=1$, $\epsilon=0.1$.

## References

  * [1] E. Ostrom, _Governing the commons: The evolution of institutions for collective action_. Cambridge university press, 1990.
  * [2] J. S. Weitz, C. Eksin, K. Paarporn, S. P. Brown, and W. C. Ratcliff, “An oscillating tragedy of the commons in replicator dynamics with game-environment feedback,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 47, pp. E7518–E7525, 2016.
  * [3] A. R. Tilman, J. B. Plotkin, and E. Akçay, “Evolutionary games with environmental feedbacks,” _Nature communications_ , vol. 11, no. 1, p. 915, 2020.
  * [4] L. Gong, W. Yao, J. Gao, and M. Cao, “Limit cycles analysis and control of evolutionary game dynamics with environmental feedback,” _Automatica_ , vol. 145, p. 110536, 2022.
  * [5] X. Chen and A. Szolnoki, “Punishment and inspection for governing the commons in a feedback-evolving game,” _PLoS computational biology_ , vol. 14, no. 7, p. e1006347, 2018.
  * [6] K. Paarporn, C. Eksin, J. S. Weitz, and Y. Wardi, “Optimal control policies for evolutionary dynamics with environmental feedback,” in _2018 IEEE conference on decision and control (CDC)_. IEEE, 2018, pp. 1905–1910.
  * [7] X. Wang, Z. Zheng, and F. Fu, “Steering eco-evolutionary game dynamics with manifold control,” _Proceedings of the Royal Society A_ , vol. 476, no. 2233, p. 20190643, 2020.
  * [8] A. Satapathi, N. K. Dhar, A. R. Hota, and V. Srivastava, “Coupled evolutionary behavioral and disease dynamics under reinfection risk,” _IEEE Transactions on Control of Network Systems_ , 2023.
  * [9] K. Frieswijk, L. Zino, M. Cao, and A. Morse, “Modeling the co-evolution of climate impact and population behavior: A mean-field analysis,” _IFAC-PapersOnLine_ , vol. 56, no. 2, pp. 7381–7386, 2023.
  * [10] H. Khazaei, K. Paarporn, A. Garcia, and C. Eksin, “Disease spread coupled with evolutionary social distancing dynamics can lead to growing oscillations,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 4280–4286.
  * [11] M. R. Arefin and J. Tanimoto, “Imitation and aspiration dynamics bring different evolutionary outcomes in feedback-evolving games,” _Proceedings of the Royal Society A_ , vol. 477, no. 2251, p. 20210240, 2021\.
  * [12] L. Stella and D. Bauso, “The impact of irrational behaviors in the optional prisoner’s dilemma with game-environment feedback,” _International Journal of Robust and Nonlinear Control_ , 2021.
  * [13] L. Stella, W. Baar, and D. Bauso, “Lower network degrees promote cooperation in the prisoner’s dilemma with environmental feedback,” _IEEE Control Systems Letters_ , vol. 6, pp. 2725–2730, 2022.
  * [14] Y. Kawano, L. Gong, B. D. Anderson, and M. Cao, “Evolutionary dynamics of two communities under environmental feedback,” _IEEE Control Systems Letters_ , vol. 3, no. 2, pp. 254–259, 2018.
  * [15] A. Govaert, L. Zino, and E. Tegling, “Population games on dynamic community networks,” _IEEE Control Systems Letters_ , vol. 6, pp. 2695–2700, 2022\.
  * [16] J. Certório, R. J. La, and N. C. Martins, “Epidemic population games for policy design: two populations with viral reservoir case study,” in _2023 62nd IEEE Conference on Decision and Control (CDC)_. IEEE, 2023, pp. 7667–7674.
  * [17] G. Hardin, “The tragedy of the commons,” _Science_ , vol. 162, no. 3859, pp. 1243–1248, 1968.

Throughout the Appendix, we make use of the following notations. The payoff
difference function $g_{i}(x_{i},n)$ ($i=1,2$) is bilinear in its two
arguments, which we can write as

$g_{i}(x_{i},n)=a_{i}x_{i}n+b_{i}x_{i}+c_{i}n+d_{i}$

where we define

$\displaystyle a_{i}$
$\displaystyle:=\delta_{SP0}^{i}-\delta_{RT0}^{i}+\delta_{PS1}^{i}-\delta_{TR1}^{i}$
(19) $\displaystyle b_{i}$ $\displaystyle:=\delta_{RT0}^{i}-\delta_{SP0}^{i}$
$\displaystyle c_{i}$ $\displaystyle:=-(\delta_{PS1}^{i}+\delta_{SP0}^{i})$
$\displaystyle d_{i}$ $\displaystyle:=\delta_{SP0}^{i}.$

Consequently, we can write $\frac{\partial g_{i}}{\partial
n}(x_{i})=a_{i}x_{i}+c_{i}<0$ (Assumption 2) and $\frac{\partial
g_{i}}{\partial x_{i}}(n)=a_{i}n+b_{i}$. Moreover, we will make use of the
variable

$Y:=b_{1}c_{1}-a_{1}d_{1}=\delta_{TR1}^{1}\delta_{SP0}^{1}-\delta_{RT0}^{1}\delta_{PS1}^{1}$
(20)

which is positive if and only if
$\delta_{RT0}^{1}<\frac{\delta_{TR1}}{\delta_{PS1}}\delta_{SP0}^{1}$, and zero
if it holds with equality.

###  -A Analysis of dynamical system

To establish stability properties of system (7), we need to analyze its
Jacobian matrix, which we compute below:

$J(\boldsymbol{z}):=\begin{bmatrix}J_{11}&0&J_{13}\\\ 0&J_{22}&J_{23}\\\
\epsilon J_{31}&\epsilon J_{32}&\epsilon J_{33}\\\ \end{bmatrix}$ (21)

for any $\boldsymbol{z}\in[0,1]^{3}$, where

$\displaystyle J_{11}$ $\displaystyle=x_{1}(1-x_{1})\frac{\partial
g_{1}}{\partial x_{1}}(n)+(1-2x_{1})g_{1}(x_{1},n)$ $\displaystyle J_{13}$
$\displaystyle=x_{1}(1-x_{1})\frac{\partial g_{1}}{\partial n}(x_{1})$
$\displaystyle J_{22}$ $\displaystyle=x_{2}(1-x_{2})\frac{\partial
g_{2}}{\partial x_{2}}(n)+(1-2x_{2})g_{2}(x_{2},n)$ $\displaystyle J_{23}$
$\displaystyle=x_{2}(1-x_{2})\frac{\partial g_{2}}{\partial n}(x_{2})$
$\displaystyle J_{31}$ $\displaystyle=n(1-n)(\theta_{1}+\alpha_{1})$
$\displaystyle J_{32}$ $\displaystyle=n(1-n)(\theta_{2}+\alpha_{2})$
$\displaystyle J_{33}$
$\displaystyle=(1-2n)(\theta_{1}x_{1}+\theta_{2}x_{2}-\alpha_{1}(1-x_{1})-\alpha_{2}(1-x_{2}))$

By Assumption 3, no fixed point with $x_{2}>0$ can be stable because
$g_{2}(x_{2},n)<0$ for all $x_{2},n$. Hence, we only need to analyze stability
properties of fixed points of the form $(x_{1},0,n)$.

###### Proof of Theorem 3.1.

We first observe that by Assumption 3, it must be the case that
$x_{2}(t)\rightarrow 0$ because $g_{2}(x_{2},n)<0$ for all $x_{2},n$. In other
words, for any $\varepsilon>0$ and for any initial condition in $(0,1)^{3}$,
there is some finite time $t_{\varepsilon}>0$ such that $x_{2}(t)<\varepsilon$
for all $t\geq t_{\varepsilon}$.

(1) $\underline{\alpha_{2}>\theta_{1}}$. We can analyze the limiting behavior
by considering the function $V(\boldsymbol{z}):=n$ defined for
$\boldsymbol{z}\in[0,1]\times[0,\varepsilon]\times[0,1]$. Its time derivative
is simply
$\dot{n}=n(1-n)(_{1}(\theta_{1}+\alpha_{1})-(\alpha_{1}+\alpha_{2}))$ (here,
we neglect $\varepsilon$ since it is arbitrarily small). Thus, $\dot{V}<0$ for
any $x_{1}\in(0,1)$. By LaSalle’s invariance principle, the set $n=0$ is
globally asymptotically stable.

(2) $\underline{\alpha_{2}\leq\theta_{1}}$.

The fixed point of the form $(x_{1},0,n)$ with $x_{1},n\in(0,1)$ is determined
to be $x_{1}=\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}}$ (from
$h(x_{1},n)=0$) and $n=-\frac{b_{1}x_{1}+d_{1}}{a_{1}x_{1}+c_{1}}$ (from
$g_{1}(x_{1},n)=0$). It holds that $n\in(0,1)$ if and only if
$\frac{\alpha_{2}-\theta_{1}}{\alpha_{1}+\alpha_{2}}\delta_{SP0}^{1}\leq\delta_{RT0}^{1}$.
The eigenvalues of $J$ at this point are given by

$\lambda_{1}=-(\delta_{PS1}^{1}+\delta_{SP0}^{1})n+\delta_{SP0}^{2}<0$

and

$\displaystyle\lambda_{2,3}$
$\displaystyle=\frac{1}{2}\left[\frac{}{}(a_{1}n+b_{1})x_{1}(1-x_{1})\pm\right.$
$\displaystyle\left.\sqrt{(a_{1}n+b_{1})^{2}x_{1}^{2}(1-x_{1})^{2}-4\epsilon
K}\right]$

where we denote

$K=n(1-n)(\theta_{1}+\alpha_{1})x_{1}(1-x_{1})|\frac{\partial g_{1}}{\partial
n}(x_{1})|>0.$

The sign of the real parts of $\lambda_{2,3}$ is thus given by the sign of
$a_{1}n+b_{1}=\frac{-Y}{|\frac{\partial g_{1}}{\partial n}(x_{1})|}<0$.

The fixed point of the form $(x_{1},0,0)$ with $x_{1}\in(0,1)$ is determined
to be $x_{1}=\frac{\delta_{SP0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}$ (from
$g_{1}(x_{1},0)=0$). This implies $\delta_{RT0}^{1}<0$. The eigenvalues of $J$
at this point are
$\\{x_{1}(1-x_{1})b_{1},\delta_{SP0}^{2},\epsilon(x_{1}(\theta_{1}+\alpha_{1})-(\alpha_{1}+\alpha_{2}))\\}$.
The first and second are negative, and the third is also negative if and only
if
$\frac{\alpha_{2}-\theta_{1}}{\alpha_{1}+\alpha_{2}}\delta_{SP0}^{1}>\delta_{RT0}^{1}$.

We now verify that all other fixed points cannot be attractive. To see why, we
can split the $(x_{1},n)$ plane into four disjoint sectors:

$\displaystyle Q_{1}$
$\displaystyle:=\\{(x_{1},n):x_{1}>\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},n<\frac{b_{1}x_{1}+d_{1}}{-(a_{1}x_{1}+c_{1})}\\}$
(22) $\displaystyle Q_{2}$
$\displaystyle:=\\{(x_{1},n):x_{1}>\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},n\geq\frac{b_{1}x_{1}+d_{1}}{-(a_{1}x_{1}+c_{1})}\\}$
$\displaystyle Q_{3}$
$\displaystyle:=\\{(x_{1},n):x_{1}\leq\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},n\geq\frac{b_{1}x_{1}+d_{1}}{-(a_{1}x_{1}+c_{1})}\\}$
$\displaystyle Q_{4}$
$\displaystyle:=\\{(x_{1},n):x_{1}\leq\frac{\alpha_{1}+\alpha_{2}}{\alpha_{1}+\theta_{1}},n<\frac{b_{1}x_{1}+d_{1}}{-(a_{1}x_{1}+c_{1})}\\}$

The vector field of the dynamics satisfy $\dot{x}_{1},\dot{n}>0$ in $Q_{1}$,
$\dot{x}_{1}\leq 0,\dot{n}>0$ in $Q_{2}$, $\dot{x}_{1}\leq 0,\dot{n}\leq 0$ in
$Q_{3}$, and $\dot{x}_{1}>0,\dot{n}\leq 0$ in $Q_{4}$. The equalities hold
only on the shared boundaries between the sectors. The unique interior fixed
point lies at the single point of intersection of all four sectors. Moreover,
we also observe that no point on the borders of $[0,1]^{2}$ can be attractive,
since the vector field is oriented such that none of the four quadrants are
positively invariant. This establishes the claim. ∎

###  -B Analysis of optimal consumption

The proof of Theorem 4.1 relies on establishing certain properties about the
utility function $U(\alpha_{2})$. Let us define the _support_ of $U$ as

$\text{supp }U:=\\{\alpha_{2}\geq 0:U(\alpha_{2})>0\\}.$ (23)

Throughout the proof, we will use the shorthand notations

$\bar{\alpha}_{1}:=\frac{\alpha_{1}}{\alpha_{1}+\theta_{1}},\quad\bar{\alpha}_{2}:=\frac{\alpha_{2}}{\alpha_{1}+\theta_{1}}$
(24)

###### Proof of Theorem 4.1.

First, we notice that $U(\alpha_{2})=0$ for all $\alpha_{2}>\theta_{1}$, so
$U(\alpha_{2})>0$ only if $\alpha_{2}\leq\theta_{1}$. The support of $U$ is
the interval $[0,\theta_{1}]$ in the case that $\delta_{RT0}^{1}>0$, and is
the interval
$[0,\frac{\delta_{SP0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\theta_{1}+\frac{\delta_{RT0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\alpha_{1}]$
in the case that
$\frac{-\theta_{1}}{\alpha_{1}}\delta_{SP0}^{1}\leq\delta_{RT0}^{1}\leq 0$. At
any point in the support, the first derivative of $U$ is

$\displaystyle U^{\prime}(\alpha_{2})$
$\displaystyle=R(\alpha_{2})+\alpha_{2}\frac{\partial R}{\partial\alpha_{2}}$
(25)
$\displaystyle=R(\alpha_{2})-\alpha_{2}\frac{Y}{(\alpha_{1}+\theta_{1})(\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1}+\bar{\alpha}_{2}))^{2}}$
$\displaystyle=\frac{b_{1}(\bar{\alpha}_{1}+\bar{\alpha}_{2})+d_{1}}{-(a_{1}(\bar{\alpha}_{1}+\bar{\alpha}_{2})+c_{1})}-\bar{\alpha}_{2}\frac{Y}{(a_{1}(\bar{\alpha}_{1}+\bar{\alpha}_{2})+c_{1})^{2}}$

We note that $Y\geq 0$ is equivalent to
$\delta_{RT0}^{1}\leq\frac{\delta_{TR1}^{1}}{\delta_{PS1}^{1}}\delta_{SP0}^{1}$.
Evaluating the second derivative, we obtain

$\displaystyle U^{\prime\prime}(\alpha_{2})$
$\displaystyle=-2\frac{Y}{(\alpha_{1}+\theta_{1})(\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1}+\bar{\alpha}_{2}))^{2}}$
$\displaystyle\quad+\frac{\alpha_{2}Y}{(\alpha_{1}+\theta_{1})}\frac{(-2a_{1})}{(\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1}+\bar{\alpha}_{2}))^{3}}$
$\displaystyle\propto-\left(1+\frac{\alpha_{2}}{\alpha_{1}+\theta_{1}}\frac{a_{1}}{|\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1}+\bar{\alpha}_{2})|}\right)$

The sign of $U^{\prime\prime}$ then depends only on the factor in parentheses.
This factor is monotonic in $\alpha_{2}\in\text{supp }U$. Indeed, the sign of
its derivative is the sign of ${-a_{1}^{2}\bar{\alpha}_{1}-a_{1}c_{1}}$, which
is independent of $\alpha_{2}$. If it is positive, then $U^{\prime\prime}<0\
\forall\alpha_{2}\in\text{supp }U$. If it is negative, we need to verify that
$1+(1-\bar{\alpha}_{1})\frac{a_{1}}{|\frac{\partial g_{1}}{\partial
n}(1)|}>0$. Indeed, this is equivalent to the inequality

$(1-\bar{\alpha}_{1})\frac{\delta_{SP0}^{1}+\delta_{PS1}}{\delta_{RT0}^{1}+\delta_{TR1}}>-\bar{\alpha}_{1}$
(26)

This is always satisfied because the LHS is positive due to Assumption 2
($\delta_{RT0}^{1}>-\delta_{TR1}$), and the RHS is negative. This establishes
concavity of $U$ in its interval of support.

(a) In this sub-regime, $\text{supp }U=[0,\theta_{1}]$. We need to establish
that $U$ is strictly increasing. Since we know that $U^{\prime}(0)=R(0)>0$ and
$U$ is concave, we just need to verify $U^{\prime}(\theta_{1})\geq 0$. We have

$\displaystyle U^{\prime}(\theta_{1})$
$\displaystyle=R(\theta_{1})-(1-\bar{\alpha}_{1})\frac{Y}{\frac{\partial
g_{1}}{\partial n}(1)^{2}}$
$\displaystyle=\frac{b_{1}+d_{1}}{-(a_{1}+c_{1})}-(1-\bar{\alpha}_{1})\frac{Y}{(a_{1}+c_{1})^{2}}$

Multiplying by $(a_{1}+c_{1})^{2}$, we have $U^{\prime}(\theta_{1})\geq 0$ if
and only if

$\displaystyle(\delta_{RT0}^{1})^{2}+((1-\bar{\alpha}_{1})\delta_{PS1}+\delta_{TR1})\delta_{RT0}^{1}$
$\displaystyle\quad-(1-\bar{\alpha}_{1})\delta_{TR1}\delta_{SP0}^{1}\geq 0$

Observe the LHS above is quadratic in $\delta_{RT0}^{1}$, and has a negative
and a positive root. Since we are considering $\delta_{RT0}^{1}>0$ in this
regime, the condition above is satisfied if and only if $\delta_{RT0}^{1}$ is
greater than the positive root, i.e.

$\delta_{RT0}^{1}\geq C(\delta_{SP0}^{1})$

(b) Case 1: First we consider $0<\delta_{RT0}^{1}<C(\delta_{SP0}^{1})$. This
implies $\text{supp }U=[0,\theta_{1}]$ and $U^{\prime}(\theta_{1})<0$.
Therefore, $U(\alpha_{2})$ attains its maximum for some
$\alpha_{2}^{*}\in(0,\theta_{1})$. To find $\alpha_{2}^{*}$, we need to solve
$U^{\prime}(\alpha_{2})=0$, which yields the quadratic equation

$\displaystyle-a_{1}b_{1}\bar{\alpha}_{2}^{2}-2b_{1}\frac{\partial
g_{1}}{\partial
n}(\bar{\alpha}_{1})\bar{\alpha}_{2}+g_{1}(\bar{\alpha}_{1})\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})=0$ (27)

The “minus” solution is given by

$\displaystyle\alpha_{2}^{*}$
$\displaystyle=\frac{\theta_{1}+\alpha_{1}}{-a_{1}b_{1}}\left[b_{1}\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})-\right.$
$\displaystyle\quad\left.\sqrt{b_{1}^{2}\frac{\partial g_{1}}{\partial
n}(\bar{\alpha}_{1})^{2}-a_{1}b_{1}g_{1}(\bar{\alpha}_{1},0)\frac{\partial
g_{1}}{\partial n}(\bar{\alpha}_{1})}\right]$

which reduces to $\alpha_{2}^{*}$ in the Lemma statement. We deduce the
“minus” solution must correspond to the critical point by showing it is
positive. The sign of $1-\sqrt{\frac{Y}{b_{1}\frac{\partial g_{1}}{\partial
n}(\bar{\alpha}_{1})}}$ coincides with the sign of $a_{1}$:

$1-\sqrt{\frac{Y}{b_{1}\frac{\partial g_{1}}{\partial
n}(\bar{\alpha}_{1})}}<0\iff a_{1}(b_{1}\bar{\alpha}_{1}+d_{1})<0$

Observe $b_{1}\bar{\alpha}_{1}+d_{1}=g_{1}(\bar{\alpha}_{1},0)>0$ since
$\delta_{SP0}^{1},\delta_{RT0}^{1}>0$, and therefore $\alpha_{2}^{*}$ is
positive. On the other hand, the “plus” solution is negative since the sign of
$1+\sqrt{\frac{Y}{b_{1}\frac{\partial g_{1}}{\partial n}(\bar{\alpha}_{1})}}$
is the opposite of $a_{1}$.

Case 2: For
$\frac{-\theta_{1}}{\alpha_{1}}\delta_{SP0}^{1}\leq\delta_{RT0}^{1}\leq 0$,
$\text{supp
}U=[0,\frac{\delta_{SP0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\theta_{1}+\frac{\delta_{RT0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\alpha_{1}]$.
We can express the utility as

$U(\alpha_{2})=\alpha_{2}\cdot\begin{cases}-\frac{(\delta_{RT0}^{1}-\delta_{SP0}^{1})({\bar{\alpha}}_{1}+{\bar{\alpha}}_{2})+\delta_{SP0}^{1}}{a_{1}({\bar{\alpha}}_{1}+{\bar{\alpha}}_{2})+c_{1}},&\text{if
}\alpha_{2}\in\text{supp }U\\\ 0,&\text{else }\end{cases}$

It holds that $U^{\prime}(0)>0$,
$U^{\prime}(\frac{\delta_{SP0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\theta_{1}+\frac{\delta_{RT0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\alpha_{1})<0$.
Since $U$ is concave, it must attain its maximum in the interval
$\alpha_{2}\in(0,\frac{\delta_{SP0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\theta_{1}+\frac{\delta_{RT0}^{1}}{\delta_{SP0}^{1}-\delta_{RT0}^{1}}\alpha_{1})$.
The expression for this point is identical to $\alpha_{2}^{*}$ from Case 1 due
to similar arguments. ∎


# Achievable Rate Maximization for Underlay Spectrum Sharing MIMO System with
Intelligent Reflecting Surface

Vaibhav Kumar, Member, IEEE, Mark F. Flanagan, Senior Member, IEEE, Rui Zhang,
Fellow, IEEE,
and Le-Nam Tran, Senior Member, IEEE This publication has emanated from
research conducted with the financial support of Science Foundation Ireland
(SFI) and is co-funded under the European Regional Development Fund under
Grant Number 17/CDA/4786.
V. Kumar, M. F. Flanagan, and L.-N. Tran are with School of Electrical and
Electronic Engineering, University College Dublin, D04 V1W8 Dublin, Ireland
(email<EMAIL_ADDRESS><EMAIL_ADDRESS>nam.tran@ucd.ie).
R. Zhang is with the Department of Electrical and Computer Engineering,
National University of Singapore, Singapore (e-mail: elezhang@nus.edu.sg).

###### Abstract

In this letter, the achievable rate maximization problem is considered for
intelligent reflecting surface (IRS) assisted multiple-input multiple-output
(MIMO) systems in an underlay spectrum sharing scenario, subject to
interference power constraints at the primary users. The formulated non-convex
optimization problem is challenging to solve due to its non-convexity as well
as coupling design variables in the constraints. Different from existing works
that are mostly based on alternating optimization (AO), we propose a penalty
dual decomposition based gradient projection (PDDGP) algorithm to solve this
problem. We also provide a convergence proof and a complexity analysis for the
proposed algorithm. We benchmark the proposed algorithm against two known
solutions, namely a minimum mean-square error based AO algorithm and an inner
approximation method with block coordinate descent. Specifically, the
complexity of the proposed algorithm grows linearly with respect to the number
of reflecting elements at the IRS, while that of the two benchmark methods
grows with the third power of the number of IRS elements. Moreover, numerical
results show that the proposed PDDGP algorithm yields considerably higher
achievable rate than the benchmark solutions.

###### Index Terms:

Intelligent reflecting surface, underlay spectrum sharing, MIMO, penalty dual
decomposition, gradient projection.

##  I Introduction

With the recent development of software-controlled metasurface technology,
intelligent reflecting surfaces (IRSs) are being regarded as a promising
technology to enhance the performance of next-generation wireless
communications [1]. By virtue of flexible configuration of the radio
propagation environment, IRSs are capable of increasing the spectral as well
as the energy efficiency of wireless communication systems. On the other hand,
spectrum sharing has great potential to support massive connectivity in the
congested radio spectrum for beyond-fifth-generation (B5G) communications
standards [2, 3]. Therefore, the design and analysis of IRS-assisted spectrum
sharing to enhance the achievable rate of wireless communication systems have
gained considerable attention in the recent years.

In this context, the problem of optimal active and passive beamforming for
secondary user (SU) rate maximization in an IRS-assisted single-input single-
output (SISO) and multiple-input single-output (MISO) spectrum sharing systems
were considered in [4, 5, 6, 7, 8]. However, studies on SU rate maximization
for the more challenging IRS-assisted multiple-input multiple-output (MIMO)
underlay spectrum sharing systems is still limited. In order to maximize the
weighted sum rate of the secondary network in an IRS-assisted MIMO system in
an underlay spectrum sharing scenario with multiple secondary receivers (SRs)
and a single primary receiver (PR), a block coordinate descent algorithm along
with an inner approximation method (BCD-IAM) was proposed in [9]. Similarly,
in [10], the authors considered the problem of weighted sum rate maximization
of the secondary network for an IRS-assisted MIMO spectrum sharing system
consisting of multiple SRs as well as multiple PRs, where they adopted a
weighted minimum mean-square error (WMMSE) method and an alternating
optimization (AO)-based algorithm to find a suboptimal solution. On the other
hand, a suboptimal solution to the problem of achievable secrecy rate
maximization for an IRS-assisted MIMO system using AO, barrier method, SDR and
minorization-maximization (MM), was presented in [11].

It is noteworthy that the methods proposed to maximize the SU rate for MISO
systems in [4, 5, 6, 7, 8] are not applicable to MIMO systems, unless each
receive antenna in the MIMO system is treated as a separate user. However,
this reduces the beamforming/multiplexing gains which are major benefits of
MIMO systems. The method presented in [9] for MIMO systems is proposed only
for the case of a single PR. Similarly, the system considered in [11], when
omitting the eavesdropper, reduces to the conventional IRS-assisted MIMO
spectrum sharing system. However, the solution proposed in [11] is limited to
the scenario where the direct links between the (secondary) transmitter and
(primary and secondary) receivers are absent. On the other hand, the AO-based
solution proposed in [10] can be inefficient due to the possibly high coupling
of design variables, especially when the SU rate performance is dominated by
the interference constraints at the PRs. Motivated by the drawbacks of
existing methods, in this letter we propose a low-complexity but efficient
algorithm to maximize the achievable rate of the SU in an IRS-assisted MIMO
underlay spectrum sharing system consisting of multiple PRs, where direct
links between the (secondary) transmitter and (primary and secondary)
receivers are also present. Our main contributions in this letter are as
follows: (i) we propose a penalty dual decomposition based gradient projection
(PDDGP) algorithm to maximize the achievable rate; (ii) we provide a
convergence proof and complexity analysis of the proposed algorithm; and (iii)
we provide extensive numerical results to demonstrate the superiority of the
proposed solution over the benchmark solutions provided in [9] and [10], and
to show the effect of different system parameters on the achievable rate.

#### Notations

Bold lowercase and uppercase letters denote vectors and matrices,
respectively. $\mathbf{X}^{*}$, $\mathbf{X}^{\mathsf{T}}$,
$\mathbf{X}^{{\dagger}}$, $\left|\mathbf{X}\right|$,
$\left\|\mathbf{X}\right\|$ and $\operatorname{Tr}(\mathbf{X})$ denote the
complex conjugate, (ordinary) transpose, Hermitian transpose, determinant,
Frobenius norm and trace of a matrix $\mathbf{X}$, respectively. The vector
space of complex matrices of size $M\times N$ is denoted by
$\mathbb{C}^{M\times N}$. $\operatorname{diag}(\mathbf{x})$ denotes a square
diagonal matrix where the elements of $\mathbf{x}$ comprise the main diagonal,
$\operatorname{\mathrm{vec_{d}}}(\mathbf{X})$ denotes the column vector
containing the elements of the main diagonal of the matrix $\mathbf{X}$, and
$\nabla_{\mathbf{X}}f(\cdot)$ denotes the complex-valued gradient of
$f(\cdot)$ with respect to (w.r.t.) $\mathbf{X}^{*}$. $\mathbb{E}\\{\cdot\\}$
denotes the expectation operator. $\mathbf{X}\succeq\mathbf{Y}$ represents
that $\mathbf{X}-\mathbf{Y}$ is positive semidefinite, and
$\Pi_{\mathcal{X}}(\mathbf{x})$ denotes the Euclidean projection of the point
$\mathbf{x}$ onto the set $\mathcal{X}$, i.e.,
$\Pi_{\mathcal{X}}(\mathbf{x})\triangleq\operatorname{argmin}_{\hat{\mathbf{x}}\in\mathcal{X}}\left\|\hat{\mathbf{x}}-\mathbf{x}\right\|$.

Figure 1: System model for IRS-assisted MIMO underlay spectrum sharing.

1

2

Input: $\boldsymbol{\theta}_{0}$, $\mathbf{X}_{0}$, $\mathbf{s}_{0}$,
$\boldsymbol{\upsilon}$, $\mu_{0}$, $\alpha_{0}$, $\rho$

3

Output: $\boldsymbol{\theta}_{n},\mathbf{X}_{n}$

4

5$n\leftarrow 1$

6repeat

7

8
$\\!\\!\boldsymbol{\theta}_{n}=\Pi_{\vartheta}(\boldsymbol{\theta}_{n\\!-\\!1}\\!+\\!\mu_{n}\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n\\!-\\!1},\boldsymbol{\theta}_{n\\!-\\!1},\mathbf{s}_{n\\!-\\!1}))$;

9

10
$\\!\\!\mathbf{X}_{n}=\Pi_{\mathcal{X}}(\mathbf{X}_{n\\!-\\!1}+\alpha_{n}\nabla_{\mathbf{X}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n\\!-\\!1},\boldsymbol{\theta}_{n},\mathbf{s}_{n\\!-\\!1}))$;

11

12
$\\!\\!s_{n,k}=\max\\{0,P_{k}\\!-\\!\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}_{n}\mathbf{Z}_{k}^{{\dagger}})\\}$,$\forall
k\in\mathcal{K}$

13 where
$\mathbf{Z}_{k}=\mathbf{H}_{\mathrm{I}k}\boldsymbol{\Theta}_{n}\mathbf{H}_{\mathrm{TI}}+\mathbf{H}_{\mathrm{T}k}$;

14

15 $n\leftarrow n+1$;

16

17until  _convergence_ ;

Algorithm 1 Gradient projection algorithm to solve (3) for fixed
$\boldsymbol{\upsilon}$ and $\rho$.

##  II System Model and Problem Formulation

Let us first consider an IRS-assisted MIMO system in an underlay spectrum
sharing scenario shown in Fig. 1, consisting of one secondary transmitter
(ST), one SR, and $K$ PRs. The number of antennas at the ST and SR are denoted
by $N_{\mathrm{T}}$ and $N_{\mathrm{R}}$, respectively, while the $K$ PRs are
each equipped with $N_{\mathrm{P}}$ antennas. The communication between the ST
and SR is assisted by an IRS consisting of $N_{\mathrm{I}}$ low-cost passive
reflecting elements. The ST-SR, ST-IRS, and IRS-SR channels are denoted by
$\mathbf{H}_{\mathrm{TR}}\in\mathbb{C}^{N_{\mathrm{R}}\times N_{\mathrm{T}}}$,
$\mathbf{H}_{\mathrm{TI}}\in\mathbb{C}^{N_{\mathrm{I}}\times N_{\mathrm{T}}}$,
and $\mathbf{H}_{\mathrm{IR}}\in\mathbb{C}^{N_{\mathrm{R}}\times
N_{\mathrm{I}}}$, respectively. On the other hand, the channel between the ST
and the $k$-th PR ($k\in\mathcal{K}\triangleq\left\\{1,2,\ldots,K\right\\}$)
is denoted by $\mathbf{H}_{\mathrm{T}k}\in\mathbb{C}^{N_{\mathrm{P}}\times
N_{\mathrm{T}}}$, while the channel between the IRS and the $k$-th PR is
denoted by $\mathbf{H}_{\mathrm{I}k}\in\mathbb{C}^{N_{\mathrm{P}}\times
N_{\mathrm{I}}}$. Let $\boldsymbol{\theta}\triangleq[\theta_{1}\ \theta_{2}\
\ldots\
\theta_{N_{\mathrm{I}}}]^{\mathsf{T}}\in\mathbb{C}^{N_{\mathrm{I}}\times 1}$
denote the phase-shift vector at the IRS, where $\theta_{l}=e^{j\phi_{l}}$,
$l\in\mathcal{N}_{\mathrm{I}}\triangleq\left\\{1,2,\ldots,N_{\mathrm{I}}\right\\}$,
and $\phi_{l}\in[0,2\pi)$. Here the ST, SR and IRS constitute the secondary
network, and in order to communicate with the SR, the ST uses the spectrum
band whose license is owned by a primary network. However, with the help of
properly designed beamformers at the ST and IRS, it is ensured that the
interference inflicted by the secondary network on the PRs is below a
predefined threshold.111Following the arguments in the related works of [6, 7,
9, 11, 12], we assume that the secondary network is deployed outside the
coverage of the primary transmitter, and that the distance between the primary
transmitter and secondary receiver is large enough for the primary-to-
secondary interference to be negligible. However, it is noteworthy that the
proposed PDDGP algorithm is also applicable when the secondary network is
deployed completely or partially inside the coverage of the primary
transmitter.

The signal received at the SR is given by

$\displaystyle\mathbf{y}_{\mathrm{R}}=(\mathrm{\mathbf{H}_{\mathrm{IR}}}\boldsymbol{\Theta}\mathbf{H}_{\mathrm{TI}}+\mathbf{H}_{\mathrm{TR}})\mathbf{x}+\mathbf{n}_{\mathrm{R}},$
(1)

where $\boldsymbol{\Theta}\triangleq\operatorname{diag}(\boldsymbol{\theta})$,
$\mathbf{x}\in\mathbb{C}^{N_{\mathrm{T}}\times 1}$ is the transmitted signal
vector from the ST (which is intended for the SR),
$\mathbf{n}_{\mathrm{R}}\sim\mathcal{CN}(\boldsymbol{0},\sigma^{2}\mathbf{I})$
denotes the complex additive white Gaussian noise (AWGN) vector at the SR, and
$\sigma^{2}$ is the noise power. Without loss of generality, we assume that
the noise at all of the PRs is also distributed as
$\mathcal{CN}(\boldsymbol{0},\sigma^{2}\mathbf{I})$. With a slight abuse of
notation, in the rest of letter, we make the replacements
$\mathbf{H}_{\mathrm{TI}}\leftarrow\mathbf{H}_{\mathrm{TI}}/\sigma$,
$\mathbf{H}_{\mathrm{TR}}\leftarrow\mathbf{H}_{\mathrm{TR}}/\sigma$ and
$\mathbf{H}_{\mathrm{T}k}\leftarrow\mathbf{H}_{\mathrm{T}k}/\sigma$. This
normalization step will also mitigate potential numerical issues since we
avoid dealing with extremely small quantities.

#### Problem Formulation

Denote the transmit covariance matrix at the ST by
$\mathbf{X}\triangleq\mathbb{E}\\{\mathbf{x}\mathbf{x}^{{\dagger}}\\}$, the
maximum transmit power budget at the ST by $P_{\max}$, and the maximum
tolerable interference power (normalized to the noise power) at the $k$-th PR
by $P_{k}$. Then, the achievable rate maximization problem for the secondary
network can be formulated as follows:

$\displaystyle\underset{\mathbf{\mathbf{X}\succeq}\boldsymbol{0},\boldsymbol{\theta}}{\operatorname{maximize}}\
$ $\displaystyle
R(\mathbf{X},\boldsymbol{\theta})=\ln|\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}}|,$
(2a) $\displaystyle\operatorname{subject~{}to}\ $
$\displaystyle\operatorname{Tr}(\mathbf{X})\leq P_{\max},$ (2b)
$\displaystyle\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}\mathbf{Z}_{k}^{{\dagger}})\leq
P_{k},\forall k\in\mathcal{K},$ (2c) $\displaystyle|\theta_{l}|=1,\forall
l\in\mathcal{N}_{\mathrm{I}},$ (2d)

where
$\mathbf{Z}_{\mathrm{R}}\triangleq\mathrm{\mathbf{H}_{\mathrm{IR}}}\boldsymbol{\Theta}\mathbf{H}_{\mathrm{TI}}+\mathbf{H}_{\mathrm{TR}}$,
and
$\mathbf{Z}_{k}\triangleq\mathbf{H}_{\mathrm{I}k}\boldsymbol{\Theta}\mathbf{H}_{\mathrm{TI}}+\mathbf{H}_{\mathrm{T}k}$.
Note that the constraints in (2b)–(2d) represent the transmit power constraint
at the ST, the interference power constraints (IPCs) at the PRs and the unit-
modulus constraints (UMCs) at the IRS, respectively. Furthermore, the feasible
sets for $\mathbf{X}$ and $\boldsymbol{\theta}$ are defined as
$\mathcal{X}\triangleq\\{\mathbf{X}\in\mathbb{C}^{N_{\mathrm{T}}\times
N_{\mathrm{T}}}:\mathrm{Tr}(\mathbf{X})\leq
P_{\max};\mathbf{X}\succeq\boldsymbol{0}\\}$ and
$\vartheta\triangleq\\{\boldsymbol{\theta}\in\mathbb{C}^{N_{\mathrm{I}}\times
1}:\left|\theta_{l}\right|=1,l\in\mathcal{N}_{\mathrm{I}}\\}$, respectively.

##  III Proposed Solution

###  III-A PDDGP Algorithm

It is easy to see that problem (2) is non-convex due to the non-convexity of
the objective as well as constraints (2c) and (2d). If the IPCs were not
present, (2) would reduce to the conventional MIMO capacity problem for which
the gradient projection algorithm in [13] can be applied to obtain a
stationary solution. In fact, it is the coupling of $\mathbf{X}$ and
$\boldsymbol{\theta}$ due to the IPCs in (2c) that makes the problem much more
difficult to solve and the gradient projection algorithm proposed in [13] is
no longer applicable. The main idea of the existing optimization methods in
[10, 9] for solving (2) is to alternately optimize $\mathbf{X}$ and
$\boldsymbol{\theta}$. However, due to the coupling of $\mathbf{X}$ and
$\boldsymbol{\theta}$, these methods are not guaranteed to return a stationary
solution. Indeed we numerically demonstrate in Section IV that these methods
are inferior to the PDDGP algorithm that we propose in the following.

Specifically, in order to cope with the IPCs, we apply the penalty dual
decomposition (PDD) method recently presented in [14]. First, let
$g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})\triangleq\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}\mathbf{Z}_{k}^{{\dagger}})+s_{k}-P_{k}$.
Then it is easy to see that (2c) is equivalent to
$g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})=0,s_{k}\geq 0,\forall
k\in\mathcal{K}$. Next, the augmented Lagrangian function can be written as
follows:

$\displaystyle\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})\triangleq\ln|\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}}|$
$\displaystyle\qquad-\sum_{k\in\mathcal{K}}\upsilon_{k}g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})-\frac{1}{2\rho}\sum_{k\in\mathcal{K}}g_{k}^{2}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k}),$

where $\mathbf{s}\triangleq[s_{1},s_{2},\ldots,s_{K}]^{\mathsf{T}}$,
$\boldsymbol{\upsilon}\triangleq[\upsilon_{1},\upsilon_{2},\ldots,\upsilon_{K}]^{\mathsf{T}}$,
$\upsilon_{k}$ is the Lagrangian multiplier corresponding to the constraint
$g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})=0$ and $\rho$ is the
_penalty parameter_. Hence, for a given $(\boldsymbol{\upsilon},\rho)$ an
equivalent optimization problem can be formulated as

$\displaystyle\underset{\mathbf{X},\boldsymbol{\theta},\mathbf{s}}{\operatorname{maximize}}$
$\displaystyle\
\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s}),$
(3a) $\displaystyle\operatorname{subject~{}to}$ $\displaystyle\
\mathbf{s}\geq\boldsymbol{0},\eqref{eq:SPC},\eqref{eq:UMC}.$ (3b)

Since the IPCs (and thus the coupling of $\mathbf{X}$ and
$\boldsymbol{\theta}$) is now included in the augmented objective, the
optimization variables in (3) are decoupled. As a result, a simple but
efficient method can be applied to find a stationary solution to (3) which is
based on the alternating projected gradient method (APGM). The APGM is
motivated by the fact that the feasible set with respect to individual
variables is simple in the sense that the projection onto this set can be
expressed in closed form.

First, a projected gradient step for $\boldsymbol{\theta}$ is performed while
other variables are held fixed. In this regard, the gradient of
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
w.r.t. $\boldsymbol{\theta}$ is given by

$\displaystyle\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})=\nabla_{\boldsymbol{\theta}}\ln|\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}}|$
$\displaystyle-\sum_{k\in\mathcal{K}}\big{[}\upsilon_{k}+\frac{1}{\rho}g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})\big{]}\nabla_{\boldsymbol{\theta}}g_{k}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k}).$
(4)

Next using (4), [13, eqn. (17a)] and [15, Table 4.3 and eqn. (6.153)],
$\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
is given by (5) shown on the next page. Suppose that, from the current point
$\boldsymbol{\theta}$, after moving along
$\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
with some step size one arrives at a point
$\hat{\boldsymbol{\theta}}=[\hat{\theta}_{1},\hat{\theta}_{2},\ldots,\hat{\theta}_{N_{\mathrm{I}}}]^{\mathsf{T}}\in\mathbb{C}^{N_{\mathrm{I}}\times
1}$. To obtain the next iterate, the projection of $\hat{\boldsymbol{\theta}}$
onto the set $\vartheta$ is needed.

$\displaystyle\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})\\!=\\!\operatorname{\mathrm{vec_{d}}}\big{(}\mathbf{H}_{\mathrm{IR}}^{{\dagger}}\big{(}\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}}\big{)}^{-1}\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{H}_{\mathrm{TI}}^{{\dagger}}\big{)}\\!-\\!\sum_{k\in\mathcal{K}}\big{[}\upsilon_{k}\\!+\\!\tfrac{1}{\rho}\\{\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}\mathbf{Z}_{k}^{{\dagger}})\\!+\\!s_{k}\\!-\\!P_{k}\\}\big{]}\operatorname{\mathrm{vec_{d}}}(\mathbf{H}_{\mathrm{I}k}^{{\dagger}}\mathbf{Z}_{k}\mathbf{X}\mathbf{H}_{\mathrm{TI}}^{{\dagger}}).$
(5)

It is easy to see that this projection is given by
$\Pi_{\vartheta}(\hat{\boldsymbol{\theta}})=[\bar{\theta}_{1},\bar{\theta}_{2},\ldots,\bar{\theta}_{N_{\mathrm{I}}}]^{\mathsf{T}}$,
where for each $l\in\mathcal{N}_{\mathrm{I}}$

$\bar{\theta}_{l}=\left\\{\begin{array}[]{cc}\hat{\theta}_{l}/\big{|}\hat{\theta}_{l}\big{|},&\text{if
}\hat{\theta}_{l}\neq 0\\\
\exp(j\phi),\phi\in[0,2\pi),&\text{otherwise}.\end{array}\right.$ (6)

Next, a projected gradient step for $\mathbf{X}$ is carried out. Following a
similar line of argument and using [15, eqns. (6.207), (6.195) and Table 4.3],
a closed-form expression for the gradient of
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
w.r.t. $\mathbf{X}$ can be obtained in (7), shown on the next page.

$\nabla_{\mathbf{X}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})=\mathbf{Z}_{\mathrm{R}}^{{\dagger}}(\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}})^{-1}\mathbf{Z}_{\mathrm{R}}-\sum_{k\in\mathcal{K}}\big{[}\upsilon_{k}+\tfrac{1}{\rho}\\{\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}\mathbf{Z}_{k}^{{\dagger}})+s_{k}-P_{k}\\}\big{]}\mathbf{Z}_{k}^{{\dagger}}\mathbf{Z}_{k}.$
(7)

Moreover, the projection of a given point
$\mathbf{W}\in\mathbb{C}^{N_{\mathrm{T}}\times N_{\mathrm{T}}}$ onto the
feasible set $\mathcal{X}$ can be shown to admit a water-filling solution. Due
to space constraints, the details regarding the projection of $\mathbf{W}$
onto $\mathcal{X}$ are omitted, while the interested reader may refer to [13,
Sec. III-C] for details.

1

2

Input: $\boldsymbol{\theta}_{0}$, $\mathbf{X}_{0}$, $\mathbf{s}_{0}$,
$\boldsymbol{\upsilon}_{0}$, $\mu_{0}$, $\alpha_{0}$, $\rho$, $\kappa<1$

3

Output: $\boldsymbol{\theta}^{\star},\mathbf{X}^{\star}$

4

5repeat

6

7 Solve problem (3) using Algorithm 1;

8

9 $\mathbf{X}^{\star}\leftarrow\mathbf{X}_{n}$,
$\boldsymbol{\theta}^{\star}\leftarrow\boldsymbol{\theta}_{n}$,
$\mathbf{s}^{\star}\leftarrow\mathbf{s}_{n}$;

10

11
$\upsilon_{k}\leftarrow\upsilon_{k}+\frac{1}{\rho}g_{k}(\mathbf{X}^{\star},\boldsymbol{\theta}^{\star},P_{k},s_{k}^{\star})$,
$\forall k\in\mathcal{K}$;

12

13 $\rho\leftarrow\kappa\rho$;

14

15until  _convergence_ ;

Algorithm 2 The PDDGP algorithm.

We now turn to the optimization over $\mathbf{s}$. It is easy to see that the
optimal solution to (3), when $\boldsymbol{\theta}$ and $\mathbf{X}$ are
fixed, is simply given by

$s_{n,k}=\max\\{0,P_{k}\\!-\\!\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}_{n}\mathbf{Z}_{k}^{{\dagger}})\\},\forall
k\in\mathcal{K}$

and thus a projection step for $\mathbf{s}$ is not necessary. The proposed
APGM to find a stationary solution to (3) for fixed $\boldsymbol{\upsilon}$
and $\rho$ is summarized in Algorithm 1. Note that in Algorithm 1, $\mu_{n}$
and $\alpha_{n}$ denote the step size corresponding to the gradient step for
the $\boldsymbol{\theta}$ and $\mathbf{X}$ update in iteration $n$,
respectively. Appropriate values of $\mu_{n}$ and $\alpha_{n}$ can be obtained
by a backtracking ___line search. Specifically,_ we can set
$\mu_{n}=\mu_{n-1}\gamma^{i_{n}}$ where $i_{n}$ is the smallest positive
integer such that

$\displaystyle\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n-1},\boldsymbol{\theta}_{n},\mathbf{s}_{n-1})\geq\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n-1},\boldsymbol{\theta}_{n-1},\mathbf{s}_{n-1})$
$\displaystyle\qquad\qquad\qquad+\langle\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n-1},\boldsymbol{\theta}_{n-1},\mathbf{s}_{n}),\boldsymbol{\theta}_{n}\\!-\\!\boldsymbol{\theta}_{n-1}\rangle$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\dfrac{1}{\mu_{n-1}\gamma^{i_{n}}}\left\|\boldsymbol{\theta}_{n}-\boldsymbol{\theta}_{n-1}\right\|^{2},$
(8)

where $\langle\mathbf{x},\mathbf{y}\rangle\triangleq
2\Re(\mathbf{x}^{{\dagger}}\mathbf{y})$ and $\gamma<1$. A similar routine can
be followed to obtain $\alpha_{n}$.

After solving (3) for given $(\boldsymbol{\upsilon},\rho)$, the penalty
parameter, $\rho$, is decreased and the Lagrange multipliers are updated. The
overall description of the proposed PDDGP algorithm to find a stationary
solution to (3) is outlined in Algorithm 2.

###  III-B Extension to Multiple SRs

We have presented the system model where there is a single SR for the sake of
simplicity. However, we remark that the extension to deal with multiple SRs is
straightforward. In this case, we can consider the maximization of the
weighted sum rate (WSR) of all SRs as in [9, 10] and the transmit power
constraint (2b) becomes the sum power constraint. It is easy to check that the
projected gradient step for $\boldsymbol{\theta}$ requires minimal
modifications. Also, the projection onto the set defined by the sum power
constraint still admits a water-filling solution. Thus, it is trivial to see
that Algorithm 2 can be easily modified to solve the resulting WSR
maximization problem. The detailed derivations are skipped here due to space
limitation and only the numerical results are presented for the multiple SR
case in the next section.

###  III-C Proof of Convergence

The convergence proof of Algorithm 2 follows the same arguments as given in
[14], and thus we only summarize the main points here. First, for given
$(\boldsymbol{\upsilon},\rho)$, Algorithm 1 generates a strictly increasing
objective sequence,
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$.
Since the feasible set is bounded, the iterates returned by Algorithm 1
converge to a limit point of (3a) to some accuracy $\epsilon_{k}$. Next, it
can be shown that the sequence
$\left\|\boldsymbol{\upsilon}_{k+1}-\boldsymbol{\upsilon}_{k}\right\|$ is
bounded. Thus, from the dual update it holds that
$\sqrt{\sum_{k\in\mathcal{K}}g_{k}^{2}(\mathbf{X},\boldsymbol{\theta},P_{k},s_{k})}=\rho\left\|\boldsymbol{\upsilon}_{k+1}-\boldsymbol{\upsilon}_{k}\right\|\to
0$ as $\rho\to 0$.

###  III-D Complexity Analysis

Figure 2: Convergence results for $P_{\max}=20$ dBm.

Figure 3: Average achievable rate for a single PR.

Figure 4: Average achievable sum rate versus $P_{\max}$.

Figure 5: Average achievable rate versus $N_{\mathrm{I}}$.

It is obvious that the complexity of the proposed PDDGP algorithm is dominated
by that of Algorithm 1. Thus the complexity of Algorithm 1 is obtained by
counting the required number of complex multiplications using the big-O
notation. First, the computational complexity associated with
$\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
in (5) is as follows. The complexity for computing the terms
$\mathbf{Z}_{\mathrm{R}}$,
$\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}}$,
$(\mathbf{I}+\mathbf{Z}_{\mathrm{R}}\mathbf{X}\mathbf{Z}_{\mathrm{R}}^{{\dagger}})^{-1}\mathbf{Z}_{\mathrm{R}}\
(\triangleq\boldsymbol{\Xi})$, $\boldsymbol{\Xi}\mathbf{X}$ and
$\operatorname{\mathrm{vec_{d}}}(\hat{\mathbf{H}}_{\mathrm{IR}}^{{\dagger}}\boldsymbol{\Xi}\mathbf{X}\mathbf{H}_{\mathrm{TI}}^{{\dagger}})$
is respectively $O(N_{\mathrm{R}}N_{\mathrm{I}}N_{\mathrm{T}})$,
$O(N_{\mathrm{R}}N_{\mathrm{T}}^{2}+N_{\mathrm{R}}^{2}N_{\mathrm{T}})$,
$O(N_{\mathrm{R}}^{3}+N_{\mathrm{R}}^{2}N_{\mathrm{T}})$,
$O(N_{\mathrm{R}}N_{\mathrm{T}}^{2})$ and
$O(N_{\mathrm{I}}N_{\mathrm{T}}N_{\mathrm{R}})$. Similarly, the complexity for
obtaining
$\operatorname{Tr}(\mathbf{Z}_{k}\mathbf{X}\mathbf{Z}_{k}^{{\dagger}})$ and
$\operatorname{\mathrm{vec_{d}}}(\mathbf{H}_{\mathrm{I}k}^{{\dagger}}\mathbf{Z}_{k}\mathbf{X}\mathbf{H}_{\mathrm{TI}}^{{\dagger}})$
is respectively
$O(N_{\mathrm{P}}N_{\mathrm{T}}(N_{\mathrm{I}}+N_{\mathrm{T}}))$ and
$O(N_{\mathrm{T}}^{2}N_{\mathrm{I}}+N_{\mathrm{P}}N_{\mathrm{T}}N_{\mathrm{I}})$.
Therefore, the total complexity for calculating
$\nabla_{\boldsymbol{\theta}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
is
$O(2N_{\mathrm{R}}N_{\mathrm{I}}N_{\mathrm{T}}+2N_{\mathrm{R}}N_{\mathrm{T}}^{2}+2N_{\mathrm{R}}^{2}N_{\mathrm{T}}+N_{\mathrm{R}}^{3}+(K+1)N_{\mathrm{P}}N_{\mathrm{T}}N_{\mathrm{I}}+KN_{\mathrm{T}}^{2}N_{\mathrm{I}}+KN_{\mathrm{P}}N_{\mathrm{T}}^{2}).$
We remark that, since the projection onto $\vartheta$ has negligible
complexity, the additional complexity to find an appropriate value of
$\mu_{n}$ using the backtracking line search in (8) can be ignored.

Next, the complexity to compute
$\nabla_{\mathbf{X}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
is $O(N_{\mathrm{T}}^{2}(N_{\mathrm{R}}+KN_{\mathrm{P}}))$. Additionally, the
complexity for projecting
$\mathbf{X}_{n\\!-\\!1}+\alpha_{n\\!-\\!1}\nabla_{\mathbf{X}}\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X}_{n\\!-\\!1},\boldsymbol{\theta}_{n},\mathbf{s}_{n\\!-\\!1})$
onto the feasible set $\mathcal{X}$ is $O(N_{\mathrm{T}}^{3})$. Therefore, the
computational complexity for the projected gradient step for $\mathbf{X}$ is
$O(N_{\mathrm{R}}N_{\mathrm{T}}^{2}+KN_{\mathrm{P}}N_{\mathrm{T}}^{2}+N_{\mathrm{T}}^{3}).$

Moreover, it is straightforward to observe that computing $\mathbf{s}_{n}$
requires
$O(K(N_{\mathrm{P}}N_{\mathrm{I}}N_{\mathrm{T}}+N_{\mathrm{P}}N_{\mathrm{T}}^{2}))$
complex multiplications. Therefore, the per-iteration complexity of Algorithm
1 is
$O(N_{\mathrm{R}}^{3}+N_{\mathrm{T}}^{3}+2N_{\mathrm{R}}^{2}N_{\mathrm{T}}+N_{\mathrm{T}}^{2}(3N_{\mathrm{R}}+KN_{\mathrm{I}}+3KN_{\mathrm{P}})+N_{\mathrm{T}}N_{\mathrm{I}}(2N_{\mathrm{R}}+2KN_{\mathrm{P}})).$

It is noteworthy that in a practical IRS-assisted MIMO underlay spectrum
sharing system,
$N_{\mathrm{I}}\gg\max\\{N_{\mathrm{T}},N_{\mathrm{R}},N_{\mathrm{P}}\\}$, and
therefore, the per-iteration complexity of Algorithm 1 can be approximated by
$O(N_{\mathrm{I}}(KN_{\mathrm{T}}^{2}+2N_{\mathrm{R}}N_{\mathrm{T}}+2KN_{\mathrm{P}}))$,
which is linear w.r.t. $N_{\mathrm{I}}$. Note that the complexity of the
algorithms proposed in both [9] and [10] grows as
$O(N_{\mathrm{I}}^{3}+KN_{\mathrm{I}}^{2}\max\\{N_{\mathrm{P}},N_{\mathrm{T}}\\})$
which is significantly higher than that of our proposed method.

##  IV Results and Discussion

In this section, the results of numerical experiments are presented to
evaluate the performance of the system under consideration. We consider $K=4$
PRs (unless stated otherwise), and the location of the ST, SR, IRS and the
four PRs in Cartesian coordinates are respectively given by
$(300\,\mathrm{m},0\,\mathrm{m})$, $(600\,\mathrm{m},0\,\mathrm{m})$,
$(300\,\mathrm{m},30\,\mathrm{m})$ and
$\\{(0\,\mathrm{m},0\,\mathrm{m}),\,(0\,\mathrm{m},5\,\mathrm{m}),\,(0\,\mathrm{m},10\,\mathrm{m}),\,(0\,\mathrm{m},15\,\mathrm{m})\\}$.
The small-scale fading is assumed to be Rayleigh distributed, and that the
path loss between two nodes is modeled as
$\mathrm{PL}=(-30-10\xi\log_{10}(d/d_{0}))$ dB, where $\xi$ is the path loss
exponent, $d$ is the distance between the nodes, and $d_{0}$
$(=1\,\mathrm{m})$ is the reference distance. We assume that the path loss
exponent for the ST-SR and ST-PR links is 3.75, whereas that for the ST-IRS,
IRS-SR and IRS-PR links is 2.2 [9]. Furthermore, it is assumed that
$N_{\mathrm{R}}=4$, $N_{\mathrm{P}}=4$, $P_{k}=10^{-13}$ W $(\forall
k\in\mathcal{K})$, the noise power spectral density $\mathscr{N}_{0}=-174$
dBm/Hz and the system bandwidth is $B=10$ MHz. In Figs. 5-5, the average
achievable rate is computed over $100$ channel realizations. To run Algorithm
2, we start with $\rho=10$ and decrease it as $\rho\leftarrow\kappa\rho$ where
$\kappa=0.1$ when the relative objective progress of
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
[14, Eqn. (48)] in Algorithm 1 is less than $10^{-5}$. For Algorithm 1 we set
$\gamma=0.5$ in the backtracking line search. Algorithm 2 terminates when the
difference between $R(\mathbf{X},\boldsymbol{\theta})$ and
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
is less than $10^{-5}$.

In Fig. 5, the convergence of Algorithms 1 and 2 are shown for
$N_{\mathrm{I}}=64$ and different values of $N_{\mathrm{T}}$. In the figure,
each iteration is comprised of one update of $\mathrm{\mathbf{X}}$,
$\boldsymbol{\theta}$ and $\mathrm{\mathbf{s}}$. Each color in Fig. 5 in fact
represents the convergence of Algorithm 1 for a fixed $\rho$ whose value is
clearly indicated in the figure. We recall that Algorithm 1 aims to maximize
the augmented objective
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$,
and thus, for a given $\rho$, monotonic increase is only guaranteed for
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$,
but not for the original objective $R(\mathbf{X},\boldsymbol{\theta})$. This
point is clearly seen in Fig. 5. We also observe that as the number of
iterations increases, the difference between
$R(\mathbf{X},\boldsymbol{\theta})$ and the augmented objective
$\hat{R}_{\boldsymbol{\upsilon},\rho}(\mathbf{X},\boldsymbol{\theta},\mathbf{s})$
decreases and approaches zero for sufficiently small $\rho$, in accordance
with the principle of a penalty method. In the figure, we also show the run
time for Algorithm 2, where the algorithm converges before 150 ms for
$N_{\mathrm{T}}=8$. Note that we have implemented the algorithm using Python
v3.9.7 on a Linux PC with 7.5 GiB memory and Intel Core i5-7200U CPU.

In Fig. 5, we demonstrate the superiority of the proposed PDDGP algorithm over
the BCD-IAM and WMMSE approaches proposed in [9] and [10], respectively, for
$N_{\mathrm{I}}=64$. Note that here only the first PR located at
$(0\,\mathrm{m},0\,\mathrm{m})$ is considered since the algorithm in [9] was
proposed for a single PR. It is evident from the figure that all the schemes
result in a similar performance for small values of $P_{\max}$; however, the
proposed PDDGP algorithm offers a significantly higher average achievable rate
for large values of $P_{\max}$ compared to that achieved via the solutions
proposed in [9] and [10]. Also, it is expected that the SU rate performance
achieved via the BCD-IAM and WMMSE schemes is the same because both of these
solutions are derived from similar principles. The main reason for the
inferior performance of the BCD-IAM and WMMSE for large $P_{\max}$ is that as
$P_{\max}$ increases, the SU rate performance is mostly dominated by the IPCs
at the PRs. In other words, the IPCs become binding for large $P_{\max}$.
Consequently, the coupling between $\mathbf{X}$ and $\boldsymbol{\theta}$
becomes stronger, and it is known that AO-based optimization cannot guarantee
a stationary solution under such a scenario.

Next, in Fig. 5, the average achievable sum rate versus the maximum transmit
power budget at the ST ($P_{\max}$) is shown for $N_{\mathrm{I}}=64$. In
particular, we compare the proposed PDDGP algorithm to the WMMSE-based AO
algorithm introduced in [10] for the case of two SRs located at
$(600\,\mathrm{m},0\,\mathrm{m})$ and $(600\,\mathrm{m},5\,\mathrm{m})$. In
Fig. 5, it can be clearly seen that the proposed PDDGP always outperforms the
WMMSE-based algorithm. Note that the BCD-IAM and the WMMSE-based AO are
similar in principle and thus share the same drawbacks. This explains why the
WMMSE-based AO in [10] is inferior to our proposed algorithm. However, note
that different from the case of a single PR in Fig. 5, the PDDGP algorithm
outperforms the WMMSE algorithm even in the lower $P_{\max}$ regime for the
multiple PRs case, and the gains are more significant for larger $P_{\max}$.
The benefit of an IRS-assisted MIMO system over that of the system without IRS
can also be observed from Fig. 5, as the former results in a considerably
higher achievable rate.

To investigate the effect of the number of reflecting elements at the IRS, in
Fig. 5 we plot the average achievable rate achieved by the proposed PDDGP
algorithm for one SR as a function of $N_{\mathrm{I}}$ for
$N_{\mathrm{T}}=16$, $N_{\mathrm{R}}=4$, $N_{\mathrm{P}}=2$ and $K=4$.
Interestingly, for an exponential increase in $N_{\mathrm{I}}$, the achievable
rate of the MIMO system also increases near-exponentially. This occurs because
for large values of $N_{\mathrm{I}}$, the IRS can create a highly-focused beam
to increase the signal-to-noise ratio (SNR) at the SR while satisfying the
interference constraints at the PRs, resulting in a higher achievable rate at
the SR.

##  V Conclusion

In this letter, we have presented the optimization problem for achievable rate
maximization in an IRS-assisted SU MIMO system in underlay spectrum sharing
with multiple PRs, subject to a transmit power constraint at the ST, unit
modulus constraints at the IRS and multiple interference power constraints at
the PRs. A PDDGP algorithm has been proposed to jointly optimize the transmit
beamforming and the IRS phase shifts. The numerical results established the
superiority of the proposed algorithm over two existing approaches in terms of
the achievable rate. We also showed that the complexity of the proposed PDDGP
algorithm is considerably lower than that of the existing approaches.

## References

  * [1] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: A tutorial,” IEEE Trans. Commun., vol. 69, no. 5, pp. 3313–3351, 2021.
  * [2] P. Yang, L. Kong, and G. Chen, “Spectrum sharing for 5G/6G URLLC: Research frontiers and standards,” IEEE Commun. Stand. Mag, vol. 5, no. 2, pp. 120–125, 2021.
  * [3] V. Kumar, Z. Ding, and M. F. Flanagan, “On the performance of downlink NOMA in underlay spectrum sharing,” IEEE Trans. Veh. Technol., vol. 70, no. 5, pp. 4523–4540, 2021.
  * [4] X. Guan, Q. Wu, and R. Zhang, “Joint power control and passive beamforming in IRS-assisted spectrum sharing,” IEEE Commun. Lett., vol. 24, no. 7, pp. 1553–1557, 2020.
  * [5] D. Xu, X. Yu, and R. Schober, “Resource allocation for intelligent reflecting surface-assisted cognitive radio networks,” in IEEE SPAWC, pp. 1–5, IEEE, 2020.
  * [6] J. Yuan, Y.-C. Liang, J. Joung, G. Feng, and E. G. Larsson, “Intelligent reflecting surface-assisted cognitive radio system,” IEEE Trans. Commun., vol. 69, no. 1, pp. 675–687, 2021.
  * [7] L. Zhang, C. Pan, Y. Wang, H. Ren, and K. Wang, “Robust beamforming design for intelligent reflecting surface aided cognitive radio systems with imperfect cascaded CSI,” IEEE Trans. on Cogn. Commun. Netw., to appear.
  * [8] S. F. Zamanian, S. M. Razavizadeh, and Q. Wu, “Vertical beamforming in intelligent reflecting surface-aided cognitive radio networks,” IEEE Wireless Communications Letters, vol. 10, no. 9, pp. 1919–1923, 2021.
  * [9] L. Zhang, Y. Wang, W. Tao, Z. Jia, T. Song, and C. Pan, “Intelligent reflecting surface aided MIMO cognitive radio systems,” IEEE Trans. Veh. Technol., vol. 69, no. 10, pp. 11445–11457, 2020.
  * [10] W. Jiang, Y. Zhang, J. Zhao, Z. Xiong, and Z. Ding, “Joint transmit precoding and reflect beamforming design for IRS-assisted MIMO cognitive radio systems,” IEEE Trans. Wireless Commun., to appear.
  * [11] L. Dong, H.-M. Wang, H. Xiao, and J. Bai, “Secure intelligent reflecting surface assisted MIMO cognitive radio transmission,” in IEEE WCNC, pp. 1–6, 2021.
  * [12] V. Kumar, M. F. Flanagan, D. B. Da Costa, and L.-N. Tran, “On the secrecy rate of downlink NOMA in underlay spectrum sharing with imperfect CSI,” in 28th International Conference on Telecommunications (ICT), pp. 167–173, 2021.
  * [13] N. S. Perović, L.-N. Tran, M. Di Renzo, and M. F. Flanagan, “Achievable rate optimization for MIMO systems with reconfigurable intelligent surfaces,” IEEE Trans. Wireless Commun., vol. 20, no. 6, pp. 3865–3882, 2021.
  * [14] Q. Shi, M. Hong, X. Fu, and T.-H. Chang, “Penalty dual decomposition method for nonsmooth nonconvex optimization–Part I: Algorithms and convergence analysis,” IEEE Trans. Signal Process., vol. 68, pp. 4108–4122, 2020.
  * [15] A. Hjørungnes, Complex-valued matrix derivatives: With applications in signal processing and communications.  Cambridge University Press, 2011.