# Solving high-dimensional Fokker-Planck equation with functional hierarchical
tensor

Xun Tang  Lexing Ying

###### Abstract

This work is concerned with solving high-dimensional Fokker-Planck equations
with the novel perspective that solving the PDE can be reduced to independent
instances of density estimation tasks based on the trajectories sampled from
its associated particle dynamics. With this approach, one sidesteps error
accumulation occurring from integrating the PDE dynamics on a parameterized
function class. This approach significantly simplifies deployment, as one is
free of the challenges of implementing loss terms based on the differential
equation. In particular, we introduce a novel class of high-dimensional
functions called the functional hierarchical tensor (FHT). The FHT ansatz
leverages a hierarchical low-rank structure, offering the advantage of
linearly scalable runtime and memory complexity relative to the dimension
count. We introduce a sketching-based technique that performs density
estimation over particles simulated from the particle dynamics associated with
the equation, thereby obtaining a representation of the Fokker-Planck solution
in terms of our ansatz. We apply the proposed approach successfully to three
challenging time-dependent Ginzburg-Landau models with hundreds of variables.

###### keywords:

High-dimensional Fokker-Planck equation; Functional tensor network;
Hierarchical tensor network; Curse of dimensionality.

††journal:

[stanford]organization=Stanford University,city=Stanford, state=California,
postcode=94305, country=USA

##  1 Introduction

Solving high-dimensional partial differential equations (PDE) remains one of
the main challenges in scientific computing, with applications in molecular
dynamics, quantum mechanics, and high-dimensional control theory. A key
element underlying the challenges in solving high-dimensional PDE is the fact
that the computational resources needed to carry out traditional numerical
techniques often scale exponentially in the dimension of the problem [1]. This
paper focuses on solving the Fokker-Planck equation:

$\partial_{t}p=\nabla\cdot(p\nabla V)+\frac{1}{\beta}\Delta p,\quad
x\in\Omega\subset\mathbb{R}^{d},t\in[0,T],$ (1)

which governs the distribution function $p(t,x)$ of the particle system
governed by the Langevin dynamics $dX_{t}=-\nabla
V(X_{t})dt+\sqrt{2\beta^{-1}}dB_{t}$ over a potential function $V$ with
inverse temperature $\beta$.

Solving Fokker-Planck equations in low dimensions is a well-developed topic.
This paper is concerned with the case of hundreds or even thousands of
dimensions. One prevalent type of high-dimensional Fokker-Planck equation
comes from the discretization of an infinite-dimensional functional equation
into finite but very large dimensions. Throughout this text, we use the
motivating example of a time-dependent Ginzburg-Landau (G-L) model used for
studying the phenomenological theory of superconductivity [2, 3, 4, 5]. In the
2D G-L model, one can intuitively think of a particle taking the form of a
field $x(a)\colon[0,1]^{2}\to\mathbb{R}$, where the potential function is a
functional $V(x)$ defined as follows:

$V(x)=\frac{\lambda}{2}\int_{[0,1]^{2}}|\nabla_{a}x(a)|^{2}\,da+\frac{1}{4\lambda}\int_{[0,1]^{2}}|1-x(a)^{2}|^{2}\,da,$
(2)

where $\lambda$ is a parameter balancing the relative importance of the first
term and the second term in (2). In the numerical treatment, one considers the
discretization of the unit square $[0,1]^{2}$ into a grid of $d=m^{2}$ points
$\\{(ih,jh)\\}$, for $h=\frac{1}{m+1}$ and $1\leq i,j\leq m$. The
discretization of the field at this grid is the $d$-dimensional vector
$x=(x_{(i,j)})_{1\leq i,j\leq m}$, where $x_{(i,j)}=x(ih,jh)$. Under the
discretization scheme, the potential energy is defined as

$V(x)=V(x_{(1,1)},\ldots,x_{(m,m)}):=\frac{\lambda}{2}\sum_{v\sim
w}\left(\frac{x_{v}-x_{w}}{h}\right)^{2}+\frac{1}{4\lambda}\sum_{v}\left(1-x_{v}^{2}\right)^{2},$
(3)

where $v$ and $w$ are Cartesian grid points and $v\sim w$ if and only if they
are adjacent. The evolution of the distribution over discretized fields
satisfies the Fokker-Planck equation (1) with the discretized potential
function $V(x)$. Even a moderate grid resolution results in a model very large
in dimension, which is a recurring theme for discretizing infinite-dimensional
functional differential equations. The associated functional differential
equation is the equation for the functional $p(t,x)$ encoding the probability
density at the function $x$ in time $t$.

###  1.1 Background and contribution

#### Solving Fokker-Planck equations without time-stepping

Traditionally, numerical solvers for the high-dimensional Fokker-Planck
equations rely on performing time integration of a given initial condition
$p|_{t=0}$ along a chosen parametric function class. The time integration can
be done either by directly performing the time stepping of the solution or it
can be done by including the PDE loss in the regression target. A common
feature of the time integration scheme is error accumulation, as errors from
previous time steps tend to build up. Moreover, compressing the differential
operator $\nabla\cdot(p\nabla V)$ in (1) is often performed with heuristics,
leading to quite significant implementation challenges both in the design of
the compression scheme and in its development into code.

The key insight in this paper is that solving the high-dimensional time-
dependent Fokker-Planck equation can be reduced to a series of density
estimation tasks. The reasoning is rather simple: one can approximate the
Fokker-Planck solution $p(t,x)$ at any time $t\in[0,T]$ through density
estimation over the sampled trajectory at time $t$. Thus, one can approximate
the solution $p(t,x)$ at a collection of time steps $0=t_{0}\leq\cdots\leq
t_{K}=T$, and the continuous-in-time solution $p(t,x)$ for an arbitrary
$t\in[0,T]$ can be approximated through appropriate interpolation. More
specifically, we take advantage of the associated particle dynamics to the
Fokker-Planck equation:

$dX_{t}=-\nabla V(X_{t})dt+\sqrt{2\beta^{-1}}dB_{t}.$ (4)

By simulating the dynamics of (4) through stochastic dynamic equation (SDE)
simulation, we obtain a particle-based empirical approximation of $p(t,x)$ by
independently sampled trajectory data $\\{X_{i}(t)\\}_{t\in[0,T]}$ for
$i=1,\ldots,N$, which constitutes the input to the density estimation task.
Thus, if $p(t,x)$ belongs to a function class with an efficient denoising or
density estimation procedure, we can approximate the Fokker-Planck solution
via denoising this empirical distribution. In effect, this approach sidesteps
both the error accumulation and the coding challenges of a time integration
scheme.

Here, we choose the functional hierarchical tensor as the ansatz of the PDE
solution $p(t,x)$. The rationale for this choice is that the diffusive nature
of the Fokker-Planck equations tends to keep the tensor rank under control.
Compared to most neural network algorithms, the proposed ansatz can
efficiently calculate the normalization constant (see Section 1.2 for detailed
discussions). Moreover, while directly modeling the density through
nonparametric density estimation suffers from the curse of dimensionality [6],
our proposed approach effectively performs by denoising the empirical
distribution through restriction to the proposed parametric function class,
thus avoiding the curse of dimensionality. Among the possible choices for
functional tensor networks, while an alternative is to perform density
estimation with a functional tensor train [7, 8, 9, 10, 11], we choose a
functional hierarchical tensor to better account for low-rank structures of
more general high-dimensional Fokker-Planck equations.

#### Functional hierarchical tensor

Functional hierarchical tensor is understood in the literature [12, 13] to be
a suitable tool for reduced-order modeling due to its versatility for
multidimensional systems. In particular, the ansatz uses a modeling assumption
natural for modeling discretized functional differential equations. The
hierarchical tensor network assumes a low-rank structure along a hierarchical
bipartition of the variables, which readily applies to discretized infinite-
dimensional models. For example, the 2D G-L model has the natural geometry of
a 2D grid, and it is natural to construct hierarchical bipartition of the
variables through alternatively partitioning the variables along the two axes
of the 2D grid, as shown in Figure 1. In comparison, it is known that the
tensor train ansatz assumes the variables $x=(x_{1},\ldots,x_{d})$ satisfy
low-rankness according to the 1D structure set by the tensor network [14],
which makes it unsuitable for the 2D Ginzburg-Landau model, and the same
argument applies to models for which the variables do not obey a 1D
topological structure.

Our choice of functional hierarchical tensor over functional tensor train is
partially motivated by parallel developments in quantum physics. For the task
of modeling quantum systems, the tensor train ansatz is more commonly known as
matrix product states [15], where it is well-known that it is more suitable
for modeling the ground state of a 1D system [16]. For modeling systems with
nonlocal interaction, a hierarchical tensor network is often more advantageous
[17]. Through the use of a functional tensor network, our work extends the
capability of the hierarchical tensor network into modeling continuous
distributions.

Figure 1: Illustration of the hierarchical bipartition of a 2D Ginzburg-Landau
model discretized to an $8\times 8$ grid.

#### Density estimation through hierarchical sketching

After the trajectories are sampled through particle simulations, we process
these trajectories through a sketching-based density estimation algorithm,
obtaining a hierarchical tensor representation of the Fokker-Planck equation
solution. The main contribution of this work is an end-to-end procedure that
produces a Fokker-Planck solution from sampled trajectories.

To ensure that an accurate approximation can be obtained with a small
parameter size, we have worked out practical choices over key design
questions. The important ingredients are the structure of the hierarchical
bipartition, the sketch function, and the basis to be used for the
hierarchical tensor. Using sophisticated numerical techniques, we are able to
obtain solutions to Ginzburg-Landau models of very high dimensions, including
solving the multidimensional Ginzburg-Landau model in (3) with $d$ in the
hundreds.

###  1.2 Related work

#### Neural network for solving high-dimensional PDEs

Neural network has shown tremendous success in recent years in solving high-
dimensional PDEs. The most common approach to solve PDEs with neural networks
primarily relies on penalizing the partial-differential equation either in
terms of the strong form (as in the physics-informed neural network [18] and
the deep Galerkin method [19]) or the variational form (as in the deep Ritz
method [20]). Another separate direction is to solve the PDE by creating
regression targets special to the problem itself, as can be seen in the deep
BSDE method [1] and neural operators [21]. While these algorithms perform
well, further improving their effectiveness in high-dimensional PDEs is
challenging.

#### Neural network methods for density estimation

As our work primarily focuses on using density estimation to solve the Fokker-
Planck equation, our work is related to generative learning algorithms for
which density estimation is possible. Examples include normalizing flows [22,
23], energy-based models [24, 25], and the more recent diffusion models [26,
27]. Only normalizing flows allow direct density evaluation. In contrast, the
energy-based model and diffusion model only allow for approximating the
likelihood function, which renders them inapplicable in our case, as the
objective is to obtain the model’s associated normalized likelihood function.
Moreover, training several neural networks across many time steps might be
inefficient, while our sketching-based approach relies on quick linear
algebraic subroutines and does not involve iterative optimization.

#### Monte-Carlo method

As one has access to trajectories of the particle dynamics due to Monte-Carlo
sampling, one might hope to perform density estimation through nonparametric
density estimation directly. However, particle-based density estimations such
as kernel density estimation suffer from the curse of dimensionality [6],
making it ill-suited to our problem. If the goal is to estimate the
statistical moments of the solved Fokker-Planck equation, then indeed the
traditional Monte-Carlo approach is well-suited to estimate such quantity with
an $O(1/\sqrt{N})$ rate [28], $N$ being the number of samples. We remark that
in practical cases, the variance is often too high for the estimation to have
a satisfactory estimation accuracy.

#### Tensor network for high-dimensional PDE

There are several works using tensor networks to solve high-dimensional PDEs.
One approach is to use grid discretization of the PDE solution $p$, resulting
in the task of solving the PDE for a $d$-dimensional tensor, as can be seen in
[29]. The functional tensor network approach employs a mesh-free
representation of the PDE solution based on a tensor network. Along this line
of work, the considered tensor network structures include canonical-polyadic
(CP) [30] and tensor train [12, 13, 11, 9, 31]. The use of hierarchical tensor
in solving high-dimensional PDE has also been discussed [32], though quite
different from our proposed methodology. The work in [33, 7] bears the most
resemblance to our work, as they combine particle methods with a functional
tensor network structure. The main difference is the fact that our proposed
algorithm does not involve a time-stepping component, and therefore, our
method is less prone to error accumulation. A similar sketching-based density
estimation algorithm for solving the Fokker-Planck equation with a functional
tensor train can be done by combining our approach with the sketching methods
outlined in [14].

###  1.3 Contents and notations

We outline the structure of the remainder of the manuscript. Section 2 gives a
detailed introduction to functional hierarchical tensor and its sketching-
based density estimation algorithm. Section 3 goes through the algorithm for
solving the Fokker-Planck equation with the functional hierarchical tensor.
Section 4 details the numerical implementations and results from solving
Ginzburg-Landau models in multiple dimensions.

For notational compactness, we introduce several shorthand notations for
simple derivation. For $n\in\mathbb{N}$, let $[n]:=\\{1,\ldots,n\\}$. For an
index set $S\subset[d]$, we let $x_{S}$ stand for the subvector with entries
from index set $S$. We use $\bar{S}$ to denote the set-theoretic complement of
$S$, i.e. $\bar{S}=[d]-S$.

##  2 Hierarchical functional tensor network

To address the Ginzburg-Landau and general models coming from discretizing
high-dimensional functional equations, we propose to use the function class of
a functional hierarchical tensor based on a binary-tree-based low-rank
structure. In this section, we go over the functional hierarchical tensor and
its sketching algorithm. The symbol $d$ is reserved for the dimension of the
state-space variable, and the symbol $N$ is reserved for the number of
samples. Without loss of generality, we shall go over the procedure to perform
density estimation over a collection of $N$ samples
$\\{y^{(i)}\in\mathbb{R}^{d}\\}_{i=1}^{N}$. One can think of $y^{(i)}$ as the
data collected from the $i$-th trajectory at a fixed time $t\in[0,T]$.

###  2.1 Hierarchical bipartition

We first describe the functional tensor network representation of a
$d$-dimensional function $p\colon\mathbb{R}^{d}\to\mathbb{R}$. In general, let
$\\{\psi_{i}\\}_{i=1}^{n}$ denote a collection of orthonormal function basis
over a single variable, and let $C\in\mathbb{R}^{n^{d}}$ be the tensor
represented by a tensor network. The _functional tensor network_ is the
$d$-dimensional function defined by the following equation:

$p(x)\equiv
p(x_{1},\ldots,x_{d})=\sum_{i_{1},\ldots,i_{d}=0}^{n-1}C_{i_{1},\ldots,i_{d}}\psi_{i_{1}}(x_{1})\cdots\psi_{i_{d}}(x_{d})=\left<C,\,\bigotimes_{j=1}^{d}\vec{\Psi}(x_{j})\right>,$
(5)

where $\vec{\Psi}(x_{j})=\left[\psi_{1}(x_{j}),\ldots,\psi_{n}(x_{j})\right]$
is an $n$-vector encoding the evaluation of any $x_{j}$ over the entire single
variable function basis set.

In this work, we choose $C$ to be represented by a hierarchical tensor network
ansatz so that the associated $p$ takes the form of a functional hierarchical
tensor. In addition to density evaluation, this ansatz allows for efficient
evaluation of moments and efficient sampling (see Section 3 for more details
on potential application).

Crucial to the notion of the hierarchical tensor network is a hierarchical
bipartition of the variable set. Without loss of generality, let $d=2^{L}$ so
that the variable set admits exactly $L$ levels of variable bipartition. As
illustrated in Figure 2(a), at the $l$-th level, the variable index set is
partitioned according to

$[d]=\bigcup_{k=1}^{2^{l}}I_{k}^{(l)},\quad
I_{k}^{(l)}:=\\{2^{L-l}(k-1)+1,\ldots,2^{L-l}k\\},$ (6)

which in particular implies the recursive relation that
$I_{k}^{(l)}=I_{2k-1}^{(l+1)}\cup I_{2k}^{(l+1)}$.

(a) Binary decomposition of variables

(b) Functional hierarchical tensor diagram

Figure 2: Illustrations of functional hierarchical tensor with $d=8$.

###  2.2 Tensor network structure

We introduce the fundamental structure of a functional hierarchical tensor.
Importantly, we focus on a single node $q$, say at level index $l$ and block
index $k\in[2^{l}]$. We take the shorthand
$a(q):=I_{2k-1}^{(l+1)},b(q):=I_{2k}^{(l+1)}$ and $f(q)=[d]-a(q)\cup b(q)$, to
be its left branch, right branch, and the rest, respectively, as illustrated
in Figure 2(b). While these index sets depend on the node $q$, in order to
simplify the formulae, we refer to them just as $a$, $b$, and $f$.

Intrinsically, a hierarchical tensor network represents a tensor
$C\in\mathbb{R}^{n^{d}}$ with a low-rank structure along every hierarchical
bipartition. The bipartition at the node $q$ is defined to be $[d]=\bar{f}\cup
f$, where $\bar{f}$ is the complement of $f$ as a subset of $[d]$. With the
partition $[d]=\bar{f}\cup f$, there exists $r_{f}\in\mathbb{N}$ so that the
$n^{|\bar{f}|}\times n^{|f|}$ unfolding matrix $C(i_{\bar{f}};i_{f})$ is of
rank $r_{f}$. In other words, there exist
$C_{\bar{f}}\colon\mathbb{R}^{n^{|\bar{f}|}}\times\mathbb{R}^{r_{f}}$ and
$C_{f}\colon\mathbb{R}^{r_{f}}\times\mathbb{R}^{n^{|f|}}$ that form a low-rank
decomposition of $C$ (illustrated in Figure 3(a)):

$C(i_{1},\ldots,i_{d})=\sum_{\theta=1}^{r_{f}}C_{\bar{f}}(i_{\bar{f}},\theta)C_{f}(\theta,i_{f}).$
(7)

If the low-rankness property in (7) holds for any node $q$, the exponential-
sized tensor $C$ can be described by a mere $O(dr^{3})$ number of parameters
by a hierarchical tensor network, where $r:=\max_{q}(r_{f(q)})$. The tensor
network is parameterized by a collection of tensor cores denoted by
$\\{G_{q}\\}_{q}$. As Figure 2 illustrates, the tensor cores have the same
tree structure as in the variable bipartition. Each internal bond of the
tensor network is associated with a bond in the binary decomposition of the
variables. Each physical bond of the tensor network connects a leaf node with
a $\vec{\Psi}(x_{j})$ node. In summary, the graphical structure of the
hierarchical tensor network is formed by a binary tree, and the physical index
is at the leaf node. The fact that such cores exist can be proved by induction
using (7), and a proof can be found in [34]. In terms of runtime complexity
and memory complexity, the functional hierarchical tensor enjoys the same
$O(d)$ scaling as that of a hierarchical tensor [34].

###  2.3 Sketching algorithm

We shall go over how to use a sketch-based method to obtain a functional
hierarchical tensor representation from the given collection of samples
$\left\\{y^{(i)}:=\left(y_{1}^{(i)},\ldots,y_{d}^{(i)}\right)\right\\}_{i=1}^{N}$.
For the reader’s convenience, important equations in this subsection are
included in Figure 3 in terms of the tensor diagram. Below, we review the main
equation behind hierarchical sketching in the functional case. The goal is to
solve for the tensor core $G_{q}$, where $q$ is the node of the $k$-th block
at level $l$.

#### Equations

For simplicity, we first assume $0<l<L$ so that $q$ is neither the root nor
the leaf node. Similar as before, we have
$a:=I_{2k-1}^{(l+1)},b:=I_{2k}^{(l+1)}$ and $f=[d]-a\cup b$. Then, the
structural low-rankness property of (7) implies that there exist
$r_{a},r_{b},r_{f}\in\mathbb{N}$,
$C_{a}\colon[n^{|a|}]\times[r_{a}]\to\mathbb{R}$,
$C_{b}\colon[n^{|b|}]\times[r_{b}]\to\mathbb{R}$,
$C_{f}\colon[r_{f}]\times[n^{|f|}]\to\mathbb{R}$ such that the following
equation holds (illustrated in Figure 3(a)):

$C(i_{1},\ldots,i_{d})=\sum_{\alpha,\beta,\theta}C_{a}(i_{a},\alpha)C_{b}(i_{b},\beta)G_{q}(\alpha,\beta,\theta)C_{f}(\theta,i_{f}).$
(8)

Figure 3: Tensor diagram representation of main equations used in Section 2.3.
Equation (7) and (8) are summarized in (a). Equation (9) is shown in (b). The
coefficient term in (9) is shown in (c). The definition of $B_{q}$ is shown in
(d). Equation (10) for $Z_{a;\bar{a}}$ is shown in (e), where the scissor
symbol indicates the rank-$r_{a}$ compression through SVD.

Equation (8) is exponential-sized, and therefore, one would not solve this
equation directly. The hierarchical sketching algorithm essentially solves the
over-determined linear system (8) for $G_{q}$ with the use of sketch
functions. Let $\tilde{r}_{a},\tilde{r}_{b},\tilde{r}_{f}$ be integers so that
$\tilde{r}_{a}>r_{a},\tilde{r}_{b}>r_{b},\tilde{r}_{f}>r_{f}$. Through the
sketch functions
$S_{a}\colon[n^{|a|}]\times[\tilde{r}_{a}]\to\mathbb{R},S_{b}\colon[n^{|b|}]\times[\tilde{r}_{b}]\to\mathbb{R},S_{f}\colon[n^{|f|}]\times[\tilde{r}_{f}]\to\mathbb{R}$,
we can contract the linear system in (8) with $S_{a},S_{b},S_{f}$ at the
variables $i_{a},i_{b},i_{f}$ respectively, which leads to the following
sketched linear system (illustrated in Figure 3(b) and Figure 3(d)):

$B_{q}(\mu,\nu,\zeta)=\sum_{\alpha,\beta,\theta}A_{a}(\mu,\alpha)A_{b}(\nu,\beta)A_{f}(\zeta,\theta)G_{q}(\alpha,\beta,\theta),$
(9)

where $A_{a},A_{b},A_{f}$ are respectively the contraction of
$C_{a},C_{b},C_{f}$ by $S_{a},S_{b},S_{f}$ and $B_{q}$ is the contraction of
$C$ by $S_{a}\otimes S_{b}\otimes S_{f}$, as illustrated by the first equality
of Figure 3(d). The equation (9) can be seen as the following linear system of
$\tilde{r}_{a}\tilde{r}_{b}\tilde{r}_{f}$ equations for the unknown $G_{q}$
(illustrated in Figure 3(d)):

$(A_{a}\otimes A_{b}\otimes A_{f})G_{q}=B_{q}.$

While obtaining $B_{q}$ in (9) is done by contracting known tensors, obtaining
terms such as $A_{a}$ is only possible after specifying the gauge degree of
freedom for $C_{a}$. However, we show that one can obtain $A_{a}$ without
having access to $C_{a}$. The structural low-rankness of the ansatz implies
the existence of $C_{\bar{a}}$ such that there exists a linear system
$C(i_{1},\ldots,i_{d})=\sum_{\alpha}C_{a}(i_{a},\alpha)C_{\bar{a}}(\alpha,i_{{\bar{a}}})$.
One can sketch this linear system by contracting it with
$S_{a},S_{{\bar{a}}}$, which leads to the following linear system (illustrated
in Figure 3(e)):

$Z_{a;{\bar{a}}}(\mu,\psi)=\sum_{\alpha}A_{a}(\mu,\alpha)A_{{\bar{a}}}(\alpha,\psi),$
(10)

where $A_{{\bar{a}}}$ is the contraction of $C_{{\bar{a}}}$ with
$S_{{\bar{a}}}$. Then, due to a gauge degree of freedom in choosing $C_{a}$,
there also exists a gauge degree of freedom in choosing $A_{a}$. Thus, one can
perform singular value decomposition (SVD) on $Z_{a;{\bar{a}}}$, obtaining
$Z_{a;{\bar{a}}}(\mu,\psi)=\sum_{\alpha}U(\mu,\alpha)V(\alpha,\psi)$, and the
choice of $A_{a}=U$ and $A_{{\bar{a}}}=V$ forms a consistent choice of gauge
between the pair $(C_{a},C_{{\bar{a}}})$.

Likewise, one can obtain $Z_{b;\bar{b}}$ and the SVD of $Z_{b;\bar{b}}$
results in $A_{b}$ and $A_{\bar{b}}$. The same holds for $Z_{\bar{f};f}$,
$A_{f}$, and $A_{\bar{f}}$. To solve for $G_{q}$ in the linear system (9), one
simply contracts $B_{q}$ with the pseudo-inverse of $A_{a},A_{b},A_{f}$. Thus,
one can use (9) to solve for $G_{q}$ for any $l\not=0,L$.

Our construction likewise gives rise to the equation of $G_{q}$ for $l=0,L$ as
special cases. For $l=0$, one can go through with the same calculation by
setting $C_{f}=S_{f}=A_{f}=1$ in (8) and (9). For $l=L$, one can sketch the
linear system in (8) for $f=[d]-\\{k\\}$ by contraction with $S_{f}$. The
detail for both cases is described in Algorithm 1.

#### Approximation via samples

We now describe how one can obtain a consistent estimation of $G_{q}$ through
samples. As before, suppose $\\{\psi_{i}\\}_{i=1}^{n}$ is the collection of
univariate orthonormal function basis which defines the functional tensor
network basis. To the sketch tensor
$S_{a}\colon[n^{|a|}]\times[\tilde{r}_{a}]\to\mathbb{R}$, one associates a
continuous sketch function
$s_{a}\colon\mathbb{R}^{|a|}\times[\tilde{r}_{a}]\to\mathbb{R}$ through the
following construction:

$s_{a}(x_{a},\mu)=\sum_{i_{j},j\in a}S_{a}(i_{a},\mu)\prod_{j\in
a}\psi_{i_{j}}(x_{j}),$ (11)

and one can likewise obtain $s_{\bar{a}}$ through $S_{\bar{a}}$. By
orthonormality of the function basis, it is easy to check that

$Z_{a;\bar{a}}(\mu,\psi)=\int_{\mathbb{R}^{d}}p(x_{1},\ldots,x_{d})s_{a}(x_{a},\mu)s_{\bar{a}}(x_{\bar{a}},\psi)\,d\,x_{1}\ldots
x_{d}=\mathbb{E}_{X\sim
p}\left[s_{a}(X_{a},\mu)s_{\bar{a}}(X_{\bar{a}},\psi)\right],$

and the formulae for $(Z_{b;\bar{b}}$ and $Z_{\bar{f};f}$ follow likewise.
Similarly, one can obtain $B_{q}$ in (9) by

$B_{q}(\mu,\nu,\zeta)=\mathbb{E}_{X\sim
p}\left[s_{a}(X_{a},\mu)s_{b}(X_{b},\nu)s_{f}(X_{f},\zeta)\right].$

Since $\\{y^{(i)}\in\mathbb{R}^{d}\\}_{i=1}^{N}$ provides an empirical
approximation to $p$, one obtains a finite sample approximation of
$Z_{a;\bar{a}}$ in the following way:

$Z_{a;\bar{a}}(\mu,\psi)\approx\frac{1}{N}\sum_{i=1}^{N}s_{a}(y^{(i)}_{a},\mu)s_{\bar{a}}(y^{(i)}_{\bar{a}},\psi)$
(12)

and likewise one can also obtain a sample estimation of $B_{q}$ and thus
approximately solve for $G_{q}$. It is easy to check that a consistent
estimator of $Z_{a;\bar{a}}$ leads to a consistent estimator of $A_{a}$, and
therefore one can likewise consistently estimate $A_{b},A_{f}$, and thus the
approximated solution is a consistent estimator for $G_{q}$.

#### Algorithm summary

The sketching algorithm is summarized as Algorithm 1, which demonstrates the
whole procedure to carry out the functional hierarchical tensor sketching
algorithm given a sample collection, including the details for the omitted
edge cases of $l=0,L$.

Algorithm 1 Functional hierarchical tensor sketching.

0: Sample $\\{y^{(i)}\\}_{i=1}^{N}$.

0: Chosen function basis $\\{\psi_{i}\\}_{i=1}^{n}$.

0: Collection of sketch functions
$\\{s_{I_{k}^{(l)}}\\},\\{s_{[d]-I_{k}^{(l)}}\\}$ and target internal ranks
$\\{r_{I_{k}^{(l)}}\\}$ for each level index $l$ and block index $k$.

1: for each node $q$ on the hierarchical tree do

2: Set $l$ as the level index of $q$. Set $k$ as the block index of $q$.

3: if $q$ is not leaf node then

4: $(a,b)\leftarrow(I_{2k-1}^{(l+1)},I_{2k}^{(l+1)})$

5: Obtain $Z_{a;\bar{a}},Z_{b;\bar{b}}$ by (12)

6: Obtain $A_{a}$ as the left factor of the best rank $r_{a}$ factorization of
$Z_{a;\bar{a}}$

7: Obtain $A_{b}$ as the left factor of the best rank $r_{b}$ factorization of
$Z_{b;\bar{b}}$.

8: end if

9: if $q$ is not root node then

10: $f\leftarrow[d]-I_{k}^{(l)}$

11: Obtain $Z_{\bar{f};f}$ by (12).

12: Obtain $A_{f}$ as the right factor of the best rank $r_{f}$ factorization
of $Z_{\bar{f};f}$.

13: end if

14: if $q$ is leaf node then

15: Obtain
$B_{q}(\mu,\nu)=\frac{1}{N}\sum_{i}s_{a}(y^{(i)}_{a},\mu)s_{b}(y^{(i)}_{b},\nu).$

16: Obtain $G_{q}$ by solving the over-determined linear system $(A_{a}\otimes
A_{b})G_{q}=B_{q}$.

17: else if $q$ is root node then

18: Obtain
$B_{q}(j,\zeta)=\frac{1}{N}\sum_{i}\psi_{j}(y^{(i)}_{k})s_{f}(y^{(i)}_{f},\zeta).$

19: Obtain $G_{q}$ by solving the over-determined linear system
$(A_{f})G_{q}=B_{q}$.

20: else

21: Obtain
$B_{q}(\mu,\nu,\zeta)=\frac{1}{N}\sum_{i}s_{a}(y^{(i)}_{a},\mu)s_{b}(y^{(i)}_{b},\nu)s_{f}(y^{(i)}_{f},\zeta).$

22: Obtain $G_{q}$ by solving the over-determined linear system $(A_{a}\otimes
A_{b}\otimes A_{f})G_{q}=B_{q}$.

23: end if

24: end for

##  3 Solving for Fokker-Planck equation

In this section, we detail our approach to solving the Fokker-Planck equation.
The given information regarding the PDE includes the potential function
$V\colon\mathbb{R}^{d}\to\mathbb{R}$, terminal time $T$, and initial
distribution $p_{0}$. The main workflow is summarized in the following simple
steps:

  1. 1.

Select time steps $0=t_{0}\leq t_{1}\leq\ldots\leq t_{K}=T$. Then, with a
chosen numerical scheme (e.,g. Euler-Maruyama), run $N$ independent instance
of SDE simulations on the stochastic dynamic equation (4) to time $T$ with the
initial condition $p_{0}$. Obtain $N$ trajectories and record the trajectory
data $\\{X^{(i)}(t_{j})\\}_{j\in[K],i\in[N]}$.

  2. 2.

For $j=1,\ldots,K$, obtain an approximated functional hierarchical solution
$p_{j}(x)\approx p(x,t_{j})$ by calling Algorithm 1 on the sample
$\\{X^{(i)}(t_{j})\\}_{i\in[N]}$.

  3. 3.

Output the continuous-in-time solution for $p(t,x)$, defined as follows:
Suppose that $t\in[0,T]$ satisfies $t\in(t_{j},t_{j+1})$, then $p(t,x)$ can be
approximated with $p_{j}(x)$ and $p_{j+1}(x)$ via an appropriate
interpolation.

Once the approximated solution $p$ is obtained, one can perform sampling or
compute observables based on the ansatz. As $p(t,x)$ interpolates between
$p_{j}$ and $p_{j+1}$ for $t\in(t_{j},t_{j+1})$, the task of sampling and
computing observables reduces to performing such tasks on snapshot solutions
$p_{j}$ and $p_{j+1}$. Therefore, it suffices to illustrate how to perform
such tasks on an obtained snapshot solution $p_{j}$.

#### Application I: Observable estimation

For an observable function $O\colon\mathbb{R}^{d}\to\mathbb{R}$, note that one
has $\mathbb{E}_{X\sim
p_{j}}[O(X)]=\left<O(x),p_{j}(x)\right>_{L^{2}(\mathbb{R}^{d})}$. When $O$ can
be written in terms of a finite sum of a monomial over the basis functions,
one can compute the inner product efficiently through a tensor contraction, as
the orthonormality of the function basis implies the function space inner
product will coincide with the inner product between the coefficient tensor.
In general, when $O$ is a functional hierarchical tensor, one can compute
$\left<O(x),p_{j}(x)\right>_{L^{2}(\mathbb{R}^{d})}$ through tensor diagram
contraction, as is shown in Figure 4. The tensor contraction diagram in Figure
4 admits efficient evaluation through iterative tensor contractions starting
from the middle two layers.

Figure 4: Tensor diagram representation for the measuring observable
$\mathbb{E}_{X\sim p_{j}}[O(X)]$. The top half of the tensor cores represents
the tensor network for $p_{j}$, and the bottom half represents the tensor
network for $O$.

#### Application II: Sample generation

The counting-to-sampling [35] perspective reduces the sampling task to the
previous one for observable estimation. Using conditional distributions, we
have the following decomposition of $p_{j}$ at time $t_{j}$

$p_{j}(x_{1},\ldots,x_{d})=p_{j}(x_{1})\cdot p_{j}(x_{2}|x_{1})\ldots
p_{j}(x_{d}|x_{1},\ldots,x_{d-1}).$

One then uses sequential Monte Carlo (see [28] for details) to generate i.i.d.
samples: For $k\in[d]$, suppose one has sampled the distribution up to index
$k-1$ with $x_{i}=a_{i}$ for $i=1,\ldots,k-1$. One can sample the $k$-th
variable through $X_{k}\sim p_{j}(x_{k}=x|x_{[k-1]}=a_{[k-1]})$, where
$p_{j}(x_{k}=x|x_{k-1}=a_{k-1},\ldots,x_{1}=a_{1})=\mathbb{E}_{X\sim
p_{j}}[\prod_{i\in[k-1]}\delta(X_{i}-a_{i})\delta(X_{k}-x))]$.

###  3.1 Extensions

Though we have been focused on Langevin dynamics for simplicity, our method
can be generalized to a wider range of settings. We shall remark a few notable
extensions to the current approach.

#### Fokker-Planck beyond Langevin dynamics

First, one can readily see that the approach is agnostic to the simulated SDE
structure considered. In particular, one can consider the generalized Fokker-
Planck equation,

$\partial_{t}p=\nabla\cdot(pf)+\sum_{i,j,k=1}^{d}\partial_{x_{i}x_{j}}(\sigma_{ik}\sigma_{jk}p),\quad
x\in\Omega\subset\mathbb{R}^{d},t\in[0,T],$

for a general particle system with the following form:

$dX_{t}=-f(X_{t},t)dt+\sigma(X_{t},t)dB_{t},$ (13)

where $f\colon\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d}$ and
$\sigma\colon\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d\times d}$ are general
advection and diffusion terms. While Algorithm 1 does not depend on the exact
method used to generate the sample trajectory, it might be advisable to use a
higher-order SDE simulation to guarantee sample quality in simulating (13).

#### Non-linear Fokker-Planck equations

Secondly, one can also extend our method to nonlinear Fokker-Planck equations.
Using the perspective that the Fokker-Planck equation is the Wasserstein
gradient flow on an Energy potential [36], one can consider a functional

$E(p)+\beta^{-1}\int p(x)\ln p(x)dx$

with a non-linear $E(p)$. The associated PDE is

$\partial_{t}p(t,x)=\nabla\cdot(p\nabla\frac{\delta E}{\delta
p})+\beta^{-1}\Delta p,$

where $\frac{\delta E}{\delta p}$ is the Frechet derivative of $E$ with
respect to $p$. The particle dynamics is

$dX=-(\nabla\frac{\delta E}{\delta p})(X)+\sqrt{2\beta^{-1}}dB_{t}$

with the drift term now depends on $p$. Our approach can be readily extended
to such nonlinear Fokker-Planck equations if one simulates multiple
trajectories jointly with an ensemble method (for detail, see [28]).

##  4 Application to Ginzburg-Landau model

Our main numerical treatment is dedicated to a discretized Ginzburg-Landau
model in varying physical dimensions. The high dimensionality here comes from
the discretization of an infinite-dimensional functional equation [37]. For
the sake of simplicity, we apply functional hierarchical tensor density
estimation only at time $T=1$. Generally, one runs the algorithm at multiple
time slots to solve for the solution across $[0,T]$.

#### The choice of basis function

For the Ginzburg-Landau model, the $\int_{a}|1-(x(a))^{2}|da$ term ensures
that it is highly unlikely for $\lVert x\rVert_{\infty}$ to be significantly
bigger than $1$. For the discretized model, each variable $x_{j}$ is
essentially bounded within an interval $[-B,B]$ with $B=2.5$. As a result, the
samples can be effectively seen as a distribution compactly supported in the
domain $[-B,B]^{d}$. Therefore, we can use the Fourier representation choosing
a maximal degree parameter $q$ and picking the first $n=2q+1$ Fourier basis
functions in the sine-cosine form in $[-B,B]$, denoted
$\\{\psi_{i}\\}_{i=-q}^{q}$. In effect, if one denotes $p^{\star}$ as the
target density at time $t\in[0,T]$, then our target is to approximate
$\mathcal{P}_{q}p^{\star}$, which denotes the truncated Fourier series
approximation of $p^{\star}$ restricted to degree $q$. This relationship is
expressed as follows:

$(\mathcal{P}_{q}p^{\star})(x_{1},\ldots,x_{d})=\sum_{i_{1},\ldots,i_{d}=-q}^{q}C^{\star}_{i_{1},\ldots,i_{d}}\psi_{i_{1}}(x_{1})\cdots\psi_{i_{d}}(x_{d}),$
(14)

where $C^{\star}\in\mathbb{R}^{n^{d}}$ is the coefficient tensor. When the
diffusion term in the Fokker-Planck equation has a strong influence, the
spectral approximation in (14) provides an accurate approximation with a small
$q$. The goal is then to approximate $p^{\star}$ through a compressed
hierarchical tensor $C\in\mathbb{R}^{n^{d}}$:

$p^{\star}\approx\mathcal{P}_{q}p^{\star}\approx p_{\mathrm{FHT}},\quad
p_{\mathrm{FHT}}(x_{1},\ldots,x_{d})=\sum_{i_{1},\ldots,i_{d}=-q}^{q}C_{i_{1},\ldots,i_{d}}\psi_{i_{1}}(x_{1})\cdots\psi_{i_{d}}(x_{d}).$
(15)

#### Bipartition structure

Similar to the proposed hierarchical bipartition for the 2D case, we
recursively bipartitions the Cartesian grid along its axes in an alternating
fashion. The structure for the two-dimensional case is shown in Figure 1 and
can be done likewise to other dimensions. Explicitly, assume for simplicity
that a function $u$ in $\Delta$-dimensions is discretized to the
$d=m^{\Delta}$ points denoted
$\\{u_{i_{1},\ldots,i_{\Delta}}\\}_{i_{1},\ldots,i_{\Delta}\in[m]}$ with
$m=2^{\mu}$. For each index $i_{\delta}$, one performs a length $\mu$ binary
expansion $i_{\delta}=a_{\delta 1}\ldots a_{\delta\mu}$. The variable
$u_{i_{1},\ldots,i_{\Delta}}$ is given index $k$ with the length $\mu\Delta$
binary expansion $k=a_{11}a_{21}\ldots a_{\Delta-1,\mu}a_{\Delta,\mu}$, and
the associated binary decomposition structure follows from the given index
according to (6). One can check that the construction exactly corresponds to
performing binary partition by sweeping along the $\Delta$ dimensions. For
intuition on this construction, one can check that the bipartition in Figure 1
corresponds to the indexing $k=a_{11}a_{21}a_{12}a_{22}a_{13}a_{23}$.

#### The choice of sketch function

In practice, one does not have ready access to a pre-determined sketch
function, and a higher quality in the sketch function would mean that one
would be better able to capture the spatial correlation of the density
function. As can be seen in (12), one needs to choose relatively smooth sketch
functions to ensure that an accurate approximation of terms is possible with
sample estimation. Moreover, the sketch functions should be designed to
capture important features so that terms such as $Z_{a;\bar{a}}$ could provide
meaningful linear systems to solve for $G_{q}$.

In our work, the sketch functions we choose are of two kinds. The first kind
of sketch function comprises monomials of the Fourier basis of low degree. For
$a=I_{k}^{(l)}$ or $a=[d]-I_{k}^{(l)}$ (see definition in (6)), the chosen
sketch functions are of the following form:

$f(x_{a})=\prod_{j\in a}{\psi_{i_{j}}(x_{j})},$

for some choice of $-q\leq i_{j}\leq q$ for each $j$. The total Fourier degree
can be characterized by $\mathrm{deg}(f)=\sum_{j\in a}|i_{j}|$, and one
typically choose $f$ with small $\mathrm{deg}(f)$.

The second kind of sketch function is motivated by the renormalization group,
in particular by the coarse-graining of the variables. In particular, one
chooses a cluster $h\subset a$ and a specific Fourier mode index $i$. The
resultant function $g_{h,i}$ is an averaging of the $i$-th Fourier mode over
variables in $h$:

$g_{h,i}(x_{a})\equiv\frac{1}{|h|}\sum_{j\in h}{\psi_{i}(x_{j})}.$

In our particular implementation, for $a=I_{k}^{(l)},$ a typical example of
chosen sketch function would be $h=I_{2k}^{(l+1)}$ and $i=1$. Then, the
resultant sketch function would be

$g_{h,1}(x_{a})\equiv\frac{1}{|h|}\sum_{j\in h}{\sin(\pi x_{j}/B)},$

which can be thought of as coarse-grained information on half of the variables
in $a$. One can similarly define $g_{a-h,1}(x_{a})$. One can augment the
sketch function set by including terms such as $g_{h,1}^{2}$ and
$g_{a-h,1}^{2}$. Empirically, including such coarse-grained low-order Fourier
modes greatly increases model performance and stability.

###  4.1 1D Ginzburg-Landau potential

The first numerical example is a 1D Ginzburg-Landau model. The potential
energy is defined as

$V(x_{1},\ldots,x_{m}):=\frac{\lambda}{2}\sum_{i=1}^{m+1}\left(\frac{x_{i}-x_{i-1}}{h}\right)^{2}+\frac{1}{4\lambda}\sum_{i=1}^{m}\left(1-x_{i}^{2}\right)^{2},$
(16)

where $h=\frac{1}{m+1}$ and $x_{0}=x_{m+1}=0$. In particular, we fix $m=256$,
$\lambda=0.01$ and $\beta=\frac{1}{8}$. The initial condition is
$p_{0}(x)=\delta(x-(0,\ldots,0))$ and the dimension $d$ is $256$.

Figure 5: 1D Ginzburg-Landau model. The plots of the marginal distribution at
$(x_{19},x_{59})$.

We perform $N=6000$ SDE simulations with $T=1$ and $\delta t=\frac{1}{2000}$,
all starting from the initial location $\left(0,\ldots,0\right)$. The maximal
Fourier degree $q$ is taken to be $q=15$ to account for the localization of
the samples, and the maximal internal bond dimension is $r=6$.

Figure 5 compares the marginal distribution of $(x_{15},x_{59})$ between
empirical distributions and the model obtained by hierarchical tensor
sketching. One can see that hierarchical tensor sketching is successful at
capturing the correlation between two variables at relatively faraway points.

###  4.2 2D Ginzburg-Landau potential

In the second numerical example, we consider the usual 2D Ginzburg-Landau
model. The potential energy is defined as

$V(x_{(1,1)},\ldots,x_{(m,m)}):=\frac{\lambda}{2}\sum_{v\sim
w}\left(\frac{x_{v}-x_{w}}{h}\right)^{2}+\frac{1}{4\lambda}\sum_{v}\left(1-x_{v}^{2}\right)^{2},$
(17)

where $h=\frac{1}{m+1}$ and $x_{(l,0)}=x_{(l,m+1)}=x_{(0,l)}=x_{(m+1,l)}=0$
for $l=1,\ldots,m$. The parameters are $m=16$, $\lambda=0.03$,
$\beta=\frac{1}{5}$, and the initial condition is
$p_{0}(x)=\delta(x-(0,\ldots,0))$. The dimension $d$ is $256$.

Figure 6: 2D Ginzburg-Landau model. The plots of marginal distribution at
$(x_{(3,5)},x_{(6,8)})$. Figure 7: Two-point correlation function
$f(i,j)=\mathrm{Corr}\left(x_{(8,8)},x_{(i,j)}\right)$ of the sketched
functional hierarchical tensor, compared with the ground truth from a Monte-
Carlo estimation from 60,000 samples.

We perform $N=6000$ SDE simulation with $T=1$ and $\delta t=0.003$, starting
from the initial condition $\left(0,\ldots,0\right)$. We take $q=10$ and a
maximal internal bond of $r=20$. To ensure numerical stability, the internal
bond dimension is dynamically chosen according to the singular values
information obtained during sketching.

Figure 6 shows that the obtained marginal distribution of
$(x_{(3,5)},x_{(6,8)})$ closely matches the empirical distributions. To test
the accuracy of observable estimation, Figure 7 plots the predicted two-point
correlation function

$f(i,j):=\mathrm{Corr}\left(x_{(8,8)},x_{(i,j)}\right)\equiv\frac{\mathrm{Cov}\left[x_{(8,8)},x_{(i,j)}\right]}{\sigma_{x_{(8,8)}}\sigma_{x_{(i,j)}}},$

where the result closely matches the ground truth obtained from a Monte-Carlo
estimation from 60,000 samples. The accuracy is quite remarkable as the
internal bond $r$ is kept at a relatively small level.

###  4.3 3D Ginzburg-Landau potential

Here, we consider a 3D Ginzburg-Landau model. Similarly, define a 3D Cartesian
grid $D=(Q,E)$, where $Q:=\left\\{(i,j,k)\mid i,j,k=0,\ldots,m+1\right\\}$.
The potential energy is defined as

$V(x_{(1,1,1)},\ldots,x_{(m,m,m)}):=\frac{\lambda}{2}\sum_{v\sim
w}\left(\frac{x_{v}-x_{w}}{h}\right)^{2}+\frac{1}{4\lambda}\sum_{v}\left(1-x_{v}^{2}\right)^{2},$
(18)

where $h=\frac{1}{m+1}$, and the boundary condition is such that
$x_{(i,j,k)}=0$ if one of $i,j,k$ equals to $0$ or $m+1$. For parameters, we
fix $m=8$, $\lambda=0.01$ and $\beta=\frac{1}{10}$. The initial condition is
$p_{0}(x)=\delta(x-(0,\ldots,0))$ and the dimension $d$ is $512$.

Figure 8: 3D Ginzburg-Landau model. The plots of marginal distribution at
$(x_{(1,1,1)},x_{(3,3,3)})$.

We perform $N=6000$ SDE simulation with $T=1$ and $\delta t=0.002$, starting
from the initial condition $\left(0,\ldots,0\right)$. We take $q=10$ and a
maximal internal bond of $r=20$, where the internal bond is likewise chosen
dynamically. Figure 8 shows the performance on the marginal distribution of
$(x_{(1,1,1)},x_{(3,3,3)})$ between empirical distributions and the model
obtained by hierarchical tensor sketching, and similarly the result shows
quite good performance.

##  5 Conclusion

We introduce a novel density estimation approach that combines particle
methods with a functional hierarchical tensor network in solving the high-
dimensional Fokker-Planck equations. The algorithm is applied to the
discretized Ginzburg-Landau model in 1D, 2D, and 3D. This method points to a
new direction for tackling general high-dimensional Fokker-Planck equations.
Combined with sensible numerical treatment, this approach has the potential
for modeling high-dimensional particle dynamics with high fidelity. An open
question is whether one can extend a tensor network approach to solving
similar equations, such as the forward Kolmogorov equation and the Hamilton-
Jacobi equation.

#### Acknowledgement

We thank Yuehaw Khoo for constructive discussions.

## References

  * Han et al. [2018] J. Han, A. Jentzen, W. E,  Solving high-dimensional partial differential equations using deep learning,  Proceedings of the National Academy of Sciences 115 (2018) 8505–8510.
  * Ginzburg et al. [2009] V. L. Ginzburg, V. L. Ginzburg, L. Landau, On the theory of superconductivity, Springer, 2009.
  * Hoffmann and Tang [2012] K.-H. Hoffmann, Q. Tang, Ginzburg-Landau phase transition theory and superconductivity, volume 134, Birkhäuser, 2012.
  * Hohenberg and Krekhov [2015] P. C. Hohenberg, A. P. Krekhov,  An introduction to the ginzburg–landau theory of phase transitions and nonequilibrium patterns,  Physics Reports 572 (2015) 1–42.
  * E et al. [2004] W. E, W. Ren, E. Vanden-Eijnden,  Minimum action method for the study of rare events,  Communications on pure and applied mathematics 57 (2004) 637–656.
  * Hastie et al. [2009] T. Hastie, R. Tibshirani, J. H. Friedman, J. H. Friedman, The elements of statistical learning: data mining, inference, and prediction, volume 2, Springer, 2009.
  * Chen and Khoo [2023] Y. Chen, Y. Khoo,  Combining particle and tensor-network methods for partial differential equations via sketching,  arXiv preprint arXiv:2305.17884 (2023).
  * Ren et al. [2023] Y. Ren, H. Zhao, Y. Khoo, L. Ying,  High-dimensional density estimation with tensorizing flow,  Research in the Mathematical Sciences 10 (2023) 30.
  * Dektor and Venturi [2021] A. Dektor, D. Venturi,  Dynamic tensor approximation of high-dimensional nonlinear pdes,  Journal of Computational Physics 437 (2021) 110295.
  * Dektor et al. [2021] A. Dektor, A. Rodgers, D. Venturi,  Rank-adaptive tensor methods for high-dimensional nonlinear pdes,  Journal of Scientific Computing 88 (2021) 36.
  * Soley et al. [2021] M. B. Soley, P. Bergold, A. A. Gorodetsky, V. S. Batista,  Functional tensor-train chebyshev method for multidimensional quantum dynamics simulations,  Journal of Chemical Theory and Computation 18 (2021) 25–36.
  * Bigoni et al. [2016] D. Bigoni, A. P. Engsig-Karup, Y. M. Marzouk,  Spectral tensor-train decomposition,  SIAM Journal on Scientific Computing 38 (2016) A2405–A2439.
  * Gorodetsky et al. [2019] A. Gorodetsky, S. Karaman, Y. Marzouk,  A continuous analogue of the tensor-train decomposition,  Computer methods in applied mechanics and engineering 347 (2019) 59–84.
  * Hur et al. [2023] Y. Hur, J. G. Hoskins, M. Lindsey, E. M. Stoudenmire, Y. Khoo,  Generative modeling via tensor train sketching,  Applied and Computational Harmonic Analysis 67 (2023) 101575.
  * Östlund and Rommer [1995] S. Östlund, S. Rommer,  Thermodynamic limit of density matrix renormalization,  Physical review letters 75 (1995) 3537.
  * Cirac et al. [2021] J. I. Cirac, D. Perez-Garcia, N. Schuch, F. Verstraete,  Matrix product states and projected entangled pair states: Concepts, symmetries, theorems,  Reviews of Modern Physics 93 (2021) 045003.
  * Vidal [2007] G. Vidal,  Entanglement renormalization,  Physical review letters 99 (2007) 220405.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational physics 378 (2019) 686–707.
  * Sirignano and Spiliopoulos [2018] J. Sirignano, K. Spiliopoulos,  Dgm: A deep learning algorithm for solving partial differential equations,  Journal of computational physics 375 (2018) 1339–1364.
  * Yu et al. [2018] B. Yu, et al.,  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems,  Communications in Mathematics and Statistics 6 (2018) 1–12.
  * Kovachki et al. [2021] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar,  Neural operator: Learning maps between function spaces,  arXiv preprint arXiv:2108.08481 (2021).
  * Tabak and Vanden-Eijnden [2010] E. G. Tabak, E. Vanden-Eijnden,  Density estimation by dual ascent of the log-likelihood,  Communications in Mathematical Sciences 8 (2010) 217–233.
  * Rezende and Mohamed [2015] D. Rezende, S. Mohamed,  Variational inference with normalizing flows,  in: International conference on machine learning, PMLR, 2015, pp. 1530–1538.
  * Hinton [2002] G. E. Hinton,  Training products of experts by minimizing contrastive divergence,  Neural computation 14 (2002) 1771–1800.
  * LeCun et al. [2006] Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, F. Huang,  A tutorial on energy-based learning,  Predicting structured data 1 (2006).
  * Song and Ermon [2019] Y. Song, S. Ermon,  Generative modeling by estimating gradients of the data distribution,  Advances in Neural Information Processing Systems 32 (2019).
  * Song et al. [2021] Y. Song, C. Durkan, I. Murray, S. Ermon,  Maximum likelihood training of score-based diffusion models,  Advances in Neural Information Processing Systems 34 (2021) 1415–1428.
  * Liu [2001] J. Liu, Monte Carlo strategies in scientific computing, volume 75, Springer, 2001.
  * Chen et al. [2023] Y. Chen, J. Hoskins, Y. Khoo, M. Lindsey,  Committor functions via tensor networks,  Journal of Computational Physics 472 (2023) 111646.
  * Khoromskij and Schwab [2011] B. N. Khoromskij, C. Schwab,  Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes,  SIAM journal on scientific computing 33 (2011) 364–385.
  * Eigel et al. [2019] M. Eigel, J. Neumann, R. Schneider, S. Wolf,  Non-intrusive tensor reconstruction for high-dimensional random pdes,  Computational Methods in Applied Mathematics 19 (2019) 39–53.
  * Bachmayr et al. [2016] M. Bachmayr, R. Schneider, A. Uschmajew,  Tensor networks and hierarchical tensors for the solution of high-dimensional partial differential equations,  Foundations of Computational Mathematics 16 (2016) 1423–1472.
  * Richter et al. [2021] L. Richter, L. Sallandt, N. Nüsken,  Solving high-dimensional parabolic pdes using the tensor train format,  in: International Conference on Machine Learning, PMLR, 2021, pp. 8998–9009.
  * Peng et al. [2023] Y. Peng, Y. Chen, E. M. Stoudenmire, Y. Khoo,  Generative modeling via hierarchical tensor sketching,  arXiv preprint arXiv:2304.05305 (2023).
  * Jerrum [2003] M. Jerrum, Counting, sampling and integrating: algorithms and complexity, Springer Science & Business Media, 2003.
  * Jordan et al. [1998] R. Jordan, D. Kinderlehrer, F. Otto,  The variational formulation of the fokker–planck equation,  SIAM journal on mathematical analysis 29 (1998) 1–17.
  * Rosen [1971] G. Rosen,  Functional calculus theory for incompressible fluid turbulence,  Journal of Mathematical Physics 12 (1971) 812–820.


# Multivariate polynomial splines on generalized oranges

Maritza Sirvent<EMAIL_ADDRESS>Tatyana Sorokina<EMAIL_ADDRESS>Nelly Villamizar<EMAIL_ADDRESS>Beihui Yuan
<EMAIL_ADDRESS>
###### Abstract

We consider spaces of multivariate splines defined on a particular type of
simplicial partitions that we call (generalized) oranges. Such partitions are
composed of a finite number of maximal faces with exactly one shared medial
face. We reduce the problem of finding the dimension of splines on oranges to
computing dimensions of splines on simpler, lower-dimensional partitions that
we call projected oranges. We use both algebraic and Bernstein-Bézier tools.

###### keywords:

Multivariate spline functions , dimension of spline spaces , Bernstein-Bézier
methods, cofactor criterion.

###### MSC:

[2020] 41A15 , 13D40 , 65D07.

[ohio]organization=Department of Mathematics, The Ohio State
University,addressline=231 West 18th Avenue, city=Columbus, postcode=OH 43210,
state=Ohio, country=United States

[towson]organization=Department of Mathematics, Towson
University,addressline=8000 York Road, city=Towson, postcode=MD 21252,
state=Maryland, country=United States

[swansea]organization=Department of Mathematics, Swansea
University,addressline=Fabian Way, city=Swansea, postcode=SA1 8EN,
country=United Kingdom

##  1 Introduction and preliminaries

Let $\Delta$ be a simplicial partition of an underlying domain $\Omega$ in
$\mathbb{R}^{k}$ formed as a geometric realization of a $k$-dimensional
simplicial complex. For integers $0\leqslant r\leqslant d$, the space of
splines $S^{r}_{d}(\Delta)$ is defined as the set of $C^{r}$-smooth piecewise
polynomial functions of degree at most $d$ on the partition $\Delta$. The
space $S^{r}_{d}(\Delta)$ is a real vector space. These spline spaces are used
throughout numerical analysis and approximation theory to solve diverse
problems such as the interpolation and approximation of data, design of curves
and surfaces, and, in the finite element method, for the solution of
differential equations, among others. A fundamental question is to determine
the dimension of such spaces. The problem has proven to be difficult due to
the dependency of the dimension on the specific geometry of $\Delta$.

In 1973, Strang [12], published a famous conjecture on the generic dimension
of $S^{1}_{d}(\Delta)$ for the bivariate case. A simplicial partition $\Delta$
is said to be generic provided that for all sufficiently small perturbations
of the location of the vertices of $\Delta$, the resulting decomposition
$\tilde{\Delta}$ satisfies $\dim S^{r}_{d}(\tilde{\Delta})=\dim
S^{r}_{d}(\Delta)$, see [6] where this definition was explicitly coined.
Although the term generic dimension was not explicitly mentioned by Strang, he
was aware of the fact that the dimension might increase for some particular
configurations. In [3], Alfeld and Schumaker prove a dimension formula for
$S^{r}_{d}(\Delta)$, when the polynomial degree $d\geqslant 4r+1$, and
$\Delta$ is a general triangulation. This result was later extended to
$d\geqslant 3r+2$ by Hong in [9], and to $d\geqslant 3r+1$ for generic
triangulations by Alfeld and Schumaker in [4]. For $r<d<3r+2$, the only
dimension known explicitly is the case when $r=1$ and $d=4$, see [2].

In 1988, Billera, see [7], pioneered the use of algebraic homological methods
to study multivariate splines, and proved the conjecture given by Strang.

For the trivariate case, Alfeld, Schumaker, and Whiteley, see [6], use the
Bernstein-Bézier analysis and algebraic homology and rigidity theory to give
an explicit expression for the generic dimension of $S^{1}_{d}(\Delta)$, for
$d\geqslant 8$.

For the spatial dimension $k\geqslant 3$, no explicit dimension formula is
known for general triangulations. There are results on upper and lower bounds
and exact dimensions for some specific triangulations. For example, in [8],
DiPasquale and Villamizar proved a lower bound for the dimension for splines
on tetrahedral stars of a vertex. A natural generalization of a star of a
vertex is a star of a simplex. A star of a simplex $\tau\in\Delta$ is the set
of simplices in $\Delta$ that contain $\tau$. If a simplicial complex $\Delta$
equals the closure of the star of $\tau$, we just call $\Delta$ a star of
$\tau$. In this paper, we study stars of a simplex of arbitrary dimension. We
call such stars _generalized oranges_ or _oranges_.

Recall that a $k$-dimensional simplicial complex is _pure_ is all its maximal
faces are of dimension $k$.

######  Definition 1.1 (Generalized orange).

For integers $0\leqslant i\leqslant k$, a _$(k,i)$ -orange_ is a pure
$k$-dimensional simplicial complex $\mathcal{O}$ composed by $n$ maximal faces
and with exactly one face $\tau$ of dimension $k-i$ such that every maximal
face of $\mathcal{O}$ contains $\tau$. We say that $\tau$ is a _medial
simplex_.

Note that the two extremal cases are $i=0$ and $i=k$. For $i=0$, $\Delta=\tau$
is a simplex. For $i=k$, the medial face $\tau$ is a vertex, and $\Delta$ is a
star of a vertex. Examples of $(3,2),(2,1)$, and $(3,1)$-oranges are shown in
Figs. 1, 2 and 3, respectively.

Figure 1: Example of a $(3,2)$-orange, it is a 3-dimensional simplicial
complex with medial simplex $\tau$, which is the $1$-face common to all
tetrahedra.

Historically, the term “orange” was introduced in [5] to describe a
tetrahedral partition with exactly one interior edge which is common to all
tetrahedra in the partition, see Fig. 1. This matches our definition of a
$(3,2)$-orange except that the medial edge does not have to be interior. The
dimension for the spline spaces on a $(3,2)$-orange with interior medial edge
was found in [5] using the trivariate cofactor method, see [10] for the proof
as well. In [11], it was noted that the dimension for the spline spaces on a
$(3,2)$-orange with interior medial edge can be computed by adding up the
dimensions of bivariate splines on a planar vertex star obtained by projecting
the orange along the common interior edge.

In this article we generalize the idea of the projection to an $i$-dimensional
star of a vertex, and relate the dimension of the spline spaces over a
$(k,i)$-orange $\mathcal{O}$ to those over a projected orange $\mathcal{C}$,
see Section 2 for the exact definition, by proving the following result.

######  Theorem 1.1.

Let $\mathcal{C}\subseteq\mathbb{R}^{i}$ be the projected orange of a
$(k,i)$-orange $\mathcal{O}\subseteq\mathbb{R}^{k}$. Then

$\dim S^{r}_{d}(\mathcal{O})=\sum_{j=0}^{d}\binom{d+k-i-j-1}{k-i-1}\dim
S^{r}_{i}(\mathcal{C}).$ (1)

For the two extremal cases, $i=0$ and $i=k$, Theorem 1.1 trivially holds. In
the case $i=0$, the orange $\mathcal{O}=\tau$ is a simplex, $\mathcal{C}$ is a
point, and

$\dim
S^{r}_{d}(\mathcal{O})=\sum_{j=0}^{d}\binom{d+k-j-1}{k-1}=\binom{d+k}{k}.$

For $i=k$, we know that $\mathcal{O}=\mathcal{C}$. In this case, formula (1)
is just the trivial identity $\dim S^{r}_{d}(\mathcal{O})=\dim
S^{r}_{d}(\mathcal{C})$.

The paper is organized as follows. In Section 2, we construct a projection
$\mathcal{C}$ of a $(k,i)$-orange $\mathcal{O}$ along the medial face. We
prove that the projection $\mathcal{C}$ is also a simplicial complex. In
Section 3, we analyze the relation between the space of splines on
$\mathcal{O}$ and that on $\mathcal{C}$ by using the so-called cofactor
criterion. Finally, we prove our main result, Theorem 1.1, by two methods: the
algebraic method in Section 4, and the Bernstein-Bézier method in Section 5.
The algebraic method is based on Hilbert series of spline spaces viewed as
filtered vector spaces, while the Bernstein-Bézier method uses minimal
determining sets and smoothness conditions in Bernstein-Bézier form.

##  2 Geometry

We embed a $(k,i)$-orange $\mathcal{O}$ in $\mathbb{R}^{k}$, and we assume
that the medial face $\tau$ in $\mathcal{O}$ of dimension $(k-i)$ is in the
$(k-i)$-dimensional coordinate subspace of $\mathbb{R}^{k}$ defined by
$x_{1}=\cdots=x_{i}=0$. As a simplifying assumption, we identify $\mathcal{O}$
with its embedding in $\mathbb{R}^{k}$. For this choice of coordinates, we
define the following projection:

$\begin{split}\pi\colon\mathbb{R}^{k}&\to\mathbb{R}^{i},\\\
(x_{1},\dots,x_{k})&\mapsto(x_{1},\dots,x_{i}).\end{split}$ (2)

Additionally, we write $\mathcal{C}$ for the collection of sets
$\pi(\sigma)\subseteq\mathbb{R}^{i}$ of images under the projection $\pi$ of
the faces $\sigma\in\mathcal{O}$. Namely,

${\mathcal{C}}:=\pi(\mathcal{O})=\bigl{\\{}\pi(\sigma)\colon\sigma\in\mathcal{O}\bigr{\\}}.$
(3)

We shall refer to ${\mathcal{C}}$ as a _projected orange_. Note that the
projection of the medial face $\pi(\tau)$ is the origin $O$ of
$\mathbb{R}^{i}$. For any $X,Y\subseteq\mathbb{R}^{k}$, denote by $X*Y$ the
_join_ of $X$ and $Y$ which is the union of all line segments joining the
points in $X$ to the points in $Y$. Note that

$\pi(X*Y)=\pi(X)*\pi(Y).$ (4)

If $X=\\{v\\}$ and $Y=\\{w\\}$, we simply write $v*w$ instead of $X*Y$. Recall
that the set of points $\\{p_{1},\dots,p_{j}\\}\subseteq\mathbb{R}^{k}$ is
said to be in _general position_ if its span is a $(j-1)$-subspace of
$\mathbb{R}^{k}$. Note that every $m$-simplex $\sigma\in\mathcal{O}$ can be
written as the join of $m+1$ points in general position.

######  Proposition 2.1.

If $\mathcal{O}$ is a $(k,i)$-orange, then the projected orange $\mathcal{C}$
defined in (3) is a simplicial complex. Moreover, $\mathcal{C}$ is a star of a
vertex.

Before we prove Proposition 2.1, we need some preliminary lemmas that will
help us verify that every element in $\mathcal{C}$ is a simplex, every face of
a simplex in $\mathcal{C}$ is also in $\mathcal{C}$, and that the intersection
of any two elements in $\mathcal{C}$ is a face of each of them.

######  Lemma 2.2.

Let $\mathcal{O}$ be a $(k,i)$-orange with the medial simplex $\tau$. The
image $\pi(\sigma)$ of a simplex $\sigma\in\mathcal{O}$ is a simplex.
Moreover, if $\sigma$ and $\tau$ do not intersect, then $\pi(\sigma)$ and
$\sigma$ have the same dimension. Otherwise,
$\dim(\pi(\sigma))=\dim(\sigma)-\dim(\sigma\cap\tau)$.

###### Proof.

For a simplex $\sigma\in\mathcal{O}$, the vertices of $\sigma$ can always be
divided into two sets, one of which contains all vertices that are in $\tau$,
the other contains those not in $\tau$. This means $\sigma$ can be written in
the form

$\sigma=\alpha*\beta,$ (5)

where $\alpha,\beta$ are faces of $\sigma$ such that
$\alpha\cap\tau=\emptyset$ and $\beta=\sigma\cap\tau$. Then by (4), we have

$\pi(\sigma)=\begin{cases}\pi(\alpha),&\mbox{if}~{}\beta=\emptyset,\\\
\pi(\alpha)*O,&\mbox{if}~{}\beta\neq\emptyset.\end{cases}$ (6)

Let $z_{i+1},\dots,z_{k+1}$ be the vertices of $\tau$. Without loss of
generality assume that $z_{k+1}=O$. Let $p_{1},\dots,p_{j}$ be the vertices of
$\alpha$, where $j\leqslant i$. Then $p_{1},\dots,p_{j},z_{i+1},\dots,z_{k+1}$
are in general position and $p_{1},\dots,p_{j},z_{i+1},\dots,z_{k}$ are
linearly independent in $\mathbb{R}^{k}$. We claim that then
$\bigl{\\{}\pi(p_{1}),\dots,\pi(p_{j})\bigr{\\}}$ is also linearly
independent. Note that in general the projection of a linearly independent set
of vectors is not necessarily linearly independent. Let
$V=\mathrm{span}\\{p_{1},\dots,p_{j},z_{i+1},\dots,z_{k}\\}$, and let
$U=\mathrm{span}\\{z_{i+1},\dots,z_{k}\\}$. Then $\dim V=k+j-i$, and $\dim
U=k-i$. Since $U$ is the $(k-i)$-subspace containing $\tau$, we have
$\mathrm{ker}(\pi)=U$, and $\dim(\pi(V))=j$. Therefore, the set
$\bigl{\\{}\pi(p_{1}),\dots,\pi(p_{j})\bigr{\\}}$ spans $\pi(V)$, and it is
linearly independent.

This shows that both $\pi(\alpha)$ and $\pi(\alpha)*O$ are simplices.
Moreover, if $\beta=\emptyset$, then $\dim\pi(\sigma)=\dim\pi(\alpha)$.
Otherwise, $\dim(\pi(\sigma))=\dim(\sigma)-\dim(\beta)$, and the proof is
complete. ∎

######  Corollary 2.3.

The projected orange $\mathcal{C}$ is a collection of simplices. A face of an
element in $\mathcal{C}$ is also in $\mathcal{C}$.

###### Proof.

It immediately follows from Lemma 2.2 that every element of $\mathcal{C}$ is a
simplex. Now if a simplex $\omega\in\mathcal{C}$, then by Definition 3, there
is $\sigma\in\mathcal{O}$ such that $\pi(\sigma)=\omega$. We may assume
$\sigma=v_{1}*\dots*v_{j}$, for vertices $v_{s}\in\mathbb{R}^{k}$ and
$j\leqslant k$. Then
$\pi(v_{1}*\dots*v_{j})=\pi(v_{1})*\cdots*\pi(v_{j})=\omega$. Thus any face of
$\omega$ can be written as
$\pi(v_{j_{1}})*\cdots*\pi(v_{j_{\ell}})=\pi(v_{j_{1}}*\cdots*v_{j_{\ell}})$
and, hence, by Definition 3, this face is in $\mathcal{C}$. ∎

######  Lemma 2.4.

Let $v_{1},v_{2}$ be two vertices in $\mathcal{O}$. If neither $v_{1}$ nor
$v_{2}$ belongs to $\tau$, then $\pi(v_{1})=\pi(v_{2})$ implies $v_{1}=v_{2}$.

###### Proof.

Let $\alpha_{1}=v_{1}*\tau$ and $\alpha_{2}=v_{2}*\tau$. Both $\alpha_{1}$ and
$\alpha_{2}$ are $(k-i+1)$-faces in $\mathcal{O}$ and
$\tau\subseteq\alpha_{1}\cap\alpha_{2}$. Since $\mathcal{O}$ is a simplicial
complex and $\dim\tau=k-i$, then there are only two possibilities:

$\tau=\alpha_{1}\cap\alpha_{2}\quad\text{or}\quad\alpha_{1}=\alpha_{2}.$ (7)

Let $z_{i+1},\dots,z_{k+1}$ be the vertices of $\tau$, where we again assume
that $z_{k+1}$ is the origin. Since $\pi(v_{1})=\pi(v_{2})$, we note that
$v_{1}-v_{2}$ is in the $(k-i)$-dimensional coordinate space spanned by
$\tau$. Thus

$v_{1}-v_{2}=\beta_{1}z_{i+1}+\dots+\beta_{k-i}z_{k}.$ (8)

Without loss of generality, we assume that the first $j$ coefficients
$\bigl{\\{}\beta_{\ell}\bigr{\\}}_{\ell=1}^{j}$ are negative, while the
remaining $k-i-j$ coefficients
$\bigr{\\{}\beta_{\ell}\bigr{\\}}_{\ell=j+1}^{k-i}$ are non-negative. We next
rewrite (8) as follows:

$v_{1}-\beta_{1}z_{i+1}-\dots-\beta_{j}z_{i+j}=v_{2}+\beta_{j+1}z_{i+j+1}+\dots+\beta_{k-i}z_{k}=v.$
(9)

Let $m=\left(1+\sum_{\ell=1}^{k-i}|\beta_{\ell}|\right)^{-1}$, and let $u=mv$.

Since $u$ can be written as two different convex combinations, one in the
vertices of $\alpha_{1}$ and a second in the vertices of $\alpha_{2}$, then
$u\in\alpha_{1}\cap\alpha_{2}$. From (9), it follows that $u\notin\tau$, so by
(7) we must have $\alpha_{1}=\alpha_{2}$, and hence $v_{1}=v_{2}$. ∎

######  Corollary 2.5.

Let $\omega\in\mathcal{C}$ such that the origin $O$ is not contained in
$\omega$. Then there exists a unique $\alpha\in\mathcal{O}$ such that
$\pi(\alpha)=\omega$.

###### Proof.

The existence is immediate from the definition of $\mathcal{C}$. Assume there
are two simplices $\alpha$ and $\alpha^{\prime}$ in $\mathcal{O}$ such that
$\pi(\alpha)=\pi(\alpha^{\prime})=\omega$. Since $O\notin\omega$, we have
$\alpha\cap\tau=\emptyset$, and $\alpha^{\prime}\cap\tau=\emptyset$. Then by
Lemma 2.2, $\dim(\alpha)=\dim(\alpha^{\prime})=\dim(\omega)=:j$. We may assume
$\alpha=p_{1}*\dots*p_{j}$ and
$\alpha^{\prime}=p^{\prime}_{1}*\dots*p^{\prime}_{j}$ for vertices $p_{s}$’s
and $p^{\prime}_{t}$’s in $\mathcal{O}$. Hence
$\pi(p_{1})*\dots*\pi(p_{j})=\pi(p^{\prime}_{1})*\dots*\pi(p^{\prime}_{j})=\omega$.
Therefore, possibly after a permutation,
$\pi(p_{\ell})=\pi(p^{\prime}_{\ell})$ for all $\ell=1,\dots,j$. By Lemma 2.4,
this implies $p_{\ell}=p^{\prime}_{\ell}$ for all $\ell=1,\dots,j$, and
$\alpha=\alpha^{\prime}$. ∎

######  Corollary 2.6.

Let $\omega\in\mathcal{C}$, and let the origin $O$ be contained in $\omega$.
Then there is a unique face $\sigma\in\mathcal{O}$ such that
$\pi(\sigma)=\omega$ and $\tau\subseteq\sigma$.

###### Proof.

By Corollary 2.3, every face $\omega$ of $\mathcal{C}$ is a simplex, so we may
assume $\omega=\theta*O$, where $O\not\in\theta$ and $\theta$ is also a face
of $\mathcal{C}$. By Corollary 2.5, there is a face $\alpha\in\mathcal{O}$
such that $\pi(\alpha)=\theta$. Hence, $\pi(\alpha*\tau)=\omega$. So we can
take $\sigma=\alpha*\tau$. To prove the uniqueness, assume there is a face
$\sigma^{\prime}$ such that $\pi(\sigma^{\prime})=\omega$ and
$\tau\subseteq\sigma^{\prime}$. Then we can write
$\sigma^{\prime}=\alpha^{\prime}*\tau$ for some
$\alpha^{\prime}\cap\tau=\emptyset$. Then $\pi(\alpha^{\prime})=\theta$.
Corollary 2.5 implies $\alpha=\alpha^{\prime}$. Hence,
$\sigma=\sigma^{\prime}$. ∎

Next corollary is an immediate consequence of Lemma 2.2 and Corollary 2.6. We
are going to use it in Section 3.

######  Corollary 2.7.

There is a one-to-one correspondence between the maximal faces of an orange
$\mathcal{O}$ and its projection $\mathcal{C}$.

###### Proof.

If $\sigma\in\mathcal{O}$ is a maximal face then the medial simplex
$\tau\subseteq\sigma$, and so Lemma 2.2 implies that $\pi(\sigma)$ is a
simplex of dimension $i$. The existence and uniqueness of
$\sigma\in\mathcal{O}$ for every maximal face of $\mathcal{C}$ follow from
Corollary 2.6. ∎

######  Lemma 2.8.

For any two faces $\sigma_{1},\sigma_{2}\in\mathcal{O}$ containing $\tau$, we
have

$\pi(\sigma_{1})\cap\pi(\sigma_{2})=\pi(\sigma_{1}\cap\sigma_{2}).$ (10)

In particular, (10) implies that
$\pi(\sigma_{1})\cap\pi(\sigma_{2})\in\mathcal{C}$ in this case.

###### Proof.

We only need to prove
$\pi(\sigma_{1})\cap\pi(\sigma_{2})\subseteq\pi(\sigma_{1}\cap\sigma_{2})$.
Let $b\in\pi(\sigma_{1})\cap\pi(\sigma_{2})$, and let $b\neq O$. Then there
exist $a_{1}\in\sigma_{1}$ and $a_{2}\in\sigma_{2}$ such that
$b=\pi(a_{1})=\pi(a_{2})$. If either $a_{1}\in\sigma_{1}\cap\sigma_{2}$ or
$a_{2}\in\sigma_{1}\cap\sigma_{2}$ then $b\in\pi(\sigma_{1}\cap\sigma_{2})$,
and we are done. Suppose that neither $a_{1}$ nor $a_{2}$ is in
$\sigma_{1}\cap\sigma_{2}$. Then $(a_{1}*\tau)\cap\sigma_{2}=\tau$ and
$(a_{2}*\tau)\cap\sigma_{1}=\tau$. Since
$(a_{1}*\tau)\cap(a_{2}*\tau)\subseteq\sigma_{1}\cap\sigma_{2}$, we have

$(a_{1}*\tau)\cap(a_{2}*\tau)=(a_{1}*\tau)\cap(a_{2}*\tau)\cap\sigma_{1}=\tau.$

Without loss of generality, we may assume that $x_{1}(b)=1$ and
$x_{2}(b)=\dots=x_{i}(b)=0$. Note that $a_{1}*\tau$ and $a_{2}*\tau$ are in
the same $(k-i+1)$-subspace $V\subseteq\mathbb{R}^{k}$ defined by
$x_{2}=\dots=x_{i}=0$. By a similar reasoning as in the proof of Lemma 2.4, we
get $\tau\neq(a_{1}*\tau)\cap(a_{2}*\tau)$, which is a contradiction.
Therefore, the hypothesis that neither $a_{1}$ nor $a_{2}$ is in
$\sigma_{1}\cap\sigma_{2}$ is false. ∎

######  Corollary 2.9.

For any two faces $\alpha,\sigma\in\mathcal{O}$ such that
$\alpha\cap\tau=\emptyset$, we have

$\pi(\alpha)\cap\pi(\sigma)=\pi(\alpha\cap\sigma).$

In particular, $\pi(\alpha)\cap\pi(\sigma)\in\mathcal{C}$ in this case.

###### Proof.

We only need to show that
$\pi(\alpha)\cap\pi(\sigma)\subseteq\pi(\alpha\cap\sigma)$. Note that
$\sigma=\alpha^{\prime}*\beta^{\prime}$ such that
$\alpha^{\prime}\cap\tau=\emptyset$ and $\beta^{\prime}\subseteq\tau$. Since
$\mathcal{O}$ is a simplicial complex, then
$\alpha\cap\sigma=\alpha\cap\alpha^{\prime}$, and
$(\alpha*\tau)\cap(\alpha^{\prime}*\tau)=(\alpha\cap\alpha)*\tau$. By Lemma
2.8, we have that
$\pi(\alpha\cap\alpha^{\prime})*\pi(\tau)=(\pi(\alpha)*\pi(\tau))\cap(\pi(\alpha^{\prime})*\pi(\tau))$.
This means that for any $b\in\pi(\alpha)\cap\pi(\sigma)$, it holds
$b\in\pi(\alpha\cap\alpha^{\prime})*O$. So we can assume that $b$ is on the
line segment $\overline{b^{\prime}O}$, where
$b^{\prime}\in\pi(\alpha\cap\alpha^{\prime})$. This implies that both $b$ and
$b^{\prime}$ are in $\pi(\alpha)$. Since $O\not\in\pi(\alpha)$, then
$b=b^{\prime}$. This shows $\pi(\alpha)\cap\pi(\sigma)=\pi(\alpha\cap\sigma)$.
∎

######  Lemma 2.10.

For any $\omega_{1},\omega_{2}\in\mathcal{C}$, we have
$\omega_{1}\cap\omega_{2}\in\mathcal{C}$.

###### Proof.

If $O\not\in\omega_{1}$, then by Corollary 2.5, there exists a unique
$\alpha\in\mathcal{O}$ such that $\pi(\alpha)=\omega_{1}$. Let $\sigma$ be a
face in $\mathcal{O}$ such that $\pi(\sigma)=\omega_{2}$. By Corollary 2.9 we
have that $\omega_{1}\cap\omega_{2}=\pi(\alpha\cap\sigma)$, and so
$\omega_{1}\cap\omega_{2}\in\mathcal{C}$.

If $O$ is in both $\omega_{1}$ and $\omega_{2}$, then by Corollary 2.6, there
exists $\sigma_{l}\in\mathcal{O}$ such that $\pi(\sigma_{l})=\omega_{l}$ and
$\tau\subseteq\sigma_{l}$, for $l=1,2$ respectively. By Lemma 2.8, it follows
$\omega_{1}\cap\omega_{2}=\pi(\sigma_{1}\cap\sigma_{2})\in\mathcal{C}$. ∎

Now we are ready to prove Proposition 2.1.

###### Proof of Proposition 2.1.

The fact that $\mathcal{C}$ is a simplicial complex follows from Corollary 2.3
and Lemma 2.10. Note that $O$ is the only interior vertex of $\mathcal{C}$ and
every facet of $\mathcal{C}$ contains $O$, so $\mathcal{C}$ must be the
closure of the star of the vertex $O$. ∎

##  3 Cofactors

In this section and thereafter, we use the notion of filtered vector spaces.
Recall that a _filtered vector space_ is a vector space $V$ with a nested
sequence of subspaces $\bigl{\\{}V_{\leqslant d}\subseteq V\colon
d=0,1,2,\dots\bigr{\\}}$ such that

$\\{0\\}\subseteq V_{\leqslant 0}\subseteq V_{\leqslant 1}\subseteq
V_{\leqslant 2}\subseteq\cdots$

and that

$V=\bigcup_{d\geqslant 0}V_{\leqslant d}.$

In this article, we always assume $V_{\leqslant d}$ is finite dimensional for
each $d$. Let $\Delta$ be a fixed finite $k$-dimensional triangulation of a
domain $\Omega\subseteq\mathbb{R}^{k}$. If $r\geqslant 0$ is an integer, the
_total spline space_ over $\Delta$ is defined as

$S^{r}(\Delta)=\bigcup_{d\geqslant 0}S^{r}_{d}(\Delta).$ (11)

In particular, $S^{r}(\Delta)$ is a filtered vector space. Denote by
$\Delta_{j}$ the set of $j$-dimensional faces (or $j$-faces) of $\Delta$, for
$0\leqslant j\leqslant k$. We say that two $k$-faces
$\sigma,\sigma^{\prime}\in\Delta_{k}$ are _maximal adjacent faces_ , or simply
_adjacent_ , if their intersection is a $(k-1)$-face of $\Delta$ i.e., if
$\sigma$ and $\sigma^{\prime}$ share (or have in common) a $(k-1)$-face of
$\Delta$.

We recall from [7] that $f=(f_{\sigma}\colon\sigma\in\Delta_{k})\in
S^{r}(\Delta)$ if and only if for every pair
$\sigma,\sigma^{\prime}\in\Delta_{k}$ such that
$\sigma\cap\sigma^{\prime}=\varepsilon\in\Delta_{k-1}$ there exists a
polynomial $c_{\varepsilon}$ such that
$f_{\sigma}-f_{\sigma^{\prime}}=c_{\varepsilon}\cdot\ell_{\varepsilon}^{r+1}$,
where $\ell_{\varepsilon}$ is the linear polynomial defining the hyperplane
containing $\varepsilon$. This result is called the _cofactor criterion_ for
$C^{r}$-splines.

Throughout this section, we follow the notation introduced in Section 2, and
denote by $\mathcal{O}$ a $(k,i)$-orange. We embedded $\mathcal{O}$ in
$\mathbb{R}^{k}$, and up to a change of coordinates assume that the medial
simplex $\tau$ of $\mathcal{O}$ satisfies $x_{1}=\cdots=x_{i}=0$. In these
coordinates, we consider the projection $\pi$ along $\tau$ given in (2). We
denote by $\mathcal{C}$ the image of $\mathcal{O}$ by $\pi$, as defined in
(3).

The orange $\mathcal{O}$ is by definition a $k$-dimensional simplicial complex
(see Definition 1.1), and Proposition 2.1 implies that $\mathcal{C}$ is an
$i$-dimensional simplicial complex. Following the notation above, we denote by
$\mathcal{O}_{k}$ the set of $k$-dimensional (or maximal) faces of
$\mathcal{O}$.

The following proposition shows how the total $C^{r}$-spline spaces over
$\mathcal{O}$ and $\mathcal{C}$ are related.

######  Proposition 3.1.

For any integer $r\geqslant 0$, the total spline space over $\mathcal{O}$
satisfies

$S^{r}(\mathcal{O})\simeq
S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}]$ (12)

as filtered vector spaces, that is, for each degree $d$

$S^{r}_{d}(\mathcal{O})\simeq(S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}])_{\leqslant
d},$

where
$(S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}])_{\leqslant
d}$ is the space spanned by

$\bigl{\\{}(g,h)\in
S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}]\colon\deg
g+\deg h\leqslant d\bigr{\\}}.$

###### Proof.

By Corollary 2.7, the projection $\pi$ establishes a bijective correspondence
between the maximal faces of $\mathcal{O}$ and those of $\mathcal{C}$. Hence,
every spline $g\in S^{r}(\mathcal{C})$ can be written as the tuple
$g=\bigl{(}g_{\pi(\sigma)}\colon\sigma\in\mathcal{O}_{k}\bigr{)}$, where
$g_{\pi(\sigma)}=g|_{\pi(\sigma)}\in\mathbb{R}[x_{1},\dots,x_{i}]$. Let
$R=\mathbb{R}[x_{1},\dots,x_{k}]$, and consider the map

$\displaystyle\varphi\colon(S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}])_{\leqslant
d}$ $\displaystyle\to\bigoplus_{\sigma\in\mathcal{O}_{k}}R,$
$\displaystyle\left(g,h\right)$
$\displaystyle\mapsto\bigl{(}hg_{\pi(\sigma)}\colon\sigma\in\mathcal{O}_{k}\bigr{)}.$

Notice that
$S_{c}^{r}(\mathcal{O})\subseteq\bigoplus_{\sigma\in\mathcal{O}_{k}}R$. We
want to show that $\varphi$ is an isomorphishm of vector spaces and
$\mathrm{Im}\,\varphi=S_{d}^{r}(\mathcal{O})$.

It is clear that $\varphi$ is $\mathbb{R}$-linear and it is injective because
$S^{r}(\mathcal{C})\subseteq\mathbb{R}[x_{1},\dots,x_{i}]$, so none of the
variables $x_{t}$ for $t>i$ is involved in the polynomials $g_{\pi}(\sigma)$.
First we show that $\mathrm{Im}\,\varphi\subseteq S^{r}_{d}(\mathcal{O})$.

Let $g\in S^{r}(\mathcal{C})$. Then, for any
$h\in\mathbb{R}[x_{i+1},\dots,x_{k}]$, and any pair of adjacent faces
$\sigma,\sigma^{\prime}\in\mathcal{O}_{k}$ such that
$\sigma\cap\sigma^{\prime}=\varepsilon\in\mathcal{O}_{k}$, the cofactor
criterion implies
$hg_{\pi(\sigma)}-hg_{\pi(\sigma^{\prime})}=hc_{\pi(\varepsilon)}\cdot\ell_{\pi(\varepsilon)}^{r+1}$,
for some polynomial $c_{\pi(\varepsilon)}\in\mathbb{R}[x_{1},\dots,x_{i}]$,
where $\ell_{\pi(\varepsilon)}$ is the linear polynomial vanishing on
$\pi(\varepsilon)$. Notice that by construction, the $(k-1)$-face
$\varepsilon$ contains the medial simplex $\tau$ of $\mathcal{O}$.
Consequently, the linear polynomial $\ell_{\varepsilon}$ vanishing on
$\varepsilon$ is in $\mathbb{R}[x_{1},\dots,x_{i}]$, and hence
$\ell_{\pi(\varepsilon)}=\ell_{\varepsilon}$. Hence, $\varphi(g,h)\in
S^{r}(\mathcal{O})$. By filtering on the degree we have that
$\deg\varphi(g,h)\leqslant d$. This shows that
$\mathrm{Im}\,(\varphi)\subseteq S^{r}_{d}(\mathcal{O})$.

We now prove that $\mathrm{im}\varphi\supseteq S^{r}_{d}(\mathcal{O})$. Let
$f=(f_{\sigma}\colon\sigma\in\mathcal{O}_{k})\in S^{r}_{d}(\mathcal{O})$,
where $f|_{\sigma}=f_{\sigma}\in R$ and $\deg f_{\sigma}\leqslant d$. Notice
that we can rewrite each $f_{\sigma}$ as a polynomial in $x_{i+1},\dots,x_{k}$
with coefficients in $\mathbb{R}[x_{1},\dots,x_{k}]$. More precisely, we have

$f_{\sigma}=\sum_{{\bm{j}}\in\mathbb{Z}_{\geqslant
0}^{k-i},\;|\bm{j}|\leqslant
d}a_{\bm{j},\sigma}\,\bm{y}^{\bm{j}},\quad\text{for
$a_{\bm{j},\sigma}\in\mathbb{R}[x_{1},\dots,x_{i}]$},$

where $\bm{y}=(x_{i+1},\dots,x_{k})$, and for a tuple of non-negative integers
$\bm{j}=(j_{i+s})_{s=1}^{k-i}$ we write $|\bm{j}|=\sum_{s=1}^{k-i}j_{s}$ and
$\bm{y}^{\bm{j}}=x_{i+1}^{j_{i+1}}\cdots x_{k}^{j_{k}}$. Then, if
$a_{\bm{j}}=\bigl{(}a_{\bm{j},\sigma}\colon\sigma\in\mathcal{O}_{k}\bigr{)}$,
we can write

$f=\sum_{{\bm{j}}\in\mathbb{Z}_{\geqslant 0}^{k-i},\;|\bm{j}|\leqslant
d}{a}_{\bm{j}}\,\bm{y}^{\bm{j}},$ (13)

where $\deg{a}_{\bm{j}}=\max\\{\deg
a_{\bm{j},\sigma}\colon\sigma\in\mathcal{O}_{k}\\}$. We prove that
$a_{\bm{j}}\in S^{r}(\mathcal{C})$ for every tuple
$\bm{j}\in\mathbb{Z}_{\geqslant 0}$ in (13) as follows. Notice that for
adjacent faces $\sigma$ and $\sigma^{\prime}$ as above, we have that
$f_{\sigma}-f_{\sigma^{\prime}}=c_{\varepsilon}\cdot\ell_{\varepsilon}^{r+1}$,
and we can write

$c_{\varepsilon}=\sum_{{\bm{j}}\in\mathbb{Z}_{\geqslant
0}^{k-i},\;|\bm{j}|\leqslant
d}b_{\bm{j},\varepsilon}\cdot\bm{\bm{y}}^{\bm{j}},\quad\text{where
}b_{\bm{j},\varepsilon}\in\mathbb{R}[x_{1},\dots,x_{i}].$

Thus,

$f_{\sigma}-f_{\sigma^{\prime}}=\sum_{{\bm{j}}\in\mathbb{Z}_{\geqslant
0}^{k-i},\;|\bm{j}|\leqslant
d}(b_{\bm{j},\varepsilon}\ell_{\varepsilon}^{r+1})\cdot{\bm{y}}^{\bm{j}}=\sum_{{\bm{j}}\in\mathbb{Z}_{\geqslant
0}^{k-i},\;|\bm{j}|\leqslant
d}(a_{\bm{j},\sigma}-a_{\bm{j},\sigma^{\prime}})\bm{y}^{\bm{j}}.$

Hence, for every $\bm{j}\in\mathbb{Z}_{\geqslant 0}$, we have
$a_{\bm{j},\sigma}-a_{\bm{j},\sigma^{\prime}}=b_{\bm{j},\varepsilon}\ell_{\varepsilon}^{r+1}$,
which implies that $a_{\bm{j}}\in S^{r}(\mathcal{C})$. Note that $\deg
f_{\sigma}\leqslant d$, so $\deg a_{\bm{j}}+|\bm{j}|\leqslant d$. This shows
that $\mathrm{Im}\,\varphi\supseteq S^{r}_{d}(\mathcal{O})$, and it completes
the proof. ∎

One immediate consequence of Proposition 3.1 is that the algebraic structure
of $S^{r}_{d}(\mathcal{O})$ only depends on that of $S^{r}(\mathcal{C})$ and
the $(k-i)$-subspace in which the medial simplex $\tau$ lies. Note that the
projection map $\pi$ only depends on the $(k-i)$-subspace containing $\tau$,
we have the following corollary.

######  Corollary 3.2.

Let $\mathcal{O}$ and $\mathcal{O}^{\prime}$ be two $(k,i)$-oranges with
medial simplices $\tau$ and $\tau^{\prime}$, respectively. If
$\tau,\tau^{\prime}$ lie in the same $(k-i)$-subspace and the projected
oranges $\mathcal{C}=\mathcal{C}^{\prime}$, then for any degree $d$,

$S^{r}_{d}(\mathcal{O})\simeq S^{r}_{d}(\mathcal{O}^{\prime}).$

In particular,

$\dim S^{r}_{d}(\mathcal{O})=\dim S^{r}_{d}(\mathcal{O}^{\prime}).$

###### Proof.

By Proposition 3.1 it is clear that both $S^{r}_{d}(\mathcal{O})$ and
$S^{r}_{d}(\mathcal{O}^{\prime})$ are isomorphic to
$(S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}])_{\leqslant
d}$. Therefore, they are isomorphic and have the same dimension. ∎

##  4 Hilbert series

In this section we prove the main result in this paper which is Theorem 1.1.

Note that for any order simplicial complex $\Delta$ and any order of
continuity $r\geqslant 0$, we can see $S^{r}(\Delta)$ as a filtered vector
space for the sequence of spline spaces $\\{0\\}\subseteq
S_{0}^{r}(\Delta)\subseteq S_{1}^{r}(\Delta)\subseteq\cdots$. The _Hilbert
series_ $\mathrm{Hilb}\bigl{(}S^{r}(\Delta),t\bigr{)}$ of $S^{r}(\Delta)$ is
by definition

$\mathrm{Hilb}\bigl{(}S^{r}(\Delta),t\bigr{)}=\sum_{d=0}^{\infty}\dim
S^{r}_{d}(\Delta)\,t^{d}.$

Following the notation in previous sections, if $0\leqslant i\leqslant k$ we
denote by $\mathcal{O}$ a $(k,i)$-orange and by $\mathcal{C}$ the image of
$\mathcal{O}$ by the projection $\pi\colon\mathbb{R}^{k}\to\mathbb{R}^{i}$ in
(2). The following result is a consequence of Proposition 3.1 and relates the
Hilbert series of the spline spaces $S^{r}(\mathcal{O})$ and
$S^{r}(\mathcal{C})$.

######  Corollary 4.1.

If $\mathcal{O}$ is a $(k,i)$-orange and $\mathcal{C}=\pi(\mathcal{O})$ as
defined in (3), then

$\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{O}),t\bigr{)}=\frac{1}{(1-t)^{k-i}}\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{C}),t\bigr{)}.$
(14)

###### Proof.

By Proposition 3.1, we know that
$S^{r}(\mathcal{O})=S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}]$
and
$S^{r}_{d}(\mathcal{O})\cong(S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}])_{\leqslant
d}$ for every $d\geqslant 0$. We will prove (14) by induction on the variables
$x_{i+j}$ for $1\leqslant j\leqslant k-i$. Define
$\phi\colon\bigl{(}S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}]\bigr{)}_{d}\to
S^{r}_{d}(\mathcal{C})$ by taking $\phi\bigl{(}x_{i+1}\bigr{)}=0$. Then,
$\ker_{d}(\phi)=\bigl{(}S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}]\bigr{)}_{d-1}\cdot
x_{i+1}\cong\bigl{(}S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}]\bigr{)}_{d-1}$.
This holds for every degree $d\geqslant 0$, and therefore it implies

$\mathrm{Hilb}(\ker(\phi),t)=\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{C},t)\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}],t\bigr{)}\cdot
t,$ (15)

where $\ker(\phi)$ is the filtered vector of $\ker_{d}(\phi)$. But $\phi$ is
surjective and linear, hence $\dim\ker_{d}(\phi)+\dim
S^{r}(\mathcal{C})_{d}=\dim\bigl{(}S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}]\bigr{)}_{d}$,
and this together with (15) yields

$\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1}],t\bigr{)}(1-t)=\mathrm{Hilb}(S^{r}(\mathcal{C})).$

Since
$\mathbb{R}[x_{i+1},\dots,x_{k}]\cong\mathbb{R}[x_{i+1}]\otimes_{\mathbb{R}}\dots\otimes_{\mathbb{R}}\mathbb{R}[x_{k}]$,
and
$S^{r}(\mathcal{O})=S^{r}(\mathcal{C})\otimes_{\mathbb{R}}\mathbb{R}[x_{i+1},\dots,x_{k}]$,
we easily see by induction that
$\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{O}),t\bigr{)}(1-t)^{k-i}=\mathrm{Hilb}\bigl{(}S^{r}(\mathcal{C}),t\bigr{)}$,
which proves (14). ∎

Theorem 1.1 can be proved from Corollary 4.1 as follows.

###### Proof Theorem 1.1.

Notice that we can rewrite

$\frac{1}{(1-t)^{k-i}}=\sum_{j=0}^{\infty}\binom{k-i-1+j}{j}t^{j}.$

Thus, by Corollary 4.1, we have

$\sum_{d=0}^{\infty}\dim
S_{d}^{r}(\mathcal{O})t^{d}=\left(\sum_{j=0}^{\infty}\binom{k-i-1+j}{j}t^{j}\right)\cdot\left(\sum_{j=0}^{\infty}\dim
S_{j}^{r}(\mathcal{C})t^{j}\right),$

which implies

$\displaystyle\dim S_{d}^{r}(\mathcal{O})$
$\displaystyle=\sum_{\begin{subarray}{c}j+l=d\\\ j,\,l\geqslant
0\end{subarray}}\binom{k-(i+1)+l}{l}\dim S^{r}_{j}(\mathcal{C})$
$\displaystyle=\sum_{j=0}^{d}\binom{d+k-i-j-1}{k-i-1}\dim
S^{r}_{j}(\mathcal{C}),$

for every degree $d\geqslant 0$. ∎

##  5 Bernstein-Bézier techniques

The idea is to use Corollary 3.2 to transform $\mathcal{O}$ into a special
$\mathcal{O}^{\prime}$ constructed from the projection $\mathcal{C}$, so that
its geometric structure allows to “lift” the Bernstein-Bézier basis on the
projected orange $\mathcal{C}$ to a Bernstein-Bézier basis in
$\mathcal{O}^{\prime}$. The main advantage of Bernstein-Bézier techniques is
the tool called domain points that essentially replaces basis functions with
actual points located in specified positions in each simplex. We assume
familiarity of the reader with basic Bernstein-Bézier concepts, and refer to
[10] for a comprehensive treatment of bivariate and trivariate case as well as
to a survey paper [1].

Recall that by assumption the medial face $\tau$ of $\mathcal{O}$ is in the
plane $x_{1}=\cdots=x_{i}=0$. Without loss of generality assume that the first
vertex $v_{0}$ of $\tau$ is at the origin $O$, and the remaining vertices are
$v_{1},\dots v_{k-i}$. We also assume that $\mathcal{C}$ is embedded in the
subspace $x_{i+1}=\dots=x_{k}=0$, and its only interior vertex is located at
$O$. We construct $\mathcal{O}^{\prime}$ as follows:

$\mathcal{O}^{\prime}:=\mathcal{C}*\tau=\bigl{\\{}\omega*v_{1}*\dots*v_{k-i}\colon\omega\in\mathcal{C}\bigr{\\}}.$

It is easy to verify that $\mathcal{O}^{\prime}$ is also a $(k,i)$-orange. It
has the same medial simplex $\tau$ as $\mathcal{O}$, and its projected orange
$\mathcal{C}^{\prime}$ equals $\mathcal{C}$.

The set of all domain points in $\tau$ for a polynomial of degree at most $d$
in $(k-i)$ variables is given by:

${\mathcal{D}}^{\tau}:=\bigl{\\{}\xi^{\tau}_{j_{0},j_{1}\dots,j_{k-i}}=(j_{1}v_{1}+\cdots+j_{k-i}v_{j_{k-i}})/d\colon
j_{0}+j_{1}+\cdots+j_{k-i}=d\bigr{\\}}.$

Next we define a family of scaled versions of $\mathcal{C}$ as follows

$\mathcal{C}_{j}:=(j/d)\,\mathcal{C},\;\text{for }j=0,\dots,d,$

where if $j=0$, the projected orange $\mathcal{C}$ scales to the origin. Let
$\omega$ be a maximal simplex in $\mathcal{C}$ with the first vertex $u_{0}$
at the origin $O$, and the remaining vertices $u_{1},\dots u_{i}$. We also
define a a family of scaled versions of $\omega$ as follows

$\omega_{j}:=(j/d)\,\omega,\;\text{for }j=0,\dots,d,$

where if $j=0$, the simplex $\omega$ scales to the origin.

$\mathcal{C}$$\mathcal{C}_{2}$$\mathcal{C}_{1}$$\mathcal{C}_{0}$ Figure 2:
Four parallel layers of domain points forming $\mathcal{D}^{\prime}$ for
$S_{3}^{1}$ on a $(2,1)$-orange with the associated spaces of univariate
splines $S_{j}^{r}(\mathcal{C}_{j})$, for $j=0,1,2,3$, and
$\mathcal{C}=[-1,0]\cup[0,1]$.

The set of all domain points in $\omega_{j}$ for a polynomial of degree at
most $j$ in $i$ variables is given by:

${\mathcal{D}}_{j}^{\omega}:=\bigl{\\{}\xi^{\omega,j}_{\ell_{0},\ell_{1}\dots,\ell_{i}}=\frac{j}{d}(\ell_{1}u_{1}+\cdots+\ell_{i}u_{i})/j=(\ell_{1}u_{1}+\cdots+\ell_{i}u_{i})/d\colon\\\
\ell_{0}+\ell_{1}+\cdots+\ell_{i}=j\bigr{\\}}.$

The purpose of the following lemma is to show that the domain points in
$\mathcal{O}^{\prime}$ can be split into parallel layers orthogonal to $\tau$,
and, thus, the smoothness conditions for a spline in
$S_{d}^{r}(\mathcal{O}^{\prime})$ can be also split into independent blocks.

######  Lemma 5.1.

Let $\mathcal{D}^{\prime}$ be the set of all domain points for a spline in
$S_{d}^{0}(\mathcal{O}^{\prime})$. Then $\mathcal{D}^{\prime}$ can be lifted
from the domain points in $\mathcal{C}$ as follows:

$\mathcal{D}^{\prime}=\bigcup_{j=0}^{d}\bigl{\\{}\cup_{\omega\in\mathcal{C}}\mathcal{D}^{\omega}_{j}+\xi_{j,j_{1},j_{2},\dots,j_{k-i}}\colon\xi_{j,j_{1},j_{2},\dots,j_{k-i}}\in\mathcal{D}^{\tau}\bigr{\\}}.$

###### Proof.

Let $\sigma$ be a maximal simplex in $\mathcal{O}^{\prime}$. Therefore,
$\sigma=\omega*\tau$ for some maximal simplex $\omega\in\mathcal{C}$. Let
$\eta_{m_{0},\cdots,m_{k}}$, $m_{0}+\cdots m_{k}=d$ be an arbitrary domain
point in $\mathcal{D}^{\prime}\cap\sigma$. Without loss of generality, assume
that the first vertex $u_{0}$ of $\omega$ is at the origin, the next $i$
vertices $u_{1},\dots,u_{i}$ are in $\mathcal{C}$, and the last $(k-i)$
vertices $v_{1},\dots,v_{k-i}$ are in $\tau$. Then

$\displaystyle\eta_{m_{0},\cdots,m_{k}}=$
$\displaystyle(m_{1}u_{1}+\cdots+m_{i}u_{i}+m_{i+1}v_{1}+\cdots+m_{k}v_{k-i})/d$
$\displaystyle=(m_{1}u_{1}+\cdots
m_{i}u_{i})/d+(m_{i+1}v_{1}+\cdots+m_{k}v_{k-i})/d$
$\displaystyle=:\xi^{\omega,j}_{m_{0},m_{1},\dots,m_{i}}+\xi^{\tau}_{j,m_{i+1},\dots,m_{k}},$

where for $j=d-m_{i+1}-\cdots-m_{k}$, we have
$\xi^{\omega,j}_{m_{0},m_{1},\dots,m_{i}}\in\mathcal{D}^{\omega}_{j}$, and
$\xi^{\tau}_{j,m_{i+1},\dots,m_{k}}\in\mathcal{D}^{\tau}$, since
$j+m_{i+1}+\cdots m_{k}=d.$

Conversely, for $0\leqslant j\leqslant d$,
$\ell_{0}+\ell_{1}+\cdots+\ell_{i}=j$, and $j+j_{1}+\cdots+j_{k-i}=d$,
consider

$\xi^{\omega,j}_{\ell_{0},\ell_{1},\dots,\ell_{i}}+\xi^{\tau}_{j,j_{1},\dots,j_{k-i}}=(\ell_{1}u_{1}+\cdots+\ell_{i}u_{i})/d+(j_{1}v_{1}+\cdots+j_{k-1}v_{k-i})/d\\\
=\ell_{1}u_{1}+\cdots+\ell_{i}u_{i}+j_{1}v_{1}+\cdots+j_{k-1}v_{k-i})/d\;=\eta_{\ell_{0},\ell_{1},\dots,\ell_{i},j_{1},\dots,j_{k-i}},$

since $\ell_{0}+\ell_{1}+\cdots+\ell_{i}+j_{1}+\cdots+j_{k-i}=j+d-j=d$, and
the simplex in $\mathcal{O}^{\prime}$ in $\mathbb{R}^{k}$ has $(k+1)$ vertices
$O,u_{1},\cdots,u_{i},v_{1},\dots,v_{k-i}$. ∎

1 long solid segment with 7 points is $\mathcal{C}$2 medium dashed segments
with 5 points are shifts of $\mathcal{C}_{2}$3 short dashed segments with 3
points each are shifts of $\mathcal{C}_{1}$$v_{1}$4 points on $[v_{1},v_{2}]$
are shifts of $\mathcal{C}_{0}$$v_{2}$ Figure 3: Domain points for
$S^{0}_{3}(\mathcal{O}^{\prime})$: large dots are the domain points in $\tau$;
small dots are the domain points on shifts of $\mathcal{C}_{j}$,
$j=0,\dots,3$.

Before proving our next result, we demonstrate the approach on two examples.
The first one is a $(2,1)$-orange $\mathcal{O}^{\prime}$ depicted in Fig. 2.
The smoothness conditions for a spline in $S_{3}^{r}(\mathcal{O}^{\prime})$
are easily seen to be the same as for a univariate spline of smoothness $r$ on
the $\mathcal{C}$: a partition of an interval with one node. Moreover, $\dim
S^{r}_{3}(\mathcal{O}^{\prime})=\sum_{j=0}^{3}\dim S^{r}_{j}(\mathcal{C})$.
Our next example is a $(3,1)$-orange $\mathcal{O}^{\prime}$ depicted in Fig.
3. Here again, the smoothness conditions for a spline in
$S_{3}^{r}(\mathcal{O}^{\prime})$ are the same as for a univariate spline of
smoothness $r$ on the $\mathcal{C}$: a partition of an interval with one node.
Moreover,

$\dim S^{r}_{3}(\mathcal{O}^{\prime})=\sum_{j=0}^{3}(4-j)\dim
S^{r}_{j}(\mathcal{C}).$

######  Theorem 5.2.

For each $0\leqslant j\leqslant d$, let ${\mathcal{M}}_{j}$ be a minimal
determining set for a spline in $S_{j}^{r}(\mathcal{C}_{j})$. Then the minimal
determining set ${\mathcal{M}}^{\prime}$ for a spline in
$S_{d}^{r}(\mathcal{O}^{\prime})$ is given by

${\mathcal{M}}^{\prime}=\bigcup_{j=0}^{d}\bigl{\\{}{\mathcal{M}}_{j}+\xi_{j,j_{1},\dots,j_{k-i}}\colon\xi_{j,j_{1},j_{2},\dots,j_{k-i}}\in\mathcal{D}^{\tau}\bigr{\\}}.$

Moreover, the dimension of $S_{j}^{r}(\mathcal{C}_{j})$ can be computed by
adding the corresponding cardinalities as follows:

$\sum_{j=0}^{d}\binom{d+k-j-i-1}{k-i-1}|\mathcal{M}_{j}|.$ (16)

###### Proof.

Lemma 5.1 implies that the domain points in $\mathcal{D}^{\prime}$ lie on
parallel layers that are shifts of $\mathcal{C}_{j}$ by a vector
$\xi_{j,j_{1},j_{2},\dots,j_{k-i}}$ that is orthogonal to $\mathcal{C}_{j}$.
Thus, all smoothness conditions for a spline in
$S_{d}^{r}(\mathcal{O}^{\prime})$ are essentially $i$-dimensional, not
$k$-dimensional, see Fig. 2 and 3. In fact, they are exactly the same as the
ones for $S_{j}^{r}(\mathcal{C})$ for all $j=0,\dots,d$. We now assume that we
know Bernstein-Bézier bases for $S_{j}^{r}(\mathcal{C})$ for all
$j=0,\dots,d$. If we scale $\mathcal{C}$, these Bernstein-Bézier bases do not
change since they are affine invariant. Thus, we know the corresponding
minimal determining sets ${\mathcal{M}}_{j}$ for a spline in
$S_{j}^{r}(\mathcal{C}_{j})$, for $j=0,\dots,d$. We assume that if $j=0$, we
have one point in the corresponding MDS. Thus, we only need to count the
number of $\mathcal{C}_{j}$ layers to complete the proof. Since
$j+j_{1}+\cdots+j_{k-i}=d$, for each $0\leqslant j\leqslant d$, we have
$\binom{d+k-j-i-1}{k-i-1}$ layers of $\mathcal{C}_{j}$, and (16) follows. ∎

Note that (16) is equal to (1). Although it is tempting to claim that the
minimal determining set $\mathcal{M}^{\prime}\subseteq\mathcal{O}^{\prime}$ in
Theorem 5.2 gives rise to a corresponding minimal determining set in
$\mathcal{O}$, this is, of course, not true in general, since the associated
change of basis matrix is not guaranteed to be non-singular. The approximation
order of $S_{d}^{r}(\mathcal{O})$ is hard to define since this partition is
not easily refinable into oranges. Assuming that we can create macro-elements
based on oranges, the only conclusion that can be derived is that the
approximation order of $S_{d}^{r}(\mathcal{O})$ cannot be better than that of
$S_{d}^{r}(\mathcal{C})$.

## Acknowledgements

N. Villamizar and Beihui Yuan were supported by the UK Engineering and
Physical Sciences Research Council (EPSRC) New Investigator Award
EP/V012835/1.

## References

  * [1] P. Alfeld.  Multivariate splines and the Bernstein-Bézier form of a polynomial.  Comput. Aided Geom. Design, 45:2–13, 2016.
  * [2] P. Alfeld, Bruce Piper, and L. Schumaker.  An explicit basis for $C^{1}$ quartic bivariate splines.  SIAM J. Numer. Anal., 24(4):891–911, 1987.
  * [3] P. Alfeld and L. Schumaker.  The dimension of bivariate spline spaces of smoothness $r$ for degree $d\geqslant 4r+1$.  Constr. Approx., 3(2):189–197, 1987.
  * [4] P. Alfeld and L. Schumaker.  On the dimension of bivariate spline spaces of smoothness $r$ and degree $d=3r+1$.  Numer. Math., 57(6-7):651–661, 1990.
  * [5] P. Alfeld, L. Schumaker, and M. Sirvent.  On dimension and existence of local bases for multivariate spline spaces.  J. Approx. Theory, 70(2):243–264, 1992.
  * [6] P. Alfeld, L. Schumaker, and W. Whiteley.  The generic dimension of the space of $C^{1}$ splines of degree $d\geqslant 8$ on tetrahedral decompositions.  SIAM J. Numer. Anal., 30(3):889–920, 1993.
  * [7] L. Billera.  Homology of smooth splines: generic triangulations and a conjecture of Strang.  Trans. Amer. Math. Soc., 310(1):325–340, 1988.
  * [8] M. DiPasquale and N. Villamizar.  A lower bound for splines on tetrahedral vertex stars.  SIAM J. Appl. Algebra Geom., 5(2):250–277, 2021.
  * [9] D. Hong.  Spaces of bivariate spline functions over triangulation.  Approx. Theory Appl., 7(1):56–75, 1991.
  * [10] M.-J. Lai and L. Schumaker.  Spline functions on triangulations, volume 110 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2007.
  * [11] M. Sirvent.  The dimension of the space of multivariate splines.  PhD thesis, Department of Mathematics, University of Utah, 1990.
  * [12] G. Strang.  Piecewise polynomials and the finite element method.  Bull. Amer. Math. Soc., 79:1128–1137, 1973.