# Role of Bootstrap Averaging in Generalized Approximate Message Passing

Takashi Takahashi  Institute for Physics of Intelligence, The University of
Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan, takashi-
<EMAIL_ADDRESS>
###### Abstract

Generalized approximate message passing (GAMP) is a computationally efficient
algorithm for estimating an unknown signal
$\mathbf{\boldsymbol{w}}_{0}\in\mathbb{R}^{N}$ from a random linear
measurement
$\mathbf{\boldsymbol{y}}=X\mathbf{\boldsymbol{w}}_{0}+\mathbf{\boldsymbol{\epsilon}}\in\mathbb{R}^{M}$,
where $X\in\mathbb{R}^{M\times N}$ is a known measurement matrix and
$\mathbf{\boldsymbol{\epsilon}}$ is the noise vector. The salient feature of
GAMP is that it can provide an unbiased estimator
$\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}\sim\mathcal{N}(\mathbf{\boldsymbol{w}}_{0},\hat{s}^{2}I_{N})$, which can be
used for various hypothesis-testing methods. In this study, we consider the
bootstrap average of an unbiased estimator of GAMP for the elastic net. By
numerically analyzing the state evolution of _approximate message passing with
resampling_ , which has been proposed for computing bootstrap statistics of
the elastic net estimator, we investigate when the bootstrap averaging reduces
the variance of the unbiased estimator and the effect of optimizing the size
of each bootstrap sample and hyperparameter of the elastic net regularization
in the asymptotic setting $M,N\to\infty,M/N\to\alpha\in(0,\infty)$. The
results indicate that bootstrap averaging effectively reduces the variance of
the unbiased estimator when the actual data generation process is inconsistent
with the sparsity assumption of the regularization and the sample size is
small. Furthermore, we find that when $\mathbf{\boldsymbol{w}}_{0}$ is less
sparse, and the data size is small, the system undergoes a phase transition.
The phase transition indicates the existence of the region where the ensemble
average of unbiased estimators of GAMP for the elastic net norm minimization
problem yields the unbiased estimator with the minimum variance.

##  I Introduction

Consider estimating an unknown signal
$\mathbf{\boldsymbol{w}}_{0}\in\mathbb{R}^{N}$ from a random linear
measurement $\mathbf{\boldsymbol{y}}\in\mathbb{R}^{M}$ in the form,

$\mathbf{\boldsymbol{y}}=X\mathbf{\boldsymbol{w}}_{0}+\mathbf{\boldsymbol{\epsilon}}.$
(1)

In (1), $X\in\mathbb{R}^{M\times N}$ is a known measurement matrix whose
elements are independent and identically distributed (i.i.d.) standard
Gaussian variables, and
$\mathbf{\boldsymbol{\epsilon}}\sim\mathcal{N}(\mathbf{\boldsymbol{0}},\Delta
I_{M})$ is the measurement noise. We also assume that each element of the
unknown signal $\mathbf{\boldsymbol{w}}_{0}$ is i.i.d., according to a
distribution $q_{0}$.

Generalized approximate message passing (GAMP) [1, 2] is a computationally
efficient algorithm for solving this problem. A striking feature of GAMP is
its applicability to various hypothesis testing [3]. Specifically, GAMP can
provide an unbiased estimator $\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}\sim\mathcal{N}(\mathbf{\boldsymbol{w}}_{0},\hat{s}^{2}I_{M})$ [2] in a
high-dimensional asymptotic setting with
$M,N\to\infty,M/N\to\alpha\in(0,\infty)$, where the variance $\hat{s}^{2}$
depends on the quality of the measurement $\mathbf{\boldsymbol{y}}$ and
denoising function used in GAMP. This unbiased estimator has been used to test
the significance of estimated signals [4, 5, 3, 6].

The statistical power of these tests depends on the variance $\hat{s}^{2}$,
with a lower variance leading to a higher statistical power. However, reducing
the variance $\hat{s}^{2}$ is not a trivial task. Replacing the denoising
function used in GAMP with a powerful one based on nonconvex regularization,
for example, can worsen the convergence of GAMP [7] or, even if convergence is
achieved, the improvement may be insignificant [8]. This study aims to find an
alternative way to reduce the variance without nonconvex regularization.

To reduce the variance, we use the _bootstrap averaging_ [9] of computational
statistics (commonly known as _ensemble learning_ in machine learning [10,
11]). Specifically, we consider averaging the unbiased estimators of GAMP for
multiple bootstrap samples with arbitrary size
$M\mu_{B},\,\mu_{B}\in(0,\infty]$. For denoising function of GAMP, we consider
the one for the elastic net [12].

However, the efficient computation and theoretical analysis of an averaged
unbiased estimator remains unresolved. To resolve this problem, we use _AMP
with resampling_ (AMPR) [13, 14] that has been proposed for computing
bootstrap statistics of the elastic net estimator by running a variant of GAMP
once. We will argue that AMPR is actually computing the bootstrap average of
the unbiased estimators of GAMP. That is, the averaged unbiased estimator can
be computed efficiently, and its variance can be analyzed using the state
evolution (SE) of AMPR, which has been developed to analyze the performance of
AMPR in [14, 13]. We then conduct a thorough numerical analysis of the SE of
the AMPR to investigate when bootstrap averaging reduces the variance of the
unbiased estimator, and what phenomena occur when optimizing the bootstrap
sample size and the hyperparameter of the elastic net regularization.

The findings of this study are summarized as follows:

  * •

As mentioned above, the averaged unbiased estimator is obtained by AMPR.
Furthermore, its variance can be estimated using the output of AMPR without
knowing the actual signal $\mathbf{\boldsymbol{w}}_{0}$. Thus, we can minimize
the variance by adjusting the bootstrap sample size and the hyperparameters of
the elastic net (see Section III and IV)

  * •

The variance of the averaged unbiased estimator can be reduced via
bootstrapping, especially when the true data generation process is
inconsistent with the sparsity assumption of the regularization and the data
size is insufficient (see Section V-B)

  * •

When $\mathbf{\boldsymbol{w}}_{0}$ is less sparse, and the data size is small,
a phase transition occurs. This phase transition indicates the existence of
the region where the value of the regularization parameter is infinitesimally
small, and the number of unique data points in each bootstrap sample is less
than the dimension of $\mathbf{\boldsymbol{w}}_{0}$. That is, in this region,
the ensemble average of the unbiased estimators of GAMP for the elastic norm
minimization problem (also known as the minimum norm _interpolation_ in
machine learning [15, 16, 17]) yields the best averaged unbiased estimator
(see Section V-C)

###  I-A Notation

$\mathcal{N}(\mu,\sigma^{2})$ denotes a Gaussian distribution with mean $\mu$
and variance $\sigma^{2}$ and ${\rm Poisson}(\mu_{B})$ denotes a Poisson
distribution with mean $\mu_{B}$. For a random variable $X\sim p_{X}$, we
denote by $\mathbb{E}_{X}[\dots]$ an average $\int(\dots)p_{X}(x)dx$. Given a
vector $\mathbf{\boldsymbol{x}}\in\mathbb{R}^{N}$ and scalar function
$f:\mathbb{R}\to\mathbb{R}$, we write $f(\mathbf{\boldsymbol{x}})$ for the
vector obtained by applying $f$ componentwise. For a vector
$\mathbf{\boldsymbol{x}}=(x_{1},x_{2},\dots,x_{N})^{\top}\in\mathbb{R}^{N}$,
we denote by
$\mathbf{\boldsymbol{x}}^{2}=(x_{1}^{2},x_{2}^{2},\dots,x_{N}^{2})^{\top}$ the
componentwise operations and by
$\langle\mathbf{\boldsymbol{x}}\rangle=N^{-1}\sum_{i=1}^{N}x_{i}$ the
empirical average.

##  II Background on AMPR

Algorithm 1 AMPR

1:Measurement matrix $X\in\mathbb{R}^{M\times N}$, measurement
$\mathbf{\boldsymbol{y}}\in\mathbb{R}^{M}$, denoising function $g$ and its
derivative $g^{\prime}$ from (4) and (7), the bootstrap sample size
$\mu_{B}\in(0,\infty]$, and the number of iterations $T_{\rm it}$.

2:Select initial
$\mathbf{\boldsymbol{h}}_{0}\in\mathbb{R}^{N},\hat{Q}_{0},\hat{v}_{0}\in(0,\infty)$.

3:Set
$\mathbf{\boldsymbol{z}}_{0}=\mathbf{\boldsymbol{0}}_{M},\alpha_{N}=M/N$.

4:Let $\mathbf{\boldsymbol{\xi}}$ and $c$ be random variables distributed as
$\mathbf{\boldsymbol{\eta}}\sim\mathcal{N}(0,I_{N})$ and $c\sim{\rm
Poisson}(\mu_{B})$.

5:for $t=0,1,\dots,T_{\rm it}-1$ do

6:
$\hat{\mathbf{\boldsymbol{w}}}_{t}=\mathbb{E}_{\mathbf{\boldsymbol{\eta}}}[g(\mathbf{\boldsymbol{h}}_{t}+\sqrt{\hat{v}_{t}}\mathbf{\boldsymbol{\eta}},\hat{Q}_{t})]$.

7:
$\chi_{t}=\langle\mathbb{E}_{\mathbf{\boldsymbol{\eta}}}[g^{\prime}(\mathbf{\boldsymbol{h}}_{t}+\sqrt{\hat{v}_{t}}\mathbf{\boldsymbol{\eta}},\hat{Q}_{t})]\rangle$.

8:
$v_{t}=\langle\mathbb{E}_{\mathbf{\boldsymbol{\eta}}}[g(\mathbf{\boldsymbol{h}}_{t}+\sqrt{\hat{v}_{t}}\mathbf{\boldsymbol{\eta}},\hat{Q}_{t})^{2}]-\hat{\mathbf{\boldsymbol{w}}}_{t}^{2}\rangle$.

9:
$f^{(1)}_{t+1}=\mathbb{E}_{c}[\frac{c/\mu_{B}}{1+c/\mu_{B}\chi_{t}}],\,f^{(2)}_{t+1}=\mathbb{E}_{c}[(\frac{c/\mu_{B}}{1+c/\mu_{B}\chi_{t}})^{2}]$.

10: $\hat{Q}_{t+1}=\alpha_{N}f^{(1)}_{t+1}$.

11:
$\mathbf{\boldsymbol{a}}_{t+1}=f_{t+1}^{(1)}(\mathbf{\boldsymbol{y}}-X\hat{\mathbf{\boldsymbol{w}}_{t}}+\chi_{t}\mathbf{\boldsymbol{a}}_{t})$.

12:
$\mathbf{\boldsymbol{h}}_{t+1}=X^{\top}\mathbf{\boldsymbol{a}}_{t}+\hat{Q}_{t+1}\hat{\mathbf{\boldsymbol{w}}}_{t}$.

13:
$\hat{v}_{t+1}=\alpha_{N}(f_{t+1}^{(2)}v_{t}+\frac{f_{t+1}^{(2)}-(f_{t+1}^{(1)})^{2}}{(f_{t+1}^{(t)})^{2}}\langle\mathbf{\boldsymbol{a}}_{t+1}^{2}\rangle)$
.

14:end for

15:Return $(\mathbf{\boldsymbol{h}}_{T_{\rm it}},\hat{Q}_{T_{\rm
it}},\hat{v}_{T_{\rm it}})$.

###  II-A AMPR

Algorithm 1 shows AMPR [13] with an elastic net denoising function. Function
$g:\mathbb{R}\times(0,\infty)\to\mathbb{R}$ is the elastic net denoising
function and $g^{\prime}$ is a derivative of $g$ with respect to the first
argument:

$\displaystyle g(h,\hat{Q})$ $\displaystyle=$
$\displaystyle\left\\{\,\begin{IEEEeqnarraybox}[]{[}][c]{l?s}\IEEEstrut
0&if$|h| \leq\lambda\gamma$,\\\ \frac{h-{\rm
sgn}(h)\gamma\lambda}{\hat{Q}+\lambda(1-\gamma)}&otherwise,\IEEEstrut\end{IEEEeqnarraybox}\right.$
(4) $\displaystyle g^{\prime}(h,\hat{Q})$ $\displaystyle=$
$\displaystyle\left\\{\,\begin{IEEEeqnarraybox}[]{[}][c]{l?s}\IEEEstrut
0&if$|h| \leq\lambda\gamma$,\\\
\frac{1}{\hat{Q}+\lambda(1-\gamma)}&otherwise.\IEEEstrut\end{IEEEeqnarraybox}\right.$
(7)

where $\lambda>0$ represents the regularization strength and $\gamma\in[0,1]$
is the $\ell_{1}$-ratio. At a fixed point, AMPR offers bootstrap statistics of
the elastic net estimator as follows:

######  Proposition 1 (bootstrap statistics based on AMPR [13])

Let $\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ be the elastic net estimator for a
bootstrap sample $D^{\ast}$ of size $\mu_{B}M$.

$\hat{\mathbf{\boldsymbol{w}}}^{\ast}=\mathop{\rm
argmin}_{\mathbf{\boldsymbol{w}}\in\mathbb{R}^{N}}\sum_{\mu=1}^{M}\frac{c_{\mu}}{2\mu_{B}}(y_{\mu}-\mathbf{\boldsymbol{x}}_{\mu}^{\top}\mathbf{\boldsymbol{w}})^{2}+\sum_{i=1}^{N}\lambda(\gamma|w_{i}|+\frac{1-\gamma}{2}w_{i}^{2}),$
(8)

where $c_{\mu}\sim_{\rm i.i.d.}{\rm Poisson}(\mu_{B})$ represents the number
of times the data point $(\mathbf{\boldsymbol{x}}_{\mu},y_{\mu})$,
$\mathbf{\boldsymbol{x}}_{\mu}$ being the $\mu$-th row of $X$, appears in the
bootstrap sample $D^{\ast}$. Then once the AMPR reaches its fixed point at
sufficiently large $T_{\rm it}$, the bootstrap statistics of
$\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ can be computed as

$\mathbb{E}_{\mathbf{\boldsymbol{c}}}\left[\psi(\hat{w}_{i}^{\ast})\right]=\mathbb{E}_{\eta}[\psi(g(h_{i}+\sqrt{\hat{v}}\eta,\hat{Q}))],\quad\eta\sim\mathcal{N}(0,1),$
(9)

where $\psi:\mathbb{R}\to\mathbb{R}$ is such that the expression in (9) is
well-defined and otherwise arbitrary. The variables without iteration indexes
$(h_{i},\hat{v},\hat{Q})$ are the output of AMPR at a fixed point.

Note that the averages
$\mathbb{E}_{\eta}[g(h+\sqrt{\hat{v}}\eta;\hat{Q})],\mathbb{E}_{\eta}[g(h+\sqrt{\hat{v}}\eta;\hat{Q})^{2}]$
and $\mathbb{E}_{\eta}[g^{\prime}(h+\sqrt{\hat{v}}\eta;\hat{Q})]$ can be
written in closed-form by the error function, thus computing these quantities
is computationally easy. Also, $\mathbb{E}_{c}[\dots]$ is an average over a
one-dimensional discrete random variable. It can be computed numerically
without much computational overhead. Hence the computational complexity of
computing the RHS of (9) is dominated by the matrix-vector product operations
in lines 11-12 of Algorithm 1 instead of repeatedly computing
$\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ for numerous realizations of
$\mathbf{\boldsymbol{c}}$, making AMPR a computationally efficient algorithm
for computing the bootstrap statistics.

###  II-B SE of AMPR

Algorithm 2 SE of AMPR

1:Initial state of AMPR $(\mathcal{E}_{0},\hat{Q}_{1},\hat{v}_{1})$, variance
of the measurement noise $\Delta$, distribution $q_{0}$ of the signal
$\mathbf{\boldsymbol{w}}_{0}$, measurement ratio $\alpha$, denoising function
$g$ and its derivative $g^{\prime}$ from (4) and (7), the bootstrap sample
size $\mu_{B}\in(0,\infty]$, and the number of iterations $T_{\rm it}$.

2:Set $\hat{\chi}_{1}=\alpha^{-1}\hat{Q}_{1}^{2}(\mathcal{E}_{0}+\Delta)$

3:Let $\xi,\eta,w_{0}$ and $c$ be random variables distributed as
$\xi,\eta\sim\mathcal{N}(0,1),w_{0}\sim q_{0}$, and $c\sim{\rm
Poisson}(\mu_{B})$.

4:for $t=1,\dots,T_{\rm it}-1$ do

5:
$\hat{w}_{t}=g(\hat{Q}_{t}w_{0}+\sqrt{\hat{\chi}_{t}}\xi+\sqrt{\hat{v}_{t}}\eta;\hat{Q}_{t})$.

6:
$\mathcal{E}_{t}=\mathbb{E}_{w_{0},\xi}[(\mathbb{E}_{\eta}[\hat{w}_{t}]-w_{0})^{2}]$

7:
$\chi_{t}=\mathbb{E}_{w_{0},\xi,\eta}[g^{\prime}(\hat{Q}_{t}w_{0}+\sqrt{\hat{\chi}_{t}}\xi+\sqrt{\hat{v}_{t}}\eta,\hat{Q}_{t})]$.

8:
$v_{t}=\mathbb{E}_{w_{0},\xi}[\mathbb{E}_{\eta}[\hat{w}_{t}^{2}]-\mathbb{E}_{\eta}[\hat{w}_{t}]^{2}]$.

9:
$f^{(1)}_{t+1}=\mathbb{E}_{c}[\frac{c/\mu_{B}}{1+c/\mu_{B}\chi_{t}}],\,f^{(2)}_{t+1}=\mathbb{E}_{c}[(\frac{c/\mu_{B}}{1+c/\mu_{B}\chi_{t}})^{2}]$.

10: $\hat{Q}_{t+1}=\alpha f^{(1)}_{t+1}$.

11: $\hat{\chi}_{t+1}=\alpha(f_{t+1}^{(1)})^{2}(\mathcal{E}_{t}+\Delta)$

12:
$\hat{v}_{t+1}=\alpha(f_{t+1}^{(2)}v_{t}+(f_{t+1}^{(2)}-(f_{t+1}^{(1)})^{2})\mathcal{E}_{t})$
.

13:end for

AMPR displays remarkable behavior. Let
$\hat{\mathbf{\boldsymbol{r}}}_{t}=\mathbf{\boldsymbol{h}}_{t}/\hat{Q}_{t},t=1,2,\dots,T_{\rm
it}.$ Then $\hat{\mathbf{\boldsymbol{r}}}_{t}$ behaves like a white Gaussian
noise-corrupted version of the true signal $\mathbf{\boldsymbol{w}}_{0}$ [13].
Furthermore, the variance can be estimated using SE.

######  Proposition 2 (SE of AMPR [13])

$\hat{\mathbf{\boldsymbol{r}}}_{t}$ behaves as a white Gaussian noise-
corrupted version of the actual signal $\mathbf{\boldsymbol{w}}_{0}$:

$\hat{\mathbf{\boldsymbol{r}}}_{t}\sim\mathcal{N}(\mathbf{\boldsymbol{w}}_{0},\hat{\sigma}_{t}^{2}),\quad\hat{\sigma}^{2}_{t}=\hat{\chi}_{t}/\hat{Q}_{t}^{2},$
(10)

for some positive value $\hat{\chi}_{t}$, indicating that
$\hat{\mathbf{\boldsymbol{r}}}_{t}$ can be used as an unbiased estimate of
$\mathbf{\boldsymbol{w}}_{0}$. The variance is predicted in the asymptotic
setting $M,N\to\infty,M/N\to\alpha\in(0,\infty)$ using the scalar SE specified
in Algorithm 2. There, $\mathcal{E}_{t},t=1,2,\dots$ corresponds to the mean
squared error (MSE) of the AMPR estimate $\hat{\mathbf{\boldsymbol{w}}}_{t}$:
$\mathcal{E}_{t}=N^{-1}\|\hat{\mathbf{\boldsymbol{w}}}_{t}-\mathbf{\boldsymbol{w}}_{0}\|_{2}^{2}$.
To track the performance of AMPR, $\mathcal{E}_{0}$ should be inputted as the
MSE of $\hat{\mathbf{\boldsymbol{w}}}_{0}$.

Using the SE of the AMPR, we can predict the variance of the unbiased
estimator in the asymptotic setting $M,N\to\infty,M/N\to\alpha\in(0,\infty)$
for each value of $\mu_{B},\lambda$, and $\gamma$. Hence variance can be
minimized by tuning $(\mu_{B},\lambda,\gamma)$ using versatile black-box
optimization methods implemented in various optimization libraries [18, 19].

##  III Averaged unbiased estimator

Here, we explain the meaning of the SE of AMPR. Subsequently, we argue that
$\hat{\mathbf{\boldsymbol{r}}}_{t}$ of AMPR is the bootstrap average of the
unbiased estimators of GAMP.

The first point is the meaning of $\xi$ and $\eta$ that appear in SE.
Propositions 1 and 2 indicate that, once the AMPR reaches its fixed point at
sufficiently large $T_{\rm it}$, for any functions
$\phi,\psi:\mathbb{R}\to\mathbb{R}$ such that the following expression is well
defined, the following holds in the asymptotic setting
$M,N\to\infty,M/N\to\alpha\in(0,\infty)$:

$\frac{1}{N}\sum_{i=1}^{N}\phi\left(\mathbb{E}_{\mathbf{\boldsymbol{c}}}[\psi(\hat{w}_{i}^{\ast})]\right)\to\mathbb{E}_{w_{0},\xi}\left[\phi\left(\mathbb{E}_{\eta}[\psi(g(\hat{Q}w_{0}\right.\right.\\\
\left.\left.+\sqrt{\hat{\chi}}\xi+\sqrt{\hat{v}}\eta;\hat{Q}))]\right)\right],$
(11)

where the variables $\hat{Q},\hat{\chi},\hat{v}$ are the outputs of SE at a
fixed point. In (11), one can interpret that $\xi$ and $\eta$ effectively
represents the randomness originating from $(X,\mathbf{\boldsymbol{y}})$ and
$\mathbf{\boldsymbol{c}}$, respectively.

The next point is the relationship between AMPR and GAMP. For this, we define
$\bar{q}_{t}\equiv\mathcal{E}_{t}+v_{t}$ and
$\bar{\chi}_{t}=\hat{\chi}_{t}+\hat{v}_{t}$. Then, the evolution of
$(\hat{Q}_{t},\bar{\chi}_{t},\bar{q}_{t},\chi_{t})$ is equivalent to the SE of
GAMP for computing $\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ for one realization
of $\mathbf{\boldsymbol{c}}$ [1, 2]. Recall that
$\hat{\mathbf{\boldsymbol{w}}}^{\ast}$ is specified in (8). In addition,
$\bar{\chi}_{t}/\hat{Q}_{t}^{2}\,(=(\mathcal{E}_{t}+\hat{v}_{t})/\hat{Q}_{t}^{2})$
represents the variance of the unbiased estimate of GAMP
$\hat{\mathbf{\boldsymbol{r}}}_{t}^{\rm G}(\mathbf{\boldsymbol{c}})$ at each
iteration. These observations indicate that the unbiased estimator of GAMP at
each iteration $t$ can be decomposed into two Gaussian variables:

$\displaystyle\hat{\mathbf{\boldsymbol{r}}}_{t}^{\rm
G}(\mathbf{\boldsymbol{c}})$ $\displaystyle=$
$\displaystyle\mathbf{\boldsymbol{w}}_{0}+\frac{\sqrt{\bar{\chi}_{t}}}{\hat{Q}_{t}}\mathbf{\boldsymbol{\xi_{0}}}=\mathbf{\boldsymbol{w}}_{0}+\frac{\sqrt{\hat{\chi}_{t}}}{\hat{Q}_{t}}\mathbf{\boldsymbol{\xi}}+\frac{\sqrt{\hat{v}_{t}}}{\hat{Q}_{t}}\mathbf{\boldsymbol{\eta}},$
(12)

where
$\mathbf{\boldsymbol{\xi}}_{0},\mathbf{\boldsymbol{\xi}},\mathbf{\boldsymbol{\eta}}\sim\mathcal{N}(0,I_{N})$,
and $\mathbf{\boldsymbol{\xi}}$ and $\mathbf{\boldsymbol{\eta}}$ represent the
randomness originating from $(X,\mathbf{\boldsymbol{y}})$ and
$\mathbf{\boldsymbol{c}}$, respectively. We can interpret that AMPR computes
the bootstrap statistics by decomposing $\mathbf{\boldsymbol{\xi}}_{0}$ into
$\mathbf{\boldsymbol{\xi}}$ and $\mathbf{\boldsymbol{\eta}}$, and averaging
over $\mathbf{\boldsymbol{\eta}}$. This yields the unbiased estimator
$\hat{\mathbf{\boldsymbol{r}}}_{t}$ as in (10) that is the bootstrap average
of $\hat{\mathbf{\boldsymbol{r}}}^{\rm G}(\mathbf{\boldsymbol{cka}})$. We
verify this point numerically in section V-A (see Fig. 1 and 2).

Figure 1: Q-Q plot of
$\hat{\mathbf{\boldsymbol{r}}}-\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})$ for one realization of
$X,\mathbf{\boldsymbol{y}}$, and $\mathbf{\boldsymbol{c}}$. The parameters are
set as
$(N,\alpha,\Delta,\lambda,\gamma,\mu_{B},\rho)=(4096,0.25,0.1,0.5,0.5,0.1)$.
Figure 2: Scatter plot for $\hat{\mathbf{\boldsymbol{r}}}$ and
$\mathbb{E}_{\mathbf{\boldsymbol{c}}}[\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})]$ for one realization of
$X,\mathbf{\boldsymbol{y}}$. The parameters are set as
$(N,\alpha,\Delta,\lambda,\gamma,\mu_{B},\rho)=(4096,0.8,0.25,0.1,0.5,0.5,0.1)$.
$\mathbb{E}_{\mathbf{\boldsymbol{c}}}[\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})]$ was computed from $65536$ realizations of
$\mathbf{\boldsymbol{c}}$.

(a)

(b)

(c)

Figure 3:  Ratio of optimal variances $\hat{\sigma}^{2}/\hat{s}^{2}$ is
plotted against the sparsity of the true signal $\rho$ and the measurement
ratio $\alpha$ as heat maps. In panels (a) and (b), the $\ell_{1}$-ratio of
the elastic net regularization is fixed. In panel (c), the $\ell_{1}$-ratio is
also optimized. The measurement noise is set as $\Delta=0.15$.

##  IV Variance optimization

In this section, we describe the properties of the variance of the unbiased
estimate of AMPR.

First, the following proposition states that bootstrap averaging for the
optimal choice of $(\mu_{B},\lambda,\gamma)$ does not increase the variance of
the unbiased estimator.

###### Proposition 3

Let $\hat{s}^{2}$ be the variance of the unbiased estimator of GAMP at a fixed
point to compute the elastic net estimator $\hat{\mathbf{\boldsymbol{w}}}$ of
$D=(X,\mathbf{\boldsymbol{y}})$: $\hat{\mathbf{\boldsymbol{w}}}=\mathop{\rm
argmin}_{\mathbf{\boldsymbol{w}}\in\mathbb{R}^{N}}\frac{1}{2}\sum_{\mu=1}^{M}(y_{\mu}-\mathbf{\boldsymbol{a}}_{\mu}^{\top}\mathbf{\boldsymbol{w}})^{2}+\sum_{i=1}^{N}\lambda(\gamma|w_{i}|+\frac{1-\gamma}{2})w_{i}^{2}$.
We can select the parameter $(\mu_{B}^{\star},\lambda^{\star},\gamma^{\star})$
such that the variance of AMPR’s unbiased estimator $\hat{\sigma}^{2}$ at a
fixed point does not exceed $\hat{s}^{2}$.

###### Proof:

We denote by $\hat{v},v$ the values of $\hat{v}_{t}$ and $v_{t}$ at the fixed
point of the SE of AMPR. We consider the limit $\mu_{B}\to\infty$. In this
limit, the random variable $c/\mu_{B},c\sim{\rm Poisson}(\mu_{B})$ that
appears in (8) and Algorithm 1 behaves deterministically; the mean and
variance of $c/\mu_{B}$ converge to $1$ and $0$. Then proposition 1 indicates
that $\hat{v}=0$. Moreover, from line 12 in Algorithm 2, $\hat{v}=0$ implies
$v=0$, which yields the GAMP algorithm for computing the elastic net estimator
of $D$ [1]. ∎

Although SE prediction of the variance of $\hat{\mathbf{\boldsymbol{r}}}_{t}$
requires information on the unknown signal $\mathbf{\boldsymbol{w}}_{0}$, we
can predict the variance from the data. In other words, estimating the
variance of $\hat{\mathbf{\boldsymbol{r}}}_{t}$ does not require explicit
knowledge on $\mathbf{\boldsymbol{w}}_{0}$.

######  Proposition 4 (Variance estimation from data)

In the asymptotic setting $M,N\to\infty,M/N\to\alpha\in(0,\infty)$, the
variance $\hat{\sigma}_{t}^{2}$ of the unbiased estimate
$\hat{\mathbf{\boldsymbol{r}}}_{t}$ can be estimated as

$\hat{\sigma}_{t}^{2}=\alpha\langle\mathbf{\boldsymbol{a}}_{t}^{2}\rangle/\hat{Q}_{t}^{2}.$
(13)

###### Proof:

In SE of AMPR, $\hat{\chi}_{t}/\hat{Q}_{t}^{2}$ is determined by the MSE
$\mathcal{E}_{t-1}$ and variance of measurement noise $\Delta$ as
$\hat{\chi}_{t}/\hat{Q}_{t}^{2}=\alpha^{-1}(\mathcal{E}_{t-1}+\Delta)$. For
linear models, $\mathcal{E}_{t-1}+\Delta$ corresponds to the prediction error
for a new sample and can be estimated using the leave-one-out error (LOOE)
[20]. LOOE can be estimated from $\mathbf{\boldsymbol{a}}_{t}$ because it is
proportional to the leave-one-out estimate for the data point $\mu$:
$a_{t,\mu}=\hat{Q}_{t}\alpha^{-1}(y_{\mu}-\mathbf{\boldsymbol{x}}_{\mu}^{\top}\hat{\mathbf{\boldsymbol{w}}}_{t}^{\backslash\mu})$,
where $\hat{\mathbf{\boldsymbol{w}}}_{t}^{\backslash\mu}$ is the AMPR’s
estimate of $\mathbf{\boldsymbol{w}}_{0}$ without the sample $\mu$ (Equation
(19) of reference [13]). ∎

Propositions 3 and 4 indicate that the variance $\hat{\sigma}^{2}$ can be
minimized even if the signal $\mathbf{\boldsymbol{w}}_{0}$ is unknown.
However, we use SE for the theoretical assessment in the next section for
convenience.

##  V Numerical analysis

In the sequel, by numerically minimizing the variance using SE, we investigate
when bootstrapping reduces the variance of the unbiased estimator and the
phenomena that occur when optimizing the bootstrap sample size $\mu_{B}$ and
the hyperparameter of the elastic net regularization $(\lambda,\gamma)$.

For this, we searched for the optimal parameter
$(\mu_{B}^{\star},\lambda^{\star},\gamma^{\star})$ that yielded the minimum
variance using the SE of AMPR and the Nelder-Mead algorithm in the _Optim.jl_
library [18]. We obtained the fixed point of the SE by iterating the SE a
sufficient number of times. For comparison, the same optimization was
performed for the non-bootstrap case. For the signal distribution $q_{0}$, we
consider the Gauss-Bernoulli model:
$q_{0}=\rho\delta_{0}+(1-\rho)\mathcal{N}(0,1)$, with sparsity $\rho\in(0,1)$.

In this section, we denote the outputs of the AMPR or GAMP at fixed points by
unindexed variables.

###  V-A Distribution of the unbiased estimator

We verified the interpretation of AMPR’s output
$\hat{\mathbf{\boldsymbol{r}}}$ described in Section III. For this, we
compared the output of GAMP $\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})$ for each realization of $\mathbf{\boldsymbol{c}}$
as in (8), and the output of AMPR. The parameters used to produce the figure
were set as
$(N,\alpha,\Delta,\lambda,\gamma,\mu_{B},\rho)=(4096,0.8,0.25,0.1,0.5,0.5,0.1)$.

Fig. 1 shows the sample quantiles of
$\hat{\mathbf{\boldsymbol{r}}}-\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})$ versus the normal distribution with variance
$\hat{v}/\hat{Q}^{2}$. The scattered points are approximately aligned with a
line with a slope of $1$ and an intercept of $0$. Fig. 2 shows the scatter
plot of $\hat{\mathbf{\boldsymbol{r}}}$ versus
$\mathbb{E}_{\mathbf{\boldsymbol{c}}}[\hat{\mathbf{\boldsymbol{r}}}^{\rm
G}(\mathbf{\boldsymbol{c}})]$. Again, the scattered points are approximately
aligned with a line with a slope of $1$ and an intercept of $0$. This is
consistent with the decomposition of the Gaussian noise (12). Thus, AMPR
computes the bootstrap average of the unbiased estimator of GAMP.

###  V-B Variance reduction

We quantitatively compare the variance $\hat{\sigma}^{2}$ of the averaged
unbiased estimator with that of $\hat{s}^{2}$ without bootstrapping. Fig. 3
shows the ratio of the optimal $\hat{\sigma}^{2}$ and the optimal
$\hat{s}^{2}$. Panels (a) and (b) show the results when the $\ell_{1}$-ratio
$\gamma$ is fixed and panel (c) shows when the $\ell_{1}$-ratio is optimized.
As expected from Proposition 3, the variance is reduced compared to the case
without bootstrapping. The magnitude of the reduction is larger when the
$\ell_{1}$-ratio is fixed and close to $1$ (LASSO [21] case). In particular,
the largest improvement is obtained when $\rho$ is large (less sparse) and the
measurement ratio $\alpha$ is small. The improvement is minor when
$\ell_{1}$-ratio is optimized. This suggests that using the bootstrap average
is effective when the actual data generation process is inconsistent with the
sparsity assumption of the regularization and the data size is insufficient.
However, when the data size is too small, meaningful improvement cannot be
obtained.

(a)

(b)

(c)

Figure 4: The optimal regularization parameter $\log_{10}\lambda$ for the
bootstrap averaged unbiased estimator are plotted against the sparsity of the
true signal $\rho$ and the measurement ratio $\alpha=M/N$ as a heat map. The
measurement noise is set as $\Delta=0.15$.

(a)

(b)

(c)

Figure 5: The typical number of unique samples in each bootstrap sample
$\alpha(1-e^{-\mu_{B}^{\star}})$ scaled by $N$ is visualized. In the purple
region, $\alpha(1-e^{-\mu_{B}^{\star}})<1$, and in the yellow region,
$\alpha(1-e^{-\mu_{B}^{\star}})\geq 1$. The red dashed line shows $\alpha=1$.
$\Delta=0.15$.

###  V-C Phase transition to an ensemble of interpolators

Fig. 4 shows the optimal regularization strength $\lambda^{\star}$ for
$\hat{\sigma}^{2}$. Panels (a) and (b) show the results when the
$\ell_{1}$-ratio is fixed, and panel (c) shows when the $\ell_{1}$-ratio is
optimized. In all cases, it is clear that a phase transition has occurred in
which $\lambda^{\ast}$ drops to an infinitesimally small value (although for
numerical reasons, $\lambda$ is constrained to exceed $10^{-7}$) as the
measurement ratio $\alpha$ decreases or $\rho$ increases. Moreover, Fig. 5
shows $\alpha(1-\mu_{B}^{\star})$, the typical number of unique data points in
each bootstrap sample scaled by $N$: $\lim_{M,N\to\infty}(M\mu_{B})/N$. From
Fig. 5, it is clear that the typical number of unique data points is always
smaller than $1$ in the region where $\lambda^{\star}\simeq+0$. This holds,
even if the measurement ratio $\alpha>1$. Thus, when
$\lambda^{\star}\simeq+0$, the elastic net estimator in (8) becomes the
minimum elastic net norm estimator as

$\hat{\mathbf{\boldsymbol{w}}}^{\ast}=\min_{\mathbf{\boldsymbol{w}}\in\mathbb{R}^{N}}\sum_{i=1}^{N}\gamma|w_{i}|+\frac{1-\gamma}{2}w_{i}^{2},\,{\rm
subject\,to}\\\
\quad\mbox{1}\mbox{l}[c_{\mu}>0](y_{\mu}-\mathbf{\boldsymbol{x}}_{\mu}^{\top}{w})=0,\,\mu=1,2,\dots,M,$
(14)

which is commonly known as minimum elastic net norm _interpolator_ in machine
learning [15, 16, 17]. These suggest that when elastic net regularization
cannot determine an appropriate sparse structure of
$\mathbf{\boldsymbol{w}}_{0}$, it is better to use an over-parameterized
setting in which the number of unique data points of each bootstrap data is
smaller than the dimension of $\mathbf{\boldsymbol{w}}_{0}$ and use an
ensemble of interpolators.

##  VI Summary and discussion

In this study, we investigated the behavior of the bootstrap-averaged unbiased
estimator of GAMP using AMPR and its SE. We found that the bootstrap averaging
procedure can effectively reduce the variance of the unbiased estimator when
the actual data generation process is inconsistent with the sparsity
assumption of the regularization and the data size is insufficient. We also
found a phase transition where the regularization strength drops to
infinitesimally small values by decreasing the measurement ratio $\alpha$ or
increasing $\rho$.

Although increasing the variance of weak learners is key to the success of
ensemble learning [22, 23, 24], the phase transition to an ensemble of
interpolators may be unexpected. Investigating whether similar phase
transitions occur in other more sophisticated machine-learning models, such as
neural networks, would be an interesting future direction.

We also remark that the landscape in $(\mu_{B},\lambda)$-space may not be
convex in general. Thus, the uniqueness of the optimizer should be
investigated in the future, paying attention to the relationship with the
implicit regularization [25, 26, 27].

On the technical side, the key to this study was the precise performance
characterization of the averaged estimator by AMPR, which was developed by
combining GAMP and the replica method of statistical physics [28, 29]. Such a
combination of approximate inference algorithms and the replica method has
been used to develop approximate computation algorithms [13, 30, 14, 31, 32,
33, 34] and has not been applied to precise performance analysis of ensemble
methods. It would be an interesting direction to try similar performance
analysis for other bootstrap methods or ensemble learning.

## Acknowledgement

This study was supported by JSPS KAKENHI Grant Number 21K21310 and 17H00764.

## References

  * [1] S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” in _2011 IEEE International Symposium on Information Theory Proceedings_. IEEE, 2011, pp. 2168–2172.
  * [2] A. Javanmard and A. Montanari, “State evolution for general approximate message passing algorithms, with applications to spatial coupling,” _Information and Inference: A Journal of the IMA_ , vol. 2, no. 2, pp. 115–144, 2013.
  * [3] P. Sur and E. J. Candès, “A modern maximum-likelihood theory for high-dimensional logistic regression,” _Proceedings of the National Academy of Sciences_ , vol. 116, no. 29, pp. 14 516–14 525, 2019.
  * [4] A. Javanmard and A. Montanari, “Hypothesis testing in high-dimensional regression under the gaussian random design model: Asymptotic theory,” _IEEE Transactions on Information Theory_ , vol. 60, no. 10, pp. 6522–6554, 2014.
  * [5] T. Takahashi and Y. Kabashima, “A statistical mechanics approach to de-biasing and uncertainty estimation in lasso for random measurements,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2018, no. 7, p. 073405, 2018.
  * [6] P. Sur, Y. Chen, and E. J. Candès, “The likelihood ratio test in high-dimensional logistic regression is asymptotically a rescaled chi-square,” _Probability theory and related fields_ , vol. 175, no. 1, pp. 487–558, 2019.
  * [7] A. Sakata and T. Obuchi, “Perfect reconstruction of sparse signals with piecewise continuous nonconvex penalties and nonconvexity control,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2021, no. 9, p. 093401, 2021.
  * [8] L. Zheng, A. Maleki, H. Weng, X. Wang, and T. Long, “Does $\ell_{p}$-minimization outperform $\ell_{1}$-minimization?” _IEEE Transactions on Information Theory_ , vol. 63, no. 11, pp. 6896–6935, 2017.
  * [9] B. Efron, “Bootstrap Methods: Another Look at the Jackknife,” _The Annals of Statistics_ , vol. 7, no. 1, pp. 1 – 26, 1979. [Online]. Available: https://doi.org/10.1214/aos/1176344552
  * [10] Z.-H. Zhou, _Ensemble methods: foundations and algorithms_. CRC press, 2012.
  * [11] K. P. Murphy, _Probabilistic machine learning: an introduction_. MIT press, 2022.
  * [12] H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” _Journal of the royal statistical society: series B (statistical methodology)_ , vol. 67, no. 2, pp. 301–320, 2005.
  * [13] T. Obuchi and Y. Kabashima, “Semi-analytic resampling in lasso,” _The Journal of Machine Learning Research_ , vol. 20, no. 1, pp. 2544–2576, 2019.
  * [14] T. Takahashi and Y. Kabashima, “Semi-analytic approximate stability selection for correlated data in generalized linear models,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2020, no. 9, p. 093402, 2020\.
  * [15] V. Muthukumar, K. Vodrahalli, V. Subramanian, and A. Sahai, “Harmless interpolation of noisy data in regression,” _IEEE Journal on Selected Areas in Information Theory_ , vol. 1, no. 1, pp. 67–83, 2020.
  * [16] P. L. Bartlett, P. M. Long, G. Lugosi, and A. Tsigler, “Benign overfitting in linear regression,” _Proceedings of the National Academy of Sciences_ , vol. 117, no. 48, pp. 30 063–30 070, 2020.
  * [17] T. Hastie, A. Montanari, S. Rosset, and R. J. Tibshirani, “Surprises in high-dimensional ridgeless least squares interpolation,” _The Annals of Statistics_ , vol. 50, no. 2, pp. 949–986, 2022.
  * [18] P. K. Mogensen and A. N. Riseth, “Optim: A mathematical optimization package for Julia,” _Journal of Open Source Software_ , vol. 3, no. 24, p. 615, 2018.
  * [19] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [20] T. Obuchi and Y. Kabashima, “Cross validation in lasso and its acceleration,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2016, no. 5, p. 053304, 2016.
  * [21] R. Tibshirani, “Regression shrinkage and selection via the lasso: a retrospective,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 73, no. 3, pp. 273–282, 2011.
  * [22] A. Krogh and P. Sollich, “Statistical mechanics of ensemble learning,” _Physical Review E_ , vol. 55, no. 1, p. 811, 1997.
  * [23] P. Sollich and A. Krogh, “Learning with ensembles: How overfitting can be useful,” _Advances in neural information processing systems_ , vol. 8, 1995\.
  * [24] L. I. Kuncheva and C. J. Whitaker, “Measures of diversity in classifier ensembles and their relationship with the ensemble accuracy,” _Machine learning_ , vol. 51, no. 2, pp. 181–207, 2003.
  * [25] D. LeJeune, H. Javadi, and R. Baraniuk, “The implicit regularization of ordinary least squares ensembles,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2020, pp. 3525–3535.
  * [26] J.-H. Du, P. Patil, and A. K. Kuchibhotla, “Subsample ridge ensembles: Equivalences and generalized cross-validation,” _arXiv preprint arXiv:2304.13016_ , 2023.
  * [27] T. Yao, D. LeJeune, H. Javadi, R. G. Baraniuk, and G. I. Allen, “Minipatch learning as implicit ridge-like regularization,” in _2021 IEEE International Conference on Big Data and Smart Computing (BigComp)_. IEEE, 2021, pp. 65–68.
  * [28] M. Mézard and A. Montanari, _Information, physics, and computation_. Oxford University Press, 2009\.
  * [29] A. Montanari and S. Sen, “A short tutorial on mean-field spin glass techniques for non-physicists,” _arXiv preprint arXiv:2204.02909_ , 2022.
  * [30] T. Takahashi and Y. Kabashima, “Replicated vector approximate message passing for resampling problem,” _arXiv preprint arXiv:1905.09545_ , 2019.
  * [31] D. Malzahn and M. Opper, “An approximate analytical approach to resampling averages,” _Journal of Machine Learning Research_ , vol. 4, no. Dec, pp. 1151–1173, 2003.
  * [32] ——, “A statistical mechanics approach to approximate analytical bootstrap averages,” in _Advances in Neural Information Processing Systems_ , 2003, pp. 343–350.
  * [33] ——, “Approximate analytical bootstrap averages for support vector classifiers,” in _Advances in Neural Information Processing Systems_ , S. Thrun, L. Saul, and B. Schölkopf, Eds., vol. 16. MIT Press, 2003. [Online]. Available: https://proceedings.neurips.cc/paper/2003/file/2c6ae45a3e88aee548c0714fad7f8269-Paper.pdf
  * [34] M. Opper, “An approximate inference approach for the pca reconstruction error,” _Advances in Neural Information Processing Systems_ , vol. 18, 2005\.


YITP-20-147, IPMU20-0119

# Path Integral Optimization from Hartle-Hawking Wave Function

Jan Borucha, Pawel Caputaa and Tadashi Takayanagib,c,d aFaculty of Physics,
University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland.  bYukawa
Institute for Theoretical Physics, Kyoto University,
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan  cInamori Research
Institute for Science, 620 Suiginya-cho, Shimogyo-ku, Kyoto 600-8411 Japan
dKavli Institute for the Physics and Mathematics of the Universe,
University of Tokyo, Kashiwa, Chiba 277-8582, Japan

###### Abstract

We propose a gravity dual description of the path integral optimization in
conformal field theories arXiv:1703.00456, using Hartle-Hawking wave functions
in anti-de Sitter spacetime. We show that the maximization of the Hartle-
Hawking wave function is equivalent to the path integral optimization
procedure. Namely, the variation of the wave function leads to a constraint,
equivalent to the Neumann boundary condition on a bulk slice, whose classical
solutions reproduce metrics from the path integral optimization in conformal
field theories. After taking the boundary limit of the semi-classical Hartle-
Hawking wave function, we reproduce the path integral complexity action in two
dimensions as well as its higher- and lower-dimensional generalizations. We
also discuss an emergence of holographic time from conformal field theory path
integrals.

Introduction$-$ The AdS/CFT correspondence Ma  provides us with a surprising
relation between gravity and quantum many-body systems. Nevertheless, the
fundamental mechanism of how it works so well is still not understood. This
problem is one of the main obstacles when we try to extend the holographic
duality to more general spacetimes, including realistic universe. One
interesting idea, pioneered in Swingle , toward uncovering the basic mechanism
behind AdS/CFT is to relate the emergent AdS geometries to tensor networks
such as Multi-scale Entanglement Renormalization Ansatz MERA ; TNR ; cMERA  or
more general ones Qi:2013caa ; HAPPY ; HQ , realizing emergence of spacetimes
from quantum entanglement Ra . In particular, this tensor network
interpretation beautifully explains the geometric calculation of entanglement
entropy in AdS/CFT RT ; HRT . Refer to e.g. Beny ; NRT ; MT ; Cz ;
Czech:2015kbp ; Milsted:2018yur ; Milsted:2018san ; Jahn:2020ukq  for further
developments in this direction. However, these tensor network approaches have
been limited to toy models on discrete lattices and precise relations between
them and genuine AdS/CFT is not clear. See also recent attempts directly from
AdS/CFT Takayanagi:2018pml ; VanRaamsdonk:2018zws ; Bao:2018pvs ; Bao:2019fpq
.

On the other hand, the path integral optimization Caputa:2017urj ; MTW , that
we will now review, provides a useful framework that describes tensor networks
for conformal field theories (CFTs) in terms of path integrals. We take the
Euclidean $\mathbb{R}^{2}$ coordinates $(\tau,x)$ and denote all fields in the
CFT by $\Phi(\tau,x)$. The ground state wave functional
$\Psi_{\text{CFT}}[\Phi(x)]$ at the time slice $\tau=0$ is defined by the path
integral

$\displaystyle\Psi_{\text{CFT}}[\Phi(x)]$ (1)
$\displaystyle=\int\prod_{-\infty<\tau\leq
0,x}[D\tilde{\Phi}(\tau,x)]e^{-S_{\text{CFT}}[\tilde{\Phi}]}\delta(\tilde{\Phi}(0,x)-\Phi(x)),$

where $S_{\text{CFT}}$ is the action of the 2D CFT.

In the path integral optimization, we deform the metric of our 2D space on
which we perform the path integral as follows:

$\displaystyle ds^{2}=e^{2\phi(\tau,x)}(d\tau^{2}+dx^{2}).$ (2)

We take $e^{2\phi(\tau,x)}=1/\epsilon^{2}$ for the flat metric of
$\mathbb{R}^{2}$ used in the original path integral that computes
$\Psi_{\text{CFT}}[\Phi(x)]$, where $\epsilon$ is a UV regularization scale
(i.e., lattice constant) when we discretize path integrals of quantum fields
into those on a lattice. The curved space metric is interpreted as a choice of
discretization such that there is a single lattice site per a unit area.

Let us write the wave functional obtained from the path integral on the curved
space (2) as $\Psi_{\text{CFT}}^{\phi}[\Phi(x)]$. If we impose the boundary
condition

$e^{2\phi(0,x)}={1\over\epsilon^{2}}\equiv e^{2\phi_{0}},$ (3)

this wave functional is proportional to the one
$\Psi_{\text{CFT}}[\Phi(x)](=\Psi_{\text{CFT}}^{\phi_{0}}[\Phi(x)])$ for the
flat space $\mathbb{R}^{2}$ since the CFT is invariant up to the Weyl anomaly

$\displaystyle\Psi_{\text{CFT}}^{\phi}[\Phi(x)]=e^{S_{L}[\phi]-S_{L}[\phi_{0}]}\cdot\Psi_{\text{CFT}}^{\phi_{0}}[\Phi(x)].$
(4)

Here $S_{L}$ is the Liouville action Po

$\displaystyle S_{L}[\phi]\\!=\\!{c\over
24\pi}\int^{\infty}_{-\infty}\\!dx\\!\int^{0}_{-\infty}\\!d\tau\\!\left[(\partial_{x}\phi)^{2}+(\partial_{\tau}\phi)^{2}+\mu
e^{2\phi}\right],$ (5)

and $c$ is the central charge of the 2D CFT. The assumption of the
discretization, that one unit area of the metric (2) has a single lattice
site, fixes the values of $\mu$ to $\mu=1$ Caputa:2017urj . Nevertheless, it
is useful to keep this cosmological constant parameter for later purpose.

Relation (4) guarantees that the quantum state is still the same CFT vacuum
$|0\rangle$ for any choice of the metric (2) as long as the boundary condition
(3) is satisfied. Since the potential term in (5) originates from the UV
divergence and we consider $S_{L}$ as a bare action, it should dominate over
the kinetic term when we take the UV cutoff to infinity. This is realized when

$\displaystyle(\partial_{i}\phi)^{2}\ll e^{2\phi}\ \ \ (i=x,\tau).$ (6)

The idea of path integral optimization is to coarse grain the discretization
as much as possible, which makes computational costs minimal, while keeping
the correct answer to the final wave functional. This path integral
optimization is performed by minimizing the functional $S_{L}[\phi]$
Caputa:2017urj . This is because we want to minimize the overall factor of the
wave functional, which is proportional to $e^{S_{L}[\phi]}$ as in (4). Even
though the overall factor does not affect physical quantities in quantum
mechanics, this estimates the number of repetitions of numerical integrals
when we discretizes the required path integral into lattice calculations whose
regularization is specified by the metric (2). Therefore, the Liouville action
$S_{L}$ (at $\mu=1$) was identified with a measure of computational complexity
Susskind , called the path integral complexity Caputa:2017urj  (refer to
Czech:2017ryf ; Caputa:2018kdj ; Camargo:2019isp ; Erdmenger:2020sup  for
connections to circuit complexity). The minimization procedure picks up the
most efficient discretization of path integral which leads to the correct
vacuum state. This method was generalized to various CFT setups in
Sato:2019kik ; Caputa:2020mgb ; Bhattacharyya:2018wym  and used to compute
entanglement of purification in 2D CFTs Caputa:2018xuf , which was recently
verified numerically in Camargo:2020yfv .

For the vacuum, the minimization is performed by solving the Liouville
equation $(\partial_{x}^{2}+\partial^{2}_{\tau})\phi=\mu e^{2\phi}$ with the
boundary condition (3), leading to the solution

$\displaystyle e^{2\phi(\tau,x)}={1\over(\sqrt{\mu}\tau-\epsilon)^{2}},\ \ \
(-\infty<\tau\leq 0).$ (7)

The solution at $\mu=1$ is the genuine optimized one, which means the
minimization of the original Liouville action. The choice $\mu<1$ may be
regarded as partially optimized solution where the UV cutoff length scale is
taken to be larger, while the choice $\mu>1$ is not allowed as this
corresponds to fine-graining the cutoff scale more than the current lattice
spacing.

The observation that (7) coincides with the time slice of a three dimensional
AdS (AdS3), leads to the main implication that the path integral optimization
can explain an emergent AdS geometry purely from CFT Caputa:2017urj . The
discretized path integral takes the form of a (nonunitary) tensor network and
its relation to AdS can be regarded as a path integral version of the
conjectured interpretation of AdS/CFT as tensor networks. However, a direct
connection between the path integral optimization and AdS/CFT has remained an
open problem.

Another subtle issue is that in the solution (7), we find
$(\partial_{i}\phi)^{2}$ and $e^{2\phi}$ are of the same order, which does not
satisfy the criterion (6). This suggests that the path integral optimization
using the Liouville action is qualitative and can have finite cutoff
corrections.

Moreover, it has not been clear how to promote the classical Liouville theory
equivalent to the above path integral optimization to a quantum Liouville
theory. Indeed, it was found in Caputa:2017urj  that to properly reproduce the
correct gravity metric dual to a primary state with $1/c$ corrections, we need
to replace a classical Liouville theoretic result with a quantum Liouville
theoretic one. This is because the path-integration
$\int[D\phi]e^{S_{L}[\phi]}$ in (4) does not make sense as it is not bounded
from below. Instead the quantum Liouville theory is defined by the path
integral $\int[D\phi]e^{-S_{L}[\phi]}$. In other words, we cannot get the
minimization as a saddle point approximation of path integrals and a
derivation of path integral optimization from AdS/CFT remained a challenge.

In this paper, we would like to resolve these important issues by introducing
a gravity dual description in terms of a Hartle-Hawking wave function which
evolves from the AdS boundary. This corresponds to the gravity action in the
shaded region in Fig.1, assuming the Euclidean Poincare AdSd+1 geometry

$\displaystyle ds^{2}={dz^{2}+d\tau^{2}+\sum_{i=1}^{d-1}dx_{i}^{2}\over
z^{2}}.$ (8)

Figure 1: The on-shell action in the colored region $M$ provides a gravity
evaluation of the Hartle-Hawking wave function.

Hartle-Hawking wave function with Boundary$-$ Consider a Hartle-Hawking wave
function Hartle:1983ai  in an AdSd+1, denoted by $\Psi_{HH}[g_{ab}]$, which is
a functional of the metric $g_{ab}$ on a surface $Q$. Respecting the timelike
boundary in AdS, we impose an initial condition on the AdS boundary $\Sigma$
given by $z=\epsilon$ and $\tau<0$. Then we consider a path integral of
Euclidean gravity from this asymptotic boundary to the surface $Q$ which
extends from $z=\epsilon$ and $\tau=0$ toward the bulk as depicted in Fig.1.
Focusing on translational invariant setups for simplicity, we assume the
diagonal form metric on $Q$

$\displaystyle ds^{2}=e^{2\phi}(dw^{2}+\sum_{i=1}^{d-1}dx_{i}^{2}),$ (9)

where $w$ is a function of $\tau$. This way, we can write the Hartle-Hawking
wave function as $\Psi_{HH}[\phi(w,x)]$, defined by

$\displaystyle\Psi_{HH}[\phi]=\int[Dg_{\mu\nu}]e^{-I_{G}[g]}\delta(g_{ab}|_{Q}-e^{2\phi}\delta_{ab}),$
(10)

where $I_{G}$ is the $d+1$-dimensional gravity action

$\displaystyle I_{G}\\!=\\!-{1\over 16\pi
G_{N}}\\!\int_{M}\\!\sqrt{g}(R\\!-\\!2\Lambda)\\!-\\!{1\over 8\pi
G_{N}}\\!\int_{Q\cup\Sigma}\\!\sqrt{h}K.~{}$ (11)

We implicitly imposed a boundary condition on $\Sigma$. Even though we choose
that of pure AdS dual to the CFT vacuum, in principle, we can consider
$\Psi_{HH}[\phi(w,x)]$ corresponding to excited states of a CFT (see below).

Finally, we propose to identify the metric (9) with (2) (in $d=2$ case) and,
after setting $w=\tau$, we argue that the optimization procedure corresponds
to the maximization of $\Psi_{HH}[\phi]$ with respect to $\phi$. This
maximization can be understood naturally when we consider an evaluation of
correlation function as

$\displaystyle\langle
O_{1}O_{2}\ldots\rangle=\int[D\phi]|\Psi_{HH}[\phi]|^{2}O_{1}O_{2}\ldots,$
(12)

by applying the saddle point approximation. In this way, the maximization of
Hartle-Hawking wave function works well even in the presence of quantum
fluctuations of $\phi$, as opposed to the minimization of $e^{S_{L}[\phi]}$.
Indeed, below, we will show that $\Psi_{HH}[\phi]$ is proportional to
$e^{-S_{L}[\phi]}$ up to finite cutoff corrections.

It will also be useful to add a tension term on the brane $Q$ as in the
AdS/BCFT Ta  (we assume $T<0$ below)

$\displaystyle I_{T}[h]={T\over 8\pi G_{N}}\int_{Q}\sqrt{h}$ (13)

and define one parameter family of deformed Hartle-Hawking wave functions as
follows:

$\displaystyle\Psi^{(T)}_{HH}[\phi]\\!=\\!\int\\![Dg_{\mu\nu}\\!]e^{-I_{G}[g]-I_{T}[e^{2\phi}]}\,\delta(g_{ab}|_{Q}\\!-\\!e^{2\phi}\delta_{ab}).~{}~{}$
(14)

The standard Hartle-Hawking wave functions are obtained by setting $T=0$. We
can regard $T$ as a chemical potential or Legendre transformation for the area
of $Q$, acting on Hartle-Hawking wave functions. As we will see, the tension
term plays the role of the cosmological constant term in the Liouville action.
More importantly, since the maximum of $\Psi^{(T)}_{HH}[\phi]$ corresponds to
a family of surfaces $Q$ in AdS parameterized by the tension $T$, we will
observe that $T$ plays the role of an emergent holographic time.

Evaluation of $\Psi^{(T)}_{HH}[\phi]$ $-$ Let us evaluate
$\Psi^{(T)}_{HH}[\phi]$ using semi-classical approximation i.e., as an on-
shell gravity action

$\displaystyle\Psi^{(T)}_{HH}[\phi]\simeq e^{-I_{G}[g]-I_{T}[h]}|_{\mbox{on-
shell}}.$ (15)

For simplicity we assume that the metric (9) has a translational invariance in
the $x$ direction. For a given choice of such a metric $e^{2\phi}$, we find a
surface $Q$ specified by the profile

$\displaystyle z=f(\tau).$ (16)

Then the semi-classical evaluation of (10) will give the value of Hartle-
Hawking wave function.

Note that the above construction assumes that a gravity solution which
calculates the Hartle-Hawking wave function, is given by a subreagion in a
Poincare AdS. The metric (9) obtained in this way covers all possible metrics
on $Q$ for $d=2$ with the condition (3) because all solutions to the vacuum
Einstein equation are locally equivalent to the Poincare AdS3. However, for
$d>2$, the above construction covers only a part of the metrics on the
$d$-dimensional surface $Q$. Nevertheless, this ansatz includes a class of
metric we want, as we will see below. It is also useful to note that we can
extend our targets to general metrics by directly solving the Einstein
equation.

The $w$ coordinate and the metric in (9) is found as

$\displaystyle e^{\phi}=f^{-1},\ \ \ {d\tau\over dw}=\sqrt{1-\dot{f}^{2}},$
(17)

where $\dot{f}$ means $\partial_{w}f$. We set $w=0$ at $\tau=0$ and thus we
have $e^{-\phi}=\epsilon$ at $w=0$. Then, the on-shell action on $M$, defined
by $\epsilon<z<f(\tau)$ (see $M$ in Fig.1) is evaluated in terms of the
coordinate $w$ as follows:

$\displaystyle I_{G}+I_{T}$ $\displaystyle=$
$\displaystyle-{(d-1)V_{x}L_{\tau}\over 8\pi
G_{N}\epsilon^{d}}+{(d-1)V_{x}\over 8\pi G_{N}}\int dwe^{d\phi}G(\dot{\phi})$
(18) $\displaystyle-{V_{x}\over 8\pi
G_{N}}\left[e^{(d-1)\phi}\arcsin(\dot{\phi}e^{-\phi})\right]^{0}_{-\infty},$

where $V_{x}$ and $L_{\tau}$ are infinite volumes in $x$ and $\tau$ directions
and $G$ is the following function bounded from below

$\displaystyle
G(\dot{\phi})\\!=\\!\sqrt{1-e^{-2\phi}{\dot{\phi}}^{2}}\\!+\\!\dot{\phi}e^{-\phi}\arcsin(\dot{\phi}e^{-\phi})\\!+\\!{T\over
d-1}.$ (19)

When we neglect the finite cutoff corrections assuming (6), we can approximate
the on-shell action (18) as a quadratic action of $\phi$. In $d=2$, this
expansion yields

$\displaystyle I_{G}+I_{T}\simeq$ $\displaystyle{c\over 12\pi}\int
dxdw\left[{1\over 2}\dot{\phi}^{2}+(1+T)e^{2\phi}-\epsilon^{-2}\right]$ (20)
$\displaystyle-{c\over 12\pi}\int dx{\theta_{0}\over\epsilon},$

where $\theta_{0}$ is the value of $\arcsin(\dot{\phi}e^{-\phi})$ at $w=0$. We
can cancel $\theta_{0}$ dependence by adding the corner term Hayward:1993my
localized on $\Sigma\cap Q$. This reproduces the ”Path Integral Complexity”
action Caputa:2017urj  $I[\phi,\phi_{0}]=S_{L}[\phi]-S_{L}[\phi_{0}]$ with the
correct coefficient of the kinetic term (remember we assumed
$\partial_{x}\phi=0$). Note that the above gravity computation (18) gives the
correct finite cutoff corrections to the Liouville action and thus provides a
full answer to the path integral optimization. The same is true in higher
dimensions. Notice that, unusual from the perspective of complexity,
properties $I[g_{1},g_{2}]=-I[g_{2},g_{1}]$ and
$I[g_{1},g_{2}]+I[g_{2},g_{3}]=I[g_{1},g_{3}]$ become manifest from the
gravity action with boundaries.

Solutions$-$ Now we would like to maximize the Hartle-Hawking wave function
(15). This is performed by taking a variation of the on-shell action (18) with
respect to $\phi$, leading to

$\displaystyle{e^{-2\phi}(\ddot{\phi}+(d-1)\dot{\phi}^{2})-d\over\sqrt{1-e^{-2\phi}{\dot{\phi}}^{2}}}={d\over
d-1}T.$ (21)

By imposing the boundary condition

$\displaystyle e^{2\phi}|_{w=0}={1\over\epsilon^{2}},$ (22)

we obtain the solution to (21) when $T<0$ as follows:

$\displaystyle
e^{2\phi}={1\over\left(\sqrt{1-{T^{2}\over(d-1)^{2}}}w-\epsilon\right)^{2}}.$
(23)

This corresponds to the following surface in (8):

$\displaystyle z=\epsilon+\sqrt{1-{T^{2}\over(d-1)^{2}}}{(d-1)\over T}\tau.$
(24)

For the solution (23), the on-shell action is evaluated as (see more in
Boruch:2021hqs )

$\displaystyle I_{G}+I_{T}=-{(d-1)V_{x}L_{\tau}\over 8\pi G_{N}\epsilon^{d}}.$
(25)

Previously we observed that the maximization of the Hartle-Hawking wave
function is given by minimizing the integral of (19) i.e., the Liouville
action plus finite cutoff corrections. On the other hand, the original path
integral optimization is proposed to be simply given by the minimization of
Liouville action. However, it is possible that this apparent deviation arises
because the regularization scheme is different between the CFT and gravity
formation, though they are actually equivalent. Though we do not have any
definite argument for this, the fact that both give the same profile of
optimized solution may imply this equivalence. Indeed if we set

$\mu=1-{T^{2}\over(d-1)^{2}},$ (26)

then the solution (23) is equal to that from the path integral optimization
(7). Remember that changing $\mu$ from $\mu=0$ to $\mu=1$ means that we
gradually increase the amount of optimization. In the gravity dual, this
corresponds to changing the tension from $T=-(d-1)$ to $T=0$ which tilts the
surface $Q$ from the asymptotic boundary $\Sigma$ to the time slice $\tau=0$.

Note also that the equation (21) for $\phi$ is equivalent to imposing the
Neumann boundary condition on $Q$

$\displaystyle K_{ab}-Kh_{ab}=-Th_{ab},$ (27)

which is imposed in the AdS/BCFT construction Ta .

Let us finally stress that, in higher dimensions $d>2$ there has not been a
complete formulation of path integral optimization known till now as we do not
know a higher-dimensional counterpart of the Liouville action. Remarkably, our
approach using Hartle-Hawking wave function gives the full answer to this
question for CFTs with gravity duals.

Excited states$-$ Furthermore, we consider a family of Euclidean BTZ-type
metrics in three dimensions

$\displaystyle ds^{2}=(r^{2}-r^{2}_{h})d\tau^{2}+{dr^{2}\over
r^{2}-r^{2}_{h}}+r^{2}dx^{2},$ (28)

where $r^{2}_{h}=M-1$ can be positive or negative depending on the mass of the
excitation. We can repeat our analysis of $\Psi^{T}_{HH}$ for region ${1\over
f(\tau)}\leq r\leq{1\over\epsilon}$, where $r=1/f(\tau)$ and $r=1/\epsilon$
describe the surface $Q$ and the asymptotic boundary $\Sigma$, respectively.
Our action has the same form as (18) for $d=2$ and

$\displaystyle G(\dot{\phi})$ $\displaystyle=$
$\displaystyle\sqrt{1-e^{-2\phi}\left(\dot{\phi}^{2}+r^{2}_{h}\right)}+T$ (29)
$\displaystyle+$
$\displaystyle\dot{\phi}e^{-\phi}\arcsin\left({\dot{\phi}e^{-\phi}\over\sqrt{1-r^{2}_{h}e^{-2\phi}}}\right).$

Variation with respect to $\phi$ yields $K_{Q}=2T$ and is again equivalent to
the Neumann condition (27). For negative tension, we can solve it by

$\displaystyle
e^{2\phi}={r^{2}_{h}\over\left(1-T^{2}\right)\sin^{2}\left(r_{h}(w-c_{1})\right)}.$
(30)

This family of solutions precisely matches those in the path integral
optimization Caputa:2017urj  via the identification (26). For
$r_{h}=(2\pi)/\beta$ and $w$ on the strip, we reproduce our optimal geometry
for the thermofield double (TFD) state dual to the time slice of Einstein-
Rosen bridge MaE . For $r^{2}_{h}=-(1-M)$ we reproduce excited states from the
optimization for primary operators in 2D CFT i.e., conical singularity
geometries, including the finite size vacuum. In all these examples we choose
$c_{1}$ such that (3) is fulfilled at each boundary. Moreover, we can verify
(either by explicit computation or using the Wheeler-DeWitt equation) that our
solutions (7), (23) and (30) have constant negative curvature that depends on
$\mu$ (or $T$ via (26)).

Last but not least, we can test our prescription in the context of JT gravity
dual to the SYK model SYKa ; SYKb ; SYKc ; SYKd ; Jensen:2016pah . In this
case, it turns out that it is advantageous to introduce the tension on Q by
coupling to the dilaton $\Phi$

$\displaystyle I_{JT}+I_{T_{\Phi}}$ $\displaystyle=$
$\displaystyle-\left[\int_{M}\sqrt{g}\Phi(R+2)+2\int_{\partial M}\sqrt{h}\Phi
K\right]$ (31) $\displaystyle-$
$\displaystyle\Phi_{0}\chi(M)+2T_{\Phi}\int_{Q}\sqrt{h}\Phi,$

with $\chi(M)$ being the Euler characteristic of our region $M$. As an
example, we can take the analogous 2D solution

$ds^{2}=(r^{2}-r^{2}_{h})d\tau^{2}+{dr^{2}\over
r^{2}-r^{2}_{h}},\quad\Phi=\Phi_{b}r.$ (32)

Defining M bounded by $\Sigma$ at $r=r_{0}\to\infty$ and Q by $r=1/f(\tau)$,
with induced metric

$ds^{2}=e^{2\phi}dw^{2},\qquad e^{\phi}=f^{-1},$ (33)

we can derive a JT analog of (29), show that the saddle point equation is the
Neumann b.c. and find that its solution is given by (30) with $T\to T_{\Phi}$.
Moreover, in the UV limit of small $\dot{\phi}$ our action reproduces the
effective Schwarzian description of the SYK with the symmetry breaking term.
This confirms the validity of our approach in all dimensions. We performed
analogous studies for higher dimensional black holes as well as examples of
Lorentzian spacetimes and details will be presented in Boruch:2021hqs .

Conclusions and Discussion$-$ In this paper, we showed that the path integral
optimization corresponds to the maximization of Hartle-Hawking wave function,
which is a functional of the metric (9) on a surface $Q$:
$\mbox{Max}_{\phi}[\Psi_{HH}[\phi]]$. This Hartle-Hawking wave function with
the boundary condition (22), describes an evolution from an initial condition
set by the AdS boundary, dual to the target CFT state for which we perform the
path integral optimization. An important requirement is that the surface $Q$
ends on the AdS boundary and this gives the boundary condition (22). Owing to
this requirement, we can calculate CFT quantities such as correlations
functions from an inner product of Hartle-Hawking wave functional as in (12),
whose saddle point approximation gives the maximization of $\Psi_{HH}[\phi]$.

Furthermore, we generalize our correspondence to non-trivial parameter: $\mu$
in Liouville theory and tension $T$ in gravity, related by (26). In the path
integral optimization $\mu$ controls the scale up to which we perform the
coarse-graining and this optimization procedure is maximized at $\mu=1$. On
the gravity side, this scale is fixed by the tension term (13), which plays
the role of a chemical potential to the area of the surface $Q$. Even though
the original Hartle-Hawking wave function does not have any time-dependence,
”time” emerges by considering the $T$-dependent one: $\Psi^{(T)}_{HH}[\phi]$,
where the tension plays a role of external field. From the CFT side, time
emerges as the scale $\mu$ of the optimization, related to $T$ via (26).
Indeed, using the optimized solution (23) or (30), we can write the full
AdSd+1 space as follows:

$\displaystyle ds^{2}={d\mu^{2}\over
4\mu^{2}(1-\mu)}+e^{2\phi}(dw^{2}+dx^{2}_{i}).$ (34)

Note that this foliation is a special case of the York time York:1972sj
(refer to Belin:2018bpg  for an interesting interpretation of York time from
AdS/CFT). It will be a very important future direction to derive the genuine
AdS/CFT itself by starting from the purely CFT analysis of path integral
optimization. We believe that this emergent time observation provides us with
an important clue in this direction.

Acknowledgements $-$ We are grateful to Sumit Das, Diptarka Das, Dongsheng Ge,
Jacek Jezierski, Masamichi Miyaji, Onkar Parrikar for useful discussions. TT
is supported by Grant-in-Aid for JSPS Fellows No. 19F19813. TT is supported by
the Simons Foundation through the “It from Qubit” collaboration. TT is
supported by Inamori Research Institute for Science and World Premier
International Research Center Initiative (WPI Initiative) from the Japan
Ministry of Education, Culture, Sports, Science and Technology (MEXT). TT is
supported by JSPS Grant-in-Aid for Scientific Research (A) No. 16H02182. TT is
also supported by JSPS Grant-in-Aid for Challenging Research (Exploratory)
18K18766. PC and JB are supported by NAWA “Polish Returns 2019” and NCN Sonata
Bis 9 grants.

## References

  * (1) J. M. Maldacena, Adv. Theor. Math. Phys. 2 (1998) 231 [Int. J. Theor. Phys. 38 (1999) 1113] [arXiv:hep-th/9711200];
  * (2) B. Swingle, “Entanglement Renormalization and Holography,” Phys. Rev. D 86, 065007 (2012), arXiv:0905.1317 [cond-mat.str-el].
  * (3) G. Vidal,“A class of quantum many-body states that can be efficiently simulated,” Phys. Rev. Lett. 101, 110501 (2008) ,“Entanglement renormalization,” Phys. Rev. Lett. 99, 220405 (2007) ,
  * (4) G. Evenbly and G. Vidal, “ Tensor Network Renormalization ,” arXiv:1412.0732 [cond-mat.str-el], Phys. Rev. Lett. 115, 180405 (2015); “ Tensor network renormalization yields the multi-scale entanglement renormalization ansatz, ” arXiv:1502.05385 [cond-mat.str-el], Phys. Rev. Lett. 115, 200401 (2015).
  * (5) J. Haegeman, T. J. Osborne, H. Verschelde and F. Verstraete, “Entanglement renormalization for quantum fields,” arXiv:1102.5524 [hep-th].
  * (6) X. L. Qi, “Exact holographic mapping and emergent space-time geometry,” [arXiv:1309.6282 [hep-th]].
  * (7) F. Pastawski, B. Yoshida, D. Harlow and J. Preskill, “Holographic quantum error-correcting codes: Toy models for the bulk/boundary correspondence,” JHEP 1506 (2015) 149 [arXiv:1503.06237 [hep-th]].
  * (8) P. Hayden, S. Nezami, X. L. Qi, N. Thomas, M. Walter and Z. Yang, “Holographic duality from random tensor networks,” arXiv:1601.01694 [hep-th].
  * (9) M. Van Raamsdonk, “Building up spacetime with quantum entanglement,” Gen. Rel. Grav. 42 (2010), 2323-2329 [arXiv:1005.3035 [hep-th]].
  * (10) S. Ryu and T. Takayanagi, “Holographic derivation of entanglement entropy from AdS/CFT,” Phys. Rev. Lett. 96 (2006) 181602 [hep-th/0603001].
  * (11) V. E. Hubeny, M. Rangamani and T. Takayanagi, “A Covariant holographic entanglement entropy proposal,” ]JHEP 0707 (2007) 062 [arXiv:0705.0016 [hep-th]].
  * (12) C. Beny, “Causal structure of the entanglement renormalization ansatz,” New J. Phys. 15 (2013) 023020 [arXiv:1110.4872 [quant-ph]].
  * (13) M. Nozaki, S. Ryu and T. Takayanagi, “Holographic Geometry of Entanglement Renormalization in Quantum Field Theories,” JHEP 1210 (2012) 193 [arXiv:1208.3469 [hep-th]].
  * (14) M. Miyaji and T. Takayanagi, “Surface/State Correspondence as a Generalized Holography,” PTEP 2015 (2015) 7, 073B03 [arXiv:1503.03542 [hep-th]].
  * (15) B. Czech, L. Lamprou, S. McCandlish and J. Sully, “Integral Geometry and Holography,” JHEP 1510 (2015) 175 [arXiv:1505.05515 [hep-th]].
  * (16) B. Czech, L. Lamprou, S. McCandlish and J. Sully, “Tensor Networks from Kinematic Space,” JHEP 07 (2016), 100 [arXiv:1512.01548 [hep-th]].
  * (17) A. Milsted and G. Vidal, “Tensor networks as path integral geometry,” [arXiv:1807.02501 [cond-mat.str-el]].
  * (18) A. Milsted and G. Vidal, “Geometric interpretation of the multi-scale entanglement renormalization ansatz,” [arXiv:1812.00529 [hep-th]].
  * (19) A. Jahn, Z. Zimborás and J. Eisert, “Tensor network models of AdS/qCFT,” [arXiv:2004.04173 [quant-ph]].
  * (20) T. Takayanagi, “Holographic Spacetimes as Quantum Circuits of Path-Integrations,” JHEP 12 (2018), 048 [arXiv:1808.09072 [hep-th]].
  * (21) M. Van Raamsdonk, “Building up spacetime with quantum entanglement II: It from BC-bit,” [arXiv:1809.01197 [hep-th]].
  * (22) N. Bao, G. Penington, J. Sorce and A. C. Wall, “Beyond Toy Models: Distilling Tensor Networks in Full AdS/CFT,” JHEP 19 (2020), 069 [arXiv:1812.01171 [hep-th]].
  * (23) N. Bao, G. Penington, J. Sorce and A. C. Wall, “Holographic Tensor Networks in Full AdS/CFT,” [arXiv:1902.10157 [hep-th]].
  * (24) P. Caputa, N. Kundu, M. Miyaji, T. Takayanagi and K. Watanabe, “Anti-de Sitter Space from Optimization of Path Integrals in Conformal Field Theories,” Phys. Rev. Lett. 119 (2017) no.7, 071602 [arXiv:1703.00456 [hep-th]]; “Liouville Action as Path-Integral Complexity: From Continuous Tensor Networks to AdS/CFT,” JHEP 11 (2017), 097 [arXiv:1706.07056 [hep-th]].
  * (25) M. Miyaji, T. Takayanagi and K. Watanabe, “From Path Integrals to Tensor Networks for AdS/CFT,” Phys. Rev. D 95 (2017) no.6, 066004 arXiv:1609.04645 [hep-th].
  * (26) A. M. Polyakov, “Quantum Geometry of Bosonic Strings,” Phys. Lett. 103B (1981) 207.
  * (27) L. Susskind, “Computational Complexity and Black Hole Horizons,” Fortsch. Phys. 64 (2016) 24 [arXiv:1403.5695 [hep-th], arXiv:1402.5674 [hep-th]]; A. R. Brown, D. A. Roberts, L. Susskind, B. Swingle and Y. Zhao, “Holographic Complexity Equals Bulk Action?,” Phys. Rev. Lett. 116 (2016) no.19, 191301 [arXiv:1509.07876 [hep-th]].
  * (28) B. Czech, “Einstein Equations from Varying Complexity,” Phys. Rev. Lett. 120 (2018) no.3, 031601 [arXiv:1706.00965 [hep-th]].
  * (29) P. Caputa and J. M. Magan, “Quantum Computation as Gravity,” Phys. Rev. Lett. 122 (2019) no.23, 231302 [arXiv:1807.04422 [hep-th]].
  * (30) H. A. Camargo, M. P. Heller, R. Jefferson and J. Knaute, “Path integral optimization as circuit complexity,” Phys. Rev. Lett. 123 (2019) no.1, 011601 [arXiv:1904.02713 [hep-th]].
  * (31) J. Erdmenger, M. Gerbershagen and A. L. Weigel, “Complexity measures from geometric actions on Virasoro and Kac-Moody orbits,” [arXiv:2004.03619 [hep-th]].
  * (32) Y. Sato and K. Watanabe, “Does Boundary Distinguish Complexities?,” JHEP 11 (2019), 132
  * (33) P. Caputa and I. MacCormack, “Geometry and Complexity of Path Integrals in Inhomogeneous CFTs,” [arXiv:2004.04698 [hep-th]].
  * (34) A. Bhattacharyya, P. Caputa, S. R. Das, N. Kundu, M. Miyaji and T. Takayanagi, “Path-Integral Complexity for Perturbed CFTs,” JHEP 07 (2018), 086 [arXiv:1804.01999 [hep-th]].
  * (35) P. Caputa, M. Miyaji, T. Takayanagi and K. Umemoto, “Holographic Entanglement of Purification from Conformal Field Theories,” Phys. Rev. Lett. 122 (2019) no.11, 111601 [arXiv:1812.05268 [hep-th]].
  * (36) H. A. Camargo, L. Hackl, M. P. Heller, A. Jahn, T. Takayanagi and B. Windt, “Entanglement and Complexity of Purification in (1+1)-dimensional free Conformal Field Theories,” [arXiv:2009.11881 [hep-th]].
  * (37) J. B. Hartle and S. W. Hawking, “Wave Function of the Universe,” Adv. Ser. Astrophys. Cosmol. 3 (1987), 174-189
  * (38) T. Takayanagi, “Holographic Dual of BCFT,” Phys. Rev. Lett. 107 (2011) 101602 [arXiv:1105.5165 [hep-th]]; M. Fujita, T. Takayanagi and E. Tonni, “Aspects of AdS/BCFT,” JHEP 11 (2011), 043 [arXiv:1108.5152 [hep-th]].
  * (39) G. Hayward, “Gravitational action for space-times with nonsmooth boundaries,” Phys. Rev. D 47 (1993), 3275-3280
  * (40) J. Boruch, P. Caputa, D. Ge and T. Takayanagi, “Holographic Path-Integral Optimization,” [arXiv:2104.00010 [hep-th]].
  * (41) J. M. Maldacena, “Eternal black holes in anti-de Sitter,” JHEP 0304 (2003) 021 [hep-th/0106112].
  * (42) A. Kitaev, “A simple model of quantum holography.” Talks at KITP, April 7, 2015 and May 27, 2015.
  * (43) S. Sachdev and J.-w. Ye, “Gapless spin fluid ground state in a random, quantum Heisenberg magnet,” Phys. Rev. Lett. 70 (1993) 3339,
  * (44) J. Maldacena and D. Stanford, “Remarks on the Sachdev-Ye-Kitaev model,” Phys. Rev. D 94 (2016) no.10, 106002 [arXiv:1604.07818 [hep-th]].
  * (45) J. Maldacena, D. Stanford and Z. Yang, “Conformal symmetry and its breaking in two dimensional Nearly Anti-de-Sitter space,” arXiv:1606.01857 [hep-th].
  * (46) K. Jensen, “Chaos in AdS2 Holography,” Phys. Rev. Lett. 117 (2016) no.11, 111601 [arXiv:1605.06098 [hep-th]].
  * (47) J. W. York, Jr., “Role of conformal three geometry in the dynamics of gravitation,” Phys. Rev. Lett. 28 (1972), 1082-1085
  * (48) A. Belin, A. Lewkowycz and G. Sárosi, “Complexity and the bulk volume, a new York time story,” JHEP 03 (2019), 044 [arXiv:1811.03097 [hep-th]].

##  Appendix A Appendix: Details of Evaluations of the On-Shell Action

(i) AdS${}_{d+1}/$CFTd ($d\geq 2$)

The free energy in classical gravity is computed as the on-shell action on the
region $M$ as in (11) and we proceed with this computation in detail. The
induced metric on $Q$ defined by $z=f(\tau)$ is given by

$\displaystyle ds^{2}={dx^{2}_{i}+(1+f^{\prime 2})d\tau^{2}\over f^{2}}\equiv
e^{2\phi(w)}\left(dw^{2}+dx^{2}_{i}\right),$ (35)

where we introduced coordinate $w$

${dw\over d\tau}=\sqrt{1+f^{\prime 2}}=\left(1-\dot{f}^{2}\right)^{-1/2},$
(36)

to make the metric diagonal and denoted
$f^{\prime}(\tau)=\partial_{\tau}f(\tau)$, $\partial_{w}f=\dot{f}$ as well as
$\phi=-\log(f)$. The (extrinsic) curvatures read $R-2\Lambda=-2d$,
$K|_{\Sigma}=d$ and

$\displaystyle K|_{Q}=-{ff^{\prime\prime}+d(1+f^{\prime 2})\over(1+f^{\prime
2})^{3/2}}={e^{-2\phi}\left(\ddot{\phi}+(d-1)\dot{\phi}^{2}\right)-d\over\sqrt{1-\dot{\phi}^{2}e^{-2\phi}}},$

where in the second equality we changed to $w$ coordinate and $\phi(w)$. In
the following, we also write the infinite lengths in the $x$ direction and
$\tau$ direction as $\int^{\infty}_{-\infty}d^{d-1}x=V_{x}$ and
$\int^{0}_{-\infty}d\tau=L_{\tau}$. We also denote $\kappa^{2}=8\pi G_{N}$.
The gravity action (Einstein-Hilbert with Gibbons-Hawking term) in region $M$
is evaluated as follows:

$\displaystyle I_{G}$ $\displaystyle=$
$\displaystyle{d\over\kappa^{2}}\int_{M}\sqrt{g}-{d\over\kappa^{2}}\int_{\Sigma}\sqrt{h}-{1\over\kappa^{2}}\int_{Q}\sqrt{h}K|_{Q}$
(37) $\displaystyle=$ $\displaystyle{d\,V_{x}\over\kappa^{2}}\int
d\tau\int^{f(\tau)}_{\epsilon}{dz\over
z^{d+1}}-{d\over\kappa^{2}}{V_{x}L_{\tau}\over\epsilon^{d}}$ $\displaystyle+$
$\displaystyle{V_{x}\over\kappa^{2}}\int{d\tau\over
f^{d}}\left({ff^{\prime\prime}\over 1+f^{\prime 2}}+d\right).$

After performing the $z$-integral, using coordinate $w$ and field $\phi$, we
can rewrite (37) as

$\displaystyle I_{G}=-{(d-1)\over\kappa^{2}}{V_{x}L_{\tau}\over\epsilon^{d}}+$
$\displaystyle{V_{x}\over\kappa^{2}}\int^{0}_{-\infty}dw{(d-1)e^{d\,\phi}-e^{(d-2)\phi}\left(\ddot{\phi}+(d-2)\dot{\phi}^{2}\right)\over\sqrt{1-\dot{\phi}^{2}e^{-2\phi}}}.$

Finally, we can perform partial integration to rewrite this action in the
first derivative form

$\displaystyle I_{G}$ $\displaystyle=$
$\displaystyle-{d-1\over\kappa^{2}}{V_{x}L_{\tau}\over\epsilon^{d}}+{(d-1)V_{x}\over\kappa^{2}}\times$
(38)
$\displaystyle\int^{0}_{-\infty}dwe^{d\phi}\left[\sqrt{1-\dot{\phi}^{2}e^{-2\phi}}+\dot{\phi}e^{-\phi}\arcsin\left(\dot{\phi}e^{-\phi}\right)\right]$
$\displaystyle-$
$\displaystyle{V_{x}\over\kappa^{2}}\left[e^{(d-1)\phi}\arcsin\left(\dot{\phi}e^{-\phi}\right)\right]^{0}_{-\infty}.$

The very last (co-dimension 2) boundary term can be written in terms of the
angle $\theta_{0}$ between $\Sigma$ and $Q$ as

$\displaystyle-{V_{x}\over\kappa^{2}}\left[e^{(d-1)\phi}\arcsin\left(\dot{\phi}e^{-\phi}\right)\right]^{0}_{-\infty}=-{V_{x}\over\kappa^{2}}{\theta_{0}\over\epsilon^{d-1}}.$
(39)

We can confirm that function $G(x)$ defined by
$G(\dot{\phi}e^{-\phi})=\sqrt{1-\dot{\phi}^{2}e^{-2\phi}}+\dot{\phi}e^{-\phi}\arcsin(\dot{\phi}e^{-\phi})$
is monotonically increasing. When

$\displaystyle\dot{\phi}e^{-\phi}\ll 1,$ (40)

we can approximate the finite term of the gravity action (38) as

$\displaystyle\int^{0}_{-\infty}dwe^{d\phi}\left[\sqrt{1-\dot{\phi}^{2}e^{-2\phi}}+\dot{\phi}e^{-\phi}\arcsin\left(\dot{\phi}e^{-\phi}\right)\right]$
$\displaystyle\simeq{1\over
2}\int^{0}_{-\infty}dw\left[e^{(d-2)}\dot{\phi}^{2}+2e^{d\phi}\right].$ (41)

Indeed, e.g. for $d=2$, using ${3\over 2G_{N}}=c$, leads to

$\displaystyle I_{G}\simeq{c\over 24\pi}\int
dxd\tau\left[\dot{\phi}^{2}+2(e^{2\phi}-\epsilon^{-2})\right]-{L_{x}\over 8\pi
G_{N}}\cdot{\theta_{0}\over\epsilon}+\cdot\cdot\cdot.$

agreeing with the kinetic terms of the Liouville action (5). Moreover, for
higher dimensions, we also get a perfect agreement with our boundary action
(see formula (8.2) in Caputa:2017urj ). Notice also that, if we add the
Hayward term Hayward:1993my , which is localized at the corner between $Q$ and
$\Sigma$:

$\displaystyle
I^{\prime}_{H}={1\over\kappa^{2}}\int_{Q\cap\Sigma}\sqrt{\gamma}\theta_{0}={L_{x}\theta_{0}\over
8\pi G_{N}\epsilon},$ (42)

we can eliminate the co-dimension 2 boundary term. See more in Boruch:2021hqs
.

(ii) JT gravity

Our solutions based on Einstein-Hilbert action appear to break down for $d=1$
($AdS_{2}$). For this case, we perform the analysis using JT gravity with
tension term coupled to the dilaton. Indeed we can show that starting from

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\Phi_{0}\left[\int_{M}\sqrt{g}R+2\int_{\partial
M}\sqrt{h}K\right]+2T_{\Phi}\int_{Q}\sqrt{h}\Phi$ (43) $\displaystyle-$
$\displaystyle\left[\int_{M}\sqrt{g}\Phi(R+2)+2\int_{\partial M}\sqrt{h}\Phi
K\right]$

we can analyze e.g. excited states from the metric and dilaton solutions

$ds^{2}=(r^{2}-r^{2}_{h})d\tau^{2}+{dr^{2}\over
r^{2}-r^{2}_{h}},\quad\Phi=\Phi_{b}r.$ (44)

Defining M bounded by $\Sigma$ at $r=r_{0}\to\infty$ and Q by $r=1/f(\tau)$,
with induced metric

$ds^{2}=e^{2\phi}dw^{2},\qquad e^{\phi}=f^{-1},$ (45)

we derive an action

$\displaystyle I$ $\displaystyle=$
$\displaystyle-2L_{\tau}\Phi_{b}r^{2}_{0}+2\Phi_{b}\int_{Q}e^{2\phi}G(\dot{\phi})$
$\displaystyle-$ $\displaystyle
2\left[\left(\Phi_{0}+\Phi_{b}e^{\phi}\right)\arcsin\left({\dot{\phi}e^{-\phi}\over\sqrt{1-r^{2}_{h}e^{-2\phi}}}\right)\right]^{0}_{-\infty},$

with

$\displaystyle G(\dot{\phi})$ $\displaystyle=$
$\displaystyle{\sqrt{1-e^{-2\phi}(\dot{\phi}^{2}+r^{2}_{h})}\over
1-r^{2}_{h}e^{-2\phi}}+T_{\Phi}$ (46) $\displaystyle+$
$\displaystyle\dot{\phi}e^{-\phi}\arcsin\left({\dot{\phi}e^{-\phi}\over\sqrt{1-r^{2}_{h}e^{-2\phi}}}\right).$

The equation of motion from this action is again equivalent to the Neumann
b.c. in $AdS_{2}$ (CMC slice)

$K_{Q}={\ddot{\phi}e^{-2\phi}-1\over\sqrt{1-e^{-2\phi}(\dot{\phi}^{2}+r^{2}_{h})}}=T_{\Phi},$
(47)

that can be solved for $T_{\Phi}<0$ by

$e^{2\phi}={r^{2}_{h}\over(1-T^{2}_{\Phi})\sin^{2}\left(r_{h}(w-c_{1})\right)}.$
(48)

Finally, we can check that expanding the finite contribution for small
$\dot{\phi}$ and small $r_{h}$, we get

$\Phi_{b}\int\left(\dot{\phi}^{2}+2(1+T_{\Phi})e^{2\phi}\right)+...$ (49)

With $e^{\phi}\sim\tilde{w}^{\prime}(w)$, this describes the Schwarzian theory
see e.g. SYKd ; Jensen:2016pah .