# Piecewise-constant Neural ODEs

Sam Greydanus  Stefan Lee  Alan Fern

###### Abstract

Neural networks are a popular tool for modeling sequential data but they
generally do not treat time as a continuous variable. Neural ODEs represent an
important exception: they parameterize the time derivative of a hidden state
with a neural network and then integrate over arbitrary amounts of time. But
these parameterizations, which have arbitrary curvature, can be hard to
integrate and thus train and evaluate. In this paper, we propose making a
piecewise-constant approximation to Neural ODEs to mitigate these issues. Our
model can be integrated exactly via Euler integration and can generate
autoregressive samples in 3-20 times fewer steps than comparable RNN and ODE-
RNN models. We evaluate our model on several synthetic physics tasks and a
planning task inspired by the game of billiards. We find that it matches the
performance of baseline approaches while requiring less time to train and
evaluate.

RNNs, sequence models, neural networks, neural ODEs, autoregressive models

##  1 Introduction

Figure 1: Predicting the dynamics of two billiards balls using an RNN (left),
an ODE-RNN (center) and a PC-ODE (right). The RNN produces a hidden state
$h_{t}$ at each time step and the ODE-RNN allows it to change continuously but
is hard to integrate. Our model strikes a balance between the two by using a
continuous-time hidden state that changes linearly over each interval $\Delta
t_{i}$ and is thus easy to integrate. This allows it to skip over long spans
of predictable motion and focus on key events such as collisions.

It is said that change happens slowly and then all at once. Billiards balls
move across a table before colliding and changing trajectories; water
molecules cool slowly and then undergo a rapid phase transition into ice; and
economic systems enjoy periods of stability interspersed with abrupt market
downturns. That is to say, many time series exhibit periods of relatively
homogeneous change divided by important events. Despite this, recurrent neural
networks (RNNs), popular for time series modeling, treat time in uniform
intervals – potentially wasting prediction resources on long intervals of
relatively constant change.

A recent family of models called Neural ODEs (Chen et al., 2018) has attracted
interest as a means of mitigating these problems. They parameterize the time
derivative of a hidden state with a neural network and then integrate it over
arbitrary amounts of time. This allows them to treat time as a continuous
variable. Integration can even be performed using adaptive integrators like
Runge-Kutta (Runge, 1895), thus allocating more compute to difficult state
transitions.

Adaptive integration is especially attractive in scenarios where “key events”
are separated by variable amounts of time. In the game of billiards, these key
events may consist of collisions between balls, walls, and pockets. Between
these events, the balls simply undergo linear motion. That motion is not
difficult to predict, but it is non-trivial for a model to learn to skip over
it so as to focus on the more chaotic dynamics of collisions; this requires a
model to employ some notion of temporal abstraction (Sutton et al., 1999;
Neitz et al., 2018). This problem is not unique to billiards. The same
challenge occurs in robotics, where a robot arm occasionally interacts with
external objects at varying intervals. It may also occur in financial markets,
scientific timeseries, and other environments where change happens at a
variable rate.

One of the challenges of using adaptive integration to train neural networks
is that the integrator must perform several function evaluations in order to
estimate local curvature when performing an integration step. The curvature
information determines how far the integrator can step forward in time,
subject to a constant error budget. This is a particularly serious issue in
the context of neural networks, which may have very irregular local curvatures
at initialization. A single Neural ODE training step can take up to five times
longer to evaluate than a comparable RNN architecture, making it challenging
to scale these models (Rubanova et al., 2019). The curvature problem has, in
fact, already motivated some work on regularizing the curvature of Neural ODEs
so as to train them more efficiently (Finlay et al., 2020). But even with
regularization, these models are more difficult to train than RNNs.
Furthermore, there are many tasks where regularizing curvature is
counterproductive, for example, modeling elastic collisions between two bodies
(Jia and Benson, 2019; Chen et al., 2020).

In this work, we propose a more direct approach to constraining curvature. We
assume that our model’s hidden state dynamics can be described by a sequence
of piecewise-linear functions over variable temporal extents. Correspondingly,
our Piecewise-constant ODE (PW-ODE) model defines the time derivative as a
piece-wise constant function. This formulation enables trivial integration
over arbitrary timespans. It also eliminates the need for complicated ODE
solvers, since a simple Euler update of the form $h(\Delta
t)=h_{0}+\frac{dh_{0}}{dt}\Delta t$ is exact. Even though the resulting
dynamics are piecewise-linear in latent space, they can represent highly
nonlinear observation dynamics when paired with nonlinear encoders and
decoders. These simplifications reduce training time, simplify implementation,
and allow us to scale our model to non-linear time series with complex
observations and chaotic dynamics.

We begin by describing the connection between Neural ODEs and our model. Then
we perform a series of experiments, beginning with proof-of-concept tasks and
proceeding to more difficult ones. At each stage, we compare our model to two
baseline models: an ODE-RNN, which uses a Neural ODE to represent continuous
hidden state dynamics, and a vanilla RNN, which has piecewise-constant hidden
state dynamics. All code for this project is available
online.111github.com/greydanus/piecewise_node Our core contributions are:

  * •

Introducing a piece-wise constant Neural ODE that performs well and is simple
to implement and train

  * •

Evaluating it in comparison to ODE-RNNs and RNNs on several physics tasks,
some of them chaotic and perceptually complex

  * •

Showcasing its adaptive integration capabilities on a model-based planning
task.

##  2 Piecewise-constant Neural ODEs

Neural ODEs. Neural ODEs are a class of models that parameterize the time
derivative of the hidden state using a neural network. Letting
$\dot{h}(t)=f_{\theta}(x,t)$ where $f_{\theta}$ is the neural network and $x$
is an initial observation, we can solve for the dynamics of the hidden state
by integration:

$h(t_{1})=h(t_{0})+\int_{t_{0}}^{t_{1}}\dot{h}(t)~{}dt$ (1)

This approach is rather general: we can use it to model any continuous-time
dynamics we choose. But this same generality is also a weakness, for when we
allow $f_{\theta}$ to have arbitrary, nonlinear curvature, it becomes
difficult to integrate and train. In order to mitigate these issues, we would
like to parameterize the hidden-state dynamics with a more constrained class
of functions that are easier to integrate.

RNN hidden state dynamics. In contrast to Neural ODEs, RNNS entirely avoid
integration by applying a potentially complex update function to the hidden
state at fixed time intervals. So for RNNs, hidden state dynamics are
piecewise constant over time. Letting $z(i)$ denote an encoding of an
observation at time $i$, we can write these dynamics as:

$h(t)=h_{i}=RNN(h_{i-1},z(i))\mbox{~{}~{}for~{}~{}}t\in[i,i{+}1)$ (2)

While simple to evaluate, piecewise-constant dynamics prevent RNN hidden
states from changing between cell updates. This makes selecting an appropriate
step size an important hyperparameter. Too large and fine detail will be lost.
Too small and commputational costs increase for long time spans. Regardless,
the step sizes must be of uniform duration.

Recent works including “Neural Event Functions” in Chen et al. (2020) and ODE-
RNNs in Rubanova et al. (2019) have sought to make RNNs more expressive by
replacing these piecewise-constant dynamics with Neural ODEs between time
steps. However, these models still must contend with complex and expensive
integration of arbitrary curvature.

Piecewise-constant ODEs (PC-ODEs). We follow a similar strategy to these
works, but restrict ourselves to a simple model class where integration
between time steps can be computed efficiently. We make a locally-linear
approximation of hidden-state dynamics, writing the Taylor series expansion
around $h(t_{0})$,

$\displaystyle h(t)$
$\displaystyle=h(t_{0})+\dot{h}(t_{0})(t-t_{0})+\tfrac{1}{2}\ddot{h}(t_{0})(t-t_{0})^{2}\dots$
(3) $\displaystyle\approx h(t_{0})+\dot{h}(t_{0})\Delta t$ (4)

This simple, first-order approximation assumes that changes in the hidden
state $h$ are close to linear with respect to $\Delta t$. Unlike Neural ODEs
with arbitrary curvature, integrating this linear form is trivial; a step of
Euler integration is exact.

While local-linearity may seem like a strong assumption, we note that this is
with respect to the _hidden state_ and joint training with complex, non-linear
encoders and decoders can learn mappings that are linear to achieve low loss.
This may hold even in cases where the observation space is perceptually
complex. Consider, for example, a video of a ball undergoing linear motion.
While the dynamics are simple when projected into the right coordinate system,
they vary in a nonlinear way in pixel space. Sufficiently complex encoders and
decoders may learn this transformation if constrained to linear dynamics
within the hidden state.

The error induced by this locally-linear approximation depends on the
curvature of the underlying function. Analogously to adaptive integration
methods, a good piecewise-linear approximation should produce many short
linear segments in areas of high curvature and longer ones when change in the
hidden state is more constant.

To formalize this setup, we consider a sequence of non-uniform time intervals
$\Delta t_{i}$ with timestamps $\tau_{i}$. Put another way,
$\tau_{i}=\sum_{j=0}^{i-1}\Delta t_{j}$. This lets us define our hidden state
as:

$\displaystyle h(t)=$ $\displaystyle
h(\tau_{i})+\dot{h}(\tau_{i})(t-\tau_{i})$ (5)

where the index $i$ above denotes the most recent timestep in the sequence.
This model allows $h$ to change in a continous linear fashion between
timestamps. And at each timestamp $\tau_{i}$, $h$ can undergo arbitrary
discontinuous updates (computed by an RNN cell update).

Implementing PC-ODEs as RNNs. To implement PC-ODEs, we need to predict a
sequence of locally-constant ODEs, parameterized by $h(\tau_{i})$,
$\dot{h}(\tau_{i})$, and the duration of each segment, $\Delta t_{i}$. As
$h(\tau_{i})$ represents a discontinuous change in $h$ at $\tau_{i}$, we
predict it directly using the RNN cell update.

We build on a standard RNN cell – predicting $h_{\tau_{i}}$,
$\dot{h}_{\tau_{i}}$, and $\Delta t_{i}$ at each cell update. Letting $[]$
denote concatenation and $z(\tau_{i})$ the encoding of an observation at time
$\tau_{i}$, we write

$\displaystyle h_{\tau_{i}},\dot{h}_{\tau_{i}}=$
$\displaystyle\mbox{RNNCell}\left(\left[h_{\tau_{i-1}},\dot{h}_{\tau_{i-1}}\Delta
t_{i-1}\right],z(\tau_{i})\right)$ (6) $\displaystyle\Delta t_{i}=$
$\displaystyle 1+\mbox{LeakyReLU}(Wh_{\tau_{i}}+b)$ (7)

To compute the hidden state at some time $t$, we simply use Equation 5. This
hidden state can be decoded to an output $\hat{x}(t)$ by a decoder. Note that
this system has two notions of time – continuous time denoted $t$ and a
sequence $i=0,1,...$ of RNN cell updates occurring at continuous times
$\tau_{i}$ depending on predicted durations $\Delta t_{i}$. In contrast to a
standard RNN, which operates at a fixed time intervals, our underlying RNN
model jumps forward in time dynamically. In Equation 7, we bias our model to
always stepping forward at least one unit of time.

One challenge associated with this model is that the optimal value of $\Delta
t_{i}$ changes depending on the predicted dynamics. The dynamics, in turn,
depend on all previous values of $\Delta t$. Given that these variables are in
a dynamic relationship, we cannot use standard teacher forcing – where $\Delta
t$ is always one – for training. Rather, we must use a variant of the approach
which dynamically solves for the optimal value of $\Delta t$. We will discuss
how to do this in the next section.

##  3 Training our model

Even though our model is defined over continuous time, training data for time-
series models typically consists of observation sequences of the form
$x_{0},x_{1},x_{2},\dots,x_{T}$. For the purposes of this work, we will assume
that these observations are separated by a constant $\Delta t=1$.

When training a standard RNN on a sequence modeling task, standard practice is
to directly minimize the error of producing the next output in the ground
truth sequence given the preceding ones. This loss can be written
$\mathcal{L}_{x}=\sum_{t=0}^{T}\ell(x_{t},\hat{x}(t))$, where
$\ell(\cdot,\cdot)$ is a loss between predicted and actual values. However, in
order to learn an adaptive timestepping behavior, we must also encourage our
model to take large time steps (predict large $\Delta t_{i}$) when possible.

One way to strike a balance between accuracy and computational efficiency is
to limit the amount of error associated with each step. This is analogous to
setting the error tolerance of an ODE solver, but here the size of the time
step is determined by $\Delta t_{i}$. Let us define this requirement as
$\ell(x_{t},\hat{x}(t))<\epsilon$ where $\epsilon$ is a hyperparameter that
determines the model’s error tolerance.222We recognize the importance of
minimizing the number of hyperparameters in a model. This hyperparameter is
inevitable, though, as one cannot perform adaptive integration without setting
an error tolerance. Now our goal is to predict $\Delta t_{i}$ to be as large
as possible while satisfying this requirement. In other words, we seek to find
$\Delta t_{i}^{*}$ which satisfies

$\displaystyle\max_{\Delta t_{i}\geq 1}~{}~{}\Delta
t_{i}~{}~{}~{}~{}\mbox{s.t.}~{}~{}~{}~{}~{}~{}\ell\left(x_{t},\hat{x}(t)\right)<\epsilon$
(8) $\displaystyle\forall\left\\{x_{t}\mid t\in[\tau_{i},\tau_{i}+\Delta
t_{i}]\right\\}.$

This can be solved using a simple forward line search that starts at
$\tau_{i}$ and moves forward through the samples until our model’s prediction
loss under the current dynamics exceeds $\epsilon$. If this occurs on the
first step (implying no time passes between ticks), we set $\Delta t_{i}^{*}$
to 1. In Appendix A.1 we show this can be efficiently vectorized for an entire
batch as part of a single teacher forcing pass.333This is a practical
improvement over Neural ODEs, which are not easy to minibatch and integrate
adaptively; see discussions in Chen et al. (2018) and Rubanova et al. (2019).
We compute the optimal $\Delta t_{i}^{*}$ on-the-fly during training and step
forward accordingly – effectively ignoring the predicted $\Delta t_{i}$, as is
also done in standard teacher forcing training. To update the model to predict
the correct $\Delta t_{i}^{*}$, we augment the prediction loss $L_{x}$ with a
step size loss $\mathcal{L}_{\Delta t}=\sum_{i}||\Delta t_{i}-\Delta
t_{i}^{*}||_{2}^{2}$. At inference, the predicted $\Delta t_{i}$ sets the step
size.

Interpreting the error threshold. Examining Equation 8, $\epsilon$ acts as a
trade-off between step size and approximation error. In the extreme of setting
$\epsilon=0$, the model predictions will likely never be below $\epsilon$,
$\Delta t_{i}$ will remain at 1, and our model will perform cell updates at
the same rate as a standard RNN. At the other extreme, choosing
$\epsilon=\infty$ forces our model to predict entire time series with a
constant ODE. In practice, we found that setting $\epsilon$ to the final
training loss of the baseline RNN model444We use the baseline RNN model
because it performs better than the baseline Neural ODE model, thus offering a
more challenging target for prediction accuracy. yields models that update
3x-20x less than a standard RNN while maintaining the same test error.

##  4 Experiments

Figure 2: A comparison between PC-ODE and two baseline models trained the
circles task. Columns 1 and 2 show samples from the models. The PC-ODE model
summarizes almost the entire trajectory with a single step. Column 3 shows how
step sizes increase throughout training. Column 4 shows that our model is able
to produce autoregressive samples of comparable quality to the baseline while
using twenty times fewer cell updates.

We evaluate our model on a suite of tasks including both coordinate and video
representations of a small-scale billiards simulation. Our model attains test
errors and autoregressive sample errors comparable to those of baseline models
while performing 3-20x fewer steps. We also show how to leverage our model for
model-based planning.

Implementation Details. We use the same architecture and training regime
across all experiments unless otherwise specified. We use a GRU cell with
$128$ hidden units and our encoder and decoder networks are four-layer
residual MLPs with $128$ hidden units and ReLU activations. We train our model
with the Adam optimizer for $10,000$ steps with a batch size of $256$, a
learning rate of $1\times 10^{-3}$ and a decay schedule of $\gamma=0.9$ every
$5000$ steps. We used a mean square error loss function to train both teacher
forcing and $\Delta t$ predictions, and scaled the latter by $1\times 10^{-5}$
relative to the former. For model selection, we used early stopping.

The initial transition to stepping more than one step is challenging as it
requires latent dynamics to extrapolate essentially without supervision to do
so. To bootstrap, we force a step larger than one even when error is greater
than $\epsilon$ in 1% of cases. This increases stability; for more discussion
see Appendix A.2.

Baselines. As baselines, we consider a vanilla RNN trained at a fixed time
step ($\Delta t=1$) and an ODE-RNN (Rubanova et al., 2019) which also includes
an RNN operating at a fixed time step but predicts the parameters of an ODE at
each step in order to integrate continuous hidden states between steps. The
same encoders and decoders are used for all models. In early experiments, we
tried using a Neural ODE to solve for hidden state dynamics without
intervening RNN cell updates; such a model might be called a latent Neural
ODE. In nearly all experiments, it was worse than the ODE-RNN by at least an
order of magnitude.

###  4.1 Predicting linear motion

To begin, we consider two simple settings that involve an object undergoing
constant linear motion. These settings demonstrate the limitations of standard
RNNs and ODEs and showcase our model’s ability to translate linear latent
dynamics into non-linear motion. The first of these settings is the Lines
task, which involves motion along a line.

We generate 10,000 trajectories corresponding to linear motion along one
dimension in $(x,y)$ Cartesian coordinates. Specifically,
$\\{(x_{0},y_{0}),\dots,(x_{20},y_{20})\\}$ where $x_{t}=t$ and $y_{t}=c$
where $c$ is a random constant in $(0,1)$. This corresponds to a flat,
straight line starting at a random height. Unsurprisingly, our model quickly
learns to jump to the end of the sequence (mean $\Delta t_{i}$ of 20) while
the RNN operates at a fixed rate. Although it is theoretically possible for
the ODE-RNN to summarize this trajectory in one integration step, in practice
it required more function evaluations (and walltime) to train than the RNN. We
trained on $90\%$ of the data and evaluated on the remainder and all models
trivially achieved very low error ($10^{-8}$). See Appendix B.1 for samples.

###  4.2 Predicting circular motion

Our second task also involves constant motion, but here the motion is
circular. Even though this motion is non-linear in the Cartesian observation
space, it is linear in polar coordinates. This experiment examines our model’s
ability to encode non-linear observations with linear latent dynamics –
effectively learning a Cartesian-to-Polar conversion. As before, our dataset
contained $10,000$ examples, this time with $25$ time steps each. Each
sequence started from a random angle in $(0,2\pi)$ and a random radius in
$(1,2)$. We fixed the tangential velocity to a constant value across all
trajectories based on the intuition of a particle undergoing constant velocity
along 1D manifolds of various curvatures.

As in the lines task, our model learned to step through significant spans of
time (mean $\Delta t_{i}$ of 17). Figure 2 provides our primary visualization
of these results. The first two columns show representative samples from the
three models. Dots correspond to RNN updates with black dots corresponding to
priming observations. Recall that our model can make a prediction at any time
point between updates. We plot these intermediate evaluations with a dashed
line and they follow the ground truth trajectory closely. The third column
shows the average optimal step size $\Delta t_{i}^{*}$ during training. As we
can see, it rises during training as the model learns.

The fourth column presents a key summary plot of our experiment – a comparison
of mean squared error vs. efficiency in terms of step size. In this chart,
further lower-right is better, and corresponds to low error with large step
sizes. In green and red, we show the progression of error for the ODE-RNN and
RNN over training (light to dark). Likewise, we show our model’s progress in
blue. Unlike the RNN which is at a fixed $\Delta t$ or the ODE-RNN, which
performs even more function evaluations than the RNN, our approach changes
both in error and time step size. We see that our approach is significantly
jumpier while maintaining a similar error rate.

The gray curve represents the performance of the our model when it is forced
to take constant step sizes of size $\Delta t$. In this task, our model
converges to a point only slightly beyond this curve, indicating limited
utility in dynamic $\Delta t_{i}$ prediction compared to the optimal fixed
time step. This is unsurprising because the dynamics are uniform throughout
each sequence. We will see in subsequent experiments that dynamic predictions
of $\Delta t_{i}$ can be impactful in more complex settings. Finally, in
Figure 3, we vary the dimensionality of our model’s hidden state and retrain;
we find that model capacity and step size are positively correlated.

Figure 3: Mean step size increases as we increase model capacity (hidden state
size). We observed this effect across all datasets and experiments.

Figure 4: Comparing PC-ODE to baseline models trained on four billiards tasks.
Columns 1 and 2 illustrate qualitatively that our model’s predictions are at
least as good as those of the baselines. RNN samples are shown for the pixel
tasks, as they improve over the ODE-RNN. Columns 3 and 4 illustrate that our
our model’s time stepping varied adaptively, giving larger time steps on easy
tasks and shorter time steps on more difficult ones. Meanwhile, it matched or
improved over the predictions of the baseline models while being more
efficient to sample from. Forcing our model to use a constant $\Delta t$ (gray
line) gives much worse performance. This suggests that our model’s adaptive
time stepping skips over the most predictable time spans and concentrates on
the less predictable ones. Supplementary Materials contains videos of
predicted trajectories.

|  | Test MSE | Sample MSE | Mean $\Delta t_{i}$
---|---|---|---|---
Task | MSE Scale | Baselines | PC-ODE | Baselines | PC-ODE | Baselines | PC-ODE
1: Lines | $\times 10^{-8}$ | $1.77/1.77$ | $44.2$ | $17.9/\mathbf{13.9}$ | $21.9$ | $1/1$ | $\mathbf{20}$
2: Circles | $\times 10^{-4}$ | $0.03/0.03$ | $3.64$ | $5.47/\mathbf{3.98}$ | $5.25$ | $1/1$ | $\mathbf{17}$
3: Billiards1D | $\times 10^{-5}$ | $2.44/2.43$ | $2.93$ | $7.29/\mathbf{1.77}$ | $4.37$ | $1/1$ | $\mathbf{5.1}$
4: Billiards2D | $\times 10^{-5}$ | $5.13/2.19$ | $2.37$ | $2700/2360$ | $\mathbf{1570}$ | $1/1$ | $\mathbf{2.71}$
5: PixBill1D | $\times 10^{-4}$ | $1.57/1.57$ | $1.74$ | $3.87/\mathbf{2.37}$ | $2.64$ | $1/1$ | $\mathbf{4.1}$
6: PixBill2D | $\times 10^{-3}$ | $1.99/2.08$ | $3.90$ | $64.1/64.8$ | $\mathbf{37.1}$ | $1/1$ | $\mathbf{3.2}$
Table 1: Quantitative comparison between PC-ODE and baseline models, denoted
here as Baselines collectively for brevity. We present statistics from the two
baseline models in the format ODE-RNN/RNN.

###  4.3 Billiards from coordinates

To examine more complex dynamics, we consider a billiards-like simulation in
which two balls bounce around within a walled enclosure. We examine this
setting for 1- and 2-dimensional versions from both coordinate and image-based
inputs. Ball trajectories in this setting exhibit periods of essentially
linear motion punctuated by collisions. As such, we expect our PC-ODE model to
adjust its step sizes so as to have cell updates coincide with collision
events.

In each setting, we generate $10,000$ trajectories of 45 time steps each. The
environment lacks friction, collisions are perfectly elastic, and walls are
placed along 0 and 1 for each coordinate. For coordinate settings, we record
trajectories as $x,y$ coordinates for each ball (or just $y$ in the 1D case).
For image settings, we render $28\times 28$ grayscale images depicting a top-
down view of the balls.

The first two rows of Figure 4 visualize the result of 1D and 2D billiards
respectively. To visualize samples in the 1D case, we treat time as the
x-axis. For 2D, we fade from light to dark to indicate time. In samples from
both cases, we see the model learns to concentrate cell updates near
collisions. However, we do see instances where collisions are skipped entirely
and filled in with linear latent interpolation alone as in the 1D example. The
‘Step size vs. error’ plots (last column) also reflect the impact of this
adaptive timestepping. Our approach is able to improve in both accuracy and
average step size compared to using a fixed step size (darkest blue triangle
vs. gray curve). For 1D, the baseline models achieve lower error but require
roughly 6x the compute. In the 2D case, our model outperforms the baselines
while executing 4x fewer updates.

###  4.4 Billiards from pixels

Learning from pixels directly is a challenging setting. While the underlying
motion of the balls is the same as before, the way pixels change as the balls
move is extremely non-linear as individual pixels turn on or off abruptly. For
these experiments, we increase the hidden layer size to 512 units and train
for $50,000$ steps. The third and forth rows of Figure 4 show result summaries
for 1D and 2D pixel-based billiards. In the 1D setting, we find that both the
baselines and our model achieve similar error, but our model does so with 4x
fewer steps. In the 2D setting however, the baseline RNN yields a 50% higher
error rate and generates unrealistic trajectories where balls bend along their
trajectory like the sample shown. The ODE-RNN is even worse. Reviewing
training curves, we noted that both baseline models overfit more severely.

Quantitative results. All experimental settings are summarized with dataset-
level metrics at convergence in Table 1. We find our PC-ODE model performs at
or near the same sample error rate while making 3 to 20 times fewer RNN steps.
Note that absolute scale of error terms varies by orders of magnitude across
tasks and is specified only in the MSE Scale column for compactness.

###  4.5 Model-based Planning

Our final experiment demonstrates our PW-ODE model in a simple model-based
planning application in the 2D billiards domain. We consider a simplified
billiards game consisting of a $1\times 1$ table with a single corner pocket
at (1,0) and two balls – a cue ball and a target ball. Both balls start at
rest. The planning task is to set initial $x,y$ velocities of the cue ball
(i.e. take a shot) such that the target ball enters the pocket (< 0.17 units
from (1,0)) within 35 time steps. As our goal is to examine the usefulness of
our learned dynamics models (and not to make an advance in planning), we use a
simple random search procedure.

We consider 250 random initial configurations of the balls. For each, we
generate random initial $x,y$ velocities for the cue ball as candidate
actions. For each candidate, we simulate the resulting trajectory using each
of our learned models. We select the action that results in the minimum
distance between the target ball and the pocket over the course of the
rollout. We evalate this action by executing it in the ground-truth simulator.

Figure 5 shows the success rates of our models as a function of the number of
planning rollouts they perform. We compare our model with the baselines, an
upper bound using the actual simulator to rank actions, and a random action
lower bound. We find our PC-ODE model outperforms the other learned methods –
achieving a 61% success rate. The RNN performs nearly as well but executes
almost four times more cell updates than our model. The ODE-RNN performs
significantly worse because of its poor predictive accuracy. The wallclock
time of the ODE-RNN was also far worse, 3 times worse than the others,
especially when evaluated without a GPU.

Figure 5: We compare PC-ODE to our two baselines (RNN and ODE-RNN) and the
ground-truth simulator on a billiards planning task. We find our model is
comparable to an RNN-based planner and significantly better than the ODE-RNN
planner. Further, our model requires 3x fewer function evaluations than the
RNN and >6x fewer than ODE-RNN.

##  5 Related Work

Neural ODEs. The concept of neural ODEs has existed since the 1990’s (LeCun et
al., 1988; Pearlmutter, 1995), but has recently seen a resurgence of interest
due to (Chen et al., 2018). As already discussed, while these models and their
piecewise variants Rubanova et al. (2019); Chen et al. (2020) have shown
promise, computational efficiency is still a challenge. Our work addresses
this by imposing a restricted model that allows for trivially efficient
integration. An important point of future work is to better understand if and
when the additional computational overhead of general neural ODEs pays off in
practice, especially when coupled with powerful encoders and decoders.

RNNs for Irregular Timeseries. RNN models have been studied for handling time
series with irregular or missing time steps, e.g. including the time gap
between observations as an input to the RNN cell. In practice, though,
preprocessing the data by interpolating (Che et al., 2018) or masking (Lipton
et al., 2016) and then training a regular RNN tends to work better. Several
works have proposed applying an exponential decay function to the RNN hidden
state in proportion to the duration of time between observations (Che et al.,
2018; Cao et al., 2018; Rajkomar et al., 2018; Mozer et al., 2017). This
approach makes hidden state dynamics continuous in time in a limited sense:
hidden states are constrained to decay towards zero at a particular rate. In
contrast, our approach for adaptively predicting piecewise linear dynamics of
arbitrary length offers more flexibility.

Related RL and Planning Models. Recent work by Gregor et al. (2018) in the
context of reinforcement learning develops a jumpy planning model which does
not use an RNN cell or perform continuous interpolation of latent states.
Another relevant work is Embed to Control by Watter et al. (2015) which, like
our model, assumes that dynamics is linear in a latent space. As with our
work, they train their model on several dynamics problems, some of them from
pixels. Unlike our work, their model performs inference over uniform time
steps and does not learn an adaptive timestepping behavior.

Related work such as Adaptive Skip Intervals (Neitz et al., 2018) and Time-
Agnostic Prediction (Jayaraman et al., 2019) use a “minimum-over-time” loss
across a horizon of future states. This allows them to predict the next
predictable state, regardless of how distant in time it might be to support
temporally-abstract planning. However, they cannot reconstruct intervening
states, estimate the “rate of change” of the system, or treat time as a
continuous quantity. Infobot (Goyal et al., 2019a) includes a similar notion
of temporal abstraction but parts of Infobot still tick at every time step and
the model treats time as a discrete variable.

Temporal Abstraction in RNNs. Koutnik et al. (2014) proposed dividing an RNN
internal state into groups and only performing cell updates on the $i^{th}$
group after $2^{i-1}$ time steps. More recent works have aimed to make this
hierarchical structure more adaptive, either by data-specific rules (Ling et
al., 2015) or by a learning mechanism (Chung et al., 2019). More recently,
(Goyal et al., 2019b) have proposed Recurrent Independent Mechanisms, which
can learn to make dynamics predictions semi-independently at different time
scales. All these models exhibit temporal abstraction, but parts of them – the
encoders at the very least – must update at each time step. Additionally, they
do not treat time as a continuous variable.

##  6 Conclusions

Neural networks are already a widely used tool, but they still have
fundamental limitations. In this paper, we reckoned with the fact that they
struggle at adaptive timestepping and the computational expense of
integration. In order to make RNNs and Neural ODEs more useful in more
contexts, it is essential to find solutions to such restrictions. With this in
mind, we proposed a PC-ODE model which can skip over long durations of
comparatively homogeneous change and focus on pivotal events as the need
arises. We hope that this line of work will lead to models that can represent
time more efficiently and flexibly.

## References

  * Cao et al. [2018] W. Cao, D. Wang, J. Li, H. Zhou, L. Li, and Y. Li.  Brits: Bidirectional recurrent imputation for time series.  In _Advances in Neural Information Processing Systems_ , pages 6775–6785, 2018.
  * Che et al. [2018] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu.  Recurrent neural networks for multivariate time series with missing values.  _Scientific reports_ , 8(1):1–12, 2018.
  * Chen et al. [2018] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In _Advances in neural information processing systems_ , pages 6571–6583, 2018.
  * Chen et al. [2020] R. T. Chen, B. Amos, and M. Nickel.  Learning neural event functions for ordinary differential equations.  _arXiv preprint arXiv:2011.03902_ , 2020.
  * Chung et al. [2019] J. Chung, S. Ahn, and Y. Bengio.  Hierarchical multiscale recurrent neural networks.  In _5th International Conference on Learning Representations, ICLR 2017_ , 2019.
  * Finlay et al. [2020] C. Finlay, J.-H. Jacobsen, L. Nurbekyan, and A. M. Oberman.  How to train your neural ode: the world of jacobian and kinetic regularization.  2020\.
  * Goyal et al. [2019a] A. Goyal, R. Islam, D. Strouse, Z. Ahmed, M. Botvinick, H. Larochelle, Y. Bengio, and S. Levine.  Infobot: Transfer and exploration via the information bottleneck.  _International Conference on Learning Representations_ , 2019a.
  * Goyal et al. [2019b] A. Goyal, A. Lamb, J. Hoffmann, S. Sodhani, S. Levine, Y. Bengio, and B. Schölkopf.  Recurrent independent mechanisms.  _arXiv preprint arXiv:1909.10893_ , 2019b.
  * Gregor et al. [2018] K. Gregor, G. Papamakarios, F. Besse, L. Buesing, and T. Weber.  Temporal difference variational auto-encoder.  In _International Conference on Learning Representations_ , 2018.
  * Jayaraman et al. [2019] D. Jayaraman, F. Ebert, A. A. Efros, and S. Levine.  Time-agnostic prediction: Predicting predictable video frames.  _International Conference on Learning Representations_ , 2019.
  * Jia and Benson [2019] J. Jia and A. R. Benson.  Neural jump stochastic differential equations.  _arXiv preprint arXiv:1905.10403_ , 2019.
  * Koutnik et al. [2014] J. Koutnik, K. Greff, F. Gomez, and J. Schmidhuber.  A clockwork rnn.  In _International Conference on Machine Learning_ , pages 1863–1871, 2014.
  * LeCun et al. [1988] Y. LeCun, D. Touresky, G. Hinton, and T. Sejnowski.  A theoretical framework for back-propagation.  In _Proceedings of the 1988 connectionist models summer school_ , volume 1, pages 21–28. CMU, Pittsburgh, Pa: Morgan Kaufmann, 1988.
  * Ling et al. [2015] W. Ling, I. Trancoso, C. Dyer, and A. W. Black.  Character-based neural machine translation.  _arXiv preprint arXiv:1511.04586_ , 2015.
  * Lipton et al. [2016] Z. C. Lipton, D. Kale, and R. Wetzel.  Directly modeling missing data in sequences with rnns: Improved classification of clinical time series.  In _Machine Learning for Healthcare Conference_ , pages 253–270, 2016\.
  * Mozer et al. [2017] M. C. Mozer, D. Kazakov, and R. V. Lindsey.  Discrete event, continuous time rnns.  _arXiv preprint arXiv:1710.04110_ , 2017.
  * Neitz et al. [2018] A. Neitz, G. Parascandolo, S. Bauer, and B. Schölkopf.  Adaptive skip intervals: Temporal abstraction for recurrent dynamical models.  In _Advances in Neural Information Processing Systems_ , pages 9816–9826, 2018.
  * Pearlmutter [1995] B. A. Pearlmutter.  Gradient calculations for dynamic recurrent neural networks: A survey.  _IEEE Transactions on Neural networks_ , 6(5):1212–1228, 1995.
  * Rajkomar et al. [2018] A. Rajkomar, E. Oren, K. Chen, A. M. Dai, N. Hajaj, M. Hardt, P. J. Liu, X. Liu, J. Marcus, M. Sun, et al.  Scalable and accurate deep learning with electronic health records.  _NPJ Digital Medicine_ , 1(1):18, 2018.
  * Rubanova et al. [2019] Y. Rubanova, R. T. Chen, and D. Duvenaud.  Latent odes for irregularly-sampled time series.  _Advances in Neural Information Processing Systems_ , 2019.
  * Runge [1895] C. Runge.  Über die numerische auflösung von differentialgleichungen.  _Mathematische Annalen_ , 46(2):167–178, 1895\.
  * Sutton et al. [1999] R. S. Sutton, D. Precup, and S. Singh.  Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning.  _Artificial intelligence_ , 112(1-2):181–211, 1999.
  * Watter et al. [2015] M. Watter, J. Springenberg, J. Boedecker, and M. Riedmiller.  Embed to control: A locally linear latent dynamics model for control from raw images.  In _Advances in neural information processing systems_ , pages 2746–2754, 2015.

##  Appendix A Additional Methods

###  A.1 Efficient Batch-wise Line Search with Masking.

In Section 3 we defined the optimal jump width, $\Delta t_{i}^{*}$, and
described how to find it for a given step $i$ using the constrained
optimization objective in Equation 8. Then we described how one might solve
that equation by starting from $\tau_{i}$ and iteratively moving through the
sample sequence, decoding our model’s prediction under the current dynamics,
until the loss exceeds $\epsilon$. In this section, we describe how to
vectorize this process efficiently over a batch of examples.

At a high level, we will compute the RNN cell update batch at each time step,
but only update the parameters defining the linear region for instances that
exceed the threshold. At each point in time, the model is within a linear
segment defined by its starting point $h_{\tau_{i}}$ and its state velocity
$\dot{h}_{\tau_{i}}$. Let the $H_{t}$ be a $B\times d$ matrix of starting
points for the linear regions active at time $t$ for a batch of $B$ instances
(i.e. $h_{\tau_{i}}$ if $t\in[\tau_{i},\tau_{i}+\Delta t_{i}$). Likewise let
${M}_{i}$ be the slopes for the active linear regions. Let $\mathcal{M}$ be a
$B$ dimensional, binary vector where the $b^{\mbox{th}}$ entry denotes whether
the loss for the $b^{\mbox{th}}$ instance is greater than $\epsilon$ at the
current step. With these definitions, we can write the batched update of
$H_{t}$ and $M_{t}$ via a dynamic masking operation as follows:

$\displaystyle\hat{H}_{t},\hat{\dot{H}}_{t}$
$\displaystyle=\mbox{RNNCell}\left([H_{i},~{}\dot{H}_{i}\Delta
t_{i-1}],v\left(\hat{X}(t)\right)\right)$ $\displaystyle H_{t}$
$\displaystyle=\mathcal{M}\odot H_{t}+(1-\mathcal{M})\odot\hat{H}_{t}$
$\displaystyle\dot{H}_{t}$
$\displaystyle=\mathcal{M}\odot\dot{H}_{t}+(1-\mathcal{M})\odot\hat{\dot{H}}_{t}$
$\displaystyle\qquad\qquad\text{where}\quad\mathcal{M}=\ell\left(x_{t},\hat{x}(t)\right)<\epsilon$

And the matrix of hidden state vectors at $t$ can then be computed as:

$H(t)=H_{t}+(t-\tau_{i})\dot{H}_{t}$ (9)

Since we must compute $\mathcal{M}$ dynamically as we progress over the time
dimension, we cannot fuse the RNN cell updates over time, a common technique
used to speed up teacher forcing. This, along with the additional cost of
decoding $\hat{H}_{t}$ and evaluating prediction error, increases the runtime
of a training step of our model by a factor of three over a baseline RNN.
Unlike reference RNN implementations, however, our code is not optimized for
speed. Furthermore, we are more interested in our model’s performance during
evaluation than during training, and thus are willing to tolerate this
slowdown (which scales as a constant factor).

###  A.2 Discussion of Bootstrapping

In early experiments, we noticed that our model had a difficult time
transitioning from no jumpiness (selecting $\Delta t=1$ at every time step) to
some jumpiness (selecting $\Delta t=2$ at rare intervals). Upon further
investigation, we found that, having never needed to linearly extrapolate
dynamics for more than one time step, our model was very bad at doing so early
in training. In fact, the prediction errors at $\Delta t=2$ were consistently
much higher than $\epsilon$. We considered two methods of solving this
problem. One was to anneal $\epsilon$ starting from a large value to promote
larger step sizes early in training. The other was to make the model to take
an extra linear step 1% of the time, regardless of the value of $\epsilon$.
The second option gave the best and most stable results so we used it in our
experiments.

##  Appendix B Additional experimental results

###  B.1 Visualizing the Lines Task

Figure 6: A comparison between our model and the two baselines trained on the
Lines task. Columns 1 and 2 show samples from the models. Here we see that the
PC-ODE model can summarize an entire trajectory with a single linear function,
whereas the baseline models cannot. Columns 3 and 4 report key quantitative
statistics of the two models. In particular, Column 3 shows how step size
increases throughout training. Column 4 shows that the PC-ODE model is able to
produce autoregressive samples of comparable quality to the baselines while
using twenty times fewer cell updates.


# Tensor network simulations for non-orientable surfaces

Haruki Shimizu<EMAIL_ADDRESS>Atsushi Ueda  Institute for
Solid State Physics, University of Tokyo, Kashiwa 277-8581, Japan

###### Abstract

In this study, we explore the geometric construction of the Klein bottle and
the real projective plane ($\mathrm{RP}^{2}$) within the framework of tensor
networks, focusing on the implementation of crosscap and rainbow boundaries.
Previous investigations have applied boundary matrix product state (BMPS)
techniques to study these boundaries. We introduce a novel approach that
incorporates such boundaries into the tensor renormalization group (TRG)
methodology, facilitated by an efficient representation of a spatial
reflection operator. This advancement enables us to compute the crosscap and
rainbow free energy terms and the one-point function on $\mathrm{RP}^{2}$ with
enhanced efficiency and for larger system sizes. Additionally, our method is
capable of calculating the partition function under isotropic conditions of
space and imaginary time. The versatility of this approach is further
underscored by its applicability to constructing other (non-)orientable
surfaces of higher genus.

††preprint: APS/123-QED

##  I Introduction

Classifying the phases of matter represents a fundamental challenge and goal
within the field of physics. This classification enables the tailored design
of models to exhibit desired phenomena, navigating the intricate landscape of
physical properties. However, the complexity inherent in realistic models
presents significant challenges. To address these, the concept of universality
in critical phenomena occurring at the phase boundaries becomes instrumental.
Universal quantities, such as critical exponents, offer insights into the
mechanisms underpinning phase transitions. These insights are invaluable, as
the fundamental processes of phase transitions are often not apparent through
the microscopic details of the systems involved.

The utility of universality is underscored in the study of critical phenomena
within both one-dimensional quantum systems and two-dimensional classical
systems, where conformal field theory (CFT) serves to describe these universal
properties. Recent advancements have further highlighted the importance of
CFT-derived universal quantities in analyzing edge states within
(2+1)-dimensional symmetry-protected topological (SPT) phases, making the
acquisition of CFT data crucial for understanding the underlying mechanisms of
phase transitions [1, 2, 3].

An intriguing development in this context is that universal quantities can be
discerned through the partition function on non-orientable surfaces. This
approach indicates that a spatial reflection introduces a universal term to
the partition function. In particular, the ratio of the partition function on
the Klein bottle to that on the torus is given by a universal ratio [4, 5]:

$\displaystyle\frac{Z^{\mathcal{K}}(\frac{\beta}{2},2L)}{Z^{\mathcal{T}}(\beta,L)}\simeq
g,$ (1)

when the system length is significantly larger than its inverse temperature
$L\gg\beta$. Here, $g$ is a universal constant unique to its universality
class. More precisely, $g$ is $\sum_{a}\frac{d_{a}}{\mathcal{D}}$ for diagonal
CFTs, with $d_{a}$ being the quantum dimension of the primary field $a$, and
$\mathcal{D}$ representing the total quantum dimension. As $g$ takes a unique
value depending on the criticality, it can be used as a tool to detect phase
transitions ($g=1$ for trivially gapped phases). This theory has been
corroborated through numerical studies, including Monte-Carlo and tensor-
network simulations employing boundary matrix product state (BMPS) techniques
[6, 7, 8, 9, 10, 11, 12, 13].

This concept was further expanded to include the non-orientable surface
$\mathrm{RP^{2}}$, leading to the discovery of another universal ratio [14,
10]:

$\displaystyle\frac{Z^{\mathrm{RP^{2}}}(\frac{\beta}{2},2L)}{Z^{\mathcal{T}}(\beta,L)}\sim\beta^{c/4},$
(2)

where $c$ denotes the central charge, a critical universal parameter within
CFT. The identification of these universal ratios is pivotal for elucidating
the universal critical and topological characteristics of systems. Yet, direct
simulation of geometrically twisted boundaries, crucial for the realization of
both the Klein bottle and $\mathrm{RP}^{2}$, presents considerable challenges.
As a workaround, previous efforts have resorted to approximations, such as
computing crosscap and rainbow states via BMPS techniques, applicable
predominantly in the anisotropic limit where $L\gg\beta$. Given the
foundational role of $\mathrm{RP^{2}}$ in constructing generalized non-
orientable surfaces, its accurate realization across a broader range of cases
remains a critical objective.

Addressing this gap, we introduce tensor-network strategies for simulating the
Klein bottle and $\mathrm{RP^{2}}$ configurations directly. By integrating the
higher-order tensor renormalization group (HOTRG) [15] with spatial
reflections, we enable the computation of $Z^{\mathcal{K}}$ and
$Z^{\mathrm{RP^{2}}}$ for extensive system sizes, thereby facilitating a
deeper understanding of universal critical and topological features.

In addition to this, we demonstrate that HOTRG can be utilized to extract
another piece of universal data. The partition function of CFT on a torus is
determined by its operator content and operator product expansion (OPE)
coefficients, $C_{ijk}$. Prior research has illustrated the effectiveness of
both tensor renormalization group (TRG) and tensor network renormalization
(TNR) methods in efficiently retrieving these CFT data [16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27]. For non-orientable surfaces, however, a different
type of universal data is required—the one-point function of the primary
operators [4]:

$\displaystyle\langle\phi_{k}\rangle_{\mathrm{RP^{2}}}=\frac{\Gamma_{k}}{(1+z\bar{z})^{2h}},$
(3)

with $\mathrm{RP}^{2}$ conceptualized as the complex plane $z$, under the
involution $I(z)=\frac{-1}{\bar{z}}$ and $h$ being a conformal weight. This
conceptualization stems from viewing $\mathrm{RP}^{2}$ as a projective space
in which antipodal points on a Riemann sphere of radius $\frac{1}{2}$ are
identified. Our findings indicate that the parameter $\Gamma_{k}$, pivotal for
understanding the behavior of physical quantities on non-orientable Riemann
surfaces, can be directly calculated using both TRG and TNR approaches.

This paper is organized as follows: Initially, we provide an overview of
tensor network methodologies and the concept of CFT on the Klein bottle and
$\mathrm{RP}^{2}$ in Section II. Subsequently, The numerical outcomes derived
from our algorithm are showcased in Section III. Finally, we explore the
potential implications and applicability of our findings in Section IV.

Our numerical calculations were performed using ITensors [28]. Sample Codes to
reproduce the figures in Section III can be found at the following repository:
https://github.com/elle-et-noire/HOTRG-nonorientablesurface

##  II Method

In this study, we focus on statistical models defined on a square lattice with
nearest neighbor interactions. Although our demonstrations primarily utilize
the HOTRG framework, the methodologies we introduce are versatile enough to be
applied across a range of TRG and TNR algorithms.

###  II.1 Review on HOTRG

In statistical mechanics, the probability of a specific configuration is
captured by its Boltzmann weight. The aggregate of these weights across all
possible configurations yields the partition function, a cornerstone for
calculating various physical properties. Within the framework of tensor
network representations, this comprehensive summation process is facilitated
by tensor contractions. Contracting tensors effectively sums over the indices
being contracted, which correspond to the degrees of freedom in the system,
such as a spin in lattice models.

Let $T^{(0)}=T$ be an initial tensor of HOTRG which represents the local
Boltzmann weights of the partition function. The torus partition function is
obtained by contracting all the tensors in the system with bonds connecting
the nearest neighboring sites:

$Z=\ \
\smashoperator[]{\sum_{\quantity{x_{i},x_{i}^{\prime},y_{i},y_{i}^{\prime}}}^{}}\
\ \prod T_{x_{i}x_{i}^{\prime}y_{i}y_{i}^{\prime}}$

with periodic boundary condition.

At the $n$-th step of coarse graining, first we contract the previously
obtained tensor $T^{(n-1)}$ in vertical direction:

$M^{(n)}_{x_{1}x_{2}x_{1}^{\prime}x_{2}^{\prime}yy^{\prime}}=\sum_{i}T^{(n-1)}_{x_{1}x_{1}^{\prime}yi}T^{(n-1)}_{x_{2}x_{2}^{\prime}iy^{\prime}}.$
(4)

Then we obtain the isometry by singular value decomposition (SVD):

$\sum_{xx^{\prime}yy^{\prime}}M^{(n)}_{x_{1}x_{2}xx^{\prime}yy^{\prime}}M^{(n)}_{x_{1}^{\prime}x_{2}^{\prime}xx^{\prime}yy^{\prime}}\approx\sum_{ii^{\prime}}U^{(n)}_{ix_{1}x_{2}}\Lambda^{(n)}_{ii^{\prime}}U^{(n)}_{i^{\prime}x_{1}^{\prime}x_{2}^{\prime}}$

where $i,i^{\prime}$ is truncated up to $\chi$ dimension. Contracting the
isometry $U^{(n)}$, we get vertically renomalized tensor:

$W^{(n)}_{xx^{\prime}yy^{\prime}}=\
\smashoperator[]{\sum_{x_{1}x_{1}^{\prime}y_{1}y_{1}^{\prime}}^{}}\
U^{(n)}_{xx_{1}x_{2}}M^{(n)}_{x_{1}x_{2}x_{1}^{\prime}x_{2}^{\prime}yy^{\prime}}U^{(n)}_{x^{\prime}x_{1}^{\prime}x_{2}^{\prime}}.$

Then, we repeat the same procedure in horizontal direction:

$T^{(n)}_{xx^{\prime}yy^{\prime}}=\
\smashoperator[]{\sum_{iy_{1}y_{1}^{\prime}y_{2}y_{2}^{\prime}}^{}}\
V^{(n)}_{yy_{1}y_{2}}W^{(n)}_{xiy_{1}y_{1}^{\prime}}W^{(n)}_{ix^{\prime}y_{2}y_{2}^{\prime}}V^{(n)}_{y^{\prime}y_{1}^{\prime}y_{2}^{\prime}}$

with an appropriate isometry $V^{(n)}$ obtained by SVD. Thus, the height and
width of the system renormalized into $T^{(n)}$ is equal.

###  II.2 Renormalization of a spatial reflection operator

Figure 1: Illustrations of the partition functions on $(a)$ the Klein bottle
and $(b)$ $\mathrm{RP}^{2}$. These partition functions can be cut in the
middle as shown in the blue arrows(left panel), and then sewn again after
flipping the latter half in $x$-direction. The resulting boundaries are
crosscap states for the Klein bottle and rainbow and crosscap states for
$\mathrm{RP}^{2}$. The schematic picture of the crosscap and rainbow
boundaries are shown in $(c)$ and $(d)$, respectively. $(e)$ In the limit
$L\gg\beta\gg 1$, the bulk part of the partition function
$Z^{\mathrm{RP^{2}}}(\beta,L)$ is dominated by the leading eigenvector
$|i_{0}\rangle$ of the transfer matrix on a cylinder. The boundary
contributions are then represented as the overlaps with $|i_{0}\rangle$ as
denoted as $\langle\mathcal{R}|i_{0}\rangle$ and $\langle
i_{0}|\mathcal{C}\rangle$.

Consider a system characterized by a width $L$ along the $x$-direction and a
height $\beta$ along the $y$-direction, adhering to the thermodynamic limit
$L\gg\beta\gg 1$. To compute the partition functions for the Klein bottle and
$\mathrm{RP}^{2}$ within a tensor network framework, we employ the cut-and-sew
method as described in Refs. [4, 10]. We briefly review these concepts for
clarity.

The partition function on the Klein bottle is periodic in the spatial
direction and undergoes a geometric twist in the $\beta$-direction, as
depicted in Fig. 1$(a)$. To evaluate this quantity, we bisect the partition
function medially, then overlay the right half atop the left, subsequent to a
spatial inversion. This process is visualized in the right panel of Fig.
1$(a)$. The resultant partition function has a width of $L/2$ and a height of
$2\beta$. Notably, the $\beta$-direction retains its periodicity, whereas the
boundaries at $x=0$ and $L/2$ are no longer periodic. Instead, they adopt
crosscap boundary conditions, as illustrated in $(c)$.

A similar “cut and sew” process applied to $\mathrm{RP^{2}}$ results in the
imposition of rainbow and crosscap boundaries at $x=0$ and $L/2$,
respectively, as depicted in $(b)$. The rainbow boundary is further
illustrated in $(d)$. Notably, the $\beta$-direction maintains periodicity in
$\mathrm{RP^{2}}$ as well. This periodicity is instrumental in facilitating
the computation of partition functions on both the Klein bottle and
$\mathrm{RP^{2}}$.

Let $M_{\beta}$ be a transfer matrix in $x$-direction, which has height
$\beta$. In the limit $L\gg\beta\gg 1$, the partition function
$Z(\beta,L)=\matrixelement*{\mathcal{B}_{L}}{M_{\beta}^{L}}{\mathcal{B}_{R}}$
is approximated using the largest eigenvector and its eigenvalue, respectively
denoted as $|i_{0}\rangle_{\beta}$ and $\lambda_{0}(\beta)$:

$Z(\beta,L)\stackrel{{\scriptstyle
L\gg\beta}}{{\approx}}\innerproduct{\mathcal{B}_{L}}{i_{0}}_{\beta}[\lambda_{0}(\beta)]^{L}\leftidx{{}_{\beta}}{\innerproduct{i_{0}}{\mathcal{B}_{R}}}{}.$

In the case of the torus, the Klein bottle and $\mathrm{RP^{2}}$, the
corresponding partition functions are given as

$\displaystyle Z^{\mathcal{T}}(\beta,L)$
$\displaystyle\approx[\lambda_{0}(\beta)]^{L},$ (5) $\displaystyle
Z^{\mathcal{K}}(\beta,L)$
$\displaystyle\approx\innerproduct{\mathcal{C}}{i_{0}}_{2\beta}[\lambda_{0}(2\beta)]^{L/2}{\phantom{\lvert}}_{2\beta}\innerproduct{i_{0}}{\mathcal{C}},$
(6) $\displaystyle Z^{\mathrm{RP^{2}}}(\beta,L)$
$\displaystyle\approx\innerproduct{\mathcal{R}}{i_{0}}_{2\beta}[\lambda_{0}(2\beta)]^{L/2}{\phantom{\lvert}}_{2\beta}\innerproduct{i_{0}}{\mathcal{C}},$
(7)

where $\ket{\mathcal{C}}$ ($\ket{\mathcal{R}}$) represents a crosscap
(rainbow) state. Equation (7) is pictorially shown in Fig. 1 $(e)$.

Given (5-7), $Z^{\mathcal{K}}$ and $Z^{\mathcal{\mathrm{RP^{2}}}}$ contain
extra terms arising from the boundaries compared to $Z^{\mathcal{T}}$ as
follows:

$\displaystyle\ln\dfrac{Z^{\mathcal{K}}(\frac{\beta}{2},2L)}{Z^{\mathcal{T}}(\beta,L)}$
$\displaystyle=2\ln\absolutevalue{\innerproduct{\mathcal{C}}{i_{0}}},$
$\displaystyle\ln\dfrac{Z^{\mathrm{RP^{2}}}(\frac{\beta}{2},2L)}{Z^{\mathcal{T}}(\beta,L)}$
$\displaystyle=\ln\absolutevalue{\innerproduct{\mathcal{C}}{i_{0}}}+\ln\absolutevalue*{\innerproduct{\mathcal{R}}{i_{0}}_{\beta}}$

This boundary free energy
$F_{\mathcal{C}}\coloneqq\ln\absolutevalue{\innerproduct{\mathcal{C}}{i_{0}}}$
and
$F_{\mathcal{R}}\coloneqq\ln\absolutevalue*{\innerproduct{\mathcal{R}}{i_{0}}_{\beta}}$
is respectively referred to as crosscap and rainbow free energy. In CFT,
$F_{\mathcal{C}}$ and $F_{\mathcal{R}}$ are expressed with universal terms as
[6, 10]

$\displaystyle F_{\mathcal{C}}$ $\displaystyle=\frac{1}{2}\ln g,$ (8)
$\displaystyle F_{\mathcal{R}}$ $\displaystyle=\frac{c}{4}\ln\beta+b,$ (9)

where $c$ and $b$ are the central charge and non-universal constant. Thus,
computing the crosscap and rainbow free energy allows us to directly read out
the universal data.

To compute $F_{\mathcal{C}}$ or $F_{\mathcal{R}}$ via tensor network
techniques, consider the scenario at the $n$-th step of HOTRG, specifically
when $\beta=2^{n}$. Here, the matrix $M_{\beta}$ is essentially derived from
$T^{(n)}$ by contracting its indices along the $y$-direction. The leading
eigenvector of $M_{\beta}$, denoted as $\ket{i_{0}}_{\beta}$, is required to
possess two tensor indices to facilitate contraction with the states
$\ket{\mathcal{C}}$ or $\ket{\mathcal{R}}$, each of which also manifests in a
dual-index tensor format as elucidated subsequently. This necessitates the
vertical duplication of $T^{(n-1)}$, leading to the formulation:

$(M_{\beta})_{(x_{1}x_{2})(x_{1}^{\prime}x_{2}^{\prime})}=\sum_{y}M^{(n)}_{x_{1}x_{2}x_{1}^{\prime}x_{2}^{\prime}yy}\
{},$

where $(M_{\beta})_{(x_{1}x_{2})(x_{1}^{\prime}x_{2}^{\prime})}$ represents
the transfer matrix along the $x$-direction. In this context, $(x_{1}x_{2})$
and $(x_{1}^{\prime}x_{2}^{\prime})$ serve as the composite matrix indices,
with the notation $(ab)$ signifying the amalgamation of indices $a$ and $b$.

Let ${L}_{x_{1}x_{2}}^{(n)}$ be the leading eigenvector of $M_{\beta}$:

$\sum_{x_{1}x_{2}}L_{x_{1}x_{2}}^{(n)}(M_{\beta})_{(x_{1}x_{2})(x_{1}^{\prime}x_{2}^{\prime})}=\lambda_{0}(\beta)L_{x_{1}^{\prime}x_{2}^{\prime}}^{(n)}.$

Then we can assume $L^{(n)}_{x_{1}x_{2}}$ as a tensor representation of
$\ket{i_{0}}_{\beta}$ with two indices.

$\approx$$\approx$

Figure 2: Tensor network representations of a crosscap boundary condition. The
six triangles are the left-hand side isometries to renormalize the bulk tensor
$T^{(0)}$ into $T^{(2)}$, with four smaller ones being $U^{(1)}$ and two
larger ones $U^{(2)}$. The arrows indicate the order of two indices before
amalgamation and truncation, which correspond to the order of $x_{1}$ and
$x_{2}$ in $U^{(n)}_{xx_{1}x_{2}}$.

####  II.2.1 Crosscap free energy term

To construct a crosscap within the coarse-grained tensor framework, we proceed
by contracting the indices $x_{1}$ and $x_{2}$ of $L^{(n)}_{x_{1}x_{2}}$, as
illustrated in Fig. 2. This operation effectively encapsulates the crosscap
structure at the $n$-th level of coarse-graining. The crosscap free energy,
${F}_{\mathcal{C}}^{(n)}$, can subsequently be computed as follows:

${F}_{\mathcal{C}}^{(n)}=\ln\quantity\Bigg(\sum_{x}{L}^{(n)}_{xx}).$ (10)

This equation provides a direct method for evaluating the free energy
associated with the crosscap configuration, leveraging the simplicity of
tensor contractions to elucidate the energetic contributions of topologically
non-trivial structures within the tensor network formalism.

####  II.2.2 Rainbow free energy term

$O$$\approx$$O$$O$$\approx$

Figure 3: Tensor network representations of a rainbow boundary condition. With
the spatial reflection operator $O^{(n)}$ in between the isometries, the lower
half of the indices are in the reverse order.

$\approx$$\approx$

Figure 4: Alternative way to construct a rainbow boundary condition in tensor
network. Instead of contracting a renormalized spatial reflection operator, we
just reflect the vertically copied bulk tensor $T^{(n)}$ in vertical
direction, which means that the included isometries are also reflected.

In order to construct a rainbow state, we renormalize a spatial reflection
operator which inverts the direction of a edge of the bulk tensor $T^{(n)}$.

Let the initial spatial reflection operator be an unit matrix:

$O^{(0)}_{ab}=\delta_{ab}.$

At the $n$-th step, we renormalize the spatial reflection operator reusing the
isometry used to renormalize $T^{(n)}$:

$O^{(n)}_{ab}=\sum_{ii^{\prime}jj^{\prime}}U^{(n)}_{aij}O^{(n-1)}_{ii^{\prime}}O^{(n-1)}_{jj^{\prime}}U^{(n)}_{bj^{\prime}i^{\prime}}.$

This can be graphically represented as below:

$\scalebox{0.6}{ \leavevmode\hbox to233.21pt{\vbox
to59.37pt{\pgfpicture\makeatletter\hbox{\hskip
61.68889pt\lower-29.68738pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{1.6pt}\pgfsys@invoke{ }
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {{}}{}{{}}{}{{}}
{}{}{{}} {}{} {}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{.75,1,0}\pgfsys@color@rgb@fill{.75}{1}{0}\pgfsys@invoke{
}{}\pgfsys@moveto{-28.45276pt}{28.45276pt}\pgfsys@lineto{-28.45276pt}{-28.45276pt}\pgfsys@lineto{-46.6625pt}{0.0pt}\pgfsys@closepath\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{}}{}{{}}{}{{}} {}{}{{}}
{}{} {}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{.75,1,0}\pgfsys@color@rgb@fill{.75}{1}{0}\pgfsys@invoke{
}{}\pgfsys@moveto{28.45276pt}{28.45276pt}\pgfsys@lineto{28.45276pt}{-28.45276pt}\pgfsys@lineto{46.6625pt}{0.0pt}\pgfsys@closepath\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{}}{}{}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}{}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{}\pgfsys@moveto{18.4568pt}{18.96648pt}\pgfsys@curveto{18.4568pt}{24.4457pt}{14.01512pt}{28.88737pt}{8.5359pt}{28.88737pt}\pgfsys@curveto{3.05669pt}{28.88737pt}{-1.385pt}{24.4457pt}{-1.385pt}{18.96648pt}\pgfsys@curveto{-1.385pt}{13.48726pt}{3.05669pt}{9.04558pt}{8.5359pt}{9.04558pt}\pgfsys@curveto{14.01512pt}{9.04558pt}{18.4568pt}{13.48726pt}{18.4568pt}{18.96648pt}\pgfsys@closepath\pgfsys@moveto{8.5359pt}{18.96648pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{4.58313pt}{15.54982pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$O$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{}\pgfsys@moveto{18.4568pt}{-18.96648pt}\pgfsys@curveto{18.4568pt}{-13.48726pt}{14.01512pt}{-9.04558pt}{8.5359pt}{-9.04558pt}\pgfsys@curveto{3.05669pt}{-9.04558pt}{-1.385pt}{-13.48726pt}{-1.385pt}{-18.96648pt}\pgfsys@curveto{-1.385pt}{-24.4457pt}{3.05669pt}{-28.88737pt}{8.5359pt}{-28.88737pt}\pgfsys@curveto{14.01512pt}{-28.88737pt}{18.4568pt}{-24.4457pt}{18.4568pt}{-18.96648pt}\pgfsys@closepath\pgfsys@moveto{8.5359pt}{-18.96648pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{4.58313pt}{-22.38313pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$O$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {{}}{}{{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{}{} {{}{}{{}}{}}{{}{}}{{}{}}{{}}
{{}{}{{}}{}}{{{}}{{}}}{{}}{{}{}}{{{}}{{}}}{{}}{}{{}}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}\pgfsys@moveto{28.45276pt}{18.96648pt}\pgfsys@lineto{19.2568pt}{18.96648pt}\pgfsys@moveto{-2.02222pt}{17.10486pt}\pgfsys@curveto{-19.22157pt}{14.07224pt}{-14.14665pt}{-8.94922pt}{-28.45276pt}{-18.96648pt}\pgfsys@stroke\pgfsys@invoke{
} {{}}{}{{}}{}{{}} {{{{{}}{}{}{}{}{{}}}}}{}{}{} {{}{}{{}}{}}{{}{}}{{}{}}{{}}
{{}{}{{}}{}}{{{}}{{}}}{{}}{{}{}}{{{}}{{}}}{{}}{}{{}}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}\pgfsys@moveto{28.45276pt}{-18.96648pt}\pgfsys@lineto{19.2568pt}{-18.96648pt}\pgfsys@moveto{-2.02222pt}{-17.10486pt}\pgfsys@curveto{-19.22157pt}{-14.07224pt}{-14.14665pt}{8.94922pt}{-28.45276pt}{18.96648pt}\pgfsys@stroke\pgfsys@invoke{
} {{}}{}{{}}{}{{}}
{}{}{}\pgfsys@moveto{-28.45276pt}{28.45276pt}\pgfsys@lineto{-28.45276pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
}{\pgfsys@beginscope\pgfsys@invoke{ }
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}}
}}{{}}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}\pgfsys@beginscope\pgfsys@invoke{
} {\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0} {{}}
{{}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{-1.0}{1.0}{0.0}{-28.45276pt}{5.55005pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}
{{{}}} }{{}{}}{{}{}}{{{{}{}{{}} }}{{}} {{{}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {{}}{}{{}}{}{{}}
{}{}{}\pgfsys@moveto{28.45276pt}{28.45276pt}\pgfsys@lineto{28.45276pt}{-28.45276pt}\pgfsys@stroke\pgfsys@invoke{
}{\pgfsys@beginscope\pgfsys@invoke{ }
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}}
}}{{}}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}\pgfsys@beginscope\pgfsys@invoke{
} {\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0} {{}}
{{}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{-1.0}{1.0}{0.0}{28.45276pt}{5.55005pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}}{{{{}{}{{}} }}{{}}
{{{}}} }{{}{}}{{}{}}{{{{}{}{{}} }}{{}} {{{}}} }
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} {{}}{}{{}}{}
{}{}{}\pgfsys@moveto{46.6625pt}{0.0pt}\pgfsys@lineto{60.88889pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} {{}}{}{{}}{}
{}{}{}\pgfsys@moveto{-46.6625pt}{0.0pt}\pgfsys@lineto{-60.88889pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{} {}{}{}{{ {\pgfsys@beginscope{} {} {} {}
\pgfsys@moveto{3.79997pt}{0.0pt}\pgfsys@lineto{-2.27998pt}{3.03998pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@lineto{-2.27998pt}{-3.03998pt}\pgfsys@fill\pgfsys@endscope}}
}{}{}{{}}\pgfsys@moveto{75.11526pt}{0.0pt}\pgfsys@lineto{95.78468pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{95.78468pt}{0.0pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{1,0.625,0.25}\pgfsys@color@rgb@fill{1}{0.625}{0.25}\pgfsys@invoke{
}{}\pgfsys@moveto{157.2638pt}{0.0pt}\pgfsys@curveto{157.2638pt}{8.28438pt}{150.54817pt}{15.0pt}{142.2638pt}{15.0pt}\pgfsys@curveto{133.97942pt}{15.0pt}{127.2638pt}{8.28438pt}{127.2638pt}{0.0pt}\pgfsys@curveto{127.2638pt}{-8.28438pt}{133.97942pt}{-15.0pt}{142.2638pt}{-15.0pt}\pgfsys@curveto{150.54817pt}{-15.0pt}{157.2638pt}{-8.28438pt}{157.2638pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{142.2638pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{138.31102pt}{-3.41666pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$O$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{}{}
{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}\pgfsys@moveto{113.81104pt}{0.0pt}\pgfsys@lineto{126.46379pt}{0.0pt}\pgfsys@moveto{158.0638pt}{0.0pt}\pgfsys@lineto{170.71655pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}} }.$

Using the spatial reflection operator and the leading eigenvector
$L^{(n)}_{x_{1}x_{2}}$, we can construct a rainbow boundary in coarse grained
tensor and obtain the $n$-th rainbow term, $F_{\mathcal{R}}^{(n)}$, as follows
111Note that the order $n$ of $O^{n}$ is crucial. The indices $x_{1}$ and
$x_{2}$ of $L^{(n)}_{x_{1}x_{2}}$ are coming from the left-hand side index of
$T^{(n-1)}$, which is originally the truncated index of $U^{(n-1)}$. Because
$O^{(n-1)}$ also has indices of $U^{(n-1)}$, $L^{(n)}$ has to be contracted
with $O^{(n-1)}$.:

${F}_{\mathcal{R}}^{(n)}=\ln\quantity\Bigg(\,\sum_{x_{1}x_{2}}{L}^{(n)}_{x_{1}x_{2}}O^{(n-1)}_{x_{1}x_{2}}).$
(11)

This process is depicted in Fig. 3.

An alternative approach to constructing a rainbow state involves inverting the
direction of the copied tensor at the vertical contraction step in (4), which
effectively reflects the lower half of the left vertical boundary without
necessitating the renormalization of the spatial reflection operator. This is
expressed as:

$\widetilde{M}^{(n)}_{x_{1}x_{2}x_{1}^{\prime}x_{2}^{\prime}yy^{\prime}}=\sum_{i}{T}^{(n-1)}_{x_{1}x_{1}^{\prime}yi}{T}^{(n-1)}_{x_{2}x_{2}^{\prime}y^{\prime}i}.$
(12)

Following this adjustment, the subsequent HOTRG steps proceed by treating
$\widetilde{M}^{(n)}$ as if it were $M^{(n)}$. This method provides a
straightforward mechanism for embodying the geometric nuances of a rainbow
state within the tensor network formalism, leveraging the inherent symmetries
and structural modifications to simulate complex boundary conditions
effectively.

In this construction, we get rainbow boundary state by just contracting the
indices of the largest eigenvector as shown in Fig. 4:

$\widetilde{F}_{\mathcal{R}}^{(n)}=\ln\quantity\Bigg(\sum_{x}\widetilde{L}_{xx}^{(n)}),$
(13)

where $\widetilde{L}^{(n)}$ is the largest eigenvector obtained from
diagonalizing $\widetilde{M}^{(n)}$ in (12) 222One might expect that the
crosscap free energy can be calculated also in this case using the spatial
reflection operator, but it does not work. The reason would be that the
isometry ${\widetilde{U}}^{(n)}_{xx_{1}x_{2}}$ becomes highly symmetric in
$y$-direction (to the exchange between $x_{1}\leftrightarrow x_{2}$) and the
information of reflection is tend to drop by truncation..

##  III Results

41664$0.25$$0.50$$0.75$$1.00$$\beta$$F_{\mathcal{R}}$Ising ($c=0.5$)Potts
($c=0.8$)Fitting Formula:$F_{\mathcal{R}}=\frac{c}{4}\ln\beta+b$

Figure 5: The rainbow free energy term $F_{\mathcal{R}}$ of the Ising
(“$\times$”) and three-state Potts model (“$+$”) obtained from HOTRG with the
bond dimension $\chi=45$. For each model, these values are obtained by both
the renormalization of the spatial reflection operator as in (11), and by
reflecting the left-hand side isometries as in (13), which are consistent and
indistinguishable in the figure. They are compared with the theoretical
predictions denoted by the orange and purple dotted lines.

41664$0.20$$0.25$$0.30$$0.35$$0.40$$0.45$$\beta$$F_{\mathcal{C}}$$\frac{1}{2}\ln
g_{\mathrm{Ising}}=\frac{1}{2}{\ln}\hskip
1.20007pt(1+1/\sqrt{2})$$\frac{1}{2}\ln
g_{\mathrm{Potts}}=\frac{1}{2}\ln\sqrt{3+6/\sqrt{5}}$

Figure 6: The crosscap free energy term $F_{\mathcal{C}}$ of the Ising
(“$\times$”) and three-state Potts model (“$+$”) obtained from HOTRG with the
bond dimension $\chi=45$. For each model, these values are obtained by the
renormalization of the spatial reflection operator as in (10).

We implemented the crosscap and rainbow boundary conditions on both the Ising
and three-state Potts model to calculate the free energy terms
$F_{\mathcal{C}}$ and $F_{\mathcal{R}}$. As shown in Fig. 5 and Fig. 6, both
terms obtained from tensor network calculations fit well with (1) or (2), up
to $\beta\sim 128$ and $\beta\sim 32$ for the Ising model and three-state
Potts model, respectively. The rainbow free energy term is also calculated by
(13), with resulting $\widetilde{F}_{\mathcal{R}}$ in perfect agreement with
$F_{\mathcal{R}}$, the value obtained by using spatial reflection operator.
The deviation from the theoretical value for the three-state Potts model is
attributed to a finite-bond dimension effect. Notably, an enhancement in
accuracy was achieved by increasing the bond dimension $\chi$ (Not shown
here). These findings underscore the effectiveness and precision of our
proposed methodologies. Finally, we proceed to calculate another set of CFT
data, specifically $\Gamma_{k}$, as outlined in Eq. (3). This equation can
also be interpreted as the overlap $\langle\phi_{k}|\mathcal{C}\rangle$, which
allows us to compute $\Gamma_{k}$ by evaluating the overlap between the
normalized primary state $|\phi_{k}\rangle$ and the crosscap state
$|\mathcal{C}\rangle$. Recalling that the transfer matrix on a cylinder is
given by $e^{-2\pi(L_{0}+\bar{L}_{0})}$, where $L_{0}$ and $\bar{L}_{0}$
represent the holomorphic and anti-holomorphic components of the dilatation
operators, respectively, its eigenvectors correspond to the primary operators
and their descendants. Importantly, the leading eigenvector $|i_{0}\rangle$ in
the previous sections is associated with the identity operator $\mathbb{I}$.
Consequently, we derive the following relation:

$\displaystyle\Gamma_{\mathbb{I}}^{2}$ $\displaystyle=|\langle
0|\mathcal{C}\rangle|^{2},$ $\displaystyle=g.$

To compute $\Gamma_{k}$ for a generic $k$, we utilize the methodology outlined
in Sec.II.2.1 and apply it to additional eigenvectors. Specifically, each
eigenvector of the transfer matrix, constructed following the procedure
described in Sec.II.2.1, possesses two indices. The overlap with the crosscap
state is determined by contracting these indices similar to (10). To
illustrate this process, we apply our method to the Ising model. Within this
context, the second leading and third leading eigenvectors correspond to the
magnetic operator $\sigma$ and the energy operator $\epsilon$, respectively,
in the Ising CFT. The numerical results are summarized in Fig. 7.

41664$0.0$$0.5$$1.0$$1.5$$2.0$$\beta$$\Gamma_{k}^{2}$$\Gamma_{\mathbb{I}}^{2}=1+1/\sqrt{2}$$\Gamma_{\sigma}^{2}=0$$\Gamma_{\epsilon}^{2}=1-1/\sqrt{2}$

Figure 7: $\Gamma_{k}$ obtained from the Ising model by HOTRG with the bond
dimension $\chi=45$. The crosscap boundary state is constructed in the way
described in Fig. 2.

$\Gamma_{k}$ of the Ising model has analytical solutions:
$\Gamma_{\mathbb{I}}^{2}=\frac{2+\sqrt{2}}{2}$, $\Gamma_{\sigma}^{2}=0$, and
$\Gamma_{\epsilon}^{2}=\frac{2-\sqrt{2}}{2}$. Our numerical results are highly
consistent with the theoretical values, denoted with dotted lines. Similarly,
computing the same quantities on other lattice models enables direct access to
the universal information about non-orientable manifolds $\Gamma_{k}$, which
highlights the the power of our scheme.

##  IV Conclusion and Discussion

In our work, we introduced an novel method for calculating the partition
functions on non-orientable surfaces, fundamentally advancing tensor network
simulations. The core of our approach hinges on maintaining data on boundary
configurations during the coarse-graining steps by tracking a spatial
reflection operator, denoted as $\mathbf{\it O}$. While a standard leg
contraction reflects periodic boundary conditions, additionally incorporating
$\mathbf{\it O}$ enables the theoretical construction of partition functions
across diverse manifolds. For instance, the orientable genus-$n$ surface can
be realized by contracting (sewing) $2n$ renormalized tensors, and we can
extend this idea to non-orientable ones by inserting $\mathbf{\it O}$ when
sewing the tensors 333This will be discussed elsewhere.. This development not
only broadens the applicability of tensor network simulations but also opens
new avenues for exploring the physics of complex geometries and topologies. As
a practical application, we showcase the calculation of the one-point function
on the non-orientable surface $\mathrm{RP^{2}}$, addressing a previously
unexplored aspect of universal information within tensor network research. By
synthesizing our method with insights from prior studies, we establish a
comprehensive framework capable of deriving all critical data pertinent to CFT
from tensor network analyses. This achievement underscores the potential of
our approach to contribute significantly to the understanding of CFT data and
its implications for phase transitions and critical phenomena in statistical
physics and condensed matter theory.

## Acknowledgements

A. U. is supported by the MERIT-WINGS Program at the University of Tokyo, the
JSPS fellowship (DC1). He was supported in part by MEXT/JSPS KAKENHI Grants
No. JP21J20523. A part of the computation in this work has been done using the
facilities of the Supercomputer Center, the Institute for Solid State Physics,
the University of Tokyo.

## References

  * Ryu and Zhang [2012] S. Ryu and S.-C. Zhang, Phys. Rev. B 85, 245132 (2012).
  * Sule  _et al._ [2013] O. M. Sule, X. Chen, and S. Ryu, Phys. Rev. B 88, 075125 (2013).
  * Han  _et al._ [2017] B. Han, A. Tiwari, C.-T. Hsieh, and S. Ryu, Phys. Rev. B 96, 125105 (2017).
  * Fioravanti  _et al._ [1994] D. Fioravanti, G. Pradisi, and A. Sagnotti, Physics Letters B 321, 349–354 (1994).
  * Tu [2017] H.-H. Tu, Phys. Rev. Lett. 119, 261603 (2017).
  * Tang  _et al._ [2017] W. Tang, L. Chen, W. Li, X. C. Xie, H.-H. Tu, and L. Wang, Phys. Rev. B 96, 115136 (2017).
  * Kaneda and Okabe [2001] K. Kaneda and Y. Okabe, Phys. Rev. Lett. 86, 2134 (2001).
  * Lu and Wu [2001] W. T. Lu and F. Y. Wu, Phys. Rev. E 63, 026107 (2001).
  * Tang  _et al._ [2019] W. Tang, X. C. Xie, L. Wang, and H.-H. Tu, Phys. Rev. B 99, 115105 (2019).
  * Wang  _et al._ [2018] H.-X. Wang, L. Chen, H. Lin, and W. Li, Phys. Rev. B 97, 220407 (2018).
  * Li  _et al._ [2020] Z.-Q. Li, L.-P. Yang, Z. Y. Xie, H.-H. Tu, H.-J. Liao, and T. Xiang, Phys. Rev. E 101, 060105 (2020).
  * Chen  _et al._ [2017] L. Chen, H.-X. Wang, L. Wang, and W. Li, Phys. Rev. B 96, 174429 (2017).
  * Li and Yang [2021] H. Li and L.-P. Yang, Phys. Rev. E 104, 024118 (2021).
  * Maloney and Ross [2016] A. Maloney and S. F. Ross, Classical and Quantum Gravity 33, 185006 (2016).
  * Xie  _et al._ [2012] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Physical Review B 86, 10.1103/physrevb.86.045139 (2012).
  * Levin and Nave [2007] M. Levin and C. P. Nave, Phys. Rev. Lett. 99, 120601 (2007).
  * Gu and Wen [2009] Z.-C. Gu and X.-G. Wen, Phys. Rev. B 80, 155131 (2009).
  * Evenbly and Vidal [2015] G. Evenbly and G. Vidal, Phys. Rev. Lett. 115, 180405 (2015).
  * Evenbly [2017] G. Evenbly, Phys. Rev. B 95, 045117 (2017).
  * Yang  _et al._ [2017] S. Yang, Z.-C. Gu, and X.-G. Wen, Phys. Rev. Lett. 118, 110504 (2017).
  * Bal  _et al._ [2017] M. Bal, M. Mariën, J. Haegeman, and F. Verstraete, Phys. Rev. Lett. 118, 250602 (2017).
  * Hauru  _et al._ [2018] M. Hauru, C. Delcamp, and S. Mizera, Phys. Rev. B 97, 045111 (2018).
  * Homma and Kawashima [2023] K. Homma and N. Kawashima, Nuclear norm regularized loop optimization for tensor network (2023), arXiv:2306.17479 [cond-mat.stat-mech] .
  * Ueda and Oshikawa [2023] A. Ueda and M. Oshikawa, Physical Review B 108, 10.1103/physrevb.108.024413 (2023).
  * Evenbly and Vidal [2016] G. Evenbly and G. Vidal, Phys. Rev. Lett. 116, 040401 (2016).
  * Li  _et al._ [2022] G. Li, K. H. Pai, and Z.-C. Gu, Phys. Rev. Res. 4, 023159 (2022).
  * Ueda and Yamazaki [2023] A. Ueda and M. Yamazaki, Fixed-point tensor is a four-point function (2023), arXiv:2307.02523 [hep-th] .
  * Fishman  _et al._ [2020] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020), arXiv:2007.14822 [cs.MS] .
  * Note [1] Note that the order $n$ of $O^{n}$ is crucial. The indices $x_{1}$ and $x_{2}$ of $L^{(n)}_{x_{1}x_{2}}$ are coming from the left-hand side index of $T^{(n-1)}$, which is originally the truncated index of $U^{(n-1)}$. Because $O^{(n-1)}$ also has indices of $U^{(n-1)}$, $L^{(n)}$ has to be contracted with $O^{(n-1)}$.
  * Note [2] One might expect that the crosscap free energy can be calculated also in this case using the spatial reflection operator, but it does not work. The reason would be that the isometry ${\mathaccent 869{U}}^{(n)}_{xx_{1}x_{2}}$ becomes highly symmetric in $y$-direction (to the exchange between $x_{1}\leftrightarrow x_{2}$) and the information of reflection is tend to drop by truncation.
  * Note [3] This will be discussed elsewhere.