\section{Introduction}

In the recent study of string theory, one of the most important issues is to understand
a holographic principle \cite{Hooft, Susskind} at the full quantum level
(For a review see \cite{Bousso}).
The AdS/CFT correspondence \cite{M, GKP, Witten} is a realization of the holography.
This is, however, a conjectured relation and there is no rigorous proof so far.
The integrable structure behind the correspondence at the planar level
has played an important role in checking conjectured relations in non-BPS regions
(For a comprehensive review, see \cite{review}).
But the proof is still far from the completion and furthermore
it does not seem likely that the integrability would work in the presence of a black hole.

\medskip

Towards the complete understanding of holography,
it is significant to try to construct a simple toy model of quantum holography.
In fact, Kitaev proposed such a model \cite{Kitaev}, which is a variant of the Sachdev-Ye (SY) model
\cite{SY}.  More concretely, this model is a one-dimensional system composed of $N \gg 1$
fermions with a random, all-to-all quartic interaction.
This model is now called the Sachdev-Ye-Kitaev (SYK) model\footnote{For the recent progress
on the SYK model, see \cite{PR,MS,GR,FGMS,Witten2,many,Klebanov,Nishinaka}. }.
It should be remarked that the Lyapunov exponent computed from an out-of-time-order
four-point function \cite{LO, Kitaev-corr} in the SYK model saturates the bound presented
in \cite{MSS}. This is the onset to open a window to a toy model of holography
because the Lyapunov exponent of black hole in Einstein gravity is $2\pi/\beta$ \cite{SS1,SS2},
where $\beta$ is the inverse of the Hawking temperature
(For a black hole S-matrix approach, see \cite{P}.).

\medskip

A promising candidate of the gravity dual for the SYK model is a 1+1 D dilaton gravity.
This system was originally introduced by Jackiw \cite{Jackiw} and Teitelboim \cite{Teitelboim}
(For a nice reviw of the 1+1 D dilaton gravity system, see \cite{V}).
From a renewed interest, the dilaton gravity with a certain dilaton potential
was intensively studied in the recent work \cite{AP},
and this model is called the Almheiri-Polchinski (AP) model. A black hole solution exists
as a vacuum solution of the AP model. They studied its various properties
like the RG flow structure at zero temperature, the Bekenstein-Hawking entropy, the renormalized
boundary stress tensor, and the contribution of conformal matter fields to the entropy.
For the recent progress on the AP model, see \cite{MSY,EMV}.

\medskip

In this paper, we are concerned with deformations of the AP model.
Why is it so interesting to study the deformations?
There are some observation and motivation based on the recent progress.
The first is an observation that the SY model is constructed
by performing a disordered quench to an {\it isotropic} quantum Heisenberg magnet \cite{SY}.
The Heisenberg model itself is integrable. Hence, supposing that the conjectured duality is true,
it is natural to expect that integrable deformations of it lead to the associated deformations
of the AP model.
The second is a motive to understand the holographic duals of deformed AdS$_2$ geometries.
Recently, a systematic way to perform integrable deformations, which is called the Yang-Baxter
deformation \cite{Klimcik,DMV,MY-YBE}\footnote{For the (affine) symmetry algebras,
see \cite{symmetry} and \cite{Sch-alg}).},
has been intensively
studied\footnote{For Yang-Baxter deformations of type IIB superstring on
AdS$_5\times$S$^5$\,, see \cite{DMV2,ABF,ABF2,HRT,LRT,HT-sol,scale,
KMY-Jordanian-typeIIB,LM-MY,MR-MY,Sch-MY,KMY-SUGRA,Stijn1,Stijn2,MY-duality,KKSY,KY-super,
HvT,ORSY,BW,OvT,HT-NA, BW-NA,Stijn3}}.
However, the holographic duals of the deformed geometries
have been poorly understood. In particular, even the location of the holographic screen has not been
clarified, though there is a proposal \cite{Kameyama}.
Hence it is important to get much deeper understanding of the simplest
case like AdS$_2$\,. Furthermore, the Yang-Baxter deformation is not
applicable to black hole geometries in general, because those cannot be described usually
as a coset, homogeneous space. However, it is not the case for a 1+1 D black hole
presented in this paper.

\medskip

Based on the observation and motivation described above,
we will study deformations of the AP model by employing the Yang-Baxter deformation technique.
The general deformed AdS$_2$ metric
becomes a solution of a deformed AP model. In particular, the dilaton potential is deformed
from a simple quadratic form to a hyperbolic function-type potential
similarly to integrable deformations. A specific solution is a deformed black hole solution.
Because the deformation makes the spacetime structure
around the boundary change drastically and a new naked singularity appears,
the holographic interpretation is far from trivial. The Hawking temperature is the same
as the undeformed case but the Bekenstein-Hawking entropy is modified due to the deformation.
This entropy can also be reproduced by evaluating the renormalized stress tensor
with an appropriate counter term on the regularized screen close to the singularity.

\medskip

This paper is organized as follows.
In section 2 we study the most general Yang-Baxter deformation of AdS$_2$\,.
Section 3 introduces the classical action of 1+1 D dilaton gravity
and the AP model as a special case.
In section 4 we study deformations of the AP model.
The most general deformed metric constructed in section 2
becomes a solution with a deformed dilaton potential. This deformed
system allows a black hole solution as a specific solution like in the AP model.
The Bekenstein-Hawking entropy is also computed.
In section 5, we revisit the black hole entropy from the viewpoint of the renormalized
boundary stress tensor. Putting a regularized screen close to a singularity, we
evaluate the renormalized boundary stress tensor with an appropriate counter term.
The resulting entropy nicely agrees with the Bekenstein-Hawking entropy
computed in section 4.
Section 6 is devoted to conclusion and discussion.

\section{Yang-Baxter deformations of AdS$_2$}

In this section, we consider the most general Yang-Baxter deformation of the AdS$_2$ metric.
First of all, we briefly describe a coset construction of the Poincar\'e AdS$_2$
Then we study the most general Yang-Baxter deformation of Poincar\'e AdS$_2$\,.
As a result, we obtain a three-parameter family of deformed AdS$_2$ spaces.

\subsection{Coset construction of AdS$_2$}

Let us recall a coset construction of the Poincar\'e AdS$_2$ metric (For the detail of the coset construction, for example, see \cite{BKLSY}).

\medskip

The starting point is that the AdS$_2$ geometry is represented by a coset
\begin{equation}
\mbox{AdS}_2 ~~= ~~ SL(2)/U(1)\,. \label{coset}
\end{equation}
By using the coordinates $t$ and $z$\,, a coset representative $g$ is parametrized as
\begin{equation}
g = {\rm exp}\left[ t H \right]
{\rm exp}\left[ (\log z) D \right]\,,
\label{group para}
\end{equation}
where $H$ and $D$ are the time translation and dilatation generators, respectively.
By involving the special conformal generator $C$\,,
the $\mathfrak{sl}(2)$ algebra in the conformal basis is spanned as
\begin{eqnarray}
[D,H] = H\,, \qquad [C,H] = 2 D\,, \qquad [D,C] = -C\,.
\end{eqnarray}
These generators can be represented by the $\mathfrak{so}(1,2)$
ones $T_I~~(I=0,1,2)$ like
\begin{eqnarray}
H \equiv T_0 + T_2\,, \qquad C \equiv T_0 - T_2\,, \qquad
D \equiv T_1\,,
\end{eqnarray}
where $T_I$'s satisfy the commutation relations:
\begin{equation}
[T_0,T_1] = - T_2\,, \qquad  [T_1,T_2] = T_0\,, \qquad [T_2,T_0] = -T_1\,.
\label{symmetric}
\end{equation}
In the following, we will work with $T_I$'s in the fundamental representation,
\begin{eqnarray}
T_0=\frac{i}{2}\sigma_1\,,\qquad T_1=\frac{1}{2}\sigma_2\,,\qquad T_2=\frac{1}{2}\sigma_3\,,
\end{eqnarray}
where $\sigma_i~(i=1,2,3)$ are the standard Pauli matrices.

\medskip

Note here that the coset (\ref{coset}) is symmetric as one can readily understand
from (\ref{symmetric})\,. When vector spaces $\mathfrak{h}$ and $\mathfrak{m}$ are spanned
as
\[
\mathfrak{h} = \mbox{span}_{\mathbb{R}}\langle\, T_2\, \rangle\,,
\qquad \mathfrak{m} = \mbox{span}_{\mathbb{R}}\langle T_0,T_1 \rangle\,,
\]
the $\mathbb{Z}_2$-grading structure is expressed as
\begin{eqnarray}
[\mathfrak{h},\mathfrak{h}] \subset \mathfrak{h}\,, \qquad
[\mathfrak{m},\mathfrak{h}] \subset \mathfrak{m}\,, \qquad
[\mathfrak{m},\mathfrak{m}] \subset \mathfrak{h}\,.
\end{eqnarray}
When representing the $\mathfrak{sl}(2)$ algebra by
a direct product (as vector spaces)\,,
\[
\mathfrak{sl}(2) = \mathfrak{h} \oplus \mathfrak{m}\,,
\]
the projection operator $P :~\mathfrak{sl}(2) \,\to\, \mathfrak{m}$ can be defined as
\begin{eqnarray}
P(X) \equiv \frac{ {\rm Tr}(X~T_0)}{{\rm Tr} (T_0~ T_0) } T_0
+  \frac{ {\rm Tr}(X~T_1)}{{\rm Tr} (T_1~T_1) }T_1\,, \qquad X \in \mathfrak{sl}(2)\,.
\label{projection}
\end{eqnarray}

\medskip

Now the Poincar\'e AdS$_2$ metric can be computed by performing coset construction.
The left invariant one-form $J = g^{-1}{\rm d} g$ is expanded as
\[
J = e^{0}\, T_0 + e^1\, T_1 + \frac{1}{2}\omega^{01}\,T_2\,.
\]
Here $e^0$ and $e^1$ are the zweibeins, and $\omega^{01}$ is the spin connection.
With the parametrization (\ref{group para})\,, the zweibeins are given by
\[
e^0 = \frac{{\rm d} t} {z}\,, \qquad e^1 = \frac{{\rm d} z}{z}\,.
\]
By using the projection operator $P$ in (\ref{projection}) and the explicit expressions of
the zweibeins $e^0$ and $e^1$\,,
the resulting metric is obtained as
\begin{eqnarray}
{\rm d} s^2 &=& 2 \mbox{Tr} \left[ J P(J) \right]  = - e^0 e^0 + e^1 e^1\nonumber \\
&=& \frac{-{\rm d} t^2 +{\rm d} z^2}{z^2}\,.
\label{AdS2}
\end{eqnarray}
This is nothing but the AdS$_2$ metric in the Poincar\'e coordinates.

\medskip

Hereafter, it is often convenient to use the the light-cone coordinates defined as
\begin{equation}
x^{\pm} \equiv t \pm z\,.
\label{lc}
\end{equation}
Then the metric is rewritten as
\begin{eqnarray}
{\rm d} s^2 = -{\rm e}^{2\omega (x^+,x^-)}\,{\rm d} x^+ {\rm d} x^-
= -\frac{4 {\rm d} x^+ {\rm d} x^-}{(x^+-x^-)^2}\,.
\end{eqnarray}
The exponential factor will play an important role in later discussion.

\subsection{The general Yang-Baxter deformation}

Let us next consider Yang-Baxter deformations of the AdS$_2$ metric (\ref{AdS2})\,.
In the usual discussion, Yang-Baxter deformations \cite{Klimcik, DMV,MY-YBE}
are performed for 2D non-linear sigma models.
Then the anti-symmetric two-form is also involved as well as the metric.
Here we will concentrate on the metric part only.

\medskip

The prescription of the deformation is very simple. It is just to insert a factor as follows:
\begin{eqnarray}
{\rm d} s^2 = 2 \mbox{Tr} \left[J\, \frac{1}{1-2\eta R_g \circ P}\, P(J)\right]\,.
\label{YB-metric}
\end{eqnarray}
Here $\eta$ is a constant parameter which measures the deformation.
Then $R_g$ is defined as a chain of operation like
\begin{eqnarray}
R_g(X) \equiv g^{-1} \circ R(g X g^{-1}) \circ g\,,
\end{eqnarray}
where $g$ is the group element in (\ref{group para})\,.
The key ingredient is a linear operator $R:~\mathfrak{sl}(2) \to \mathfrak{sl}(2)$\,,
and satisfy the (modified) classical Yang-Baxter equation [(m)CYBE]:
\begin{eqnarray}
[R(X),R(Y)]-R([R(X),Y]+[X,R(Y)])=c \cdot [X,Y]
\qquad (X\,,Y\in \mathfrak{sl}(2))\,. \label{mCYBE}
\end{eqnarray}
Here $c$ is a real constant parameter.
The case with $c \neq 0$ is the mCYBE and
the case with $c=0$ is the homogeneous CYBE.

\medskip

We consider the most general deformations with the following $R$-operator
\begin{eqnarray}
R(T_I) &=& \widetilde{\Omega}_{IJ}M^{JK}\,T_K\,, \label{general}
\end{eqnarray}
where $\widetilde{\Omega}_{IJ}$ and $M^{IJ}$ are defined as\footnote{A similar study
was done for the Nappi-Witten model \cite{NW} in \cite{KY-NW}.
Note that the present definition of $M^{IJ}$ is slightly different
from the one in \cite{KY-NW}.}
\begin{eqnarray}
\widetilde{\Omega}_{IJ} &\equiv& \text{Tr}(T_IT_J)=\frac{1}{2}\eta_{IJ}\,,  \qquad
M^{IJ} \equiv
\begin{pmatrix}
0&~m_1~&~m_2\;\\
-m_1&0&m_3\\
-m_2&-m_3&0\\
\end{pmatrix}\,,
\label{ansatz}
\end{eqnarray}
Putting the ansatz (\ref{general}) into the (m)CYBE (\ref{mCYBE}) leads to an algebraic relation,
\begin{eqnarray}
-m_1^2-m_2^2+m_3^2=4c\,.
\label{constraint}
\end{eqnarray}

\medskip

After evaluating the expression (\ref{YB-metric}) with the general ansatz (\ref{general})\,,
one can obtain the following metric:
\begin{eqnarray}
{\rm d} s^2 = \frac{-{\rm d} t^2 + {\rm d} z^2}{z^2-\eta^2\left(\alpha + \beta t+\gamma (-t^2+z^2) \right)^2}\,.
\label{3-metric}
\end{eqnarray}
Here $\alpha$\,, $\beta$ and $\gamma$ are defined
as linear combinations of $m_p~~(p=1,2,3)$ as follows:
 \begin{eqnarray}
\alpha \equiv \frac{1}{2}(m_1+m_3)\,,\qquad
\beta \equiv -m_2\,,\qquad
\gamma \equiv \frac{1}{2}(m_1-m_3)\,.
\end{eqnarray}
When $\eta=0$\,, the undeformed metric (\ref{AdS2}) is reproduced.
Note here that the four constant parameters $m_p~(p=1,2,3)$ and $c$ appear
in our discussion. Then a constraint (\ref{constraint})\,, which comes from the (m)CYBE,
is imposed. Hence, three of them are independent each other.

\medskip

The Ricci scalar curvature of the metric (\ref{3-metric}) is
\begin{eqnarray}
\label{Ricci}
R=-2\left(1-\widetilde{\omega}\eta^2\right)\,
\frac{z^2+\eta^2\left(\alpha+\beta t +\gamma(-t^2+z^2)\right)^2}{
z^2-\eta^2\left(\alpha+\beta t +\gamma(-t^2+z^2)\right)^2}\,,
\end{eqnarray}
where we have introduced a new quantity,
\begin{eqnarray}
\widetilde{\omega} &\equiv&  \beta^2+4\alpha\gamma=m_1^2+m_2^2-m_3^2 \nonumber \\
&=& -4c\,.
\end{eqnarray}
At the last equality, the (m)CYBE (\ref{mCYBE}) has been utilized.
The scalar curvature (\ref{Ricci}) changes (even its sign) depending on the values of parameters
and coordinates, while it becomes a constant $-2$ in the undeformed limit $\eta \to 0$\,.
The expression (\ref{Ricci}) indicates that the deformed geometry contains both AdS and dS
in general.

\section{A brief review of the AP model}

In this section, we shall introduce the classical action of 1+1 D dilaton gravity
system. Then we briefly describe the AP model and its properties related to our later discussion.

\subsection{1+1 D dilaton gravity system}

The dilaton gravity system in 1+1 dimensions is composed of the metric $g_{ab}~(a,b=0,1)$
and the dilaton $\Phi$\,. The coordinates are parametrized as $x^a = (x^0,x^1) = (t,z)$\,.

\medskip

The classical action $S$ is given by
\begin{eqnarray}
\label{action}
&&S=S_{g,\Phi}+S_{\rm matter}\,,\no\\
&&S_{g,\Phi}=\frac{1}{16\pi G}\int\!d^2 x\, \sqrt{-g}\left(\Phi^2 R - U(\Phi)\right)\,,\no\\
&&S_{\rm matter}=\frac{1}{32\pi G}\int\!d^2 x\, \sqrt{-g}\, \Omega(\Phi)\, (\nabla f)^2\,.
\end{eqnarray}
Here $G$ is the Newton constant in 1+1 dimensions and $U(\Phi)$ is the dilaton potential.

\medskip

In the following, we will work with the metric in the conformal gauge,
\begin{eqnarray}
{\rm d} s^2 &=& -{\rm e}^{2 \omega(x^+,x^-)}{\rm d} x^+ {\rm d} x^-\,,
\end{eqnarray}
where the light-cone coordinates are defined in (\ref{lc})\,.

\medskip

Then the equations of motion are given by
\begin{eqnarray}
\partial_+(\Omega\partial_- f)+\partial_-(\Omega\partial_+ f)&=&0\,,\no\\
4\partial_+\partial_-\Phi^2 - {\rm e}^{2\omega}U(\Phi)&=&0\,,\no\\
2\partial_+({\rm e}^{-2\omega}\partial_- {\rm e}^{2\omega})
-\frac{1}{2}e^{2\omega}\partial_{\Phi^2}U(\Phi)&=&(\partial_{\Phi^2}\Omega)\partial_
+f \partial_-f\,,\no\\
- {\rm e}^{2\omega}\partial_+({\rm e}^{-2\omega}\partial_+\Phi^2)&=&
\frac{\Omega}{2}\partial_+f\partial_+f\,,\no\\
- {\rm e}^{2\omega}\partial_-({\rm e}^{-2\omega}\partial_-\Phi^2)&=&
\frac{\Omega}{2}\partial_-f\partial_-f\,.
\label{2D-eom}
\end{eqnarray}
The energy-momentum tenser for the matter field $f$ is normalized as
\begin{eqnarray}
(T_{\rm matter})_{ab} &\equiv&
-\frac{2}{\sqrt{-g}}\frac{\delta S_{\rm matter}}{\delta g^{ab}}\no\\
 &=&-\frac{\Omega(\Phi)}{16\pi G}\left(\partial_af \,\partial_b f -\frac{1}{2}g_{ab}\,
\partial^c f \partial_c f\right)\,.
\label{tm}
\end{eqnarray}
This expression (\ref{tm}) is valid for the general form of $\Omega(\Phi)$\,.

\subsection{The AP model}

The AP model corresponds to a special case of 1+1 D dilaton gravity
specified by the following condition:
\begin{eqnarray}
U (\Phi) = 2 - 2 \Phi^2, \qquad \Omega(\Phi)=1\,. \label{AP}
\end{eqnarray}
This model exhibits nice properties. Among them, we are concerned with
the vacuum solution of this model. For our later convenience, we shall give a brief review
of the work \cite{AP} by focusing upon the vacuum solution in the following.

\medskip

The general vacuum solution is given by
\begin{eqnarray}
{\rm d} s^2 &=&\frac{1}{z^2}(-{\rm d} t^2 + {\rm d} z^2) \,,\\
\Phi^2 &=& 1+\frac{a+b \,t +c \,(-t^2+z^2)}{z}\,,
\label{dilaton-AP}
\end{eqnarray}
and depends on three real constants $a$\,, $b$ and $c$\,.
This three-parameter family contains interesting solutions as specific examples.
For example, the case with $a=1/2$\,, $b=0$ and $c=0$ corresponds to
a renormalization group flow solution from a conformal Lifshitz spacetime
to AdS$_2$ \cite{AP}, with an appropriate lift-up to higher dimensions.

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.3]{ads-undeform.eps}
\vspace*{-0.5cm}
\caption{\label{fig:AdS penrose} \footnotesize Penrose diagram of the black hole solution.
The global AdS$_2$ is parametrized by $\tau$ and $\theta$\,. The largest red triangle describes
the Poincar\'e patch of AdS$_2$ as usual. The coordinate system (\ref{undeformed BH1}) covers
the inside of smaller triangle bounded by green lines. The right vertex corresponds to
the black hole horizon which is specified as the point that $T$ is finite but $Z$ is infinity.}
\end{center}
\end{figure}

\medskip

Another intriguing example is a black hole solution specified with
$a=1/2$\,, $b=0$ and $c=\mu/2$\,, where $\mu$ is a real positive constant.
Then by performing a coordinate transformation,
\begin{eqnarray}
x^\pm=\frac{1}{\sqrt{\mu}}
\tanh\left(\sqrt{\mu}\, (T\pm Z)\right),
\label{to-BH}
\end{eqnarray}
the solution is rewritten into the following form:
\begin{eqnarray}
{\rm d} s^2 &=& \frac{4\, \mu}{\sin (2 \sqrt{\mu}\, Z)} (-{\rm d} T^2+{\rm d} Z^2)\,,\no\\
\Phi^2 &=& 1+ \sqrt{\mu}\,\text{coth}(2\sqrt{\mu}\, Z)\,.
\label{undeformed BH1}
\end{eqnarray}
The new coordinates $T$ and $Z$ cover a smaller region which is in the inside of
the entire Poincar\'e AdS$_2$\,, as depicted in Fig.\,\ref{fig:AdS penrose}.

\medskip

The background (\ref{undeformed BH1}) indeed describes a black hole geometry,
but it may not be so manifest.
To figure out the black hole geometry, it is nice to move to the Schwarzschild coordinates
by performing a further coordinate transformation,
\begin{eqnarray}
Z=\frac{1}{2\sqrt{\mu}}\text{arccoth}\left(\frac{\rho}{\sqrt{\mu}}\right)\,.
\end{eqnarray}
Then the background (\ref{undeformed BH1}) can be rewritten as\footnote{The factor 4 is included so that
the Bekenstein-Hawking entropy should match with the holographic computation.
This normalization guarantees the matching of the bulk and boundary times (or temperatures).
We are grateful to Ahmed Almheiri for this point. }
\begin{eqnarray}
{\rm d} s^2 = -4 (\rho^2-\mu){\rm d} t^2 + \frac{{\rm d} \rho^2}{\rho^2-\mu}\,, \qquad
\Phi^2 = 1+\rho\,.
\end{eqnarray}
In this metric, the black hole horizon is located at $\rho=\sqrt{\mu}$\,,
and the Hawking temperature $T_{\rm H}$ can be evaluated in the standard manner as
\begin{eqnarray}
T_{\rm H} =\left.\frac{1}{4\,\pi}\partial_\rho
\sqrt{\frac{-g_{tt}}{g_{\rho\rho}}}\right|_{\rho=\sqrt{\mu}} =\frac{\sqrt{\mu}}{\pi}\,.
\end{eqnarray}
Thus one can see that the background (\ref{undeformed BH1}) describes a black hole
whose horizon is located at $Z=\infty$\,.

\medskip

The Bekenstein-Hawking entropy can also be computed as
\begin{eqnarray}
S_{\rm BH} = \left.\frac{A}{4 G_{\rm eff}}\right|_{Z\to \infty}
= \left.\frac{\Phi^2}{4 G}\right|_{\sqrt{\mu}=\pi \, T_H}= \frac{1+\pi\,  T_{\rm H}}{4 G}\,.
\end{eqnarray}
Here the area $A$ is taken as $A=1$ because the horizon is just a point, and
the effective Newton constant $G_{\rm eff}$ can be read off
from the classical action as
\begin{eqnarray}
\frac{1}{G_{\rm eff}} = \frac{\Phi^2}{G} \,.
\label{Geff}
\end{eqnarray}

\medskip

On the other hand, the holographic entropy can be computed by using the renormalized
boundary stress tensor. For the detailed computation like the regularization and the counter term,
see \cite{AP}. As a result, the renormalized boundary stress tensor is evaluated as
 \begin{eqnarray}
\langle \hat{T}_{tt} \rangle = \frac{\mu}{8\pi G} \equiv E\,.
\end{eqnarray}
Then by using the thermodynamic relation
\begin{eqnarray}
{\rm d} S = \frac{{\rm d} E}{T_{\rm H}}\,, \label{thermo}
\end{eqnarray}
the entropy is obtained as
\begin{eqnarray}
S = \frac{\pi T_{\rm H}}{4 G} + S_{T_{\rm H}=0}\,,
\end{eqnarray}
where $S_{T_{\rm H}=0}$ is an integration constant.
Thus the holographic entropy agrees with the Bekenstein-Hawking entropy,
up to the temperature-independent constant.

\medskip

The main goal of this paper is to realize this correspondence of the entropies
for a deformed black hole solution introduced in the next section.

\section{Deforming the AP model}

In this section, we consider deforming the AP model so that
the deformed AdS$_2$ metric (\ref{3-metric}) is supported as a solution.
For simplicity, the matter fields are turned off hereafter. Along this line, as well as the dilaton itself,
the dilaton potential also has to be deformed from a simple quadratic one (\ref{AP})
to a hyperbolic function, similarly to integrable deformations.

\subsection{The deformed AP model}

\subsubsection*{The deformed metric}

Before discussing the dilaton and the dilaton potential,
it is helpful to rewrite the deformed metric (\ref{3-metric}) as
\begin{eqnarray}
{\rm d} s^2 &=& \frac{-{\rm d} t^2 + {\rm d} z^2}{z^2-\eta^2\left(\alpha + \beta t+\gamma (-t^2+z^2) \right)^2}
\nonumber \\
&=& \frac{1}{1-\eta^2 (X\cdot P)^2}\frac{-{\rm d} t^2+ {\rm d} z^2}{z^2}\,. \label{d-met}
\end{eqnarray}
Here we have introduced new quantities:  a coordinate vector $X^I$ and
a parameter vector $P_I$ defined as
\begin{eqnarray}
X^I &\equiv& \frac{1}{z}\left(t~,~\frac{1}{2}(1+t^2-z^2)~,~\frac{1}{2}(1-t^2+z^2)\right)\,,\no\\
P_I &\equiv&(\beta~,~\alpha-\gamma~,~\alpha+\gamma)\,\qquad (I,J=1,2,3)\,.
 \end{eqnarray}
The metric of the embedding space $\mathbb{M}^{2,1}$ is taken as $\eta_{IJ}=\text{diag}(-1,+1,-1)$\,.
The inner products are defined as
\begin{eqnarray}
&& X \cdot P \equiv X^I P_I = \frac{\alpha+\beta \, t+\gamma\,(-t^2+z^2)}{z}\,,
 \label{def-para}
\\ && X \cdot X \equiv \eta_{IJ}X^I X^J = -1\,, \qquad
P\cdot P \equiv \eta^{IJ}P_I P_J = -\tilde{\omega}\,.
\end{eqnarray}
These three products $X \cdot X$, $P \cdot P$ and $X\cdot P$ are transformed as scalars
under the $SL(2,\mathbb{R})$ transformation\footnote{
This $SL(2,\mathbb{R})$ transformation is the usual one generated by
three transformations,
1) time translation, 2) dilatation and 3) special conformal transformation.}.
For example, $X\cdot P$ is transformed as $X\cdot P= \tilde{X}\cdot \tilde{P}$,
where $\tilde{X}$ and $\tilde{P}$ are new coordinate and parameter vectors, respectively.

\medskip

Using the $SL(2,\mathbb{R})$ transformation, we can choose the vector $\tilde{P}$ freely
 as long as it satisfies the relation $\tilde{P}\cdot \tilde{P}=P\cdot P=-\tilde{\omega}$\,.
 Note that only the warped factor of the metric changes like
 \begin{eqnarray}
 {\rm d} s^2 &=&\frac{1}{1-\eta^2 (\tilde{X}\cdot \tilde{P})^2}
 \frac{-{\rm d}\tilde{t}^2+{\rm d}\tilde{z}^2}{\tilde{z}^2}\,
 \end{eqnarray}
because the rigid AdS$_2$ part is invariant under the $SL(2,\mathbb{R})$ transformation.

\subsubsection*{The dilaton sector}

Given the deformed metric (\ref{3-metric}) [or equivalently (\ref{d-met})]\,,
by solving the equations of motion (\ref{2D-eom}) without the matter fields,
the dilaton $\Phi$ is determined as
\begin{eqnarray}
\Phi^2 &=& \frac{c_1}{2\eta}\log\left|\frac{z+\eta\left(\alpha + \beta t
+ \gamma (-t^2+z^2) \right)}{z-\eta\left(\alpha + \beta t
+\gamma (-t^2+z^2) \right)}\right| +c_2 \nonumber \\
&\equiv &
\frac{c_1}{2\eta}\log\left|\frac{1+\eta (X\cdot P)}{1-\eta (X\cdot P)}\right| +c_2\,, \label{deformed-dilaton}
\end{eqnarray}
when the dilaton potential is deformed as
\begin{eqnarray}
U(\Phi)&=&
\left\{ \begin{array}{ll}\displaystyle
~-(1-\tilde{\omega}\, \eta^2)\frac{c_1}{\eta}\sinh\left[\frac{2\eta}{c_1}(\Phi^2-c_2)\right]
& ~~~\quad ({\text{for}}~~1>|\eta (X \cdot P)|)
\vspace*{0.3cm} \\ \displaystyle
~+(1-\tilde{\omega}\, \eta^2)\frac{c_1}{\eta}\sinh\left[\frac{2\eta}{c_1}(\Phi^2-c_2)\right]
& ~~~\quad  ({\text{for}}~~1<|\eta (X \cdot P)|)  \\
\end{array} \right.\,.
\label{d-pot}
\end{eqnarray}
Here $c_1$ and $c_2$ are arbitrary constants

\medskip

In the undeformed limit $\eta \to 0$\,, the dilaton (\ref{deformed-dilaton}) is reduced to
\[
\Phi^2 = c_2 + c_1\, \frac{\alpha + \beta\, t + \gamma (-t^2+z^2)}{z}\,,
\]
and thus the dilaton (\ref{dilaton-AP}) in the AP model has been reproduced
when $c_1=1$ and $c_2 =1$\,.  Remarkably, the three parameters $\alpha$\,, $\beta$ and $\gamma$
correspond to $a$\,, $b$ and $c$ in (\ref{dilaton-AP})\,, respectively.
Similarly, as $\eta \to 0$\,, the upper branch of the potential (\ref{d-pot}) reduces to
\[
U(\Phi) = 2(c_2 - \Phi^2)\,,
\]
while the lower branch vanishes.
Thus the dilaton potential of the AP model is reproduced when $c_2 =1$\,.
In total, the case with $c_1 = c_2 =1$ is associated with the AP model
and hence we will work with $c_1=c_2=1$ hereafter.

\subsubsection*{The vacuum solution in the deformed AP model}

In summary, the deformed AP model is specified by the deformed dilaton potential,
\begin{eqnarray}
U(\Phi)&=&
\left\{ \begin{array}{ll}\displaystyle
~-(1-\tilde{\omega}\, \eta^2)\frac{1}{\eta}\sinh\left[2\eta(\Phi^2-1)\right]
& ~~~\quad ({\text{for}}~~1>|\eta (X \cdot P)|)
\vspace*{0.3cm} \\ \displaystyle
~+(1-\tilde{\omega}\, \eta^2)\frac{1}{\eta}\sinh\left[2\eta(\Phi^2-1)\right]
& ~~~\quad  ({\text{for}}~~1<|\eta (X \cdot P)|)  \\
\end{array} \right.\,, \nonumber
\end{eqnarray}
and the vacuum solution is given by
\begin{eqnarray}
 {\rm d} s^2 = \frac{1}{1-\eta^2 (X\cdot P)^2}\frac{-{\rm d} t^2+ {\rm d} z^2}{z^2}\,, \qquad
 \Phi^2 = \frac{1}{2\eta}\log\left|\frac{1+\eta (X\cdot P)}{1-\eta (X\cdot P)}\right| +1\,,
\end{eqnarray}
where
\begin{eqnarray}
 X \cdot P = \frac{\alpha+\beta \, t+\gamma\,(-t^2+z^2)}{z}\,. \nonumber
\end{eqnarray}

\subsection{A deformed black hole solution}

In this subsection, we study a deformed black hole solution contained as a special case
of the general vacuum solution (obtained in the previous subsection).
This solution can be regarded as a deformation of the black hole solution presented
in \cite{AP}.

\medskip

In the following,  instead of $\tilde{\omega}$\,,
we use a new parameter $\mu$ defined as
\[
\mu \equiv - \tilde{P}\cdot \tilde{P} =\tilde{\omega} = -4c\,,
\]
so as to make our notation the same as that of the AP model. Here it may be worth noting
that the black hole temperature is related to the modification of the CYBE.
The zero temperature case corresponds to the homogeneous CYBE
and the temperature is measured by negative values of $c$\,. Solutions of the mCYBE
with negative (positive) $c$ are called the split (non-split) type. The well-known example
of the non-split type is the $q$-deformation of AdS$_5$ \cite{DMV2}, while the split type
has gotten little attention.  For the recent progress on the split type, see \cite{split}.
It may be interesting to seek some connection between black hole geometries and
solutions of split type.

\medskip

By performing the same coordinate transformation as in the undeformed case like
\begin{eqnarray}
x^\pm = \frac{1}{\sqrt{\mu}} \tanh(\sqrt{\mu}\, (T \pm Z))\,,
\end{eqnarray}
the deformed black hole solution is obtained as
\begin{eqnarray}
{\rm d} s^2 &=&\frac{4\mu}{-\eta^2\mu+(1-\eta^2\mu)\sinh^2(2\sqrt{\mu}\, Z)}
\bigl(-{\rm d} T^2+{\rm d} Z^2 \bigr)\,,\no\\
\Phi^2&=&
1+\frac{1}{2\eta}{\text{log}}\left|
\frac{1+\eta\sqrt{\mu}\, \coth(\sqrt{\mu}\,Z)}{1-\eta\sqrt{\mu}\,
\coth(\sqrt{\mu}\, Z)}\right|\,.
\label{d-BH}
\end{eqnarray}

\medskip

In this coordinate, the Ricci scalar (\ref{Ricci}) is rewritten as
\begin{eqnarray}
R=-(1-\eta^2 \mu)\frac{1-\eta^2 \mu -(1+\eta^2 \mu)\text{cosh}(4\sqrt{\mu} Z) }{\eta^2 \mu-(1-\eta^2 \mu)\,\text{sinh}^2 (2\sqrt{\mu} Z) }\,.
\end{eqnarray}
In the following, we impose that
\begin{eqnarray}
\eta^2 < \frac{1}{\mu}
\label{condition}
\end{eqnarray}
so as to ensure the existence of  the undeformed limit\footnote{
Otherwise, it is not posible to take the undeformed limit
$\eta\rightarrow 0$ because $\eta^2 >1/\mu$\,.}.
Note here that this background has a naked singularity at $Z=Z_0$\,, where
\begin{eqnarray}
Z_{0} \equiv \frac{1}{2\sqrt{\mu} }\, \text{arctanh}(\eta \sqrt{\mu})\,.
\end{eqnarray}
This is a peculiar feature of the Yang-Baxter deformed geometry
based on the modified CYBE like the $\eta$-deformation of AdS$_5$ \cite{ABF}.
From (\ref{Ricci}), in the region with $Z>Z_0$ the Ricci scalar takes negative values,
while for $0<Z<Z_0$, it has positive values (See Fig.\,\ref{fig:deformed AdS penrose})\,.
In the undeformed limit $\eta \rightarrow 0$\,, $Z_0$ is sent to zero
and the singularity disappears because the undeformed spacetime is just AdS$_2$\,.
In the following discussion, we focus upon the negative-curvature region ($Z>Z_0$)\,.
Therefore, we are concerned with only the upper branch of the potential (\ref{d-pot})\,.

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.3]{ads-deform.eps}
\vspace*{-0.5cm}
\caption{\label{fig:deformed AdS penrose} \footnotesize Penrose diagram of
the deformed black hole. In this diagram, a curvature singularity is depicted
in the global AdS$_2$ coordinates with $\alpha=1/2$, $\beta=0$ and $\gamma=\mu/2$
in (\ref{def-para}), where $\tau$ and $\theta$ are the same global coordinates as
in the undeformed AdS$_2$\,. The black curves represent the curvature singularities
of the deformed spacetime. In the blue region, the scalar curvature is positive,
while in the red and orange regions, it takes negative values.
The black hole coordinates in (\ref{d-BH}) covers the interior bounded by the green lines.
By employing a Schwarzschild-like coordinate system (\ref{d-Sch})\,,
we focus on the orange region in order to evaluate the black hole entropy.}
\end{center}
\end{figure}

\medskip

By performing the following coordinate transformation,
\begin{eqnarray}
r = \frac{1}{\eta} \text{arctanh}{ \bigl(\eta \sqrt{\mu}\, \coth \left(2 \sqrt{\mu}\, Z\right)
\bigr) }\,, \label{nice-coord}
\end{eqnarray}
the metric takes a Schwarzschild-like form\footnote{The factor 4
is included so as to reproduce the result of \cite{AP}. }
\begin{eqnarray}
{\rm d} s^2 &=& -4 F(r)\,{\rm d} T^2+\frac{ {\rm d} r^2}{ F(r) }\,,
\label{d-Sch}
\end{eqnarray}
where the scalar function $F(r)$ is defined as
\begin{eqnarray}
F(r)&\equiv& \frac{-1-\eta^2 \mu +(1-\eta^2 \mu)\text{cosh}(2\, \eta \, r)}{2\eta^2}\,.
\label{func-F}
\end{eqnarray}
In this coordinate system, the dilaton takes the simplest form,
\begin{eqnarray}
\Phi^2=1+r\,.
\end{eqnarray}
The locations of the boundary and black hole horizon are
\begin{eqnarray}
\text{boundary}\,: ~~r= \infty\,,\qquad  \text{BH horizon}\,: ~~
r=r^*\equiv \frac{1}{\eta}\text{arctanh}(\eta\sqrt{\mu})\,.
\end{eqnarray}

\subsubsection*{Bekenstein-Hawking entropy}

Let us compute the Bekenstein-Hawking entropy of the deformed black hole
by utilizing the coordinate system (\ref{d-Sch})\,.

\medskip

The Hawking temperature $T_{\rm H}$ is given by the standard formula:
\begin{eqnarray}
T_{\rm H} &=& \left. \frac{1}{4 \pi} \partial_{r} \sqrt{ -\frac{ g_{tt} }{ g_{rr} } }\right|_{r=r^*}
= \frac{ \sqrt{\mu} }{\pi}\,. \label{HT}
\end{eqnarray}
This is the same result as the undeformed case.
By assuming that the horizon area A is $1$ and using the effecting Newton constant
$G_{\rm eff}$ in (\ref{Geff})\,,
the Bekenstein-Hawking entropy $S_{\rm BH}$ can be computed as
\begin{eqnarray}
S_{\rm BH} = \left. \frac{A}{4 G_{\rm eff}}\right|_{r=r^*}
= \frac{\text{arctanh}{ (\pi \,T_{\rm H} \, \eta)} }{4 G \eta} + \frac{1}{4G}\,.
\label{BH}
\end{eqnarray}
In the undeformed limit $\eta \to 0$\,, the entropy is reduced to
\[
S_{\rm BH}^{(\eta=0)} = \frac{\pi T_{\rm H}}{4G} + \frac{1}{4G}\,,
\]
and thus the result of AP model has been reproduced.

\section{The boundary computation of entropy}

In this section, we compute the entropy of the deformed black hole
by evaluating the renormalized boundary stress tensor.
Now that the boundary structure is drastically changed, the first thing is to
determine the location of the holographic screen. In the following, we take the screen
on the singularity by following the proposal of \cite{Kameyama}. More precisely,
by introducing a UV cut-off $\epsilon$\,,  the boundary is taken just before
the singularity ($Z=Z_0+\epsilon$)\,.

\medskip

In the conformal gauge, the total action including the Gibbons-Hawking term
can be rewritten as
\begin{eqnarray}
S_{g,\Phi} &=&\frac{1}{16\pi G}\int\!{\rm d}^2 x\, \sqrt{-g}\left[\Phi^2 R - U(\Phi)\right] +\frac{1}{ 8 \pi G}
\int\!{\rm d} t\, \sqrt{-\gamma}\,\Phi^2 K\,\no\\
&=&\frac{1}{8 \pi G}\int\!{\rm d}^2 x \left[
-4 \partial_{(+} \Phi^2 \partial_{-)} \omega  - \frac{1}{2}U(\Phi)\, {\rm e}^{2 \omega}
\right]\,.
\end{eqnarray}
$K$ is the extrinsic curvature and $\gamma$ is the extrinsic metric.
By using the explicit expression of the deformed black hole solution in (\ref{d-BH}),
the on-shell bulk action can be evaluated on the boundary,
\begin{eqnarray}
S_{g,\Phi} &=& \int \! {\rm d} t\,  \left. \frac{ -\mu }{2\pi G ( 1+\eta^2 \mu +(-1+\eta^2 \mu )
 \cosh(4\mu^{ \frac{1}{2} } Z) ) } \right|_{Z\rightarrow Z_0}\,.
\end{eqnarray}
Recall that the regulator $\epsilon$ is introduced such that $Z-Z_{0}=\epsilon$\,,
the on-shell action can be expanded as
\begin{eqnarray}
S_{g,\Phi} &=& \int \! {\rm d} t\, \left[
\frac{1}{ 16 \pi G\, \eta\, \epsilon  } - \frac{1+\eta^2 \mu }{16 \pi G\, \eta^2}\,\epsilon^0
+ \frac{ (3+\eta^2 \mu (-2+3\eta^2 \mu) )\,\epsilon }{48 \pi G\, \eta^3} +O(\epsilon^2) \right]\,.
\end{eqnarray}

\medskip

To cancel the divergence that occurs as the bulk action approaches the boundary,
it is appropriate to add the following counter term:\footnote{The dual-theory interpretation
of it is not so clear
because it cotains an infinite number of polynomials and also depends
on the temperature explicitly. Another counter term may be allowed and it would be nice
to seek for it by following the procedure in \cite{Cvetic}.
We are grateful to Ioannis Papadimitriou for this point. }
\begin{eqnarray}
S_{\rm ct} &=& -\frac{1}{8 \pi G}\int\!{\rm d} t\, \sqrt{-\gamma_{tt}}\,
\sqrt{ F(\Phi^2-1) - \frac{1}{ \eta^2 } \log(1-\eta^2 \mu)}\,.
\label{counter}
\end{eqnarray}
Here the scalar function $F$ is already given in (\ref{func-F}) and hence
\begin{eqnarray}
 F(\Phi^2-1) &=& \frac{ -1-\eta^2 \mu +(1-\eta^2 \mu) \cosh( 2\eta (\Phi^2-1) )}{ 2 \eta^2 }
\end{eqnarray}
Note that the inside of the root of (\ref{counter}) is positive
due to the condition (\ref{condition})\,.
The extrinsic metric $\gamma_{tt}$ on the boundary is obtained as
\[
\gamma_{tt} =\left. - {\rm e}^{2 \omega} \right|_{Z \to Z_{0}} \,.
\]
In the undeformed limit $\eta \to 0$, this counter term reduces to
\begin{eqnarray}
S_{\rm ct}^{(\eta =0)} &=& \frac{1}{8 \pi G}\int\!{\rm d} t\, \sqrt{-\gamma_{tt}}\, (1-\Phi^2)\,,
\end{eqnarray}
because $\Phi^2-1>0$. This is nothing but the counter term utilized in the AP model \cite{AP}.

\medskip

It is straightforward to check that the sum $S = S_{g,\Phi}+S_{\rm ct}$ becomes finite
on the boundary by using the expanded form of the counter term (\ref{counter}):
\begin{eqnarray}
S_{\rm ct} &=& \int \! {\rm d} t\, \left[
 \frac{-1}{16 \pi G \eta\, \epsilon } + \frac{1+\eta^2 \mu
 +2 \log(1-\eta^2 \mu) }{ 16 \pi G\, \eta^2}\, \epsilon^0 +O(\epsilon) \right]\,.
\end{eqnarray}

\medskip

Around the boundary, the warped factor of the metric in (\ref{d-BH}) can be expanded as
\begin{eqnarray}
 {\rm e}^{2 \omega} &=& \frac{1}{ \eta\, \epsilon } - \left[ \frac{1}{ \eta^2} +\mu \right]\,\epsilon^0
 + O(\epsilon)\,.
\end{eqnarray}
Hence, by normalizing the boundary metric as
\[
\hat{\gamma}_{tt}  = \eta\, \epsilon\, \gamma_{tt}\,,
\]
the boundary stress tensor can be defined as
\begin{eqnarray}
\langle \hat{T}_{tt} \rangle
&\equiv& \frac{-2}{\sqrt{-\hat{\gamma}_{tt}}}\, \frac{\delta S}{\delta \hat{\gamma}^{tt} }
= \lim_{ \epsilon \to 0 } \sqrt{\eta\, \epsilon}\, \frac{-2}{\sqrt{-\gamma_{tt} } }\,
\frac{\delta S}{\delta \gamma^{tt} }\,.
\end{eqnarray}
After all, $\langle \hat{T}_{tt} \rangle$ has been evaluated as
\begin{eqnarray}
\langle \hat{T}_{tt} \rangle &=&-\frac{ \log(1-\eta^2 \mu) }{8 \pi G\, \eta^2}\,.
\end{eqnarray}

\medskip

To compute the associated entropy,
$\langle \hat{T}_{tt} \rangle$ should be identified with energy $E$ like
\begin{eqnarray}
E &=& -\frac{ \log(1-\pi^2 T_{\rm H}^2 \eta^2 ) }{8 \pi G\, \eta^2}\,,
\end{eqnarray}
where we have used the expression of the Hawking temperature (\ref{HT})\,.
Then by solving the thermodynamic relation (\ref{thermo}) again,
the entropy is obtained as
\begin{eqnarray}
S &=& \frac{ \text{arctanh}{(\pi T_{\rm H} \eta)} }{ 4 G \eta } + S_{T_{\rm H}=0}\,.
\end{eqnarray}
Here $S_{T_{\rm H}=0}$ has appeared as an integration constant that measures
the entropy at zero temperature.
Thus the resulting entropy precisely agrees with the Bekenstein-Hawking entropy
(\ref{BH})\,, up to the temperature-independent constant.

\medskip

Finally, it should be remarked that this agreement is quite non-trivial because
the deformation changes the UV region of the geometry drastically.
Hence the location of the holographic screen and the choice of the counter term
are far from trivial. Although the holographic screen was supposed to be the singularity,
inversely speaking, this agreement of the entropies here supports that the proposal
in \cite{Kameyama} would work well. As for the geometrical meaning of the counter term
(\ref{counter}), we have no definite idea. It is significant to figure out a systematic prescription
to produce the counter term (\ref{counter})\,.

\section{Conclusion and discussion}

In this paper, we have discussed deformations of the AP model
by following the Yang-Baxter deformation technique. To support the deformed
AdS$_2$ metric, the dilaton itself is deformed and the dilaton potential is also modified
from the polynomial to the hyperbolic function-type potential, similarly to integrable
deformations. We have obtained the general vacuum solution for the deformed potential.

\medskip

A particularly interesting example is a deformed black hole solution.
The deformation makes the spacetime structure
around the boundary change drastically and a new naked singularity appears.
The Hawking temperature is the same as in the undeformed case,
but the Bekenstein-Hawking entropy is modified due to the deformation.
This entropy has also been reproduced by evaluating the renormalized stress tensor
with an appropriate counter term on the regularized screen close to the singularity.

\medskip

There are some open problems. A possible generalization is to include matter fields,
though it has not succeeded yet. The matter contribution would not be so simple
in comparison to the AP model. It is also interesting to consider lifting up our results to
higher dimensional setups. Possibly, the most intriguing issue is to clarify the dual quantum mechanics
for the deformed black hole presented here. A candidate would be a deformed
SYK model which would be constructed by performing a disordered quench for a $q$-deformed
Heisenberg magnet. When an infinitesimal deformation of the deformed AdS$_2$ geometry is considered,
one would encounter a deformed Schwarzian derivative, though it seems difficult to determine
what it is because there is no $SL(2)$ invariance on the boundary
in comparison to the standard setup studied in \cite{Jensen,MSY,EMV}. It is also interesting to study
Yang-Baxter deformations of the Callan-Giddings-Harvey-Strominger (CGHS) model \cite{CGHS}
by following \cite{BKLSY,MORSY}.

\medskip

There are some future directions associated with Yang-Baxter deformations as well.
Now that we know the classical $r$-matrix which leads to the black hole geometry,
it would be interesting to consider a Yang-Baxter deformation of higher-dimensional AdS
with this $r$-matrix. In the study of Yang-Baxter deformations, it has been a long standing problem
to determine where the holographic screen is, while there was a proposal for the $\eta$-deformed
AdS$_5$ \cite{Kameyama} but it has not been supported by concrete evidence before this paper.
It is significant to find out more supports to clarify the holographic interpretation
for general Yang-Baxter deformations.

\medskip

We hope that the deformed AP model would provide a new arena to study
the correspondence between nearly AdS$_2$ geometries and 1D quantum mechanical system
like the SYK model or its cousins.

\subsection*{Acknowledgments}

We are very grateful to Ahmed Almheiri,
Heng-Yu Chen, Koji Hashimoto, Keiju Murata, Ioannis Papadimitriou,
Jun-ichi Sakamoto and Yuki Yokokura for useful comments and discussions.
The work of H.K. is supported by the Japan Society for the Promotion of Science (JSPS).
The work of K.Y. is supported by the Supporting Program for Interaction-based Initiative Team Studies
(SPIRITS) from Kyoto University and by a JSPS Grant-in-Aid for Scientific Research (C) No.\,15K05051.
This work is also supported in part by the JSPS Japan-Russia Research Cooperative Program
and the JSPS Japan-Hungary Research Cooperative Program.


\section{Introduction}\label{sec:intro}

In recent decades, functional data analysis (FDA) has received substantial attention and become increasingly important especially as the advent of ``Big Data'' era.
Representative monographs on FDA
include
\cite{Ramsay-Silverman05}, \cite{FerrV06}, \cite{HorvK12}, and \cite{HsinE15}.
Typically functional data are collected from $n$ curves $\{X_i: i=1, \ldots, n\}$ that are regarded as independent copies of a real-valued $L^2$
stochastic process $X$ defined on a compact domain $\mathcal{T}$ with mean function $\mu_0(t) = \mathbb{E} \{X(t)\}, t \in \mathcal{T},$ and covariance function $C_0(s, t) = \text{cov}\{X(s), X(t)\}, s, t \in \mathcal{T}$. In reality, due to discrete recording and the presence of noise, the data are often represented by
$\{(T_{ij},Y_{ij}): i=1,\dots, n; j=1,\dots m_i\}$,
where $m_i$ is the number of observations from the $i$-th curve $X_i$, and $Y_{ij} $ is
the noisy  observation from $X_i$ measured at the discrete time point
$T_{ij}$,
i.e., $Y_{ij} =
X_i(T_{ij}) +\varepsilon_{ij}$. Here $\{\varepsilon_{ij}: i=1,\dots, n;
j=1,\dots m_i\}$ are independent errors with zero mean and finite variance.
For simplicity and without loss of generality we assume $m_i=m$ for all $i$.

Among various population quantities, the covariance function $C_0$ is fundamental in FDA. Generally $C_0$ has two major roles. It is not only an important quantity that characterizes the temporal dependency \citep{Yao-Muller-Wang05, ZhanC07, Li-Hsing10, ZhanW16},
but also a building block for more advanced approaches in FDA
such as functional principal component analysis (FPCA) and functional linear regression \citep{YaoMW05, Hall-Muller-Wang06}.
Hence the FDA literature that involves covariance function estimation may accordingly be categorized into two types depending on the role of $C_0$.
As for the estimation of $C_0$, a variety of nonparametric methods have been proposed, such as local polynomial smoothing \citep{Li-Hsing10, ZhanW16},  B-splines \citep{JameHS00, RiceW01, Paul-Peng09}, penalized splines \citep{Goldsmith-Bobb-Crainiceanu11, XiaoLR13}, and smoothing splines \citep{Rice-Silverman91, Cai-Yuan10}.

Positive semi-definiteness is an essential characteristic of covariance functions.  Therefore, a valid covariance estimator is usually desired to be positive semi-definite, especially when this estimator is involved in subsequent analyses. See \cite{YaoMW05} and Section~\ref{sec:data} for examples. Meanwhile, it is also appealing if a covariance estimator is of low rank since this will encourage dimension reduction, alleviate computational and storage burdens, and facilitate simple interpretations.
In addition, low rank is often needed in trajectory prediction and some other advanced FDA methods \citep[e.g.,][]{Yao-Muller-Wang05, DelaH12,  LiWC13, JianAW16}.
Unfortunately, a majority of existing methods in FDA cannot directly produce a covariance estimator that is positive semi-definite or of low rank.
Hence a two-step procedure is typically performed in order to achieve at least one property, where a constraint-free covariance estimator is first obtained, then followed by
{a reconstruction step (e.g., via FPCA and truncation).}
See
\cite{Hall-Vial06} and \cite{Poskitt-Sengarapillai13}
for instances.
This two-step procedure, however, is unfavorable since it not only complicates the theoretical analysis of the final estimator, but also makes computation unstable due to the non-smooth truncation.

In this paper, we utilize
a reproducing kernel Hilbert space (RKHS) framework
to achieve a coherent ``one-step"
covariance estimation procedure such that the resulted estimator is automatically both positive semi-definite and of low rank.
The application of RKHS has gained popularity recently in FDA \citep[e.g.,][]{YuanC10, ZhuYZ14, WangR15}. In the same vein as penalized splines \citep[e.g.,][]{Pearce-Wand06} and smoothing splines
\citep[e.g.,][]{Wahba90, EggeL09, Gu13},
we suppose that the sample path of $X$ belongs to a RKHS $\mathcal{H}(K)$ defined on $\mathcal{T}$,
with a continuous and square-integrable reproducing kernel $K(\cdot, \cdot)$
defined on $\mathcal{T} \times \mathcal{T}$. A key property of $K$
is the so-called reproducing property:
\[
\langle K(t, \cdot), f(\cdot) \rangle_{\mathcal{H}(K)} = f(t), \quad \text{for any $t \in \mathcal{T}$ and $f \in \mathcal{H}(K)$}.
\]
Moreover, $K$ also uniquely determines the inner product and norm of $\mathcal{H}(K)$, denoted by $\langle \cdot, \cdot \rangle_{\mathcal{H}(K)}$ and $\|\cdot\|_{\mathcal{H}(K)}$ respectively.
A canonical example of RKHS is the $r$-th order Sobolev-Hilbert space on $\mathcal{T}=[0, 1]$:
\[
\mathcal{W}^r= \{g: g^{(v)}, v=0, \ldots, r-1, \,\text{are absolutely continuous};\, g^{(r)} \in L^2([0, 1])\},
\]
equipped with the squared norm
\[
\|g\|^2 = \sum^{r-1}_{v=0} \left\{ \int^1_0 g^{(v)}(t)dt\right\}^2 + \int^1_0 \left\{g^{(r)}(t)\right\}^2 dt.
\]
In this paper, we use $\mathcal{W}^2$ in all numerical implementations, but establish theoretical results for $\mathcal{W}^r, r \geq 2$, with general equivalent norms.

The RKHS framework was also used by \citet{Cai-Yuan10} for covariance function estimation, which is perhaps the most related
work to ours.
This will be made clear that their estimator is a non-positive semi-definite version of a special case in our general spectral regularization framework.
Under the assumption $\mathbb{E}\|X\|^2_{\mathcal{H}(K)}<\infty$,
they showed that
$C_0\in\mathcal{H}(K\otimes K)$, where $\mathcal{H}(K \otimes K)$ is the tensor
product RKHS
equipped with the norm $\| \cdot \|_{\mathcal{H}(K\otimes K)}$ and the reproducing kernel
\[
K \otimes K ((s_1, t_1), (s_2, t_2)) = K(s_1, s_2) K(t_1, t_2), \quad s_1,s_2,t_1, t_2\in\mathcal{T}.
\]
This suggests a tensor product RKHS modeling of $C_0$,
which we also adopt in this paper. With slight abuse of notation, we hereafter also use the notation $\otimes$ to denote the tensor product of functions, i.e., $f\otimes g(s, t) = f(s) g(t)$.

\citet{Cai-Yuan10} proposed
to estimate the covariance function $C_0$ by solving
\begin{equation}
  \underset{C\in\mathcal{H}(K\otimes K)}{\min} \left\{ \ell(C) + \lambda \| C\|^2_{\mathcal{H}(K\otimes K)} \right\},  \label{eqn:CY}
\end{equation}
where $\ell$ is a convex and smooth loss function characterizing the fidelity to the data,
and $\lambda > 0$ is a tuning parameter for the penalty term. Unfortunately, this approach cannot ensure the covariance estimator to be positive semi-definite or of low rank, so the aforementioned two-step procedure must be performed to improve the estimator.

We propose a new class of tensor product RKHS covariance estimators via a variety of spectral regularizations of an operator on $\mathcal{H}(K)$. The spectral regularizations generalize the penalty in (\ref{eqn:CY}), and can easily enable low-rank modeling, e.g., when the trace-norm penalty is used. The estimation framework respects the semi-positivity structure of covariance functions by imposing a constraint, so the resulted estimator automatically inherits this characteristic. Given any penalty, the covariance estimator is obtained by one step, which can reduce the computational and theoretical complexities of the two-step method. We establish a representer theorem to provide a finite dimensional representation for this class of covariance estimators, which makes the estimation procedure practically computable. Compared with its classical counterparts \citep[e.g.,][]{Wahba90, Cai-Yuan10}, the representer theorem is unconventional due to the semi-positivity constraint and a wide range of regularizations (e.g., trace-norm regularization). As a byproduct of the representer theorem, a closed form of the $L^2$ eigen-decomposition admitted by the covariance estimator can be easily obtained, without any numerical approximations needed in common FPCA approaches.

To promote dimension reduction, we particularly focus on trace-norm regularization
to additionally encourage low-rank estimation.
The corresponding objective function involved in the estimation framework is convex but non-differentiable. {An efficient algorithm is developed for this optimization problem} based on the representer theorem and the accelerated proximal gradient method \citep{Beck-Teboulle09}.
Note that, asymptotically, the use of trace-norm regularization does not rule out the
cases when $C_0$ is of high or infinite rank.
Irrespective of the true rank, our
estimator is consistent
with the optimal convergence rate, up to some order of $\log n$, as implied by the theoretical results below.

Despite the lack of a closed-form solution due to the
semi-positivity constraint and possibly non-differentiable penalties,
we develop the empirical $L^2$ rate of convergence for covariance estimators in the tensor product Sobolev-Hilbert spaces.
This result is broad since it allows for a variety of spectral regularizations, including the trace-norm regularization and others, and incorporates both fixed and random designs.
Generally, the rate is comparable to
the optimal rate of standard two-dimensional nonparametric smoothers. If $X$ is additionally
periodic, we can improve our results significantly
  such that the optimal one-dimensional nonparametric rate, up to some order of $\log n$, is attained.
For periodic functional spaces, when the data
are sparse, i.e., $m < \infty$, the rate of convergence is comparable to the
minimax rate obtained by \citet{Cai-Yuan10} and the $L^2$ rate achieved by
\citet{Paul-Peng09}.
Different from these two pioneer works, our objective function is not necessarily
  differentiable, which thus requires separate theoretical treatments.
Our theoretical results are established in terms of empirical processes techniques.
The success of the relevant proofs depends on the upper bound of the entropy for tensor
product Sobolev-Hilbert spaces, which is the first appearance in the FDA
literature to our best knowledge.

The rest of the paper is organized as follows.
The proposed methodology for covariance function estimation is presented in Section~\ref{sec:spec}. Computational issues and theoretical results are given in Sections~\ref{sec:computation} and~\ref{sec:asym} respectively. The empirical performance of the proposed approach is evaluated by a simulation study in Section~\ref{sec:sim} and a real data application in Section~\ref{sec:data}.
Additional materials, including technical details and further algorithmic descriptions, are provided in a separate supplemental document.

\section{Methodology}\label{sec:spec}

\subsection{Spectral decomposition on RKHS}\label{subsec:RKHS}

We first introduce spectral decomposition on RKHS and then define a variety of spectral regularizations which we will use to obtain a class of covariance function estimators.

For a bivariate function $C(\cdot, \cdot)$ on $\mathcal{T}\times \mathcal{T}$, define its transpose, denoted by $C^\top$, as $C^\top(s,t)=C(t,s)$ for any $s,t\in\mathcal{T}$.
Due to the symmetry of covariance functions, we focus on the space $\mathcal{S}(K)=\{C\in\mathcal{H}(K\otimes K): C = C^\top\}$.
For any $C\in\mathcal{S}(K)$,
define its self-adjoint operator $\mathcal{C}_C: \mathcal{H}(K)\rightarrow\mathcal{H}(K)$ by
\begin{equation}
  (\mathcal{C}_C f) (s) = \langle C(s,\cdot), f(\cdot) \rangle_{\mathcal{H}(K)}, \quad\quad \text{for any $f \in \mathcal{H}(K)$ and $s\in\mathcal{T}$.} \label{eqn:RKHSOp}
\end{equation}
Note that $\|C\|_{\mathcal{H}(K\otimes K)}<\infty$ since $C \in \mathcal{S}(K)$ and that the Hilbert-Schmidt norm of $\mathcal{C}_C$ coincides with $\|C\|_{\mathcal{H}(K\otimes K)}$. Therefore,
$\mathcal{C}_C$ is a Hilbert-Schmidt operator and hence
admits a spectral decomposition. In Section~\ref{subsec:spec}, we will define a penalty function based on this spectral decomposition.

In the FDA literature, the spectral analysis is often based on the Hilbert-Schmidt integral operator
$\mathcal{L}_C: L^2(\mathcal{T}) \rightarrow L^2(\mathcal{T})$ defined by
\begin{equation}
(\mathcal{L}_C f) (s) = \langle C(s,\cdot), f(\cdot)\rangle_{L^2(\mathcal{T})} =
\int_\mathcal{T} C(s,t) f(t)\, dt, \quad \text{for any $f \in L^2(\mathcal{T})$ and $s\in\mathcal{T}$.} \label{eqn:L2Op}
\end{equation}

There are two reasons why we adopt $\mathcal{C}_C$ instead of $\mathcal{L}_C$.
First, $\mathcal{C}_C$ is more aligned with the RKHS modeling of $X$, especially when the inner product of $\mathcal{H}(K)$ is chosen to mimic the physical reality.
Many examples can be found in the work on $L$-splines, e.g., Chapter~4.5 in \cite{Gu13} and Chapter~21 of \cite{Ramsay-Silverman05}.
Second, using $\mathcal{C}_C$ enables a finite dimensional representation
of our proposed covariance estimators as in Theorem \ref{thm:represent} below, and thus simplifies its practical computation.

\subsection{Spectrally regularized covariance estimator} \label{subsec:spec}

For any $C\in\mathcal{S}(K)$, let $\tau_1(C), \tau_2(C), \dots$ be the eigenvalues corresponding to the spectral decomposition of $\mathcal{C}_C$ such that $|\tau_1(C)|\ge|\tau_2(C)|\ge\cdots$.
We propose the following covariance estimator:
\begin{equation}
  \qquad \hat{C} = \argmin_{C\in\mathcal{S}^{+}(K)} \left\{ \ell(C) + \lambda \Psi(C) \right\},
  \label{eqn:ggoal}
  \end{equation}
where $\mathcal{S}^{+}(K)=\{ C\in\mathcal{S}(K): \langle\mathcal{C}_C f, f\rangle_{\mathcal{H}(K)} \ge 0,\ \text{for all $f\in\mathcal{H}(K)$}\}$,
$\ell$ is a convex and smooth loss function, $\lambda>0$ is a tuning parameter, and $\Psi(C) = \sum_{k\ge1} \psi(|\tau_k(C)|)$ with
$\psi$ being a non-decreasing penalty function satisfying $\psi(0)=0$ \citep{Abernethy-Bach-Evgeniou09}.
We assume that $\ell$ depends on $C$ through $\{C(T_{ij}, T_{ik}):i=1,\ldots,n; j, k =1, \ldots, m\}$.
The choice of $\psi$, and thus $\Psi$, is broad. In below we list a few interesting forms and briefly discuss their effects on the corresponding estimator.

\begin{ex}[Rank regularization] \label{ex:rank}
If $\psi(\tau) = I(\tau\neq 0)$ where $I(\cdot)$
    is the indicator function, $\Psi(C)$ is the rank of the operator $\mathcal{C}_C$. This penalty obviously encourages a low-rank solution.
    However, the minimization (\ref{eqn:ggoal}) is now difficult owing to its non-convexity, and
over-fitting may occur since no regularizations are imposed on non-zero eigenvalues.
\end{ex}
\begin{ex}[Hilbert-Schmidt-norm regularization] \label{ex:HS}
If $\psi(\tau)=\tau^2$, $\Psi(C)$ becomes the squared Hilbert-Schmidt norm of the operator $\mathcal{C}_C$, which equals $\|C\|_{\mathcal{H}(K\otimes K)}^2$.
Similar to the $\ell_2$-norm regularization for vectors, the Hilbert-Schmidt-norm regularization ensures the convexity of the objective function in (\ref{eqn:ggoal}), but does not encourage sparsity in eigenvalues, so the resulted covariance estimator is usually of high rank. \cite{Cai-Yuan10} used this regularization as in (\ref{eqn:CY}), but did not impose the constraint $C\in\mathcal{S}^{+}(K)$, so a positive semi-definite covariance estimator was not guaranteed.
\end{ex}

\begin{ex}[Trace-norm regularization] \label{ex:trace}
If $\psi(\tau) = \tau$,
    $\Psi(C)$ is the trace norm of $\mathcal{C}_C$,
    which is a convex relaxation of its rank.
    Similar to the celebrated $\ell_1$-regularization for vectors
    and the trace-norm regularization for matrices
    the trace-norm penalty $\Psi$ for operators
    not only
    promotes the sparsity of eigenvalues and hence
    low-rank solutions, but also
   regularizes non-zero eigenvalues.
    The minimization (\ref{eqn:ggoal}) now becomes a convex optimization which allows leveraging
    many recent developments in non-smooth convex optimizations
    \citep[e.g.,][]{Beck-Teboulle09}
    to achieve feasible computations. See Section \ref{sec:computation}
    for more
    details.
\end{ex}

Obviously the penalty $\Psi(C)$ generalizes the regularization as in (\ref{eqn:CY}). Regardless of the form of $\Psi(C)$, the covariance estimator obtained by (\ref{eqn:ggoal}) is always positive semi-definite since the solution to the minimization (\ref{eqn:ggoal}) is searched only within $\mathcal{S}^{+}(K)$.

\subsection{Representer theorem}\label{subsec:representer}

Since commonly used $\mathcal{H}(K)$, including $\mathcal{W}^r$, are infinite dimensional, solving (\ref{eqn:ggoal}) is typically an infinite dimensional optimization problem.
Therefore, $\hat{C}$ is of little
practical value if a finite dimensional representation, based on data, is unavailable.
To address this,
we provide a representer theorem which holds for the entire class
of estimators defined in (\ref{eqn:ggoal}).

Write $N=nm$ and $(\tilde{T}_1,\dots,\tilde{T}_N) =
(T_{11},\dots,T_{1m},T_{21}, \dots, T_{2m}, \dots, T_{n1},\dots, T_{nm})$.

\begin{theorem}[Representer theorem]\label{thm:represent}
If the solution set of (\ref{eqn:ggoal}) is not empty, then there always exists a solution
lying in the space
  $\mathcal{K}\otimes\mathcal{K}=\mathrm{span}\{  K(\cdot, \tilde{T}_i)\otimes K(\cdot,
    \tilde{T}_{j}): i,j=1,\dots,N \}$,
where $\mathcal{K} = \mathrm{span}\{ K(\cdot, \tilde{T}_i): i=1,\dots,N \}$.
Moreover, the solution takes the form:
\begin{equation}
C(s,t) = z(s)^\top {A}\, z(t),\label{eqn:represent2}
\end{equation}
where $A$ is a $N\times N$ symmetric matrix and $z(\cdot) = (K(\cdot,\tilde{T}_1),\dots, K(\cdot, \tilde{T}_N))^\top$.
\end{theorem}
\noindent

Classical representer theorems \citep[e.g.,][]{Wahba90}, as adopted in
\citet{Cai-Yuan10}, do not
cover the scenario addressed by Theorem \ref{thm:represent} due to
the semi-positivity constraint and
a wide choice of regularizations, e.g., the trace-norm regularization.
To show this theorem, we significantly utilize the fact that the spectral analysis is based on
the RKHS geometry.
This is the main reason for using the operator $\mathcal{C}_C$
in (\ref{eqn:RKHSOp}) instead of
$\mathcal{L}_C$ in (\ref{eqn:L2Op}).
We remark that the conclusion of Theorem \ref{thm:represent} also holds when the semi-positivity is not imposed, i.e., $\mathcal{S}^{+}(K)$ is replaced by
$\mathcal{S}(K)$ in (\ref{eqn:ggoal}).
In Section \ref{sec:sim}, this fact will be used to compute unconstrained estimators
for comparison.

At a first glance, a significant number of scalar parameters ($(N+1)N/2$) is involved in (\ref{eqn:represent2}). However, if a low-rank inducing penalty, such as the trace-norm regularization, is used,
the resulted estimator is often of low rank, which will benefit computation and storage in its estimation, and subsequent uses.
In Section~\ref{sec:computation}, we particularly focus on developing an efficient computational tool when the trace-norm regularization is imposed in (\ref{eqn:ggoal}).

\subsection{Parametrization}\label{subsec:para}
By Theorem \ref{thm:represent}, we are able to parametrize the solution to (\ref{eqn:ggoal}) in terms of a finite dimensional representation since
it suffices to merely focus on covariance functions of the form $C(\cdot, \cdot) = \sum^N_{i=1}\sum^N_{j=1} {A}_{ij} K(\cdot, \tilde{T}_i)
\otimes K(\cdot, \tilde{T}_j)$.
The eigenvalues of the operator $\mathcal{C}_C$, $\{\tau_j(C): j \ge 1\}$,
are the eigenvalues of the matrix
$B=M^\top{A} M$,
where
$M$ is any $N\times q$ matrix such that $MM^\top=\tilde{K}
=[K(\tilde{T}_i,
\tilde{T}_j)]_{1\le i,j\le N}$ with $q=\mathrm{rank}(\tilde{K})$.
The matrix $M$ provides a representation based on an orthonormal basis (of $\mathcal{K}$)
$\{v_1,\dots, v_q\}$:
\[
  K(\cdot, \tilde{T}_i) \otimes K(\cdot, \tilde{T}_j) = \sum^q_{k=1}\sum^q_{l=1} M_{ik}M_{jl}
  v_k \otimes v_l, \quad\quad 1\le i,j\le N.
\]
Therefore $C=\sum^q_{k=1}\sum^q_{l=1} B_{kl} v_k \otimes v_l$.
As $C(s,t)=\langle K(\cdot,s), \mathcal{C}_C K(\cdot, t)\rangle_{\mathcal{H}(K)}$,
we have
  $[C(\tilde{T}_i,\tilde{T}_j)]_{1\le i, j\le N} = MBM^\top$.
Moreover, write
$M=(M_1^\top, \dots, M_n^\top)^\top$,
where $\{M_i: i=1, \ldots, n\}$ are $m\times q$ matrices.
Then the loss function depends on $C$ through $[C(T_{ij}, T_{ik})]_{1\le j, k\le
  m}=M_i B M_i^\top$ for $i=1,\dots,n$.
Compared with $A$, the new parametrization $B$ is unique even
when $\{T_{ij}: i=1, \ldots, n; j=1, \ldots, m\}$ are not all unique.
Now (\ref{eqn:ggoal}) can be rewritten as
\begin{align}
\argmin_{B\in\mathcal{S}_q^{+}} \left\{\tilde{\ell}(B) + \lambda \tilde{\Psi}(B)\right\}, \label{eqn:goalB}
\end{align}
where
$\mathcal{S}_q^+$ is the set of all $q\times q$ positive semi-definite matrices,
$\tilde{\ell}(B) =
\ell(\sum^q_{k=1}\sum^q_{l=1} B_{kl} v_k \otimes v_l)$,
and $\tilde{\Psi}(B)=\sum_{k =1}^q\psi(|\xi_k(B)|)$ with $\xi_1(B),\dots, \xi_q(B)$
being the eigenvalues of the matrix $B$ such that $|\xi_1(B)|\ge\dots\ge|\xi_q(B)|$.
Conversely, with the new parametrization $B$, we can represent $C(s,t)=z(s)^\top (M^{+})^\top B M^{+} z(t)$ where
$M^+$ is the Moore-Penrose pseudoinverse of $M$.
Consequently, solving (\ref{eqn:ggoal}) is equivalent to solving (\ref{eqn:goalB}),
a finite-dimensional optimization.

\subsection[for the toc]{Closed-form expression of $L^2$ eigen-decomposition}\label{subsec:eigfun}
By Mercer's theorem, we can represent an arbitrary covariance function $C\in\mathcal{S}^{+}(K)$ in terms of the typical spectral
decomposition via the $L^2$ inner product, i.e.,
$C(s,t) = \sum_{k\ge1} \zeta_k \phi_k(s)\phi_k(t)$,
where $\{\phi_k:k\ge 1\}$ are the $L^2$ eigenfunctions and
$\{\zeta_k: k\ge 1\}$ are the corresponding $L^2$ eigenvalues.
This eigen-decomposition is a key component of FPCA and other advanced FDA methods. In the literature \citep[e.g.,][]{Rice-Silverman91}, approximate computations are commonly involved where the eigen-decomposition is obtained based on the discretized covariance function estimator.  In contrast, due to Theorem~\ref{thm:represent}, our covariance estimator posseses a closed-form expression of this eigen-decomposition
so that such computational complication can be avoided.

Following the notations in Sections~\ref{subsec:representer} and \ref{subsec:para}, let $Q = [\int_\mathcal{T} K(s, \tilde{T}_i)K(s,\tilde{T}_j )\, ds]_{1\le i, j\le N}=M R M^\top$ where
$R = [\int_\mathcal{T} v_k(s)v_l(s)\, ds]_{1\le k,l \le q}$.
Note that $R = M^+ Q (M^+)^\top$.
Similar to Lemma~3 of \citet{Cai-Yuan10}, once $\hat{B}$ is obtained from (\ref{eqn:goalB}),
the $L^2$ eigenfunctions of
$\hat{C}=\sum^q_{k=1}\sum^q_{l=1} \hat{B}_{kl} v_k \otimes v_l$
can be expressed as
  $\hat{\phi}_k(\cdot) = U^\top_k z(\cdot)$, for $k=1,\dots, n$,
where $U_k$ is the $k$-th column of $U=(M^+)^\top R^{-1/2} V$ and $V$ is the eigenvectors of
$R^{1/2} \hat{B} R^{1/2}$.
The $L^2$ eigenvalues of $\hat{C}$ coincide with those of $R^{1/2} \hat{B} R^{1/2}$, and the number of nonzero eigenvalues is the same as the rank of $\hat{B}$.

\section{Computational issues for trace-norm regularization}\label{sec:computation}

To achieve a desirable low-rank covariance estimator, we develop an algorithm when the trace-norm regularization is used.

\subsection{Algorithm}\label{subsec:alg}
With the trace-norm regularization in (\ref{eqn:ggoal}), it is equivalent to solving the convex optimization
\begin{equation}
\argmin_{B\in\mathcal{S}_q^+} \left\{\tilde{\ell}(B) + \lambda \|B\|_*\right\}, \label{eqn:goalnuclear}
\end{equation}
where $\|\cdot\|_*$ represents the typical trace norm for matrices.
We can also rewrite (\ref{eqn:goalnuclear}) as
\begin{equation}
\argmin_{B\in\mathcal{S}_q} \left\{\tilde{\ell}(B) + \lambda h(B)\right\},
\quad\mbox{where}\quad
h(B) =
\begin{cases} \|B\|_*, & B\in\mathcal{S}_q^+ \\ \infty, &
  B\not\in\mathcal{S}_q^+ \end{cases}.
\label{eqn:goalnuclear2}
\end{equation}
Here $\mathcal{S}_q$ represents the set of all $q\times q$ matrices.

The objective function in (\ref{eqn:goalnuclear2}) is the sum of a smooth function
and a non-smooth function.
A popular approach to such optimizations is the accelerated proximal gradient (APG) method \citep{Beck-Teboulle09}.
To apply this method, define an operator $\mathrm{svec} : \mathcal{S}_q \rightarrow \mathbb{R}^{q(q+1)/2}$ by
\[
  \mathrm{svec}(B) = [B_{11}, \sqrt{2}B_{21}, \dots, \sqrt{2}B_{q 1},
  B_{22}, \sqrt{2}B_{32}, \dots, \sqrt{2}B_{q 2},\dots,
  B_{q\ra}]^\top,
\]
for any $B=[B_{ij}]_{1\le i,j\le q}\in\mathcal{S}_q$.
This operator provides an isometry between $\mathcal{S}_q$ and $\mathbb{R}^{q(q+1)/2}$.
Denote its inverse by $\mathrm{svec}^{-1}$.
We write $\check{\ell}(b) = \tilde{\ell}(\mathrm{svec}^{-1}(b))$ for
any $b\in\mathbb{R}^{q(q+1)/2}$.
The APG algorithm of our case involves the proximal operator
$\mathrm{prox}_\nu: \mathbb{R}^{q(q+1)/2} \rightarrow \mathbb{R}^{q(q+1)/2}$ defined by
\begin{align*}
  \mathrm{prox}_{\nu} (b)
  &=\argmin_{d\in\mathbb{R}^{q(q+1)/2}} \left\{\frac{1}{2} \|d - b\|_E^2 + \nu h(\mathrm{svec}^{-1}(d))\right\}\nonumber\\
  &=\mathrm{svec}\left[\argmin_{D\in\mathcal{S}^+_q}\left\{ \frac{1}{2} \|D - B\|_F^2 + \nu\|D\|_*\right\}\right],
  \label{eqn:prox}
\end{align*}
for any $b\in\mathbb{R}^{q(q+1)/2}$ and $\nu>0$. Here $\|\cdot\|_E$ and $\|\cdot\|_F$
represent the Euclidean norm and the Frobenius norm respectively.
The following proposition states the closed-form solution of this proximal operator.

\begin{proposition}\label{prop:prox}
  For any $\nu>0$ and $b\in\mathbb{R}^{q(q+1)/2}$ with
  eigen-decomposition
  $\mathrm{svec}^{-1}(b)=P\mathrm{diag}(\tilde{b})P^\top$,
  \[\mathrm{prox}_\nu (b) = \mathrm{svec}(P\mathrm{diag}(\tilde{c}) P^\top),\]
where
$\tilde{c}=( g_\nu(\tilde{b}_{1}),\dots,g_\nu(\tilde{b}_{q}))^\top$ and
$\tilde{b}=(\tilde{b}_{1},\dots,\tilde{b}_{q})^\top$. Here $g_\nu(x) = (x-\nu)_{+}$ for any $x\in\mathbb{R}$.
\end{proposition}
Due to this closed-form solution,
we can avoid the application of an inner numerical optimization within every
iteration of the APG algorithm.
The proof uses the same technique as in the proof of Lemma~1 in
\citet{Mazumder-Hastie-Tibshirani10}, and is thus omitted.

The standard APG method requires the knowledge of the Lipschitz
constant of $\nabla\check\ell$,
which directly relates to the step size in each iteration of
the algorithm.
For many choices of the loss function $\ell$, the corresponding Lipschitz
constant of $\nabla\check{\ell}$ is difficult to obtain.
Moreover, even when the Lipschitz constant is known (e.g., the choice of $\ell$ described
in Section \ref{sec:ell}),
the algorithm usually suffers from conservative
step sizes \citep{Becker-Candes-Grant11}.
Hence the APG with backtracking steps is usually preferred.
Following the suggestions of \cite{Becker-Candes-Grant11}, we adopt a modified version of the APG
method. See details in
Algorithm \ref{alg:trace_back}.
Modifying Steps 7--8 will result in other variants of
proximal gradient methods. See Section 5.2 of \citet{Becker-Candes-Grant11}
for more discussions. For convergence properties of the APG method, we refer interested readers to \cite{Beck-Teboulle09}.

  \renewcommand{\baselinestretch}{1}
\begin{algorithm}[h]
  \SetKwInOut{Input}{Input}

   \caption{The APG algorithm with backtracking for trace-norm-regularized covariance estimation }\label{alg:trace_back}
   \Input{$B_0\in\mathcal{S}_q^+$, $\hat{L}>0$, $\eta> 1$, $\alpha<1$}
   \BlankLine
    $b_0\leftarrow\mathrm{svec}(B_0)$, $\bar{b}_0 \leftarrow b_0$, $\theta_{-1}\leftarrow +\infty$, $L_{-1} \leftarrow \hat{L}$\\
    \For{$k=0,1,2,\dots$}{
     $L_k \leftarrow \alpha L_{k-1}$\\

     \Repeat{convergence}{
     $\theta_{k}\leftarrow 2/[1+ \{1+ 4L_k/ (L_{k-1}\theta_{k-1}^2)\}^{1/2}]$\\
     $e_k \leftarrow (1-\theta_k) b_k + \theta_k \bar{b}_k$\\
     ${b}_{k+1} \leftarrow \mathrm{prox}_{\lambda/L_k}(e_k - \nabla \check\ell(e_k) / L_k)$\\
     $\bar{b}_{k+1}\leftarrow \{b_{k+1} - (1-\theta_k)b_{k}\}/\theta_k$\\
     $\hat{L}\leftarrow 2|(e_k-b_{k+1})^\top \{\nabla\check\ell (b_{k+1})-\nabla\check\ell (e_{k})\}|/\|b_{k+1}-e_k\|_E^2$\\
     \If{$L_k\ge \hat{L}$}{break}
   $L_k \leftarrow \max\{ \eta L_k, \hat{L} \}$
 }
    }
\end{algorithm}
  \renewcommand{\baselinestretch}{1.81}

A similar algorithm can be obtained for (\ref{eqn:ggoal}) coupled
with the Hilbert-Schmidt-norm regularization, which
will be implemented in Section~\ref{sec:sim} for a direct comparison.
Details are given in Section S1 of the supplemental document.

\subsection[for the toc]{A choice of $\ell$}\label{sec:ell}
Hereafter, we adopt the following {quadratic loss function}
\begin{equation}
  \ell(C) = \frac{1}{nm(m-1)} \sum^n_{i=1}\sum_{1\le j\neq k \le m}
  \left\{Z_{ijk} - C(T_{ij}, T_{ik})\right\}^2,  \label{eqn:lchoice}
\end{equation}
where
$Z_{ijk}=\{Y_{ij}-\hat{\mu}(T_{ij})\}\{Y_{ik}-\hat{\mu}(T_{ik})\}$ and
$\hat{\mu}$ is an estimator of the mean function $\mu_0$.
Since $[C(T_{ij},T_{ik})]_{1\le j,k\le m} = M_i B M_i^\top$ as shown in Section~\ref{subsec:para}, $\ell(C)$ becomes
\begin{align*}
  \tilde{\ell}(B) &= \frac{1}{nm(m-1)} \sum^n_{i=1}\| \rho(Z_i - M_i B M_i^\top)\|^2_F
  = \frac{1}{2} \mathrm{vec} (B)^\top \nabla^2 \tilde{\ell}(B) \mathrm{vec}(B) -
  \left(\sum^n_{i=1} M_i^\top Z_i M_i\right)^\top\mathrm{vec}(B),
\end{align*}
up to an additive constant independent of $B$.
Here $Z_i=[Z_{ijk}]_{1\le j,k\le m}$, $\rho$ is an operator setting the
diagonal entries of its input to zero,
and
\[
  \nabla^2 \tilde{\ell}(B) = \frac{2}{nm(m-1)} \sum^n_{i=1} (M_i^\top \otimes M_i^\top)
  \mathrm{diag}\{\mathrm{vec}(\tilde{I})\} (M_i \otimes M_i),
\]
with $\tilde{I}\in \mathbb{R}^{q\times q}$ consisting of elements $\tilde{I}_{ij} = I(i\neq j)$.
With straightforward derivations, one can obtain the closed-form expressions
of $\check{\ell}$ and $\nabla\check\ell$ as required in Algorithm \ref{alg:trace_back}, which we omit here.

\section{Asymptotic properties}\label{sec:asym}
In this section, we develop the empirical $L^2$ rate of convergence for a variety of spectrally regularized covariance estimators in the tensor product Sobolev-Hilbert spaces. This result is broad since it incorporates both fixed and random designs, and also allows for a variety of spectral regularizations, including the trace-norm, Hilbert-Schmidt regularizations and others.
\subsection{Assumptions} \label{subsec:assump}
Without loss of generality, we take $\mathcal{T} = [0, 1]$. Below we establish
the asymptotic properties for the $r$-th order Sobolev-Hilbert space on $[0, 1]$ where $r \ge 2$, i.e.,
\[
\mathcal{H}(K) = \{g: g^{(v)}, v=0, \ldots, r-1, \,\text{are absolutely continuous};\, g^{(r)} \in L^2([0, 1])\}.
\]
The space $\mathcal{H}(K)$ is equipped with squared norm $\|g\|^2= \sum^{r}_{v=0} \int^1_0 \{g^{(v)}(t)\}^2\,dt $.
The asymptotic results also hold for its equivalent norms, e.g., $\|g\|^2= \int^1_0 \{g(t)\}^2\, dt + \int^1_0 \{g^{(r)}(t)\}^2\, dt$, $\|g\|^2=( [\int^1_0 \{g(t)\}^2\, dt]^{1/2} + [\int^1_0 \{g^{(r)}(t)\}^2\, dt]^{1/2} )^2$, and $\|g\|^2= \sum^{r-1}_{v=0} \{ \int^1_0 g^{(v)}(t)\,dt\}^2 + \int^1_0 g^{(r)}(t)^2\, dt$.

We list the assumptions needed for the asymptotic properties as follows.

 \begin{assumption} \label{assum:C0}
 $C_0 \neq 0$ and $C_0 \in\mathcal{F}
 \subseteq \mathcal{H}(K\otimes K)$ where $\mathcal{F}$ is the hypothesis space for estimation.
 \end{assumption}
  \begin{assumption} \label{assum:T}
The time points $\{T_{ij}: i=1,\ldots, n; j=1,\ldots, m\}$ are either fixed or random, and are independent of $\{X_i: i=1, \ldots, n\}$. The errors $\{\varepsilon_{ij}: i=1,\dots, n; j=1,\dots m\}$ are independent of both $\{T_{ij}: i=1,\ldots, n; j=1,\ldots, m\}$ and $\{X_i: i=1, \ldots, n\}$.
 \end{assumption}
 \begin{assumption} \label{assum:X}
 For each $t\in [0, 1]$, $X(t)$ is sub-Gaussian with a parameter $b_X>0$ which does not depend on $t$, i.e., $\mathbb{E}(\exp\{\beta X(t)\})\le \exp\{b_X^2 \beta^2/2\}$ for all $\beta>0$ and $t\in [0, 1]$.
 \end{assumption}
 \begin{assumption} \label{assum:varesp}
   For each $i, j$, $\varepsilon_{ij}$ is sub-Gaussian with a parameter $b_\varepsilon$ independent of $i$ and $j$.
 \end{assumption}

As shown in \citet{Cai-Yuan10},
$C_0\in\mathcal{H}(K\otimes K)$
under the assumptions that $X\in\mathcal{H}(K)$ almost surely and
$\mathbb{E}\|X\|^2_{\mathcal{H}(K)}<\infty$.
Assumption~\ref{assum:T} is standard in FDA modeling.
Assumptions~\ref{assum:X} and~\ref{assum:varesp} are sub-gaussian conditions
of the stochastic process and the measurement error.

\subsection{Rate of convergence}\label{subsec:rate}
We investigate the asymptotic property of a class of covariance estimators given by
\begin{equation}\label{eqn:ggoalnew}
\hat{C}_{\lambda} = \argmin_{C\in \mathcal{F}} \left\{ \ell(C) + \lambda \Psi(C) \right\},
\end{equation}
where $\Psi(C) =\sum_{k\ge1} |\tau_k(C)|^p$ for $1 \leq p \leq 2$
and the loss function $\ell$ is chosen as (\ref{eqn:lchoice}).
Apparently, the penalty term $\Psi$ generalizes both trace-norm ($p=1$) and Hilbert-Schmidt-norm ($p=2$) regularizations. Moreover, $\hat{C}_{\lambda}$ becomes the estimator by \citet{Cai-Yuan10} if $\mathcal{F} = \mathcal{H}(K\otimes K)$ and $p=2$.
For simplicity, we assume known $\mu_0=0$ so we let $\hat{\mu}=0$ and accordingly $Z_{ijk}=Y_{ij} Y_{ik}$.

For arbitrary bivariate functions $g_1$ and $g_2$, define an empirical inner product and the corresponding empirical norm as follows:
\[
\ip{g_1, g_2}_n = \frac{1}{nm(m-1)}\sum_{i=1}^n \sum_{1 \leq j \neq k \leq m} g_1(T_{ij}, T_{ik})g_2(T_{ij}, T_{ik}) \quad\mbox{and}\quad \norm{g_1}_n^2 = \ip{g_1, g_1}_n.
\]
Recall that we say a random variable $S_n=\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_{p}(k_n)$ if
    \[
    \lim_{L\rightarrow\infty} \limsup_{n \rightarrow \infty}\, \mathrm{Pr}(S_n\ge L k_n) =0.
    \]
  To accommodate the flexibility of the design $\mathbb{T}=\{T_{ij}: i=1, \ldots, n; j=1, \ldots, m\}\in\mathcal{T}^{nm}$, we denote $S_n=\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_{p}^T(k_n)$ if
  \[
    \lim_{L\rightarrow\infty} \limsup_{n \rightarrow \infty}
    \sup_{\mathbb{T}\in\mathcal{T}^{nm}}\, \mathrm{Pr}(S_n\ge L k_n \mid \mathbb{T}) =0.
  \]

We first provide the empirical $L^2$ rate of convergence for $\hat{C}_{\lambda}$.
\begin{theorem}\label{thm:2drate}
Under Assumptions~\ref{assum:C0}--\ref{assum:varesp},
if $\Psi(C_0) > 0$ and $\lambda^{-1} = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p \left\{ n^{r/(1+r)}\right\}$,
we have
$\norm{\hat{C}_{\lambda} - C_0}_n =\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p( \lambda^{1/2} )$.
Further, if
$\lambda^{-1} = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p^T \left\{ n^{r/(1+r)}\right\}$,
we have
$\norm{\hat{C}_{\lambda} - C_0}_n =\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p^T( \lambda^{1/2} )$.
\end{theorem}
In Theorem \ref{thm:2drate}, the asymptotic accuracy of $\hat{C}_\lambda$ is guaranteed for both fixed and random designs. In particular, both independent and dependent designs are also allowed if the design is random. Furthermore, Theorem \ref{thm:2drate} provides a uniform result over all designs under a stronger condition of $\lambda$. For instance, such $\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p^T$-condition degenerates to the weaker $\ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p$-condition if the choice of $\lambda$ is nonrandom or independent of the design.

Theorem~\ref{thm:2drate} incorporates a variety of regularizations as long as
$1 \leq p \leq 2$, where the commonly used trace-norm ($p=1$) and
Hilbert-Schmidt-norm ($p=2$) penalties are both special cases. It shows that
the empirical $L^2$ rate of convergence of $\hat{C}_{\lambda}$ is comparable to
that of standard two-dimensional nonparametric smoothers. For example, the rate
of convergence is $n^{1/3}$ for the second order Sobolev-Hilbert space, i.e.,
$r=2$.
The conclusion in Theorem~\ref{thm:2drate} is generally true for all two-dimensional Sobolev spaces, but the rate is sub-optimal within the scope of tensor product Sobolev-Hilbert spaces.
For periodic functions, however, we are able to significantly improve this rate by utilizing appropriate and specific entropy results for tensor product Sobolev-Hilbert spaces.

\begin{theorem}\label{thm:rate}
Suppose that $\mathcal{F} \subseteq \{ C \in \mathcal{H}(K\otimes K): \text{$C$ is a periodic function}\}$. Under Assumptions~\ref{assum:C0}--\ref{assum:varesp},
if $\Psi(C_0) > 0$, and
$\lambda^{-1} = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p \{ {n^{2r/(1+2r)}}/{\log n}\}$,
we have
$\norm{\hat{C}_{\lambda} - C_0}_n = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p( \lambda^{1/2} )$.
Further, if
$\lambda^{-1} = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p^T \{ {n^{2r/(1+2r)}}/{\log n}\}$,
we have
$\norm{\hat{C}_{\lambda} - C_0}_n = \ensuremath{\mathop{}\mathopen{}\mathcal{O}\mathopen{}}_p^T( \lambda^{1/2} )$.
\end{theorem}

Similar to Theorem~\ref{thm:2drate}, Theorem~\ref{thm:rate} also allows for both fixed and random designs.
Theorem~\ref{thm:rate} demonstrates that $\hat{C}_{\lambda}$ can achieve the empirical $L^2$ rate of convergence for one-dimensional nonparametric estimation, up to some order of $\log n$, although the target function $C_0$ is two-dimensional.
For instance, if we let $r=2$, the rate of $\hat{C}_{\lambda}$ is $(\log n)^{-1/2} n^{2/5}$, which is much faster than the two-dimensional nonparametric rate $n^{1/3}$.
For sparse functional data, i.e., $m < \infty$, up to some order of $\log n$, the rate of $\hat{C}_{\lambda}$ is comparable to the minimax rate obtained by \citet{Cai-Yuan10} and the $L^2$ rate achieved by \citet{Paul-Peng09} for $r=4$. However, the rates in both theorems are sub-optimal for functional data that are not sparse \citep{ZhanW16}.

The covariance estimator $\hat{C}_{\lambda}$ defined in (\ref{eqn:ggoalnew}) does not have a closed form due to the possible non-differentiability of the penalty term (e.g., when $p=1$), and the flexibility of $\mathcal{F}$. This explains the technical challenges and highlights the novelties of the proofs for Theorems~\ref{thm:2drate} and~\ref{thm:rate}.
In Theorem~\ref{thm:rate}, the particular structure of the tensor product RKHS accounts for the appealing rate of convergence of $\hat{C}_{\lambda}$. The upper bound of the entropy for tensor product Sobolev-Hilbert spaces, as given in Lemma~1 of the supplemental document, is a crucial component for the technical success. To our best knowledge, this paper is the first one in the FDA literature that achieves this result.

\section{Simulation experiments}\label{sec:sim}

Numerical experiments were conducted to illustrate the practical performance of
the proposed methodology.
We generated $\{X_i: i=1, \ldots, n\}$ where $n=200$ from a Gaussian process with $\mu_0(t) = 3\sin\{3\pi (t + 0.5)\} + 2 t^3$
and $C_0(s, t) = \sum_{k=1}^L (k+1)^{-2} \phi_k(s) \phi_k(t)$,
with $\phi_1(t) = 2^{1/2} \cos(2 \pi t)$, $\phi_2(t) = 2^{1/2} \sin(2 \pi t)$, $\phi_3(t) = 2^{1/2} \cos(4 \pi t)$, and $\phi_4(t) = 2^{1/2} \sin(4 \pi t)$.
We also sampled $\{T_{ij}: i=1,\dots, n; j=1,\dots m\}$ independently from the uniform distribution on $[0, 1]$ and $\{\varepsilon_{ij}: i=1,\dots, n; j=1,\dots m\}$ independently from $N(0, 0.01)$ to produce $Y_{ij} = X_i(T_{ij}) + \varepsilon_{ij}$.
We studied six settings in total, where $L=2$ or $4$, and $m=5$, $10$ or $20$.
In each setting, we simulated 300 datasets
where we compared various covariance function estimators.
Other than our proposed estimators, we also included popular alternatives, including the
covariance smoothing estimators by local polynomial regression \citep{Yao-Muller-Wang05} and bivariate P-splines \citep{Goldsmith-Bobb-Crainiceanu11} respectively.
In these common alternatives, a raw smoothed estimate is first computed. Then
a truncation step via FPCA is often applied to reconstruct a covariance function
that is both positive semi-definite and of low rank.
That means, the reconstructed covariance estimator, which we refer to as a two-step estimator below,
takes the form: $\sum_{k=1}^J \hat{\zeta}_k\hat{\phi}(s)\hat{\phi}_k(t)$ where $\hat{\zeta}_k$'s and $\hat{\phi}_k$'s are the largest $J$ positive estimated eigenvalues and corresponding estimated eigenfunctions based on the raw smoothed estimate.

Altogether we compared the following ten methods, {of which the first five are based on our proposed framework while the rest are popular alternatives}:
1) $\hat{C}^{+}_{\mathsf{trace}}$: {obtained from} (\ref{eqn:ggoal}) with trace-norm regularization;
2) $\hat{C}_{\mathsf{trace}}$: {obtained from} (\ref{eqn:ggoal}) with trace-norm regularization but {without}
the semi-positivity constraint, i.e., $\mathcal{S}^+(K)$ replaced by $\mathcal{S}(K)$;
3) $\hat{C}^{+}_{\mathsf{HS}}$: {obtained from} (\ref{eqn:ggoal}) with Hilbert-Schmidt norm regularization;
4) $\hat{C}_{\mathsf{HS}}$: {obtained from} (\ref{eqn:ggoal}) with
    Hilbert-Schmidt norm regularization but {without}
    the semi-positivity constraint, which is exactly (\ref{eqn:CY});
5) $\hat{C}_{\mathsf{CY}}$: the estimator proposed by \cite{Cai-Yuan10} as defined by (\ref{eqn:CY}) and
    implemented in their R package\footnote{Downloaded from \url{http://stat.wharton.upenn.edu/~tcai/paper/html/Covariance-Function.html}};
6) $\hat{C}_{\mathsf{PACE}}$: the raw smoothed covariance estimator using local polynomial regression
\citep{Yao-Muller-Wang05} implemented in the R package \textsf{fdapace};
7) $\hat{C}_{\mathsf{PACE,FVE}}^{+}$: the two-step estimator based on $\hat{C}_{\mathsf{PACE}}$ with the
number of components $J$ selected by fraction of variation explained (\textsc{fve})\footnote{With package default threshold 0.9999}, implemented in \textsf{fdapace};
8) $\hat{C}_{\mathsf{PACE,AIC}}^{+}$: the two-step estimator based on $\hat{C}_{\mathsf{PACE}}$ with
$J$ selected by Akaike Information Criterion (\textsc{aic}), implemented in \textsf{fdapace};
9) $\hat{C}_{\mathsf{PACE,BIC}}^{+}$: the two-step estimator based on $\hat{C}_{\mathsf{PACE}}$ with
$J$ selected by Baysian Information Criterion (\textsc{bic}), implemented in \textsf{fdapace};
10) $\hat{C}_{\mathsf{SC}}^{+}$: the two-step estimator based on tensor product
bivariate P-spline smoothing \citep{Goldsmith-Bobb-Crainiceanu11}
with $J$ selected by \textsc{fve}\footnote{With package default threshold 0.99},
implemented in the R package \textsf{refund}.
Note that \textsf{refund} does not output the raw smoothed estimate of $\hat{C}_{\mathsf{SC}}^{+}$ and therefore we do not compare such raw estimator.

{To obtain $\hat{C}^{+}_\mathsf{trace},\hat{C}_\mathsf{trace},\hat{C}^{+}_\mathsf{HS},\hat{C}_\mathsf{HS}$
and $\hat{C}_\mathsf{CY}$,}
smoothing spline was first applied to estimate $\mu$, where its smoothing parameter was selected by generalized cross-validation (\textsc{gcv}).
To obtain each of these five covariance estimators, we always used the loss function $\ell$ as in Section~\ref{sec:ell}, and $\mathcal{H}(K)$ was chosen as the second order Sobolev-Hilbert space on
$[0,1]$ with the squared norm
$\|g\|^2 = \sum^{1}_{v=0} \{ \int^1_0 g^{(v)}(t)dt\}^2 + \int^1_0 \{g^{(2)}(t)\}^2 dt$.
The tuning parameter $\lambda$ of the first four methods were chosen by
five-fold cross-validation.
The computations of $\hat{C}_{\mathsf{trace}}$, $\hat{C}^+_{\mathsf{HS}}$ and $\hat{C}_{\mathsf{HS}}$
were achieved by Algorithm \ref{alg:trace_back}
with different proximal operators (Step 7) due to the change of penalty and
the utility of positivity constraint.
For the remaining five methods, $\mu$ is estimated by the corresponding computational packages. {See their documents for further implementation details.}

Table \ref{tab:rkhscov} shows the average integrated squared errors (\textsc{aise})
and
average ranks of these covariance estimators over 300 simulated data sets.
First, we restrict our attention to the first five methods which can all be regarded as variations of (\ref{eqn:ggoal}).
Although $\hat{C}_{\mathsf{CY}}$ and $\hat{C}_{\mathsf{HS}}$ share the same definition (\ref{eqn:CY}), they differ in various implementation details
and hence
the practical performance.
However, their differences in \textsc{aise}
{are too small to}
affect the subsequent {comparisons in this section, so hereafter it suffices to include only $\hat{C}_{\mathsf{HS}}$, rather than both of them, to study rank reduction and the effect of the semi-positivity constraint. }

When we compare the two pairs, $\hat{C}^{+}_\mathsf{trace}$ versus $\hat{C}_\mathsf{trace}$, and $\hat{C}^{+}_\mathsf{HS}$ versus $\hat{C}_\mathsf{HS}$, obviously the covariance estimators with the positivity constraint always achieve smaller \textsc{aise} values than their counterparts.
This suggests that not only can imposing the semi-positivity constraint produce a valid estimator, but also improve estimation accuracy.
Notice in Table \ref{tab:rkhscov} that rank reduction can also be observed for
$\hat{C}^{+}_\mathsf{HS}$ because the semi-positivity constraint often results
in truncations of eigenvalues at zero.
When $\hat{C}^{+}_\mathsf{trace}$ is compared with $\hat{C}^{+}_\mathsf{HS}$, the former performs slightly worse in \textsc{aise}, but significantly better in rank reduction. The average ranks of $\hat{C}^{+}_\mathsf{trace}$ are in fact the closest to the true rank among the three rank-reduced estimators.
This highlights the benefits of trace-norm regularizations in computation, storage and subsequent uses as mentioned in Section~\ref{sec:intro}. {Next we only compare $\hat{C}^{+}_\mathsf{trace}$ with the five popular alternatives.}
{The three two-step estimators based on local polynomial regression, $\hat{C}_{\mathsf{PACE,FVE}}^{+}$, $\hat{C}_{\mathsf{PACE,AIC}}^{+}$ and $\hat{C}_{\mathsf{PACE,BIC}}^{+}$, have similar average ranks, which explains why their \textsc{aise} values are almost identical. When compared with the raw smooth estimator $\hat{C}_{\mathsf{PACE}}$, they can all reduce rank and achieve smaller \textsc{aise} values, which illustrates the benefit of the reconstruction step of the classical two-step procedure.
However, their average ranks and  \textsc{aise} values are both much higher than those of $\hat{C}_{\textsf{trace}}^{+}$ in all settings. Thus our proposed $\hat{C}_{\textsf{trace}}^{+}$ uniformly and significantly outperforms all these estimators with respect to both rank reduction and estimation accuracy.}

{When compared with $\hat{C}_{\mathsf{SC}}^{+}$, $\hat{C}_{\textsf{trace}}^{+}$ achieves similar \textsc{aise} values and performs slightly but uniformly better in rank reduction in all six settings. Table \ref{tab:rkhscov} also shows that $\hat{C}_{\textsf{trace}}^{+}$ is numerically more stable than $\hat{C}_{\mathsf{SC}}^{+}$. Computational errors occurred in some simulation runs to obtain $\hat{C}_{\mathsf{SC}}^{+}$, but none appeared for $\hat{C}_{\textsf{trace}}^{+}$. Therefore, $\hat{C}_{\textsf{trace}}^{+}$ is a more desirable covariance estimator than $\hat{C}_{\mathsf{SC}}^{+}$.}

{In summary, the overall performance of $\hat{C}_{\textsf{trace}}^{+}$ is the best among the ten covariance estimators regarding rank reduction, estimation accuracy, and computational stability. This motivates us to use $\hat{C}^{+}_\mathsf{trace}$ in the following real data application.}

\renewcommand{\baselinestretch}{1}
\begin{table}[ht]
\centering
{\small  \caption{\textsc{aise} ($\times 10^3$) values with standard errors ($\times 10^3$) in parentheses for the {ten} covariance estimators,
and average ranks for those estimators with rank reduction. {For $\hat{C}^{+}_{\mathsf{SC}}$, whenever computational errors occurred in a setting, its statistics are computed only based on successful runs, with the proportion of success runs additionally shown in square brackets.}
  }
  \label{tab:rkhscov}
\begin{tabular}{rrr|rrrrr}
    \hline\hline
$m$ & $L$ &   & $\hat{C}^{+}_\mathsf{trace}$ & $\hat{C}_\mathsf{trace}$ &
    $\hat{C}^{+}_\mathsf{HS}$ & $\hat{C}_\mathsf{HS}$ & $\hat{C}_\mathsf{CY}$ \\\hline
5 & 2 & \textsc{aise} &  4.86 (0.168) &  7.64 (0.159) &  4.77 (0.159) &  8.22 (0.171) &  6.43 (0.146) \\
  && rank &2.77 & 6.14 & 14.8 & - & - \\
  & 4 & \textsc{aise} &  9.66 (0.262) & 14.98 (0.244) &  9.15 (0.237) & 15.54 (0.220) & 12.26 (0.293) \\
  && rank  & 3.90 & 7.28 & 15.3 & - & - \\
  \hline
10 & 2 & \textsc{aise} &  2.85 (0.0918) &  3.56 (0.0953) &  2.81 (0.0888) &  4.08 (0.0904) &  4.30 (0.0825) \\
 && rank  & 2.74 & 8.00 & 14.5 & - & - \\
& 4 & \textsc{aise} &  4.94 (0.107) &  7.97 (0.123) &  4.78 (0.0991) &  8.25 (0.125) &  7.25 (0.104) \\
&&rank  & 4.44 & 12.9 & 14.9 & - & - \\
  \hline
20 & 2& \textsc{aise} &  2.07 (0.0777) &  2.27 (0.0782) &  2.06 (0.0772) &  2.54 (0.0782) &  3.46 (0.0787) \\
 && rank  & 2.68 & 8.47 & 14.3 & - & - \\
 & 4 & \textsc{aise} &  3.27 (0.0809) &  4.69 (0.0846) &  3.24 (0.0808) &  4.69 (0.0853) &  6.04 (0.0764) \\
&& rank  & 4.52 & 22.0& 15.0 & - & - \\
   \hline\hline
$m$ & $L$ &   & $\hat{C}_{\mathsf{PACE}}$ & $\hat{C}_{\mathsf{PACE,FVE}}^{+}$ & $\hat{C}_{\mathsf{PACE,AIC}}^{+}$ & $\hat{C}_{\mathsf{PACE,BIC}}^{+}$ & $\hat{C}^{+}_{\mathsf{SC}}$\\
  \hline
5 & 2 & \textsc{aise}  &  8.99 (0.1940) &  8.29 (0.1976) &  8.29 (0.1976) &  8.29 (0.1976) &  5.20 (0.1715) [90.0\%] \\
  &   & rank & - & 6.57 & 5.58 & 5.57 & 4.15 [90.0\%] \\
 & 4 & \textsc{aise} & 11.54 (0.2016) & 10.71 (0.2030) & 10.71 (0.2030) & 10.71 (0.2030) &  7.98 (0.3392) [77.7\%] \\
  &   & rank & - & 7.16 & 6.16 & 6.16 & 5.21 [77.7\%] \\\hline
  10 & 2 & \textsc{aise} &  6.69 (0.1295) &  6.37 (0.1264) &  6.37 (0.1264) &  6.37 (0.1264) &  3.18 (0.0974) \\
    &   & rank & - & 5.92 & 4.95 & 4.93 & 4.14 \\
 & 4 & \textsc{aise} & 9.01 (0.1393) &  8.60 (0.1360) &  8.60 (0.1360) &  8.60 (0.1360) &  4.65 (0.1035) \\
    &   & rank & - & 6.46 & 5.50 & 5.49 & 5.35 \\\hline
  20 & 2 & \textsc{aise} &  5.81 (0.1159) &  5.56 (0.1118) &  5.56 (0.1118) &  5.56 (0.1118) &  2.23 (0.0804) \\
    &   & rank & - & 5.35 & 4.40 & 4.38 & 3.86 \\
 & 4 & \textsc{aise} & 7.88 (0.1205) &  7.58 (0.1160) &  7.58 (0.1160) &  7.58 (0.1160) &  3.14 (0.0812) \\
    &   & rank & - & 5.99 & 5.08 & 5.04 & 5.09 \\
   \hline\hline
\end{tabular}}

\end{table}

\renewcommand{\baselinestretch}{1.81}

\section{Real data application}\label{sec:data}

We apply the proposed method to a loop sensor dataset which
contains vehicle counts recorded every five minutes at an on-ramp on the 101
North freeway in Los Angeles, U.S.A.. This on-ramp is located near Dodger Stadium, the home field of the Los Angeles Dodgers
baseball team,
so unusual traffic is expected after a Dodgers home game.
These measurements were collected by the Freeway Performance Measurement System
(PeMS)\footnote{\url{http://pems.dot.ca.gov}}
and can be obtained from the UCI Machine Learning Repository\footnote{\url{https://archive.ics.uci.edu/ml/datasets/Dodgers+Loop+Sensor}}.
We focus on the after-game traffic measurements of 78 games between April 2005 and October 2005 available in this dataset.
For each game, we have 31 measurements that cover the time interval from 30 minutes
before the
end of the game, to 120 minutes after the end of the game.
This time interval is presented as $[-30, 120]$, where zero
marks the end of a game.

\renewcommand{\baselinestretch}{1}
\begin{figure}[ht]
  \centerline{
    \includegraphics[height=7cm]{data2.pdf}
  }
  \caption{Vehicle counts over a time interval from 30 minutes before the end
    of a game, to 120 minutes after the end of the game.
    The black line represents a  smoothing spline estimate of the mean function.}\label{fig:real:data}
\end{figure}
\renewcommand{\baselinestretch}{1.81}

The vehicle counts of the 78 games are displayed
in Figure \ref{fig:real:data}, where the mean function was estimated by smoothing splines with its tuning parameter determined by \textsc{gcv}.
The estimated mean curve demonstrates a traffic peak that
emerges at around 20 minutes after the end of a game. This characteristic is consistent with the finding of \citet{Zhang-Wang15} and conforms to common sense.

\renewcommand{\baselinestretch}{1}
\begin{figure}[ht]
  \centerline{
    \includegraphics[height=6cm]{Cfit2.pdf}
    \includegraphics[height=6cm]{Corfit2.pdf}
  }
  \caption{Contour plots of the estimated covariance (left) and correlation (right) functions respectively for $\hat{C}^{+}_\mathsf{trace}$.}\label{fig:real:contour}
\end{figure}
\renewcommand{\baselinestretch}{1.81}

We provided the covariance estimator $\hat{C}^{+}_\mathsf{trace}$ of the vehicle counts, as described in Section \ref{sec:sim},
and constructed the corresponding correlation function estimate by the simple transformation: $\hat{C}(s,t)/\{\hat{C}(s,s) \hat{C}(t,t)\}^{1/2}$ for any covariance estimate $\hat{C}$ with $\hat{C}(s,s)>0$ for all $s$.
Note that positive semi-definiteness guarantees the validity of
the correlation function estimate obtained by the above simple transformation.
Namely, it has value in-between -1 and 1.
However, this property could be violated for non-positive semi-definite
estimators such as $\hat{C}_\mathsf{CY}$.
The covariance and correlation estimate for $\hat{C}^{+}_\mathsf{trace}$ are depicted in Figure \ref{fig:real:contour}.
One intriguing feature with respect to the temporal dependency of the vehicle counts
is the high correlations of traffic between time $0$ and
time points after around time $30$.
When compared with adjacent time points such as $-20$ and $20$,
this feature is so distinctive that a ridge is formed at time $0$.

To provide further insights of such phenomenon, we investigate the $L^2$ eigen-decomposition
of $\hat{C}^{+}_\mathsf{trace}$.
Due to the built-in low-rank estimation, $\hat{C}^{+}_\mathsf{trace}$ is automatically of rank $5$ without further truncation of eigenvalues. Its corresponding five $L^2$ eigenfunctions, as described in Section \ref{subsec:eigfun}, are shown in Figure \ref{fig:real:eig_score} (Left).
The first eigenfunction explains over 80\% of the total variance, i.e., the
first eigenvalue is greater than 80\% of the sum of all five eigenvalues.
Therefore, the first eigenfunction plays a major role in the variation
of the traffic profile.
Of interest is that this eigenfunction possesses two peaks located near times $0$ and $50$, where
the second peak is spanning over the time interval roughly between 30 and 120.
This eigenfunction characterizes the high correlation we have observed between
time $0$ and the time interval between $30$ and $120$.
Since a positive variation along this eigenfunction will add traffic to
these two peaks, this implies that
some audiences may choose to leave shortly after the game or even earlier, while
some others take longer than usual to leave.
As suggested by \citet{Zhang-Wang15}, one possible explanation for this phenomenon is high game attendance.
For games with high attendance, one may choose
to leave earlier than usual to avoid traffic. Meanwhile,
heavy traffic would also last longer due to high attendance.
To further verify this explanation,
we produced the functional principal component (FPC) scores
by pre-smoothing individual vehicle count curves and then projecting them
onto the first eigenfunction.
Smoothing spline with \textsc{gcv} was used to implement the pre-smoothing.
The scatter plot between FPC scores and game attendance as shown in
Figure \ref{fig:real:eig_score} (Right), together with the fact that their Pearson correlation is $0.57$,
evidently indicates a positive association.

\renewcommand{\baselinestretch}{1}
\begin{figure}[ht]
  \centerline{
    \includegraphics[height=6cm]{eigenfun2.pdf}
    \includegraphics[height=6cm]{scores.pdf}
  }
  \caption{Left: $L^2$ eigenfunctions of $\hat{C}^{+}_\mathsf{trace}$. Right: Scatterplot of game attendance versus functional principal component scores (with respect to the first eigenfunction).}\label{fig:real:eig_score}
\end{figure}
\renewcommand{\baselinestretch}{1.81}

\section*{Acknowledgements}
The research of Raymond K.~W.~Wong is partially supported by National Science Foundation grant DMS-1612985. The research of Xiaoke Zhang is partially supported by National Science Foundation grant DMS-1613018.

\bibliographystyle{agsm}