[1]AnnaDe Magistris

These authors contributed equally to this work. [1]Elvira Romano These authors
contributed equally to this work.

These authors contributed equally to this work.

[1]Dipartimento di Matematica e Fisica, University of Camapania ”Luigi
Vanvitelli”, Viale Abramo Lincoln 5, Caserta, 81100, Italy , Caserta

# Roughness regularization for functional data analysis with free knots spline
estimation

<EMAIL_ADDRESS>ValentinaDe Simone
<EMAIL_ADDRESS><EMAIL_ADDRESS>Gerardo Toraldo
<EMAIL_ADDRESS>*

###### Abstract

In the era of big data, an ever-growing volume of information is recorded,
either continuously over time or sporadically, at distinct time intervals.
Functional Data Analysis (FDA) stands at the cutting edge of this data
revolution, offering a powerful framework for handling and extracting
meaningful insights from such complex datasets. The currently proposed FDA
methods can often encounter challenges, especially when dealing with curves of
varying shapes. This can largely be attributed to the method’s strong
dependence on data approximation as a key aspect of the analysis process. In
this work, we propose a free knots spline estimation method for functional
data with two penalty terms and demonstrate its performance by comparing the
results of several clustering methods on simulated and real data.

###### keywords:

functional data analysis, free knots splines, double penalty terms,
clustering.

##  1 Introduction

In recent years, there has been a rapid increase in the availability of large
data resulting from the observation of a phenomenon along a continuous domain,
such as time, space, or frequency domains. The growing need to analyze data of
this type, characterized by an intrinsic functional nature, is the basis for
the development of the field of Functional Data Analysis (FDA) [1],[2].
Functional data analysis typically involves two key steps: first, representing
discrete observations as functions and then applying statistical methods. The
step of functional representation helps to reduce noise and fluctuations in
the data and plays a key rule in the FDA because it sets the foundation for
subsequent statistical analyses. There are several methods and techniques that
can be employed to perform functional representation. Some common approaches
include smoothing spline, wavelet analysis, Fourier analysis, principal
component analysis, smoothing techniques, and interpolation methods. The
choice of the functional representation method depends on the nature of the
data, the specific goals of the analysis, and the assumptions made about the
underlying structure of the functions. Each method has its advantages and
limitations, and the selection should be based on the characteristics of the
data and the objectives of the analysis. To address the efficient
representation and analysis of functional data with various shapes, we focus
on the problem of representing functional data using free knots spline
estimation. In the last years, the problem of representing functional data
with free knots spline has been addressed by the same authors. Gervini [3]
proposes free knots regression spline estimators for the mean and the variance
components of a sample of curves, showing that free knots splines estimate
salient features of the functions (such as sharp peaks) more accurately than
smoothing splines. Inspired by a Bayesian model [4] propose estimating various
curves using piecewise polynomials. The first one establishes a joint
distribution over the number and positions of knots defining the piecewise
polynomials and employs reversible jump Markov Chain Monte Carlo methods for
posterior computation. The second, based on a marginalised chain on the knot
number and locations, provides a method for inference about the regression
coefficients and functions of them in both normal and non normal models. Among
others, a novel knot selection algorithm is introduced by [1] for
nonparametric regression using regression splines. Recently [5], building upon
a previous proposal of [6], introduces a data-driven approach that uses a
machine learning-style algorithm to select knots efficiently and employs
orthogonal spline bases (splinets) to represent sparse functional data. In
this work, generalizing ’the balanced-discrepancy principle’ proposed by [7]
in the functional data analysis framework, we propose a simple double
penalization criterion to improve the smoothing process with free knots spline
approximation for various types of curve shapes. This approach is designed to
improve the precision and overall effectiveness of fitting spline curves to
diverse datasets. The paper is organized as follows. Section $2$ introduces
mathematical foundations of smoothing with roughness penalty. In section $3$,
we illustrate free knots spline estimation and their construction with a
double penalization criterion. In section $4$ we illustrate performance of the
proposed penalization criterion on simulated data. An extensive comparative
analysis to evaluate which approximation method is most suitable for
functional data in terms of clustering performance is conducted. Section $5$
shows the performance of a clustering method on COVID-19 pandemic data for 30
countries by comparing free knots spline and free knots spline with double
penalization criterion.

##  2 Smoothing Functional Data

Let $\\{x_{i}(t),\;i=1,...,n\\}$ be a sample of real-valued functions on a
compact interval $T$ related to a functional variable $X$. The sample curves
can be considered realizations of a stochastic process $X=\\{X(t):t\in T\\}$
whose sample functions belong to the Hilbert space $L^{2}(T)$ of square
integrable functions with the usual inner product $f,g=\int_{T}f(s)g(s)ds$.

In real applications, $x_{i}(t)$ often cannot be observed directly, but may be
collected discretely over time points $\\{t_{1},...,t_{h}\\}\subset T$. FDA
aims to reconstruct the true functions from the observed data using the basis
function approach, that reduces the estimation of the function $x_{i}(t)$ to
the estimation of a finite linear combination of known basis functions
$\phi_{k}(t),\;k=1\ldots n_{\mathcal{B}}$

${x}_{j}(t)=\sum_{k=1}^{n_{\mathcal{B}}}\textbf{c}_{k}^{j}{\phi}_{k}({t}),\>\>\>j=1,...,n$

where $\mathbf{c}^{j}\in\mathbb{R}^{n_{\mathcal{B}}}$ is a vector of
coefficients. The choice of the base and the dimension $n_{\mathcal{B}}$
depend strongly on the nature of the data. The most known bases are: the
Fourier base, B-spline, polynomial, exponential, wavelets. The Fourier base is
used when the data has a cyclical nature, the exponential base when the data
shows exponential growth, and the B-spline base are the most used when the
data does not have a strong cyclical trend. There are different ways of
obtaining the basis coefficients depending on the kind of observations.
Generally observed data are contaminated by random noise that can be viewed as
random fluctuations around a smooth trajectory, or as actual errors in the
measurements.

The observed data is collected in a matrix $\textbf{Y}\in\mathbb{R}^{h\times
n}$ whose elements are:

$y_{i,j}={x}_{j}({t}_{i})+\bm{\epsilon}_{i}^{j},\;\;\;i=1,...,h;\;j=1,...,n$

where $\bm{\epsilon}_{i}$ is unobserved error term independent and identically
distributed random variable with zero mean and constant variance $\sigma^{2}$.
In this case an appropriate way to estimate the basis coefficients from the
data is by using a least squares approximation

$\hat{\textbf{C}}=\underset{{\textbf{C}}}{\operatorname{argmin}}\left\|\textbf{Y}-\bm{\Phi}^{T}\textbf{C}\right\|^{2}_{2},$

with $\bm{\Phi}\in\mathbb{R}^{{n_{\mathcal{B}}\times h}}$,
$\bm{C}\in\mathbb{R}^{{n_{\mathcal{B}}}\times n}$, where
$\phi_{ij}=\bm{\phi}_{k}(\mathbf{t}_{j})$ and $c_{ij}=\mathbf{c}_{j}^{i}$.

The smoothness is implicitly controlled by the number of basis functions,
${n_{\mathcal{B}}}$. If we assume that $n_{\mathcal{B}}\leq h$ and
$rank(\bm{\Phi})=n_{\mathcal{B}}$ then

$\hat{\textbf{C}}={(\bm{\Phi}\bm{\Phi}^{T})}^{-1}\bm{\Phi}\textbf{Y}.$ (1)

and the matrix $\hat{\textbf{Y}}$ of the approximation values is:

$\hat{\textbf{Y}}=\bm{\Phi}^{T}\hat{\textbf{C}}=\bm{\Phi}^{T}{(\bm{\Phi}\bm{\Phi}^{T})}^{-1}\bm{\Phi}\textbf{Y}.$
(2)

Increasing ${n_{\mathcal{B}}}$ leads to overfitting and generates a $x$ curve
that is overly ”wiggly. One way to overcome this drawback is to use a
roughness penalty term.

###  2.1 Smoothing functional data with a roughness penalty

The roughness penalty approach defines a measure of the roughness of the
fitted function $x$ using the derivatives of some order $m\geq 1$

$PEN_{m}(x)=\int{[D^{m}x(s)]}^{2}ds.$ (3)

This allows to measure the closeness of $x(t)$ to a polynomial of order $m$.
The most commonly used roughness penalty is $m=2$ that permits to keep under
control the curvature of the curve, that is we control the variability of the
slope of the curve [8].

Often, there is a need for a wider class of measures of deviation. Especially,
when there is periodicity in the data or an exponential trend, it would not be
sufficient to use the integrated squared mth derivative because it can only
penalize deviations from polynomials. More generally, a measure of roughness
is then given by

$PEN_{L}(x)=\int{[Lx(s)]}^{2}ds,$ (4)

where $L$ is the linear differential operator defined as [2]

$\displaystyle Lx(t)=\alpha_{0}(t)x(t)+\alpha_{1}(t)D^{1}x(t)+....+$
$\displaystyle+\alpha_{m-1}(t)D^{m-1}x(t)+D^{m}x(t).$

Obviously $PEN_{m}(x)$ is a special case of $PEN_{L}(x)$ with
$\alpha_{0}(t)=\alpha_{1}(t)=...=\alpha_{m-1}(t)=0$ and $\alpha_{m}(t)=1$. In
the following we will assume the $\alpha_{i}(t)$ to be constant. Then, the
penalized least squares approximation is given by

$\hat{\textbf{C}}=\underset{{\textbf{C}}}{\operatorname{argmin}}\left\|\textbf{Y}-\bm{\Phi}^{T}\textbf{C}\right\|^{2}_{2}+\lambda\textbf{C}^{T}\mathbf{R}_{L}\textbf{C},$
(5)

with
$\mathbf{R}_{L}=\alpha_{0}\mathbf{R}_{0}+\alpha_{1}\mathbf{R}_{1}+...\alpha_{m}\mathbf{R}_{m}$
discretization of $PEN_{L}(x)$, where
$\mathbf{R}_{l}\in\mathbb{R}^{{n_{\mathcal{B}}}\times{n_{\mathcal{B}}}}$
($0\leq l\leq m)$ :

$(r_{l})_{ij}=\int_{I}D^{l}\phi_{i}(s)D^{l}\phi_{j}(s)^{T}ds,$ (6)

where $I$ is a suitable interval containing the data. The regularization
parameter $\lambda>0$ manages the compromise between the fitting to the data
and the roughness of the function: the smaller it is $\lambda$, the closer the
estimate is to the estimate of least squares and tends to interpolate the
observed points; the higher is $\lambda$, the flatter the smooth function
tends to be. This parameter can be calibrated through a generalised cross
validation (GCV)[2]. In this paper we consider the first and second-order
roughness in order to control both the course and the variability of slope of
the curve and to avoid placing an excessive burden on the cost of the penalty
matrix. Then the penalized least square problem becomes

$\displaystyle\hat{\textbf{C}}=\underset{{\textbf{C}}}{\operatorname{argmin}}\left\|\textbf{Y}-\bm{\Phi}^{T}\textbf{C}\right\|^{2}_{2}+\lambda_{1}\textbf{C}^{T}\mathbf{R}_{1}\textbf{C}+$
(7) $\displaystyle+\lambda_{2}\textbf{C}^{T}\mathbf{R}_{2}\textbf{C}\ \ \ \ \
\ \ \lambda_{2}>0,\lambda_{1}>0.$

The solution of (7) can be computed by solving the system arising by first
order optimality conditions

$\hat{\textbf{C}}={(\bm{\Phi}\bm{\Phi}^{T}+\lambda_{1}\textbf{R}_{1}+\lambda_{2}\textbf{R}_{2})}^{-1}\bm{\Phi}\textbf{Y}.$
(8)

The expression for the data-fitting $\hat{\textbf{Y}}$ is:

$\hat{\textbf{Y}}=\bm{\Phi}^{T}{(\bm{\Phi}\bm{\Phi}^{T}+\lambda_{1}\textbf{R}_{1}+\lambda_{2}\textbf{R}_{2})}^{-1}\bm{\Phi}\textbf{Y}$
(9)

##  3 Free knots spline

For their compact support and fast computation, as well as the ability to
create smooth approximations of non periodic data, B-splines are a common
choice in the functional data framework. The approximations by splines can be
significantly improved if knots are allowed to be free rather than at a priori
fixed locations [3].

It is well-known that the primary advantage of free knots spline over
smoothing splines is their greater flexibility in modelling data, as they
allow for a better adaptation of the curve’s shape to the specific
characteristics of the data.

The free knots spline is a spline in which the positions of the knots are
considered parameters to be estimated by the data. Adjusting the position of
the nodes allows you to adapt the shape of the function to the target
function.

This section considers smooth estimators for knots selection for the
approximation of a given curve. Existing approaches are based on individual
levelling of sample curves, followed by the mean of the cross-section and the
calculation of covariance [9]. These methods, however, do not take strength
from the dataset that we have available in the leveling phase. A further
drawback is that analytical expressions for optimal knots locations, or even
for the general characteristics of optimal knots distributions, are not easy
to derive. We introduce free knots spline estimators that avoid individual
levelling. We show that this approach applied to the methods seen in the
previous section (smoothing spline with one parameter and with two parameters)
often produces better estimators than sanding splines [2] in which knots are
chosen randomly and equally at the cost of a modest increase in computational
complexity. In this section we introduce the algorithm of optimal knots
placement. Given a vector of nodes $\bm{\tau}\in\Re^{p}$ with
$a<\tau_{1}<\tau_{2}..<\tau_{p}<b$, the Jupp[10] transformation k of
$\bm{\tau}$, $\textbf{k}=J(\bm{\tau})$, is defined componentwise as

${k}_{i}=\log\frac{(\bm{\tau}_{i+1}-\bm{\tau}_{i})}{(\bm{\tau}_{i}-\bm{\tau}_{i-1})}\>\>\>\>\>\>\>\>\>i=1,...,p$

where ${\tau}_{0}=a$ and ${\tau}_{p+1}=b$. This one-to-one transformation maps
constrained and ascending knots vectors $\bm{\tau}$ on unconstrained and
unclassified carriers k, which has some practical and theoretical advantages.
Note that for each fixed set of knots, the class of such splines is a linear
space of functions with $(p+r)={n_{\mathcal{B}}}$ free parameters where $r$ is
the order of spline and $p$ is the length of knots.

Let $\bm{\phi}(\textbf{t},\textbf{k})\in\mathbb{R}^{{n_{\mathcal{B}}}}$ be the
vector of the basic functions B-spline corresponding to a set of nodes k and
$\bm{\Phi}(\textbf{k})$ the matrix $h\times{n_{\mathcal{B}}}$ whose j-row is
$\bm{\phi}(\textbf{t}_{j},\textbf{k})^{T}$. We find the coefficients of linear
expansion and the vector of the optimal nodes minimizing the penalized least
squares problem111Note that for $\lambda_{1}=\lambda_{2}=0$ (10) becomes the
classic free knots spline problem.:

$\displaystyle(\hat{\textbf{C}},\hat{\textbf{k}})=\underset{({\textbf{C}},{\textbf{k}})\in\mathbb{R}^{n_{\mathcal{B}}}\times\mathbb{R}^{p}}{\operatorname{argmin}}\left\|\textbf{Y}-\bm{\Phi}^{T}\textbf{C}\right\|^{2}_{2}+\lambda_{1}\textbf{C}^{T}\textbf{R}_{1}\textbf{C}+$
(10) $\displaystyle+\lambda_{2}\textbf{C}^{T}\textbf{R}_{2}\textbf{C}\ \ \ \ \
\lambda_{1}\geq 0,\lambda_{2}\geq 0$

Alternating minimization algorithm can be used to solve problem (10); at each
iteration the problem is splitted in two optimization ones. Closed form
solution can be obtained for minimization with respect to C fixing k. The
optimal value can be obtained by solving the system with coefficient matrix

$\mathbf{H}(\mathbf{k})=\bm{\Phi}(\textbf{k})\bm{\Phi}(\textbf{k})^{T}+\lambda_{1}\textbf{R}_{1}+\lambda_{2}\textbf{R}_{2}.$

######  Proposition 1.

If the $\bm{\Phi}(\textbf{k})$ is full rank, the matrix H(k) is SPD and let
$\sigma(\textbf{H})$ be the set of eigenvalues of matrix H(k) we have
$\sigma(\textbf{H})\subseteq[\sigma^{-},\sigma^{+}]$, with

$\sigma^{-}=\sigma_{min}(\bm{\Phi}(\textbf{k})\bm{\Phi}(\textbf{k})^{T})+\lambda_{1}\sigma_{min}(\textbf{R}_{1})+\lambda_{2}\sigma_{min}(\textbf{R}_{2})$
$\sigma^{+}=\sigma_{max}(\bm{\Phi}(\textbf{k})\bm{\Phi}(\textbf{k})^{T})+\lambda_{1}\sigma_{max}(\textbf{R}_{1})+\lambda_{2}\sigma_{max}(\textbf{R}_{2})$

.

######  Proof 3.2.

H(k) is the sum of $\bm{\Phi}(\textbf{k})\bm{\Phi}(\textbf{k})^{T}$ that is
SPD, and of other matrices that are non negative definite [8]; then the first
statement follows. For the second statement we recall that if A and B are
real, symmetric matrices, then A + B has real eigenvalues, and the following
inequalities hold

$\displaystyle\sigma_{min}(\textbf{A})+\sigma_{min}(\textbf{B})\leq\sigma_{min}(\textbf{A
+ B})\leq$ $\displaystyle\leq\sigma_{max}(\textbf{A +
B})\leq\sigma_{max}(\textbf{A})+\sigma_{max}(\textbf{B}).$

Proposition 1 suggests to use Cholesky factorization to compute

$\hat{\textbf{C}}(\textbf{k})=\mathbf{H}(\mathbf{k})^{-1}\bm{\Phi}(\textbf{k})\textbf{Y}.$

The minimization with respect to k fixing C gives the following nonlinear
optimization problem:

$\hat{\textbf{k}}=\underset{\textbf{k}\in\mathbb{R}^{p}}{\operatorname{argmin}}\left\|\textbf{Y}-\bm{\Phi}(\textbf{k})^{T}\textbf{H}(\textbf{k})^{-1}\bm{\Phi}(\textbf{k})\textbf{Y}\right\|^{2}$
(11)

To solve (11), we apply the knots addition algorithm that produces knots
sequences of increasing dimensions. We define the functional
$f_{j}:\mathbb{R}^{j}\rightarrow\mathbb{R}$ as follows:

$f_{j}(\textbf{k})=\left\|\textbf{Y}-\bm{\Phi}(\textbf{k})^{T}\textbf{H}(\textbf{k})^{-1}\bm{\Phi}(\textbf{k})\textbf{Y}\right\|^{2}\>\>\>1\leq
j\leq p.$ (12)

Now, we will present the procedure of the gradual node addition algorithm[3].

Gradual node addition algorithm Initialization Choose an ordered grid
$F_{1}=\\{s_{1}^{1},...,s_{N}^{1}\\}\subset(a,b)$. Compute
$J_{1}=\\{J(\\{s_{1}^{1}\\}),...,\ J(\\{s_{N}^{1}\\})\\}$. Find
$\tilde{\textbf{k}}_{1}=\underset{\textbf{k}\in
J_{1}}{\operatorname{argmin}}f_{1}(\textbf{k})$ Compute $\hat{\textbf{k}}_{1}$
solution of (11) with $p=1$ using the Gauss-Newton algorithm with
$\tilde{\textbf{k}}_{1}$ as the starting point.
$\hat{\bm{\tau}}_{1}=J^{-1}(\hat{\textbf{k}}_{1})$ Forward addition For
$i=2,...,p$ Choose an ordered grid
$F_{i}=\\{s_{1}^{i},...,s_{N}^{i}\\}\subset(a,b)$. Compute
$J_{i}=\\{J\bigl{(}\hat{{\bm{\tau}}}_{j-1}\bigcup\\{s_{1}^{i}\\}\bigr{)},...,\bigl{(}\hat{{\bm{\tau}}}_{j-1}\bigcup\\{s_{N}^{i}\\}\bigr{)}\\}$.
Find $\tilde{\textbf{k}}_{i}=\underset{\textbf{k}\in
J_{i}}{\operatorname{argmin}}f_{i}(\textbf{k})$

Compute $\hat{\textbf{k}}_{i}$ solution of (11) with $p=i$ using the Gauss-
Newton algorithm with $\tilde{\textbf{k}}_{i}$ as the starting point.
$\hat{\bm{\tau}}_{i}=J^{-1}(\hat{\textbf{k}}_{i})$ . EndFor

Although there is no guarantee that this (or any other) algorithm will find
the global minimizer of (11), we have found that it works well in practice. In
our simulations and examples, knots have been added in the ”right” order
[3],[10]. This is important for the selection of the model, since the optimal
number of knots $p$ is never known in practice and will be chosen on the basis
of sequences of intermediate nodes.

##  4 Computational experiments

In this section, we present some computational results using our algorithm on
functional data coming both on synthetic data and on data from a real-world
application. More precisely, we present applications of three different
clustering methods to evaluate the benefits of detecting clusters when using
free knots splines estimation with two penalty terms with respect to free
knots splines and free knots splines with one penalty term. To facilitate the
notation, we will use FS0 to indicate the traditional free knots spline
approximation, FS1 to indicate the free knots spline approximation with a
single penalty term on the second derivative and FS2 to indicate the free
knots spline approximation with a double penalty term.

Figure 1: Simulated data.

We will consider the classical k-means method for functional data [11], a
model based clustering method [12] and a hierarchical agglomerative clustering
methods [13],[14], which we will refer to respectively R packages as
kmeans.fd, funFEM, fdahclust.

###  4.1 Synthetic datasets

In the simulated scenario, four groups of 50 functions each were generated as
in [15] using the following functions as average functions

  * •

$y(t)=-2\sin(t-1)\log(t+0.5)$,

  * •

$y(t)=2\cos(t)\log(t+0.5)$,

  * •

$y(t)=-0.5-0.2\cos(0.5(t-1))t^{1.5}\sqrt{5\sqrt{t}+0.5}$,

  * •

$y(t)=1.2\cos(t)\log(t+0.5)\sqrt{t+0.5}$,

i.e. functions that resemble a negative sine (or cosine) wave with a regular
oscillation period, multiplied either by the natural logarithm of $(t+0.5)$
which modulates the amplitude of the curve, causing it to progressively
decrease as t increases, or the product of a combination of log, power and
square root functions which makes the curve strongly nonlinear and intricate.

Errors determined by a normal distribution with zero mean and standard
deviation are added to each curve. Random errors are randomly generated and
modeled to incorporate heterogeneity in the variance along the functional
curves. Each function was sampled in a common set of 50 randomly chosen points
in the interval $[0,5]$.

Each group, as can be seen in Figure 1, is characterized by differences in
terms of amplitude, shape, and complexity.

The first step in approximating the data was to determine the regularization
parameters through cross validation: we chose the values $\lambda_{1}=10^{-7}$
and $\lambda_{2}=10^{-5}$ for FS2 obtained by minimizing with respect to
$\lambda_{1}$ and $\lambda_{2}$ the GCV on a $L\times L$ grid, with
$L=\\{10^{l},\;l=-8,-6\ldots 3,4\\}$ and $\lambda_{2}=10^{-5}$ for FS1.

The number of basis for the approximation, selected via cross validation, is
set to 12. This simulation compares FS0, FS1 and FS2. The graphical results
are shown in Figure 2, and the numerical results are presented in Table 1.

Performances of the two different approximation methods were measured using
the classic Integrated Sum of Square Errors (ISSE) and its local version
defined on the tails of the curves. The ISSE is a statistical metric used to
assess the goodness of fit of a regression model or interpolation model to
observed data. It is particularly useful in scenarios where the dependent
variable is a continuous function of an independent variable, such as in
functional data analysis. The ISSE is calculated by summing the squared
differences between the observed values and the predicted values over the
entire data range.

Figure 2: Comparison of FS0,FS1 and FS2 for data approximation in four case.

Similarly to the traditional ISSE, we define a Local ISSE (ISSEinf,ISSEsup)
with the aim to quantify the fit in the initial and final portions of the
curves.

The expression of ISSE remains the same, but the boundaries of the integration
interval need to be adjusted to reflect the specific regions of interest. This
allows you to focus on the initial and final tails of the curve rather than
the entire curve. The choice of regions for calculating the Integrated Sum of
Squared Errors depends on the analysis objectives and the nature of the
functional data or curves under examination. In our context, cross validation
was used to test the model on different parts of the curves and identify the
regions in the tails. Degrees of freedom ($df$), Integrated Sum of Square
Errors, and the GCV scores are used to evaluate the quality of both the
overall and local model fit. The results in Table 1 show that the lowest ISSE
values were obtained using the free knots splines incorporating two
regularization terms, both on the entire interval and on the tails.

Table 1: Table of comparison between FS0, FS1 and FS2. | FS0 | FS1 | FS2
---|---|---|---
$df$ | $0.120e+2$ | $0.119e+2$ | $0.119e+2$
$ISSE$ | $0.154e+0$ | $0.154e+0$ | $0.966e-1$
$ISSE_{inf}$ | $0.977e-1$ | $0.977e-1$ | $0.667e-1$
$ISSE_{sup}$ | $0.237e-1$ | $0.240e-1$ | $0.227e-1$
$GCV$ | $0.506e+1$ | $0.509e+1$ | $0.198e+1$

The advantages obtained by using the dual penalty approach compared to free
knots spline are even more evident if we look at the clustering.

Figure 3: Evaluation of the number of clusters in the simulation dataset. The
optimal number is indicated from the dotted line.

To determine the number of clusters, we used the Elbow method. In cluster
analysis, the Elbow method is a heuristic used to identify the optimal number
of clusters in a given dataset. This method involves plotting the explained
variation as a function of the number of clusters and selecting the ’elbow
point’ on the curve as the optimal number of clusters. The results in Figure 3
show that according to the index values, the optimal cluster number is c = 4.
By applying kmeans.fd clustering to functional data, we can ascertain the
number of curves assigned to each cluster for the methods under consideration.

Table 2: Table of comparison between clustering methods with FS0, FS1 and FS2: cluster cardinality and false positives ($FP$). | FS0 | FP | FS1 | FP | FS2 | FP
---|---|---|---|---|---|---
cluster 1 | 34 | 0 | 34 | 0 | 39 | 0
cluster 2 | 48 | 1 | 48 | 1 | 49 | 0
cluster 3 | 66 | 17 | 66 | 17 | 61 | 11
cluster 4 | 52 | 2 | 52 | 2 | 51 | 1
Figure 4: Clustering structure: Cluster 1 in red, Cluster 2 in blue, Cluster 3
in green, and Cluster 4 in orange. First column represents synthetic
datasets,the second column represents the data approximated using FS0
segmented into clusters, the third column shows the data approximated using
FS1 segmented into clusters, the last column illustrates the data approximated
using FS2 segmented into clusters.

Results reported in Table 2 and in Figure 4, with respect to FS2, many more
misclassifications are observed for both FS0 and FS1.

The approximation that yields the best clusters is the one using FS2.

In conclusion, it can be seen that introducing two regularization terms into
the free knots spline approximation method led to have better approximation
results also for further analyses such as clustering. Indeed, when data
exhibit very similar shapes, clustering typically works regardless of
approximation; the problem arises when dealing with data of highly dissimilar
shapes, as seen in this case with simulated data, which a conventional method
fails to capture.

To consolidate the results obtained from clustering the curves approximated by
the free knot spline, we demonstrate that no method of functional data
clustering can yield the same or better results. From now on we no longer
consider $FS1$ since just adding a regularization term does not improve the
analysis through the use of clustering.

Table 3: Table of comparison between clustering methods with FS0 and FS2: cluster cardinality and false positives ($FP$). Method | FS0 | FP | FS2 | FP
---|---|---|---|---
$kmeans.fd$ | cluster 1 = 34 | 0 | cluster 1 = 39 | 0
| cluster 2 = 48 | 1 | cluster 2 = 49 | 0
| cluster 3 = 66 | 17 | cluster 3 = 61 | 11
| cluster 4 = 52 | 2 | cluster 4 = 51 | 1
$funFEM$ | cluster 1 = 56 | 41 | cluster 1 = 26 | 0
| cluster 2 = 53 | 38 | cluster 2 = 17 | 0
| cluster 3 = 35 | 40 | cluster 3 = 100 | 50
| cluster 4 = 56 | 43 | cluster 4 = 57 | 7
$fdahclust$ | cluster 1 = 69 | 19 | cluster 1 = 62 | 12
| cluster 2 = 16 | 0 | cluster 2 = 8 | 0
| cluster 3 = 111 | 61 | cluster 3 = 87 | 37
| cluster 4 = 4 | 24 | cluster 4 = 43 | 0

Results resented in Table 3 for all the clustering methods show how the best
classification is achieved by applying clustering to curves approximated with
$FS2$. Nevertheless, none of the classifications results are comparable to the
classification obtained by the kmeans.fd method.

We can draw the same conclusions by evaluating the results obtained using the
previous clustering methods by performing $150$ simulations. In order to
compare the performances, we use the Adjusted Rand index $ARI$ measuring the
percentage of times a clustering algorithm correctly determines the number of
clusters out of a total of 150 simulations. Clustering can be thought of as a
series of decisions where the goal is to group two individuals into the same
cluster if and only if they are similar. A ”true positive” ($TP$) decision
correctly assigns two similar individuals to the same cluster, while a ”true
negative” ($TN$) decision correctly assigns two dissimilar individuals to
different clusters. On the other hand, a ”false positive” ($FP$) decision
incorrectly assigns two dissimilar individuals to the same cluster, and a
”false negative” ($FN$) decision incorrectly assigns two similar individuals
to different clusters.

The Rand index (RI) is defined as:

$RI=\frac{TP+TN}{TP+FP+FN+TN}$

and falls within the range between 0 and 1. A value of 0 indicates that the
two data clusterings do not agree on any pair of points, while a value of 1
means that the data clusterings are identical. However, RI’s value may not be
close to 0 when category labels are randomly assigned, leading to potential
issues.

Table 4: A comparison table of clustering methods, including FS0 and FS2, based on 150 simulations.  Method | ARI FS0 | ARI FS2
---|---|---
$kmeans.fd$ | 0.78 | 0.83
$funFEM$ | 0.54 | 0.65
$fdahclust$ | 0.52 | 0.64

To address this problem, the Adjusted Rand index (ARI) is introduced, defined
as:

$ARI=\frac{RI-E[RI]}{\max(RI)-E[RI]}.$

ARI’s range is between $-1$ and $1$, with values closer to $1$ indicating
better clustering results. From Table 4 we can observe that we can observe
that for all the clustering methods, the use in the functional approximation
of $FS2$ leads to higher ARI compared $FS0$.

###  4.2 Application: Real Datasets

In this section, we aim to provide further confirmation of the validity and
competitiveness of the proposed method by examining a real-world dataset.

We will examine the application of clustering to the COVID-19 dataset on “New
cases in different countries from 2020 to 2021, specifically COVID-19 pandemic
data for 30 countries. In this context, our goal is to demonstrate the
feasibility of analyzing pandemic models through time series clustering using
knot-free splines with two regularization terms for approximating functional
data and finding taxonomies within the data. The data used in this paper were
sourced from the World Health Organization (WHO) and the COVID-19 Data
Hub[15],[16]. These databases are known for their transparency, open
accessibility, and high credibility, ensuring a high level of accuracy.

As a significant number of countries experienced substantial outbreaks in
March 2020, we selected this month as the starting point for our sequential
data in this research. The endpoint of our dataset corresponds to the date of
November 30st, 2021. As in [17], to mitigate the impact of factors such as
varying population size, land area, population density, and population
mobility across different countries, which can lead to significant differences
in the magnitude of daily new COVID-19 cases, we at first proceeded to
standardize the raw data as follows:

$y_{i,t}^{*}=\frac{y_{it}-\bar{y}_{i}}{s_{i}}\;\;\;\;t=1,2,...,T\;\;\;\;i=1,2,...,N$

where $y_{i,t}^{*}$ represents the normalized number of daily new cases in the
country $i$ on day $t$, $y_{it}$ stands for the number of daily new cases in
the country $i$ on day $t$, $\bar{y}_{i}$ represents the mean of daily new
cases in the country $i$ during the whole study period, $s_{i}$ is the
standard deviation of daily new cases in the country $i$ during the whole
study period.

Figure 5: New Cases of Covid19 in different Country from 2020 to 2021.

The data obtained can be seen in Figure 5. The next step was to take the
values $\lambda_{1}=10^{-3}$ and $\lambda_{2}=10^{-7}$. for the regularization
parameters, obtained through a general cross-validation on a pre-established
grid.

Table 5: Numerical results. | Free knots spline | Free knots spline
---|---|---
| $\lambda_{1}=\lambda_{2}=0$ | $\lambda_{1}>0,\lambda_{2}>0$
$df$ | $0.230e+2$ | $0.228e+2$
$SSE$ | $0.330e+0$ | $0.312e+0$
$SSE_{inf}$ | $0.722e-2$ | $0.674e-2$
$SSE_{sup}$ | $0.845e-1$ | $0.735e-1$
$GCV$ | $0.338e+1$ | $0.301e+1$

From Figure 6 and Table 5, which provide comparisons between the two
approximation methods, we can conclude that, as in previous cases, the double
penalization approximation performs better along the extreme regions, ensuring
that valuable information is not lost.

Figure 6: Covid19 Data smoothed by FS0 and FS2. Figure 7: Evaluation of the
number of clusters in the simulation dataset. The optimal number is indicated
from the dotted line.

The advantages we gain from using the double penalization approximation method
are observed in the application of clustering. To select the number of
clusters, we apply the Elbow method. From the Figure 7, we can determine that
the appropriate number of clusters is 4.

The clustering method we will use is kmeans.fd. In the first analysis we
perform is clustering on data approximated using FS2.

Figure 8 displays the time series of daily new cases for each country based on
the clusters. We can visually observe that the model varies significantly for
countries in different clusters.

We can observe that the pattern undergoes significant changes across countries
within distinct clusters. The countries in Cluster 1 exhibit a consistent
trend, followed by a sudden surge in cases from January 2022 to March 2022.
Conversely, countries in Cluster 3 demonstrate a decline in new cases during
the same period.

The countries in Cluster 2 show an increase in cases from the beginning of
January 2021, followed by a decrease starting in February. During the same
period, countries in Cluster 3 exhibited a consistent trend, followed by a
rise in cases in March 2021. It can be observed that the pattern of pandemic
development in Cluster 3 is often at the opposite pace to that of Cluster 4.

In general, each cluster of countries exhibits distinct characteristics that
set their pandemic patterns apart from those of the other clusters.

Table 6: Specific clustering results for daily new cases in 30 countries, in the case of FS2.  Cluster | Countries
---|---
Cluster 1 | Japan, Russian, Chile, Romania
| Ukraine, Germany, Czech, Netherlands
Cluster 2 | Indonesia, South Africa, Brazil, India, Colombia
Cluster 3 | Spain, Mexico, The U.S., The U.K.
| Argentina, Turkiye, France, Canada
| Italy, Poland, Peru
Cluster 4 | Philippines, Malaysia, Thailand
| The Islamic Republic of Iran, Bangladesh, Iraq
Table 7: Specific clustering results for daily new cases in 30 countries, in the case of FS0.  Cluster | Countries
---|---
Cluster 1 | Russia, Romania, Ukraine
| Netherlands, Germany, Czech
Cluster 2 | Japan, Brazil, Chile, Poland
Cluster 3 | Mexico, South Africa, The U.S., The U.K.
| Peru, India, Argentina, Colombia
| Canada, Turkiye, Italy, France, Spain
Cluster 4 | Malaysia, Thailand, The Islamic Republic Iran,
| Indonesia, Bangladesh, Iraq, Philippines

Table 6 shows the specific clustering results for the 30 countries. We can see
that all Cluster 1 and Cluster 4 countries are located in regions of Asia,
while most of Cluster 3 is located in different regions in the European
region.

For this reason we might think that geographical location can influence
clustering. However, geographical proximity is not a decisive factor for
clustering; just think of Cluster 2.

All the considerations we have made so far concern clustering on data
approximated using FS2. From the Figure 8 (upper figure) and the Table 6, we
can analyze the results of clustering on data approximated using FS0. The
Figure 8 (down figure) and Table 7 displays the different clusters for FS0:
the clusters are not well-defined, and the curves, not being correctly
allocated to the clusters, making it challenging to read and interpret the
data. The considerations that were previously made in this case are not as
clear: in Clusters 1 and 2, we observe a similar trend, especially towards the
end of 2021. Similar observations can be made for Clusters 3 and 4, which
exhibit similar trends. Therefore, in the case of $FS0$, specific
characteristics to characterize each cluster are not found.

Figure 8: K-means clustering results by using FS2 (upper figure) and FS0 (down
figure)

##  5 Conclusions

In this work, we considered the use of free knots spline in the context of
functional data estimation, analyzing, in particular, the impact of the
roughness regularization term. More specifically, we compared two different
penalty regularization schemes, namely a standard regularization scheme with a
one parameter roughness term handling the boundedness of the variability of
the function [8] and a two-penalty regularization scheme which controls
monotonicity and smoothness. Our numerical experiments seem to demonstrate
that compared to the free knots spline without any penalty, things do not
change much when using the one-parameter scheme, while the two-parameter
scheme shows notable improvements. However, the most significant advantages
appear in the data analysis phase based on the functional data approximation
obtained. In particular, when applied to simulated data, our method shows a
higher improved ability to detect ties, highlighting a clearer clustering
structure. A promising direction emerges from the non-linear effects linked to
the selection of knots in the initial basis. The incorporation of machine
learning method to enhance the analytical approach, open a different
perspective for continuous advancements in the domain of functional data
analysis.

Acknowledgements This work was partially supported by Italian Ministry of
University and Research (MIUR), PRIN 2022 Project: _Numerical Optimization
with Adaptive Accuracy and Applications to Machine Learning_ , grant n.
2022N3ZNAX, PNRR PRIN 2022 Project: _A multidisciplinary approach to evaluate
ecosystems resilience under climate change_ , grant n. P2022WC2ZZ, and by
Istituto Nazionale di Alta Matematica - Gruppo Nazionale per il Calcolo
Scientifico (INdAM-GNCS), by MIUR.

## References

  * Basna et al. [2022] R. Basna, H. Nassar, and K. Podgórski.  Data driven orthogonal basis selection for functional data analysis.  _Journal of Multivariate Analysis_ , 189:1048–1068, 2022\.  ISSN 0047-259X.
  * Ramsay and Silverman [1997] J.O. Ramsay and B. Silverman.  _Functional Data Analysis_.  Springer, New York, 1997.
  * Gervini [2006] D. Gervini.  Free-knot spline smoothing for functional data.  _Journal of the Royal Statistical Society_ , 68(4):671–687, 2006.
  * Denison et al. [2002] D.G.T. Denison, B.K. Mallick, and A.F.M. Smith.  Automatic bayesian curve fitting.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 60(2):333–350, 01 2002.  ISSN 1369-7412.  10.1111/1467-9868.00128.
  * Guidotti and Ardia [2020] E. Guidotti and D. Ardia.  Covid-19 data hub.  _The Open Journal_ , 5(51):2376, 2020.
  * Ferraty and Vieu [2006] F. Ferraty and P. Vieu.  _Nonparametric Functional Data Analysis: Theory and Practice_.  2006\.
  * Sivananthan [2016] S. Sivananthan.  Multi-penalty regularization in learning theory.  _Journal of Complexity_ , 36:141–165, 2016.
  * O’Sullivan [1986] F. O’Sullivan.  A statistical perspective on ill-posed inverse problems.  _Statistical Science_ , 1(4):502–518, 1986.
  * Rice and Silverman [1991] J.A. Rice and B.W. Silverman.  Estimating the mean and covariance structure nonparametrically when the data are curves.  _Journal of the Royal Statistical Society_ , 53(1):233–243, 1991.
  * Jupp [1978] D.L.B. Jupp.  Approximation to data by splines with free knots.  _SIAM Journal Numerical Analysis_ , 15:328–343, 1978.
  * Hartigan and Wong [1979] J.A. Hartigan and M.A. Wong.  A k-means clustering algorithm.  _Applied Statistics_ , 28:100–108, 1979.
  * Bouveyron et al. [2014] C. Bouveyron, E. Côme, and J. Jacques.  The discriminative functional mixture model for the analysis of bike sharing systems.  _HAL n.01024186_ , 2014.
  * Ramsay et al. [2015] J.O. Ramsay, J.S. Marron, L.M. Sangalli, and A. Srivastava.  Functional data analysis of amplitude and phase variation.  _Statistical Science_ , pages 468–484, 2015.
  * Tucker et al. [2013] J.D. Tucker, W. Wu, and A. Srivastava.  Generative models for functional data using phase and amplitude separation.  _Computational Statistics & Data Analysis_, 61:50–66, 2013.
  * Chen et al. [2014] H. Chen, P.T. Reiss, and T. Tarpey.  Optimally weighted l2 distance for functional data.  _Biometrics_ , 70:516–525, 2014.
  * Guidotti [2022] E. Guidotti.  A worldwide epidemiological database for covid-19 at fine-grained spatial resolution.  _Scientific Data_ , 9(1):1–7, 2022.
  * Luo et al. [2023] Z. Luo, L. Zhang, N. Liu, and Y. Wu.  Time series clustering of covid-19 pandemic-related data.  _Data Science and Management_ , 6:79–87, 2023.


,$ (276a) $\displaystyle(\mathcal{L}_{0}G_{0})_{\epsilon}$
$\displaystyle=8\pi^{2}R\,c_{0}\Bigg{[}\delta_{0}[f]\left(\sin
2\theta-2\theta\right)-\left(\cos
2\theta-1\right)\int_{\epsilon}^{\infty}\frac{dr}{r}\,\partial_{\theta}f\Bigg{]}_{\theta=\alpha}^{\theta=\beta}+O(\epsilon)\
,$ (276b) $\displaystyle(\mathcal{L}_{-1}G_{-1}[f])_{\epsilon}$
$\displaystyle=4\pi^{2}R\,c_{-1}\Bigg{[}e^{-2i\theta}\bigg{(}i\delta_{0}[f]\left(e^{-2i\theta}-2\right)-2\left(\cos
2\theta-1\right)\int_{\epsilon}^{\infty}\frac{dr}{r}\,\partial_{\theta}f\bigg{)}\Bigg{]}_{\theta=\alpha}^{\theta=\beta}+O(\epsilon)\
.$ (276c)

For generic values of $k$, we can take $\alpha=0$ and $\beta=\pi$ to obtain
(after taking the $\epsilon\to 0$ limit)

$\mathcal{L}_{k}G_{k}[f]=\frac{8\pi^{2}R\left(e^{2\pi
ik}-1\right)}{k}c_{k}\delta_{0}[f]\ ,$ (277)

so the Green’s function equation (235) is satisfied if we take

$c_{k}=\frac{k}{16\pi^{3}R^{2}\left(e^{2\pi ik}-1\right)}\ .$ (278)

However, like the 3d case this is not possible when $k\in\mathbb{Z}$, which is
the case of interest.

Let us follow the same steps as before to obtain the correct equations for the
tower of 5d two-point functions from the 6d Green’s function equation. The
time-ordered 6d two-point function is given in terms of the modes by

$\bra{\Omega}T\\{\hat{\phi}(\hat{x}^{+},\hat{x}^{-},\hat{x})\bar{\hat{\phi}}(0)\\}\ket{\Omega}=\cos^{2}\left(\frac{x^{+}}{2R}\right)\sum_{k\in\mathbb{Z}}e^{-\frac{ikx^{+}}{R}}G_{k}(x^{-},x)\
.$ (279)

The modes of $\hat{\phi}$ are related to those of $\phi$ by

$\hat{\phi}_{k}=\frac{1}{4}\left(\phi_{k-1}+2\phi_{k}+\phi_{k+1}\right)\ ,$
(280)

so the conditions (249) become

$\phi_{k}\ket{\Omega}=\bar{\phi}_{-k}\ket{\Omega}=0\ \forall\,k>-1\ .$ (281)

This means that the mode two-point functions are of the form

$G_{k}(x^{-}.x)=\begin{cases}\Theta(x^{-})\bra{\Omega}\phi_{k}(x^{-},x)\bar{\phi}_{k}(0)\ket{\Omega}&k>0\\\
0&k=0\\\
\Theta(-x^{-})\bra{\Omega}\bar{\phi}_{k}(0)\phi_{k}(x^{-},x)\ket{\Omega}&k<0\end{cases}\
.$ (282)

In particular, we immediately have the condition $c_{0}$ = 0.

Acting with the 6d d’Alembertian

$\hat{\square}=\cos^{2}\left(\frac{x^{+}}{2R}\right)\left(-2\left(\partial_{+}+\frac{1}{R}\tan\left(\frac{x^{+}}{2R}\right)\right)+D^{2}\right)$
(283)

and rearranging, we see that the Green’s function equation

$\hat{\square}\bra{\Omega}T\\{\hat{\phi}(\hat{x}^{+},\hat{x}^{-},\hat{x})\bar{\hat{\phi}}(0)\\}\ket{\Omega}=i\delta^{(6)}(\hat{x})$
(284)

reduces to the tower of equations

$\mathcal{L}_{k-1}G_{k-1}+\mathcal{L}_{k+1}G_{k+1}+2\mathcal{L}_{k}G_{k}=\frac{2i}{\pi
R}\delta(x^{-})\delta^{(4)}(x)$ (285)

for $k\in\mathbb{Z}$. Substituting in $\alpha=0$, $\beta=\frac{\pi}{2}$ for
$k>0$ and $\alpha=\pi$, $\beta=\frac{\pi}{2}$ for $k<0$, the distributional
derivatives of the two-point functions are (taking $\epsilon\to 0$)

$\mathcal{L}_{k}G_{k}[f]=\frac{8\pi^{2}iRc_{k}}{k}\left((-1)^{k}-1\right)\delta_{0}[f]+16\pi^{2}R(-1)^{k}c_{k}\int_{0}^{\infty}\frac{dr}{r}\left(\partial_{\theta}f\right)_{\theta=\frac{\pi}{2}}$
(286a) for $k>0$, and
$\mathcal{L}_{k}G_{k}[f]=-\frac{8\pi^{2}iRc_{k}}{k}\left((-1)^{k}-1\right)\delta_{0}[f]-16\pi^{2}R(-1)^{k}c_{k}\int_{0}^{\infty}\frac{dr}{r}\left(\partial_{\theta}f\right)_{\theta=\frac{\pi}{2}}$
(286b)

for $k<0$. We see that the equations (285) are satisfied if we take

$c_{k}=\begin{cases}-\frac{k}{16\pi^{3}R^{2}}&k>0\\\ 0&k=0\\\
\frac{k}{16\pi^{3}R^{2}}&k<0\end{cases}\ .$ (287)

A similar calculation to the resummed 4d two-point function then
gives131313Suppressing the terms arising from the $i\epsilon$ prescription.

$\bra{\Omega}T\\{\hat{\phi}(\hat{x}^{+},\hat{x}^{-},\hat{x})\bar{\hat{\phi}}(0)\\}\ket{\Omega}=\frac{1}{4\pi^{3}\left(-2\hat{x}^{+}\hat{x}^{-}+\hat{x}^{i}\hat{x}^{i}\right)^{2}}\
,$ (288)

which we recognise as the two-point function of a canonically normalised 6d
scalar field.

###  C.3 The DLCQ Limit of $SU(1,n)$ Scalar Fields

Let us see how we recover the familiar DLCQ picture from this. As discussed in
section 4.1, we can think of the DLCQ limit by writing $k$ and $R$ as

$\displaystyle k$ $\displaystyle=Nk_{+}\ ,$ (289a) $\displaystyle R$
$\displaystyle=NR_{+}\ ,$ (289b)

and taking the limit $N\to\infty$ keeping $k_{+}$ and $R_{+}$ fixed. To get a
finite action after this limit we must also rescale the mode fields by
$\frac{1}{\sqrt{N}}$. Since we have

$\frac{k+c}{R}=\frac{k_{+}}{R_{+}}+O\Big{(}\frac{1}{N}\Big{)}$ (290)

for any $c$ that doesn’t scale with $N$, sums over $k\in\mathbb{Z}+\alpha$
become sums over $k_{+}\in\mathbb{Z}$ in the $N\to\infty$ limit. From here
onwards we will drop pluses for ease of notation. The action for the DLCQ
reduction of a complex scalar field in $2n$-dimensions can then be seen to be

$S=2\pi R\sum_{k\in\mathbb{Z}}\int
dx^{-}d^{2n-2}x\left(\frac{2ik}{R}\bar{\phi}_{k}\partial_{-}\phi_{k}+\bar{\phi}_{k}\partial^{2}\phi_{k}\right)\
.$ (291)

We observe that in the DLCQ limit the Green’s function-like equations (256)
and (285) both reduce to the standard Schrödinger field Green’s function
equation

$\left(\frac{2ik}{R}\partial_{-}+\partial^{2}\right)\bra{\Omega}T\\{\phi_{k}(x^{-},x)\bar{\phi}_{k}(0)\\}\ket{\Omega}=\frac{i}{2\pi
R}\delta(x^{-})\delta^{2n-2}(x)\ ,$ (292)

which is the obvious set of equations obtained from the action (291). The
solution to this is well known, and we briefly review it. Using the fact that
the two-point function satisfies the classical equation of motion away from
the origin and the DLCQ vacuum condition

$\phi_{k}\ket{\Omega}=\bar{\phi}_{-k}\ket{\Omega}=0$ (293)

for $k>0$, we find

$G_{k}\equiv\bra{\Omega}T\\{\phi_{k}(x^{-},x)\bar{\phi}_{k}(0)\\}\ket{\Omega}=\begin{cases}c_{k}\Theta(x^{-})(x^{-})^{1-n}\exp\left(\frac{ikx^{2}}{2Rx^{-}}\right)&k>0\\\
0&k=0\\\
c_{k}\Theta(-x^{-})(x^{-})^{1-n}\exp\left(\frac{ikx^{2}}{2Rx^{-}}\right)&k<0\end{cases}\
.$ (294)

Let us focus on the $k>0$ sector. Since $G_{k}$ is only non-differentiable at
$x^{-}=0$, we define the distributional derivative as

$(\mathcal{L}_{k}G_{k})_{\epsilon}=\left(\int_{-\infty}^{-\epsilon}dx^{-}+\int_{\epsilon}^{\infty}dx^{-}\right)\int_{\mathbb{R}^{2n-2}}d^{2n-2}x\left(-\frac{2ik}{R}\partial_{-}f+\partial^{2}f\right)G_{k}$
(295)

where we use $\mathcal{L}_{k}$ to denote the differential operator in (292).
Substituting in our expression for $G_{k}$ and formally treating our spatial
integrals as if they are convergent we find

$\displaystyle(\mathcal{L}_{k}G_{k})_{\epsilon}$
$\displaystyle=\frac{2ik}{R}\epsilon^{1-n}c_{k}\,\delta_{0}[f]\left(\int_{\mathbb{R}}dx\,e^{\frac{ikx^{2}}{2R\epsilon}}\right)^{2n-2}+O(\epsilon)$
$\displaystyle=\frac{2ik}{R}\left(-\frac{2\pi
R}{ik}\right)^{n-1}c_{k}\,\delta_{0}[f]+O(\epsilon)\ .$ (296)

We can now take $\epsilon\to 0$ and compare with (292) to find the coefficient

$c_{k}=\frac{1}{4\pi k}\left(-\frac{ik}{2\pi R}\right)^{n-1}$ (297)

for $k>0$. Similarly, working through the same calculation for the $k<0$
sector gives

$c_{k}=-\frac{1}{4\pi k}\left(-\frac{ik}{2\pi R}\right)^{n-1}\ .$ (298)

We now consider directly taking the DLCQ limit of the $SU(1,n)$ two-point
functions. We saw above that for a complex scalar field these are of the form

$G_{k}=\frac{c_{k}F_{k}(\theta)e^{2ik\theta}}{r^{n-1}}\ ,$ (299)

where

$\displaystyle r$
$\displaystyle=\sqrt{(x^{-})^{2}+\left(\frac{x^{2}}{4R}\right)^{2}}\ ,$ (300a)
$\displaystyle\theta$
$\displaystyle=-\frac{i}{2}\ln\left(\frac{1+\frac{ix^{2}}{4Rx^{-}}}{1-\frac{ix^{2}}{4Rx^{-}}}\right)\
.$ (300b)

In the DLCQ limit we have

$\displaystyle r$ $\displaystyle=x^{-}+O\left(\frac{1}{R}\right)\ ,$ (301a)
$\displaystyle\theta$
$\displaystyle=\frac{x^{2}}{4Rx^{-}}+O\left(\frac{1}{R^{2}}\right)\ ,$ (301b)

so the limit of the two-point function is (assuming that no limit needs to be
taken in the coefficients)

$G_{k}\to\begin{cases}c_{k}\Theta(x^{-})(x^{-})^{1-n}\exp\left(\frac{ikx^{2}}{2Rx^{-}}\right)&k>0\\\
0&k=0\\\
c_{k}\Theta(-x^{-})(x^{-})^{1-n}\exp\left(\frac{ikx^{2}}{2Rx^{-}}\right)&k<0\end{cases}\
,$ (302)

where we add the $k=0$ case in by hand for $n=2$. As hoped, this is identical
to the DLCQ correlation functions in (294) as long as the coefficients agree
in the $SU(1,n)$ and DLCQ calculations. Comparing (297) and (298) with (262)
and (266) for $n=2$ and (287) for $n=3$, we see that this is indeed true.

It will be interesting to resum the modes and obtain the higher-dimensional
DLCQ two-point function. Doing this gives

$\displaystyle\bra{\Omega}T\\{\phi(x^{+},x^{-},x)\bar{\phi}(0)\\}\ket{\Omega}$
$\displaystyle=\sum_{k=1}^{\infty}\left(e^{-\frac{ikx^{+}}{R}}G_{k}+e^{\frac{ikx^{+}}{R}}G_{-k}\right)$
$\displaystyle=\frac{1}{4\pi}\left(\frac{-i}{2\pi
Rx^{-}}\right)^{n-1}\left(\Theta(x^{-})\mathrm{Li}_{2-n}(q)+(-1)^{n-1}\Theta(-x^{-})\mathrm{Li}_{2-n}(q^{-1})\right)\
,$ (303)

where we have defined $q$ as in (120) and used the usual $i\epsilon$
prescription (86) to ensure convergence of the sum. We can then use the
identity141414Note that this means we are restricting our attention to the
DLCQ of scalar field theories in even spacetime dimensions. (121) to simplify
the two-point function to

$\bra{\Omega}T\\{\phi(x^{+},x^{-},x)\bar{\phi}(0)\\}\ket{\Omega}=\frac{1}{4\pi}\left(\frac{-i}{2\pi
Rx^{-}}\right)^{n-1}\mathrm{Li}_{2-n}(q)\ .$ (304)

For the exceptional $n=2$ case we instead have

$\bra{\Omega}T\\{\phi(x^{+},x^{-},x)\bar{\phi}(0)\\}\ket{\Omega}=\frac{1}{16\pi^{2}R}\left(\frac{1}{x^{-}}\cot\left(\frac{x_{\mu}x^{\mu}}{4R\,x^{-}}\right)+\frac{i}{\absolutevalue{x^{-}}}\right)\
,$ (305)

where as before

$x_{\mu}x^{\mu}\equiv-2x^{+}x^{-}+x^{i}x^{i}\ .$ (306)

In order to recover the Minkowski two-point function we need to take the
$R\to\infty$ limit to decompactify the null direction; using the large-$R$
expansions (128) we see that in both cases we get

$\lim_{R\to\infty}\Big{(}\bra{\Omega}T\\{\phi(x^{+},x^{-},x)\bar{\phi}(0)\\}\ket{\Omega}\Big{)}=\frac{(n-2)!}{4\pi^{n}\left(x_{\mu}x^{\mu}\right)^{n-1}}\
,$ (307)

which we recognise as the canonically-normalised two-point function of a
massless scalar field in $2n$-dimensional Minkowski spacetime.

## References

  * (1) J. Gomis and H. Ooguri, Nonrelativistic closed string theory, J. Math. Phys. 42 (2001) 3127–3151, [hep-th/0009181].
  * (2) J. Gomis, J. Gomis, and K. Kamimura, Non-relativistic superstrings: A New soluble sector of AdS(5) x S**5, JHEP 12 (2005) 024, [hep-th/0507036].
  * (3) E. Bergshoeff, K. T. Grosvenor, J. Lahnsteiner, Z. Yan, and U. Zorba, Non-Lorentzian IIB Supergravity from a Polynomial Realization of SL(2,R), arXiv:2306.04741.
  * (4) G. Oling and Z. Yan, Aspects of Nonrelativistic Strings, Front. in Phys. 10 (2022) 832271, [arXiv:2202.12698].
  * (5) T. Harmark, J. Hartong, L. Menculini, N. A. Obers, and G. Oling, Relating non-relativistic string theories, JHEP 11 (2019) 071, [arXiv:1907.01663].
  * (6) T. Harmark and M. Orselli, Spin Matrix Theory: A quantum mechanical model of the AdS/CFT correspondence, JHEP 11 (2014) 134, [arXiv:1409.4417].
  * (7) S. Baiguera, T. Harmark, and N. Wintergerst, Nonrelativistic near-BPS corners of $\mathcal{N}=4$ super-Yang-Mills with $SU(1,1)$ symmetry, JHEP 02 (2021) 188, [arXiv:2009.03799].
  * (8) S. Baiguera, T. Harmark, and Y. Lei, The Panorama of Spin Matrix theory, JHEP 04 (2023) 075, [arXiv:2211.16519].
  * (9) L. Donnay, A. Fiorucci, Y. Herfray, and R. Ruzziconi, Carrollian Perspective on Celestial Holography, Phys. Rev. Lett. 129 (2022), no. 7 071602, [arXiv:2202.04702].
  * (10) L. Donnay, A. Fiorucci, Y. Herfray, and R. Ruzziconi, Bridging Carrollian and celestial holography, Phys. Rev. D 107 (2023), no. 12 126027, [arXiv:2212.12553].
  * (11) A. Bagchi, S. Banerjee, R. Basu, and S. Dutta, Scattering Amplitudes: Celestial and Carrollian, Phys. Rev. Lett. 128 (2022), no. 24 241601, [arXiv:2202.08438].
  * (12) C. Duval, G. W. Gibbons, and P. A. Horvathy, Conformal Carroll groups, J. Phys. A 47 (2014), no. 33 335204, [arXiv:1403.4213].
  * (13) C. Duval, G. W. Gibbons, and P. A. Horvathy, Conformal Carroll groups and BMS symmetry, Class. Quant. Grav. 31 (2014) 092001, [arXiv:1402.5894].
  * (14) D. T. Son, Toward an AdS/cold atoms correspondence: A Geometric realization of the Schrodinger symmetry, Phys. Rev. D 78 (2008) 046003, [arXiv:0804.3972].
  * (15) K. Yamawaki, Zero mode problem on the light front, in 10th Summer School and Symposium on Nuclear Physics: QCD, Light cone Physics and Hadron Phenomenology (NuSS 97), pp. 116–199, 1, 1998.  hep-th/9802037.
  * (16) N. Lambert, A. Lipstein, R. Mouland, and P. Richmond, Five-dimensional non-Lorentzian conformal field theories and their relation to six-dimensions, JHEP 03 (2021) 053, [arXiv:2012.00626].
  * (17) N. A. Nekrasov, Seiberg-Witten prepotential from instanton counting, Adv. Theor. Math. Phys. 7 (2003), no. 5 831–864, [hep-th/0206161].
  * (18) N. Lambert, R. Mouland, and T. Orchard, Non-Lorentzian SU(1, n) Spacetime Symmetry In Various Dimensions, Front. in Phys. 10 (2022) 864800, [arXiv:2112.14860].
  * (19) J. J. Heckman and T. Rudelius, Top Down Approach to 6D SCFTs, J. Phys. A 52 (2019), no. 9 093001, [arXiv:1805.06467].
  * (20) N. Lambert, A. Lipstein, and P. Richmond, Non-Lorentzian M5-brane Theories from Holography, JHEP 08 (2019) 060, [arXiv:1904.07547].
  * (21) N. Lambert and T. Orchard, Non-Lorentzian Avatars of (1,0) Theories, arXiv:2011.06968.
  * (22) I. L. Shapiro, Covariant Derivative of Fermions and All That, Universe 8 (2022), no. 11 586, [arXiv:1611.02263].
  * (23) J. de Boer, J. Hartong, N. A. Obers, W. Sybesma, and S. Vandoren, Carroll stories, arXiv:2307.06827.
  * (24) S. Baiguera, G. Oling, W. Sybesma, and B. T. Søgaard, Conformal Carroll scalars with boosts, SciPost Phys. 14 (2023), no. 4 086, [arXiv:2207.03468].
  * (25) V. Balasubramanian, J. de Boer, M. M. Sheikh-Jabbari, and J. Simon, What is a chiral 2d CFT? And what does it have to do with extremal black holes?, JHEP 02 (2010) 017, [arXiv:0906.3272].
  * (26) G. ’t Hooft, Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle, Phys. Rev. D 14 (1976) 3432–3450. [Erratum: Phys.Rev.D 18, 2199 (1978)].
  * (27) L. Bhardwaj and Y. Tachikawa, Classification of 4d N=2 gauge theories, JHEP 12 (2013) 100, [arXiv:1309.5160].
  * (28) D. Gaiotto, N=2 dualities, JHEP 08 (2012) 034, [arXiv:0904.2715].
  * (29) D. Gaiotto, G. W. Moore, and A. Neitzke, Wall-crossing, Hitchin systems, and the WKB approximation, Adv. Math. 234 (2013) 239–403, [arXiv:0907.3987].
  * (30) D. Z. Freedman and A. Van Proeyen, Supergravity.  Cambridge Univ. Press, Cambridge, UK, 5, 2012.
  * (31) F. G. Friedlander and M. S. Joshi, Introduction to the theory of distributions.  Cambridge Univ. Press, 2003.