\section{Introduction}

Machine learning is increasingly applied in control approaches, where first principle
models are not available or expensive to obtain due to the complexity of systems. In
order to apply learning based control approaches to real world applications, it is
crucial to certify their safety using theoretically rigorous methods. Although there exists
a wide variety of machine learning methods, this theoretical analysis is often difficult,
such that safe control approaches typically rely on a few, well-understood
learning methods.\looseness=-1

Gaussian process (GP) regression is such a machine learning method which bases on
solid Bayesian foundations. Due to its inherent bias-variance trade-off it allows efficient
generalization from few training data, which makes it an appealing method both for control
practitioners and theoreticians, and has lead to its increasing use in control. For the control
theoretic analysis of Gaussian processes, several approaches have been developed. In
\citep{Beckers2016a, Beckers2016} stochastic stability and equilibria of Gaussian processes with
certain covariance kernels are investigated analytically. Employing numerical quadrature, an
efficient method to analyze stability of the posterior mean function of a Gaussian process state
space model with squared exponential kernel is developed in \citep{Vinogradska2017a}. A statistical
learning error analysis of learned MPC control laws is proposed in \citep{Hertneck2018} to ensure
closed-loop
stability. By extending the Lyapunov stability conditions on a grid to a continuous state space
using Lipschitz continuity, the region of attraction of a Gaussian process state space model
is learned in \citep{Berkenkamp2016a}. A similar approach is employed to determine the region of
attraction of a nonlinear system based on the learned infinite horizon cost function in
\citep{Lederer2019b}. The Lipschitz constants of the posterior mean function used in these approaches
are also of interest themselves, e.g., for computing uniform regression error bounds \citep{Lederer2019}
or in batch parallelization of Bayesian optimization in \citep{Gonzalez2016}.

Despite of the similarity of many of these problems, each one is based on a separate analysis method,
which is typically not optimized for computational efficiency. In order to overcome this issue,
we propose a novel framework for the analysis of Gaussian process mean functions called GP3:
\textbf{G}eneral \textbf{P}urpose computation on \textbf{G}raphics \textbf{P}rocessing units for
\textbf{G}aussian \textbf{P}rocesses. By defining a common problem formulation for many problems,
we can employ interval analysis to derive local Lipschitz constants on hyperrectangles
covering the region of interest. These local Lipschitz constants allow to extend verified
properties on a discrete grid to the whole region of interest such that multi-resolution
sampling can be used for efficient analysis. As the Lipschitz constants can be computed
independently for each hyperrectangle, the method is parallelized using general purpose
graphics processing units in order to exploit the full computational power of modern hardware.
We demonstrate the flexibility and efficiency of the GP3 framework by applying it to a region
of attraction estimation and a Lipschitz constant bounding problem.

The remaining paper is structured as follows. In Section~\ref{sec:prob} we define the general
property analysis problem. The theoretical background on Gaussian process regression and
interval analysis is presented in Section~\ref{sec:theoBack}. The theoretical foundations of the
GP3 framework are explained in Section~\ref{sec:GP_Ana}, before it is evaluated in simulations
in Section~\ref{sec:numeval}.

\section{Problem Statement}
\label{sec:prob}

Consider the posterior mean function\footnote{Lower/upper case bold symbols denote vectors/matrices,
	$\mathbb{R}_+$ denotes all real positive numbers,~$\bm{I}_n$ the~$n\times n$ identity matrix and
	$\|\cdot\|$ the Euclidean norm
	}
$\mu:\mathbb{R}^d\rightarrow\mathbb{R}$ of a Gaussian process with continuous covariance
kernel \mbox{$k:\mathbb{R}^d\times\mathbb{R}^d\rightarrow\mathbb{R}$}, a continuous function
$g:\mathbb{R}^d\rightarrow\mathbb{R}$ for comparison and a continuous state transformation
$\bm{f}:\mathbb{R}^d\rightarrow\mathbb{R}^d$. An abstract problem that finds many practical
applications is to find bounds~$\epsilon_1,\epsilon_2\in\mathbb{R}_+$ on the difference
\begin{align}
\label{eq:prob}
	-\epsilon_1\leq g(\bm{f}(\bm{x}))-\mu(\bm{x}) \leq \epsilon_2 \quad\forall\bm{x}\in\mathbb{X}
\end{align}
on a compact set~$\mathbb{X}\!\subset\!\mathbb{R}^d$. Various choices for the functions~$\bm{f}(\cdot), g(\cdot)$ immediately come to mind: if we choose~$g(\cdot)$ as the function generating the training data of the Gaussian process and set~$\bm{f}(\bm{x})\!=\!\bm{x}$ we can determine the maximum learning error. Furthermore,~$\bm{g}(\cdot)$ can be the Gaussian process mean~$\mu(\cdot)$ itself and~$\bm{f}(\cdot)$ can be defined such that it returns the closest point in a discrete set. Thereby the variation of the mean function with respect to the discrete set is analyzed. Finally, we can consider~$\bm{f}(\cdot)$ as autonomous, discrete-time dynamics and choose~$g(\cdot)\!=\!\mu(\cdot)$, such that we can immediately investigate if the mean function~$\mu(\cdot)$ satisfies the second condition of Lyapunov's theorem. Although this problem formulation offers such a high flexibility, it allows a straight forward, uniform treatment through interval analysis, which we exploit in the following sections.\looseness=-1

\section{Theoretical Background}
\label{sec:theoBack}

\subsection{Gaussian Process Regression}

Gaussian process regression is a supervised machine learning method, which is frequently applied
in control and system identification due to its Bayesian foundations. The Gaussian process
distribution assigns to any finite subset~$\{ \bm{x}_i,\ldots,\bm{x}_N \}\subset\mathbb{X}$ from
a continuous input domain~$\mathbb{X}\subset\mathbb{R}^d$ a joint Gaussian distribution~\citep{Rasmussen2006}.
It is completely defined through the mean function~$m:\mathbb{R}^d\rightarrow\mathbb{R}$ and
the covariance function~$k:\mathbb{R}^d\times\mathbb{R}^d\rightarrow\mathbb{R}$. Although prior
information like approximate models can be incorporated as mean function~$m(\cdot)$ into GP
regression, such knowledge is often not available, such that the mean function~$m(\cdot)$ is
usually set to zero. We also assume this in the following. In contrast, a wide variety of
different covariance functions~$k(\cdot,\cdot)$ is applied in GP regression to encode
prior information such as smoothness, periodicity and stationarity. Frequently used covariance
functions are the squared exponential kernel
\begin{align}
k_{\mathrm{SE}}(r)=\sigma_f^2\exp\left( -\frac{1}{2}r^2 \right)
\end{align}
and the Mat\'ern class kernels
\begin{align}
&k_{m+\frac{1}{2}}(r)=\sigma_f^2p_m(r)\exp\left( -\sqrt{(2m+1)}r \right)
\end{align}
with signal variance~$\sigma_f^2\in\mathbb{R}_+$,~$p_m(\cdot)$ a polynomial of order
$m\in\mathbb{N}$ and automatic relevance determination distance~\citep{Neal1996}
\begin{align}
\label{eq:autorel}
	r=\sqrt{\sum\limits_{i=1}^d\frac{(x_i-x_i')^2}{l_i^2}}
\end{align}
with length scales~$l_i\in\mathbb{R}_+$. Mat\'ern class kernels are commonly used with
$m=1$ or~$m=2$ which results in the polynomials
\begin{align}
	p_1(r)&=1+\sqrt{3}r\\
	p_2(r)&=\left(1+\sqrt{5}r+\frac{5}{3}r^2\right).
\end{align}
A major reason for their popularity is the fact that a GP with either of these covariance
functions is a universal approximator~\citep{Steinwart2001}, i.e., any continuous function
can be approximated with arbitrary precision.

GP regression is based on the assumption that the training data set
$\mathbb{D}=\{ (\bm{x}^{(n)},y^{(n)}) \}_{n=1}^N$ is generated through
noisy observations of a function~$h:\mathbb{R}^d\rightarrow\mathbb{R}$, i.e.,
\begin{align}
y^{(i)}= h(\bm{x}^{(i)})+\omega^{(i)},
\end{align}
where~$\omega^{(i)}\sim\mathcal{N}(0,\sigma_n^2)$ are i.i.d. random variables.
By conditioning the prior joint Gaussian distribution of a prediction~$h(\bm{x})$
and the training outputs~$y^{(i)}$ on the training data~$\mathbb{D}$, we obtain
the predictive mean
\begin{align}
\label{eq:mean}
\mu(\bm{x})=\mathbb{E}[h(\bm{x})|\mathbb{D},\bm{x}]=\bm{k}^T(\bm{x})\bm{\lambda},
\end{align}
where
\begin{align}
\bm{\lambda}=(\bm{K}+\sigma_n^2\bm{I}_N)^{-1}\bm{y}
\end{align}
and the data covariance matrix~$\bm{K}\in\mathbb{R}^{N\times N}$ and the
covariance vector~$\bm{k}(\bm{x})\in\mathbb{R}^{N}$
are defined through~$K_{ij}=k(\bm{x}^{(i)},\bm{x}^{(j)})$ and
$k_i(\bm{x})=k(\bm{x}^{(i)},\bm{x})$. The observation noise variance~$\sigma_n^2$, the
signal variance~$\sigma_f^2$ and the length scales~$l_i$ are considered hyperparameters
of the GP regression and can be determined by maximization of the
log-likelihood~\citep{Rasmussen2006}.

\subsection{Interval Analysis for Property Analysis}

Interval analysis is a method to approach the problem of calculating bounds of functions. Instead of operating on exact values, interval analysis uses real compact intervals~$[a] = [\underline{a},\bar{a}] = \{x \in \mathbb{R}|\underline{a} \leq x \leq \bar{a}\}$. Basic mathematical operations for intervals are defined as
\begin{align}
	[a] + [b] &= [\underline{a} + \underline{b}, \overline{a} + \overline{b}], \label{eq:IntervalAdd}\\
	[a] - [b] &= [\underline{a} - \overline{b}, \overline{a} - \underline{b}], \label{eq:IntervalSub}\\
	\label{eq:IntervalMul}
	[a] \cdot [b] &= [\mathrm{min}\{\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b}\}, \mathrm{max}\{\underline{a}\underline{b},\underline{a}\overline{b},\overline{a}\underline{b},\overline{a}\overline{b}\}].
\end{align}
A more thorough introduction into operations on intervals can be found, e.g., in~\citep{Alefeld2000}.

Based on this interval arithmetic, it is possible to use intervals as inputs to functions~$f: \mathbb{R} \rightarrow \mathbb{R}$ and calculate output intervals~$[y] = f([x])$, where~$\underline{y}$ serves as a lower bound and~$\bar{y}$ as an upper bound of the function values on the input interval. It is straight forward to adapt this approach to higher dimensional functions~$f: \mathbb{R}^d \rightarrow \mathbb{R}$ by considering so called \textit{hyperrectangles} instead of intervals. A hyperrectangle is completely defined by its center~$\bm{c}$ and length parameter~$\bm{b}$, such that it defines the multi-dimensional interval~$[\bm{c}-\bm{b},\bm{c}+\bm{b}]$ with edge widths~$2\bm{b}$. Using a grid of hyperrectangles, interval analysis allows to efficiently expand the validity of \eqref{eq:prob} on the hyperrectangle centers to the area covered by the hyperrectangles. Using multi-resolution grids as proposed, e.g., in \citep{Bobiti2018}, this enables efficient determination of valid constants~$\epsilon_1$ and~$\epsilon_2$ in \eqref{eq:prob}.

\section{Sampling-Based Analysis of Gaussian Processes}
\label{sec:GP_Ana}

Due to the strong nonlinearity of typical Gaussian process mean functions,
standard interval operations are not directly applicable to their analysis.
Therefore, we develop an efficient multi-resolution sampling algorithm for
the analysis of Gaussian processes in Section~\ref{subsec:multi-res}. As this
algorithm requires upper and lower bounds for the derivative of covariance
kernels depending on the training data, we investigate such bounds for squared exponential
and Mat\'ern class kernels in Section~\ref{subsec:deriv int}.

\subsection{Multi-resolution Analysis of Gaussian Processes}
\label{subsec:multi-res}

Gaussian processes exhibit a strongly nonlinear mean function in general. In order to efficiently analyze
their mean functions we develop a multi-resolution sampling algorithm in this section.
Exploiting Lipschitz continuity we derive a theorem to calculate the bounds~$\epsilon_1$ and~$\epsilon_2$ in \eqref{eq:prob} inside hyperrectangles. In order to evaluate the bound, we determine a local Lipschitz constant of the GP on a hyperrectangle. Due to the sampling structure, the analysis of hyperrectangles can be efficiently parallelized using GPUs in each grid refinement iteration.

The basis of our approach lies in the independent analysis of the bounds~$\epsilon_1$ and~$\epsilon_2$ on hyperrectangles. This analysis is founded on the following theorem, which relies on Lipschitz continuity of all involved functions.
\begin{thm}
\label{thmbounds}
	Consider a function~$\bm{f}(\cdot)$ with local Lipschitz constant~$L_f$ and a posterior mean function
	$\mu(\cdot)$ of a Gaussian process with local Lipschitz constant~$L_{\mu}$ on a hyperrectangle
	with center~$\bm{c}$ and edge widths~$2\bm{b}$. Furthermore, assume the function~$g(\cdot)$ has
	a local Lipschitz constant~$L_g$ on the hyperrectangle with center~$\bm{f}(\bm{c})$ and edge widths
	$2L_f\bm{b}$. Then, \eqref{eq:prob} holds with
	\begin{align}
	\label{eq:eps1}
		\epsilon_1&=g(\bm{f}(\bm{c}))-\mu(\bm{c})-(L_fL_g+L_{\mu})\|\bm{b}\|\\
		\epsilon_2&=g(\bm{f}(\bm{c}))-\mu(\bm{c})+(L_fL_g+L_{\mu})\|\bm{b}\|
		\label{eq:eps2}
	\end{align}
	for all~$|\bm{x}-\bm{c}|\leq \bm{b}$, where the
	absolute values and the comparison are performed element-wise.
\end{thm}
\begin{pf}
	Due to Lipschitz continuity of~$\bm{f}(\cdot)$ and~$g(\cdot)$ it follows that
	\begin{align*}
		|g(\bm{f}(\bm{x}))-g(\bm{f}(\bm{c}))|\leq L_fL_g\|\bm{b}\| \quad\forall\bm{x}:\|\bm{x}-\bm{c}\|\leq \|\bm{b}\|.
	\end{align*}
	Furthermore, Lipschitz continuity of the GP yields
	\begin{align*}
		|\mu(\bm{x})-\mu(\bm{c})|\leq L_{\mu}\|\bm{b}\|.
	\end{align*}
	Applying the triangle inequality we finally obtain the bounds \eqref{eq:eps1} and \eqref{eq:eps2}.
\end{pf}

The application of this theorem crucially relies on the Lipschitz constant of the GP mean function. The computation of this Lipschitz constant on a hyperrectangle can be performed independently for each hyperrectangle and depends merely on the parameter vector~$\bm{\lambda}$, the kernel function~$k(\cdot,\cdot)$ and the data points~$\bm{x}^{(i)}$,~$i=1,\ldots,N$, as shown in the following theorem.

\begin{thm}
\label{thmlipconst}
	Consider the posterior mean function~$\mu(\cdot)=\bm{k}^T(\cdot)\bm{\lambda}$
	of a GP with covariance kernel~$k(\cdot,\cdot)$ on a hyperrectangle with
	center~$\bm{x}$ and length~$2\bm{b}$. Let~$\bm{L}_{k,i}^{\partial j}$ denote
	a vector of partial derivative bounds of the~$i$-th element of~$\bm{k}(\cdot)$
	with respect to~$x_j$. Then, a local Lipschitz constant of the mean on the
	hyperrectangle is given by
	\begin{align}
		L_{\mu}(\bm{x},\bm{b})=\sqrt{\sum\limits_{j=1}^d \max\left\{ \left(\sum\limits_{i=1}^N\bm{R}(\lambda_i)\lambda_i\bm{L}_{k,i}^{\partial j} \right)^2 \right\}},
	\end{align}
	where
	\begin{align}
		\bm{R}(\lambda_i)=\begin{cases}
		\bm{I}_2  & \lambda_i>0\\
		\begin{bmatrix}
		0&1\\1&0
		\end{bmatrix}&\lambda_i\leq 0.
		\end{cases}
	\end{align}
\end{thm}
\begin{pf}
	Due to the scalar product of parameter vector~$\bm{\lambda}$ and the kernel vector~$\bm{k}(\bm{x})$ in \eqref{eq:mean}, the partial derivatives of the mean function are given by
	\begin{align*}
		\frac{\partial}{\partial x_j}\mu(\bm{x})=\begin{bmatrix}
		\frac{\partial}{\partial x_j} k_1(\bm{x}) & \ldots & \frac{\partial}{\partial x_j} k_N(\bm{x})
		\end{bmatrix}\bm{\lambda}.
	\end{align*}
	The vectors~$\bm{L}_{k,i}^{\partial j}$ contain upper and lower bound on the partial
	derivatives~$\frac{\partial}{\partial x_j} k_i(\bm{x})$ in the first and second row, respectively.
	By multiplying these vectors with the matrix~$\bm{R}(\lambda_i)$, the order of elements is
	changed if~$\lambda_i$ is negative, such that upper derivative bounds are multiplied with
	positive~$\lambda_i$s and lower bounds with negative ones in the first row of the resulting vector.
	The multiplication is performed in the inverse combination in the second row. By summing up each
	row, upper and lower bounds on the partial derivatives are obtained. We take the maximum squared
	value of these two rows as squared Lipschitz constant in the~$j$-th direction and finally,
	calculate the overall Lipschitz~$L_{\mu}(\bm{x},\bm{b})$ constant by taking the square root of
	the sum of squared Lipschitz constants in all directions.\looseness=-1
\end{pf}

Theorem \ref{thmbounds} and Theorem \ref{thmlipconst} allow for a straightforward implementation and integration in a multi-resolution, GPU parallelized sampling algorithm, which is depicted in Algorithm \ref{Alg1}. Initially, a sampling grid is created by dividing the analyzed region~$\mathbb{X}$ into~$M$ hyperrectangles in line 2, where the centers~$\bm{c}_i$ are concatenated in~$\bm{C}$ and edge-lengths~$\bm{b}_i$ are concatenated in~$\bm{B}$. Generally, any initial sampled grid is possible for this step, since our approach does not crucially depend on it. The grid is refined in line 6 by using an indication value~$s_i$, which divides hyperrectangles that must be further refined for~$s_i = \texttt{false}$ from those that require no further refinement. For the refinement procedure itself, several ways are possible in general. A simple method is to resample the grid with smaller hyperrectangle sizes~$\bm{B}$ and skip the calculation where~$s_i=\texttt{true}$ already. Another example for a refinement method is the~$2$-refinement used in~\citep{Bobiti2018}, where hyperrectangles with~$s_i=\texttt{false}$ are divided into~$2^d$ smaller hyperrectangles with new edge widths~$\frac{1}{2}\bm{b}_i$. \looseness=-1

The computation of the bounds~$\epsilon_1,\epsilon_2$ based on Theorems~\ref{thmbounds} and~\ref{thmlipconst} takes place inside the parallelized For (\textbf{ParFor}) loop of line 7. By shifting the execution to a highly parallel GPU, the computation time of the algorithm can be significantly reduced, which is an important feature of this approach. The termination of the refinement process of the algorithm is ensured by setting~$s_i \!=\! \texttt{true}$, once the bounds~$\epsilon_1,\epsilon_2$ inside the hyperrectangle satisfy predefined desired bounding functions~$\bar{\epsilon}_1(\bm{c}_i)$ and~$\bar{\epsilon}_2(\bm{c}_i)$ or the size~$\|\bm{b}_i\|$ of the hyperrectangle falls below a specified minimum size~$b_{\mathrm{min}}$.\looseness=-1
\begin{algorithm}[t!]
\caption{Multi-resolution bound calculation}
\label{Alg1}
\begin{algorithmic}[1]
\Function{Bounds}{$g(\cdot),\mu(\cdot),\mathbb{X},M,\bar{\epsilon}_1(\cdot),\bar{\epsilon}_2(\cdot),b_{\mathrm{min}}$}
  \State~$[\bm{C};\bm{B}] \gets$ \Call{SamplingGrid}{$\mathbb{X},M$}
  \State~$\bm{\epsilon}_1,\bm{\epsilon}_2 \gets \texttt{zeros}[\textsf{len}(\bm{C})]$
  \State~$\bm{s}\gets\texttt{false}[\textsf{len}(\bm{C})]$
  \While{$\textsf{any}(s_i) = \texttt{false}$}
  	\State~$[\bm{C};\bm{B}] \gets$ \Call{RefineGrid}{$\bm{C};\bm{B};s_i$}
  	\ParFor{$i=1$ to len$(\bm{C})$}
  		\State Calculate~$L_{\mu}(\bm{c}_i, \bm{b}_i)$ using Theorem~\ref{thmlipconst}
  		\State Calculate~$\epsilon_{1,i},\epsilon_{2,i}$ based on Theorem~\ref{thmbounds}
  		\State$s_i \gets (\epsilon_{1,i}\leq\bar{\epsilon}_1(\bm{c}_i) \land \epsilon_{2,i}\leq\bar{\epsilon}_2(\bm{c}_i)) \lor$
  		\State\hphantom{$s_i \gets\ $}$(\|\bm{b}_i\| \leq b_{\mathrm{min}})$
  		\EndParFor
  \EndWhile
\State \Return~$\bm{C},\bm{\epsilon}_1,\bm{\epsilon}_2$
\EndFunction
\end{algorithmic}
\end{algorithm}

\subsection{Derivative Intervals of Kernels}
\label{subsec:deriv int}

Typical covariance kernels are highly nonlinear functions themselves, which
complicates the interval analysis of their derivatives. However, many kernels
exhibit a structure such that derivatives are monotonous on large intervals.
We exploit this behavior for deriving upper and lower bounds for squared
exponential and Mat\'ern class kernels in the following theorem.
\begin{thm}
	Consider a multivariate squared exponential or Mat\'ern class
	kernel~$k(\cdot)$ with length scale~$\bm{l}$
	on a hyperrectangle with center~$\bm{x}$ and edge lengths~$2\bm{b}$.
	Then, the derivative bounds for a training point~$\bm{x}^{(i)}$ with
	respect to the~$j$-th component is given by
	\begin{align}
	\label{eq:boundsmult}
	\bm{L}_{k,i}^{\partial j}=\bm{T}(x_j^{(i)}-x_j)\bm{\kappa}(\bm{x}^{(i)}-\bm{x})
	\end{align}
	where
	\begin{align}
	&\bm{\kappa}(\bm{\Delta})=\\
	&\begin{cases}
	\begin{bmatrix}
	-\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta-\overline{\bm{\tau}},|\Delta_j|-b_j\right)\right)\\
	-\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta-\underline{\bm{\tau}},|\Delta_j|+b_j\right)\right)
	\end{bmatrix}& |\Delta_j|>b_j+\tilde{l}_j\\[18pt]
	\begin{bmatrix}
	-\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta-\overline{\bm{\tau}},|\Delta_j|+b_j\right)\right)\\
	-\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta-\underline{\bm{\tau}},|\Delta_j|-b_j\right)\right)
	\end{bmatrix}& |\Delta_j|<\tilde{l}_j-b_j\\[18pt]
	\begin{bmatrix}
	-\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta-\overline{\bm{\tau}},\tilde{l}_j\right)\right)\\
	\!\min\left\{\! -\frac{\partial }{\partial \Delta_j}k\left(\rho_j\left(\Delta\!-\!\underline{\bm{\tau}},|\Delta_j|\!\pm\! b_j\right)\right) \!\right\}\!
	\end{bmatrix}\!&\!||\Delta_j|\!-\!\tilde{l}_j|\leq b_j
	\end{cases}
	\label{eq:kappa}\\
	&\bm{T}(\Delta_j)=\begin{cases}
	\bm{I}_2&\Delta_j>0\\
	\begin{bmatrix}
	0&-1\\-1&0
	\end{bmatrix}&\Delta_j\leq 0
	\end{cases},
	\end{align}
	with modified weighted distance~$\rho_j(\bm{\Delta},\tilde{\Delta})$, maximum point~$\tilde{l}_j$,
	maximum distance~$\overline{\bm{\tau}}$, minimum distance~$\underline{\bm{\tau}}$ and index
	set~$\mathbb{I}_j$ defined as follows
	\allowdisplaybreaks
	\begin{align}
	\rho_j(\!\bm{\Delta},\tilde{\Delta}\!)&=\sqrt{\!\frac{\tilde{\Delta}^2}{l_j^2}\!+\!\sum\limits_{i\in\mathbb{I}_j}\!
		\frac{\Delta_i^2}{l_i^2}\!}\\
	\tilde{l}_j&=\begin{cases}
	l_j\!& \!\text{squared exponential kernel}\\
	\frac{l_j}{\sqrt{3}}\!&\!\text{Mat\'ern kernel with }m=1\\
	\frac{5\!+\!\sqrt{5}}{10}l_j\!&\!\text{Mat\'ern kernel with }m=2
	\end{cases}\\
	\label{eq:taubar}
	\overline{\bm{\tau}}&= \min\left\{ \|\bm{b}_{\neq j}\|, \|\bm{x}_{\neq j}^{(i)}\!-\!\bm{x}_{\neq j} \|\right\}
	\frac{\bm{x}^{(i)}\!-\!\bm{x}}{\|\bm{x}_{\neq j}^{(i)}\!-\!\bm{x}_{\neq j} \|}\\
	\label{eq:tauunder}
	\underline{\bm{\tau}}&=\|\bm{b}_{\neq j}\|\frac{\bm{x}\!-\!\bm{x}^{(i)}}{\|\bm{x}_{\neq j}^{(i)}\!-\!\bm{x}_{\neq j} \|}\\
	\mathbb{I}_j&=\{1,\ldots,j-1,j+1,\ldots,d\}.
	\end{align}
\end{thm}
\begin{pf}
	We start this proof by first showing that \eqref{eq:boundsmult} holds for
	the squared exponential before we highlight the differences in the proof
	for Mat\'ern kernels. Due to the exponent product rule and the automatic relevance
	determination \eqref{eq:autorel}, the derivative of
	the multivariate squared exponential kernel can be split into
	\begin{align*}
	\frac{\partial }{\partial x_j}k(\bm{x}^{(i)},\bm{x})=-\frac{\partial }{\partial \Delta_j}k\left(\frac{|\Delta_j|}{l_j}\right)k\left(\sqrt{\sum\limits_{i\in\mathbb{I}_j} \frac{\Delta_j^2}{l_i^2}}\right),
	\end{align*}
	where~$\bm{\Delta}=\bm{x}^{(i)}-\bm{x}$. This separation allows to determine interval
	bounds for both factors independently and combine the maximizers or minimizers, respectively,
	with the help of~$\rho_j(\bm{\Delta},\Delta_j)$ based on \eqref{eq:IntervalMul}. Due to odd
	symmetry of the univariate derivative, it is sufficient to derive only bounds for the positive
	real line and obtain the interval bounds for the negative real line by multiplication with
	$\bm{T}(\Delta_j)$. The kernel derivative is monotonous on the intervals~$[0,l_j]$ and~$[l_j,\infty]$
	with a maximum at~$l_j$. Therefore, we obtain three intervals with maximizers~$\nu$
	\begin{align*}
		|\Delta_j|>b_j+l_j&: \nu=|\Delta_j|-b_j\\
		|\Delta_j|<l_j-b_j&: \nu=|\Delta_j|+b_j\\
		||\Delta_j|-l_j|\leq b_j&: \nu=l_j
	\end{align*}
	which can be analogously obtained for the minimizers.
	The multivariate kernel is trivially maximized by considering the minimal distance to the
	training point, which is given by~$\bm{\Delta}-\overline{\bm{\tau}}$, while
	it is minimized by~$\bm{\Delta}-\underline{\bm{\tau}}$. Therefore, \eqref{eq:boundsmult}
	provides an upper and lower bound for the derivative of the squared exponential kernel
	on the hyperrectangle and we proceed with Mat\'ern class kernels. Although the derivatives
	of Mat\'ern class kernels cannot be separated as the squared exponential kernel, they
	exhibit a similar behavior in that they are monotonous in the non-derived directions
	$i\neq j$. Therefore, we can choose the same maximizer and minimizer in these directions
	as for the squared exponential kernel. Furthermore, the derivatives of Mat\'ern class
	kernels exhibit the same odd symmetry in the derived direction~$j$ and have two monotonous
	intervals on the positive real line with maximum at~$\tilde{l}_j=l_j/\sqrt{3}$ for
	$m=1$ and~$\tilde{l}_j=(5+\sqrt{5})l_j/10$. Therefore, we can upper and
	lower bound the derivative for Mat\'ern class kernels on the positive real line using
	\eqref{eq:kappa}, which concludes the proof.\looseness=-1
\end{pf}
Although this theorem might appear rather complicated, it allows a straightforward
implementation with few conditional operators, which is beneficial for GPU parallelization~\citep{GPUcomputing}.
Therefore, it allows efficient, parallelized analysis of Gaussian process mean functions
in combination with Theorems~\ref{thmbounds} and~\ref{thmlipconst}.

\section{Numerical Evaluation}
\label{sec:numeval}

\subsection{Efficiency of the GP3 Framework}
\label{subsec:effic}
In order to demonstrate the advantages of the GP3 framework\footnote{Code is available at \url{https://gitlab.lrz.de/alederer/gp3}}
over a global Lipschitz
constant computation and a CPU parallelized implementation, we compare these approaches
on the example proposed in~\citep{Lederer2019}. A crucial step therein is the derivation of
a Lipschitz constant for a posterior mean representing the function\looseness=-1
\begin{align}
	f(\bm{x})=1-\sin(x_1)+\frac{1}{1+\exp(-x_2)}.
\end{align}
The Gaussian process is trained with~$100$ samples which are uniformly spaced
over the analyzed region
$\mathbb{X}=[-6,4]\times[-4,4]$.
We train three GPs with a squared exponential and
Mat\'ern kernels with~$m=1$ and~$m=2$. The hyperparameters
obtained via log-likelihood maximization for each of these GPs are depicted in
Table~\ref{tab:hyp_param}.
\begin{table}[b]
	\centering
	\begin{tabular}{l| c c c}
		kernel &~$\sigma_n^2$ &~$\sigma_f^2$ &~$\bm{l}$\\ \hline
		squared exponential &~$0.1$ &~$0.956$ &~$[1.762\ 5.537]^T$\\
		Mat\'ern with~$m=1$ &~$0.1$ &~$1.274$ &~$[3.755\ 15.052]^T$\\
		Mat\'ern with~$m=2$ &~$0.1$ &~$1.012$ &~$[2.333\ 8.496]^T$
		\rule{0pt}{3ex}
	\end{tabular}
	\caption{Hyperparameters of the different covariance kernels}
	\label{tab:hyp_param}
\end{table}

We compute the Lipschitz constants of the posterior mean functions using
the GP3 framework for different numbers of hyperrectangles and the naive, global
Lipschitz constant proposed in~\citep{Lederer2019}. The Lipschitz constants obtained from
the GP3 approach are depicted in Figure~\ref{fig:Lip}. Since the naive approach yields
Lipschitz constants of~$343.62$,~$245.66$ and~$268.20$ for the squared exponential and
Mat\'ern kernel with~$m=1$ and~$m=2$, respectively,
the corresponding constant curves are not displayed in the figure. In contrast, the
Lipschitz constants obtained by the GP3 approach can provide Lipschitz constants
smaller than~$20$ with less than~$2000$ hyperrectangles. Furthermore, with
approximately~$10^6$ hyperrectangles the Lipschitz constant has almost converged to
a constant value.\looseness=-1
\begin{figure}
	\centering
	\includegraphics[]{plots/lip_val}
	\vspace{-0.3cm}
	\caption{Decrease rate of Lipschitz constants for increasing number of analyzed hyperrectangles}
	\label{fig:Lip}
\end{figure}

Additionally, we compare the time required to compute the Lipschitz constants on the
hyperrectangles with a GPU and CPU parallelization. We compute the Lipschitz constants
on a system with a NVIDIA TITAN V GPU, which has~$2\!\!:\!\!1$ single to double precision unit
ratio, two AMD EPYC~$16$-Core CPUs and~$1$TB RAM. The speedup of the GPU
compared to the CPU parallelization averaged over~$10$ runs of both
implementations are displayed for the three different covariance kernels in
Figure~\ref{fig:comptime}. Although we parallelize the Lipschitz constant computation
with~$64$ threads on the CPU, the GPU parallelization achieves a speedup
of at least~$30$ for large numbers of hyperrectangles. Merely at low numbers of hyperrectangles
the speedup can be in the single digit region due to the computational overhead of GPU computation.
Therefore, GP3 allows to exploit the advantages of modern general purpose GPU computing to achieve
low computation times, while requiring a reasonable amount of hyperrectangles for convergence of the
obtained Lipschitz constant.

\begin{figure}
	\centering
	\includegraphics[]{plots/time_val}
	\vspace{-0.3cm}
	\caption{Average speedup of GPU over CPU parallelization over the number of analyzed hyperrectangles}
	\label{fig:comptime}
\end{figure}

\subsection{Region of Attraction for Power Systems}

As an example for a problem of the form \eqref{eq:prob}, we analyze the region of attraction
of a nonlinear autonomous system. We consider the single machine infinite bus
system~\citep{Munz2013}
\begin{align}
\label{eq:smibs}
m_1\ddot{\phi}_1+d_1\dot{\phi}_1=-a_{12}(\sin(\theta_1+\phi_1)-\sin(\theta_1)),
\end{align}
which models a synchronous machine with inertia~$m_1$, damping~$d_1=20$ and steady state
phase~$\theta_1=\arcsin(0.05)$ as generator bus connected to an infinite bus with~$\theta_2=\phi_2=0$,
such that~$a_{12}=10$ is the product of the susceptance between both buses and the root mean square
voltages~$u_1$ and~$u_2$ at bus~$1$ and~$2$, respectively. For analyzing the regions, where the system
is discrete-time stable, we consider the cost of finite time trajectories proposed in \citep{Bobiti2018}
as Lyapunov function, i.e.,
\begin{align}
\label{eq:Lyap}
	V(\bm{x}_0)=\sum\limits_{k=0}^K\bm{x}^2(k\Delta t,\bm{x}_0)
\end{align}
where~$\bm{x}(t,\bm{x}_0)\!=\![\dot{\phi}(t)\ \phi(t)]$ denotes the solution of the
differential equation \eqref{eq:smibs} for initial state~$\bm{x}_0$ and
\linebreak~$\Delta t\!\in\!\mathbb{R}_+$ is the sampling time. For determining the states at
the sampling times~$k\Delta t$, we numerically
integrate the system using the Bogacki-Shampine method \citep{Bogacki1989}. In order to avoid
the computational complexity of performing this numerical integration for many test points,
we compute the value of \eqref{eq:Lyap} with~$K\!=\!1000$ and~$\Delta t\!=\!0.01$ only for~$N\!=\!1024$ initial states
$\bm{x}_0$ uniformly spread over
the rectangle~$\mathbb{X}\!=\![-5,5]^2$ and train a Gaussian process with the data to obtain a learned Lyapunov
function~$\tilde{V}(\cdot)\!=\!\mu(\cdot)$.\looseness=-1

In order to analyze the region of attraction of the system, we follow the approach proposed in
\citep{Lederer2019b}. First, we determine the regions~$\mathbb{W}$
of the state space satisfying the inequality
\begin{align}
	\Delta V(\bm{x}_0)=\tilde{V}(\bm{x}(\Delta t,\bm{x}_0))-\tilde{V}(\bm{x}_0)\leq 0,
\end{align}
which corresponds to verifying \eqref{eq:prob} with~$\epsilon_1=\infty$,~$\epsilon_2=0$,
$g(\cdot)=\mu(\cdot)$ and
$\bm{f}(\cdot)=\bm{x}(\Delta t,\cdot)$. Based on the decrease region~$\mathbb{W}$, the region of attraction
is a level set of the learned Lyapunov given by
\begin{align}
	\mathbb{V}=\{\bm{x}\in\mathbb{X}:\tilde{V}(\bm{x})\leq \min\limits_{\bm{x}\in\mathbb{X}: \Delta \tilde{V}(\bm{x})>0} \tilde{V}(\bm{x})\},
\end{align}
which can also be determined using the GP3 framework. Since the Lipschitz based analysis of
our approach does not allow the verification of the decrease condition close to the origin, we assume
stability in a ball around~$\bm{0}$ with radius~$r=0.1$ similarly as in \citep{Lederer2019b}.\looseness=-1

We apply Alg.~\ref{Alg1} with minimal size~$b_{\min}=10^{-9}$ of the hyperrectangles
and a conservative Lipschitz constant~$L_f=20$ to both problems.
The resulting decrease region~$\mathbb{W}$ as well as region of attraction~$\mathbb{V}$ of the
single machine infinite bus system are illustrated in Fig.~\ref{fig:roa}. Additionally, an estimate
of the region of attraction~$\mathbb{V}_{\mathrm{num}}$ obtained through examining the convergence of
trajectories after~$10^4$ simulation steps for~$4000$ initial states~$\bm{x}_0$ is depicted. The
verified decrease region is
different from the approximated true region of attraction~$\mathbb{V}_{\mathrm{num}}$ at its boundary.
Due to the saw tooth behavior at the boundary of the verified decrease region~$\mathbb{W}$, the resulting
region of attraction~$\mathbb{V}$ is smaller than the numerical approximation~$\mathbb{V}_{\mathrm{num}}$.
However, this underestimation results mainly from imprecision of the learning, while the GP3 approach allows to
analyze the Gaussian process very accurately, since the boundary of~$\mathbb{V}$ visually
touches the boundary of~$\mathbb{W}$ at~$\bm{x}=[1.2\ 1.7]^T$. Therefore, better estimates of
the region of attraction can easily be obtained by training the Gaussian process with more data of the
discrete-time Lyapuonv function \eqref{eq:Lyap}.

\begin{figure}
	\includegraphics[]{plots/decreaseregion}
	\vspace{-0.3cm}
	\caption{Lyapunov decrease region~$\mathbb{W}$ and region of attraction~$\mathbb{V}$
		obtained with Alg.~\ref{Alg1} as well as numerically approximated region of
		attraction~$\mathbb{V}_{\mathrm{num}}$ for the single machine infinite bus system}
	\label{fig:roa}
\end{figure}

\section{Conclusion}
\label{sec:conc}

This paper introduces a novel framework for the analysis of Gaussian process mean
functions called GP3: general purpose computation on graphics processing units for
Gaussian processes. Based on interval analysis to compute local Lipschitz constants,
the posterior mean function is analyzed using multi-resolution sampling. Due independence
of the computations for each sample, the method can be parallelized on a GPU for computational
efficiency. In order to demonstrate the computational benefits of the GP3 framework, it is
applied to a Lipschitz constant bounding and a region of attraction estimation problem.


\section{Introduction}

\label{sec1}

Random matrix theory (RMT) is a powerful mathematical tool for studying
complex quantum systems, it describes the universal properties of random
matrix that are determined only by the system's symmetry while independent
of microscopic details. For this reason, RMT has been applied to various
fields ranging from disordered nuclei to isolated quantum many-body systems%
\cite{Porter,RMP,PR}.

Among various statistical quantities of RMT, the most widely-used ones are
the distribution of nearest level spacings $\left\{
s_{i}=E_{i+1}-E_{i}\right\} $ and gap ratios $\left\{
r_{i}=s_{i+1}/s_{i}\right\} $. It is well established that in a chaotic
system, the distribution of level spacing $P\left( s\right) $ will follow a
Wigner-Dyson distribution\cite{Mehta,Haake2001} (see Eq.~(\ref{equ:WD}) in
Sec.~\ref{sec3}), which reveals the level repulsion in a direct way.
However, when accounting $P\left( s\right) $, an unfolding procedure is
required to erase the model-dependent information about local density of
states (DOS). In contrast, gap ratios distribution $P\left( r\right) $ is
independent of DOS and requires no unfolding procedure\cite{Oganesyan,Atas},
and has found various applications especially in the context of many-body
localization (MBL)\cite%
{Oganesyan,Avishai2002,Regnault16,Regnault162,Huse1,Huse2,Huse3,Sarma,Luitz,Rao182,Sarkar,Corps}%
.

Both the nearest level spacing and gap ratio account for the short-range
level correlations. However, long-range correlations are also important,
especially in the study of MBL transition phenomena. Actually, there're a
number of RMT models accounting for the intermediate level statistics
between Wigner-Dyson and Poisson ensembles\cite{Data,Jain,Auberson,SRPM,Ndawana,Mirlin}, and some of them are suggested
to describe the level statistics in MBL transition regime\cite%
{Shukla,Serbyn,Sierant19,Buijsman,Sierant20}, all of which more or less
describe the short-range level correlations well, while their difference can
only be revealed when long-range correlations are concerned. Commonly, the
long-range correlation in a random matrix ensemble is described by the
number variance $\Sigma ^{2}$ or the Dyson-Mehta $\Delta _{3}$ statistics%
\cite{Haake2001}, however, both of them are very sensitive to the concrete
unfolding strategy which has already been a source of misleading signatures%
\cite{Gomez2002}. Instead, it's more direct and numerically easier to study
the higher order level spacings and gap ratios, as have been done in a
number of recent works\cite%
{Sierant19,Sierant20,Tekur1,Tekur,Bhosale,Atas2,Chavda,Magd,Duras,Rubah,Rao20,Gong}%
.

Formally, the $n$-th order level spacing and gap ratio are defined as $%
\left\{ s_{i}^{\left( n\right) }=E_{i+n}-E_{i}\right\} $ and $\left\{
r_{i}^{\left( n\right) }=s_{i+n}^{\left( n\right) }/s_{i}^{\left( n\right)
}\right\} $ respectively. The higher-order gap ratios are first studied in
Ref.~[\onlinecite{Tekur}], where the authors provide strong numerical
evidences that $P\left( r^{\left( n\right) }\right) $ in a random matrix
ensemble with Dyson index $\beta $ bears the same form as $P\left( r\right) $
with a rescaled parameter
\begin{equation}
\gamma =\frac{n(n+1)}{2}\beta +n-1\text{,}  \label{equ:0}
\end{equation}%
where $\beta =1,2,4$ correspond to the Gaussian orthogonal ensemble (GOE),
Gaussian unitary ensemble (GUE) and Gaussian symplectic ensemble (GSE)
respectively. On the other hand, it is later proved in Ref.~[\onlinecite{Rao20}]
that $P\left( s^{\left( n\right) }\right) $ has the same form as
Wigner-Dyson distribution with a rescaled parameter $\gamma $ that is
identical to the one in Eq.~(\ref{equ:0}) for $\beta \in \left( 0,\infty
\right) $. This strongly hints a homogeneous relation between the
higher-order level spacing and gap ratio, but an explanation is still
lacking, which gives the first motivation of this work.

The second motivation comes from the recent works that encounter physical
systems that go beyond the three standard ensembles with $\beta =1,2,4$. For
example, the $\beta =3$ behavior has been found in a 2D lattice with
non-Hermitian disorder\cite{Tzortzakakis}; and the ensembles with
non-integer $\beta $ is suggested to describe the level statistics in the
whole region along the MBL transition in 1D random spin system\cite{Buijsman}%
, while their efficiency in describing long-range spectral correlations is
controversial\cite{Sierant20}. Therefore, it's beneficial to have an
expression for the higher-order spacing distributions in these ensembles,
which will also be offered in this study.

In this work, we find that the key to link the higher-order level spacing
and gap ratio is the reduced energy spectrum $\left\{ E_{i}^{\left( n\right)
}\equiv E_{in}\right\} $, i.e. the spectrum constructed by picking one level
from every $n$ levels in the original spectrum $\left\{ E_{i}\right\} $. By
this construction, the $n$-th order level spacing and gap ratio in $\left\{
E_{i}\right\} $ become the lowest-order ones in $\left\{ E_{i}^{\left(
n\right) }\right\} $. It will be verified that the joint probability
distribution of $\left\{ E_{i}^{\left( n\right) }\right\} $ (to leading
order) bears the same form as $\left\{ E_{i}\right\} $ with a rescaled
parameter $\gamma $ expressed in Eq.~(\ref{equ:0}), which hence explains
their distributions by virtue of Wigner surmise. Furthermore, this rescaling
relation holds for general $\beta $ beyond $\beta =1,2,4$, therefore, the
higher-order level spacing and gap ratio distributions in these ensembles
can be obtained accordingly, which is thus a natural extension of Ref.~[%
\onlinecite{Tekur}].

This paper is organized as follows. In Sec.~\ref{sec2} we summarize the
formulas regarding the higher-order level spacings and gap ratios, which
motivates the construction of reduced energy spectrum. In Sec.~\ref{sec3} we
focus on the cases that $\left\{ E_{i}^{\left( n\right) }\right\} $ has two
and three levels, and provide compelling numerical evidence for the scaling
of its joint probability distribution in ensembles with general $\beta $. In
Sec.~\ref{sec4} we provide numerical simulations for the distributions of
nearest level spacing and gap ratio in $\left\{ E_{i}^{\left( n\right)
}\right\} $ with large number of levels, and confirm their coincidence with
the higher-order ones in $\left\{ E_{i}\right\} $. In Sec.~\ref{sec5} we
briefly discuss the higher-order reduced energy spectrums. Conclusion and
discussion come in Sec.~\ref{sec6}.

\section{Motivating The Reduced Energy Spectrum}

\label{sec2}

The starting point to study the spectral statistics in the Gaussian ensemble
of random matrix is the joint probability distribution (JPDT) of its energy
levels, whose form is\cite{Mehta,Haake2001}%
\begin{equation}
P\left( \beta ,\left\{ E_{i}\right\} \right) =C\prod_{i<j}\left\vert
E_{i}-E_{j}\right\vert ^{\beta }e^{-A\sum_{i=1}^{N}E_{i}^{2}}\text{,}
\label{equ:Dist0}
\end{equation}%
where the Dyson index $\beta \in \left( 0,\infty \right) $ is a continuous
parameter, $C$ and $A$ are coefficients correlated by the normalization
condition $\int \prod_{i}dE_{i}P\left( \beta ,\left\{ E_{i}\right\} \right)
=1$. It's worth noting that only the subset $\beta =1,2,4$ are physically
invariant ensembles, that is, each matrix element is allowed to drawn from a
Gaussian distribution provided the matrix being invariant under orthogonal,
unitary or symplectic transformations. While for other values of $\beta $,
the JPDT stems from a special tridiagonal random matrix (see Eq.~(\ref%
{equ:Beta}) in Sec.~\ref{sec4}).

From the JPDT in Eq.~(\ref{equ:Dist0}), we can in principle calculate any
statistical quantity we want, in particular the distribution of nearest
level spacings and gap ratios. The general distributions for them in large
dimension $N$ are complicated, while for most practical purpose it is
sufficient to adopt the so-called Wigner surmise that deals with smallest
matrix that holds the quantity of interest. For example, to study nearest
level spacing it's sufficient to consider a $2\times 2$ matrix, which gives
the celebrated Wigner-Dyson distribution\cite{Haake2001}
\begin{equation}
P\left( \beta ,s\right) =C\left( \beta \right) s^{\beta }e^{-A\left( \beta
\right) s^{2}}\text{,}  \label{equ:WD}
\end{equation}%
where the coefficients $C\left( \beta \right) ,A\left( \beta \right) $ are
determined by the normalization conditions%
\begin{equation}
\int_{0}^{\infty }P\left( \beta ,s\right) ds=1\text{, }\int_{0}^{\infty
}sP\left( \beta ,s\right) ds=1\text{.}
\end{equation}%
It's easy to see $P\left( \beta ,s\rightarrow 0\right) \sim s^{\beta }$,
hence $\beta $ is the parameter that controls the strength of level
repulsion.

As for the nearest gap ratios $\left\{ r_{i}=s_{i+1}/s_{i}\right\} $, a
Wigner like surmise is applicable by studying $3\times 3$ matrices\cite{Atas}%
, which gives%
\begin{equation}
P\left( \beta ,r\right) =\frac{1}{Z_{\beta }}\frac{\left( r+r^{2}\right)
^{\beta }}{\left( 1+r+r^{2}\right) ^{1+3\beta /2}}  \label{equ:Pr}
\end{equation}%
where $Z_{\beta }$ is the normalization factor determined by requiring $%
\int_{0}^{\infty }P\left( \beta ,r\right) dr=1$. It is crucial to note that
the derivations for Eq.~(\ref{equ:WD}) and Eq.~(\ref{equ:Pr}) are purely
mathematical, that is, applicable for arbitrary positive $\beta $.

For the higher order level spacings $\left\{ s_{i}^{\left( n\right)
}=E_{i+n}-E_{i}\right\} $, its distribution is studied in Ref.~[%
\onlinecite{Rao20}] using a Wigner-like surmise that deals with $\left(
n+1\right) \times \left( n+1\right) $ matrix, and the result shows they
follow a generalized Wigner-Dyson distribution that bears the same form as
Eq.~(\ref{equ:WD}) with the parameter $\beta $ rescaled to $\gamma $ as
expressed in Eq.~(\ref{equ:0}).

On the other hand, higher order gap ratios come in two different ways, i.e.
the ``overlapping'' way\cite{Atas2} and ``non-overlapping'' way\cite{Tekur}.
In the former case we are dealing with%
\begin{equation}
\widetilde{r}_{i}^{\left( n\right) }=\frac{E_{i+n+1}-E_{i+1}}{E_{i+n}-E_{i}}=%
\frac{s_{i+n}+s_{i+n-1}+...+s_{i+1}}{s_{i+n-1}+s_{i+n-2}+...+s_{i}}\text{,}
\end{equation}%
which is termed ``overlapping'' gap ratio since there are shared spacings
between the numerator and denominator. While the $n$-th order
``non-overlapping'' gap ratio is defined as%
\begin{equation}
r_{i}^{\left( n\right) }=\frac{E_{i+2n}-E_{i+n}}{E_{i+n}-E_{i}}=\frac{%
s_{i+2n-1}+s_{i+2n-2}+...+s_{i+n}}{s_{i+n-1}+s_{i+n-2}+...+s_{i}}\text{.}
\end{equation}%
Both these two generalizations reduce to the nearest gap ratio when $n=1$,
but they are very different when studying their distributions using Wigner
surmise: for non-overlapping ratio $r^{\left( n\right) }$, the smallest
matrix dimension is $\left( 2n+1\right) \times \left( 2n+1\right) $; while
it is $\left( n+2\right) \times \left( n+2\right) $ for overlapping ratios.
Naively, it's expected that $P\left( \widetilde{r}^{\left( n\right) }\right)
$ is more complicated due to the overlapping spacings. Indeed, the
analytical form of $P\left( \widetilde{r}^{\left( 2\right) }\right) $ was
worked out in Ref.~[\onlinecite{Atas2}] and the result is quite involving.
For the non-overlapping ratios, Ref.~[\onlinecite{Tekur}] provides
compelling numerical evidence that in cases with $\beta =1,2,4$, $P\left(
r^{\left( n\right) }\right) $ bears the same form as $P\left( r\right) $
with the same rescaling parameter as higher order level spacing, that is,
Eq.~(\ref{equ:0}).

In summary, we have%
\begin{eqnarray}
P\left( \beta ,s^{\left( n\right) }\right) &=&P\left( \gamma ,s\right) \text{%
,}  \label{equ:Dists} \\
\text{\thinspace }P\left( \beta ,r^{\left( n\right) }\right) &=&P\left(
\gamma ,r\right) \text{,}  \label{equ:Dists2}
\end{eqnarray}%
where $\gamma $ is expressed in Eq.~(\ref{equ:0}). For the rest of this
paper, ``gap ratio'' will always refer to the non-overlapping one.

The identical rescaling behavior for higher level spacing and gap ratio
hints they may be attributed to one single reason, which is found to be the
reduced energy spectrum $\left\{ E_{i}^{\left( n\right) }\right\} $.
Formally, a reduced energy spectrum $\left\{ E_{i}^{\left( n\right)
}\right\} $ is constructed by picking one level from every $n$ levels in the
original spectrum $\left\{ E_{i}\right\} $, which is mathematically achieved
by tracing out every $n-1$ levels in between. This construction is very
similar to that of the reduced density matrix there we trace out the degrees
of freedom in a subsystem, hence $\left\{ E_{i}^{\left( n\right) }\right\} $
is named ``reduced energy spectrum''. By this construction, the $n$-th order
level spacing and gap ratio in $\left\{ E_{i}\right\} $ are mapped to the
lowest-order counterparts in $\left\{ E_{i}^{\left( n\right) }\right\} $.
It's then natural to conjecture that $\left\{ E_{i}\right\} $ and $\left\{
E_{i}^{\left( n\right) }\right\} $ (to leading order) bear the same form for
their probability distributions, with the Dyson index $\beta $ for the
former rescaled to $\gamma $ for the latter according to Eq.~(\ref{equ:0}).
Therefore, by applying Wigner surmise to $\left\{ E_{i}^{\left( n\right)
}\right\} $, the scaling behaviors in Eq.~(\ref{equ:Dists}) and Eq.~(\ref%
{equ:Dists2}) can be explained simultaneously. This is the main task of
current work.

Before proceeding, we want to mention the relatively trivial case of Poisson
ensemble. The reduced energy spectrum in Poisson ensemble has been studied
in Ref.~[\onlinecite{Daisy}] (which is named ``Daisy model'' by the
authors), where $n$-th order level spacing is shown to follow the
generalized semi-Poisson distribution
\begin{equation}
P\left( s^{\left( n\right) }=s\right) =\frac{n^{n}}{(n-1)!}s^{n-1}e^{-ns}%
\text{,}  \label{equ:P1}
\end{equation}%
which reduces to the conventional Poisson distribution $P\left( s\right)
=\exp \left( -s\right) $ when $n=1$. For the $n$-th order gap ratios, we can
derive their distribution from the results in Ref.~[\onlinecite{Daisy}],
that is
\begin{equation}
P\left( r^{(n)}=r\right) =\frac{r^{n-1}}{\left( 1+r\right) ^{2n}}\text{,}
\label{equ:P2}
\end{equation}%
which reduces to the one given in Ref.~[\onlinecite{Atas}] when $n=1$. The
formulas in Eq.~(\ref{equ:Dists}), Eq.~(\ref{equ:Dists2}), Eq.~(\ref{equ:P1}%
) and Eq.~(\ref{equ:P2}) will be used for later numerical simulations.

\section{Scaling of $P\left( \left\{ E_{i}^{\left( n\right) }\right\}
\right) $}

\label{sec3}

\bigskip By the construction of reduced energy spectrum $\left\{
E_{i}^{\left( n\right) }\right\} $, its formal joint probability
distribution is%
\begin{equation}
P\left( \left\{ E_{i}^{\left( n\right) }\right\} \right)
=\prod_{i}\int_{E_{in}}^{E_{\left( i+1\right) n}}\prod_{j=in+1}^{\left(
i+1\right) n-1}dE_{j}P\left( \beta ,\left\{ E_{i}\right\} \right) \text{.}
\label{equ:Dist}
\end{equation}%
For reasons described in previous section, we conjecture $P\left( \left\{
E_{i}^{\left( n\right) }\right\} \right) $ (to leading order) bear the same
form as $P\left( \left\{ E_{i}\right\} \right) $ with the rescaled parameter
$\gamma $ as in Eq.~(\ref{equ:0}), that is,%
\begin{equation}
P\left( \left\{ E_{i}^{\left( n\right) }\right\} \right) \sim
\prod_{i<j}\left\vert E_{i}^{\left( n\right) }-E_{j}^{\left( n\right)
}\right\vert ^{\gamma }e^{-A^{\prime }\sum_{i=1}^{N/n}\left( E_{i}^{\left(
n\right) }\right) ^{2}}\text{.}  \label{equ:Rescale}
\end{equation}%
An analytical derivation from Eq.~(\ref{equ:Dist}) to Eq.~(\ref{equ:Rescale}%
) for arbitrary matrix dimension $N$ is mathematically challenging, only the
special case with $\beta =2/k$ ($k$ being positive integer) has been proven
rigorously in Ref.~[\onlinecite{Forrester}]. It is not our purpose to give a
general proof, instead, to explain the behaviors of the $P\left( s^{\left(
n\right) }\right) $ and $P\left( r^{\left( n\right) }\right) $ in Eq.~(\ref%
{equ:Dists}) and Eq.~(\ref{equ:Dists2}), we only need to verify Eq.~(\ref%
{equ:Rescale}) in the sense of Wigner surmise, that is, in the cases that $%
\left\{ E_{i}^{\left( n\right) }\right\} $ has only two and three levels,
for which we will provide strong numerical evidence in the following.

First of all, the constant $A^{\prime }$ is not important since it is only a
decay rate parameter, whose value can be tuned by global rescaling of the
energy levels without affecting the distribution of level spacing or gap
ratio. Therefore, we will focus on the parameter $\gamma $ that controls the
strength of level repulsion.

The ensembles with general positive $\beta $ can be divided into four
typical categories: (i) $\beta =1,2,4$, corresponding to the three standard
Gaussian ensembles, which are of most physical interest; (ii) $\beta$ is an
integer that goes beyond the three standard ensembles, for which we choose $%
\beta =3$; (iii) $\beta$ is a fraction, for which we choose $\beta =1/3$;
(iv) $\beta$ is an irrational value, for which we choose $\beta =\left(
\sqrt{5}-1\right) /2$ (the golden ratio). We will verify the recaling
relation Eq.~(\ref{equ:Rescale}) in these cases.

We start with the case that $\left\{ E_{i}^{\left( n\right) }\right\} $ has
only two levels, where the rescaling in Eq.~(\ref{equ:Rescale}) becomes%
\begin{eqnarray}
I\left( E_{0},E_{n}\right)
&=&\int_{E_{0}}^{E_{n}}\prod_{i=1}^{n-1}dE_{i}P\left( \beta ,\left\{
E_{i}\right\} \right) \\
&\sim &\left\vert E_{0}-E_{n}\right\vert ^{\gamma }e^{-A^{\prime }\left(
E_{0}^{2}+E_{n}^{2}\right) }\text{.}  \notag
\end{eqnarray}%
Denote $E_{0}=R\cos \theta $ and $E_{n}=R\sin \theta $, and keeping $R$
constant, we reach to%
\begin{equation}
\log I\left( \theta \right) =\gamma \log \left\vert \cos \theta -\sin \theta
\right\vert +\text{const.}  \label{equ:1}
\end{equation}%
Without loss of generality, we take $A=1$ and $R=1$, and randomly generate $%
200$ sets of $\theta $ in the range $[0,2\pi )$. We then numerically
calculate $\log I\left( \theta \right) $ and $\log \left\vert \cos \theta
-\sin \theta \right\vert $, and collect the results for $n=2,3,4$, which are
presented in Fig.~\ref{fig:rescale0}. As can be seen, the $\log I\left(
\theta \right) $ and $\log \left\vert \cos \theta -\sin \theta \right\vert $
shows a perfect linear dependence in all cases, with the fitted values of $%
\gamma $ quite close to the expected ones in Eq.~(\ref{equ:0}), as
summarized in Table~\ref{tab:1}.

\begin{figure*}[htbp]
\centering
\includegraphics[width=2\columnwidth]{figure1.pdf}
\caption{The fitting of Eq.~(\protect\ref{equ:1}) for $\protect\beta=1,2,3,4,%
\frac{1}{3},\frac{\protect\sqrt{5}-1}{2}$ with $n=2,3,4$. The fitted slopes
are shown in the figure legends, where the numbers in the bracket are the
expected values according to Eq.~(\protect\ref{equ:0}). }
\label{fig:rescale0}
\end{figure*}

Next we consider the case that $\left\{ E_{i}^{\left( n\right) }\right\} $
has three levels, now the rescaling in Eq.~(\ref{equ:Rescale}) becomes
\begin{eqnarray}
&&Q\left( E_{0},E_{n},E_{2n}\right)  \notag \\
&\equiv
&\int_{E_{0}}^{E_{n}}\prod_{i=1}^{n-1}dE_{i}\int_{En}^{E_{2n}}%
\prod_{j=n+1}^{2n-1}dE_{j}P\left( \beta ,\left\{ E_{i}\right\} \right) \\
&\sim &\left\vert E_{0}-E_{n}\right\vert ^{\gamma }\left\vert
E_{n}-E_{2n}\right\vert ^{\gamma }\left\vert E_{0}-E_{2n}\right\vert
^{\gamma }e^{-A^{\prime }\left( E_{0}^{2}+E_{n}^{2}+E_{2n}^{2}\right) }\text{%
.}  \notag
\end{eqnarray}%
With the transformation to spherical coordinates%
\begin{equation}
E_{0}=R\sin \theta \cos \varphi ,E_{n}=R\sin \theta \sin \varphi
,E_{2n}=R\cos \theta
\end{equation}%
and keeping $R$ constant, we can reach to%
\begin{equation}
\log Q\left( R,\theta ,\varphi \right) =\gamma \log G\left( \theta ,\varphi
\right) +\text{const.}  \label{equ:2}
\end{equation}%
where%
\begin{eqnarray}
G\left( \theta ,\varphi \right) &=&|\left( \sin \theta \cos \varphi -\sin
\theta \sin \varphi \right)  \notag \\
&&\times \left( \sin \theta \cos \varphi -\cos \theta \right)  \notag \\
&&\times \left( \sin \theta \sin \varphi -\cos \theta \right) |\text{.}
\end{eqnarray}%
We then perform numerical checks, where we fix $R=1$ and $A=1$ as before. We
randomly generate $200$ pairs of $\left( \theta ,\varphi \right) $ and
numerically determine $\log Q\left( \theta ,\varphi \right) \equiv \log
Q\left( 1,\theta ,\varphi \right) $ and $\log G\left( \theta ,\varphi
\right) $, the results are displayed in Fig.~\ref{fig:rescale}. As can be
seen, the linear dependence between $\log Q\left( \theta ,\varphi \right) $
and $\log G\left( \theta ,\varphi \right) $ are still perfect in all cases,
and the fitted values of $\gamma $ are close to the expected ones in Eq.~(%
\ref{equ:0}).

\begin{figure*}[htbp]
\centering
\includegraphics[width=2\columnwidth]{figure2.pdf}
\caption{The fitting of Eq.~(\protect\ref{equ:2}) for $\protect\beta=1,2,3,4,%
\frac{1}{3},\frac{\protect\sqrt{5}-1}{2}$ with $n=2,3,4$. The fitted $%
\protect\gamma$ values are shown in the figure legends, with the expected
values according to Eq.~(\protect\ref{equ:0}) in the brackets.}
\label{fig:rescale}
\end{figure*}

For convenience, we collected the theoretical and numerical values of $%
\gamma $ in Table~\ref{tab:1}, where $\gamma _{\text{e}}$ refers to the
expected value according to Eq.~(\ref{equ:0}), $\gamma _{\text{num.}%
}^{\left( 1\right) }$ and $\gamma _{\text{num.}}^{\left( 2\right) }$ are the
numerical values fitted from Eq.~(\ref{equ:1}) and Eq.~(\ref{equ:2}), and
error$^{\left( 1\right) }$ and error$^{\left( 2\right) }$ are their relative
deviations from $\gamma _{\text{e}}$ respectively. In general, the numerical
errors increases with larger $\beta $, but they're controlled within a
satisfactory level in all cases.

\begin{table}
\centering%
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$\beta $ & $n$ & $\gamma _{\text{e}}$ & $\gamma _{\text{num.}}^{\left(
1\right) }$ & error$^{\left( 1\right) }$ & $\gamma _{\text{num.}}^{\left(
2\right) }$ & error$^{\left( 2\right) }$ \\ \hline
& $2$ & $4$ & $4.10$ & $2.50\%$ & $4.03$ & $0.75\%$ \\ \cline{2-7}
$1$ & $3$ & $8$ & $8.25$ & $3.13\%$ & $8.27$ & $3.38\%$ \\ \cline{2-7}
& $4$ & $13$ & $13.41$ & $3.15\%$ & $13.49$ & $3.77\%$ \\ \hline
& $2$ & $7$ & $7.14$ & $2.00\%$ & $7.20$ & $2.86\%$ \\ \cline{2-7}
$2$ & $3$ & $14$ & $14.29$ & $2.07\%$ & $15.10$ & $7.86\%$ \\ \cline{2-7}
& $4$ & $23$ & $23.37$ & $1.61\%$ & $24.75$ & $7.29\%$ \\ \hline
& $2$ & $13$ & $13.19$ & $1.46\%$ & $13.42$ & $3.23\%$ \\ \cline{2-7}
$4$ & $3$ & $26$ & $26.23$ & $0.88\%$ & $27.66$ & $6.38\%$ \\ \cline{2-7}
& $4$ & $43$ & $43.36$ & $0.84\%$ & $46.48$ & $8.09\%$ \\ \hline
& $2$ & $10$ & $10.16$ & $1.60\%$ & $10.42$ & $4.20\%$ \\ \cline{2-7}
$3$ & $3$ & $20$ & $20.24$ & $1.20\%$ & $21.35$ & $6.75\%$ \\ \cline{2-7}
& $4$ & $33$ & $33.39$ & $1.18\%$ & $35.46$ & $7.45\%$ \\ \hline
& $2$ & $2$ & $2.11$ & $5.5\%$ & $1.91$ & $-4.5\%$ \\ \cline{2-7}
$1/3$ & $3$ & $4$ & $4.20$ & $5.0\%$ & $3.90$ & $-2.5\%$ \\ \cline{2-7}
& $4$ & $19/3$ & $6.71$ & $6.0\%$ & $6.25$ & $-1.3\%$ \\ \hline
& $2$ & $\frac{3\sqrt{5}-1}{2}$ & $2.98$ & $4.4\%$ & $2.82$ & $-1.2\%$ \\
\cline{2-7}
$\frac{\sqrt{5}-1}{2}$ & $3$ & $3\sqrt{5}-1$ & $5.88$ & $3.0\%$ & $5.75$ & $%
0.7\%$ \\ \cline{2-7}
& $4$ & $5\sqrt{5}-2$ & $9.51$ & $3.6\%$ & $9.43$ & $2.7\%$ \\ \hline
\end{tabular}
\caption{The values of $\gamma$ with different $\beta$ and order $n$,
where $\gamma_{\text{e}}$ is the expected value according to Eq.~(\ref{equ:0}).
$\gamma _{\text{num.}}^{\left( 1\right)}$ and $\gamma _{\text{num.}}^{\left( 2\right) }$
are the values fitted from Eq.~(\ref{equ:1}) and Eq.~(\ref{equ:2}),
their relative errors with respect to $\gamma_{\text{e}}$ are denoted by error$^{\left( 1\right) }$
and error$^{\left( 2\right) }$ respectively.} \label{tab:1}
\end{table}

Up to now, we have verified the scaling behavior for the probability
distribution of $\left\{ E_{i}^{\left( n\right) }\right\} $ with two and
three levels, the cases with more levels can be verified with the same
method, but it will become more and more tedious when the number of
integrals increases. Nevertheless, present results are sufficient to justify
the scaling behaviors of higher order level spacing/gap ratio by applying
Wigner surmise to $\left\{ E_{i}^{\left( n\right) }\right\} $. More
importantly, we have verified the scaling behavior to hold for general $%
\beta $s that go beyond GOE,GUE,GSE ($\beta =1,2,4$), even when $\beta $ is
non-integer or irrational. This indicates the distribution for higher order
level spacing and gap ratio in Eq.~(\ref{equ:Dists}) and Eq.~(\ref%
{equ:Dists2}) also hold for general $\beta $ ensembles, for which we will
present numerical evidences in the next section.

\section{Numerical Simulations}

\label{sec4}

\bigskip In this section, we numerically check the distribution of nearest
level spacing/gap ratio in $\left\{ E_{i}^{\left( n\right) }\right\} $ with
many levels, and show they indeed coincide with the $n$-th order
counterparts in $\left\{ E_{i}\right\} $. Before that, a technical issue
needs to be pointed out, that is, the nearest level spacings in $\left\{
E_{i}^{\left( n\right) }\right\} $ is actually $\left\{ E_{i+1}^{\left(
n\right) }-E_{i}^{\left( n\right) }=E_{\left( i+1\right) n}-E_{in}\right\} $%
, whose total number is $\left[ \frac{N}{n}\right] -1$; while the $n$-th
order level spacing in original energy spectrum $\left\{ E_{i}\right\} $ are
$\left\{ E_{i+n}-E_{i}\right\} $ with total number $N-n$; therefore the
mapping does not strictly hold, the same thing happens to gap ratios.
However, since the distribution is extracted from a large number of level
spacings (gap ratios), it's natural to suspect the difference is negligible
when the number of samples and matrix dimension are large, which we will
soon justify.

\begin{figure*}[htbp]
\centering
\includegraphics[width=2\columnwidth]{figure3.pdf}
\caption{Distribution of nearest level spacing and gap ratio in the reduced
energy spectrum $\left\{ E_{i}^{\left( n\right) }\right\}$ of model Eq.~(%
\protect\ref{equ:H}) in a $L=12$ chain with $n=2$ ((a) and (d)), $n=3$ ((b)
and (e)) and $n=4$((c) and (f)). The data from $h=1$ in the orthogonal
(unitary) model represent GOE (GUE), and those from $h=5$ in orthogonal
model represent Poisson respectively. The reference curves corresponds to
the higher-order spacing distributions in $\left\{E_{i}\right\}$ according
to Eq.~(\protect\ref{equ:Dists}), Eq.~(\protect\ref{equ:Dists2}) for GOE and
GUE, and Eq.~(\protect\ref{equ:P1}), Eq.~(\protect\ref{equ:P2}) for Poisson,
the parameter $\protect\gamma$ for the former is calculated by Eq.~(\protect
\ref{equ:0}). The perfect fittings in all cases confirms the coincidence
between higher-order spacing distributions in $\left\{E_{i}\right\}$ and the
lowest-order counterparts in $\left\{E_{i}^{(n)}\right\}$.}
\label{fig:Spin}
\end{figure*}

\bigskip For GOE,GUE and Poisson ensemble, we perform simulations from a
real physical system, that is, the 1D Heisenberg chain with random external
fields, which is the canonical model to study many-body localization\cite%
{Alet}. The Hamiltonian reads%
\begin{equation}
H=J\sum_{i=1}^{L}\mathbf{S}_{i}\cdot \mathbf{S}_{i+1}+\sum_{i=1}^{L}\sum_{%
\alpha =x,y,z}h^{\alpha }\varepsilon _{i}^{\alpha }S_{i}^{\alpha }\text{,}
\label{equ:H}
\end{equation}%
where $\mathbf{S}_{i}$ is spin-$1/2$ operator. The anti-ferromagnetic
coupling strength $J$ is set to be unity, and $\varepsilon _{i}^{\alpha }$s
are random numbers within range $\left[ -1,1\right] $. The $h^{\alpha }$ is
referred as randomness strength. We consider two sets of $h^{\alpha }$: (i) $%
h^{x}=h^{z}=h\neq 0$ and $h^{y}=0$, the model is orthogonal and belongs to
GOE; (ii) $h^{x}=h^{y}=h^{z}=h\neq 0$, the model is unitary and belongs to
GUE. This model undergos a thermal-MBL transition at roughly $h_{c}\simeq 3$
($2.5$) in the orthogonal (unitary) case, where the level spacing
distribution evolves from GOE (GUE) to Poisson\cite{Regnault16,Regnault162}.

We choose a $L=12$ system to perform simulations, and prepare $500$ samples
of energy spectrum at $h=1$ in both the orthogonal and unitary model,
representing GOE and GUE respectively. We also prepare $500$ samples at $h=5$
in the orthogonal model to represent Poisson ensemble. For each sample of
energy spectrum, we manually construct the reduced energy spectrum $\left\{
E_{i}^{\left( n\right) }\right\} $ with $n=2,3,4$ and count the
corresponding nearest level spacing and gap ratio distributions, and compare
them to the formulas in Eq.~(\ref{equ:Dists})-(\ref{equ:P2}), the results
are displayed in Fig.~\ref{fig:Spin}. As can be seen, the fittings are quite
good, confirming the correspondence between nearest level spacing/gap ratio
in $\left\{ E_{i}^{\left( n\right) }\right\} $ with the $n$-th order
counterparts in $\left\{ E_{i}\right\} $.

For ensembles with general $\beta $, we perform numerical simulations from
modelling random matrices. It was proven in Ref.~[\onlinecite{Beta}] that
the eigenvalues of the following tridiagonal matrix ensemble%
\begin{equation}
M_{\beta }=\frac{1}{\sqrt{2}}\left(
\begin{array}{ccccc}
x_{1} & y_{1} &  &  &  \\
y_{1} & x_{2} & y_{2} &  &  \\
&
\begin{array}{ccc}
\text{.} &  &  \\
& \text{.} &  \\
&  & \text{.}%
\end{array}
&
\begin{array}{ccc}
\text{.} &  &  \\
& \text{.} &  \\
&  & \text{.}%
\end{array}
&
\begin{array}{ccc}
\text{.} &  &  \\
& \text{.} &  \\
&  & \text{.}%
\end{array}
&  \\
&  & y_{N-2} & x_{N-1} & y_{N-1} \\
&  &  & y_{N-1} & x_{N}%
\end{array}%
\right)  \label{equ:Beta}
\end{equation}%
will follow the distribution in Eq.~(\ref{equ:Dist0}) with continuous
parameter $\beta \in \left( 0,\infty \right) $ provided the diagonals $%
x_{i}\,$($i=1,2,...,N$) follow the normal distribution $\mathit{N}\left(
0,2\right) $, that is%
\begin{equation}
P\left( x_{i}\right) =\frac{1}{2\sqrt{2\pi }}e^{-x_{i}^{2}/8}\text{, }%
i=1,2,...,N\text{,}
\end{equation}%
and $y_{k}$ ($k=1,2,...,N-1$) follows the $\chi $ distribution with
parameter $\left( N-k\right) \beta $, that is%
\begin{equation}
P\left( y_{k}\right) =\left\{
\begin{array}{ll}
\frac{2}{2^{\left( N-k\right) \beta }\Gamma \left( \left( N-k\right) \beta
/2\right) }y_{k}^{\left( N-k\right) \beta -1}e^{-y_{k}^{2}/2}\text{,} &
y_{k}\geq 0 \\[1mm]
0\text{,} & y_{k}<0%
\end{array}%
\right. \text{.}
\end{equation}%
By virtue this remarkable construction, we can efficiently generate energy
spectrums with any $\beta $ we want\cite{note}.

In accordance with Sec.~\ref{sec3}, we choose to simulate the cases with $%
\beta =3,\frac{1}{3},\frac{\sqrt{5}-1}{2}$. For each $\beta $, we generate $%
500$ samples of energy spectrums by using Eq.~(\ref{equ:Beta}), with the
number of energy levels in the original spectrum $\left\{ E_{i}\right\} $
kept to be $600$. Then we construct the corresponding reduced energy
spectrums $\left\{ E_{i}^{\left( n\right) }\right\} $ with $n=2,3,4$ and
determine the distributions of nearest level spacings/gap ratios in them,
the results are shown in Fig.~\ref{fig:beta}. As can be seen, the fittings
are quite satisfactory in all cases.

\begin{figure*}[tbph]
\centering
\includegraphics[width=2\columnwidth]{figure4.pdf}
\caption{Distribution of nearest level spacing and gap ratio in the reduced
energy spectrum $\left\{ E_{i}^{\left( n\right) }\right\}$ with $n=2,3,4$
from ensembles with $\protect\beta=3$ ((a),(d)), $\frac{1}{3}$ ((b),(e)) and
$\frac{\protect\sqrt{5}-1}{2}$ ((c),(f)), and the reference curves with
index $\protect\gamma$ correspond to the ones in Eq.~(\protect\ref{equ:Dists}%
) and Eq.~(\protect\ref{equ:Dists2}) for level spacing and gap ratio
respectively. These fittings confirm the higher-order spacing distributions
in general $\protect\beta$ ensembles.}
\label{fig:beta}
\end{figure*}

Up to now, we have verified the coincidence between nearest level
spacing/gap ratio in $\left\{ E_{i}^{\left( n\right) }\right\} $ with the $n$%
-th order ones in $\left\{ E_{i}\right\} $\ for general Gaussian $\beta $
ensembles. As a result, we have generalized the scaling behavior for $%
P\left( s^{\left( n\right) }\right) $ and $P\left( r^{\left( n\right)
}\right) $ in Eq.~(\ref{equ:Dists}) and Eq.~(\ref{equ:Dists2}) to cases
beyond GOE,GUE,GSE. The Gaussian ensemble with non-integer $\beta $ has been
used to describe the spacing distributions along the thermal-MBL transition%
\cite{Buijsman}, while its efficiency in describing long-range level
correlations is under debate\cite{Sierant19}. Our results thus provide a
numerical criteria for such studies.

\section{Higher-order Reduced Energy Spectrums}

\label{sec5}

\bigskip Results in Sec.~\ref{sec3} verified the scaling of $\left\{
E_{i}^{\left( n\right) }\right\} $ in cases with two/three levels, which is
sufficient for applying Wigner surmise to $\left\{ E_{i}^{\left( n\right)
}\right\} $ regarding level spacing/gap ratio, and the numerical simulations
in Sec.~\ref{sec4} provide strong support for the cases where $\left\{
E_{i}^{\left( n\right) }\right\} $ has more levels. Actually, we can
continue to construct the $2$nd order reduced energy spectrums, i.e. the
``reduced energy spectrum of reduced energy spectrum'', whose nearest level
spacing/gap ratio will correspond to the higher-order ones in $\left\{
E_{i}^{\left( n\right) }\right\} $, and will rescale in a similar manner.
Denote the $m$-th order level spacing and gap ratio in $\left\{
E_{i}^{\left( n\right) }\right\} $ as $s_{i}^{\left( n,m\right) }$ and $%
r_{i}^{\left( n,m\right) }$, that is,%
\begin{equation}
s_{i}^{\left( n,m\right) }=E_{i+m}^{\left( n\right) }-E_{i}^{\left( n\right)
}\text{, }r_{i}^{\left( n,m\right) }=\frac{E_{i+2m}^{\left( n\right)
}-E_{i+m}^{\left( n\right) }}{E_{i+m}^{\left( n\right) }-E_{i}^{\left(
n\right) }}\text{,}
\end{equation}%
it's straightforward to write down their expected probability distributions%
\begin{eqnarray}
&&P\left( \beta ,s^{\left( n,m\right) }\right) =P\left( \delta ,s\right)
\text{, }P\left( \beta ,r^{\left( n,m\right) }\right) =P\left( \delta
,r\right) \\
&&\delta =\frac{m(m+1)}{2}\gamma +m-1\text{, }\gamma =\frac{n(n+1)}{2}\beta
+n-1\text{.}  \notag
\end{eqnarray}%
The same procedure can continue for higher-order reduced energy spectrums.

Of course, such a construction is artificial, meanwhile, it reveals a
hierarchy of energy spectrums can emerge from single spectrum, which (to
lowest order) bear the same form of probability distributions. Moreover, by
taking a closer look at the scaling expression $\gamma =\frac{n\left(
n+1\right) }{2}\beta +n-1$, we immediately recognize an infinite number of
coincident relations between $\left\{ E_{i}^{\left( n\right) }\right\} $
from different ensembles. For examples, $\left\{ E_{i}^{\left( 2\right)
}\right\} $ in $\beta =1$ (GOE) has the same structure as $\left\{
E_{i}\right\} $ in $\beta =4$ (GSE) -- a result known before\cite%
{GSE,Forrester}, and both of them coincide with $\left\{ E_{i}^{\left(
3\right) }\right\} $ in $\beta =\frac{1}{3}$. Actually, it's easy to verify
that for $\gamma \in (k,k+1]$ there exists $k$ different sets of $\left(
\beta ,n\right) $ that have equal $\gamma $, and their lower-order level
statistics are expected to be identical.

\section{Conclusion and Discussion}

\label{sec6}

We studied the reduced energy spectrum $\left\{ E_{i}^{\left( n\right)
}\right\} $, constructed by picking one level from every $n$ levels in
original spectrum $\left\{ E_{i}\right\} $. It is verified the distribution
of $\left\{ E_{i}^{\left( n\right) }\right\} $ (to leading order) bears the
same form as $\left\{ E_{i}\right\} $, with the Dyson index rescaled from $%
\beta $ to $\gamma =\frac{n\left( n+1\right) }{2}\beta +n-1$. It's then
demonstrated the nearest level spacing and gap ratio in $\left\{
E_{i}^{\left( n\right) }\right\} $ correspond to the $n$-th order ones in $%
\left\{ E_{i}\right\} $, which explains the distributions of the latter
found recently in Ref.~[\onlinecite{Rao20}] and Ref.~[\onlinecite{Tekur}]
simultaneously.

Moreover, we find the rescaling of reduced energy spectrum holds for
Gaussian ensembles that go beyond the standard GOE,GUE,GSE, and establish
the distributions of higher-order level spacings and gap ratios in these
ensembles. We also confirmed such correspondences in the Poisson ensemble,
and discovered the distribution of $n$-th order gap ratios, as expressed in
Eq.~(\ref{equ:P2}).

It's noted the reduced energy spectrum has been studied for the Poisson
ensemble\cite{Daisy} and some Gaussian ensembles, that is, the special ones
with $\beta =2/k$ ($k$ being positive integer)\cite{Forrester}, which
contains the well-known coincidence between $\left\{ E_{i}^{\left( 2\right)
}\right\} $ in GOE and $\left\{ E_{i}\right\} $ in GSE. Our work is thus a
natural extension of these studies.

\bigskip The significance of our work is three-folded. First, we explained
the distributions of higher-order level spacings and gap ratios -- both of
which are widely used in the study of MBL and whose distributions are found
separately in recent studies -- by a single common mechanism: the reduced
energy spectrum. Second, we generalized the higher-order spacing
distributions in Eq.~(\ref{equ:Dists}) and Eq.~(\ref{equ:Dists2}) to general
$\beta $ ensembles, which may be beneficial for studying systems that go
beyond the standard Gaussian ensembles\cite{Lucas,Tzortzakakis,Buijsman,Ahmed,Jain2}. Third, our results reveals a rich set of
structures hidden in the energy spectrum, by constructing the reduced energy
spectrums.

\bigskip Last but not least, in our numerical simulations in Sec.~\ref{sec4}%
, we were employing the modelling matrix\ of general $\beta $ ensemble as
expressed in Eq.~(\ref{equ:Beta}). It's thus natural and interesting to ask
whether this ``parent matrix'' corresponds to a real quantum system, and
what's the property of such a system if it does. Given the physically
relevant GOE,GUE,GSE are incorporated in Eq.~(\ref{equ:Beta}), we conjecture
such a ``parent Hamiltonian'' does exist, whose construction is left for a
future study.

\section*{Acknowledgements}

This work is supported by the National Natural Science Foundation of China
through Grant No.11904069 and No.11847005 and No.11804070.