# Intent-driven Intelligent Control and Orchestration in O-RAN Via
Hierarchical Reinforcement Learning

Md Arafat Habib1, Hao Zhou1, Pedro Enrique Iturria-Rivera1, Medhat Elsayed2,
Majid Bavand2,
Raimundas Gaigalas2, Yigit Ozcan2 and Melike Erol-Kantarci1, 2Ericsson Inc.,
Ottawa, Canada
Emails:{mhabi050, hzhou098, pitur008<EMAIL_ADDRESS>
{medhat.elsayed, majid.bavand, raimundas.gaigalas<EMAIL_ADDRESS>1School of Electrical Engineering and Computer Science, University of Ottawa,
Ottawa, Canada

###### Abstract

rApps and xApps need to be controlled and orchestrated well in the open radio
access network (O-RAN) so that they can deliver a guaranteed network
performance in a complex multi-vendor environment. This paper proposes a novel
intent-driven intelligent control and orchestration scheme based on
hierarchical reinforcement learning (HRL). The proposed scheme can orchestrate
multiple rApps or xApps according to the operator’s intent of optimizing
certain key performance indicators (KPIs), such as throughput, energy
efficiency, and latency. Specifically, we propose a bi-level architecture with
a meta-controller and a controller. The meta-controller provides the target
performance in terms of KPIs, while the controller performs xApp orchestration
at the lower level. Our simulation results show that the proposed HRL-based
intent-driven xApp orchestration mechanism achieves $7.5\%$ and $21.4\%$
increase in average system throughput with respect to two baselines, i.e. a
single xApp baseline and a non-machine learning-based algorithm, respectively.
Similarly, $17.3\%$ and $37.9\%$ increase in energy efficiency is observed in
comparison to the same baselines.

###### Index Terms:

O-RAN, rApps, xApp, hierarchical reinforcement learning, orchestration

##  I Introduction

Open radio access network (O-RAN) facilitates openness and intelligence to
support diverse traffic types and their requirements in 5G and beyond
networks[1] as well as, multi-vendor RAN deployments. In a multi-vendor
environment, rApps and xApps can be hosted in a non-real-time RAN intelligent
controller (non-RT-RIC) and near-real-time RAN intelligent controller (near-
RT-RIC). In the literature, xApps and rApps have been studied for resource and
power allocation, beamforming and management, cell sleeping, traffic steering
and so on[2, 3, 4]. Advanced reinforcement learning (RL) algorithms can be
used to develop intelligent network functions in O-RAN. However, a multi-rApp
or a multi-xApp scenario with a variety of AI-enabled Apps will require
intelligent control and orchestration among the Apps to avoid performance
degradation.

Note that, we focus on xApps as a case study but our work generalizes to rApps
as well. To elevate autonomy in O-RAN via xApp orchestration, intent-driven
network optimization goals can play a pivotal role. The intent is defined as
an optimization goal that is a high-level command given by the operator
usually in plain language and it determines a key performance indicator (KPI)
that the network should meet, such as “increase throughput by $10\%$” or
“increase energy efficiency by $5\%$” [5]. To better support autonomous
orchestration of the xApps in a multi-vendor environment, emphasis on
operators’ intents is crucial [6]. Intents aid in achieving agile, flexible,
and simplified configuration of the wireless networks with minimum possible
intervention. Furthermore, intelligent intent-driven management has the
ability to constantly acquire knowledge and adjust to changing network
conditions by utilizing extensive real-time network data. The inclusion of
intent-driven goals for intelligent xApp control and orchestration is a
promising yet highly complex task, since there are multiple vendors involved
with different network functions and intents may trigger conflicting
optimization goals in sub-systems. There are a few works on conflict
mitigation or xApp cohabitation in O-RAN. For instance, Han et al. propose a
conflict mitigation scheme among multiple xApps using team learning [7], and
Polese et al. propose a machine learning (ML)-based pipeline for the
cohabitation of multiple xApps in an O-RAN environment [8]. The work outlined
in [9] introduces a method for achieving automation throughout the entire life
cycle of xApps, beginning with the utilization scenario, requirements, design,
verification, and ultimately, the deployment within networks. However, the
operator intent is not involved in these works.

To this end, we propose a hierarchical reinforcement learning (HRL) method for
intent-driven xApp orchestration. Different from the previous works, the
proposed scheme has a bi-level architecture, where we can pass the intents to
the top-level hierarchy, and process it as optimization goals for the lower-
level controller to control and orchestrate xApps. Orchestration can avoid
xApp conflicts and improve performance by combining xApps with similar
performance objectives. The proposed method is compared with two baselines:
non-machine learning (non-ML) solution and a single xApp scenario. Our
simulation results show that the proposed HRL-based intent-driven xApp
orchestration mechanism achieves $7.5\%$ and $21.4\%$ increase in average
system throughput along with $17.3\%$ and $37.9\%$ increase in energy
efficiency, compared to the single xApp and non-ML baselines, respectively.

The rest of the paper is organized as follows: Section II discusses the
related works, followed by Section III which presents the system model
elaborately. The proposed HRL-based xApp orchestration in O-RAN is covered in
Section IV. Performance analysis and comparison of the proposed method along
with the baselines are presented in Section V. Lastly, we present our
conclusions in Section VI.

##  II Related work

There are a few works that investigate ML-based xApps for RAN optimization and
control. Polese et al. propose an ML pipeline for multiple xApps in an O-RAN
environment [8]. Han et al. propose a conflict mitigation scheme among deep
reinforcement learning (DRL)-based power allocation and resource allocation
xApps [7]. Polese et al. propose an Orchest-RAN scheme, in which network
operators can specify high-level control objectives in non-RT-RIC to sort out
the optimal set of data-driven algorithms to fulfill the provided intent [10].
While the work presented in [10] focuses on selecting the appropriate machine
learning models and their execution locations for given inputs from the
operator, it does not put emphasis on the network operator’s goals as
optimization objectives to select and orchestrate xApps.

An intent-driven orchestration of cognitive autonomous networks of RAN
management is presented in [11], where the authors propose a generic design of
intent-based management for controlling RAN parameters and KPIs. Zhang et al.
propose an intent conflict resolution scheme to realize conflict avoidance in
machine learning-based xApps [12]. A graph-based solution is proposed in [13]
to determine the specific network function required to fulfill an intent.

Compared with existing literature, the main contribution of this work is that
we propose an HRL scheme for intent-driven orchestration of xApps. The HRL
scheme can well fit the inherent O-RAN hierarchy with non-RT-RIC and near-RT-
RIC, and intent-based orchestration enables higher flexibility for network
control and management. The intents from the human level operator are provided
as goals for the system to achieve, which leads to the orchestration of xApps
to achieve the provided goal.

##  III System Model

###  III-A System Model

We consider an O-RAN-based downlink orthogonal frequency division multiplexing
cellular system having $B$ BSs serving $U$ users simultaneously, and multiple
small cells in the system are within the range of a macro cell. There are $K$
classes of traffic in the system and users are connected with multiple RATs
via dual connectivity. There are $Q$ classes of RATs
($q_{1},q_{2},...,q_{n}$), where $q$ represents a certain access technology
(LTE, 5G, etc.). The wireless system model considered in this work is
presented in Fig. 1. RIC platforms in the figure (non and near-RT-RIC) can
host rApps and xApps which are control and optimization applications operating
at different time scales.

Figure 1: O-RAN based system model with macro cell and small cells.

We design three xApps, namely traffic steering, cell sleeping, and beam
forming xApps. In each xApp, we apply deep reinforcement learning for
optimization within this xApp, which will be introduced in the following.

####  III-A1 Traffic Steering xApp

The traffic steering xApp aims to achieve a simultaneous balance of QoS
requirements for various traffic classes by introducing a traffic steering
scheme based on Deep Q-Network (DQN) [14]. We design the reward and state
functions to ensure satisfactory performance, focusing on two essential KPIs:
network delay and average system throughput. Traffic can be steered to a
certain BS based on load experienced, link quality, and traffic type. The
details of this xApp can be found in [2].

####  III-A2 Cell Sleeping xApp

The cell sleeping xApp is designed to reduce power consumption in the system
by turning off idle or less busy BSs. The xApp can perform cell sleeping based
on traffic load ratios and queue length of each BS. The energy consumption
model for the BS is:

$P_{in}=\left\\{\begin{array}[]{lcl}P_{0}+\delta_{p}P_{out},&&0<P_{out}\leq
P_{max},\\\ P_{sleep},&&P_{out}=0,\end{array}\right.$ (1)

where $P_{0}$ is the fixed power consumption, $\delta_{p}$ is the slope of
load-dependent power consumption, $P_{out}$ is the transmission power,
$P_{max}$ is the maximum transmission power, and $P_{sleep}$ is the constant
power consumption in sleep mode [15].

The goal of the cell sleeping xApp is to maximize energy efficiency as much as
possible without overloading the active BSs. The optimization goal is
formulated as follows:

$\begin{split}\max_{P_{b}}\quad\frac{\sum_{u\in{U_{o}}}\sum_{b\in
B}T_{u,b}}{P_{b}}-\theta b_{u},\\\
s.t.\quad\quad\quad\quad(\ref{eq2}),\quad\quad\quad\quad\quad\quad\quad\\\
\end{split}$ (2)

where $U_{o}$ is the set of the user equipments (UEs) connected to a certain
BS, $T$ represents the throughput, $\theta$ is the penalty factor to reduce
overloading, and $b_{u}$ is the number of the BSs overloaded. Turning off the
BSs can greatly decrease power consumption. It reduces the number of BSs
active that are serving the live network traffic. This poses a risk of
overloading the active BSs. Therefore, the penalty factor related to the
number of BSs has been introduced to avoid excessive overloading.

To address the formulated problem, the following MDP is formulated:

  * •

State: The set of the state consists of: $S=\\{q_{L},L_{R}\\}$. $L_{R}$
represents the traffic load ratio of a BS, $b$. The second element of the
state space is the queue length of the BSs representing the load level.

  * •

Action: Turning the BSs on and off are put into the action set for the DQN
implementation. $A=\\{ON,OFF\\}$.

  * •

Reward: The reward function is the same as eq. (2).

####  III-A3 Beamforming xApp

The third xApp is the beamforming xApp. We deploy band-switching BSs from 3.5
GHz to mmWave frequencies [16]. This allows us to support high throughput
traffic like enhanced mobile broadband (eMBB) via accurate intelligent
beamforming. This xApp can control power based on the location of the UE and
it uses minimum transmission power needed which is energy efficient. The xApp
employs analog beamforming, and a multi-antenna setup is adopted where each BS
deploys a uniform linear array (ULA) of $M$ antennas [17]. The beamforming
weights of every beamforming vector are implemented using constant modulus
phase shifters. We also assume that there is a beam steering-based codebook,
$F$, from which every beamforming vector is selected [17].

Every BS $l$ has a transmit power $P_{TX,l}\in P$, where $P$ is the set of
candidate transmit powers. We want to optimize two parameters: throughput and
energy efficiency using this xApp. To obtain such a goal, the following
optimization problem is addressed.

$\begin{split}max\sum_{l\in\\{1,2,..,L\\}}\left[c_{1}\left(\frac{T_{k,b}}{T_{QoS}}\right)+c_{2}\left(\frac{\varepsilon}{\varepsilon_{max}}\right)\right],\\\
s.t.\quad\quad\quad\quad\quad\quad P_{TX,l}[t]\in
P,\quad\quad\quad\quad\quad\quad\\\ f_{l}[t]\in
F,\quad\quad\quad\quad\quad\quad\quad\end{split}$ (3)

where $T_{k,b}$ is the throughput achieved by the system, $T_{QoS}$ is the
defined throughput requirement for a traffic type $k$, $\varepsilon$
represents the energy efficiency associated with the BS throughput and
transmission power, $\varepsilon_{max}$ is the maximum theoretical energy
efficiency, and $c_{1}$ and $c_{2}$ are the weight factors.

To solve the formulated problem, the following MDP is defined.

  * •

State: UE co-ordinates are used as set of states,
$S=\\{C_{UE1},C_{UE2},...,C_{UEN}\\}$.

  * •

Action: The action set consists of two elements:
$A=\\{\alpha(\chi_{n}),\delta_{n}\\}$. Here, $\chi$ is the steering angle, and
$\alpha(\chi_{n})$ is the array steering vector in the direction $\chi_{n}$ of
the $n$-th element in the codebook. $\delta_{n}$ accounts for the power level
change.

  * •

Reward: The reward function is the same as eq. (3) as presented before.

##  IV Proposed HRL-based xApp orchestration Scheme

RL problems can be formulated as MDPs where we have a set of states, actions,
transition probability, and a reward function ($S,A,T,R$). The RL agent in HRL
consists of two controllers: a meta-controller and a controller [18]. The MDP
for HRL has an added element which is denoted as a set of goals ($G$).
Depending on the current state, the meta-controller is responsible for
generating high-level goals ($G={g_{1},g_{2},...,g_{n}}$) for the controller.
After that, these goals are transformed into high-level policies. The
controller chooses low-level action $a$ according to the high-level policies.
This process from the controller yields an intrinsic reward ($r_{in}$).
Finally, an extrinsic reward ($r_{ex}$) is given to the meta-controller from
the environment and it will provide the controller with a new goal
($g_{\prime}$). This section will discuss the xApp orchestration scheme via
HRL.

###  IV-A xApp Coordination Using HRL

The proposed O-RAN-based system architecture is presented in Fig 2. RIC
platforms can host rApps and xApps which are applications operating at
different time scales. Three xApps have been defined in previous sections. The
rApp in the figure works as an input panel for the network operator, and it
can convert these inputs as goals to be optimized. Also, it works as the meta-
controller in the non-RT-RIC.

Let’s assume, $X$ is a set of xApps and $Y$ is the subset of $X$ having at
least one element (xApp in our case), that can optimize the network
performance based on the operator input. Let $I$ be the set of candidate KPIs
that a xApp can optimize and $Z$ be the set of QoS requirements the system has
to satisfy. Considering all these assumptions, the xApp orchestration problem
that we want to address can be formulated as follows:

$\begin{split}max\sum_{i\in I}\sum_{z\in
Z}(P_{i}-\rho\xi_{z}),\quad\quad\quad\\\
s.t.\quad\forall(X)\exists(Z):V(O)=1,\quad\quad\end{split}$ (4)

where $P$ is the intended performance metric the operator intends to improve.
$\rho$ is the penalty parameter for QoS requirement violation, and $\xi_{z}$
is the number of UEs QoS requirements violated to. Lastly, $V(O)$ is the
proposition that “An xApp can improve a performance metric”, which is either
‘0’ or ‘1’.

As presented in Fig. 2, the rApp in the system is directly connected to the
user panel where the operator may provide input to the system. The operator
input is provided as the percentage of the increase related to a certain KPI.
For example, $x\%$ for throughput increase or $y\%$ for energy efficiency
increase or any other intent stated in natural language. The rApp has a
hierarchical deep Q-learning (h-DQN) framework [18]. The meta-controller (in
non-RT-RIC) takes the increased amount of throughput or increased amount of
energy efficiency as a goal, observes the state in the environment and
provides both the goal and states to the controller in near-RT-RIC having a
bundle of xApps. This type of data passing is done via the A1 interface by
which both the non and near-RT-RIC are connected. The controller takes the
action of choosing an xApp or a set of xApps based on the provided state and
goal. Following, we define the MDP for the meta-controller and controller to
address the xApp orchestration problem formulated in eq. (4).

  * •

State: The set of states consists of traffic flow types of different users in
the network. UEs having similar traffic types are grouped together.
$S=\\{T_{voice},..,T_{urllc},...,T_{gaming},...,T_{eMBB},..\\}$. Elements in
this set stand for five different traffic types in the system. Both meta-
controller and controller share the same states.

  * •

Action: xApp or combination of xApp selection is considered as actions to be
performed by the controller which is defined as:
{$A_{xApp1},A_{xApp1,2},....,A_{xAppN}\\}$.

  * •

Intrinsic reward: The intrinsic reward function ($r_{in}$) for the controller
is: $r_{in}=P_{i}-\rho\xi_{z}$ which is similar to eq. (4).

  * •

Goal for the controller: Increased throughput or increased energy efficiency
level that can satisfy operator intent is passed to the controller as goals.
It is $G=\\{tp_{1},tp_{2},...,tp_{n}\\}$ for throughput increasing intents or
$G=\\{ee_{1},ee_{2},...,ee_{n}\\}$ for energy efficiency increasing intents.
Note that these goals can be generalized to other KPIs however for simplicity
we target throughput and energy efficiency.

  * •

Extrinsic reward: The meta-controller is responsible for the overall
performance of the system. Therefore, we have set the extrinsic reward
function for the meta-controller as the objective of the problem formulation
presented in eq. (4). The following equation is basically the summation of the
intrinsic reward over $\tau$ steps.

$r_{ex,\tau}=\frac{1}{n}\sum_{\tau=1}^{n}r_{in,\tau}\quad\forall(u)\in
U,\forall(b)\in B,$ (5)

Figure 2: Intent-based xApp orchestration with macro and micro cells having
different types of traffic.

The whole process of xApp orchestration can be summarized as follows:

  * •

Step 1: Operator’s intent is provided as input regarding which performance
metric is to be improved.

  * •

Step 2: These performance targets are provided to the controller in near-RT-
RIC by the meta controller rApp in the non-RT-RIC as goals to achieve.

  * •

Step 3: The controller selects an xApp or a combination of xApps to reach the
target performance as close as possible. The system learns based on the reward
it gets for such kind of xApp selection.

  * •

Step 4: Selected xApps with their own DRL-based functionalities optimize the
performance of the network as a response to the intent of the operator.

###  IV-B Baseline Algorithms

This section includes two baselines. The first baseline is a simulation of the
same network scenario based on the system model we have presented so far where
there is no intelligent DRL-based xApp to optimize the network. We use non-ML
algorithms. For comparing the throughput performance of the proposed HRL-based
system, we use a threshold-based traffic steering scheme proposed in [19]. It
uses a predefined threshold. The threshold is determined considering the load
at each station, channel condition, and user service type. The mean of all
these metrics is taken to obtain the threshold ($T$) values. Weighted
summation of the same parameters is taken to form a variable ($w$). Then, the
traffic is steered to another BS based on the $w$ and $T$ values. This
baseline does not include cell sleeping, therefore BSs are always on. In our
second baseline, we consider single xApp scenarios. For example, the proposed
HRL-based xApp orchestration mechanism is compared with single xApp scenarios
where only traffic steering xApp is in action, or the cell sleeping xApp is in
action.

##  V Performance Evaluation

###  V-A Simulation setup

A MATLAB-based simulation environment has been developed having one eNB and
four gNBs to serve as one macro-cell and four small cells. In total, we deploy
60 UEs with five different traffic types: voice, gaming, video, URLLC, and
eMBB. Different types of traffic in the system have variant requirements in
terms of different KPIs. QoS requirements of different traffic types have been
defined based on our previous work [20]. eMBB and URLLC traffic types have
been added additionally to test the system compatibility. For the eMBB traffic
type, we consider packet size, $T_{QoS}$, and $D_{QoS}$ to be 1500 bytes, 100
Mbps, and 15 ms, respectively [21]. Lastly, specifications related to delay
and packet size for the URLLC traffic are set to 32 bytes and 2.5 ms.

A 5G NSA mode having different types of RAT (LTE and 5G) in the simulation
environment work together. We deploy an architecture based on [22]. The
carrier frequency for LTE is set to be 800 MHz. For 5G NR small cells, band-
switching BSs are deployed at 3.5 GHz and 30 GHz. BS transmission power for
LTE and 5G NR is set to be 38 dBm and 43 dBm, respectively [23].

For the HRL implementation, the starting rate of learning is set to 0.95. In
order to maintain stable learning performance, we reduce the learning rate
periodically after a certain number of episodes. Additionally, the discount
factor used is 0.3. The simulation is conducted 10 times using MATLAB, and the
average outcomes are presented along with a 95% confidence interval.

###  V-B Simulation results

Figure 3: Impact of operator intents on throughput.

Figure 4: Impact of operator intents on energy efficiency.

Before conducting the performance evaluation of the proposed xApp
orchestration scheme, first, we present how the intent-oriented HRL-based
orchestration scheme works. Fig. 3 shows that the operator intent of “increase
throughput” leads to the selection of certain xApps. When there is a $5\%$
throughput increase intent from the operator, after a few time slots, there is
a sharp increase in throughput. This is because xApp1 (traffic steering xApp)
has been invoked. When a 5% increase is again given as an input, a combination
of xApp1 and xApp3 (intelligent beamforming xApp) is selected. When the
operator provides an intent to decrease power consumption by $5\%$, we can see
from Fig. 3 that there is a sharp decrease in throughput. This is because
xApp1 and xApp3 have been terminated at the 461-th time slot and xApp2 (cell
sleeping xApp) has been invoked.

Fig. 4 presents a similar graph to the previous one but this time it plots the
energy efficiency in the time axis. When there is an intent from the operator
to achieve “$10\%$ increase in energy efficiency”, we can see that there is an
initiation of xApp2 at the 131-th time slot. This xApp performs cell sleeping
and saves energy. For the next energy efficiency increase intent given by the
operator, it can be seen that both xApp2 and xApp3 are working together. The
proposed HRL-based algorithm has successfully orchestrated these two xApps for
the desired performance gain. Fig. 3 and 4 basically show the utility of the
proposed system. Not only it can induce operator intent as an optimization
goal, but also it can orchestrate xApps to gain desired performance output by
using the proper combination of xApps.

Fig. 5 shows the performance comparison between the proposed HRL-based xApp
orchestration scheme and the baseline scenarios in terms of average system
throughput. Results are obtained under a constant load of 6 Mbps. The proposed
orchestration scheme achieves a $21.4\%$ increase and $7.5\%$ increase in
average system throughput compared to the non-ML algorithm and single xApp
scenario (traffic steering xApp), respectively. It is because of the efficient
orchestration mechanism that involves multiple xApps that trigger the optimal
combination of xApps to reach better performance based on the operator intent.

Fig. 6 shows the performance comparison between the proposed HRL-based xApp
orchestration scheme and the baseline scenarios in terms of average energy
efficiency. The proposed orchestration scheme obtains a $17.3\%$ increase and
$37.9\%$ increase in average energy efficiency compared to the single xApp and
non-ML scenario (cell sleeping xApp), respectively. Similar to the former, it
is because of the HRL-based orchestration mechanism that incorporates multiple
xApps to achieve better performance based on the user intent. Also, note that
we use traffic steering in the former figure and cell sleeping in this
evaluation because they are specifically optimizing throughput and energy
respectively.

Figure 5: Performance comparison in terms of average system throughput.

Figure 6: Performance comparison in terms of energy efficiency.

##  VI Conclusions

In this paper, we show that the HRL-based intent-driven orchestration
mechanism is vastly effective in not only optimizing KPIs but also providing
great flexibility and control to the operator. In this study, we have
introduced a novel HRL-based xApp orchestration mechanism that can perform
xApp management and provide recommendations for the best combination of xApps
given the operator’s intent. The optimal xApp orchestration scheme has led to
a $7.5\%$ increase in average system throughput and a $17.3\%$ increase in
energy efficiency compared to single xApp usage with no orchestration. In our
future work, we plan to extend this orchestration to rApps and other xApps
with complex KPI interactions.

## Acknowledgement

This work has been supported by MITACS and Ericsson Canada, and NSERC Canada
Research Chairs and NSERC Collaborative Research and Training Experience
Program (CREATE) under Grant 497981.

## References

  * [1] L. Bonati, S. D’Oro, M. Polese, S. Basagni, and T. Melodia, “Intelligence and Learning in O-RAN for Data-Driven NextG Cellular Networks,” _IEEE Communications Magazine_ , vol. 59, no. 10, pp. 21–27, 2021.
  * [2] M. A. Habib, H. Zhou, P. E. Iturria-Rivera, M. Elsayed, M. Bavand, R. Gaigalas, S. Furr, and M. Erol-Kantarci, “Traffic Steering for 5G Multi-RAT Deployments using Deep Reinforcement Learning,” in _2023 IEEE 20th Consumer Communications & Networking Conference (CCNC)_, 2023, pp. 164–169.
  * [3] Y. Dantas, P. E. Iturria-Rivera, H. Zhou, M. Bavand, M. Elsayed, R. Gaigalas, and M. Erol-Kantarci, “Beam Selection for Energy-Efficient mmWave Network Using Advantage Actor Critic Learning,” 2023.
  * [4] H. Zhou, L. Kong, M. Elsayed, M. Bavand, R. Gaigalas, S. Furr, and M. Erol-Kantarci, “Hierarchical Reinforcement Learning for RIS-Assisted Energy-Efficient RAN,” in _GLOBECOM 2022 - 2022 IEEE Global Communications Conference_ , 2022, pp. 3326–3331.
  * [5] K. Mehmood, K. Kralevska, and D. Palma, “Intent-Driven Autonomous Network and Service Management in Future Cellular Networks: A Structured Literature Review,” _Computer Networks_ , vol. 220, p. 109477, 2023. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1389128622005114
  * [6] A. Leivadeas and M. Falkner, “A Survey on Intent-Based Networking,” _IEEE Communications Surveys & Tutorials_, vol. 25, no. 1, pp. 625–655, 2023.
  * [7] H. Zhang, H. Zhou, and M. Erol-Kantarci, “Team Learning-Based Resource Allocation for Open Radio Access Network (O-RAN),” in _ICC 2022 - IEEE International Conference on Communications_ , 2022, pp. 4938–4943.
  * [8] M. Polese, L. Bonati, S. D’Oro, S. Basagni, and T. Melodia, “ColO-RAN: Developing Machine Learning-based xApps for Open RAN Closed-loop Control on Programmable Experimental Platforms,” 2022.
  * [9] A. Kilks, M. Dryjansky, V. Ram, L. Wong, and P. Harvey, “Towards Autonomous Open Radio Access Networks,” _ITU Journal on Future and Evolving Technologies_ , vol. 4, no. 2, 2023.
  * [10] S. D’Oro, L. Bonati, M. Polese, and T. Melodia, “OrchestRAN: Network Automation Through Orchestrated Intelligence in the Open RAN,” in _IEEE INFOCOM 2022 - IEEE Conference on Computer Communications_ , 2022, pp. 270–279.
  * [11] A. Banerjee, S. S. Mwanje, and G. Carle, “An Intent-Driven Orchestration of Cognitive Autonomous Networks for RAN Management,” in _2021 17th International Conference on Network and Service Management (CNSM)_ , 2021, pp. 380–384.
  * [12] J. Zhang, J. Guo, C. Yang, X. Mi, L. Jiao, X. Zhu, L. Cao, and R. Li, “A Conflict Resolution Scheme in Intent-Driven Network,” in _2021 IEEE/CIC International Conference on Communications in China (ICCC)_ , 2021, pp. 23–28.
  * [13] E. J. Scheid, C. C. Machado, M. F. Franco, R. L. dos Santos, R. P. Pfitscher, A. E. Schaeffer-Filho, and L. Z. Granville, “INSpire: Integrated NFV-based Intent Refinement Environment,” in _2017 IFIP/IEEE Symposium on Integrated Network and Service Management (IM)_ , 2017, pp. 186–194.
  * [14] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-Level Control Through Deep Reinforcement Learning,” _nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [15] P. Ren and M. Tao, “A Decentralized Sleep Mechanism in Heterogeneous Cellular Networks with QoS Constraints,” _IEEE Wireless Communications Letters_ , vol. 3, no. 5, pp. 509–512, Oct. 2014.
  * [16] F. B. Mismar, A. Alammouri, A. Alkhateeb, J. G. Andrews, and B. L. Evans, “Deep Learning Predictive Band Switching in Wireless Networks,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 96–109, 2021.
  * [17] F. B. Mismar, B. L. Evans, and A. Alkhateeb, “Deep Reinforcement Learning for 5G Networks: Joint Beamforming, Power Control, and Interference Coordination,” _IEEE Transactions on Communications_ , vol. 68, no. 3, pp. 1581–1592, 2020.
  * [18] T. D. Kulkarni, K. Narasimhan, A. Saeedi, and J. B. Tenenbaum, “Hierarchical Deep Reinforcement Learning: Integrating Temporal Abstraction and Intrinsic Motivation,” _CoRR_ , vol. abs/1604.06057, 2016. [Online]. Available: http://arxiv.org/abs/1604.06057
  * [19] M. Khaturia, P. Jha, and A. Karandikar, “5G-Flow: A unified Multi-RAT RAN architecture for beyond 5G networks,” _Computer Networks_ , vol. 198, p. 108412, 2021.
  * [20] M. A. Habib, H. Zhou, P. E. Iturria-Rivera, M. Elsayed, M. Bavand, R. Gaigalas, Y. Ozcan, and M. Erol-Kantarci, “Hierarchical Reinforcement Learning Based Traffic Steering in Multi-RAT 5G Deployments,” 2023.
  * [21] A. Chagdali, S. E. Elayoubi, and A. M. Masucci, “Impact of Slice Function Placement on the Performance of URLLC with Redundant Coverage,” in _2020 16th International Conference on Wireless and Mobile Computing, Networking and Communications (WiMob)_ , 2020, pp. 1–6.
  * [22] P. Frenger and R. Tano. (2019) A Technical Look at 5G Energy Consumption and Performance. [Online]. Available: https://www.ericsson.com/en/blog/2019/9/energy-consumption-5G-nr
  * [23] E. Dahlman, S. Parkvall, and J. Skold, _4G, LTE-Advanced Pro and The Road to 5G, Third Edition_ , 3rd ed. USA: Academic Press, Inc., 2016.


# Elastic Functional Changepoint Detection

J. Derek Tucker and Drew Yarger  J. Derek Tucker  Statistical Sciences, Sandia
National Laboratories  Drew Yarger  Statistical Sciences, Sandia National
Laboratories

Elastic Functional Changepoint Detection

J. Derek Tucker 111Statistical Sciences, Sandia National Laboratories
222Corresponding Author - J. Derek Tucker<EMAIL_ADDRESS>Drew Yarger:
<EMAIL_ADDRESS>Drew Yarger1

Abstract

Detecting changepoints in functional data has become an important problem as
interest in monitory of climatologies and other various processing monitoring
situations has increased, where the data is functional in nature. The observed
data often contains variability in amplitude ($y$-axis) and phase ($x$-axis).
If not accounted for properly, incorrect changepoints can be detected, as well
as underlying mean functions at those changes will be incorrect. In this
paper, an elastic functional changepoint method is developed which properly
accounts for these types of variability. Additionally, the method can detect
amplitude and phase changepoints which current methods in the literature do
not, as they focus solely on the amplitude changepoint. This method can easily
be implemented using the functions directly, or to ease the computational
burden can be computed using functional principal component analysis. We apply
the method to both simulated data and real data sets to show its efficiency in
handling data with phase variation with both amplitude and phase changepoints.

Keywords: changepoint analysis, functional data analysis, functional principal
component analysis

Short title: Elastic Functional Changepoint Detection

##  1 Introduction

The statistical analysis of functional time series has become increasingly
important to many scientific fields including climatology (Shang and Hyndman,
2011), finance (Kokoszka and Zhang, 2012), geophysics (Hörmann and Kokoszka,
2012), demography (Hyndman and Booth, 2008), manufacturing (Woodall, 2007),
and environmental modeling (Fortuna et al., 2020). A functional time series is
a sequence of functions (i.e., infinite objects), observed over time.
Functional time series are analogous to univariate or multivariate time
series, except that we observe a continuous function at each point in time
(Bosq, 2012). Similar to univariate and multivariate time series, a functional
time series can experience abrupt changes in its data generating mechanism.
These abrupt changes are known as changepoints and can complicate further
analysis by invalidating stationarity assumptions. They also can lead to
interesting analysis by determining when data has patterns which are not
homogeneous over time, or when a forcing has caused the data generation
mechanism to change.

A problem with current methods in the Functional Data Analysis (FDA)
literature on changepoint detection has been solely focused on cross-sectional
analysis. This process completely ignores phase variability that may exist in
the data. This is an issue because phase variability could be just a nuisance
parameter or actually contain information. In Berkes et al. (2009) a
Cumulative Sum (CUSUM) test was proposed for independent functional data,
which was further studied in Aue et al. (2009), where its asymptotic
properties were developed. This test was then extended to weakly dependent
functional data by Hörmann et al. (2010). Zhang et al. (2011) introduced a
test for changes in the mean of weakly dependent functional data using self-
normalization to alleviate the use of asymptotic control. Later, Sharipov et
al. (2016) similarly developed a sequential block bootstrap procedure for
these methods. More recently, Gromenko et al. (2017) considered changes in
spatially-correlated functional data, and Aue et al. (2018) proposed a fully
functional method for finding a change in the mean without losing information
due to functional principal component analysis (fPCA) truncation.

Recently, Tucker et al. (2013) and Srivastava and Klassen (2016) have
demonstrated the need to account for phase variability in functional data. If
the phase variability is ignored, subsequent statistical analysis and modeling
is affected as the metric between functions is degenerate. Specifically,
within the changepoint detection problem is the computation of the mean
function. When phase variability is ignored, the mean computed will not be
representative of the underlying data generating mechanism. An example of this
is shown in Figure 1 where Panel (a) presents a set of functional data that
contains a large amount of phase variability. Panel (b) shows the computed
mean if the phase variability is ignored, and it is easily seen that the mean
is not representative of the underlying data generating mechanism. Panel (c)
shows the alignment of the data using the method in Tucker et al. (2013) and
the resulting Karcher Mean in Panel (d) where the phase variability is
accounted for. The mean in Panel (d) is representative of the underlying data
generating mechanism.

(a) Original Data

(b) Original Mean

(c) Aligned Data

(d) Karcher Mean

Figure 1: Example of using proper mean calculation on functional data.

In this paper we present a functional changepoint method that accounts for the
amplitude and phase variability. We demonstrate that we account for this
variability in two different ways, utilizing the appropriate fPCA (Tucker et
al., 2013; Lee and Jung, 2017) or utilizing a fully functional approach
similar to Aue et al. (2018) where the mean is computed using the Karcher mean
under the Fisher-Rao Metric (Srivastava and Klassen, 2016). This metric which
we will call elastic is a proper distance and invariant to random warping.
Through this we can separate out the amplitude and phase variability for
subsequent changepoint analysis. Therefore, we are able to find changepoints
where the underlying mean changes in amplitude, phase, or a combination of
both.

The rest of the paper is organized as follows: Section 2 presents the elastic
amplitude and phase changepoint methods, and provides implementation using
fully functional or functional principal component analysis. Section 3
provides analysis of the methodologies on simulated and climate monitoring
data sets. Finally, concluding remarks and future work are presented in
Section 4.

##  2 Elastic Functional Changepoint Analysis

###  2.1 The Functional Changepoint Problem

Assume we have real-valued functions $f_{1},\dots,f_{n}$ that are absolutely
continuous on the interval $[0,1]$. Using this interval is not a restriction
and can be easily extended to any interval. The standard changepoint problem
(Berkes et al., 2009; Aue et al., 2018) assumes the data is from the following
model:

$f_{i}=\mu+\delta\vmathbb{1}(i>k^{*})+\epsilon_{i}.$

Where the point $k^{*}$ labels the time of the unknown mean change,
$\epsilon_{i}$ is a zero mean random error term with may be i.i.d. or
potentially weakly dependent as in Aue et al. (2018). The changepoint problem
then becomes a hypothesis testing problem of

$H_{0}:\delta=0~{}~{}\text{vs}~{}~{}H_{A}:\delta\neq 0.$

Aue et al. (2018) considers the test statistic of $T_{n}=\max_{1\leq k\leq
n}\lVert S_{n,k}\rVert^{2}$, where

$\displaystyle S_{n,k}=\frac{1}{\sqrt{n}}\left(k\mu_{k}-k\mu_{n}\right),$

$\mu_{k}=k^{-1}\sum_{i=1}^{k}f_{i}$, and $\lVert\cdot\rVert$ is the standard
$\mathbb{L}^{2}$ norm on the functions $f_{i}$. Aue et al. (2018) demonstrate
asymptotic properties of this test statistic. Also, Berkes et al. (2009) and
subsequent work propose a functional principal component-based test statistic
that includes dimension reduction. These two approaches will motivate our test
statistics for functional data with phase and amplitude variability.

Since the test statistics of Aue et al. (2018) and Berkes et al. (2009) are
based functional statistics that do not use functional alignment, this cross-
sectional approach assumes there is no phase variability and that there is
only a random amplitude term. We know from Tucker et al. (2013) and Srivastava
and Klassen (2016) that functional data contains two types of variation,
amplitude and phase, and both are important to consider when evaluating
functional data for statistical analysis and modeling.

###  2.2 The Elastic Functional Changepoint Problem

We now formulate the changepoint problem for functional data with amplitude
and phase variation. We can include phase variation in the model in the
following form

$f_{i}=\left(\mu+\delta\vmathbb{1}(i>k^{*})+\epsilon_{i}\right)\circ\gamma_{i}^{-1},$

where $\gamma_{i}$ is an unknown phase variation. We will assume that
$\gamma_{i}$ is from the following group $\Gamma$, the phase space of the unit
interval $[0,1]$ defined as

$\Gamma=\\{\gamma:[0,1]\mapsto[0,1]\ |\ \gamma(0)=0\text{, }\gamma(1)=1\text{,
}\gamma\text{ is a diffeomorphism}\\}.$

This diffeomorphic constraint gives rise to the notion of elasticity because
the elements of $\Gamma$, i.e. phase functions, can only smoothly stretch and
contract portions of the unit interval so that it maps back to itself.
Additionally, we will use the Fisher-Rao metric and the square-root velocity
function (SRVF) for its attractive properties (see Srivastava and Klassen
(2016) for properties and appropriate proofs). The square-root velocity
function is defined as:

$q=sign(\dot{f})\sqrt{|\dot{f}|}.$

Therefore we propose the model

$q_{i}=\left((\mu_{q}+\delta\vmathbb{1}(i>k^{*})+\epsilon_{i}),\gamma_{i}^{-1}\right).$

where the action $(q,\gamma)=q\circ\gamma\sqrt{\dot{\gamma}}$. Here, we now
aim to detect a changepoint in the sequence of SRVF functions while
controlling for phase variability. We can also view the data as

$(q_{i},\gamma_{i})=\mu_{q}+\delta\vmathbb{1}(i>k^{*})+\epsilon_{i}.$

We wish to test the hypothesis

$H_{0}:\delta=0~{}~{}\text{vs}~{}~{}H_{A}:\delta\neq 0.$

One might also observe a change in phase of the functional data that would not
be detected with cross-sectional methods. In this setting, we consider a model
for $\gamma_{i}$ as

$\gamma_{i}=\mu_{\gamma}+\delta_{\gamma}\vmathbb{1}(i>k^{*})+\epsilon_{\gamma,i}$

where $\epsilon_{\gamma,i}$ is another zero-mean error term and test the
hypothesis

$H_{0}:\delta_{\gamma}=0~{}~{}\text{vs}~{}~{}H_{A}:\delta_{\gamma}\neq 0.$

In the next two subsections, we propose a test procedure for changepoints in
the amplitude mean and phase mean of functional data, respectively.

###  2.3 Amplitude changepoint

From Srivastava and Klassen (2016) the amplitude distance (which is a proper
distance) between two functions is defined as

$d_{a}(f_{1},f_{2})=\inf_{\gamma\in\Gamma}\|q_{1}-(q_{2},\gamma)\|.$

From a proper distance we can define the Karcher mean of functions as

$\mu_{q}^{k}=\operatorname*{arg\,min}_{q}\sum_{i=1}^{k}d_{a}(q,q_{i})^{2}.$

To determine if an amplitude changepoint has occurred, we can define a test
statistic using the (scaled) functional cumulative sum. This is similar to the
approach used in Aue et al. (2018). The test statistic is defined as

$S_{n,k}=\frac{1}{\sqrt{n}}\left(k\mu_{q}^{k}-k\mu_{q}^{n}\right)$ (1)

where $\mu_{q}^{n}$ is the Karcher mean of the entire data set and for a
specific value of $k$ we can compute the Karcher mean $\mu_{q}^{k}$. We can
also define a test statistic using the $\mathbb{L}^{2}$ norm to also test for
differences in means which is similar to the test statistic proposed in
Gromenko et al. (2017). The test statistic is defined as

$\Lambda_{2}=\frac{1}{n}\sum_{k=1}^{n}\|k\mu_{q}^{k}-k\mu_{q}^{n}\|^{2}.$ (2)

One can find the changepoint by using the max-point structural break detector,

$T_{n}=\max_{1\leq k\leq n}\|S_{n,k}\|^{2}$

to test the hypothesis $H_{0}$ versus $H_{A}$ and determine the estimated
changepoint $k^{*}$ when

$k^{*}=\min\\{k:\|S_{n,k}\|=\max_{1\leq k\leq n}\|S_{n,k}\|\\}.$

For the asymptotic properties of this statistic the reader is referred to Aue
et al. (2018). Under this model $T_{n}$ is distributed according to

$T_{n}\xrightarrow{d}\sup_{0\leq x\leq
1}\sum_{l=1}^{\infty}\lambda_{l}B_{l}^{2}(x)$

where $B_{l}$ are i.i.d. standard Brownian bridges defined on $[0,1]$ and
$\\{\lambda_{l}\\}$ are the eigenvalues on the long-run covariance kernel of
the error sequence $\epsilon_{i}$. Computation of critical values of $T_{n}$
can be done by the simulation of independent Brownian bridges and truncating
the infinite sum; see Section 3.2 of Aue et al. (2018) for more details on
this computation.

This tests for mean shape changes under phase variability. We can also easily
implement this using functional principal component analysis (fPCA) as in
Berkes et al. (2009). In order to properly account for the variability, we can
use the vertical fPCA and horizontal fPCA presented in Tucker et al. (2013).
These fPCA methods account for the variability by first separating the phase
and amplitude and then performing the fPCA on the spaces separately. Using
these methods one can construct the changepoint test in the amplitude space
using the SRVF, $q$, or specifically the aligned SRVF, $\tilde{q}$, and in the
phase space using the warping functions, $\gamma$.

We can use the vertical fPCA from Tucker et al. (2013) to test for an
amplitude change. We will first give a short review of this fPCA method below.

####  2.3.1 Vertical Functional Principal Components

Let $f_{1},\cdots,f_{n}$ be a given set of functions, and $q_{1},\cdots,q_{n}$
be the corresponding SRVFs, ${\mu}_{q}$ be their Karcher Mean, and let
$\tilde{q}_{i}$s be the corresponding aligned SRVFs using Algorithm 1 from
Tucker et al. (2013). In performing vertical fPCA, one also needs to include
the variability associated with the initial values, i.e., $\\{f_{i}(0)\\}$, of
the given functions. Since representing functions by their SRVFs ignores
vertical translation, this variable is treated separately. That is, a
functional variable $f$ is analyzed using the pair $(q,f(0))$ rather than just
$q$. This way, the mapping from the function space $\mathcal{F}$ to
$\mathbb{L}^{2}\times$ is a bijection. In practice, where $q$ is represented
using a finite partition of $[0,1]$, say with cardinality $T$, the combined
vector $h_{i}=[q_{i}~{}~{}f_{i}(0)]$ simply has dimension $(T+1)$ for fPCA. We
can define a sample covariance operator for the aligned combined vector
$\tilde{h}_{i}=[\tilde{q}_{i}~{}~{}f_{i}(0)]$ as

$K_{h}=\frac{1}{n-1}\sum_{i=1}^{n}(\tilde{h}_{i}-\mu_{h})(\tilde{h}_{i}-\mu_{h})^{\mathsf{T}}\in\mathbb{R}^{(T+1)\times(T+1)}\
,$

where $\mu_{h}=[\mu_{q}~{}~{}\bar{f}(0)]$. Taking the SVD,
$K_{h}=U_{h}\Lambda_{h}V_{h}^{\mathsf{T}}$, we can calculate the directions of
principal variability in the given SRVFs using the first $p\leq n$ columns of
$U_{h}$ and can then converted back to the function space $\mathcal{F}$, via
integration, for finding the principal components of the original functional
data. Moreover, we can calculate the observed principal scores as
$\eta_{i,l}=\left\langle\tilde{h}_{i},U_{h,l}\right\rangle$.

The test statistic would then follow directly from Berkes et al. (2009) with
test statistic

$T_{N}(x)=\frac{1}{N}\sum_{l=1}^{d}\lambda_{l}^{-1}\left(\sum_{1\leq i\leq
Nx}\eta_{i,l}-x\sum_{1\leq i\leq N}\eta_{i,l}\right)^{2}$ (3)

where $\lambda$ are the eigenvalues from $\Lambda_{h}$ and $\eta$ are the
scores from the fPCA. The test statistic is distributed as

$T_{N}(x)\xrightarrow{d}\sum_{1\leq l\leq d}B_{l}^{2}(x),~{}~{}0\leq x\leq 1$

where $B_{l}(\cdot)$ are independent Brownian Bridges on $[0,1]$.

###  2.4 Phase changepoint

From the calculation of the Karcher Mean, $\mu_{q}^{k}$, we have the set of
warping functions $\\{\gamma_{i}\\}^{k}$. To perform a changepoint test on the
space of warping functions directly is difficult, as $\Gamma$ is a infinite-
dimensional nonlinear manifold and computation of the mean is not
straightforward due to the unknown geometry. To facilitate computing a mean,
we can simplify the geometry of $\Gamma$.

####  2.4.1 Simplifying Geometry of $\Gamma$

The space of warping functions, $\Gamma$, is an infinite-dimensional nonlinear
manifold and therefore cannot be treated as a standard Hilbert space. To
overcome this problem, we will use tools from differential geometry to perform
statistical analyses and model the warping functions. The following framework
was previously used in various settings including (1) modeling re-
parameterizations of curves (Srivastava and Jermyn, 2009) and (2) putting
prior distributions on warping functions (Kurtek, 2017; Tucker et al., 2021),
and many others. It is also very closely related to the square-root
representation of probability density functions introduced by Bhattacharya
(1943).

We represent an element $\gamma\in\Gamma$ by the square-root of its derivative
$\psi=\sqrt{\dot{\gamma}}$. Note that this is the same as the SRVF defined
earlier, and it takes this form since $\dot{\gamma}>0$. The identity
$\gamma_{id}$ maps to a constant function with value $\psi_{id}(t)=1$. Since
$\gamma(0)=0$, the mapping from $\gamma$ to $\psi$ is a bijection and one can
reconstruct $\gamma$ from $\psi$ using $\gamma(t)=\int_{0}^{t}\psi(s)^{2}ds$.
An important advantage of this transformation is that since
$\|\psi\|^{2}=\int_{0}^{1}\psi(t)^{2}dt=\int_{0}^{1}\dot{\gamma}(t)dt=\gamma(1)-\gamma(0)=1$,
the set of all such $\psi$s is the positive orthant of the unit Hilbert sphere
in $\mathbb{L}^{2}$: $\Psi=\mathbb{S}_{\infty}^{+}$. In other words, the
square-root representation simplifies the complicated geometry of $\Gamma$ to
a (subset of a) unit sphere. The distance between any two warping functions,
i.e., the phase distance, is exactly the arc-length between their
corresponding SRVFs on $\Psi$:

$d_{p}(\gamma_{1},\gamma_{2})=d_{\psi}(\psi_{1},\psi_{2})\equiv\cos^{-1}\left(\int_{0}^{1}\psi_{1}(t)\psi_{2}(t)dt\right)\
.$ (4)

Figure 2 depicts the SRVF space of warping functions as a unit sphere (Tucker,
2014).

$\psi_{1}$$v_{i}$$\psi_{i}$ Figure 2: Depiction of the SRVF space of warping
functions as a sphere and a tangent space at $\psi_{1}$.

While the geometry of $\Psi\subset\mathbb{S}_{\infty}$ is more tractable, it
is still a nonlinear manifold and computing standard statistics remains
difficult. Instead, we use a tangent (vector) space at a certain fixed point
for further analysis. The tangent space at any point $\psi\in\Psi$ is given
by: $T_{\psi}(\Psi)=\\{v\in\mathbb{L}^{2}|\int_{0}^{1}v(t)\psi(t)dt=0\\}$. To
map between the representation space $\Psi$ and tangent spaces, one requires
the exponential and inverse-exponential mappings. The exponential map at a
point $\psi\in\Psi$ denoted by $\exp_{\psi}:T_{\psi}(\Psi)\mapsto\Psi$, is
defined as

$\exp_{\psi}(v)=\cos(\|v\|)\psi+\sin(\|v\|)\frac{v}{\|v\|},$ (5)

where $v\in T_{\psi}(\Psi)$. Thus, $\exp_{\psi}(v)$ maps points from the
tangent space at $\psi$ to the representation space $\Psi$. Similarly, the
inverse-exponential map, denoted by $\exp_{\psi}^{-1}:\Psi\mapsto
T_{\psi}(\Psi)$, is defined as

$\exp_{\psi}^{-1}(\psi_{1})=\frac{\theta}{\sin(\theta)}(\psi_{1}-\cos(\theta)\psi),$
(6)

where $\theta=d_{p}(\gamma_{1},\gamma)$. This mapping takes points from the
representation space to the tangent space at $\psi$.

The tangent space representation $v$ is sometimes referred to as a _shooting
vector_ , as depicted in Figure 2. The remaining question is which tangent
space should be used to represent the warping functions. A sensible point on
$\Psi$ to define the tangent space is at the sample Karcher mean
$\hat{\mu}_{\psi}$ (corresponding to $\hat{\mu}_{\gamma}$) of the given
warping functions. We can compute the mean of the warping functions using
Algorithm 2 from Tucker et al. (2013).

We then can perform a test if we have a changepoint in the phase space (on the
tangent space of $\mathbb{S}_{\infty}$) using the test statistic

$S_{n,k}=\frac{1}{\sqrt{n}}\left(k\bar{v}^{k}-k\bar{v}^{n}\right),$ (7)

which is similar as defined previously, where $\bar{v}^{k}$ is the shooting
vector of the Karcher mean of the warping functions, $\hat{\mu}_{\psi}$,
computed using the first $k$ warping functions. We can also define the test
using the $\mathbb{L}^{2}$ norm similar to Gromenko et al. (2017) as

$\Lambda_{2}=\frac{1}{n}\sum_{k=1}^{n}\|k\bar{v}^{k}-k\bar{v}^{n}\|^{2}.$

Under this test we are testing if there is a phase shift or changepoint. The
definition of $T_{n}$, its distribution, and finding $k^{*}$ are the same as
the amplitude changepoint in Section 2.3. As mentioned above we can also
implement this using horizontal fPCA from Tucker et al. (2013).

####  2.4.2 Horizontal Functional Principal Components

To perform horizontal fPCA we will use the tangent space at $\mu_{\psi}$
(Srivastava and Klassen, 2016) to perform analysis, where $\mu_{\psi}$ is the
mean of the transformed warping functions. In this tangent space we can define
a sample covariance function:

$K_{\psi}=\frac{1}{n-1}\sum_{i=1}^{n}v_{i}v_{i}^{\mathsf{T}}\in^{T\times T}.$

The singular value decomposition (SVD) of
$K_{\psi}=U_{\psi}\Lambda_{\psi}V_{\psi}^{\mathsf{T}}$ provides the estimated
principal components of $\\{\psi_{i}\\}$: the principal directions
$U_{\psi,l}$ and the observed scores $\eta_{i,l}=\left\langle
v_{i},U_{\psi,l}\right\rangle$. This analysis on $\mathbb{S}_{\infty}$ is
similar to the ideas presented in Srivastava et al. (2005) although one can
also use the idea of principal nested sphere for this analysis (Jung et al.,
2012).

We can then use the test defined in Equation 3 to perform the changepoint
test. The distribution to find critical point and confidence interval apply
directly.

##  3 Results

To demonstrate our proposed method, we will compare our elastic method to the
cross-section changepoint methods of Aue et al. (2009) and Aue et al. (2018).
First, we will assess our method’s performance on simulated data where an
amplitude changepoint exists. Next, we will look at a simulated data where a
phase changepoint exists. Finally, we will apply our method to a global
temperature observation dataset from ERA5 (ECMWF) and a global aerosol optical
depth (AOD) dataset from MERRA-2 (Global Modeling and Assimilation Office
(GMAO), 2015) during the time period that contains the eruption of Mt.
Pinatubo in Indonesia in June 1991. For each dataset, we will demonstrate the
use of the changepoint methodology to a location near Mt. Pinatubo; following
this, we will present the results of the changepoint methodology applied to
different locations across the globe.

###  3.1 Simulation Data

Amplitude Changepoint
Consider the following simulated data shown in Figure 3 which presents a set
of Gaussian peaks with random phase warpings. This data was simulated using
the equation $f_{i}(t)=z_{i}e^{-(t-a_{i})^{2}/2}$,
$t\in[-6,6],~{}i=1,2,\dots,75$, where $z_{i}$ is i.i.d.
$\mathcal{N}(1,(0.05)^{2})$ and $a_{i}$ is i.i.d. $\mathcal{N}(0,(1.25)^{2})$.
(Note that although the elastic framework was developed for functions on
$[0,1]$, it can easily be adapted to an arbitrary interval). The randomness of
$z_{i}$ introduces amplitude variability while $a_{i}$ introduces phase
variability in the functional data. The dataset consists of 75 functions with
an amplitude changepoint occurring after the 30th function. The changepoint
was computed such that the mean of $z_{i}$ was changed to 1.5.

Figure 3: Simulated data set of Gaussian curves with an amplitude changepoint
and random phase variation.

Panel (a) in Figure 4 presents the elastic amplitude changepoint method on the
simulated data utilizing the fully functional approach. The elastic method
finds the changepoint at function 30 with a $p$-value of .001. The mean of the
functions before the changepoint is shown in red and the corresponding
functions are shown in pink. The mean of the functions after the changepoint
is shown in blue and the corresponding functions are shown in light blue.
Panel (b) in Figure 4 presents the elastic phase changepoint method on the
simulated data and does not find a changepoint using a significant threshold
of $\alpha=0.05$. Panel (c) in Figure 4 presents the method of Aue et al.
(2018) where a changepoint is found at function 38 with a $p$-value of 0.00.
This is not the correct changepoint and due to the cross-sectional nature of
the data the corresponding means computed before and after the changepoint are
not representative of the underlying data generating mechanism. Both curves
have a much lower amplitude and vastly different shape than the original data.

(a) Elastic Amplitude Changepoint

(b) Elastic Phase Changepoint

(c) Cross-Sectional Changepoint of Aue et al. (2018)

Figure 4: Changepoint results on simulated data with an amplitude changepoint.

Phase Changepoint
Consider the following simulated data shown in Figure 5 consisting of two
peaks with random phase warpings. This dataset is similar to data previously
used in Kneip and Ramsay (2008). The individual functions are given by:
$y_{i}(t)=z_{i,1}e^{-(t-1.5)^{2}/2}+z_{i,2}e^{-(t+1.5)^{2}/2}$,
$t\in[-3,3],~{}i=1,2,\dots,75$, where $z_{i,1}$ and $z_{i,2}$ are i.i.d.
$\mathcal{N}(1,(0.25)^{2})$. Each of these functions is then warped according
to: $\gamma_{i}(t)=6\left(\frac{e^{a_{i}(t+3)/6}-1}{e^{a_{i}}-1}\right)-3$ if
$a_{i}\neq 0$, otherwise $\gamma_{i}=\gamma_{id}$ ($\gamma_{id}(t)=t$ is the
identity warping). Here $a_{i}$ are randomly generated from a uniform
distribution on $-1$ to $1$, and the observed functions are computed using
$x_{i}(t)=y_{i}\circ\gamma_{i}(t)$. The dataset consists of 75 functions with
a phase changepoint occurring after the 30th function. The mean of the warping
functions was shifted from the identity after the 30th function.

Figure 5: Simulated functional data set of two peaks with a phase changepoint
and random phase variation.

Panel (a) in Figure 6 presents the elastic amplitude changepoint method on the
simulated data utilizing the fully functional approach. The elastic amplitude
method does not find a changepoint at any significant value. Panel (b) in
Figure 6 presents the elastic phase changepoint method on the simulated data
and does find a changepoint at function 30 using the significant threshold of
$\alpha=0.05$. The mean of the warping functions before the changepoint are
shown in red and the corresponding functions in pink. The mean of the warping
functions after the changepoint are shown in blue and the corresponding
functions in light blue. Panel (c) in Figure 6 presents the method of Aue et
al. (2018) where a changepoint is found at function 41 with a $p$-value of
0.00. This is not the correct changepoint and due to the cross-sectional
nature of the data the corresponding means computed before and after the
changepoint are not representative the underlying data generating mechanism.
Both curves are have much lower amplitude than the original data.

(a) Elastic Amplitude Changepoint

(b) Elastic Phase Changepoint

(c) Cross-Sectional Changepoint of Aue et al. (2018)

Figure 6: Changepoint results on simulated data with an phase changepoint.

Sensitivity Analysis
Lastly, we examined the ability of the method versus Aue et al. (2009) over
1000 random data sets. Figure 7 presents the results of generating 75
simulated functions and inserting a changepoint (either amplitude or phase)
after function 30. The simulated data set was generated according the
following structure

$f=a_{0}\cos(2\pi t)+b_{0}\sin(2\pi t)+a_{1}\cos(4\pi t)+b_{1}\sin(4\pi t)$

where $a_{i}\sim\mathcal{N}(a^{*},.08)$ and $b_{i}\sim\mathcal{N}(b^{*},.08)$
for $i\in\\{0,1\\}$. The components $a^{*}$ and $b^{*}$ are drawn of a uniform
distribution. For each of the 1000 runs the parameters were sampled from the
corresponding distributions, the data generated and the elastic changepoint
methods utilizing the fPCA approach and the fPCA changepoint method of Aue et
al. (2009). The change occurs when $a^{*}$ and $b^{*}$ change values.

Panel (a) in Figure 7 shows a boxplot of estimated changepoint (using a
significance level of $\alpha=0.05$) over the 1000 simulations for the elastic
amplitude method versus the cross-sectional method of Aue where an amplitude
changepoint occurs. Over the 1000 replicates, the elastic method finds the
changepoint in all but 11 simulations. The cross-sectional method finds a
later changepoint and has a high distribution of the estimated changepoint.
Panel (b) in Figure 7 shows the boxplot estimated changepoint over the 1000
simulation (where a significant $p$-value for $\alpha=0.05$) for the elastic
phase method versus the cross-sectional method of Aue where a phase
changepoint occurs. Over the 1000 replicates, the elastic method finds the
changepoint more often. The cross-sectional method finds the changepoint at
certain times, but has a more variable distribution, and overall the estimated
changepoint mean is incorrect due to the cross-sectional behavior as mentioned
above.

(a) Amplitude

(b) Phase

Figure 7: Detected changepoints over 1000 replicate simulations.

###  3.2 ERA5

We rely on the ERA5 (ECMWF) global temperature dataset which has hourly
temporal coverage from 1979 to present, horizontal resolution of 0.25 degree
grid cells, and a vertical coverage from 1 to 1,000 hPA. ERA5 is the fifth
generation European Centre for Medium-Range Weather Forecasts (ECMWF)
atmospheric reanalysis of the global climate. Reanalysis assimilates
observational data with simulation data to produce a globally consistent
dataset. For this study, we focus on daily averages between 1986 and 1995 at
near-surface level (1,000) hPa across 7.5 degree grid cells. We treat each
year as independent functional observations and apply our principal-component-
based elastic functional changepoint method to identify key shifts in temporal
climate trends from year to year. We first focus on a specific location near
the Mt. Pinatubo eruption of 1991.

Panel (a) in Figure 8 presents temperature profiles smoothed using a box
filter for the years 1986 to 1995 at a location near the Mt. Pinatubo eruption
shown by the blue dot in Panel (b). The profiles are the daily average
temperature for each year. Both the elastic amplitude and elastic phase
changepoint methods were executed on the data. The elastic amplitude
changepoint found a significant changepoint ($\alpha=0.05)$ at the year of
1991. The mean before the changepoint is shown in red and the mean after the
changepoint is shown in blue in Panel (a) of Figure 9. The estimated
functional change in the mean of the SRVFs is shown in Panel (b). The
estimated changepoint is the year of the eruption and we see a cooling in the
later months.

(a)

(b)

Figure 8: ERA5 temperatures at location near Mt. Pinatubo Eruption.

(a) Elastic Amplitude Changepoint

(b)

Figure 9: ERA5 changepoint and estimated amplitude change function.

Figure 10 presents the test statistic for the elastic amplitude test. The
clear maximum of the test statistic is shown at 1991. There is a slight peak
at 1989 or 1994, and if multiple changepoint methods such as in Harris et al.
(2022) are used, a possibility of multiple changepoints can be detected.

Figure 10: ERA5 elastic amplitude changepoint test statistic.

In Figure 11, we present the results of the elastic phase changepoint test,
which also shows a maximum in the test statistic at 1991. This detected
changepoint is significant with a p-value of $0.002$, and one can see clear
separation between the means of the warping functions before and after the
estimated phase change in the data in Panel (b).

Figure 11: (Left) ERA5 elastic phase changepoint test statistic. (Right) ERA5
warping functions before and after the estimated changepoint.

We next repeat this analysis at 173 additional locations across the globe to
provide a more comprehensive analysis. The estimated changepoints for the
amplitude PCA approach are plotted in Figure 12. For the amplitude PCA
approach, each estimated changepoint was significant at the $\alpha=0.05$
level, although this does not correct for multiple testing. In general, there
is mild-to-moderate spatial coherence in the estimated changepoints. We also
plot in Figure 12 a frequency table of the number of times each year was
detected as a changepoint. There appear to be more changepoints estimated
throughout the middle range of the years studied (1990-1992). Since many
variabilities (for example, El Niño) are present in climate and reanalysis
data, filtering out known variabilities may better isolate a potential
observed change around the time or after the eruption of Mt. Pinatubo.

Figure 12: (Left) Estimated ERA5 temperature elastic amplitude changepoints
using the principal component-based test. Points were not plotted if the
estimated p-value was greater than 0.05. (Right) Histogram of detected
changepoints using the elastic amplitude changepoint method.

###  3.3 MERRA-2

We repeat this analysis on another atmospheric variable of interest, aerosol
optical depth (AOD), a unitless quantity describing the amount of light that
aerosols block in the air. Since volcanic activity directly injects aerosols
into the atmosphere, identifying changes related to Mt. Pinatubo’s 1991
eruption may be easier to discern. See, for example, Thompson et al. (2009)
for more information on variability and variables related to Mt. Pinatubo’s
eruption. We use a reanalysis product from NASA called MERRA-2 (Global
Modeling and Assimilation Office (GMAO), 2015). For this study, we focus on
daily averages of total AOD extinction at a wavelength of 550 nanometers
between 1986 and 1997.

In Figure 13, we plot smoothed AOD curves at a location near Mt. Pinatubo.
There appears to be a drastic change in this variable around the time of the
eruption, after which the AOD decreases. In Figure 14, we present results from
the amplitude change test, which detects a significant change at the $0.05$
level in 1991. The estimated change shows an increase in AOD for the second
half of the year (corresponding to the months immediately after the eruption
for 1991), as well as an increase in AOD for the first two months of the year.

(a)

(b)

Figure 13: MERRA-2 aerosol optical depths (AOD) at location near Mt. Pinatubo
Eruption.

(a) Elastic Amplitude Changepoint

(b)

Figure 14: MERRA-2 changepoint and estimated amplitude change function.

As with the ERA5 data, we repeat this analysis at different locations, plotted
in Figure 15. Overall, significant changepoints at the $0.05$ level were
found, and the estimated changepoints are more concentrated around 1991
compared to those for the temperature reanalysis. We found that the AOD
decreased quickly to its pre-Pinatubo behavior in 1992 or 1993 in some
locations, leading to estimated changepoints in these years; therefore,
considering an “epidemic” changepoint as in Aston and Kirch (2012) may more
accurately identify the change of behavior in this variable.

Figure 15: (Left) Estimated aerosol optical depth (AOD) MERRA-2 elastic
amplitude changepoints using the principal component-based test. Points were
not plotted if the estimated p-value was greater than 0.05. (Right) Histogram
of detected changepoints using the elastic amplitude changepoint method.

##  4 Conclusion

We have proposed a new flexible elastic approach to functional changepoint
detection that handles phase and amplitude variability in the data
accordingly. We have demonstrated its advantages over the state-of-the-art
methods in the functional changepoint literature using both simulated and real
datasets. Unlike the current cross-sectional methods, the elastic method
accurately estimates the mean of the underlying data generating mechanism and
also handles phase as a nuisance or as information suggesting a potential
changepoint in the functional data. Correctly estimating the underlying mean
before and after the change can have large impacts on interpretability of the
data and implementation of downstream analysis.

For future work, we look to extend the method to the multiple change-point
problem. This extension will not be trivial as scalability is a problem as
shown in Harris et al. (2022). Multiple changepoints can be detected in
$S_{n,k}$ if using a multiple peak detection scheme, but determining the
$p$-value is not straightforward, and the test statistics will have to be
defined carefully. Additionally, we can extend this from curves to
trajectories that lie on Riemannian manifolds $\mathcal{M}$. In this case, one
has to account for the non-zero curvature of the space and in particular, the
calculation of the gradient in the Fisher Rao metric and computation of the
required statistics.

## Acknowledgment

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. This work was supported by the Laboratory Directed Research and
Development program at Sandia National Laboratories; a multi-mission
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. Department of Energy’s National Nuclear
Security Administration under contract DE-NA0003525.

## References

  * (1) ECMWF reanalysis v5 (ERA5).  (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5.  Accessed: 2022-01-15.
  * Aston and Kirch (2012) Aston, J. A. and C. Kirch (2012).  Detecting and estimating changes in dependent functional data.  Journal of Multivariate Analysis 109, 204–220.
  * Aue et al. (2009) Aue, A., R. Gabrys, L. Horváth, and P. Kokoszka (2009).  Estimation of a change-point in the mean function of functional data.  Journal of Multivariate Analysis 100(10), 2254–2269.
  * Aue et al. (2018) Aue, A., G. Rice, and O. Sonmez (2018).  Detecting and dating structural breaks in functional data without dimension reduction.  Journal of the Royal Statistical Society Series B 80, 509–529.
  * Berkes et al. (2009) Berkes, I., R. Gabrys, L. Horvath, and P. Kokoszka (2009).  Detecting changes in the mean of functional observations.  Journal of the Royal Statistical Society Series B 71, 927–946.
  * Bhattacharya (1943) Bhattacharya, A. (1943).  On a measure of divergence between two statistical populations defined by their probability distributions.  Bull. Calcutta Math. Soc. 35, 99–109.
  * Bosq (2012) Bosq, D. (2012).  Linear processes in function spaces: theory and applications, Volume 149.  Springer Science & Business Media.
  * Fortuna et al. (2020) Fortuna, F., S. A. Gattone, and T. Di Battista (2020).  Functional estimation of diversity profiles.  Environmetrics 31(8), e2645.
  * Global Modeling and Assimilation Office (GMAO) (2015) Global Modeling and Assimilation Office (GMAO) (2015).  tavg1_2d_aer_nx_daily: Merra-2 3d iau state, meteorology instantaneous 3-hourly (p-coord, 0.625x0.5l42).  Accessed 2022-11-09 at doi: 10.5067/KLICLTZ8EM9D.
  * Gromenko et al. (2017) Gromenko, O., P. Kokoszka, and M. Reimherr (2017).  Detection of change in the spatiotemporal mean function.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(1), 29–50.
  * Harris et al. (2022) Harris, T., B. Li, and J. D. Tucker (2022).  Scalable multiple changepoint detection for functional data sequences.  Environmetrics 33(2), e2710.
  * Hörmann and Kokoszka (2012) Hörmann, S. and P. Kokoszka (2012).  Functional time series.  In Handbook of statistics, Volume 30, pp. 157–186. Elsevier.
  * Hörmann et al. (2010) Hörmann, S., P. Kokoszka, et al. (2010).  Weakly dependent functional data.  The Annals of Statistics 38(3), 1845–1884.
  * Hyndman and Booth (2008) Hyndman, R. J. and H. Booth (2008).  Stochastic population forecasts using functional data models for mortality, fertility and migration.  International Journal of Forecasting 24(3), 323–342.
  * Jung et al. (2012) Jung, S., I. L. Dryden, and J. S. Marron (2012).  Analysis of principal nested spheres.  Biometrika 99(3), 551–568.
  * Kneip and Ramsay (2008) Kneip, A. and J. O. Ramsay (2008).  Combining registration and fitting for functional models.  Journal of American Statistical Association 103(483), 1155–1165.
  * Kokoszka and Zhang (2012) Kokoszka, P. and X. Zhang (2012).  Functional prediction of intraday cumulative returns.  Statistical Modelling 12(4), 377–398.
  * Kurtek (2017) Kurtek, S. (2017).  A geometric approach to pairwise Bayesian alignment of functional data using importance sampling.  Electronic Journal of Statistics 11(1), 502–531.
  * Lee and Jung (2017) Lee, S. and S. Jung (2017).  Combined analysis of amplitude and phase variations in functional data.  arXiv:1603.01775 [stat.ME], 1–21.
  * Shang and Hyndman (2011) Shang, H. L. and R. J. Hyndman (2011).  Nonparametric time series forecasting with dynamic updating.  Mathematics and Computers in Simulation 81(7), 1310–1324.
  * Sharipov et al. (2016) Sharipov, O., J. Tewes, and M. Wendler (2016).  Sequential block bootstrap in a hilbert space with application to change point analysis.  Canadian Journal of Statistics 44(3), 300–322.
  * Srivastava and Jermyn (2009) Srivastava, A. and I. H. Jermyn (2009).  Looking for shapes in two-dimensional, cluttered point clouds.  IEEE Trans. Pattern Analysis and Machine Intelligence 31(9), 1616–1629.
  * Srivastava et al. (2005) Srivastava, A., S. H. Joshi, W. Mio, and X. Liu (2005).  Statistical shape analysis: Clustering, learning and testing.  IEEE Trans. Pattern Analysis and Machine Intelligence 27(4), 590–602.
  * Srivastava and Klassen (2016) Srivastava, A. and E. P. Klassen (2016).  Functional and shape data analysis.  Springer.
  * Thompson et al. (2009) Thompson, D., J. Wallace, P. Jones, and J. Kennedy (2009, 11).  Identifying signatures of natural climate variability in time series of global-mean surface temperature: Methodology and insights.  Journal of Climate - J CLIMATE 22, 6120–6141.
  * Tucker (2014) Tucker, J. D. (2014).  Functional component analysis and regression using elastic methods.  Ph. D. thesis, Florida State University.
  * Tucker et al. (2021) Tucker, J. D., L. Shand, and K. Chowdhary (2021).  Multimodal Bayesian registration of noisy functions using Hamilitonian Monte Carlo.  Computational Statistics and Data Analysis 163, 1–24.
  * Tucker et al. (2013) Tucker, J. D., W. Wu, and A. Srivastava (2013).  Generative models for functional data using phase and amplitude separation.  Computational Statistics & Data Analysis 61, 50–66.
  * Woodall (2007) Woodall, W. H. (2007).  Current research on profile monitoring.  Production 17(3), 420–425.
  * Zhang et al. (2011) Zhang, X., X. Shao, K. Hayhoe, D. J. Wuebbles, et al. (2011).  Testing the structural stability of temporally dependent functional observations and application to climate projections.  Electronic Journal of Statistics 5, 1765–1796.