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

Inference on a high-dimensional joint distribution based on a relatively small number of replicates is important in many applications, such as in ensemble-based data assimilation and in the analysis of geospatial data and functional computer-model output.
An important example is statistical climate-model emulation \citep[e.g.,][]{Castruccio2014,Nychka2018,Haugen2019}, whose goal is to learn a high-dimensional, nonstationary, and non-Gaussian spatial distribution from an ensemble of spatial fields generated by an expensive computer model (see Figure \ref{fig:precdata}).

\begin{figure}
\centering
\begin{minipage}{.94\textwidth}
\centering
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_1.png}
	\end{subfigure}%
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_2.png}
	\end{subfigure}%
	\\
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_3.png}
	\end{subfigure}%
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_4.png}
	\end{subfigure}%
\end{minipage}
\hfill
\begin{minipage}{.05\textwidth}
\includegraphics[trim=400mm 0mm 0mm 0mm, clip,width =.98\linewidth]{precip_global_legend.png}
\end{minipage}
  \caption{Four members of an ensemble of log-transformed precipitation anomalies produced by a climate model, on a global grid of size $N = 288 \times 192 = 55{,}296$. We want to infer the $N$-dimensional distribution based on an ensemble of $n<100$ training samples. See Section \ref{sec:application} for details.}
\label{fig:precdata}
\end{figure}

Existing methods for spatial inference are often geared toward a single training replicate and assume Gaussian processes (GPs) with simple parametric covariance functions \citep[e.g.,][]{Cressie1993,Banerjee2004}. Many extensions to nonstationary \citep[e.g., as reviewed by][]{Risser2016} or nonparametric covariances \citep[e.g.,][]{Huang2011,Choi2013,Porcu2019} have been proposed, but these typically still rely on implicit or explicit assumptions of Gaussianity.
For non-Gaussian spatial data, GPs can be transformed or used as latent building blocks \citep[see, e.g.,][and references therein]{Gelfand2016,Xu2017}, but relying on a GP's covariance function limits the types of dependence that can be captured. Parametric non-Gaussian Mat\'ern fields can be constructed using stochastic partial differential equations driven by non-Gaussian noise \citep{Wallin2015,Bolin2020}.
Models for non-Gaussian spatial data can also be built using copulas; for example, \citet{Graler2014} proposed vine copulas for spatial fields with extremal behavior, and the factor copula approach of \citet{Krupskii2018} assumes all locations in a homogeneous spatial region to be affected by a common latent factor.
Many existing non-Gaussian spatial methods are not scalable to large datasets without approximations.

A convenient and flexible way to model a continuous multivariate distribution with any dependence structure is via a triangular transport map \citep[see][for a review]{Marzouk2016} that transforms the target distribution to a reference distribution (e.g., standard Gaussian). For Gaussian target distributions, such a map is linear and given by the Cholesky factor of the precision matrix; non-Gaussian distributions can be obtained by allowing nonlinearities in the map. Given an invertible transport map, it is straightforward to sample from the target distribution and some of its conditionals, or to transform the non-Gaussian data to the reference space, in which simple linear operations (e.g., averaging) may be more meaningful.
Typically, the map is estimated based on training data, often by iteratively expanding a finite-dimensional parameterization of the transport map \citep[e.g.,][]{ElMoselhy2012,Bigoni2016,Marzouk2016,baptista2020adaptive}; subsequent inference is then carried out assuming that the map is known.

Here we propose Bayesian inference on a transport map that describes a multivariate continuous distribution and is learned from a limited number of samples from the distribution.
We model the map components using nonparametric, conjugate Gaussian-process priors, which probabilistically regularize the map and shrink toward linearity. The resulting procedure is flexible, naturally quantifies uncertainty, and adjusts to the amount of complexity that is discernible from the training data, thus avoiding both over- and under-fitting. The conjugacy results in simple, closed-form inference.
Instead of assuming Gaussianity for the multivariate target distribution, our approach is equivalent to a series of conditional GP regression problems that together characterize a non-Gaussian target distribution.

We then focus on the case of a spatial target distribution, for which the dimension (i.e., the number of spatial locations) is high and the number of training replicates is relatively small.
In this setting, our Bayesian transport maps impose sparsity and regularization motivated by the behavior of diffusion-type processes that are encountered in many environmental applications.
After applying a so-called maximin ordering of the spatial locations, determining the triangular transport map essentially consists of conditional spatial-prediction problems on an increasingly fine scale. We discuss how this scale decay results in conditional near-Gaussianity for a large class of non-Gaussian stochastic processes associated with quasilinear partial differential equations. Hence, our prior distributions are motivated by the behavior of Gaussian fields with Mat\'ern-type covariance, for which the so-called screening effect leads to a decay of influence that allows us to only consider nearby observations in the spatial prediction problems.
The resulting Bayesian methods scale near-linearly in the number of spatial locations, and the main computations are trivially parallel.
Hyperparameters, which are inferred from the data, determine the amount of shrinkage and sparsity.
As a special case, the hyperparameters can be set such that the transport map is forced to be linear and is essentially equivalent to a sparse inverse Cholesky factor of the precision matrix; the joint distribution then becomes Gaussian, and we obtain a model similar to that in \citet{Kidd2020}, who proposed flexible Bayesian methods for inferring the covariance matrix of a Gaussian spatial field based on a sparse Cholesky factor of the precision matrix.

We extend our approach to further increase the flexibility in the (continuous) marginal distributions by modeling the GP-regression error terms using Dirichlet process mixtures, which can be fit using a Gibbs sampler. The resulting method lets the data decide the degrees of nonlinearity, nonstationary, and non-Gaussianity, with little to no manual tuning or model-selection required.
We also discuss an extension for settings in which Euclidean distance between the locations is not meaningful or in which variables are not identified by spatial locations (e.g., multivariate spatial processes).

The remainder of this document is organized as follows. In Section \ref{sec:bayesian}, we develop Bayesian transport maps. In Section \ref{sec:spatial}, we consider the special case of high-dimensional spatial distributions. In Section \ref{sec:dpm}, we discuss extensions to non-Gaussian marginals using Dirichlet process mixtures. Sections \ref{sec:simstudy} and \ref{sec:application} provide comparisons and applications to simulated data and climate-model output, respectively. Section \ref{sec:conclusions} concludes and discusses future work. Appendices \ref{app:proofs}--\ref{app:dpmgibbs} contain proofs and further details. \texttt{R} code to run our method and reproduce all results is available at \url{https://github.com/katzfuss-group/BaTraMaSpa}.

\section{Bayesian transport maps\label{sec:bayesian}}

\subsection{Transport maps as regressions\label{sec:regression}}

Consider a continuous random vector $\mathbf{y} = (y_1,\ldots,y_N)^\top$, for example describing a spatial field at $N$ locations as in Figure \ref{fig:precdata}. For simplicity, assume that $\mathbf{y}$ has been centered to have mean zero.

For a multivariate Gaussian distribution, $\mathbf{y} \sim \mathcal{N}_N(\mathbf{0},\bm{\Sigma})$ with $\bm{\Sigma}^{-1} = \mathbf{L}^\top\mathbf{L}$, the (lower-triangular) Cholesky factor $\mathbf{L}$ represents a transformation to a standard normal: $\mathbf{z} = \mathbf{L} \mathbf{y} \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$. As a natural extension, we can characterize any continuous $N$-variate distribution $p(\mathbf{y})$ by a potentially nonlinear transport map $\mathcal{T}: \mathbb{R}^N \rightarrow \mathbb{R}^N$ \citep[][]{Villani2009}, such that $\mathbf{z} = \mathcal{T}(\mathbf{y}) \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$ for $\mathbf{y} \sim p(\mathbf{y})$.
Like $\mathbf{L}$, we can assume without loss of generality that the transport map $\mathcal{T}$ is lower-triangular \citep{Rosenblatt1952,Carlier2009},
\begin{equation}
\mathcal{T}(\mathbf{y}) = \begin{bmatrix*}[l] \mathcal{T}_1(y_1) \\ \mathcal{T}_2(y_1,y_2) \\ ~\,\vdots \\ \mathcal{T}_N(y_1,y_2,\ldots,y_N) \end{bmatrix*}, \label{eq:map}
\end{equation}
where each $\mathcal{T}_i(\mathbf{y}_{1:i})$ with $\mathbf{y}_{1:i} = (y_1,\ldots,y_i)^\top$ is an increasing function of its $i$th argument.
Letting $\mathcal{N}(x|\mu,\sigma^2)$ denote a Gaussian density with parameters $\mu$ and $\sigma^2$ evaluated at $x$, we then have
\begin{equation}
    \label{eq:mapdist}
\textstyle p(\mathbf{y}) = p_\mathbf{z}\big(\mathcal{T}(\mathbf{y})\big) \, |\text{det} \nabla \mathcal{T}| = \prod_{i=1}^N \big( \, \mathcal{N}(\mathcal{T}_i(\mathbf{y}_{1:i})|0,1) \, \left|\frac{\partial \mathcal{T}_i(\mathbf{y}_{1:i})}{\partial y_i}\right| \big),
\end{equation}
as the triangular $\mathcal{T}$ also implies a triangular $\nabla \mathcal{T} = ( \frac{\partial \mathcal{T}_i(\mathbf{y}_{1:i})}{\partial y_j} )_{i,j=1,\ldots,N}$.

Throughout, we assume each $\mathcal{T}_i$ to be additive in its $i$th argument,
\begin{equation}
    \label{eq:condmap}
\mathcal{T}_i(\mathbf{y}_{1:i}) = (y_i - f_i(\mathbf{y}_{1:i-1}))/d_i, \qquad i=1,\ldots,N,
\end{equation}
for some $d_i \in \mathbb{R}^+$, $f_i: \mathbb{R}^{i-1} \rightarrow \mathbb{R}$ for $i=2,\ldots,n$, and $f_i(\mathbf{y}_{1:i-1}) \equiv 0$ for $i=1$.
Then, $\partial_i \mathcal{T}_i(\mathbf{y}_{1:i}) = 1/d_i > 0$, as required.
Using \eqref{eq:mapdist}, it is easy to show that
\begin{equation}
\label{eq:map2norm}
   \textstyle p(\mathbf{y}) \propto  \prod_{i=1}^N \big( \exp(-\frac{1}{2d_i^2}(y_i - f_i(\mathbf{y}_{1:i-1}))^2) \, \frac{1}{d_i} \big) \propto \prod_{i=1}^N \mathcal{N}(y_i|f_i(\mathbf{y}_{1:i-1}),d_i^2).
\end{equation}
Thus, the transport-map approach has turned the difficult problem of inferring the $N$-variate distribution of $\mathbf{y}$ into $N$ independent regressions of $y_i$ on $\mathbf{y}_{1:i-1}$ of the form
\begin{equation}
y_i = f_i(\mathbf{y}_{1:i-1}) + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0,d_i^2), \qquad i=1,\ldots,N.  \label{eq:regi}
\end{equation}

\subsection{Modeling the map functions using Gaussian processes\label{sec:priors}}

In the existing transport-map literature \citep[e.g.,][]{Marzouk2016}, $f_i: \mathbb{R}^{i-1} \rightarrow \mathbb{R}$ and $d_i \in \mathbb{R}^+$ in \eqref{eq:condmap}--\eqref{eq:regi} are often assumed to have parametric form, whose parameters are estimated and then assumed known.
Instead, we here assume a flexible, nonparametric prior on the map $\mathcal{T}$ by specifying independent conjugate Gaussian-process-inverse-Gamma priors for the $f_i$ and $d_i^2$. These prior assumptions induce prior distributions on the map components $\mathcal{T}_i$ in \eqref{eq:condmap}, and thus on the entire map $\mathcal{T}$ in \eqref{eq:map}.
Our Bayesian model is most appropriate when the $y_i$ are or have been transformed to be marginally Gaussian, but this will be relaxed in Section \ref{sec:dpm}.

Specifically, for the ``noise'' variances $d_i^2$, we assume inverse-Gamma distributions,
\begin{equation}
\label{eq:dprior}
    d_i^2 \stackrel{ind.}{\sim} \mathcal{IG}(\alpha_i,\beta_i), \qquad i=1,\ldots,N.
\end{equation}
Conditional on $d_i^2$, each function $f_i$ is modeled as a Gaussian process (GP).
We implicitly assume a linear mean for the GP, but by using a Gaussian prior distribution
for the linear coefficients,
we can absorb the mean into the covariance function. This results in a zero-mean GP,
\begin{equation}
    \label{eq:gp}
f_i | d_i \stackrel{ind.}{\sim} \mathcal{GP}(0,d_i^2 K_i), \qquad i=1,\ldots,N,
\end{equation}
where we multiply the kernel function by $d_i^2$ for conjugacy and then counteract by dividing  $K_i(\cdot,\cdot) = C_i(\cdot,\cdot)/E(d_i^2)$ by $E(d_i^2) = \beta_i/(\alpha_i -1)$. Due to absorbing the GP mean into the kernel, the covariance function
\begin{equation}
    \label{eq:kernel}
C_i(\mathbf{y}_{1:i-1},\mathbf{y}_{1:i-1}') = \mathbf{y}_{1:i-1}^\top \mathbf{Q}_i \mathbf{y}_{1:i-1}' + \sigma^2_i \, \rho_i(\mathbf{y}_{1:i-1},\mathbf{y}_{1:i-1}'), \qquad i=1,\ldots,N,
\end{equation}
consists of a component corresponding to the linear mean and of a nonlinear component with a positive-definite correlation function $\rho_i$, such that $\rho_i(\mathbf{y}_{1:i-1},\mathbf{y}_{1:i-1}) =1$. The factor $\sigma_i^2$ determines the degree of nonlinearity; if $\sigma_i^2=0$, then $f_i$ becomes a linear function of $\mathbf{y}_{1:i-1}$.

The prior distributions (i.e., $\alpha_i$, $\beta_i$, $K_i$) may depend on hyperparameters $\bm{\theta}$; see Section \ref{sec:hyper} for more details.

\subsection{The posterior map\label{sec:invertiblemap}}

Now assume that we have observed $n$ independent training samples $\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(n)}$ from the distribution in Section \ref{sec:regression} conditional on $\bf=(f_1,\ldots,f_N)$ and $\mathbf{d}=(d_1,\ldots,d_N)$, such that
$\mathbf{y}^{(j)} \stackrel{iid}{\sim} p(\mathbf{y} | \bf,\mathbf{d})$ with
$\mathcal{T}(\mathbf{y}^{(j)}) \, | \, \bf,\mathbf{d} \, \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$, $j=1,\ldots,n$.
We combine the samples into a $n \times N$ data matrix $\mathbf{Y}$ whose $j$th row is given by $\mathbf{y}^{(j)}$. Then, for the regression in \eqref{eq:regi}, the responses $\mathbf{y}_i$ and the covariates $\mathbf{Y}_{1:i-1}$ are given by the $i$th and the first $i-1$ columns of $\mathbf{Y}$, respectively. Below, let $\mathbf{y}^\star$ denote a new observation sampled from the same distribution, $\mathbf{y}^\star \sim p(\mathbf{y} | \bf,\mathbf{d})$, independently of $\mathbf{Y}$.

Based on the prior distribution for $\bf$ and $\mathbf{d}$ in Section \ref{sec:priors}, we can now determine the posterior map $\widetilde{\mathcal{T}}$ learned from the training data $\mathbf{Y}$, with $\bf$ and $\mathbf{d}$ integrated out. This map is available in closed form and invertible:
\begin{proposition}
\label{prop:maps}
The transport map $\widetilde{\mathcal{T}}$ from $\mathbf{y}^\star \sim p(\mathbf{y}|\mathbf{Y})$ to
$\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star) \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$ is a triangular map with components
\begin{equation}
z_i^\star = \widetilde{\mathcal{T}}_i(y_1^\star,\ldots,y_i^\star) = \Phi^{-1}\big( F_{2\tilde\alpha_i}\big( \hat d_i^{-1} (v_i(\mathbf{y}^\star_{1:i-1})+1)^{-1/2}(y_i^\star - \hat f_i(\mathbf{y}^\star_{1:i-1})) \big)\big), \quad i=1,\ldots,N, \label{eq:singlemap}
\end{equation}
where
$\tilde\alpha_i = \alpha_i + n/2$,
$\tilde\beta_i = \beta_i + \mathbf{y}_i{}^\top \mathbf{G}_i^{-1} \mathbf{y}_i/2$,
$\hat d_i^2=\tilde\beta_i/\tilde\alpha_i$,
$\mathbf{G}_i = \mathbf{K}_i + \mathbf{I}_n$,
\begin{align}
\hat f_i(\mathbf{y}^\star_{1:i-1}) & = K_i(\mathbf{y}^\star_{1:i-1},\mathbf{Y}_{1:i-1})\mathbf{G}_i^{-1}\mathbf{y}_i, \label{eq:gppredmean}\\
v_i(\mathbf{y}^\star_{1:i-1}) & = K_i(\mathbf{y}^\star_{1:i-1},\mathbf{y}^\star_{1:i-1}) - K_i(\mathbf{y}^\star_{1:i-1},\mathbf{Y}_{1:i-1})\mathbf{G}_i^{-1} K_i(\mathbf{Y}_{1:i-1},\mathbf{y}^\star_{1:i-1}), \label{eq:gppredvar}
\end{align}
for $i=2,\ldots,N$, $\hat f_1 = v_1 = 0$ for $i=1$, and $\Phi$ and $F_{\kappa}$ denote the cumulative distribution functions of the standard normal and the $t$ distribution with $\kappa$ degrees of freedom, respectively.
The inverse map $\widetilde{\mathcal{T}}^{-1}$ can be evaluated at a given $\mathbf{z}^\star$ by solving the nonlinear triangular system $\widetilde{\mathcal{T}}(\mathbf{y}^\star) =\mathbf{z}^\star$ for $\mathbf{y}^\star$; because $\widetilde{\mathcal{T}}$ is triangular, the solution can be expressed recursively as:
\begin{equation}
y_i^\star = \hat f_i(\mathbf{y}_{1:i-1}^\star) + F_{2\tilde\alpha_i}^{-1}(\Phi(z_i^\star))\, \hat d_i (v_i(\mathbf{y}_{1:i-1}^\star)+1)^{1/2}, \quad i=1,\ldots,N. \label{eq:invmap}
\end{equation}
\end{proposition}
All proofs are provided in Appendix \ref{app:proofs}. Note that we can also write the prior map in a similar form, but the prior map will only be useful in the case of highly informative priors.

Determining $\widetilde{\mathcal{T}}_i$ requires $\mathcal{O}(n^3 + i n^2)$ time, largely for computing and decomposing the $n \times n$ matrix $\mathbf{G}_i$, for each $i=1,\ldots,N$. However, note that the $N$ rows or components of $\widetilde{\mathcal{T}}$ can be computed completely in parallel, as in the optimization-based transport-map estimation reviewed in \citet{Marzouk2016}. Each application of the transport map or its inverse then consists of the GP prediction in \eqref{eq:gppredmean}--\eqref{eq:gppredvar} and only requires $\mathcal{O}(n^2 + in)$ time for $i=1,\ldots,N$, but the inverse map is evaluated recursively (i.e., not in parallel).

In contrast to previous transport-map approaches, our approach is Bayesian and naturally quantifies uncertainty in the nonlinear transport functions. The GP priors on the $f_i$ automatically adapt to the amount of information available, only resulting in strongly nonlinear function estimates when supplied the requisite evidence by the data.
However, we can also consider some simple special cases.
If $n$ is increasing, then $\tilde\alpha_i$ increases, $F_{2\tilde\alpha_i}$ converges to $\Phi$, and $v_i(\mathbf{y}^\star_{1:i-1})$ will typically converge to zero, and so the map components will simplify to
\begin{equation}
\label{eq:detmap}
\widetilde{\mathcal{T}}_i(y_1^\star,\ldots,y_i^\star) = (y_i^\star - \hat f_i(\mathbf{y}_{1:i-1}^\star))/\hat d_i  \quad \text{and} \quad  y_i^\star = \hat f_i(\mathbf{y}_{1:i-1}^\star) + z_i \hat d_i.
\end{equation}
When employed for finite $n$, this simplified version of the map ignores posterior uncertainty in $\bf$ and $\mathbf{d}$ and simply relies on the point estimates $\hat f_i(\mathbf{y}_{1:i-1})$ and $\hat d_i^2$. For this simplified map, if we also assume that $\sigma_i=0$ in \eqref{eq:kernel} for all $i=1,\ldots,N$, then all $f_i$ and all $\widetilde{\mathcal{T}}_i$ will become linear functions; we can think of the resulting linear map $\widetilde{\mathcal{T}}(\mathbf{y}^\star) = \mathbf{L}^\top \mathbf{y}^\star$ as an inverse Cholesky factor, in the sense that $\mathbf{y}^\star |\mathbf{Y} \sim \mathcal{N}(\mathbf{0},\bm{\Lambda}^{-1})$ with $\bm{\Lambda} = \mathbf{L}\bL^\top$.

Transport maps can be used for a variety of purposes. For example, we can obtain a new sample $\mathbf{y}^\star$ from the posterior predictive distribution $p(\mathbf{y}|\mathbf{Y})$ by sampling $\mathbf{z}^\star \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$ and computing $\mathbf{y}^\star=\widetilde{\mathcal{T}}^{-1}(\mathbf{z}^\star)$ using \eqref{eq:invmap}.
The map $\widetilde{\mathcal{T}}$ in \eqref{eq:singlemap} provides a transformation from a non-Gaussian vector $\mathbf{y}^\star$ to the standard Gaussian $\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$; we call $\mathbf{z}^\star=(z_1^\star,\ldots,z_N^\star)$ the map coefficients corresponding to $\mathbf{y}^\star$.
Because the nonlinear dependencies have been removed, many operations are more meaningful on $\mathbf{z}^\star$ than on $\mathbf{y}^\star$, including linear regressions, translations using linear shifts, and quantifying similarity using inner products. We can also detect inadequacies of the map $\widetilde{\mathcal{T}}$ for describing the target distribution by examining the degree of non-Gaussianity and dependence in $\mathbf{z}^\star$.
These uses of transport maps will be considered further in Section \ref{sec:spatialinference}.

A close relative of transport maps in machine learning are normalizing flows \citep[see][for a review]{kobyzev2020normalizing}, where triangular layers are used to ensure easy evaluation and inversion of likelihood objectives.
\cite{Kovachki2020} also design generative adversarial networks with triangular generators that allow easy conditional sampling.
Our approach can be viewed as a Bayesian autoencoder \citep[e.g.,][Ch.~14]{Goodfellow2016}, with the posterior transport map and its inverse acting as the encoder and decoder, respectively.

\subsection{Hyperparameters\label{sec:hyper}}

The prior distributions on the $f_i$ and $d_i$ in Section \ref{sec:priors}
may depend on unknown hyperparameters $\bm{\theta}$. For example, by making inference on hyperparameters in the $\sigma_i$ in \eqref{eq:kernel}, we can let the data decide the degree of nonlinearity in the map and thus the non-Gaussianity in the resulting joint target distribution.
We can write in closed form the integrated likelihood $p(\mathbf{Y})$, where $\bf$ and $\mathbf{d}$ have been integrated out.

\begin{proposition}
\label{prop:lik}
The integrated likelihood is
\begin{equation}
    p(\mathbf{Y}) \textstyle \propto \prod_{i=1}^N \big( \, |\mathbf{G}_i|^{-1/2} \times ({\beta_i^{\alpha_i}}/{\tilde\beta_i^{\tilde\alpha_i}}) \times {\Gamma(\tilde\alpha_i)}/{\Gamma(\alpha_i)} \, \big),
    \label{eq:intlik}
\end{equation}
where $\Gamma(\cdot)$ denotes the gamma function, and $\tilde\alpha_i$, $\tilde\beta_i$, $\mathbf{G}_i$ are defined in Proposition \ref{prop:maps}.
\end{proposition}

Now denote by $p_{\bm{\theta}}(\mathbf{Y})$ the integrated likelihood $p(\mathbf{Y})$ computed based on a particular value $\bm{\theta}$ of the hyperparameters. There are two main possibilities for inference on $\bm{\theta}$.
First, an empirical Bayesian approach consists of estimating $\bm{\theta}$ by the value that maximizes $p_{\bm{\theta}}(\mathbf{Y})$, and then regarding $\bm{\theta}$ as fixed and known.
Second, we can carry out fully Bayesian inference by specifying a prior $p(\bm{\theta})$, and sampling $\bm{\theta}$ from its posterior distribution $p(\bm{\theta}|\mathbf{Y}) \propto p_{\bm{\theta}}(\mathbf{Y}) p(\bm{\theta})$ using Metropolis-Hastings; subsequent inference then relies on these posterior draws.

For all our numerical results below, we employed the empirical Bayesian approach, because it is faster and preserves the closed-form map properties in Section \ref{sec:invertiblemap}. In exploratory numerical experiments, we observed no significant decrease in inferential accuracy relative to the fully Bayesian approach, likely due to working with a relatively small number of hyperparameters in $\bm{\theta}$.

\section{Bayesian transport maps for large spatial fields\label{sec:spatial}}

Now assume that $\mathbf{y} = (y_1,\ldots,y_N)^\top$ consists of observations $y_i = y(\mathbf{s}_i)$ of a spatial field $y(\cdot)$ at locations $\mathbf{s}_1,\ldots,\mathbf{s}_N$ in a region or domain $\mathcal{D} \subset \mathbb{R}^\text{dim}$.
We assume Bayesian transport maps as in Section \ref{sec:regression}, with regressions of the form \eqref{eq:regi} in $(i-1)$-dimensional space for $i=1,\ldots,N$. As $N$ is very large in many relevant applications, we will specify priors distributions of the form described in Section \ref{sec:priors} that induce substantial regularization and sparsity, as a function of hyperparameters $\bm{\theta} = (\theta_{\sigma,1},\theta_{\sigma,2},\theta_{d,1},\theta_{d,2},\theta_\gamma,\theta_q)$.

\subsection{Maximin ordering and nearest neighbors\label{sec:maximin}}

\begin{figure}
\centering
	\begin{subfigure}{.24\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{maxmin13}
	\caption{$i=13$}
	\label{fig:mm1}
	\end{subfigure}%
\hfill
\centering
	\begin{subfigure}{.24\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{maxmin51}
	\caption{$i=51$}
	\label{fig:mm2}
	\end{subfigure}%
\hfill
\centering
	\begin{subfigure}{.24\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{maxmin290}
	\caption{$i=290$}
	\label{fig:mm3}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.24\textwidth}
	\centering
	\includegraphics[width =.99\linewidth]{scales}
	\caption{scale decay}
	\label{fig:scale}
	\end{subfigure}%
  \caption{Maximin ordering (Section \ref{sec:maximin}) for locations on a grid (small gray points) of size $N=60 \times 60 = 3{,}600$ on a unit square, $[0,1]^\text{dim}$ with $\text{dim}=2$. (a)--(c): The $i$th ordered location (${\color{blue}+}$), the previous $i-1$ locations (${\color{black}\circ}$), including the nearest $m=4$ neighbors (${\color{green}\mathbf{x}}$) and the distance $\ell_i$ to the nearest neighbor ({\color{red}---}). (d): For $i=1,\ldots,N$, the length scales (i.e., minimum distances) decay as $\ell_i = i^{-1/\text{dim}}$.}
\label{fig:maxmin}
\end{figure}

A triangular map $\mathcal{T}(\mathbf{y})$ as in \eqref{eq:map} depends on the ordering of the variables $y_1,\ldots,y_N$.
We assume that the entries of $\mathbf{y}$ are ordered according to a maximum-minimum-distance (maximin) ordering of the corresponding locations $\mathbf{s}_1,\ldots,\mathbf{s}_N$ (see Figure \ref{fig:maxmin}), in which we sequentially choose each location to maximize the minimum distance to all previously ordered locations.
Define $c_i(k)$ as the index of the $k$th nearest (previously ordered) neighbor of the $i$th location (and so $\mathbf{s}_{c_i(1)},\ldots,\mathbf{s}_{c_i(4)}$ are indicated by ${\color{blue}+}$ in Figure \ref{fig:maxmin}).

The maximin ordering can be interpreted as a multiresolution decomposition into coarse scales early in the ordering and fine scales later in the ordering. In particular, the minimal pairwise distance $\ell_i = \| \mathbf{s}_i - \mathbf{s}_{c_i(1)} \|$ among the first $i$ locations of the ordering decays roughly as $\ell_i \propto i^{-1/\text{dim}}$, where $\text{dim}$ here is the dimension of the spatial domain (see Figure \ref{fig:scale}).
As a result of the maximin ordering, the $i$th regression in \eqref{eq:regi} can be viewed as a spatial prediction at location $\mathbf{s}_i$ based on data at locations $\mathbf{s}_1,\ldots,\mathbf{s}_{i-1}$ that lie roughly on a regular grid with distance (i.e., scale) $\ell_i$.

When the variables $y_1,\ldots,y_N$ are not associated with spatial locations or when Euclidean distance between the locations is not meaningful (e.g., nonstationary, multivariate, spatio-temporal, or functional data), the maximin and neighbor ordering can be carried out based on other distance metrics, such as $(1 - |\text{correlation}|)^{1/2}$ based on some guess or estimate of the correlation between variables \citep{Kidd2020}.
An example of this idea can be found in Section \ref{sec:application}.

\subsection{Priors on the conditional non-Gaussianity \texorpdfstring{$\sigma_i^2$}{}\label{sec:prior_s}}

In \eqref{eq:kernel}, $\sigma_i^2$ determines the degree of nonlinearity in $f_i$; hence, $\sigma_i^2,\ldots,\sigma_N^2$ together determine the conditional non-Gaussianity in the distribution of $\mathbf{y}_{i:N}$ given $\mathbf{y}_{1:i-1}$.
A priori, we assume that the degree of nonlinearity decays polynomially with length scale $\ell_i$, namely $\sigma_i^2 = e^{\theta_{\sigma,1}} \ell_i^{\theta_{\sigma,2}}$, which allows the conditional distributions of $\mathbf{y}_{i:N}$ given $\mathbf{y}_{1:i-1}$ to be increasingly Gaussian as $i$ increases, as a function of hyperparameters $\theta_{\sigma,1},\theta_{\sigma,2}$.

This prior assumption is motivated by the behavior of stochastic processes with quasiquadratic loglikelihoods.
A quasiquadratic loglikelihood of order $r$ is the sum of a quadratic leading-order term that depends on the $r$-th derivatives of the process, and a nonquadratic term that may only depend on derivatives up to order $r - 1$.
Gaussian smoothness priors (with quadratic loglikelihoods) such as the Mat{\'e}rn model \citep{Whittle1954,whittle1963stochastic} are closely related to linear partial differential equations (PDEs); similarly, stochastic processes with quasiquadratic loglikelihoods are closely related to quasilinear PDEs, for which terms involving the highest-order derivatives are required to be linear.
A wide range of physical phenomena is governed by quasilinear PDEs.
For instance, the Cahn-Hilliard \citep{cahn1958free} and Allen-Cahn \citep{allen1972ground} equations describe phase separation in multi-component systems, the Navier-Stokes equation describes the dynamics of fluids, and the F{\"o}ppl-von K{\'a}rm{\'a}n equations describe the large deformations of thin elastic plates.

By using Poincar{\'e} inequalities \citep{adams2003sobolev} one can show, roughly speaking, that conditional on measurements on a length scale $\ell$, the $r$-th derivatives of a perturbation to the stochastic process have magnitude $\approx \ell^r$.
Thus, the quadratic leading-order term becomes the dominant contribution to the loglikelihood on smaller scales and, accordingly, the conditional distribution is approximately Gaussian.
This is illustrated in a numerical example in Figure \ref{fig:CHconditioning}. A more detailed description of these results and the illustration in Figure \ref{fig:CHconditioning} is provided in Appendix \ref{app:quasilinear}.

\begin{figure}
    \centering
    \includegraphics[width=0.24\textwidth]{plots/cahn_hilliard_heatmap.pdf}
    \includegraphics[width=0.24\textwidth]{plots/cahn_hilliard_uncond_hist.pdf}
    \includegraphics[width=0.24\textwidth]{plots/cahn_hilliard_cond_hist_1.pdf}
    \includegraphics[width=0.24\textwidth]{plots/cahn_hilliard_cond_hist_5.pdf}
    \caption{Samples from the non-Gaussian process in \eqref{eqn:CHdensity} feature regions of negative (dark) and positive (light) values (first panel). The distribution at a given location is a mixture of these two possibilities and thus non-Gaussian (second panel). By contrast, after conditioning on averages over regions of size $\ell = 2^{-1}$ (third panel) or $\ell = 2^{-5}$ (fourth panel), the conditional distribution is close to Gaussian, as these averages determine with high probability whether the location is in a positive or negative region. (See Appendix \ref{app:quasilinear} for details.)}
    \label{fig:CHconditioning}
\end{figure}

\subsection{Priors on the conditional variances \texorpdfstring{$d_i^2$}{}\label{sec:prior_d}}

As we have argued in Section \ref{sec:prior_s}, even non-Gaussian stochastic processes with quasiquadratic loglikelihoods exhibit conditional near-Gaussianity on fine scales. Thus, we will now describe prior assumptions for the $d_i^2$ and $f_i$ in \eqref{eq:regi} that are motivated by the behavior of a transport map $\mathcal{T}$ for a Gaussian target distribution with Mat{\'e}rn covariance (see Figure \ref{fig:matern}), which is a highly popular assumption in spatial statistics.
The Mat{\'e}rn covariance function is also the Green's function of an elliptic PDE \citep[][]{Whittle1954,whittle1963stochastic}.

\citet[][Thm.~2.3]{Schafer2017} show that Gaussian processes with covariance functions given by the Green's function of elliptic PDEs of order $r$ have conditional variance of order $\ell_i^{2r}$ when conditioned on the first $i$ elements of the maximin ordering (see Figure \ref{fig:matern1}).

Hence, for the noise or conditional variances $d_i^2 \sim \mathcal{IG}(\alpha_i,\beta_i)$ as in Section \ref{sec:priors}, we set $E(d_i^2) = \beta_i/(\alpha_i -1) = e^{\theta_{d,1}} \ell_i^{\theta_{d,2}}$. Assuming the prior standard deviation of $d_i^2$ to be equal to $g$ times the mean, we obtain $\alpha_i = 2+1/g^2$ and $\beta_i = e^{\theta_1} \ell_i^{\theta_2} (1+1/g^2)$. For our numerical experiments, we arbitrarily chose $g=4$ to obtain a relatively vague prior for the $d_i^2$.

\begin{figure}
\centering
	\begin{subfigure}{.37\textwidth}
	\centering
  	\includegraphics[trim=4mm 13mm 9mm 22mm, clip,width =.99\linewidth]{matern3d}
	\caption{$\mathbf{y}_i$ vs $\mathbf{y}_{c_i(1)}$, $\mathbf{y}_{c_i(5)}$ for $i=290$}
	\label{fig:maternscatter}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.27\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{matern_nugget}
	\caption{$d_{i}$ ($\circ$) \, and \,  $e^{\theta_{d,1}} \ell_i^{\theta_{d,2}}$ ({\color{red}---})}
	\label{fig:matern1}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.27\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{matern_sq_weights}
	\caption{Avg.\ $\{b_{i,k}^2: i\!=\!1,\ldots,N\}$}
	\label{fig:matern2}
	\end{subfigure
  \caption{For a Gaussian process with exponential covariance on the grid and with the ordering from Figure \ref{fig:maxmin}, expressing the joint distribution $p(\mathbf{y})$ using a transport map as in \eqref{eq:map}--\eqref{eq:condmap} results in a series of regressions as in \eqref{eq:regi} with linear predictors, $f_i(\mathbf{y}_{1:i-1}) = \sum_{k=1}^{i-1} y_{c_i(k)}b_{i,k}$, where $c_i(k)$ indicates the $k$th nearest (previously ordered) neighbor of the $i$th location.
  (For non-Gaussian $p(\mathbf{y})$, the functions $f_i$ are nonlinear.)
  (a): For $n=100$ simulations, the values of $\mathbf{y}_i$ and its 1st and 5th nearest neighbor (NN) lie on a low-dimensional manifold; the regression plane (assuming all other variables to be fixed) indicates a stronger influence of the 1st NN (see the slope of the intersection of the regression plane with the front of the box) than of the 5th NN.
  (b): The conditional standard deviations decay as a function of the length scale $\ell_i$ (see Figure \ref{fig:scale}).
  (c) The squared regression coefficients decay rapidly as a function of neighbor number $k$.
}
\label{fig:matern}
\end{figure}

\subsection{Priors on the regression functions \texorpdfstring{$f_i$}{}\label{sec:prior_f}}

The regression functions $f_i: \mathbb{R}^{i-1} \rightarrow \mathbb{R}$ in \eqref{eq:regi} were specified to be GPs in $(i-1)$-dimensional space in Section \ref{sec:priors}. For the covariance function in \eqref{eq:kernel}, we assume that $\rho_i(\mathbf{y}_{1:i-1},\mathbf{y}_{1:i-1}') = \rho(\mathbf{y}_{1:i-1}^\top\mathbf{Q}_i\mathbf{y}_{1:i-1}'/\gamma)$, where $\gamma = \exp(\theta_\gamma)$ is a range parameter and $\rho$ is an isotropic correlation function, taken to be Mat\'ern with smoothness 1.5 for all numerical experiments here.

To make this potentially high-dimensional regression feasible, we again use the example of a spatial GP with Mat\'ern covariance to motivate regularization and sparsity via the relevance matrix $\mathbf{Q}_i = diag(q_{i,1}^2,\ldots,q_{i,i-1}^2)$.
We assume that the relevance of the $k$th neighbor (see Section \ref{sec:maximin}) decays exponentially as a function of $k$, such that $q_{i,c_i(k)}$ decays as $\exp(\theta_q k)$. This type of behavior, often referred to as the screening effect \citep[e.g.,][]{stein20112010}, is illustrated in Figure \ref{fig:matern2}, and it has been exploited for covariance estimation of a Gaussian spatial field by \citet{Kidd2020}.
Recently, \citet{Schafer2017} proved exponential rates of screening for Gaussian processes derived from elliptic boundary-value problems; following the discussion in Section~\ref{sec:prior_s}, we expect similar conditional-independence phenomena to hold on the fine scales of processes with quasiquadratic loglikelihoods.

Given this exponential decay as a function of the neighbor number $k$, the relevance will be essentially zero for sufficiently large $k$, and so we achieve sparsity by setting
\begin{equation}
q_{i,c_i(k)} = \begin{cases} \exp(\theta_q k), & k \leq m,\\ 0, & k>m, \end{cases}
\label{eq:sparsity}
\end{equation}
where the sparsity parameter $m=\max \{k: \exp(\theta_q k) \geq \varepsilon\}$ is determined by the data through the hyperparameter $\theta_q$. We use $\varepsilon = 0.01$ for our numerical examples.
This assumption induces a sparse transport map, in that $f_i$ (and thus $\mathcal{T}_i$) depend on $\mathbf{y}_{1:i-1}$ only through the $m$ nearest neighbors $y_{c_i(1)},\ldots,y_{c_i(m)}$, where $\rho_i$ is isotropic as a function of the scaled inputs $y_{c_i(k)}/q_{i,c_i(k)}$.
Sparsity in the transport map is equivalent to an assumption of ordered conditional independence \citep[cf.][]{Spantini2018}, which can be seen from the expression of the conditional distributions in \eqref{eq:map2norm}.
Similar ordered-conditional-independence assumptions are also popular for fast approximations of Gaussian fields with parametric covariance functions \citep[e.g.,][]{Vecchia1988, Stein2004,Datta2016,Guinness2016a,Katzfuss2017a,Katzfuss2018,Katzfuss2020,Schafer2020}.

Identifying the regression functions $f_i$ in $m$-dimensional space is further aided by the data concentrating on a lower-dimensional manifold due to the strong dependence between most $y_{c_i(k)}$ and $y_{c_i(l)}$ for small $k,l \leq m$ (e.g., see Figure \ref{fig:maternscatter}).

\subsection{Inference\label{sec:spatialinference}}

Based on the prior distributions in Sections \ref{sec:prior_s}--\ref{sec:prior_f}, we can carry out inference and compute the transport map as in Section \ref{sec:invertiblemap}. Due to the sparsity assumption in \eqref{eq:sparsity}, the computational complexity is lower than in Section \ref{sec:invertiblemap}; specifically, determining $\widetilde{\mathcal{T}}_i$ now only requires $\mathcal{O}(n^3 + m n^2)$ time, again in parallel for each $i=1,\ldots,N$. Each application of the transport map or its inverse then requires $\mathcal{O}(N(n^2 + mn))$ time.
The maximin ordering and nearest neighbors can also be computed in quasilinear time in $N$ \citep[][Alg.~7]{Schafer2020}.

The prior distributions in Sections \ref{sec:prior_s}--\ref{sec:prior_f} depend on a vector of hyperparameters, $\bm{\theta} = (\theta_{\sigma,1},\theta_{\sigma,2},\theta_{d,1},\theta_{d,2},\theta_\gamma,\theta_q)$. When making inference on $\bm{\theta}$ as described in Section \ref{sec:hyper}, we effectively let the training data $\mathbf{Y}$ decide the degree of sparsity (through $\theta_q$ via $m$) and the degree of nonlinearity (through $\theta_{\sigma,1},\theta_{\sigma,2}$ via $\sigma_i$).

Figure \ref{fig:sine} illustrates estimation of transport-map components in a simulated example.

\begin{figure}
\centering
	\begin{subfigure}{.22\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_bw_true}
	\caption{True}
	\label{fig:sinetrue}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.22\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_bw_data}
	\caption{Data}
	\label{fig:sinedata}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.22\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_bw_lin}
	\caption{Linear}
	\label{fig:sinelin}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.22\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_bw_nonlin}
	\caption{Nonlinear}
	\label{fig:sinenonlin}
	\end{subfigure}%
\hfill
\begin{minipage}{.042\textwidth}
~

\vspace{-10mm}

\includegraphics[trim=78mm 0mm 5mm 0mm, clip,width =.93\linewidth]{sine_bw_legend}
\end{minipage}
  \caption{Simulation from a nonlinear map with sine structure in $f_i$, described as NR900 in Section \ref{sec:simstudy}. For $n=100$ and $i=80$, $y_i$ versus its 1st and 2nd nearest neighbor (NN): true $f_i$ (a), observations $\mathbf{y}_i$ (b), together with linear (c) and nonlinear (d) fit (i.e., posterior means) of $f_i$, with further variables in $\mathbf{y}_{1:i-1}$ held at their mean levels. The linear map in (c) is estimated under the restriction $\sigma_i=0$. In (d), we have a nonlinear regression in $79$-dimensional space, with $m=5$ active variables in the estimated (via $\bm{\theta}$) nonlinear model.}
\label{fig:sine}
\end{figure}

In Section \ref{sec:invertiblemap}, we discussed using $\widetilde{\mathcal{T}}$ in \eqref{eq:singlemap} to transform the non-Gaussian $\mathbf{y}^\star$ to standard Gaussian map coefficients $\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$. This concept, which is illustrated in Figure \ref{fig:trans}, is especially interesting in our spatial setting. Due to the maximin ordering (Figure \ref{fig:maxmin}), the scales $\ell_i$ are arranged in decreasing order, and in our prior the $d_i^2$ also follow a decreasing stochastic order with $E(d_i^2) = e^{\theta_{d,1}} \ell_i^{\theta_{d,2}}$ (see, e.g., Figure \ref{fig:matern1}). Thus, we can view the map components as a form of nonlinear principal components (NPCs), with the map coefficients as the corresponding component scores.
For Gaussian processes with covariance functions given by the Green's function of elliptic PDEs, similar to the Mat{\'e}rn family, it can be shown that these principal components based on the maximin ordering are approximately optimal \citep{Schafer2017}.
For example, these NPCs can be used for dimension reduction by only storing or modeling the first $k$, say, map coefficients $\mathbf{z}^\star_{1:k} = (z^\star_1,\ldots,z^\star_k)^\top$.
Note that if we set $\mathbf{z}_{k+1:N}^\star = \mathbf{0}$, we assume $y_i^\star = \hat f_i(\mathbf{y}_{1:i-1}^\star)$ for $i>k$, which overestimates dependence and underestimates variability; hence, it is preferable to draw $\mathbf{z}_{k+1:N}^* \sim \mathcal{N}(\mathbf{0},\mathbf{I})$.
In addition to reducing storage, we can also use this approach for conditional simulation \citep[][Lemma 1]{Marzouk2016}, in which we fix the large-scale features of an observed field by fixing the first $k$ map coefficients (see Figure \ref{fig:precipcond} for an illustration).
To model a time series of spatial fields, we can for example assume a linear vector autoregressive model for the NPCs, such that the map coefficients at time $t+1$, say $\mathbf{z}_{1:k}^{(t+1)}$, linearly depend on $\mathbf{z}_{1:k}^{(t)}$.
When it is of interest to regress some response on a spatial field, one could also use the first $k$ map coefficients of the field as the covariates, similar to the use of function principal component scores in regression.

\begin{figure}
\centering
	\begin{subfigure}{.19\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{matern_trans_data}
	\caption{Transformed}
	\label{fig:transdata}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.19\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{matern_trans_qq}
	\caption{QQ plot}
	\label{fig:qq}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.19\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_avg_test}
	\caption{Test data}
	\label{fig:avgtest}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.19\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_avg_latent}
	\caption{Reference avg}
	\label{fig:avglatent}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.19\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{sine_avg_orig}
	\caption{Original avg}
	\label{fig:avgorig}
	\end{subfigure
\hfill
\begin{minipage}{.036\textwidth}
~

\vspace{-9mm}

\includegraphics[trim=78mm 0mm 5mm 0mm, clip,width =.93\linewidth]{sine_avg_legend}
\end{minipage}
  \caption{Illustration of map coefficients $\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$ (see Sections \ref{sec:invertiblemap} and \ref{sec:spatialinference}) for the simulated NR900 data using $\widetilde{\mathcal{T}}$ inferred from $n=100$ training data.
  (a)--(b): The $N=900$ map coefficients corresponding to one test sample are roughly iid Gaussian.
  (c): For $1{,}000$ test samples, $y^*_{80}$ versus 1st and 2nd NNs (cf.~Figure \ref{fig:sine}).
  (d) When averaging pairs of two map-coefficient vectors in reference space and transforming back to the original space using \eqref{eq:invmap}, the sinusoidal relationship between $y^*_{80}$ and its NNs is preserved in the resulting 500 averages.
  (e) When averaging test samples directly in original space, the nonlinear structure is lost.}
\label{fig:trans}
\end{figure}

\section{Marginal non-Gaussianity\label{sec:dpm}}

So far, we have focused on nonlinear, non-Gaussian dependence structures; the model described in Sections \ref{sec:bayesian} and \ref{sec:spatial} is most appropriate if the marginal distribution of each $y_i$ is Gaussian. If this does not hold at least approximately, extensions based on additional marginal (i.e., pointwise) transformations are straightforward. For example, assume that the model from Sections \ref{sec:bayesian} and \ref{sec:spatial} holds for $\mathbf{y}$, but that we actually observe $\tilde\mathbf{y} = \mathcal{G}(\mathbf{y})$ such that $\tilde y_i = g_i(y_i)$. If the $g_i$ are one-to-one differentiable functions, the resulting posterior map is a simple extension of that in Proposition \ref{prop:maps}. The $g_i$ can be pre-determined (see Section \ref{sec:application} for an example with a log transform) or may depend on $\bm{\theta}$ and thus be inferred based on a minor modification of the integrated likelihood in Proposition \ref{prop:lik}.

For marginal distributions that are continuous but are not easily transformed to Gaussianity, we will extend our GP regressions with Gaussian error (see \eqref{eq:regi}) to modeling the GP errors $\epsilon_i^{(j)}$ using Bayesian nonparametrics for all $i=1,\ldots,N$. More precisely, we will use Dirichlet process mixtures (DPM) with Gaussian base measures. In \eqref{eq:regi}, we now assume that $f_i(\cdot) \sim \mathcal{GP}(0,C_i)$, and the $\epsilon_i^{(j)}$ are distributed according to a DPM:
\begin{equation}
    \epsilon_i^{(j)} | \mu_i,d_i \sim \mathcal{N}(\mu_i,d_i^2), \qquad (\mu_i,d_i^2) | \mathcal{F}_i \sim \mathcal{F}_i, \qquad \mathcal{F}_i \sim \mathcal{DP}(\mathcal{NIG}(\xi_i,\eta_i,\alpha_i,\beta_i),\zeta_i),
\end{equation}
where $\zeta_i$ is the concentration parameter, and the base measure $\mathcal{NIG}$ is a normal-inverse-Gamma distribution with parameters $\xi_i,\eta_i,\alpha_i,\beta_i$ with $\xi_i=0$.
The degree of non-Gaussianity allowed for the $\epsilon_i^{(j)}$ is determined by $\eta_i$ and $\zeta_i$; a small value of $\zeta_i$ concentrates the Dirichlet process near the NIG base measure, for which a large value of $\eta_i$ shrinks $\mu_i$ toward zero. Thus, in the limit as $\zeta_i \rightarrow 0$ and $\eta_i \rightarrow \infty$, we obtain a model similar to that in Section \ref{sec:priors}
(except that here $d_i^2$ does not appear in the variance of the GP $f_i$).

For the spatial setting with maximin ordering of Section \ref{sec:spatial}, we can again find a sparse parameterization in terms of hyperparameters $\bm{\theta} = (\theta_{\sigma,1},\theta_{\sigma,2},\theta_{d,1},\theta_{d,2},\theta_\gamma,\theta_q,\theta_{\zeta,1},\theta_{\zeta,2},\theta_{\eta,1},\theta_{\eta,2})$. We parameterize the $\alpha_i$, $\beta_i$, $C_i$ in terms of the first six hyperparameters as in Sections \ref{sec:prior_s}--\ref{sec:prior_f}. For the concentration parameter $\zeta_i = e^{\theta_{\zeta,1}}\ell_i^{\theta_{\zeta,2}}$, we allow increasing shrinkage toward Gaussianity for increasing $i$. We similarly set $\eta_i = e^{\theta_{\eta,1}}\ell_i^{\theta_{\eta,2}}$.
For this DPM model, we take a fully Bayesian perspective and assume an improper uniform prior for $\bm{\theta}$ over $\mathbb{R}^{10}$.

Due to the nonparametric nature of the DPM, the universal approximation property of GPs \citep{Micchelli2006}, and nonzero prior probability for the dense (non-sparse) transport map, the posterior distribution obtained using this model contracts (for $n\rightarrow \infty$ and fixed $N$) to the Kullback-Leibler (KL) projection of the actual distribution of $\mathbf{y}$ onto the space of distributions that can be described by a transport map whose components are additive in the $i$th argument as in \eqref{eq:condmap}, due to the KL optimality of the Knothe-Rosenblatt map \citep[][Sec.~4.1]{Marzouk2016}. In other words, as the number of replicates increases, the learned distribution gets as close as possible to the truth under the additivity restriction.

Inference for our DPM model cannot be carried out in closed form anymore and instead relies on a Metropolis-within-Gibbs Markov chain Monte Carlo (MCMC) sampler. We can also compute and draw samples from the posterior predictive distribution
\[
\textstyle
p(\mathbf{y}^\star|\mathbf{Y}) = \prod_{i=1}^N p(y_i^\star|\mathbf{y}_{1:i-1}^\star,\mathbf{Y}),
\]
for which each $p(y_i^\star|\mathbf{y}_{1:i-1}^\star,\mathbf{Y})$ is approximated as a Gaussian mixture based on the MCMC output.
Details for the MCMC procedure and the posterior predictive distribution are given in Appendix \ref{app:dpmgibbs}.

In the spatial setting with sparsity parameter $m$, each MCMC iteration still has time complexity $\mathcal{O}(N(n^3 + n^2m))$ and the computations within each iteration are highly parallel; however, the actual computational cost for this sampler is much higher than in Section \ref{sec:spatialinference} due to the typically large number of MCMC iterations required.

\section{Simulation study \label{sec:simstudy}}

We carried out a simulation study with data simulated based on transport maps as in Section \ref{sec:regression}, with the $f_i$ and $d_i$ as implied by various extensions of a Gaussian process (GP) with exponential covariance evaluated at locations on a unit square. Specifically, we considered four simulation scenarios, for which samples are illustrated in the top row of Figure \ref{fig:maternsine}:
\begin{description}[itemsep=1pt,topsep=2pt,parsep=1pt]
    \item[LR900:] \textbf{L}inear map (i.e., a Gaussian distribution) with components $f_i^\text{L}(\mathbf{y}_{1:i-1}) = \sum_{k=1}^{i-1} b_{i,k} y_{c_i(k)}$, where the $b_{i,k}$ are based on an exponential covariance with unit variance and range parameter 0.3 evaluated on a \textbf{R}egular grid of size $N = 30 \times 30 = \mathbf{900}$.
    \item[NR900:] \textbf{N}onlinear extension of LR900 by a sine function of a weighted sum of the nearest two neighbors: $f_i^\text{NL}(\mathbf{y}_{1:i-1}) = f_i^\text{L}(\mathbf{y}_{1:i-1}) + 2 \sin(4(b_{i,1} y_{c_i(1)}+b_{i,2} y_{c_i(2)}))$ (see Figure \ref{fig:sinetrue})
    \item[NI3600:] Same as NR900, but at $N=\mathbf{3{,}600}$ \textbf{I}rregularly spaced locations sampled uniformly at random
    \item[NR900B:] Same as NR900, but with a \textbf{B}imodal distribution for the $\epsilon_i$ in \eqref{eq:regi}: $\epsilon_i | \mu_i \sim \mathcal{N}(\mu_i,d_i^2)$ with $\mu_i$ sampled from $\{-3.5 d_i, 3.5 d_i\}$ with equal probability
\end{description}
For computational simplicity, each (true) $f_i$ was assumed to only depend on the nearest 30 previously ordered neighbors, but this gives a highly accurate approximation of a ``full'' exponential covariance in the linear MR900 case.

\begin{figure}
\centering
	\begin{subfigure}{.244\textwidth}
	\centering
	\includegraphics[trim=-22mm 0mm 7mm 0mm, clip,width =.98\linewidth]{simLR900} \\
	\includegraphics[width = 0.98\linewidth]{klcomp_sin_lin}
	\caption{LR900}
	\label{fig:ms1}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.244\textwidth}
	\centering
	\includegraphics[trim=-22mm 0mm 7mm 0mm, clip,width =.98\linewidth]{simNR900} \\
 	\includegraphics[width = 0.98\linewidth]{klcomp_sin}
 	\caption{NR900}
	\label{fig:ms2}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.244\textwidth}
	\centering
	 \includegraphics[trim=-22mm 0mm 7mm 0mm, clip,width =.98\linewidth]{simNI3600} \\
 	\includegraphics[width = 0.98\linewidth]{klcomp_sin_random}
 	\caption{NI3600}
	\label{fig:ms3}
	\end{subfigure}
	\hfill
	\begin{subfigure}{.244\textwidth}
	\centering
	\includegraphics[trim=-22mm 0mm 7mm 0mm, clip,width =.98\linewidth]{simNR900B} \\
 	\includegraphics[width = 0.98\linewidth]{klcomp_sin_DPM}
 	\caption{NR900B}
	\label{fig:ms4}
	\end{subfigure}
  \caption{Top row: Simulated spatial fields for four simulation scenarios described in Section \ref{sec:simstudy}. Bottom row: Corresponding comparisons of KL divergence as a function of ensemble size $n$ (on a log scale) for different methods. In (b)--(c), the KL divergence for \texttt{S-nonlin} was extremely high and thus is not shown. \texttt{DPM} is only included in (d).}
\label{fig:maternsine}
\end{figure}

We compared the following methods:
\begin{description}[itemsep=1pt,topsep=2pt,parsep=1pt]
\item[\texttt{nonlin}:] Our method with Bayesian uncertainty quantification described in Section \ref{sec:spatial}
\item[\texttt{S-nonlin}:] Simplified version of \texttt{nonlin} ignoring uncertainty in the $f_i$ and $d_i$, as in \eqref{eq:detmap}
\item[\texttt{linear}:] Same as \texttt{nonlin}, but forcing $\theta_{\sigma,1} = -\infty$ and hence linear $f_i$
\item[\texttt{S-linear}:] Simplified version of \texttt{linear} ignoring uncertainty in the $f_i$ and $d_i$, as in \eqref{eq:detmap}
\item[\texttt{DPM}:] The model with Dirichlet process mixture residuals described in Section \ref{sec:dpm}
\end{description}
We carried out fully Bayesian inference for \texttt{DPM} using the MCMC procedure in Appendix \ref{app:dpmgibbs}; for all other methods, we carried out empirical Bayesian inference based on a plug-in estimate of the hyperparameters $\bm{\theta}$ as described at the end of Section \ref{sec:hyper}. \texttt{S-linear} is the only method that results in a joint Gaussian posterior predictive distribution; it is similar to the approach proposed and used in numerical comparisons in \citet{Kidd2020}.

Some properties of our nonlinear approach are illustrated in Figures \ref{fig:sine} and \ref{fig:trans}. Our nonlinear method estimated $f_i$ well in the observed high-probability region.

Figure \ref{fig:maternsine} shows the results of more systematic comparisons using the Kullback-Leibler (KL) divergence,
\begin{equation}
\label{eq:kl}
E( \log p_0(\mathbf{y}) )  - E(\log p(\mathbf{y}|\mathbf{Y}) )
\end{equation}
between the true distribution $p_0(\mathbf{y})$ and the distribution $p(\mathbf{y}|\mathbf{Y})$ implied by the inferred posterior map (see \eqref{eq:tstar}), where the expectations are taken with respect to the true distribution. We approximated the expectations by the sample means over 50 simulated test fields $\mathbf{y}^\star$, and so the resulting KL divergence is the difference of the log-scores \citep[e.g.][]{Gneiting2014} of the true and inferred distributions.
Whenever nonlinear structure was not discernible from the data (because the true map was linear or because the ensemble size $n$ was too small), \texttt{nonlin} performed similarly to \texttt{linear} and hence did not suffer due to its over-flexibility. For larger ensemble size and nonlinear truths, \texttt{nonlin} at times far outperformed \texttt{linear}. \texttt{S-linear} and \texttt{S-nonlin} were generally less accurate than their counterparts with uncertainty quantification; in the linear LR900 setting, this was only an issue for small ensemble size, but \texttt{S-nonlin} performed extremely poorly when the nonlinear structure was clearly apparent in the data, likely due to overfitting without accounting for uncertainty in the fit.

\section{Climate-data application \label{sec:application}}

An important application of our methods is the analysis and emulation of output from climate models. Climate models are essentially large sets of computer code describing the behavior of the Earth system (e.g., the atmosphere) via systems of differential equations. Much time and resources have been spent on developing these models, and enormous computational power is required to produce ensembles (i.e., solve the differential equations for different starting conditions) on fine latitude-longitude grids for various scenarios of greenhouse-gas emissions. Of the large amount of data and output that have been generated, only a small fraction has been fully explored or analyzed \citep[e.g.][]{Benestad2017}.
A popular way of analyzing climate-model output is statistical emulation \citep[e.g.,][]{Castruccio2013,Castruccio2014,Nychka2018}, which aims to infer the distribution of one or more climate variables, so that relevant summaries or additional samples can be computed more cheaply than via more runs of the full climate model.

We considered log-transformed total precipitation rate (in m/s) on a roughly $1^\circ$ longitude-latitude global grid of size $N = 288 \times 192 = 55{,}296$ on July 1 in $98$ consecutive years, starting in the year 402, from the Community Earth System Model (CESM) Large Ensemble Project \citep{Kay2015}.
We obtained precipitation anomalies by standardizing the data at each grid location to mean zero and variance one, shown in Figure \ref{fig:precdata}.

For ease of comparison and illustration, we first considered a smaller of size $N = 37 \times 74 = 2{,}738$ in a subregion containing large parts of the Americas ($45^\circ$S to $45^\circ$N and $130^\circ$W to $30^\circ$W). The anomalies in this subregion exhibited a nonstationary dependence structure due to ocean, land, and mountains. Hence, we used the correlation-based ordering discussed in Section \ref{sec:maximin}, based on a correlation-matrix estimate given by the element-wise product between the sample covariance and an exponential covariance with a large range parameter.
As shown in Figure \ref{fig:pre3}, the precipitation anomalies exhibited similar features as our simulated data in Figure \ref{fig:maternscatter}, with regression data concentrating on lower-dimensional manifolds.

\begin{figure}
\centering
~
\hfill
	\begin{subfigure}{.36\textwidth}
	\centering
  	\includegraphics[trim=4mm 13mm 9mm 22mm, clip,width =.99\linewidth]{precip_data_2000_3D}
	\caption{$\mathbf{y}_{2000}$ vs 1st and 2nd NN}
	\label{fig:pre3}
	\end{subfigure
\hfill
	\begin{subfigure}{.264\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{ls_precip}
	\caption{LS (Americas)}
	\label{fig:pre4}
	\end{subfigure}%
\hfill
	\begin{subfigure}{.264\textwidth}
	\centering
 	\includegraphics[width =.99\linewidth]{ls_precip_global}
	\caption{LS (global)}
	\label{fig:global2}
	\end{subfigure}%
\hfill
~
	\caption{
	For precipitation anomalies:
	(a) Our approach can be viewed as $N$ regressions, of each $y_i$ on ordered nearest neighbors (NNs), as in Figure \ref{fig:matern}, with the regression data on low-dimensional manifolds.
(b) Comparison of log-score (LS; equal to KL divergence up to an additive constant) for estimated joint distribution for the Americas subregion as a function of ensemble size $n$; the LS for \texttt{S-nonlin}, \texttt{tapSamp}, and \texttt{expCov} were too high (all $>20{,}000$) and are thus not visible.
(c) LS for \texttt{linear} and \texttt{nonlin} for precipitation anomalies on the global grid.}
\label{fig:precip}
\end{figure}

For comparing the methods from Section \ref{sec:simstudy} on the precipitation anomalies, computing the KL divergence as in \eqref{eq:kl} was not possible, as the true distribution $p_0(\mathbf{y})$ was unknown. Hence, we compared the methods using various training data sizes $n$ in terms of log-scores, which are approximately equal to the KL divergence up to an additive constant; specifically, these log-scores consist of the second part of \eqref{eq:kl}, $- E(\log p(\mathbf{y}|\mathbf{Y}) )$, approximated by averaging over 18 test replicates and over five random training/test splits. The comparison for the Americas subregion is shown in Figure \ref{fig:pre4}; \texttt{nonlin} outperformed \texttt{linear}, and \texttt{DPM} was even more accurate than \texttt{nonlin} for large $n$, indicating that the precipitation anomalies exhibited joint and marginal non-Gaussian features.
As in Section \ref{sec:simstudy}, \texttt{S-linear} and \texttt{S-nonlin} performed poorly due to ignoring uncertainty in the estimated map.
We also compared to two baseline competitors that combine a (zero-mean) multivariate-Gaussian assumption with different covariance-matrix estimates: the sample covariance tapered (i.e., element-wise multiplied) by an exponential correlation matrix with range equal to the maximum pairwise distance among the locations (\texttt{sampTap}), and an exponential covariance with variance and range estimated from the training data (\texttt{expCov}); both baseline competitors performed poorly.

We also considered the map coefficients $\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$ discussed in Sections \ref{sec:invertiblemap} and \ref{sec:spatialinference}, using the map obtained by fitting \texttt{nonlinear} to the first $n=97$ replicates as training data.
In Figure \ref{fig:precip_coef}, the map coefficients for a held-out test field appeared roughly iid standard Gaussian, with sample autocorrelations near zero (not shown).
Figure \ref{fig:precipcoef} also illustrates that the map coefficients offer similar properties for non-Gaussian fields as principal-component scores do for Gaussian settings. For example, the medians of the posterior distributions of the $d_i$ (see \eqref{eq:nigpost}) decreased rapidly as a function of $i$, which means that the map coefficients early in the maximin ordering captured much more (nonlinear) variation than later-ordered coefficients (see, e.g., \eqref{eq:invmap} and \eqref{eq:detmap}). Further, we computed the map coefficients for all 98 replicates for July 2--30 (still based on the posterior map trained on July 1 data), and the lag-1 autocorrelation over time between map coefficients also decreased with $i$. Specifically, while most of the first 100 were greater than 0.2, many later autocorrelations were negligible; this indicates that a spatio-temporal analysis could proceed by fitting a simple (linear) autoregressive model over time to only the first $k$, say, map coefficients, while treating the remaining $N-k$ coefficients as independent over time.

\begin{figure}
\centering
~
\hfill
	\begin{subfigure}{.26\textwidth}
	\centering
  	\includegraphics[width =.95\linewidth]{precip_coef}
	\caption{Test map coefficients $\mathbf{z}^\star$}
	\label{fig:precip_coef}
	\end{subfigure
\hfill
	\begin{subfigure}{.26\textwidth}
	\centering
  	\includegraphics[width =.95\linewidth]{precip_di}
	\caption{Posterior median of $d_i$}
	\label{fig:precip_di}
	\end{subfigure
\hfill
	\begin{subfigure}{.26\textwidth}
	\centering
 	\includegraphics[width =.95\linewidth]{precip_lag1autocor}
	\caption{Lag-1 autocorrelation}
	\label{fig:precip_acf}
	\end{subfigure
\hfill
~
	\caption{Properties of the map coefficients $\mathbf{z}^\star = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$ for the precipitation anomalies on the grid of size $N=2{,}738$ in the Americas subregion.
	(a): The map coefficients corresponding to the test field in Figure \ref{fig:prec1} in the original data ordering (first by longitude, then latitude) appeared roughly iid standard Gaussian, aside from slightly heavy tails.
	(b): The posterior medians of the $d_i$ decreased rapidly as a function of $i$, meaning that the first few map coefficients captured much more variation than later-ordered coefficients.
	(c): The autocorrelation between consecutive days also decreased with $i$; while most were greater than 0.2 for $i<100$, many autocorrelations for $i>100$ were negligible.}
\label{fig:precipcoef}
\end{figure}

\begin{figure}
\centering
\begin{minipage}{.94\textwidth}
\centering
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_cond_55296.png}
 	\caption{Test field ($i=N$)}
 	\label{fig:prec1}
	\end{subfigure}%
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_cond_5000.png}
 	\caption{$i=5{,}000$}
	\end{subfigure}%
	\\
	\smallskip
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_cond_500.png}
 	\caption{$i=500$}
	\end{subfigure}%
	\begin{subfigure}{.495\linewidth}
	\centering
 	\includegraphics[width =.99\linewidth]{precip_global_cond_0.png}
 	\caption{Unconditional ($i=0$)}
 	\label{fig:prec4}
	\end{subfigure}%
\end{minipage}
\hfill
\begin{minipage}{.05\textwidth}
\includegraphics[trim=400mm 0mm 0mm 0mm, clip,width =.98\linewidth]{precip_global_legend.png}
\end{minipage}
  \caption{For the global climate data ($N=55{,}296$), we fitted a stochastic emulator using \texttt{nonlin} based on $n=97$ training replicates. Given a held-out test field $\mathbf{y}^*$ in (a), we show conditional simulations based on fixing the first $i$ map coefficients in $\mathbf{z}^* = \widetilde{\mathcal{T}}(\mathbf{y}^\star)$. (b): Only differs in some fine-scale features from (a). (c): Some large-scale features from (a) are preserved. (d): Unconditional simulation (i.e., independent from (a)).}
\label{fig:precipcond}
\end{figure}

To demonstrate scalability to large datasets, we compared \texttt{linear} and \texttt{nonlinear} on the entire global precipitation anomaly fields of size $N = 288 \times 192 = 55{,}296$. We used a maximin ordering based on chordal distance. As shown in Figure \ref{fig:global2}, \texttt{nonlin} outperformed \texttt{linear} even more decisively than for the Americas subregion. Even in the largest and most accurate setting ($n=80$), the estimated $\bm{\theta}$ for \texttt{nonlinear} implied $m=9$, meaning that the corresponding transport map was extremely sparse and hence computationally efficient. Also, all but 113 of the $N=55{,}296$ posterior medians of the $d_i$ were more than 20 times smaller than the largest posterior median (i.e., that of $d_1$), indicating that our approach could be used for massive dimension reduction without losing too much information.

Finally, the fitted map (or rather, its inverse $\widetilde{\mathcal{T}}^{-1}$) can also be viewed as a stochastic emulator of the climate model. Specifically, we can produce a new precipitation-anomaly sample by drawing $\mathbf{z}^* \sim \mathcal{N}_N(\mathbf{0},\mathbf{I}_N)$ and then computing $\mathbf{y}^* = \widetilde{\mathcal{T}}^{-1}(\mathbf{z}^*)$. One such sample (for the full global grid) is shown in Figure \ref{fig:prec4}; it can be generated in a few seconds on a basic desktop computer. Further, this approach can also be used to draw conditional samples, in which we fix the first $i$, say, map coefficients, for example at the values corresponding to a given spatial field. Such draws, which maintain the large-scale features in the held-out (98th) test field but allow for newly sampled fine-scale features, are shown in Figure \ref{fig:precipcond}.
This is related to the supervised conditional sampling ideas in \citet{Kovachki2020}, with their inputs given by our first $i$ ordered test observations.

\section{Conclusions\label{sec:conclusions}}

We have developed a Bayesian approach to inferring a non-Gaussian target distribution via a transport map from the target to a standard normal distribution. The components of the map are modeled using Gaussian processes. For the distribution of spatial fields, we have developed specific prior assumptions that result in sparse maps and thus scalability to high dimensions.
Instead of manually or iteratively expanding a finite-dimensional parameterization of the transport map \citep[e.g.,][]{ElMoselhy2012,Bigoni2016,Marzouk2016}, our Bayesian approach probabilistically regularizes the map; the resulting approach is flexible and nonparametric, but guards against overfitting and quantifies uncertainty in the estimation of the map. In our numerical comparisons, accounting for uncertainty led to vastly improved accuracy, and the flexibility did not reduce the accuracy even when the true distribution was Gaussian. Due to conjugate priors and the resulting closed-form expressions for the posterior map and its inverse, our approach also allows us to convert non-Gaussian data into iid normal map coefficients, which can be thought of as a nonlinear extension of principal components.

As our approach essentially turns estimation of a high-dimensional joint distribution into a series of GP regressions, it is straightforward to include additional covariates and examine their nonlinear, non-Gaussian effect on the distribution.
Shrinkage toward a joint Gaussian distribution with a parametric covariance function could be achieved by assuming the mean for the GP regressions to be the one implied by a Vecchia approximation of that covariance function \citep{Kidd2020}.
Extensions to more complicated input domains (e.g., space-time) could be obtained using correlation-based ordering (Section \ref{sec:maximin}). More generally, one could envision an iterative procedure: find an initial estimate of the dependence structure, carry out our transport-map procedure based on the resulting correlation-based ordering, and then use the dependence implied by the resulting map to reorder, and so on.
Finally, another major avenue of future work would be to use the inferred distribution as the prior of a latent field, which we then update to obtain a posterior given noisy observations of the field; among numerous other applications, this would enable the use of our technique to infer the forecast distribution and account for uncertainty in ensemble-based data assimilation \citep{Boyles2020}, leading to nonlinear updates for non-Gaussian applications.

\footnotesize


text:
\section{Introduction\label{Sec:Intro}}

The spin glass corresponding to the random variant of the Potts model
\cite{Wu} occured to be one of the first models where the scenario of the
transition to the nonergodic state was discovered \cite{ElSh83,GoEl85} to
be quite different from the standard one, that is of the
Sherrington--Kirkpatrick (SK) model \cite{sk}. Nevertheless, it remains to
be in the focus of investigations till now and a number of questions
still are not answered (see, e.g., the recent papers
Ref.\cite{Janis, Al} that are also dealing with the mean field Potts
glasses with few number of states).

It is well known that in the case of SK model the so callled full replica
symmetry breaking (FRSB)
\cite{Par} takes place at the instability point $T_c$ of RS solution. All
other models with two-particle interaction and with reflection symmetry
behave in the same way \cite{Full}.  The Potts glass models (with $p>2$)
are usually considered to be classical examples of models without
 reflection symmetry: FRSB does not work at $T_c$
\cite{ElSh83,GoEl85,Gross,NoSh93,Ritort,Gribova} and there is a region
where the 1RSB solution is stable (see, however,
Ref.\cite{Janis}). It is known (see, e.g., Ref.\cite{ Ritort}
and references therein) from the numerical solution, that if the number of
states $p\leq4$, then the 1RSB solution appears continuously at $T_c$ --
the bifurcation point of the RS solution. For $p>4$ the 1RSB solution
appears jumpwise. The Potts model with $p=4$ presents a boundary between
two types of phase transitions. In this case the order parameter $m$
defining the number of groups in 1RSB construction achieves its maximum
value $m=1$ just at the point $T_c=T_{1RSB}$. It is worthwhile to follow
in an analytic way, using the bifurcation theory, how the crossover from
continuous to jumpwise behavior of the order parameter happens. This paper
is devoted to this problem. To some extent we proceed in a same line as
the authors of the paper Ref.\cite{A.Crisantione} where the SK model
was considered.

It is worth to notice that the Potts symmetry leads to the cancellation of
some terms in the free energy. This by-turn
results in a trivial solution for the order parameters in  the RS
approximation. This is why the Potts models differ from other models
without reflection symmetry, in particular, from the models similar to SK
model, that is with the Hamiltonian
\begin{equation}
H=-\frac{1}{2}\sum_{i\neq j}J_{ij} U_i U_j, \label{SK}
\end{equation}
with quenched interactions $J_{ij}$  distributed with Gaussian
probability,
but without reflection symmetry:
\begin{equation}
Tr\left[{U}^{(2k+1)}\right]\neq 0 \label{odd}
\end{equation}
for some integer $k$. In this case the disorder smears out the
 phase transition: there is no trivial RS solution at high temperature
(see, e.g. \cite{Luchinskaya, Hydrogen}).

At the same time the Potts symmetry leads to the fact that the order
parameter is one and the same for all of $p$ states, that makes it related
to two--particle generalized SK--type models and allows to show that the
dependence of the free energy on the order parameters in the vicinity of
the bifurcation point $T_c$ is universal in both cases.

 We use a special representation of Potts models which originates from
 the analogy between three--state Potts spin glass
 and isotropic quadrupole glass (see below). In our representation the
 number of integrations is less than in usual approach and so it is more
 convenient for actual calculations. Moreover, this representation
 demonstrates explicitly the appearence of cubic terms in the free
 energy expansion. These terms play a crucial role for the
 order-parameter behavior.

The Potts glass model is defined by the Hamiltonian
\begin{equation}
\hat H=-\frac{p}{2}\sum_{i\neq j}J_{ij} {\delta_{\sigma_{i},\sigma_{j}}}.
\label{one1}
\end{equation}
where $p$ is the number of states and
variables $\sigma_i$, $\sigma_j$ can take the values $0,1,...,p-1$.
Here  $J_{ij}$ are random interactions
distributed with Gaussian probability
\begin{equation} P(J_{ij})=\frac{\sqrt{N}}{\sqrt{2\pi}
\tilde{J}}\exp\left[-\frac{(J_{ij})^{2}N}{ 2\tilde{J}^{2}}\right],
\label{two1}
\end{equation}
where the factor $N$ insures a sensible thermodynamic limit.

The 3-state Potts spin glass was described earlier in our papers
(see, e.g., Ref.\cite {L, Lu, Gribova}) using the representation in
terms of the operators of quadrupole momenta:
$Q=\left(3{J_{z}}^{2}-2\right)$,
$V=\sqrt{3}\left({J_{x}}^{2}-{J_{y}}^{2}\right)$, $\textbf{J}=1$,
$J_{z}=1,0,-1$. The operators $Q$ and $V$ commute and in the representation
where both of them are diagonal have the following form
$$
Q =
\left(
\begin{array}{ccc}
-2& 0& 0\\
\phantom{x}0& 1& 0\\
\phantom{x}0& 0& 1\\
\end{array}
\right)
\hspace{10mm}
V =
\left(
\begin{array}{cccc}
0& \phantom{x} 0& 0\\
0& \phantom{x} \sqrt{3}& 0\\
0& \phantom{x} 0& - \sqrt{3}
 \end{array}
\right)
$$
Throughout this paper we shall use the short-cut notations
 $Q_{k} =(-2,1,1)$ and $V_{k }=(0,\sqrt{3},-\sqrt{3})$.

The equation
\begin{equation}
\left[Q_{k
}Q_{l}+V_{k }V_{l }\right]+2=6\delta_{k ,l }.
\label{two13}
\end{equation}
provides the equivalence of the Potts model with $p=3$ to the model with
the Hamiltonian
\begin{equation}
 H=-\frac{1}{2}\sum_{i\neq j}\left[J_{ij}Q_{i}Q_{j}+G_{ij}{V}_{i}{V}_{j}\right] ,
\label{two834}
\end{equation}
(with $J_{ij}\equiv G_{ij}$).
The quenched interactions are distributed with Gaussian
probability.

Let us note that this representation is very useful for actual
calculations and allowed us to obtain, in particular, the number of
metastable states at zero temperature
\cite{Lu} and to determine the low--temperature boundary of the stability
of the 1RSB solution
\cite{Gribova}.

Using the replica method we obtain the averaged
free energy of the system in the form

\begin{eqnarray}
\label{two934}
 F/NkT=\lim_{n \rightarrow
0}\frac{1}{n}\left[-2t^{2}n+\frac{t^{2}}{2}\sum_{\alpha
\neq\beta}\left(q^{\alpha \beta }\right)^{2} \right.\nonumber
\\
\left.-\ln\Tr\exp\left[\frac{t^{2}}{2}\sum_{\alpha \neq\beta}q^{\alpha
 \beta }
\left({Q}^{\alpha}{Q}^{\beta}+{V}^{\alpha}{V}^{\beta}\right)\right]
\right],
\end{eqnarray}
where $q^{\alpha \beta }$ is the glass order parameter, $t={\tilde{J}}/kT$.
Here indices $\alpha$ and $\beta$  label replicas.
The Potts symmetry permits to describe the glass state with one
glass order parameter for both kinds of operators.

Using the standard procedure  we perform the first stage
of the replica symmetry breaking (1RSB) according to Parisi ($n$ replicas
are divided into  $n/m$ groups with  $m$ replicas in each group) and obtain
the free energy in the form (with $q^{\alpha \beta }= r_1$ if $\alpha $
and $\beta $ are from different groups and $q^{\alpha \beta }= r_1+v$ if
$\alpha $ and $\beta $ belong to the same group).

\begin{eqnarray}
\label{two137}
F_{\mathrm{1RSB}}=-NkT\times \nonumber
\\
\biggl\{t^2\left(2+(1-m)\frac{v^{2}}{2}-2v^{2}\right)+ \nonumber
\\
\frac{1}{m}\ln\int dz^G\int ds^G \left[\Psi\right]^{m}\biggr\}.
\end{eqnarray}

Here
\begin{eqnarray}
\Psi=\exp{\left(-2{\theta}_{1}\right)}+\exp{\left({\theta}_{1}\right)}\left[\exp{\left({\theta}_{2}\right)}
+\exp{\left(-{\theta}_{2}\right)}\right],\\
{\theta}_{1}=zt\sqrt{v},{\theta}_{2}=st\sqrt{3v},
\\
\int\frac{dz}{\sqrt{2\pi}}(\ldots)\exp\left(-\frac{z^2}{2}\right)\equiv \int dz^G(\ldots).
\end{eqnarray}

The Potts symmetry provides also the trivial solution
$q_{RS}=0$ and $r_{1}=0$ (see e.g.
Ref \cite{Ritort}).
Let us note that in general case of the Hamiltonian
(\ref{two834}) with $J_{ij}\neq G_{ij}$ one of the glass order parameter
can not be zero. This fact can be easily seen by analyzing the
high--temperature expansions for the RS equations for the order
parameters \cite{LuTa87}. Let us note that
the Eq.(\ref{two137}) can be obtained from the appropriate equation of
Ref. \cite{Ritort} by the integration over suitable linear
combinations of variables.

Expanding the free energy (\ref{two137}) near $T_c$ up
to $v^4$ and solving the equations that define the extremum
conditions for the free energy one can obtain~\cite{L}:
\begin{equation}
v = 8\tau -84\tau ^2,\hspace{3mm} m = \frac{1}{2} - \frac{9}{2} \tau,
\hspace{3mm}\Delta F= 16\tau ^3,
\label{mvf3}
\end{equation}
where $\Delta F=F_{1RSB}-F_{RS}$ and $ \tau = t-t_c$.

Let us consider now the case $p=4$. In this case one can use the following
operators $Q'\sim 3{J_{z}}^{2}-J(J+1)$, $V'\sim J_{z}$
for $\textbf{J}=3/2$
and the third operator $P'$ defined as to be orthogonal to $Q'$ and $V'$.
Namely, $2Q'_k= (1,-1,-1,1)$, $2\sqrt{5}V'_k=
(-3,1,-1,3)$, $2\sqrt{5}P'_k= (1,3,-3,-1).$

It is easy to see that the following equations hold:
$$Q'^2=\frac{1}{4};~~
V'^2=\frac{1}{4}+\frac{2}{5}Q';~~P'^2=\frac{1}{4}-\frac{2}{5}Q',$$
$$V'P'=-\frac{3}{10}Q';~~Q'V'= -\frac{3}{10}P'+\frac{2}{5} Q'; $$
$$Q'P'=\frac{2}{5}P'-\frac{3}{10}Q',$$
so that
\begin{equation}
Q'_kQ'_l+V'_kV'_l+P'_kP'_l +\frac{1}{4} = \delta _{k,l}.
\label{4444prs}
\end{equation}
There is no reflection symmetry and we have
\begin{equation}
\Tr P'^2Q'=-\frac{2}{5};\hspace{2mm}\Tr V'^2Q'=\frac{2}{5};\hspace{2mm}\Tr P'V'Q'=-\frac{3}{10},
\label{kuby}
\end{equation}
that leads to the appearance of cubic terms in the free energy
expansion\footnote{It is worth to notice that our representation in terms of three
matrices is not unique. For example, one can use the set: $A_k=(1,-1,-1,1)$,
$B_k=(-1,1,-1,1)$, $C_k= (1,1,-1,-1).$ We use, however, the
quadrupole--like operators in order to retain the physical meaning of
generalized anisotropic models and to preserve the analogy between even--p
and odd--p cases.}.

Now, using the replica method one can write the averaged free energy for
the Hamiltonian
\begin{equation}
 H=-\frac{1}{2}\sum_{i\neq
j}J_{ij}[Q'_{i}Q'_{j}+V'_{i}V'_{j}+P'_iP'_j] ,
\label{ham4}
\end{equation}
in the form analogous to (\ref{two934}) and perform the 1RSB.
The extremum conditions for this 1RSB free energy have the form of three
equations that look different only at the first glance.
In fact they are the
equations from Ref.\cite{Ritort}, integrated over different
variables and so have a common solution $v\equiv \langle \langle
Q'_\alpha Q'_\beta \rangle \rangle = \langle \langle V'_\alpha
V'_\beta \rangle \rangle = \langle \langle P'_\alpha P'_\beta
\rangle \rangle$. Expanding the 1RSB
free energy in the neighborhood of $t_c$ up to
$v^4$ we obtain (in terms of standard variables
(as, e.g., in \cite{Ritort})):
\begin{eqnarray}\label{df4}
\Delta F/Nkt= -\frac{1}{8}(-1+m)t^2 v^2 \times \nonumber
\\
\times [-6+6t^2+4(m-1)t^4v+(3m^2-3m-7)t^6v^2].
\end{eqnarray}
Here $t=t_c+\tau $, $m=m_0+\delta m$ and $\tau $, $v$ and $\delta m$ are small.
The extremum conditions relative to $v$ and $m$ give the system of
equations:
\begin{eqnarray}\label{min4}
-\frac{1}{2}(-1+m)t^2v\times \nonumber
\\
\times[-3 +3t^2+3(m-1)t^4v+(3m^2+3m-7)t^6v^2]=0;\nonumber
\\
\frac{1}{8}t^2v^2[6-6t^2-8(m-1)t^4v+(4+12m-9m^2)t^6v^2]=0.
\end{eqnarray}
To determine the relative orders of the small values
$\tau $, $v$ and $\delta m$, one has to use the so-called bifurcation
equation (see, e.g.
Ref.\cite{tren}) for the system (\ref{min4}). One can also check
directly that the equations
(\ref{min4}) become incompatible if
$v$ and $\delta m$  are supposed to be series in integer powers of
$\tau $. If $\delta
m=m_1 \sqrt{\tau }$ and $v=v_1\sqrt{\tau }$ (as it follows from the
bifurcation equation)
then we obtain $t_c=1$, $m_0=1$, and for $m_1$ and $v_1$ the
following system of equations can be obtained from (\ref{min4}):
$$6+3m_1v_1-7v_1^2=0;$$
$$12+8m_1v_1-7v_1^2=0.$$
The system has two solutions that differ only by the sign.
Taking into account that
$$\Delta F/NkT=\frac{1}{2}m_1^2v_1^3\tau ^{5/2},$$
and that on the physical branch $m<1$,
we choose
$$ v=\sqrt{12\tau /35}; ~~~m=1-\sqrt{21\tau /5}.$$
So, in the case of Potts spin glass with 4 states
$m=1$ at the point of bifurcation $t_c$, that is $t_{1RSB}=t_c$.

Let us consider now the Potts glass with
$p=5$ in an analogous way.
We shall use the following diagonal operators:
$P_k^{(1)}=(-4,1,1,1,1)$,
$P_k^{(2)}=(0,-3,1,-1,3)$, $P_k^{(3)}=\sqrt{5}(0,1,-1,-1,1)$,
$P_k^{(4)}=(0,1,3,-3,-1)$. The equivalence to the Potts model
follows from the relation:
$$P_k^{(1)}P_l^{(1)}+P_k^{(2)}P_l^{(2)}+P_k^{(3)}P_l^{(3)}+P_k^{(4)}P_l^{(4)}+4 =
20\delta _{k,l}.
$$
The nonzero cubic terms are
$${P^{(1)}}^3=-12;~~{P^{(2)}}^2P^{(1)}=4;~~{P^{(3)}}^2P^{(1)}=4;$$
 $${P^{(4)}}^2P^{(1)}=4;~~
{P^{(2)}}^2P^{(3)}=16/\sqrt{5};$$
$${P^{(4)}}^2P^{(3)}=-16/\sqrt{5};~~ P^{(1)}P^{(3)}P^{(4)}=-12/\sqrt{5}$$

Now the 1RSB free energy expansion near
$t_c$ has the following form (in standard variables):
\begin{eqnarray}\label{df5}
\Delta F/NkT= -\frac{1}{12}(-1+m)t^2 v^2 \times \nonumber
\\
\times [-12+12t^2+4(2m-1)t^4v+(6m^2+6m-35)t^6v^2].
\end{eqnarray}
The extremum conditions are
\begin{eqnarray}\label{min5}
-\frac{1}{3}(-1+m)t^2v\times \nonumber
\\
\times[-6 +6t^2+3(2m-1)t^4v+(6m^2+6m-35)t^6v^2]=0, \nonumber
\\
\frac{1}{12}t^2v^2[12-12t^2-4(4m-3)t^4v+(41-18m^2)t^6v^2]=0.
\end{eqnarray}

It is easy to obtain that now
$v$ and $\delta m$ can be presented as series in integer degrees of
$\tau $:
$v=v_1\tau +v_2 \tau ^2$, $m=m_0+m_1\tau $.
Here $v_1$, $v_2$, $m_0$ and $m_1$ satisfy the following
equations:
$$6-3v_1+4m_0v_1=0,$$
$$4-v_1+2m_0v_1=0$$
and
$$6v_2-12m_1-92=0,$$
$$6v_2-16m_1-43=0,$$
so that finally we obtain
$$t_c=1;\hspace{2mm}m=\frac{3}{2}+\frac{49}{4}\tau ;\hspace{2mm} v=-2\tau
+\frac{239}{6}\tau ^2$$
and the transition from RS to 1RSB
can not take place at the bifurcation point.

Here we would like to make a remark about the 5-state model without any
pretend to be rigorous. Basing on the fact that the numerically obtained
\cite{Ritort} value of $t_{1RSB}$ for $m=1$ is very close to our $t_c$, one
can hope to obtain the jump to the 1RSB solution already from the
equation (\ref{df5}) and to estimate qualitatively the characteristics in
question. In fact, if we expand (\ref{df5}) in the vicinity of $m=1$ and
let its first and second derivatives relative to $v$ be equal to zero,
we obtain $t_{1RSB}^2= 184/187$ and the jump $\Delta v = 561/8464$ (compare
with the table in \cite{Ritort}).

Qualitative behavior of the solutions for $v$ near $t_c$ for $p=3,4,5$ is
presented on the figure \ref{fig:scheme}.
\begin{figure}[t]
\begin{center}
\includegraphics[width=8 cm]{potts1.eps}
\caption{(Colour online) Schematic representation of the order parameter $v$ in the vicinity of
the bifurcation point $t_c$ for $p=3,4,5$. Solid lines for $m\leq 1$,
dotted lines for $m>1$. }
\label{fig:scheme}
\end{center}
\end {figure}

Now let us consider the 1RSB free energy expansion near the
bifurcation point $t_c$ from a more general point of view.
It is interesting
that the form of the series for $\Delta F$ in small
deviations $\delta q^{\alpha \beta}$ from $q_{RS}$ up to the third order
is one and the same for different models and coincides with the one that can be
written \cite{Hydrogen} for the random generalized Hamiltonian (\ref{SK})
where $U$ is arbitrary diagonal operator. The fact is that in the most
general case (including
(\ref{SK}) and (\ref{one1})) $\Delta F$ depends on the following (the
only nonzero) sums:
\begin{eqnarray}
&\lim_{n \rightarrow 0}\frac{1}{n}
{\sum_{\alpha,\beta}}^{'}(\delta
q^{\alpha\beta})^{2} =-[r-(m-1)v]^{2}-m(1-m)v^{2},\label{11eq:JP}\\
&\lim_{n \rightarrow 0}\frac{1}{n}
{\sum_{\alpha,\beta,\delta}}^{'}\delta
q^{\alpha\beta}\delta
q^{\alpha\delta} =[r-(m-1)v]^{2},
\label{51eq:JP}\\
&\lim_{n \rightarrow 0}\frac{1}{n}
{\sum_{\alpha,\beta}}^{'}(\delta
q^{\alpha\beta})^{3}=-[r-(m-1)v]^{3}+\nonumber\\
&3m(m-1)[r-(m-1)v]v^{2}+m(m-1)(2m-1)v^{3}, \label{61eq:JP}\\
&\lim_{n \rightarrow
0}\frac{1}{n}{\sum_{\alpha,\beta,\gamma}}^{'}\delta
q^{\alpha\beta}\delta q^{\beta\gamma}\delta
q^{\gamma\alpha}=2[r-(m-1)v]^{3} - \nonumber\\
&3m(m-1)[r-(m-1)v]v^{2}-m^{2}(m-1)v^{3}, \label{711eq:JP}\\
&\lim_{n \rightarrow
0}\frac{1}{n}{\sum_{\alpha,\beta,\gamma}}^{'}(\delta
q^{\alpha\beta})^{2}\delta
q^{\alpha\gamma}=[r-(m-1)v]^{3}- \nonumber
\\
&[r-(m-1)v]v^{2}, \label{71eq:JP}\\
&\lim_{n \rightarrow
0}\frac{1}{n}{\sum_{\alpha,\beta,\gamma,\delta}}^{'}\delta
q^{\alpha\beta}\delta q^{\alpha\gamma}\delta q^{\alpha\delta}=\nonumber\\
&\lim_{n \rightarrow
0}\frac{1}{n}{\sum_{\alpha,\beta,\gamma,\delta}}^{'}\delta
q^{\alpha\beta}\delta q^{\alpha\gamma}\delta q^{\beta\delta}=
-[r-(m-1)v]^{3},\label{81eq:JP}
\end{eqnarray}
with  $r=r_{1}-q_{RS}$.
The prime on the sum means that only the superscripts belonging to the same
 $\delta q$ are necessarily different in  $\sum'$ .

So, the deviation $\Delta F$ of the free energy $F_{1RSB}$ from its
replica symmetric part  in the most general case is

\begin{eqnarray}\label{1011frs}
\frac{\Delta F}{NkT}=\frac{t^2}{4}\left[1-t^{2}W\right]\left\{-\left[r-(m-1)v\right]^{2}-\right.
\nonumber\\
\left.v^{2}m(1-m)\right\}-\frac{t^{4}}{2}L\left[r-(m-1)v\right]^{2}- \nonumber\\
t^{6}\left\{C\left[r-(m-1)v\right]^{3}+\right. \nonumber
\\
D\left[r-(m-1)v\right]v^{2}m(m-1)-B_{3}v^{3}m^{2}(m-1)+\nonumber
\\
\left. B_{4}v^{3} m(m-1)(2m-1)\right\}+...
\end{eqnarray}
where
$t=t_c+\tau $, and  parameters
$W, L, C, D, B_3,
  B_4$ are some combinations of operators averaged over the RS solution.
 For example, the coefficient $L$ enters $\Delta F$ with the sum
 (\ref{51eq:JP}): $\frac{1}{n}
{\sum_{\alpha,\beta,\delta}}^{'}\delta q^{\alpha\beta}\delta
q^{\alpha\delta}$

Let us note that in the case of zero RS solutions for the order
parameters the expansion does not contain the terms where
some indices occur only once. In the case of reflection symmetry
there is no terms where some indices occur odd number of
times. If there is no reflection symmetry,
using
the extremum conditions for the free energy for Hamiltonian (\ref{one1})
and taking into account the fact that $L|_{t=t_c}\neq {0}$,  the
bifurcation condition gives\cite{Hydrogen}
$r-(m-1)v=0+o(\Delta t)^{2}$, i.e. the condition that there  is no linear
term for the glass order parameters. In fact, there is no other linear
term because
\begin{equation}
\left[1-t^{2}W\right]|_{t=t_c}=\lambda_{\rm RS repl}|_{t=t_c}=0
\label{40510prs}
\end{equation}
at the bifurcation point.

In the case of the Potts spin glass model the situation is quite
different, although the reflection symmetry is absent.
Now $L=0$ because it is zero RS solution that bifurcates
 so that single indices mean that $L$ is composed of zero values of
averaged
 $\Tr Q$,$\Tr V$.... Analogously the coefficients in front of
(\ref{71eq:JP}) and (\ref{81eq:JP}) are also zero. So, we have
$C=2B_{3}-B_{4}$ and $D=3(-B_{3}+B_{4})$. Here $B_{3}$ is the coefficient
in front of
(\ref{711eq:JP}), and $B_{4}$ that in front of (\ref{61eq:JP}).
Since $L|_{t=t_{0}}= {0}$, then the condition $r-(m-1)v=0$ is not fulfilled.
(If $p>2$ then $r_{1}=r=0$.\cite{Ritort})

In general case we obtain from the extremum conditions for the
free energy (\ref{1011frs}) the following equations for the order
parameters:
\begin{eqnarray}
& 2m(m-1)vZ\Delta t =
t_c^{6}m(m-1)v\left\{3\left[-B_{4}+\right.\right.\nonumber\\
&\left.\left. m(-B_{3} +2B_{4})\right]v+2D\left[r-(m-1)v\right]\right\} \label{30prs},\\
\nonumber \\
&(2m-1)v^{2}Z\Delta t =
t_c^{6}v^{2}\left\{(2m-1)\left[-B_{4}+\right.\right. \nonumber \\
&\left.m(-B_{3}+2B_{4})\right]v+m(m-1)(-B_{3}+2B_{4})v+\nonumber\\
&\left.D(2m-1)\left[r-(m-1)v\right]\right\}
\label{40prs},
\end{eqnarray}
where
\begin{equation}\label{B123eq:JP}
Z = \frac{d}{dt}{\left[\frac{t^2}{4}\left[1-t^{2}W\right]
\right]}_{\mid_{t_c}}.
\end{equation}

Hence we obtain from (\ref{30prs})-(\ref{40prs}):
\begin{equation}
m(m-1)\left[mB_{3}-B_{4}\right]=0.
\label{401prs}
\end{equation}

In fact, in the case of the Potts spin glass with $p=3$ for the Hamiltonian
 (\ref{two834})
at $t=t_{c}$ we have:
\begin{eqnarray}
W &= \frac{1}{2}\left[\left(TrQ^2/3\right)^{2}+
\left(TrV^2/3\right)^{2}\right]=4,\nonumber\\
B_{3}&=\frac{1}{6}\left[\left(TrQ^2/3\right)^{3}+
\left(TrV^2/3\right)^{3}\right]=\frac{8}{3},\nonumber\\
B_{4}&=\frac{1}{12}\left[\left(TrQ^3/3\right)^{2}+
3\left(TrQV^2/3\right)^{2}\right]=\frac{4}{3},\nonumber\\
D&=-4.
\label{0298prs}
\end{eqnarray}

In this case in the neighborhood of the bifurcation point one can get from the
equations (\ref{30prs})-(\ref{40prs}) that $v\sim \tau $,  and from the
 equation (\ref{401prs}) $m=\frac{1}{2}$.

Analogously, in the case of 4 states (for the Hamiltonian (\ref{ham4}) we
 obtain: $D=0$, $B_{3} = B_{4} = \frac{1}{128}$ so that $m=1$. The
 equation (\ref{30prs}) is the identity now. The r.h.s. of (\ref{40prs})
 becomes zero and one has to take into account higher order terms in the
 expansion of $F_{1RSB}$. This leads to $v\sim\sqrt{\tau }$. For $p=5$ we
 have $D=64, B_{3}=\frac{128}{3}, B_{4}=64$, so that $m=\frac{3}{2}$ and
 $v\sim -\tau $.

 It is worth noticing that from the free--energy expansion
 (\ref{1011frs}) it follows directly that the stability of the RS state
 is determined by the sign of
 $\lambda_{\rm RS repl}$, just as in the case of SK model
 \cite{Almeida}. To prove this fact let us consider small deviations from
 the RS solution:
$q^{\alpha  \beta}=q_{RS}+\delta q^{\alpha \beta}$
at arbitrary temperature and let us perform 1RSB. Using the notations
introduced before, (that is
  $\delta q^{\alpha \beta }= r$ if $\alpha $ and $\beta $ are
 from different groups and $\delta q^{\alpha \beta }= r+v$ if $\alpha $
and $\beta $ belong to the same group) we can write $F_{1RSB}-F_{RS}$
up to the second order in $\delta q$ in the form
(\ref{1011frs}) with arbitrary $t$. Since
$\left[1-t^{2}W\right]=\lambda_{\rm RS repl}$, it is easy to see
that for $m<1$ and $L \geq {0}$ the RS state is stable ($\Delta F<0$)
for the values of temperature such that
$\lambda_{\rm RS repl}>0$. The similar result is valid for all
subsequent stages of RSB because of
the Parisi rule
$$\lim_{n \rightarrow 0}\frac{1}{n}
{\sum_{\alpha,\beta}}^{'}(\delta
q^{\alpha\beta})^{2}<0$$

 To conclude, we introduced the representation for the Potts glass model
based on the operators of the quadrupole momenta. With the help of this
representation the expansions for the free energy functionals of the Potts spin glass models
with 3, 4 and 5 states up to the fourth order in $\delta q_{\alpha \beta }$
around the replica symmetric solution were obtained. The temperature
dependence of the 1RSB order parameters in the vicinity of the point
$T=T_c$ where the RS solution becomes unstable was derived. The crossover
from continuous to jumpwise behavior with the growing of the number of
states was traced analytically. The comparison was made of the free energy
expansions for the Potts spin glass with that for other models and the
similarity of such expansions was demonstrated.

Authors thank V.N. Ryzhov and N.M.Chtchelkatchev for
helpful discussions and valuable comments.

This work was supported in part by the Russian
Foundation for Basic Research (Grant No. 08-02-00781),
 and the Program of the Presidium of the Russian Academy of
Sciences.
\section*{References}

redpajama_set_name: RedPajamaArXiv