# A Unifying Variational Framework for
Gaussian Process Motion Planning

Lucas Cosier$*1,2$ Rares Iordan$*1$ Sicelukwanda Zwane$1$ Giovanni Franzese$6$
James T Wilson$3$ Marc Peter Deisenroth$1$ Alexander Terenin$3,4,5$ Yasemin
Bekiroğlu$1,7$

$1$University College London $2$ETH Zürich $3$Imperial College London
$4$University of Cambridge $5$Cornell University $6$TU Delft
$7$Chalmers University of Technology

###### Abstract

To control how a robot moves, motion planning algorithms must compute paths in
high-dimensional state spaces while accounting for physical constraints
related to motors and joints, generating smooth and stable motions, avoiding
obstacles, and preventing collisions. A motion planning algorithm must
therefore balance competing demands, and should ideally incorporate
uncertainty to handle noise, model errors, and facilitate deployment in
complex environments. To address these issues, we introduce a framework for
robot motion planning based on variational Gaussian Processes, which unifies
and generalizes various probabilistic-inference-based motion planning
algorithms. Our framework provides a principled and flexible way to
incorporate equality-based, inequality-based, and soft motion-planning
constraints during end-to-end training, is straightforward to implement, and
provides both interval-based and Monte-Carlo-based uncertainty estimates. We
conduct experiments using different environments and robots, comparing against
baseline approaches based on the feasibility of the planned paths, and
obstacle avoidance quality. Results show that our proposed approach yields a
good balance between success rates and path quality.

##  1 Introduction

Motion planning refers to the process of finding a motion path from a start to
a goal state. This path should avoid obstacles in the environment, satisfy
constraints, such as joint or torque limits, maintain a desired end-effector
orientation, and fulfill other task-specific requirements. The path should
also be smooth, avoiding abrupt or otherwise sharp motions. This can involve
difficulties such as working with potentially high-dimensional robot
configuration spaces and complex environment dynamics.

Modern approaches handle these challenges in different ways. Most motion
planning algorithms can be divided into three broad categories: methods based
on _sampling_ [13, 12, 2, 3], _optimization_ [31, 11, 23, 22, 18, 15], and
_probabilistic inference_ [26, 1, 28, 27]. These categories are formulated
using different mathematical principles. Sampling-based planners, such as
rapidly-exploring random trees [14], work by exploring the space of possible
trajectories via stochastic processes, such as space-filling trees.
Optimization-based planners [21, 23, 22] calculate a motion plan by solving a
constrained optimization problem defined on the space of possible paths.
Probabilistic-inference-based motion planners reinterpret the motion planning
as a Bayesian inference problem over the space of paths, enabling motion plans
to be computed using inference algorithms.

Motion planning algorithms are typically designed to find a single solution
that satisfies certain constraints. However, in many real-world applications
it is useful to reason about the uncertainty associated with the robot’s
state, environment, and actions. Incorporating uncertainty can improve
robustness and help account for sensor noise and model errors. Uncertainty can
also enable the robot to identify the most informative actions for gaining
knowledge about its environment and quantify the risk associated with the
robot’s actions, which in turn can be used to identify the most appropriate
actions to minimize risk. To be able to account for uncertainty in the motion
plan, probabilistic-inference-based planners model uncertainty by way of
producing a distribution of possible motion plans rather than a single
solution.

In this work, we introduce a unifying _variational Gaussian process motion
planning (vGPMP)_ framework, which represents motion plans as random functions
that map time into robot’s configuration space. Our framework represents the
motion planning problem as a _variational inference problem_ —an optimization
problem on a suitable space of probability distributions. This connection
enables Gaussian process motion planners to be viewed as stochastic extensions
of optimization-based planners, which therefore share many of the same
capabilities—but, owing to their stochasticity, also produce uncertainty
estimates that can be propagated upstream.

Our framework makes it possible to directly incorporate both hard and soft
constraints, including joint limits and collision constraints, into the motion
planning objective as well as task-independent properties, such as smoothness.
We enable both interval-based and Monte-Carlo-sampling-based quantification of
uncertainty, each in a scalable manner. In these ways, vGPMP increases the
technical capabilities of Gaussian-process-based motion planners compared to
prior work [6, 10]. Overall, we obtain a simple and straightforward-to-
implement way to construct motion plans with uncertainty.

##  2 Motion Planning using Gaussian Processes

Gaussian processes (GPs) are random functions
$f:\mathcal{X}\rightarrow\mathbb{R}$, where the output value of the function
at any finite set of input variables $\boldsymbol{x}\in\mathcal{X}^{n}$
follows a Gaussian distribution. The process is characterized by a mean
function $\mu:\mathcal{X}\rightarrow\mathbb{R}$ and a positive semi-definite
kernel $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, such that, if
$f\sim\mathcal{GP}(\mu,k)$, then
$\boldsymbol{f}=f(\boldsymbol{x})\sim\mathcal{N}(\boldsymbol{\mu},\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})$
with mean $\boldsymbol{\mu}=\mu(\boldsymbol{x})$ and covariance
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}=k(\boldsymbol{x},\boldsymbol{x})$.
Let $(\boldsymbol{x},\boldsymbol{y})$ be the data, and to ease notation assume
the prior mean is $\mu=0$.

If we consider the Gaussian likelihood
$p(\boldsymbol{y}\mid\boldsymbol{x})=\mathcal{N}(f(\boldsymbol{x}),\sigma^{2}\mathbf{I})$
where $\sigma^{2}>0$ is the noise variance, then the conditional distribution
of the random function $f$ given the training dataset
$(\boldsymbol{x},\boldsymbol{y})$ is also a Gaussian process. Letting
$\mathbf{K}_{(\cdot)\boldsymbol{x}}=k(\cdot,\boldsymbol{x})$, we have
$f\mid\boldsymbol{y}\sim\mathcal{GP}(\mu_{f\mid\boldsymbol{y}},k_{f\mid\boldsymbol{y}})$
where

$\displaystyle\mu_{f\mid\boldsymbol{y}}$
$\displaystyle=\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I}_{n})^{-1}\boldsymbol{y}$
$\displaystyle k_{f\mid\boldsymbol{y}}$
$\displaystyle=k(\cdot,\cdot^{\prime})-\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{K}_{\boldsymbol{x}(\cdot^{\prime})}.$
(1)

To train the GP, the hyperparameters of the kernel need to be estimated. This
is typically done by maximizing the marginal likelihood using standard
gradient-based optimization methods [20].

In the context of motion planning, the input space will always represent time,
namely $\mathcal{X}=\mathcal{T}=\mathbb{R}$, and output space will represent
unconstrained joint values and possibly their derivatives, which we write as
$\Theta=\mathbb{R}^{d}$. Our GPs will therefore be random functions
$f:\mathcal{T}\to\Theta$.

###  2.1 Gaussian Process Motion Planning

By modeling the distribution of feasible trajectories as functions that map
time to robot states, GPs provide a principled approach for formulating the
motion planning problem through the lens of probabilistic inference. To do so
[6] propose to choose a GP prior which is given by a _linear time-varying
stochastic differential equation_ , namely

$\displaystyle\mathop{}\\!\mathrm{d}\boldsymbol{f}(t)$
$\displaystyle=(\mathbf{A}(t)\boldsymbol{f}(t)+\boldsymbol{u}(t))\mathop{}\\!\mathrm{d}t+\mathbf{F}(t)\mathop{}\\!\mathrm{d}\mathbf{w}(t),$
(2)

where $\mathbf{A}(t)$ and $\mathbf{F}(t)$ are time-varying matrices,
$\boldsymbol{u}(t)$ is the control input to the system, and $\mathbf{w}(t)$ is
a Wiener process. The solution of this equation defines the GP prior. The
primary motivations for this choice are that it (i) encodes spatial
smoothness, and (ii) gives rise to structured covariance matrices, which are
block-tridiagonal because the stochastic differential equation (SDE) possesses
a Markov property. This accelerates computation and enables fast maximum a
posteriori inference [6].

The second part of the model used by [6] is given by the likelihood, which
encourages the robot to avoid collisions and constraint violations. Let
$\operatorname{k}_{\operatorname{fwd}}:\Theta\rightarrow\Psi$ be the forward
kinematics function, where $\Psi$ represents joint locations and orientations.
Given the forward kinematics output, following [21], they approximate the
robot’s body by a set of $s$ spheres, and calculate a signed distance field
for each sphere: let $\operatorname{sdf}_{s}:\Psi\rightarrow\mathbb{R}^{s}$ be
the function representing this calculation. Let
$\operatorname{h}_{\varepsilon}(x)=\max(-x+\varepsilon,0)$ be the hinge loss,
where $\varepsilon>0$ is called the _safety distance_ parameter. Let
$\operatorname{c}$ be a linear penalty function for values exceeding a given
joint limit threshold. By convention, $\operatorname{h}_{\varepsilon}$ and
$\operatorname{c}$ act component-wise on vectors. Let
$\left\lVert\boldsymbol{x}\right\rVert_{\mathbf{\Sigma}}^{2}=\boldsymbol{x}^{T}\mathbf{\Sigma}^{-1}\boldsymbol{x}$
be the Mahalonobis norm, and let $\mathbf{\Sigma}_{\operatorname{obs}}$ and
$\mathbf{\Sigma}_{\operatorname{c}}$ be diagonal matrices. With this notation,
[6] propose the likelihood

$p(\boldsymbol{e}\mid\boldsymbol{f})\propto\exp\bigg{(}-\frac{1}{2}\Big{(}\,\underbracket{\left\lVert\operatorname{h}_{\varepsilon}(\operatorname{sdf}_{s}(\operatorname{k}_{\operatorname{fwd}}(\boldsymbol{f})))\right\rVert_{\mathbf{\Sigma}_{\operatorname{obs}}}^{2}}_{\vphantom{|}\mathclap{\mathrel{\text{missing}}{collisionterm}}\vphantom{|}}+\underbracket{\left\lVert\operatorname{c}(\boldsymbol{f})\right\rVert^{2}_{\mathbf{\Sigma}_{\operatorname{c}}}}_{\vphantom{|}\mathclap{\mathrel{\text{missing}}{softconstraintterm}}\vphantom{|}}\Big{)}\bigg{)}$
(3)

where $\boldsymbol{e}$ represents the probability of _collision-free events
with motion constraints_. From here, [6] approximate the posterior
distribution as follows: they (i) apply maximum a posteriori inference,
deriving efficient algorithms which use the SDE’s Markov property to recast
the resulting optimization objective as inference in a sparse factor graph.
Once such a trajectory is obtained, they (ii) apply Laplace approximations to
approximate the posterior distribution, and (iii) interpolate joint
configuration in-between the states at times $t_{1},..,t_{T}\in\mathcal{T}$.
At run-time, the output values of $\boldsymbol{f}$ are clamped to enforce
constraints. In total, the computations performed produce an optimal
trajectory $\boldsymbol{f}^{*}$, and a Gaussian distribution representing
trajectories around it that quantifies uncertainty.

##  3 Variational Gaussian Process Motion Planning

Computing motion plans using GPs involves balancing multiple competing
factors, including obstacle avoidance, smoothness, and conformity to joint
limits. These factors depend simultaneously on the motion planner’s
hyperparameters, which must be tuned, and on the formulation of the objective
itself. To provide practitioners with precise control over what the planner
produces and enable them to modify it as needed to introduce other objectives
for specific tasks, we provide a unifying mathematical and algorithmic
formalism for Gaussian process motion planning.

###  3.1 Sparse Gaussian Processes as a Unifying Framework

From a probabilistic inference perspective, one of the main challenges in
Gaussian Process motion planning is the need to accommodate non-Gaussian
likelihoods. Due to these likelihoods, the motion planning posterior, unlike
the prior, is not a GP. This makes it impossible to perform exact
probabilistic inference and necessitates the use of numerical methods to
extract means, covariances, and Monte Carlo samples from the posterior
distribution.

_Sparse Gaussian Processes_ [8, 9, 25] are a state-of-the-art method for
applying GPs in settings with non-Gaussian likelihoods. They approximate the
exact non-Gaussian posterior with the best Gaussian approximation, in the
sense of Kullback–Leibler divergence. Originally proposed to tackle
scalability issues, which arise when working with large datasets, they provide
a broad and flexible framework for designing approximate probabilistic
inference schemes using GPs. We now apply these notions to design a unifying
framework for motion planning.

We start with an _arbitrary_ Gaussian process prior
$f\sim\mathcal{GP}(\mu,k)$, generalizing the framework of [6] beyond
stochastic differential equations. Thus, for any collection of times
$\mathbf{T}=[t_{0},t_{1},\dots,t_{n}]$, the random vector
$\boldsymbol{f}=[\boldsymbol{f}(t_{0}),\boldsymbol{f}(t_{1}),\dots,\boldsymbol{f}(t_{n})]$
follows a joint Gaussian distribution whose mean and covariance are determined
by $\mu$ and $k$. To ease exposition, suppose that we are also given a
likelihood $p(\boldsymbol{e}\mid\boldsymbol{f})$, which incorporates
constraints and collision avoidance, such as the one presented in Section
2.1—we will discuss likelihoods later in Section 3.2. These ingredients fully
specify the exact posterior distribution, which is non-Gaussian.

To handle non-Gaussianity, following the sparse GP framework [24, 19, 25], we
introduce the _variational approximation_ $f\mid\boldsymbol{u}$, which is a GP
conditioned on the event $f(\boldsymbol{z})=\boldsymbol{u}$, where
$\boldsymbol{z}\in\mathcal{T}^{m}$ are a set of locations called the _inducing
points_ , and $\boldsymbol{u}$ are called the _inducing variables_. Let
$q(\boldsymbol{u})\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{u}},\mathbf{\Sigma}_{\boldsymbol{u}})$
be the distribution of the $\boldsymbol{u}$-values, which we assume to be a
free-form multivariate Gaussian whose mean and covariance
$\boldsymbol{\mu}_{\boldsymbol{u}}$ and $\mathbf{\Sigma}_{\boldsymbol{u}}$
will be learned using optimization. The gives
$f\mid\boldsymbol{u}\sim\mathcal{GP}(\mu_{f\mid\boldsymbol{u}},k_{f\mid\boldsymbol{u}})$
with

$\displaystyle\mu_{f\mid\boldsymbol{u}}$
$\displaystyle=\mu(\cdot)+\mathbf{K}_{(\cdot)\boldsymbol{z}}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{-1}(\boldsymbol{u}-\boldsymbol{\mu}_{\boldsymbol{u}})$
$\displaystyle k_{f\mid\boldsymbol{u}}$
$\displaystyle=k(\cdot,\cdot^{{}^{\prime}})-\mathbf{K}_{(\cdot)\boldsymbol{z}}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{-1}\mathbf{K}_{\boldsymbol{z}(\cdot)}$
(4)

To learn the parameters of $q(\boldsymbol{u})$, we minimize the
Kullback–Leibler divergence between $f\mid\boldsymbol{u}$ and the true
posterior. One can show [16] that this is equivalent to maximizing the
_evidence lower bound_

$\operatorname{ELBO}=\operatorname{\mathbb{E}}_{q(\boldsymbol{u})}\operatorname{\mathbb{E}}_{p(\boldsymbol{f}\mid\boldsymbol{u})}\log
p(\boldsymbol{e}\mid\boldsymbol{f})-D_{\operatorname{KL}}(q(\boldsymbol{u})\mid\mid
p(\boldsymbol{u})),$ (5)

which we optimize by using stochastic-gradient-based techniques. This defines
our general framework for probabilistic-inference-based motion planning.
Figure 1 illustrates how inducing points affect possible trajectories of the
variational posterior. From here, the next steps are to (ii) fully specify the
likelihood $p(\boldsymbol{e}\mid\boldsymbol{f})$ to incorporate factors, such
as constraints and collision avoidance, and (ii) understand how to extract
uncertainty for use by upstream systems. We now study these aspects.

(a) Sparse Gaussian process

(b) Motion plans via GP samples

Figure 1: Intuitive illustration of representing motion paths with GPs. (1(a))
Sparse GP posterior for 1 joint with inducing variables
$\boldsymbol{u}=[\boldsymbol{u}_{\operatorname{c}},\boldsymbol{u}^{\prime}]$
denoted by ‘$\circ$’, where
$\boldsymbol{u}_{\operatorname{c}}=\left[\boldsymbol{\theta}_{0},\boldsymbol{\theta}_{1}\right]$
and $\boldsymbol{u}^{\prime}$ represent the $M=4$ inducing locations
$\boldsymbol{z}\in\mathcal{T}$ subject to the learning process, shown in-
between $t=0$ and $t=1$. (1(b)) Candidate end-effector trajectories for the
Franka Panda robot and associated uncertainty, which is illustrated as the
orange shaded region.

###  3.2 Incorporating Constraints, Smoothness, and Collision Avoidance

A motion planner needs to produce motion plans that satisfy a number of
constraints, such as the desired start and end states and collision avoidance.
Within the variational Gaussian Process motion planning framework, this
amounts to incorporating information into the posterior distribution using
conditional probability, in most cases via the likelihood
$p(\boldsymbol{e}\mid\boldsymbol{f})$. We now describe how this works for a
number of different forms of information.

#### Joint Limits and Inequality-based Constraints

To respect the robot’s joint limits, we need to incorporate an inequality
constraint into the GP model. Previous, non-variational Gaussian Process
motion planning approaches [6] generally enforce such constraints softly by
adding an extra term to the optimization objective. This does not guarantee
the constraints are satisfied, which necessitates post-processing. The most
natural alternative, which guarantees constraints to be satisfied, would be to
modify the prior using conditional probability, but this unfortunately breaks
Gaussianity of the GP. We therefore instead propose to leverage the
variational formulation and enforce constraints using a bijective non-linear
transformation. For inequality constraints of the form

$\theta_{i}^{(\min)}\leq\theta^{\vphantom{()}}_{i}\leq\theta_{i}^{(\max)}$ (6)

this can be done by passing the GP through a scaled and shifted sigmoid
function $\sigma:\Theta\rightarrow\Theta_{\operatorname{c}}$, where
$\Theta_{\operatorname{c}}$ is constrained joint space. This amounts to simply
replacing $\boldsymbol{f}$ with $\sigma(\boldsymbol{f})$ in the likelihood,
and returning motion plans based on the transformed GP sample paths, which are
guaranteed to satisfy all joint constraints. This straightforward solution
highlights the value of a variational formulation, which makes such aspects
computationally simple to add.

#### Start States, Goal States, and Equality-based Constraints

To ensure our trajectories pass through the (unconstrained) start and goal
states $\boldsymbol{\theta}_{0},\boldsymbol{\theta}_{1}$ at times $t=0$ and
$t=1$, we condition the variational process $f\mid\boldsymbol{u}$ on the
events $f(0)=\boldsymbol{\theta}_{0}$ and $f(1)=\boldsymbol{\theta}_{1}$. This
is done by partitioning the inducing variables as
$\boldsymbol{u}=[\boldsymbol{u}_{\operatorname{c}},\boldsymbol{u}^{\prime}]$,
where $\boldsymbol{u}_{\operatorname{c}}$ are the inducing variables
representing start and goal states, and $\boldsymbol{u}^{\prime}$ are the
remaining inducing variables. We then simply set
$\boldsymbol{u}_{{\operatorname{c}}}$ to the values needed and exclude them
from optimization, as shown in the posterior in Figure 1 where
$\boldsymbol{u}_{\operatorname{c}}$ were set to $[0,0]$ at $t=0,t=1$ in the 1D
case. This approach allows us to handle general equality-based constraints.

#### Collision Avoidance, Velocity Limits, and Additional Soft Constraints

To incorporate constraints that do not have a convenient mathematical form, we
include likelihood-based soft constraints in a manner that mirrors non-
variational Gaussian Process motion planning algorithms, such as the
collision-avoidance term proposed by [6]. In such cases, note that we replace
$\boldsymbol{f}$ with $\sigma(\boldsymbol{f})$ to account for inequality
constraints as described previously. This enables us to handle collision
avoidance, velocity constraints and other soft constraints, which are
implicitly enforced by penalizing violations via the optimization objective,
providing a flexible modeling tool that can be adapted to the task at hand.
For example, to enforce grasping, we can include the term

$\left\lVert\operatorname{g}_{s}(\operatorname{k}_{\operatorname{fwd}}(\sigma(\boldsymbol{f})))\right\rVert_{\mathbf{\Sigma}_{\operatorname{grasp}}}$
(7)

inside the exponential function in the log likelihood (8). Here,
$\operatorname{g}_{s}$ is the function that maps the forward kinematics to a
_good graspable pose_ , which is characterized by both the position and
orientation of the end effector. The matrix
$\mathbf{\Sigma}_{\operatorname{grasp}}$ plays a role similar to
$\mathbf{\Sigma}_{\operatorname{obs}}$, and quantifies the relative importance
of this term in the likelihood. This yields a combined overall likelihood,
which becomes part of the variational optimization problem that produces the
motion plan.

Figure 2: A sampled trajectory for grasping a can. Frames 1-2, 2-3, 3-4 depict
respectively the end-effector alignment with the final pose, the approach
stage, and the grasping stage.

As a proof-of-concept, we include a grasping experiment. Alignment for viable
grasping is achieved by adding an appropriate soft constraint term (7), which
encourages the arm to orient itself accordingly. One can take this further,
conditioning the GP at appropriate time steps to open and close the gripper of
the robot to move objects from one place to another. This reinforces the value
of the proposed framework, which allows us to treat different kinds of
constraints—here, collision-avoidance and grasping alignment—in essentially-
the-same way mathematically. We demonstrate a motion plan incorporating these
elements in Figure 2, and include a video in the supplementary material.

#### Smoothness

Smoothness can be controlled via kernel choice and its hyperparameters. In
particular, stationary kernels on $\mathbb{R}^{1}$ can be formulated to
include a _length scale_ parameter $\kappa>0$, and written
$k(x,x^{\prime})=k(\frac{|x-x^{\prime}|}{\kappa})$. This parameter determines
the characteristic distances over which the function varies [20]. Kernels on
$\mathbb{R}^{d}$ for $d>1$ can be defined similarly, with one length scale per
dimension. Larger length scales yield smoother trajectories, while smaller
length scales yield rougher, more jagged trajectories. Through kernel choice,
for instance by choosing a Matérn kernel with given smoothness [20], one can
also control properties such as the number of times a motion path is
differentiable. We illustrate this in Figure 3.

#### Computation via Variational Inference and Connections with Optimization-
based Planners

Summarizing, by combining the previous likelihood terms, we obtain the
negative ELBO

$\frac{1}{2}\operatorname{\mathbb{E}}_{q(\boldsymbol{u})}\operatorname{\mathbb{E}}_{p(\boldsymbol{f}\mid\boldsymbol{u})}\Big{(}\,\underbracket{\left\lVert\operatorname{h}_{\varepsilon}(\operatorname{sdf}_{s}(\operatorname{k}_{\operatorname{fwd}}(\sigma(\boldsymbol{f}))))\right\rVert_{\mathbf{\Sigma}_{\operatorname{obs}}}^{2}}_{\vphantom{|}\mathclap{\mathrel{\text{missing}}{collisionterm}}\vphantom{|}}+\underbracket{\left\lVert\operatorname{c}(\sigma(\boldsymbol{f}))\right\rVert^{2}_{\mathbf{\Sigma}_{\operatorname{c}}}}_{\vphantom{|}\mathclap{\mathrel{\text{missing}}{softconstraintterms}}\vphantom{|}}\Big{)}+D_{\operatorname{KL}}(q(\boldsymbol{u}^{\prime})\mid\mid
p(\boldsymbol{u}^{\prime}\mid\boldsymbol{u}_{\operatorname{c}})),$ (8)

which we minimize to obtain the optimal parameters of $q(\boldsymbol{u})$.
This enforces equality constraints through the
$\boldsymbol{u}_{\operatorname{c}}$-values, inequality constraints through
$\sigma$, and soft constraints through the collision and other soft constraint
terms. Compared to previous GP-based approaches such as for instance [6], we
do not need to introduce sparse factor graphs, Laplace approximations,
interpolation between time steps, or other scalability tricks. Instead, we
simply carry minimization out using stochastic optimization, by sampling
$\boldsymbol{u}$-values, conditionally sampling $\boldsymbol{f}$-values,
computing the objective, and taking gradient steps.

From (8), we deduce that probabilistic-inference-based motion planners based
on variational inference are stochastic extensions of traditional,
optimization-based motion planning algorithms such as STOMP and CHOMP [11,
31]. The differences are that: (i) instead of a single trajectory, we now
return a distribution over trajectories, (ii) the loss contains an extra
expectation term, and (iii) we include an additional KL-based regularizer to
ensure we do not stray too far from the prior, thereby allowing us to control
smoothness and related task-independent properties in a distributional manner.

(a) Matérn-$\nicefrac{{1}}{{2}}$

(b) Matérn-$\nicefrac{{3}}{{2}}$

(c) Matérn-$\nicefrac{{5}}{{2}}$

Figure 3: Illustration of smoothness properties of different kernels,
including the nowhere-differentiable Matérn-$\nicefrac{{1}}{{2}}$, once-
differentiable Matérn-$\nicefrac{{3}}{{2}}$, and twice-differentiable
Matérn-$\nicefrac{{5}}{{2}}$ kernel.

#### Replanning

If we need to replan, we do so by starting with the previous motion plan and
moving inducing points and inducing values around. We also condition on the
current value as initial state.

###  3.3 Uncertainty: Intervals and Monte Carlo Function Samples

Once a Gaussian process motion planning algorithm has completed optimization,
the next step is to extract an actual motion plan from it. Since, in the
variational framework, motion planning algorithms produce Gaussian processes,
the _most likely_ motion plan under the variational posterior is the mean of
the variational posterior. We can incorporate uncertainty in two ways: (a)
through posterior intervals or (b) by drawing Monte Carlo random function
samples. These work as follows:

  * (a)

Posterior intervals: for a given scaling level $\alpha>0$, we can represent
uncertainty through posterior intervals, which are
$\mu_{f\mid\boldsymbol{u}}(x)\pm\alpha\sqrt{k_{f\mid\boldsymbol{u}(x,x)}}$.
For large-enough $\alpha$, this gives a high-probability estimate on the
region through which motion plans will travel.

  * (b)

Posterior function samples: we can also draw random function
$\phi\sim\operatorname{GP}(\mu_{f\mid\boldsymbol{u}},k_{f\mid\boldsymbol{u}})$.
For stationary kernels, we do so up to a small error tolerance using the
_efficient sampling_ technique of [29, 30], which works by drawing approximate
random functions from the prior, and transforming them into the posterior.
Once drawn, we can then plug these randomly-sampled functions into upstream
systems as needed for the given task.

We emphasize that (b) is computationally tractable in part because the
variational formulation can handle general kernels. These technical primitives
can be used to propagate uncertainty into any upstream decision-making
algorithm built atop of GP. Gaussian processes are widely used to represent
uncertainty in a number of decision-making systems, including in Bayesian
optimization, active learning, and reinforcement learning. Our formulation can
therefore in principle facilitate the development of similar systems in the
context of motion planning, with further work.

###  3.4 Limitations

Compared to alternatives, such as optimization-based approaches, the vGPMP
framework includes uncertainty, which means that more effort can be required
to find algorithmic hyperparameters that produce good performance. Since this
uncertainty comes from GPs, we therefore also inherit all of the usual
limitations of GP-based techniques, including the need to solve linear
systems, which can present numerical difficulties. Our optimization objective
is a sum of multiple terms that reflect complementary aspects of motion
planning, which need to be weighted appropriately—empirically, we found it to
be important that the entries of $\mathbf{\Sigma}_{\operatorname{obs}}$ in (8)
are large enough, to ensure the motion planner does not ignore the collision
part of the objective.

Compared to other probabilistic-inference-based planners, such as GPMP2, our
formulation ensures computational tractability in a different way, namely
through sparse GPs rather than through a factor graph arising from the
stochastic differential equation prior. This means that we cannot take
advantage of fast inference algorithms for factor graphs and rely instead on
variational inference for performance, making vGPMP reliant on stochastic
optimization algorithms. Our algorithm’s complexity is
$\mathcal{O}(TM^{2}+dM^{3}+T(d+s))$ for obtaining the variational parameters
which define a motion plan, $\mathcal{O}(dM^{3})$ for computing error bars,
and $\mathcal{O}(TM^{2}+M^{3})$ for both computing the mean motion plan and
drawing one Monte Carlo sample when using efficient sampling [29, 30] and
$\mathcal{O}(T(d+s))$ to compute its likelihood, where $T$ is the number of
time points, $M$ is the number of inducing points, which is usually a dozen or
so, $d$ is the dimensionality of joint space, and $s$ is the number of
spheres. Derivations are given in Appendix B.

##  4 Experiments

(a) Start query state

(b) Goal query state

Figure 4: Example distribution of collision-free paths generated by vGPMP in
the boxes environment. A sample is shown in green.

M | vGPMP (ours) | GPMP2
---|---|---
Acc | $\boldsymbol{100\%}$ | $\boldsymbol{100\%}$
MC | $\boldsymbol{-0.90\pm 1.50}$ | $-2.58\pm 2.54$
MPL | $2.17\pm 1.00$ | $1.89\pm 0.98$

Table 1: Comparison using WAM robot in the lab environment. Metrics are
_accuracy_ , _mean clearance_ , and _mean path length_.

We now evaluate vGPMP and compare its performance and capabilities with other
motion planners, including both those based on probabilistic inference and
those based on other formalisms. To illustrate the unifying nature of our
techniques, we consider standard collision-avoidance benchmarks. In all cases,
full experiment details can be found in Appendix C.

Robot | Env | M | Planners |
---|---|---|---|---
vGPMP (ours) | RRTconnect | LBKPIECE | CHOMP |
Franka | I | Acc | $\boldsymbol{100\%}$ | $\boldsymbol{100\%}$ | $88.89\%$ | $18.06\%$ |
MPL | $2.75\pm 2.14$ | $1.16\pm 1.67$ | $1.82\pm 2.48$ | $0.96\pm 1.95$ |
MC | $\boldsymbol{-1.17\pm 1.73}$ | $-11.81\pm 3.23$ | $-14.03\pm 4.55$ | $-9.25\pm 6.94$ |
B | Acc | $81.81\%$ | $\boldsymbol{100\%}$ | $\boldsymbol{100\%}$ | $14.54\%$ |
MPL | $2.60\pm 1.86$ | $2.10\pm 2.09$ | $2.21\pm 2.19$ | $1.27\pm 2.01$ |
MC | $\boldsymbol{-6.76\pm 4.69}$ | $-84.42\pm 18.63$ | $-116.58\pm 20.05$ | $-37.47\pm 12.77$ |
UR10 | I | Acc | $\boldsymbol{100\%}$ | $97.22\%$ | $86.11\%$ | $11.11\%$ |
MPL | $1.63\pm 1.09$ | $2.42\pm 2.55$ | $4.03\pm 4.37$ | $0.20\pm 0.47$ |
MC | $\boldsymbol{-56.03\pm 11.86}$ | $-281.84\pm 22.91$ | $-479.73\pm 35.52$ | $-5940.77\pm 108.31$ |
B | Acc | $\boldsymbol{89.82\%}$ | $85.45\%$ | $16.36\%$ | $18.18\%$ |
MPL | $3.95\pm 1.91$ | $1.82\pm 2.22$ | $0.93\pm 1.56$ | $0.55\pm 0.90$ |
MC | $\boldsymbol{-7.10\pm 5.55}$ | $-42.56\pm 13.03$ | $-26.48\pm 9.97$ | $-58.53\pm 6.58$ |
Table 2: Experiments in different environments using two robots. The
environments (Env) used are named Industrial (I) and Bookshelves (B), and the
Measures (M) are accuracy (Acc), Mean Path Length (MPL), and Mean Clearance
(MC).

To understand how vGPMP and other motion planners handle collision avoidance,
we ran them on two kinds of motion planning benchmarks. For the first kind, we
compare with a state-of-the-art probabilistic-inference-based planner _GPMP2_
, using similar hyperparameters as detailed in Appendix C, on the Barrett
_WAM_ robot (7 DOF) in the _Lab_ environment [18]. For the second kind, we
compare with sampling-based and optimization-based planners, using benchmarks
originally introduced by [23], which are based on the _Franka Panda_ (7 DOF)
and _UR10_ (6 DOF) robots in the _Industrial_ and _Bookshelves_ environments
(which can be seen in Appendix C), and have 36 and 55 unique motion planning
problems, respectively.

We evaluate behavior with respect to (i) _accuracy_ , which is the percentage
of tasks within the benchmark suite correctly solved by the planner; (ii)
_clearance_ , which quantifies collision-avoidance by measuring the distance
to objects; (iii) _path length_ , which is the total distance moved by the end
effector. All measurements are repeated using different random seeds to assess
variability. Further details on measurements are given in Appendix C.

Table 1 shows that vGPMP and GPMP2 both succeed in finding obstacle-free paths
reaching the targets in all cases. Compared to GPMP2, vGPMP leads to better
clearance values and consequently slightly longer paths moving away from the
objects ensuring safety. Both methods require careful tuning of
$\mathbf{\Sigma}_{\operatorname{obs}}$ to ensure that this term is a
sufficiently-important part of the loss to achieve collision-avoidance.
Overall, in comparison to GPMP2, vGPMP performs comparably and provides a
simpler computational pipeline and more general formulation.

Next, we examine how uncertainty translates to real-world execution. For
vGPMP, we use Monte Carlo sampling to obtain 150 paths, from which we
calculate the resulting envelope around the mean trajectory. The intervals in
Figure 4 show the successful incorporation of the problem constraints as
expressed in the low uncertainty at the beginning and end of the trajectory.
This uncertainty grows as one moves further away from start and goal states,
as well as points of possible collision.

Table 2 shows that vGPMP consistently leads to high success rates in reaching
the goal states, achieving 91.42% across all experiments, which is
significantly higher than the performance of LBKPIECE (69.78%) and CHOMP
(15.69%). RRTconnect achieves a slightly higher success rate, 95.05% on
average—however, the generated paths demonstrate worse clearance values. In
most cases, vGPMP produces longer paths, since uncertainty drives paths away
from obstacles for better collision avoidance.

vGPMP consistently yields better clearance values in all environments using
the two robots, except with _Franka Panda_ in the _bookshelves_ environment,
where the performance of vGPMP drops due to collisions with the floor: this
could be alleviated by including an extra term in the loss, which we omitted
to ensure a fair comparison with other planners that also do not consider the
floor. In total, without post-processing, the uncertainty produced by vGPMP
results in safer trajectories compared to those of other planners.

We then apply our approach directly to a real robot. We demonstrate the
feasibility of the approach by executing a trajectory generated from vGPMP
given three obstacles placed around a Franka Panda robot. To execute the
generated motion plans on the physical robot, we used a low-level joint
impedance position controller [7]. We show the resulting trajectory in Figure
5, and the experiment can be seen in the supplementary video.

Figure 5: Sampled trajectory executions on the real (bottom) and the simulated
(top) robot.

##  5 Conclusion

We present a unifying framework for Gaussian-process-based motion planning
with obstacle avoidance, using a formulation based on variational inference.
The framework supports equality constraints, inequality constraints, and
general soft constraints. It also allows for fine-grained control over
smoothness through prior specification. Our approach generalizes previous
Gaussian-process-based motion planners to arbitrary kernels in a manner that
simplifies computation. At the same time, by viewing variational inference as
an optimization problem over a space of random motion paths, our framework
also generalizes optimization-based planners to incorporate uncertainty. We
evaluate the approach using a number of robots, environments, and
tasks—including obstacle avoidance and grasping—showing that it has high
accuracy in reaching target positions, while avoiding obstacles and using
uncertainty to providing better clearance compared to baselines. We
demonstrate the feasibility of the approach on a real robot, executing a
sampled trajectory while avoiding obstacles.

## References

  * [1] Hagai Attias  “Planning by probabilistic inference.”  In _Artificial Intelligence and Statistics_ , 2003
  * [2] Dmitry Berenson, Siddhartha Srinivasa and James Kuffner  “Task space regions: A framework for pose-constrained manipulation planning”  In _International Journal of Robotics Research_ , 2011
  * [3] Peng Cheng and Steven M LaValle  “Resolution complete rapidly-exploring random trees”  In _International Conference on Robotics and Automation_ , 2002
  * [4] Sachin Chitta, Ioan Sucan and Steve Cousins  “Moveit![ros topics]”  In _IEEE Robotics & Automation Magazine_ 19.1 IEEE, 2012, pp. 18–19
  * [5] Erwin Coumans and Yunfei Bai  “PyBullet, a python module for physics simulation in robotics, games and machine learning”, 2017
  * [6] Jing Dong, Mustafa Mukadam, Frank Dellaert and Byron Boots  “Motion Planning as Probabilistic Inference using Gaussian Processes and Factor Graphs.”  In _Robotics: Science and Systems_ , 2016
  * [7] Giovanni Franzese, Anna Mészáros, Luka Peternel and Jens Kober  “ILoSA: Interactive learning of stiffness and attractors”  In _Intelligent Robots and Systems_ , 2021
  * [8] James Hensman, Nicolo Fusi and Neil D. Lawrence  “Gaussian Processes for Big Data”  In _Uncertainty in Artificial Intelligence_ , 2013
  * [9] James Hensman, Alex Matthews and Zoubin Ghahramani  “Scalable Variational Gaussian Process Classification”  In _Artificial Intelligence and Statistics_ , 2015
  * [10] Eric Huang, Mustafa Mukadam, Zhen Liu and Byron Boots  “Motion planning with graph-based trajectories and Gaussian process inference”  In _International Conference on Robotics and Automation_ , 2017
  * [11] Mrinal Kalakrishnan et al.  “STOMP: Stochastic trajectory optimization for motion planning”  In _International Conference on Robotics and Automation_ , 2011
  * [12] Sertac Karaman and Emilio Frazzoli  “Sampling-based algorithms for optimal motion planning”  In _International Journal of Robotics Research_ , 2011
  * [13] James Kuffner and Steven M LaValle  “RRT-connect: An efficient approach to single-query path planning”  In _International Conference on Robotics and Automation_ , 2000
  * [14] Steven Lavalle and James Kuffner  “Rapidly-Exploring Random Trees: Progress and Prospects”  In _Algorithmic and Computational Robotics_ , 2000
  * [15] Dieter Fox Lirui Wang  “Manipulation trajectory optimization with online grasp synthesis and selection”  In _Robotics: Science and Systems_ , 2020
  * [16] Alexander G. G. Matthews, James Hensman, Richard Turner and Zoubin Ghahramani  “On Sparse Variational Methods and the Kullback-Leibler Divergence between Stochastic Processes”  In _Artificial Intelligence and Statistics_ , 2016
  * [17] Alexander G. G. Matthews et al.  “GPflow: A Gaussian process library using TensorFlow”  In _Journal of Machine Learning Research_ , 2017
  * [18] Mustafa Mukadam, Xinyan Yan and Byron Boots  “Gaussian process motion planning”  In _International Conference on Robotics and Automation_ , 2016
  * [19] Joaquin Quiñonero-Candela and Carl Rasmussen  “A Unifying View of Sparse Approximate Gaussian Process Regression”  In _Journal of Machine Learning Research_ , 2005
  * [20] Carl Rasmussen and Christopher K.. Williams  “Gaussian Processes for Machine Learning”  MIT Press, 2006
  * [21] Nathan Ratliff, David Silver and J. Bagnell  “Learning to Search: Functional Gradient Techniques for Imitation Learning”  In _Autonomous Robots_ , 2009
  * [22] John Schulman et al.  “Motion planning with sequential convex optimization and convex collision checking”  In _International Journal of Robotics Research_ , 2014
  * [23] John Schulman et al.  “Finding locally optimal, collision-free trajectories with sequential convex optimization.”  In _Robotics: Science and Systems_ , 2013
  * [24] Edward Snelson and Zoubin Ghahramani  “Sparse Gaussian Processes using Pseudo-inputs”  In _Advances in Neural Information Processing Systems_ , 2006
  * [25] Michalis Titsias  “Variational Learning of Inducing Variables in Sparse Gaussian Processes”  In _Artificial Intelligence and Statistics_ , 2009
  * [26] Marc Toussaint  “Robot trajectory optimization using approximate inference”  In _International Conference on Machine Learning_ , 2009
  * [27] Marc Toussaint and Christian Goerick  “A Bayesian View on Motor Control and Planning”  In _From Motor Learning to Interaction Learning in Robotics_ , 2010
  * [28] Marc Toussaint and Amos Storkey  “Probabilistic inference for solving discrete and continuous state Markov decision processes”  In _International Conference on Machine Learning_ , 2006
  * [29] James T. Wilson et al.  “Efficiently Sampling Functions from Gaussian Process Posteriors”  In _International Conference on Machine Learning_ , 2020
  * [30] James T. Wilson et al.  “Pathwise Conditioning of Gaussian Processes”  In _Journal of Machine Learning Research_ , 2021
  * [31] Matthew Zucker et al.  “CHOMP: Covariant Hamiltonian Optimization for Motion Planning”  In _International Journal of Robotics Research_ , 2013

##  Appendix A Incorporating Additional Constraints under The Variational
Representation

For learning the variational distribution, we choose
$q(\boldsymbol{u})=\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{u}},\mathbf{\Sigma}_{\boldsymbol{u}})$
with a particular representation. Letting
$\boldsymbol{u}=[\boldsymbol{u}^{\prime},\boldsymbol{u}_{\operatorname{c}}]$
be the partition of $\boldsymbol{u}$ into the learned and constrained parts,
the matrix $\mathbf{\Sigma}_{\boldsymbol{u}}$ contains one part for
$\boldsymbol{u}^{\prime}$, which is padded with zeros to encode the
uncertainty in the given motion planning constraints
$\boldsymbol{u}_{\operatorname{c}}$. One can additionally introduce noise in
the latent variables, but for we omit this here to ease notation.

#### Computational Considerations

To improve computational efficiency, we use a whitened representation that
only considers the learned portion of the covariance
$\mathbf{\Sigma}_{\boldsymbol{u}}$. Let $\boldsymbol{z}^{\prime}$ and
$\boldsymbol{z}_{\operatorname{c}}$ be the inducing points corresponding to
$\boldsymbol{u}^{\prime}$ and $\boldsymbol{u}_{\operatorname{c}}$,
respectively. The part of the distribution representing
$\boldsymbol{u}_{\operatorname{c}}$ is a degenerate Gaussian, making it
advantageous to omit it from the whitening process. This representation is
achieved using the block partitioned representation of
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$, because the Cholesky
decomposition of the covariance of
$p(\boldsymbol{u}^{\prime}|\boldsymbol{u}_{\operatorname{c}})$ can be written
as
$\operatorname{Chol}(\mathbf{K}_{\boldsymbol{z}^{\prime}\boldsymbol{z}^{\prime}}-\mathbf{K}_{\boldsymbol{z}^{\prime}\boldsymbol{z}_{\operatorname{c}}}\mathbf{K}_{\boldsymbol{z}_{\operatorname{c}}\boldsymbol{z}_{\operatorname{c}}}^{-1}\mathbf{K}_{\boldsymbol{z}_{\operatorname{c}}\boldsymbol{z}^{\prime}})$:
this can be seen by considering the block LU decomposition. An added benefit
of this approach is that the KL term in (8) can be computed efficiently, as
the Cholesky factor is already cached and used in the pathwise updates.

#### Velocity Modeling and Higher Order Derivatives

The inclusion of constraints regarding higher order derivatives, such as
maintaining constant velocity, is crucial in motion planning for real-world
execution. By incorporating these constraints, the resulting models can
accurately account for the dynamic limitations and constraints of physical
systems, ensuring the feasibility and reliability of planned trajectories
during execution. In our framework, one can encode velocities as prior
information by considering the inducing variables along with their time
derivatives, namely
$\boldsymbol{u}=[\boldsymbol{u}^{\prime},\dot{\boldsymbol{u}}^{\prime},\boldsymbol{u}_{\operatorname{c}},\dot{\boldsymbol{u}}_{\operatorname{c}}]$
which implies that the respective kernel matrix will have corresponding kernel
derivatives depending on the left or right argument. Typically, the start and
end velocities are set to zero. In addition, in our framework angular
velocities of the samples can easily be computed by differentiating the sample
paths. This means they can be included in the likelihood via a soft constraint
term which bring the velocity of sampled trajectories to desired values. This
gives, for instance, the term

$\big{\lVert}\sigma(\boldsymbol{\dot{f}})-\mu_{\mathbf{v}}\big{\rVert}_{\mathbf{\Sigma}_{\operatorname{velocity}}}.$
(9)

Acceleration and other higher order derivatives can also be considered
following a similar argument, assuming appropriate kernel choice to ensure
sufficient differentiability. For the Matérn-$\nicefrac{{5}}{{2}}$ kernel,
this benefit is two-fold as the GP is then Markovian, which can be used to re-
formulate the necessary computations in terms of sparse linear algebra, which
makes it possible to generalize the stochastic-differential-equation-based
techniques of [6] to GP priors with more available derivatives.

##  Appendix B Time Complexity

In this section, we examine the scalability of vGPMP. Let $s$ be the number of
spheres that model the robot body. Let $T$ be the total number of discrete
time points. Let $M$ be the number of inducing points, and $d$ be the
dimensionality of joint space. The time complexity of _one iteration_ can be
split up into

$\text{Cost}_{\text{iteration}}=C_{\text{sampling}}+C_{\text{likelihood}}+C_{\text{KL}}.$
(10)

We begin by deriving the complexity cost of computing the likelihood for a
single sample, for one single joint configuration at some timestep $t$.

The time complexity of the forward kinematics computation using the Denavit-
Hartenberg (DH) convention and homogeneous transformation matrices is
$\mathcal{O}(d)$, where $d$ is the number of degrees of freedom of the robot
arm. To compute the forward kinematics, we need to multiply the homogeneous
transformation matrices for each joint, as follows:

$\mathbf{F}_{d}=\prod_{i=0}^{d-1}\mathbf{F}_{i,i}\mathbf{F}_{i,i+1}$ (11)

Each $\mathbf{F}_{i}$ is $4\times 4$ matrix computed using the DH parameters,
and hence requires constant amount of computation. Then, an additional
transform is needed to translate the sphere offsets into world coordinates.
This can be easily achieved by considering some other homogeneous transform
$\mathbf{H}_{i,j}$ for some sphere $c_{j},j=\\{1,\cdot\cdot\cdot,s\\}$ and
some frame of reference (in world coordinates) $i$ such that we can get the
sphere placement in world coordinates, namely

$\mathbf{H}_{j}=\mathbf{F}_{i}\mathbf{H}_{i,j}.$ (12)

To do this for all joints and all spheres, the time complexity is
$\mathcal{O}(s)$, which brings the overall cost to $\mathcal{O}(d+s)$. After
the joint configuration has given the robot pose and the sphere position in
world coordinates, the likelihood can be computed. The cost can be broken up
into applying the hinge loss to the signed distances, and the cost of
computing the quadratic form. For the first half, the cost of computing the
signed distance involves a lookup in a 3D grid, incurring constant cost per
sphere. Since the hinge loss’ complexity is linear in the number of spheres,
the overall time complexity of this procedure is $\mathcal{O}(s)$.
Additionally, since the covariance matrix is diagonal, for a given vector of
spheres of size $s$, the quadratic form can be computed in linear time
$\mathcal{O}(s)$. At each iteration of the optimization there are $N_{s}$
samples drawn at $T$ timesteps, the overall cost of computing the likelihood
is $C_{\text{likelihood}}=TN_{s}(D+P)$.

The efficient sampling technique of [29, 30] involves solving a linear system
and computing the Cholesky decomposition of
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$. Therefore, drawing the samples
amounts to $C_{\text{sampling}}=TM^{2}+M^{3}$ for one sample, where we do not
consider the cost of Fourier features, as we view it as a constant
asymptotically.

Lastly, we shift our focus to the computation of the KL term. The dominant
term in the time complexity of computing the KL-divergence is the computation
of the matrix inverse and the determinant of the covariance matrices.
Additionally, we must whiten the covariance
$\mathbf{\Sigma}_{\boldsymbol{u}}$. We have (a) cost of computing $K_{uu}$ is
$\mathcal{O}(dM^{2})$, (b) the cost of computing the Cholesky decomposition
and solving the linear system takes $\mathcal{O}(dM^{3})$ time for a dense
matrix, (c) the cost of whitening the covariance is dominated by solving the
linear system, where the Cholesky factor is computed beforehand and is upper
triangular, which therefore also takes $\mathcal{O}(dM^{3})$. Overall,
$C_{\text{KL}}=dM^{3}$, which gives the overall time complexity
$\mathcal{O}(TM^{2}+dM^{3}+T(d+s))$.

##  Appendix C Experimental Details

Experiments were performed on an AMD Ryzen 7 3700x CPU. For the proposed
planner, we use a Matérn-$\nicefrac{{5}}{{2}}$ kernel, as it is a simple
choice which provides a sufficient degree of differentiability, but is
expressive enough to accommodate different ranges of motion necessary in the
motion plans. Experiments were repeated five times with different seeds. We
use the GPFlow library [17], with the complementary GPFlow-based pathwise
sampling module [29, 30] and PyBullet for simulations [5]. vGPMP uses Adam
with hyperparameters $\beta_{1}=0.8$ and $\beta_{2}=0.95$ and learning rate
$\eta=0.09$. We note that any standard optimizer can be used, as optimization
is carried out on the loss directly. Mean clearance is defined to equal the
negative of the collision term
$-\frac{1}{2}\left\lVert\operatorname{h}_{\varepsilon}(\operatorname{sdf}_{s}(\operatorname{k}_{\operatorname{fwd}}(\sigma(\boldsymbol{f}))))\right\rVert_{\mathbf{\Sigma}_{\operatorname{obs}}}^{2}$
used in the likelihood.

We include RRTConnect, LBKPIECE and CHOMP via the MoveIt! OMPL library [4].
RRTConnect, LBKPIECE were ran with default settings, with post-processing left
on to smoothen the resulting trajectory. CHOMP was ran using the configuration
detailed by [6].

We compare the proposed approach, vGPMP, with GPMP2 [6] in the _lab_
environment, originally introduced by [18], using the (7 degree of freedom)
Barrett WAM robot. We test both planners using 24 different cases with various
start and goal states. We initialize GPMP2 with parameters as outlined in [6]:
$\mathbf{\Sigma}_{\operatorname{obs}}=0.02^{2}\mathbf{I}$, $\varepsilon=0.2$
with 101 temporally equidistant states. Results are detailed in Table 1. For
vGPMP, we set $\mathbf{\Sigma}_{\operatorname{obs}}=0.005^{2}\mathbf{I}$,
$\varepsilon=0.03$, and only 14 inducing locations. For optimization, we use
14 samples per iteration for a total of 150 iterations.

Every parameter used for baselines against RRTCONNECT, LBKPIECE, CHOMP and
GPMP2 are outlined in Table 3. The number of iterations for _Franka_ and
_UR10_ is 130. As seen in Table 2, in the _industrial_ environment using UR10,
our approach finds shorter paths in general, which is in part due to higher
length scale initialization than usual.

We evaluated the performance of vGPMP and other baseline planners in two
different environments, which can be seen in Figure 6, using various start and
goal states. vGPMP generates smooth trajectories that are also further away
from the obstacles in comparison to baselines as illustrated in Figure 7. To
ensure balanced comparisons, average success rate reported was calculated as a
weighted average over the number of problems.

Robot | Env |  | Parameters
---|---|---|---
$\mathbb{E}[q(\textbf{u})]$ | length scale | kernel variance | $\sigma_{obs}$ | Ind. var. | $\gamma$
Franka | I | Zeros | $[2,2,2,2,2,2,2]$ | 0.2 | 0.005 | 10 | 100
Train | Train | Train | Not Train | Train | Not Train
B | Interp. | $[4,4,4,3,3,3,3]$ | 0.5 | 0.0005 | 24 | 100
Train | Train | Train | Not Train | Train | Not Train
UR10 | I | Interp. | $[6,6,6,6,6,6]$ | 0.1 | 0.0001 | 18 | 100
Train | Train | Train | Not Train | Train | Not Train
B | Upright | $[4,4,4,4,4,4]$ | 0.25 | 0.0005 | 12 | 100
Train | Train | Train | Not Train | Not Train | Not Train
WAM | Lab | Interp. | $[5,5,5,4,5,5,5]$ | 0.05 | 0.005 | 14 | 100
Train | Train | Train | Not Train | Not Train | Not Train
Table 3: Experiment configurations for various robots in different
environments. The environments (Env) used are named _industrial_ (I) and
_bookshelves_ (B). Interp. stands for _interpolated_. In all cases, set
$\varepsilon=0.05$, aside from UR10 on the _industrial_ environment where it
was $\varepsilon=0.08$ and WAM on the _lab_ environment where it was
$\varepsilon=0.03$.

(a) Industrial environment

(b) Bookshelves environment

Figure 6: Example distributions of collision-free paths generated by vGPMP in
two benchmark environments. Sampled trajectories are shown in green.  Figure
7: Comparison of trajectories from a baseline motion planner, RRTconnect (in
blue) and vGPMP (in green). As in this example, in general vGPMP leads to
smooth trajectories further away from the obstacles. Spheres around robot
links that are used in collision detection can also be seen.


aainstitutetext: Center for Computational Sciences, University of Tsukuba,
Tsukuba, Ibaraki 305-8577, Japanbbinstitutetext: Graduate School of Science,
The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan

# Tensor renormalization group study of (1+1)-dimensional U(1) gauge-Higgs
model at $\theta=\pi$ with Lüscher’s admissibility condition

Shinichiro Akiyama<EMAIL_ADDRESS>a Yoshinobu Kuramashi
<EMAIL_ADDRESS>
###### Abstract

We investigate the phase structure of the (1+1)-dimensional U(1) gauge-Higgs
model with a $\theta$ term, where the U(1) gauge action is constructed with
Lüscher’s admissibility condition. Using the tensor renormalization group,
both the complex action problem and topological freezing problem in the
standard Monte Carlo simulation are avoided. We find the first-order phase
transition with sufficiently large Higgs mass at $\theta=\pi$, where the
$\mathds{Z}_{2}$ charge conjugation symmetry is spontaneously broken. On the
other hand, the symmetry is restored with a sufficiently small mass. We
determine the critical endpoint as a function of the Higgs mass parameter and
show the critical behavior is in the two-dimensional Ising universality class.

††preprint: UTHEP-789, UTCCS-P-157

##  1 Introduction

At the end of the last century, Lüscher introduced an admissibility condition
for the gauge fields to be separated into disconnected subspaces corresponding
to topological charges in the continuum theory Luscher:1998du :

$\displaystyle\|1-P_{\mu\nu}(n)\|<\epsilon\;\;{\rm
for}\;\;^{\forall}n,\mu,\nu,$ (1)

where $\epsilon$ is a positive constant and $P_{\mu\nu}(n)$ is a product of
link variables $U_{\mu}(n)$ as in the standard way,

$\displaystyle
P_{\mu\nu}(n)=U_{\mu}(n)U_{\nu}(n+{\hat{\mu}})U_{\mu}^{\dagger}(n+{\hat{\nu}})U_{\nu}^{\dagger}(n).$
(2)

The link variable $U_{\mu}(n)$ lives on the link connecting the sites $n$ and
$n+\hat{\mu}$. As an example, he proposed the following gauge action to make
the link variables satisfy the above condition:

$\displaystyle\beta
S_{g}=\begin{cases}\displaystyle\beta\sum_{n,\mu>\nu}\frac{1-{\rm
Re}P_{\mu\nu}(n)}{1-\|1-P_{\mu\nu}(n)\|/\epsilon}\;\;&{\rm
if}\;\;\|1-P_{\mu\nu}(n)\|<\epsilon,\\\ \infty&{\rm otherwise},\end{cases}$
(3)

with the inverse gauge coupling $\beta$. This action should have an advantage
in investigating the topological effects of the gauge theories. However, early
numerical studies with this action revealed that the topological change is
substantially suppressed in the Monte Carlo histories Fukaya:2003ph  and it is
difficult to evaluate the contributions from different topological sectors. As
long as the Monte Carlo method is employed, the possible way is to perform
calculations in the fixed topological sectors Fukaya:2003ph ; Fukaya:2004kp ;
Bietenholz:2005rd ; Bietenholz:2016ymo  or to utilize open boundary
conditions, dismissing the translational invariance of the system
Luscher:2011kk .

However, this problem is potentially solved by the tensor renormalization
group (TRG) method. 111 In this paper, the “TRG method” or the “TRG approach”
refers to not only the original numerical algorithm proposed by Levin and Nave
Levin:2006jai  but also its extensions Gu:2010yh ; PhysRevB.86.045139 ;
Shimizu:2014uva ; Sakai:2017jwp ; Adachi:2019paf ; Kadoh:2019kqk ;
Akiyama:2020soe .  The major advantages of the TRG method over the Monte Carlo
simulation are (i) no sign problem Denbleyker:2013bea ; Shimizu:2014fsa ;
Kawauchi:2016xng ; Shimizu:2017onf ; Unmuth-Yockey:2018ugm ; Kadoh:2019ube ;
Kuramashi:2019cgs ; Butt:2019uul ; Bloch:2021uup ; Takeda:2021mnc ;
Bloch:2021mjw ; Nakayama:2021iyp ; Hirasawa:2021qvh , (ii) logarithmic
computational cost on the system size, (iii) direct manipulation of the
Grassmann variables Gu:2010yh ; Shimizu:2014uva ; Takeda:2014vwa ;
Sakai:2017jwp ; Yoshimura:2017jpk ; Kadoh:2018hqq ; Akiyama:2021xxr ;
Bloch:2022vqz ; Akiyama:2023lvr ; Yosprakob:2023tyr ; Asaduzzaman:2023pyz ,
and (iv) evaluation of the partition function or the path integral itself. The
advantage (iv) ensures that the TRG calculation automatically includes full
contributions from different topological sectors. Moreover, the TRG method
assumes the translational invariance of the system and can easily impose
periodic boundary conditions.

In this paper, we investigate the phase structure of the (1+1)-dimensional
((1+1)$d$) U(1) gauge-Higgs model with a $\theta$ term, where the topological
effects play an essential role, employing Lüscher’s gauge action of Eq. (3).
The Monte Carlo simulation of this model is extremely difficult due to a
double whammy of the complex action problem and the topological freezing.
Figure 1 illustrates the expected phase diagram Komargodski:2017dmc . The
model exhibits the first-order phase transition at $\theta=\pi$, where the
$\mathds{Z}_{2}$ charge conjugation symmetry is spontaneously broken in the
large positive Higgs mass-squared regime, including the pure gauge limit. 222
There are several earlier studies on the $2d$ pure gauge theory with a
$\theta$ term by the density of state approach Gattringer:2020mbf , complex
Langevin method Hirasawa:2020bnl , and TRG Kuramashi:2019cgs ;
Hirasawa:2021qvh .  Once the Higgs mass-squared is sufficiently reduced, the
symmetry is restored. We determine the critical endpoint as a function of the
Higgs mass-squared and show its critical behavior is in the 2$d$ Ising
universality class based on the numerical analysis of the transfer matrix and
topological charge density. We also compare our results with the previous work
employing the dual lattice simulation based on the Villain gauge action, which
is a non-compact gauge action on the lattice Gattringer:2018dlw .

This paper is organized as follows. In Sec. 2, we define the U(1) gauge-Higgs
model with a $\theta$ term on a (1+1)$d$ lattice. We also demonstrate how to
represent the path integral as a tensor network. In Sec. 3, we first present
the results for the pure U(1) gauge action which corresponds to the infinitely
heavy limit of the Higgs mass, where the first-order phase transition takes
place at $\theta=\pi$. After that, we discuss the phase transition with the
finite lattice Higgs mass and determine the critical endpoint and its
universality class. Section 4 is devoted to summary and outlook.

Figure 1:  Schematic phase diagram of (1+1)$d$ U(1) gauge-Higgs model with a
$\theta$ term. The horizontal axis denotes the Higgs mass-squared. The red
line denotes the first-order phase transition, which terminates at the
critical endpoint expressed by the red blob.

##  2 Tensor network formulation

###  2.1 The model on a lattice

The U(1) gauge-Higgs model with a $\theta$ term is defined by

$\displaystyle S=\beta S_{g}+S_{h}+S_{\theta}.$ (4)

We always consider the model on a square lattice with periodic boundary
conditions. We employ the Lüscher gauge action,

$\displaystyle S_{g}=\begin{cases}\displaystyle\sum_{n}\frac{1-{\rm
Re}P_{12}(n)}{1-[1-{\rm Re}P_{12}(n)]/\epsilon}&{\rm if~{}``admissible"}\\\
\infty&{\rm otherwise}\end{cases},$ (5)

where $P_{12}(n)$ is defined by Eq. (2) with $U_{\nu}(n)={\rm
e}^{\mathrm{i}\vartheta_{\nu}(n)}$ and $\vartheta_{\nu}(n)\in[-\pi,\pi]$. The
admissibility condition is given by

$\displaystyle 1-{\rm Re}P_{12}(n)<\epsilon.$ (6)

When this condition is satisfied, the corresponding gauge fields are called
admissible. The admissibility condition makes gauge fields smooth and
unphysical configurations are suppressed. The space of admissible gauge fields
is separated into disconnected subspaces which are labeled by the integers
corresponding to topological charges in the continuum Luscher:1998du . The
Higgs part is defined by

$\displaystyle S_{h}$
$\displaystyle=-\sum_{n}\sum_{\nu}\left[\phi^{*}(n)U_{\nu}(n)\phi(n+\hat{\nu})+\phi^{*}(n+\hat{\nu})U^{*}_{\nu}(n)\phi(n)\right]$
$\displaystyle+M\sum_{n}\left|\phi(n)\right|^{2}+\lambda\sum_{n}\left|\phi(n)\right|^{4}.$
(7)

The complex-valued Higgs fields are denoted by $\phi(n)$ and $M=m^{2}+4$ is
the lattice Higgs mass where $m$ corresponds to the Higgs mass parameter in
the continuum action. The quartic coupling constant is denoted by $\lambda$.
Finally, the $\theta$ term is defined by

$\displaystyle S_{\theta}=-\frac{\mathrm{i}\theta}{2\pi}\sum_{n}\ln
P_{12}(n).$ (8)

###  2.2 Tensor network formulation

We consider the tensor network representation of the path integral defined as

$\displaystyle
Z=\prod_{n,\nu}\int^{\pi}_{-\pi}\frac{\mathrm{d}\vartheta_{\nu}(n)}{2\pi}\prod_{n}\int_{\mathds{C}}\frac{\mathrm{d}\phi(n)}{2\pi}\exp(-S).$
(9)

Parametrizing the complex-valued Higgs field by
$\phi(n)=r(n)\mathrm{e}^{\mathrm{i}\varphi(n)}$, the integral measure is
represented as

$\displaystyle\int_{\mathds{C}}\frac{\mathrm{d}\phi(n)}{2\pi}=\int^{\infty}_{0}r(n)\mathrm{d}r(n)\int^{\pi}_{-\pi}\frac{\mathrm{d}\varphi(n)}{2\pi},$
(10)

and $S_{h}$ reads

$\displaystyle S_{h}$
$\displaystyle=-\sum_{n}\sum_{\nu}2r(n)r(n+\hat{\nu})\cos\left[\varphi(n+\hat{\nu})-\varphi(n)+\vartheta_{\nu}(n)\right]+\sum_{n}\left[Mr(n)^{2}+\lambda
r(n)^{4}\right].$ (11)

We further introduce $\ell(n)=r(n)^{2}$ and rewrite Eq. (11) as

$\displaystyle S_{h}$
$\displaystyle=-\sum_{n}\sum_{\nu}2\sqrt{\ell(n)\ell(n+\hat{\nu})}\cos\left[\varphi(n+\hat{\nu})-\varphi(n)+\vartheta_{\nu}(n)\right]+\sum_{n}\left[M\ell(n)+\lambda\ell(n)^{2}\right].$
(12)

Using the invariance of the Haar measure, or choosing the unitary gauge, we
can eliminate $\varphi(n)$ from the path integral,

$\displaystyle
Z=\prod_{n,\nu}\int^{\pi}_{-\pi}\frac{\mathrm{d}\vartheta_{\nu}(n)}{2\pi}\prod_{n}\int^{\infty}_{0}\frac{\mathrm{d}\ell(n)}{2}\exp\left[-\beta
S_{g}-S^{\prime}_{h}-S_{\theta}\right],$ (13)

with

$\displaystyle S^{\prime}_{h}$
$\displaystyle=-\sum_{n}\sum_{\nu}2\sqrt{\ell(n)\ell(n+\hat{\nu})}\cos\vartheta_{\nu}(n)+\sum_{n}\left[M\ell(n)+\lambda\ell(n)^{2}\right].$
(14)

In this study, we use the Gauss-Laguerre quadrature rule to discretize the
integral over $\ell(n)$ and the Gauss-Legendre one for $\vartheta_{\nu}(n)$.
The efficacy of these Gauss quadrature rules has been reported in the previous
TRG studies for the U(1) pure gauge theory with a $\theta$ term
Kuramashi:2019cgs  and complex $\phi^{4}$ theories Kadoh:2019ube ;
Akiyama:2020ntf . The path integral is now approximated by

$\displaystyle Z$ $\displaystyle\simeq Z(K_{g},K_{h})$
$\displaystyle=\prod_{n,\nu}\sum_{\tilde{\vartheta}_{\nu}(n)\in
D_{g}}\frac{w_{\tilde{\vartheta}_{\nu}(n)}}{2}\prod_{n}\sum_{\tilde{\ell}(n)\in
D_{h}}\frac{w_{\tilde{\ell}(n)}\mathrm{e}^{\tilde{\ell}(n)}}{2}\exp\left[-\beta\tilde{S}_{g}-\tilde{S}^{\prime}_{h}-\tilde{S}_{\theta}\right],$
(15)

where

$\displaystyle\tilde{S}_{g}=\begin{cases}\displaystyle\sum_{n}\frac{1-\cos\pi\left(\tilde{\vartheta}_{1}(n)+\tilde{\vartheta}_{2}(n+\hat{1})-\tilde{\vartheta}_{1}(n+\hat{2})-\tilde{\vartheta}_{2}(n)\right)}{1-\left[1-\cos\pi\left(\tilde{\vartheta}_{1}(n)+\tilde{\vartheta}_{2}(n+\hat{1})-\tilde{\vartheta}_{1}(n+\hat{2})-\tilde{\vartheta}_{2}(n)\right)\right]/\epsilon}&{\rm
if~{}admissible}\\\ \infty&{\rm otherwise}\end{cases},$ (16)
$\displaystyle\tilde{S}^{\prime}_{h}$
$\displaystyle=-\sum_{n}\sum_{\nu}2\sqrt{\tilde{\ell}(n)\tilde{\ell}(n+\hat{\nu})}\cos\pi\tilde{\vartheta}_{\nu}(n)+\sum_{n}\left[M\tilde{\ell}(n)+\lambda\tilde{\ell}(n)^{2}\right],$
(17)
$\displaystyle\tilde{S}_{\theta}=-\frac{\mathrm{i}\theta}{2\pi}\sum_{n}\ln\left[{\rm
e}^{\mathrm{i}\pi\left(\tilde{\vartheta}_{1}(n)+\tilde{\vartheta}_{2}(n+\hat{1})-\tilde{\vartheta}_{1}(n+\hat{2})-\tilde{\vartheta}_{2}(n)\right)}\right].$
(18)

In Eq. (2.2), $\tilde{\ell}(n)$ denotes the sampling point according to the
Gauss-Laguerre quadrature and $w_{\tilde{\ell}(n)}$ is the corresponding
weight. The number of sampling points in $D_{h}$ is denoted by $K_{h}$.
Similarly, $\tilde{\vartheta}_{\nu}(n)$ denotes the sampling point according
to the Gauss-Legendre quadrature and $w_{\tilde{\vartheta}_{\nu}(n)}$ is the
corresponding weight. The number of sampling points in $D_{g}$ is denoted by
$K_{g}$. In the limits of $K_{g}\to\infty$ and $K_{h}\to\infty$, the original
path integral is restored. Eq. (2.2) is ready to be described as a tensor
network. We introduce four-leg pure gauge tensors on each plaquette as
follows,

$\displaystyle T^{(g)}_{x_{g}y_{g}x^{\prime}_{g}y^{\prime}_{g}}$
$\displaystyle=\begin{cases}\displaystyle\frac{\sqrt{w_{x_{g}}w_{y_{g}}w_{x^{\prime}_{g}}w_{y^{\prime}_{g}}}}{2^{2}}\exp\left[-\beta\frac{1-\cos\pi\left(y^{\prime}_{g}+x_{g}-y_{g}-x^{\prime}_{g}\right)}{1-\left[1-\cos\pi\left(y^{\prime}_{g}+x_{g}-y_{g}-x^{\prime}_{g}\right)\right]/\epsilon}\right]&{\rm
if~{}admissible}\\\ 0&{\rm otherwise}\end{cases},$ (19) $\displaystyle
T^{(\theta)}_{x_{g}y_{g}x^{\prime}_{g}y^{\prime}_{g}}=\exp\left(\frac{\mathrm{i}\theta}{2\pi}\ln\left[{\rm
e}^{\mathrm{i}\pi\left(y^{\prime}_{g}+x_{g}-y_{g}-x^{\prime}_{g}\right)}\right]\right).$
(20)

For the Higgs part, we introduce the following hopping matrix,

$\displaystyle
H_{\tilde{\ell}(n)\tilde{\theta}_{\nu}(n)\tilde{\ell}(n+\hat{\nu})}$
$\displaystyle=\frac{\sqrt[4]{w_{\tilde{\ell}(n)}w_{\tilde{\ell}(n+\hat{\nu})}}\mathrm{e}^{(\tilde{\ell}(n)+\tilde{\ell}(n+\hat{\nu}))/4}}{\sqrt{2}}$
$\displaystyle\times\exp\left[2\sqrt{\tilde{\ell}(n)\tilde{\ell}(n+\hat{\nu})}\cos\pi\tilde{\theta}_{\nu}(n)-\frac{M}{4}\left(\tilde{\ell}(n)+\tilde{\ell}(n+\hat{\nu})\right)-\frac{\lambda}{4}\left(\tilde{\ell}(n)^{2}+\tilde{\ell}(n+\hat{\nu})^{2}\right)\right].$
(21)

Now, we perform the singular value decomposition (SVD) of the
$\nu$-directional hopping matrix, which gives us

$\displaystyle
H_{\tilde{\ell}(n)\tilde{\vartheta}_{\nu}(n)\tilde{\ell}(n+\hat{\nu})}\simeq\sum_{\alpha=1}^{\chi}A_{\tilde{\ell}(n)\tilde{\vartheta}_{\nu}(n)\alpha}B_{\tilde{\ell}(n+\hat{\nu})\alpha},$
(22)

where $A$ and $B$ are defined by unitary matrices multiplied by the square
root of singular values $\sigma_{\alpha}$. In this study, we choose $\chi$ in
Eq. (22) such that the singular values satisfying
$\sigma_{\alpha}/\sigma_{1}>10^{-7}$ are kept. Note that $\sigma_{1}$ is the
largest singular value and $\sigma_{\alpha}$ is in the descending order. We
are now allowed to integrate out $\tilde{\ell}(n)$ at each site $n$. As a
result, we can define a six-leg tensor at each lattice site as,

$\displaystyle
T^{(h)}_{x_{h}y_{h}x^{\prime}_{g}x^{\prime}_{h}y^{\prime}_{g}y^{\prime}_{h}}=\sum_{\tilde{\ell}(n)}A_{\tilde{\ell}(n)y^{\prime}_{g}x_{h}}A_{\tilde{\ell}(n)x^{\prime}_{g}y_{h}}B_{\tilde{\ell}(n)x^{\prime}_{h}}B_{\tilde{\ell}(n)y^{\prime}_{h}}.$
(23)

Therefore, the tensor network representation for Eq. (2.2) is obtained as

$\displaystyle Z(K_{g},K_{h})={\rm tTr}\left[\prod_{n}T_{n}\right],$ (24)

with the fundamental tensor $T_{n}$ at each site $n$,

$\displaystyle(T_{n})_{xyx^{\prime}y^{\prime}}=T^{(g)}_{x_{g}y_{g}x^{\prime}_{g}y^{\prime}_{g}}T^{(\theta)}_{x_{g}y_{g}x^{\prime}_{g}y^{\prime}_{g}}T^{(h)}_{x_{h}y_{h}x^{\prime}_{g}x^{\prime}_{h}y^{\prime}_{g}y^{\prime}_{h}},$
(25)

whose bond dimension is $K_{g}\chi$. Note that the subscripts in the left-hand
side of Eq. (25) are defined as $i=(i_{g}i_{h})$ with
$i=x,y,x^{\prime},y^{\prime}$.

###  2.3 Coarse-graining algorithm

We apply the bond-weighted TRG (BTRG) algorithm PhysRevB.105.L060402  to
approximately compute the path integral in Eq. (24). BTRG allows us to
approximately carry out the contractions among $2^{q}$ local tensors within
the $q$ times of coarse-graining. Since each local tensor in Eq. (25) is
defined on each lattice site, $q$ relates to the volume $V$ via $V=2^{q}$ and
the linear system size $L$ via $L=2^{q/2}$.

This algorithm improves the accuracy of the original Levin-Nave TRG at the
same bond dimension. Remarkably, the computational cost of BTRG is completely
the same as the Levin-Nave TRG. The essence of the BTRG is to introduce a
weight on each edge of the tensor network. These weights mimic the effect of
the environment, which is not taken into account in the original Levin-Nave
TRG. Therefore, the BTRG can be regarded as a variant of the second
renormalization group algorithms Xie:2009zzd ; PhysRevB.81.174411 . For the
algorithmic details, see Ref. PhysRevB.105.L060402 . 333 Since the model is
defined on a square lattice, we always set the hyperparameter in the BTRG
algorithm as $k=-0.5$ PhysRevB.105.L060402 ; Akiyama:2022pse .

##  3 Numerical results

In the following, we always set $\beta=3.0$ and the positive constant
$\epsilon$ in Eq. (6) as $\epsilon=1.0$. For the gauge-Higgs model, the
quartic coupling is fixed as $\lambda=0.5$. Note that the cutoff effects from
the finite lattice spacing of the Lüscher gauge action and standard Wilson
gauge action are different even at the same inverse gauge coupling $\beta$.
With the same value of $\beta$, the Lüscher gauge action is expected to be
closer to the continuum limit than the Wilson action. See Appendix A for the
comparison between the Lüscher and Wilson gauge actions.

###  3.1 Pure U(1) gauge theory

We start by studying the (1+1)$d$ pure gauge theory with a $\theta$ term to
validate our tensor network formulation. In this case, the local tensor in Eq.
(25) is defined without $T^{(h)}$ and the bond dimension in Eq. (24) is equal
to $K_{g}$.

At $\theta=\pi$, the theory is expected to undergo the first-order transition.
We calculate the topological charge density $\langle Q\rangle/V$, which is
defined by

$\displaystyle\frac{\langle
Q\rangle}{V}=-\frac{\mathrm{i}}{V}\frac{\partial\ln Z}{\partial\theta}.$ (26)

Figure 2 shows the volume dependence of $\langle Q\rangle/V$ at $\beta=3.0$ as
a function of $\theta$. We set $D=K_{g}=30$ and employed the numerical
differentiation to evaluate Eq. (26). As the volume increases, the topological
charge density gradually becomes discontinuous and we observe a clear jump at
$\theta=\pi$ with $L\geq 2^{13}$. This is a signal of the first-order
transition at $\theta=\pi$.

Figure 2:  Topological charge density for the pure U(1) gauge theory as a
function of $\theta/\pi$ at $\beta=3.0$ with $D=K_{g}=30$.

Taking advantage of the TRG method, we also compute the ground state
degeneracy by introducing

$\displaystyle X=\frac{\left({\rm Tr}A\right)^{2}}{{\rm
Tr}\left(A^{2}\right)},$ (27)

following Ref. PhysRevB.80.155131 . After sufficient times of coarse-graining,
this quantity counts the ground state degeneracy. In Eq. (27), $A$ is a
transfer matrix defined from the local tensor via

$\displaystyle A_{yy^{\prime}}=\sum_{x}T_{xyxy^{\prime}}.$ (28)

Figure 3 shows $X$ in Eq. (27) as a function of $\theta$. With sufficiently
large volume, $X=2$ is realized only at $\theta=\pi$. Therefore, the current
TRG computation successfully reproduces the spontaneous $\mathds{Z}_{2}$
symmetry breaking at $\theta=\pi$ as expected.

Figure 3:  The ground state degeneracy $X$ of Eq. (27) for the pure U(1) gauge
theory as a function of $\theta/\pi$ at $\beta=3.0$ with $D=K_{g}=30$.  Figure
4:  System size dependence of $\chi_{\rm peak}$ for the pure U(1) gauge theory
at $\beta=3.0$ with $D=K_{g}=30$.

Additionally, we apply the finite-size scaling analysis for the topological
susceptibility, which is given by

$\displaystyle\chi_{Q}=\frac{\partial}{\partial\theta}\frac{\langle
Q\rangle}{V}.$ (29)

We use the numerical differentiation with the $O(\Delta^{4})$ accuracy,
setting $\Delta=\pi\times 10^{-5}$ around $\theta=\pi$, to evaluate Eq. (29).
Since the peak of the topological susceptibility always appears at
$\theta=\pi$, we examine the system size dependence of the values of
$\chi_{Q}$ at $\theta=\pi$. Let $\chi_{\rm peak}(L)$ be the value of
$\chi_{Q}$ at $\theta=\pi$ with the linear system size $L$. Using the
following ansatz,

$\displaystyle\chi_{\rm peak}(L)=c_{0}+c_{1}L^{p},$ (30)

the data of $\chi_{\rm peak}(L)$ are fitted as shown in Figure 4. The fit is
performed on the data over $L\in[4,256\sqrt{2}]$ and we obtain $p=2.00001(6)$,
$c_{0}=-0.0196(6)$, and $c_{1}=0.0009944(3)$. Therefore, $\chi_{\rm peak}(L)$
is proportional to the volume and this is another indication of the first-
order phase transition. Although we employed the finite bond dimension $D$ and
cutoff $K_{g}$, the current computation correctly reproduces the first-order
phase transition at $\theta=\pi$. We have also tried the same analysis with
$D=K_{g}=20$, which results in $p=2.00000(5)$, $c_{0}=-0.0197(5)$, and
$c_{1}=0.0009973(3)$. Therefore, $D,K_{g}\geq 20$ seems sufficiently large to
investigate the phase transition at $\theta=\pi$ for the pure U(1) gauge
theory at $\beta=3.0$.

###  3.2 The U(1) gauge-Higgs model

We now investigate the critical behavior of the U(1) gauge-Higgs model at
$\theta=\pi$. Figure 5 shows the thermodynamic topological charge density at
$M=2.99$ and $M=3.00$, where we set $K_{g}=20$, $K_{h}=20$ which result in
$\chi=8$ according to Eq. (22). We also set $D=K_{g}\chi=160$. As we will see,
our choice of $D=K_{g}\chi=160$ is sufficiently large to determine the
universality class at the critical endpoint. The behavior of topological
charge density is highly different between these two mass parameters; $\langle
Q\rangle/V$ varies continuously around $\theta=\pi$ at $M=2.99$ but the
discontinuity appears at $\theta=\pi$ when $M=3.00$. Therefore, the critical
mass $M_{\rm c}$ should exist in the range of $2.99\leq M_{\rm c}\leq 3.00$.

Figure 5:  Topological charge density in the thermodynamic limit for the U(1)
gauge-Higgs model as a function of $\theta/\pi$ at $\beta=3.0$ and
$\lambda=0.5$ with $K_{g}=20$, $K_{h}=20$, and $D=160$. The circle and diamond
denote $M=2.99$ and $M=3.00$, respectively.  Figure 6:  The ground state
degeneracy $X$ computed on $V=2^{40}$ as a function of $M$.

To locate the critical point $M_{\rm c}$ more precisely, we use the ground
state degeneracy $X$ in Eq. (27). When $M>M_{\rm c}$, we should have $X=2$ due
to the $\mathds{Z}_{2}$ charge symmetry breaking and with $M<M_{\rm c}$, we
must have $X=1$. Figure 6 shows the degeneracy $X$ as a function of $M$, which
results in $2.99747\leq M_{\rm c}\leq 2.99748$.

To identify the universality class, we analyze the eigenvalue spectrum of the
transfer matrix $A$ defined in Eq. (28). Note that the analysis on the
transfer matrix obtained from the TRG has been employed in investigating
several spin models PhysRevB.104.165132 ; PhysRevB.80.155131 ; Huang:2023vwp ;
PhysRevB.108.024413 ; Az-zahra:2024gqr  and the lattice Schwinger model
Shimizu:2017onf . For a comprehensive overview, see Ref.
ueda2024renormalizationgroupflowfixedpoint . Although we are dealing with
continuous gauge and scalar degrees of freedom, a similar analysis should be
possible if our discretization scheme based on the Gauss quadrature rules is
working. Particularly, we employ the tensor-network-based level spectroscopy
proposed in Refs. PhysRevB.104.165132 ; PhysRevB.108.024413 . Let
$\lambda_{n}~{}(n=0,1,\cdots)$ be an $n$-th eigenvalue of the transfer matrix
$A$. Assuming that these eigenvalues are in descending order, we can obtain
the scaling dimension $x_{n}(L)$ at the finite system size $L$ via

$\displaystyle
x_{n}(L)=\frac{1}{2\pi}\ln\frac{\lambda_{0}(L)}{\lambda_{n}(L)}.$ (31)

When these quantities are computed at the criticality and on a sufficiently
large volume, they give us the scaling dimension of the conformal field theory
(CFT). Figure 7 shows $x_{1}(L)$ as a function of the Higgs mass $M$. The
volume independence can be observed at $M\sim 2.99748$ with $x_{1}=1/8$, which
is consistent with the 2$d$ Ising universality class. From now on, we assume
the critical phenomenon is in the 2$d$ Ising universality class. Following
Ref. PhysRevB.108.024413 , we consider a combined scaling dimension,

$\displaystyle x_{\rm cmb}(L)=x_{1}(L)+\frac{1}{16}x_{2}(L),$ (32)

which removes the effect of the leading irrelevant perturbation associated
with the scaling dimension 4. Since $x_{2}=1$ in the 2$d$ Ising universality
class, $x_{\rm cmb}=3/16$ is expected at the criticality. Figure 8 shows
$x_{\rm cmb}(L)$ against the Higgs mass $M$ which is in agreement with this
expectation. Note that $x_{\rm cmb}(L)$ shown in Figure 8 also agrees with the
previous estimation by $X$: $2.99747\leq M_{\rm c}\leq 2.99748$. Moreover, we
apply a scheme to determine the critical point following Ref.
PhysRevB.108.024413 , which assumes the phase transition is in the $2d$ Ising
universality class. We first choose two mass parameters $M^{(+)}$ and
$M^{(-)}$ such that $M^{(-)}\leq M_{\rm c}\leq M^{(+)}$ and compute $\delta
x_{\rm cmb}(L)=x_{\rm cmb}(L)-3/16$. When $M\leq M_{\rm c}$ ($M\geq M_{\rm
c}$), $\delta x_{\rm cmb}(L)$ should increase (decrease) by enlarging the
system size as shown in Fig. 8. Secondly, we perform linear interpolations of
$\delta x_{\rm cmb}(L)$ between $M^{(+)}$ and $M^{(-)}$, and find a crossing
point $M^{*}(L)$ of two lines with the system size $L$ and $\sqrt{2}L$. Then,
$M_{\rm c}$ is obtained via $M^{*}(L)=M_{\rm c}+aL$, where $M_{\rm c}$ and $a$
are the parameters determined by the numerical fit. The fit results in $M_{\rm
c}=2.997480(2)$, where we set $M^{(-)}=2.99747$, $M^{(+)}=2.99749$ using the
data with $L\in[2^{10},2^{14}\sqrt{2}]$. Therefore, the tensor-network-based
level spectroscopy confirms the critical point estimated by the ground state
degeneracy. We also investigate the finite-size correction for the free energy
in the thermodynamic limit via

$\displaystyle\frac{1}{L^{2}}\ln\lambda_{0}=-f_{\infty}+\frac{\pi c}{6L^{2}},$
(33)

where $f_{\infty}$ is the thermodynamic free energy and $c$ is the central
charge of the CFT at criticality. Using Eq. (33) as a fitting ansatz for
$\ln\lambda_{0}/L^{2}$, we can determine the central charge. As a
representative point, we choose $M=2.99748$ and the fit using the data with
$L\in[2^{10},2^{15}]$ results in $c=0.50(7)$, which is consistent with
$c=1/2$, as expected.

Figure 7:  The system-size dependence of the scaling dimension $x_{1}(L)$. The
dashed line denotes the theoretical value of the $2d$ Ising universality
class, $x_{1}=1/8$.  Figure 8:  The system-size dependence of the scaling
dimension $x_{\rm cmb}(L)$. The dashed line denotes the theoretical value of
the $2d$ Ising universality class, $x_{\rm cmb}=3/16$.

We finally remark on the location of $M_{\rm c}$ which slightly depends on the
algorithmic parameters, $K_{g}$, $K_{h}$, and $D=K_{g}\chi$. Table 1 shows the
critical endpoint $M_{\rm c}$ estimated by the scheme in Ref.
PhysRevB.108.024413 . In Table 1, the error originates from the fit based on
$M^{*}(L)=M_{\rm c}+aL$. We find that these estimations are comparable with
$M_{\rm c}=2.989(2)$ obtained in Ref. Gattringer:2018dlw , where the gauge-
field Boltzmann weight is given by the Villain form and the Monte Carlo
simulation is performed based on the dual representation. Note that $M_{\rm
c}=2.989(2)$ in Ref. Gattringer:2018dlw  is obtained at $\beta=3.0$ and
$\lambda=0.5$.

Table 1:  Comparison of the critical endpoint $M_{\rm c}$ against the algorithmic parameters.  $K_{g}$ | $K_{h}$ | $\chi$ | $D$ | $M_{\rm c}$
---|---|---|---|---
24 | 20 | 8 | 192 | 2.9982886(1)
22 | 20 | 8 | 176 | 2.9998263(13)
20 | 20 | 8 | 160 | 2.9974765(14)
24 | 10 | 6 | 144 | 2.9929635(1)
22 | 10 | 6 | 132 | 2.9945222(9)
20 | 10 | 7 | 140 | 2.9921698(6)

##  4 Summary and outlook

We have analyzed the phase structure of the (1+1)$d$ U(1) gauge-Higgs model
with a $\theta$ term, whose gauge action is constructed with Lüscher’s
admissibility condition. Although the model suffers from a complex action
problem and topological freezing within the Monte Carlo simulation, the TRG
approach allows us to deal with the model successfully. We have observed the
first-order phase transition at $\theta=\pi$ with sufficiently large lattice
Higgs mass $M$, including the pure gauge theory, with the finite gauge
coupling $\beta$ and the quartic coupling $\lambda$. We have determined the
critical endpoint and its universality class from the numerical analysis of
the transfer matrix, which can be directly obtained from the TRG computation.
Employing the tensor-network-based level spectroscopy, we have confirmed that
the scaling dimensions are consistent with the 2$d$ Ising universality class.

All these results show that the TRG is a promising approach to deal with the
lattice gauge theory with Lüscher’s admissibility condition. We emphasize that
one can easily combine the Lüscher gauge action with fermions and extend the
theory for the higher dimensions within the tensor network formulation. It
would be an interesting future work to investigate the Schwinger model with a
$\theta$ term under Lüscher’s admissibility condition.

###### Acknowledgements.

Numerical calculation for the present work was carried out with SQUID at the
Cybermedia Center, Osaka University (Project ID: hp240012). We also used the
supercomputer Fugaku provided by RIKEN through the HPCI System Research
Project (Project ID: hp230247) and computational resources of Wisteria/BDEC-01
and Cygnus and Pegasus under the Multidisciplinary Cooperative Research
Program of Center for Computational Sciences, University of Tsukuba. This work
is supported by Grants-in-Aid for Scientific Research from the Ministry of
Education, Culture, Sports, Science and Technology (MEXT) (Nos. 24H00214,
24H00940). SA acknowledges the support from the Endowed Project for Quantum
Software Research and Education, the University of Tokyo qsw , JSPS KAKENHI
Grant Number JP23K13096, and the Center of Innovations for Sustainable Quantum
AI (JST Grant Number JPMJPF2221).

##  Appendix A Comparison with the Wilson gauge action

Employing the pure U(1) gauge theory with a $\theta$ term, we make a
comparison between the Lüscher gauge action in Eq. (5) and the Wilson gauge
action,

$\displaystyle S_{g}=\sum_{n}\left(1-{\rm Re}P_{12}(n)\right).$ (34)

Here, we consider a $\theta$ term defined by

$\displaystyle S_{\theta}=\frac{\mathrm{i}\theta}{2\pi}\sum_{n}{\rm
Im}P_{12}(n),$ (35)

instead of Eq. (8). In the continuum limit, Eq. (35) reproduces

$\displaystyle S_{\theta}=\frac{\mathrm{i}\theta}{2\pi}\int{\rm
d}^{2}xF_{12}(x),$ (36)

which is a topological term of the continuum theory and $F_{12}(x)$ is the
field strength in two dimensions. The continuum limit is defined by
$\beta\to\infty$ and $V\to\infty$ with fixing $\beta/V$. The $\theta$ term in
Eq. (35) reproduces the $2\pi$ periodicity of observables as a function of
$\theta$ with sufficiently large $\beta$ and $V$ Gattringer:2015baa .

Since the cutoff effects in the Lüscher and Wilson gauge actions should
differ, how to approach the continuum limit should also be different. In
Figure 9, the free energy density $\ln Z/V$ is shown as a function of
$\theta/\pi$ at various $\beta$ with fixing $\beta/V=0.1$. All these results
are obtained by the BTRG with $D=K_{g}=30$. The left and right panels are
obtained by the Lüscher action with $\epsilon=1$ and Wilson action at the same
$\beta$ and $V$. Figure 10 shows the topological charge density $\langle
Q\rangle/V$ defined by Eq. (26) in the same way. With the same $\beta$, the
Lüscher action gives a result closer to the continuum limit than the Wilson
action.

For more quantitative discussion, we compute the topological susceptibility
$\chi_{Q}$ given by Eq. (29) employing the numerical differentiation with the
$O(\Delta^{4})$ accuracy, setting $\Delta=\pi\times 10^{-5}$. In particular,
we investigate the peak position of $\chi_{Q}$ varying $\beta$ and $V$ with
fixing $\beta/V$, fitting $\chi_{Q}$ via three parameters $\chi_{\rm peak}$,
$c$, and $\theta_{\rm c}$ by

$\displaystyle\chi_{Q}(\theta)=\chi_{\rm peak}+c(\theta-\theta_{\rm c})^{2}$
(37)

in the vicinity of the first peak appearing in the region where $\theta$ is
positive. Typically, we use 2000 data points around $\theta_{\rm c}$. Table 2
summarizes the resulting peak position $\theta_{\rm c}$ in the unit of $\pi$.
At the same $\beta$, the Lüscher action results in the peak position closer to
$\theta_{\rm c}=\pi$ than the Wilson action. Since we should observe the peaks
in the topological susceptibility at $\theta=n\pi~{}(n=\pm 1,\pm 3,\cdots)$ in
the limits of $\beta\to\infty$ and $V\to\infty$, we can conclude that the
Lüscher gauge action shows a faster convergence toward the continuum limit
than the Wilson gauge action.

(a) $(\beta,V)=(1.6,2^{4})$
(b) $(\beta,V)=(3.2,2^{5})$
(c) $(\beta,V)=(6.4,2^{6})$
(d) $(\beta,V)=(12.8,2^{7})$

Figure 9: Free energy density $\ln Z/V$ at various $(\beta,V)$ by the Lüscher
gauge action (left) and Wilson gauge action (right) with fixing $\beta/V=0.1$.

(a) $(\beta,V)=(1.6,2^{4})$
(b) $(\beta,V)=(3.2,2^{5})$
(c) $(\beta,V)=(6.4,2^{6})$
(d) $(\beta,V)=(12.8,2^{7})$

Figure 10: Topological charge density $\langle Q\rangle/V$ at various $(\beta,V)$ by the Lüscher gauge action (left) and Wilson gauge action (right) with fixing $\beta/V=0.1$. Table 2:  Comparison of $\theta_{\rm c}/\pi$ obtained by the Lüscher gauge action ($\epsilon=1$) and Wilson gauge action fixing $\beta/V=0.1$.  $\beta$ | Lüscher ($\epsilon=1$) | Wilson
---|---|---
$1.6$ | $1.11932(3)$ | 1.67903(2)
$3.2$ | $1.08112(2)$ | 1.26026(1)
$6.4$ | $1.05111(4)$ | 1.09604(7)
$12.8$ | $1.03070(3)$ | 1.04296(1)

## References

  * (1) M. Lüscher, _Abelian chiral gauge theories on the lattice with exact gauge invariance_ , _Nucl. Phys. B_ 549 (1999) 295–334, [hep-lat/9811032].
  * (2) H. Fukaya and T. Onogi, _Lattice study of the massive Schwinger model with theta term under Lüscher’s “admissibility” condition_ , _Phys. Rev. D_ 68 (2003) 074503, [hep-lat/0305004].
  * (3) H. Fukaya and T. Onogi, _Theta vacuum effects on the chiral condensation and the eta-prime meson correlators in the two flavor massive QED(2) on the lattice_ , _Phys. Rev. D_ 70 (2004) 054508, [hep-lat/0403024].
  * (4) W. Bietenholz, K. Jansen, K. I. Nagai, S. Necco, L. Scorzato and S. Shcheredin, _Exploring topology conserving gauge actions for lattice QCD_ , _JHEP_ 03 (2006) 017, [hep-lat/0511016].
  * (5) W. Bietenholz, C. Czaban, A. Dromard, U. Gerber, C. P. Hofmann, H. Mejía-Díaz et al., _Interpreting Numerical Measurements in Fixed Topological Sectors_ , _Phys. Rev. D_ 93 (2016) 114516, [1603.05630].
  * (6) M. Lüscher and S. Schaefer, _Lattice QCD without topology barriers_ , _JHEP_ 07 (2011) 036, [1105.4749].
  * (7) M. Levin and C. P. Nave, _Tensor renormalization group approach to two-dimensional classical lattice models_ , _Phys. Rev. Lett._ 99 (2007) 120601, [cond-mat/0611687].
  * (8) Z.-C. Gu, F. Verstraete and X.-G. Wen, _Grassmann tensor network states and its renormalization for strongly correlated fermionic and bosonic states_ , 1004.2563.
  * (9) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang and T. Xiang, _Coarse-graining renormalization by higher-order singular value decomposition_ , _Phys. Rev. B_ 86 (Jul, 2012) 045139, [1201.1144].
  * (10) Y. Shimizu and Y. Kuramashi, _Grassmann tensor renormalization group approach to one-flavor lattice Schwinger model_ , _Phys. Rev._ D90 (2014) 014508, [1403.0642].
  * (11) R. Sakai, S. Takeda and Y. Yoshimura, _Higher order tensor renormalization group for relativistic fermion systems_ , _PTEP_ 2017 (2017) 063B07, [1705.07764].
  * (12) D. Adachi, T. Okubo and S. Todo, _Anisotropic Tensor Renormalization Group_ , _Phys. Rev. B_ 102 (2020) 054432, [1906.02007].
  * (13) D. Kadoh and K. Nakayama, _Renormalization group on a triad network_ , 1912.02414.
  * (14) S. Akiyama, Y. Kuramashi, T. Yamashita and Y. Yoshimura, _Restoration of chiral symmetry in cold and dense Nambu–Jona-Lasinio model with tensor renormalization group_ , _JHEP_ 01 (2021) 121, [2009.11583].
  * (15) A. Denbleyker, Y. Liu, Y. Meurice, M. P. Qin, T. Xiang, Z. Y. Xie et al., _Controlling Sign Problems in Spin Models Using Tensor Renormalization_ , _Phys. Rev._ D89 (2014) 016008, [1309.6623].
  * (16) Y. Shimizu and Y. Kuramashi, _Critical behavior of the lattice Schwinger model with a topological term at $\theta=\pi$ using the Grassmann tensor renormalization group_, _Phys. Rev._ D90 (2014) 074503, [1408.0897].
  * (17) H. Kawauchi and S. Takeda, _Tensor renormalization group analysis of CP( $N$-1) model_, _Phys. Rev._ D93 (2016) 114503, [1603.09455].
  * (18) Y. Shimizu and Y. Kuramashi, _Berezinskii-Kosterlitz-Thouless transition in lattice Schwinger model with one flavor of Wilson fermion_ , _Phys. Rev._ D97 (2018) 034502, [1712.07808].
  * (19) J. Unmuth-Yockey, J. Zhang, A. Bazavov, Y. Meurice and S.-W. Tsai, _Universal features of the Abelian Polyakov loop in 1+1 dimensions_ , _Phys. Rev._ D98 (2018) 094511, [1807.09186].
  * (20) D. Kadoh, Y. Kuramashi, Y. Nakamura, R. Sakai, S. Takeda and Y. Yoshimura, _Investigation of complex $\phi^{4}$ theory at finite density in two dimensions using TRG_, _JHEP_ 02 (2020) 161, [1912.13092].
  * (21) Y. Kuramashi and Y. Yoshimura, _Tensor renormalization group study of two-dimensional U(1) lattice gauge theory with a $\theta$ term_, _JHEP_ 04 (2020) 089, [1911.06480].
  * (22) N. Butt, S. Catterall, Y. Meurice, R. Sakai and J. Unmuth-Yockey, _Tensor network formulation of the massless Schwinger model with staggered fermions_ , _Phys. Rev. D_ 101 (2020) 094509, [1911.01285].
  * (23) J. Bloch, R. Lohmayer, S. Schweiss and J. Unmuth-Yockey, _Effective $\mathbb{Z}_{3}$ model for finite-density QCD with tensor networks_, 10, 2021, 2110.09499.
  * (24) S. Takeda, _A novel method to evaluate real-time path integral for scalar $\phi^{4}$ theory_, _PoS_ LATTICE2021 (2022) 532, [2108.10017].
  * (25) J. Bloch, R. G. Jha, R. Lohmayer and M. Meister, _Tensor renormalization group study of the three-dimensional O(2) model_ , _Phys. Rev. D_ 104 (2021) 094517, [2105.08066].
  * (26) K. Nakayama, L. Funcke, K. Jansen, Y.-J. Kao and S. Kühn, _Phase structure of the CP(1) model in the presence of a topological $\theta$-term_, _Phys. Rev. D_ 105 (2022) 054507, [2107.14220].
  * (27) M. Hirasawa, A. Matsumoto, J. Nishimura and A. Yosprakob, _Tensor renormalization group and the volume independence in 2D U(N) and SU(N) gauge theories_ , _JHEP_ 12 (2021) 011, [2110.05800].
  * (28) S. Takeda and Y. Yoshimura, _Grassmann tensor renormalization group for the one-flavor lattice Gross-Neveu model with finite chemical potential_ , _PTEP_ 2015 (2015) 043B01, [1412.7855].
  * (29) Y. Yoshimura, Y. Kuramashi, Y. Nakamura, S. Takeda and R. Sakai, _Calculation of fermionic Green functions with Grassmann higher-order tensor renormalization group_ , _Phys. Rev._ D97 (2018) 054511, [1711.08121].
  * (30) D. Kadoh, Y. Kuramashi, Y. Nakamura, R. Sakai, S. Takeda and Y. Yoshimura, _Tensor network formulation for two-dimensional lattice $\mathcal{N}$ = 1 Wess-Zumino model_, _JHEP_ 03 (2018) 141, [1801.04183].
  * (31) S. Akiyama and Y. Kuramashi, _Tensor renormalization group approach to (1+1)-dimensional Hubbard model_ , _Phys. Rev. D_ 104 (2021) 014504, [2105.00372].
  * (32) J. Bloch and R. Lohmayer, _Grassmann higher-order tensor renormalization group approach for two-dimensional strong-coupling QCD_ , _Nucl. Phys. B_ 986 (2023) 116032, [2206.00545].
  * (33) S. Akiyama, _Matrix product decomposition for two- and three-flavor Wilson fermions: Benchmark results in the lattice Gross-Neveu model at finite density_ , _Phys. Rev. D_ 108 (2023) 034514, [2304.01473].
  * (34) A. Yosprakob, J. Nishimura and K. Okunishi, _A new technique to incorporate multiple fermion flavors in tensor renormalization group method for lattice gauge theories_ , _JHEP_ 11 (2023) 187, [2309.01422].
  * (35) M. Asaduzzaman, S. Catterall, Y. Meurice, R. Sakai and G. C. Toga, _Tensor network representation of non-abelian gauge theory coupled to reduced staggered fermions_ , _JHEP_ 05 (2024) 195, [2312.16167].
  * (36) Z. Komargodski, A. Sharon, R. Thorngren and X. Zhou, _Comments on Abelian Higgs Models and Persistent Order_ , _SciPost Phys._ 6 (2019) 003, [1705.04786].
  * (37) C. Gattringer and O. Orasch, _Density of states approach for lattice gauge theory with a $\theta$-term_, _Nucl. Phys. B_ 957 (2020) 115097, [2004.03837].
  * (38) M. Hirasawa, A. Matsumoto, J. Nishimura and A. Yosprakob, _Complex Langevin analysis of 2D U(1) gauge theory on a torus with a $\theta$ term_, _JHEP_ 09 (2020) 023, [2004.13982].
  * (39) C. Gattringer, D. Göschl and T. Sulejmanpasic, _Dual simulation of the 2d U(1) gauge Higgs model at topological angle $\theta=\pi\,$: Critical endpoint behavior_, _Nucl. Phys. B_ 935 (2018) 344–364, [1807.07793].
  * (40) S. Akiyama, D. Kadoh, Y. Kuramashi, T. Yamashita and Y. Yoshimura, _Tensor renormalization group approach to four-dimensional complex $\phi^{4}$ theory at finite density_, _JHEP_ 09 (2020) 177, [2005.04645].
  * (41) D. Adachi, T. Okubo and S. Todo, _Bond-weighted tensor renormalization group_ , _Phys. Rev. B_ 105 (Feb, 2022) L060402, [2011.01679].
  * (42) Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng and T. Xiang, _Second Renormalization of Tensor-Network States_ , _Phys. Rev. Lett._ 103 (2009) 160601, [0809.0182].
  * (43) H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai and T. Xiang, _Renormalization of tensor-network states_ , _Phys. Rev. B_ 81 (May, 2010) 174411.
  * (44) S. Akiyama, _Bond-weighting method for the Grassmann tensor renormalization group_ , _JHEP_ 11 (2022) 030, [2208.03227].
  * (45) Z.-C. Gu and X.-G. Wen, _Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order_ , _Phys. Rev. B_ 80 (Oct, 2009) 155131.
  * (46) A. Ueda and M. Oshikawa, _Resolving the berezinskii-kosterlitz-thouless transition in the two-dimensional xy model with tensor-network-based level spectroscopy_ , _Phys. Rev. B_ 104 (Oct, 2021) 165132.
  * (47) C.-Y. Huang, S.-H. Chan, Y.-J. Kao and P. Chen, _Tensor network based finite-size scaling for two-dimensional Ising model_ , _Phys. Rev. B_ 107 (2023) 205123, [2302.02585].
  * (48) A. Ueda and M. Oshikawa, _Finite-size and finite bond dimension effects of tensor network renormalization_ , _Phys. Rev. B_ 108 (Jul, 2023) 024413.
  * (49) F. I. Az-zahra, S. Takeda and T. Yamazaki, _Spectroscopy by Tensor Renormalization Group Method_ , 2404.15666.
  * (50) A. Ueda, _Renormalization group flow and fixed-point in tensor network representations_ , 2024.
  * (51) https://qsw.phys.s.u-tokyo.ac.jp/.
  * (52) C. Gattringer, T. Kloiber and M. Müller-Preussker, _Dual simulation of the two-dimensional lattice U(1) gauge-Higgs model with a topological term_ , _Phys. Rev. D_ 92 (2015) 114508, [1508.00681].