A conservative low rank tensor method for the Vlasov dynamics

Wei Guo111 Department of Mathematics and Statistics, Texas Tech University,
Lubbock, TX, 70409. E-mail<EMAIL_ADDRESS>Research is supported by NSF
grant NSF-DMS-1830838 and NSF-DMS-2111383, Air Force Office of Scientific
Research FA9550-18-1-0257.  and Jing-Mei Qiu222Department of Mathematical
Sciences, University of Delaware, Newark, DE, 19716. E-mail<EMAIL_ADDRESS>Research supported by NSF grant NSF-DMS-1818924 and 2111253, Air Force Office
of Scientific Research FA9550-18-1-0257.

Abstract. In this paper, we propose a conservative low rank tensor method to
approximate nonlinear Vlasov solutions. The low rank approach is based on our
earlier work [17]. It takes advantage of the fact that the differential
operators in the Vlasov equation are tensor friendly, based on which we
propose to dynamically and adaptively build up low rank solution basis by
adding new basis functions from discretization of the differential equation,
and removing basis from a singular value decomposition (SVD)-type truncation
procedure. For the discretization, we adopt a high order finite difference
spatial discretization together with a second order strong stability
preserving multi-step time discretization.

While the SVD truncation will remove the redundancy in representing the high
dimensional Vlasov solution, it will destroy the conservation properties of
the associated full conservative scheme. In this paper, we develop a
conservative truncation procedure with conservation of mass, momentum and
kinetic energy densities. The conservative truncation is achieved by an
orthogonal projection onto a subspace spanned by $1$, $v$ and $v^{2}$ in the
velocity space associated with a weighted inner product. Then the algorithm
performs a weighted SVD truncation of the remainder, which involves a scaling,
followed by the standard SVD truncation and rescaling back. The algorithm is
further developed in high dimensions with hierarchical Tucker tensor
decomposition of high dimensional Vlasov solutions, overcoming the curse of
dimensionality. An extensive set of nonlinear Vlasov examples are performed to
show the effectiveness and conservation property of proposed conservative low
rank approach. Comparison is performed against the non-conservative low rank
tensor approach on conservation history of mass, momentum and energy.

Key Words: Low rank; high order SVD; Conservative truncation; Hierarchical
Tuck decomposition of tensors; Mass and momentum conservation; Vlasov
Dynamics.

##  1 Introduction

The Vlasov-Poisson (VP) system is known as a fundamental model in plasma
physics which describes the dynamics of dilute charged particles due to self-
induced electrostatic forces from a statistical viewpoint. The numerical
challenges of realistic Vlasov simulations include the high-dimensionality of
the phase space, features with multiple scales in time and in phase space,
preservation of physical invariants, among many others.

The Particle-In-Cell method employs a collection of sampled macro particles to
represent the distribution function and is widely used for the plasma
simulations mainly due to its distinct dimension-independent rate of
convergence [4, 2]. On the other hand, such a method suffers the inherent
statistical noise, and hence may fail to accurately capture the physics of
interest. For example, if the resolution of the tail of the distribution
function is desired, a noise-free grid-based deterministic method is a better
choice. Over the past two decades, the development of deterministic Vlasov
solvers has attracted lots of research interest, see e.g. [12]. Meanwhile, a
conventional deterministic Vlasov simulation in a realistic and high-
dimensional setting is prohibitively expensive because of the curse of the
dimensionality and the associated huge computational and storage cost. Several
dimension reduction techniques have been developed in the literature to
alleviate the curse of dimensionality for the Vlasov equation and other high-
dimensional partial differential equations (PDEs). One such example is the
sparse grid approach [29, 32, 15], which can effectively reduce the
computational complexity and is well-suited for the problems with moderately
high dimensions. For the Vlasov simulations, we mention the sparse grid semi-
Lagrangian (SL) method [24] and the sparse grid discontinuous Galerkin method
[16, 30]. Recently, the tensor approach emerges as a promising tool for
feasible high-dimensional simulations. Such an approach aims to extract the
underlying low rank structure of the data with advanced tensor decomposition,
potentially breaking the curse of dimensionality. The popular tensor formats
include the canonical polyadic (CP) format [21, 3, 20, 22], Tucker format [31,
5], hierarchical Tucker (HT) format [19, 14], and tensor train (TT) format
[27, 25, 26]. There are several pioneering works employing low rank tensor
approach for nonlinear simulations, including the low rank SL method in the TT
format [23], a low rank method with the CP format based on the underlying
Hamiltonian formulation [7], a dynamical low rank method proposed in [9, 11]
for which the dynamical low rank approximation of the Vlasov solution is
evolved on the low rank manifold using a tangent space projection, and a
dynamical tensor approximations for high dimensional linear and nonlinear PDEs
based on functional tensor decomposition and dynamic tensor approximation [6].

With the existing understanding of the low rank solution structure for the
Vlasov dynamics, as well as the observation that the differential operator in
the Vlasov equation can be represented in the tensorized form, in [17], we
developed a low rank tensor approach to dynamically and adaptively build up
low rank solution basis, by adding new basis functions from discretization of
the PDE and then removing basis from an singular value decomposition
(SVD)-type truncation procedure. In particular, we start from a low rank
solution in a tensor format and add additional basis by applying the well-
established high order finite difference upwind method coupled with the
strong-stability-preserving (SSP) multi-step time discretizations [13]; the
solutions are being further updated by an SVD-type truncation to remove
redundant bases. We further generalize the algorithm to high-dimensional cases
with the HT decomposition, which attains a storage complexity that is linearly
scaled with the dimension, overcoming the curse of dimensionality. An issue
associated with the low rank algorithm in [17] is the loss of mass
conservation in the SVD truncation step. In fact, it is well-known that the VP
system conserves mass, momentum, and energy _locally_ by respecting a set of
macroscopic moment equations. The high order finite difference scheme with an
SSP multi-step time integrator can be shown to locally preserve the
macroscopic mass and momentum, i.e. when taking moments of the kinetic scheme,
one will arrive at a consistent and locally conservative discretization of the
corresponding macroscopic equations. Here “locally conservative” means that
the scheme can be written in a flux difference form, where the fluxes
represent the amount of macroscopic quantities transported between neighboring
computational cells. However, the SVD truncation in the low rank algorithm
would destroy local and global conservation property of the original numerical
scheme.

There exist several techniques developed to correct the conservation errors
for low rank methods in various settings. In [23], the low rank solution is
rescaled so that the total mass is conserved, and a similar mass correction
technique is proposed in [28] for a dynamical low rank method. In [1], moment
fitting is applied to the low rank solution so that the corrected moments
match those solved from the macroscopic fluid equations. In [10], a dynamical
low rank method with Lagrangian multipliers is developed to improve
conservation properties for the total mass and momentum as well as local
projected moment equations. More recently, along the same line, the truly
local conservation of mass, momentum, and energy is attained for the dynamical
low rank method [8]. The idea is to fix certain basis functions in the
dynamical low rank approximation and employ a modified Petrov–Galerkin
formulation which is compatible with the remainder of the approximation.

Inspired by the recent work in [8], in this paper, we develop a novel locally
mass, momentum and energy conservative truncation algorithm under the
framework of the low rank tensor approach [17]. In particular, we apply an
orthogonal projection operator to the low rank solution in a weighted inner
product spaces of phase variables, to extract exactly the first few moments
(i.e. the mass, momentum and kinetic energy density functions) of the low rank
solution after the “add basis” step; a weighted SVD truncation is then applied
to the remainder of the projection. Such a truncation is called “conservative
truncation” throughout the paper. We further develop the conservative
truncation algorithm for the 2D2V VP system with the HT tensor format using a
dimension tree that separates the spatial and phase variables. For high order
HT tensors, an additional projection step is needed, after the high order SVD
(HOSVD) truncation of the remaining term, to ensure conservation of
macroscopic variables from the first few moments. The low rank tensor
algorithm with the conservative projection is theoretically proved and
numerically verified to be locally mass and momentum conservative.

This paper is organized as follows. In Section 2, we introduce the kinetic
Vlasov model and the corresponding macroscopic conservation laws. In Section
3, we first review the low rank tensor approach for the 1D1V Vlasov equation
in Section 3.1; then we introduce the orthogonal projection at the continuous
level to extract exactly the first few moments of Vlasov solution in Section
3.2; we further extend such an orthogonal projection to the discrete level and
develop the main conservative truncation algorithm in Section 3.3. In Section
4, we develop the conservative truncation algorithm for the 2D2V Vlasov model,
and show that the low rank algorithm locally conserves the mass and momentum
at the macroscopic level. In Section 5, we present an extensive set of 1D1V
and 2D2V numerical results to demonstrate the effectiveness and the
conservation properties of the low rank tensor algorithm. We conclude the main
contributions of the paper and comment on future directions in Section 6.

##  2 The kinetic Vlasov model and the corresponding macroscopic systems

We consider the VP system

$\frac{\partial f}{\partial
t}+{\bf{v}}\cdot\nabla_{\bf{x}}f+{\bf{E}}({\bf{x}},t)\cdot\nabla_{\bf{v}}f=0,$
(2.1) ${\bf E}({\bf x},t)=-\nabla_{\bf x}\phi({\bf x},t),\quad-\triangle_{\bf
x}\phi({\bf x},t)={{\bf\rho}({\bf x},t)}-\rho_{0},$ (2.2)

which describes the probability distribution function $f({\bf x},{\bf v},t)$
of electrons in collisionless plasma. Here ${\bf E}$ is the electric field,
and $\phi$ is the self-consistent electrostatic potential. $f$ couples to the
long range fields via the charge density, ${\bf\rho}({\bf
x},t)=\int_{\Omega_{{\bf v}}}f({\bf x},{\bf v},t)d{\bf v}$, where we take the
limit of uniformly distributed infinitely massive ions in the background. The
VP system describes the movement of electrons due to self-induced electric
field ${\bf E}$ determined by the Poisson equation.

The Vlasov dynamics are well-known to conserve several physical invariants. In
particular, let

mass density: $\displaystyle\quad\rho({\bf x},t)=\int_{\Omega_{{\bf v}}}f({\bf
x},{\bf v},t)d{\bf v},$ (2.3) current density: $\displaystyle\quad{\bf J}({\bf
x},t)=\int_{\Omega_{{\bf v}}}f({\bf x},{\bf v},t){\bf v}d{\bf v},$ (2.4)
kinetic energy density: $\displaystyle\quad\kappa({\bf
x},t)=\frac{1}{2}\int_{\Omega_{{\bf v}}}{\bf v}^{2}f({\bf x},{\bf v},t)d{\bf
v},$ (2.5) energy density: $\displaystyle\quad e({\bf x},t)=\kappa({\bf
x},t)+\frac{1}{2}{\bf E}({\bf x})^{2}.$ (2.6)

Then, by taking a few first moments of the Vlasov equation, the following
conservation laws of mass, momentum and energy can be derived

$\displaystyle\rho_{t}+\nabla_{\bf x}\cdot{\bf J}$ $\displaystyle=0$ (2.7)
$\displaystyle\partial_{t}{\bf J}+\nabla_{{\bf x}}\cdot\mathbf{\sigma}$
$\displaystyle=\rho{\bf E}$ (2.8) $\displaystyle\partial_{t}e+\nabla_{{\bf
x}}\cdot\mathbf{Q}$ $\displaystyle=0,$ (2.9)

where $\sigma(t,{\bf x})=\int_{\Omega_{{\bf v}}}({\bf v}\otimes{\bf v})f({\bf
x},{\bf v},t)d{\bf v}$ and $\mathbf{Q}({\bf
x},t)=\frac{1}{2}\int_{\Omega_{{\bf v}}}{\bf v}{\bf v}^{2}f({\bf x},{\bf
v},t)d{\bf v}$. It is well-known that the local conservation property of
numerical schemes is critical to capture correct entropy solutions of
hyperbolic systems such as macroscopic equations (2.7)-(2.9).

##  3 A low rank tensor approach for the Vlasov dynamics with local
conservation

For simplicity of illustrating the basis idea, we only discuss a 1D1V example
in this section.

###  3.1 Review of a low rank tensor approach for Vlasov dynamics [17]

The low rank tensor approach [17] is designed based on the assumption that our
solution at time $t$ has a low-rank representation in the form of

$f(x,v,t)=\sum_{l=1}^{r}\left(C_{l}(t)\ \
U_{l}^{(1)}(x,t)U_{l}^{(2)}(v,t)\right),$ (3.1)

where $\left\\{U_{l}^{(1)}(x,t)\right\\}_{l=1}^{r}$ and
$\left\\{U_{l}^{(2)}(v,t)\right\\}_{l=1}^{r}$ are a set of time-dependent low
rank unit length orthogonal basis in $x$ and $v$ directions, respectively,
$C_{l}$ is the coefficient for the basis $U_{l}^{(1)}(x,t)U_{l}^{(2)}(v,t)$,
and $r$ is the representation rank. (3.1) can be viewed as a Schmidt
decomposition of functions in $(x,v)$ by truncating small singular values up
to rank $r$.

We assume a finite difference discretization of $f$ on a truncated 1D1V domain
of $[x_{\min},x_{\max}]\times[-v_{\max},v_{\max}]$ with uniform $N_{x}\times
N_{v}$ grid points

$x_{\text{grid}}:\quad x_{\min}=x_{1}<\cdots<x_{i}<\cdots<x_{N_{x}}=x_{\max},$
(3.2) $v_{\text{grid}}:\quad-
v_{\max}=v_{1}<\cdots<v_{j}<\cdots<v_{N_{v}}=v_{\max}.$ (3.3)

The numerical solution ${\bf f}\in\mathbb{R}^{N_{x}\times N_{v}}$, as an
approximation to point values of the solution on the grids (3.2)-(3.3), has
the corresponding low-rank counterpart to (3.1) as

${\bf f}=\sum_{l=1}^{r}\left(C_{l}\ \ {\bf U}_{l}^{(1)}\otimes{\bf
U}_{l}^{(2)}\right),\quad(\mbox{or element-wise:}\quad{\bf
f}_{ij}=\sum_{l=1}^{r}C_{l}{\bf U}_{l,i}^{(1)}{\bf U}_{l,j}^{(2)}),$ (3.4)

where ${\bf U}_{l}^{(1)}\in\mathbb{R}^{N_{x}}$ and ${\bf
U}_{l}^{(2)}\in\mathbb{R}^{N_{v}}$ can be viewed as approximations to
corresponding basis functions in (3.1). (3.4) can also be viewed as an SVD of
the matrix ${\bf f}\in\mathbb{R}^{N_{x}\times N_{v}}$. The associated storage
cost is $\mathcal{O}(rN)$, where we assume $N=N_{x}=N_{v}$.

Our low rank tensor approach adaptively updates low-rank basis and associated
coefficients by two steps: an adding basis step by conservative hyperbolic
solvers and a removing basis step via an SVD-type truncation. We consider a
simple first order forward Euler discretization of 1D1V Vlasov equation (2.1)
to illustrate the main idea. We assume the solution in the form of (3.4) with
superscript $n$ for the solution at $t^{n}$.

  1. 1.

Add basis and obtain an intermediate solution ${\bf f}^{n+1,*}$. A forward
Euler discretization of (2.1) gives

${f}^{n+1,*}={f}^{n}-\Delta
t(v\partial_{x}({f}^{n})+E^{n}\partial_{v}({f}^{n})).$ (3.5)

Here the electric field $E^{n}$ is solved by a Poisson solver. Thanks to the
tensor friendly form of the Vlasov equation, at the fully discrete level,
${\bf f}^{n+1,*}$ can be represented in the following low-rank format:

${\bf f}^{n+1,*}=\sum_{l=1}^{r^{n}}C_{l}^{n}\left[\left({\bf
U}_{l}^{(1),n}\otimes{\bf U}_{l}^{(2),n}\right)-\Delta t\left(D_{x}{\bf
U}_{l}^{(1),n}\otimes v\star{\bf U}_{l}^{(2),n}+E^{n}\star{\bf
U}_{l}^{(1),n}\otimes D_{v}{\bf U}_{l}^{(2),n}\right)\right],$ (3.6)

where $D_{x}$ and $D_{v}$ represent high order locally conservative upwind
discretization of spatial differentiation terms, and $\star$ denotes an
element-wise multiplication operation. For example the discretization of
$D_{x}{\bf U}_{l}^{(1),n}\otimes v\star{\bf U}_{l}^{(2),n}$ follows

$D^{+}_{x}{\bf U}_{j}^{(1),n}\otimes v^{+}\star{\bf
U}_{j}^{(2),n}+D^{-}_{x}{\bf U}_{j}^{(1),n}\otimes v^{-}\star{\bf
U}_{j}^{(2),n},$ (3.7)

where $D^{+}_{x}$ and $D^{-}_{x}$ are a fifth order finite difference
discretization of positive and negative velocities respectively, with
$v^{+}=\max(v,0)$ and $v^{-}=\min(v,0)$. Similarly, the discretization of
$E^{n}\star{\bf U}_{l}^{(1),n}\otimes D_{v}{\bf U}_{l}^{(2),n}$ follows

${\bf E}^{n,+}\star{\bf U}_{j}^{(1),n}\otimes D^{+}_{v}{\bf
U}_{j}^{(2),n}+{\bf E}^{n,-}\star{\bf U}_{j}^{(1),n}\otimes D^{-}_{v}{\bf
U}_{j}^{(2),n}$ (3.8)

where $D^{+}_{v}$ and $D^{-}_{v}$ are a fifth order upwind finite difference
discretization of positive and negative velocities respectively, with ${\bf
E}^{+}=\max({\bf E},0)$ and ${\bf E}^{-}=\min({\bf E},0)$.

  2. 2.

Remove basis of ${\bf f}^{n+1,*}$ to update solution ${\bf f}^{n+1}$. Since
the number of basis has increased in a single step update, we perform an SVD-
type truncation to remove redundant bases with a prescribed threshold
$\varepsilon$. The truncation step has no guarantee of any mass, momentum or
energy conservation property. The removing basis step costs
$\mathcal{O}(r^{2}N+r^{3})$, where $r$ is the SVD rank of the numerical
solution.

In this two-step process, both the basis and coefficients are updated.
Extensions to schemes with high order spatial and temporal discretizations and
to high dimensional problems, are developed in [17]. The low rank approach
[17] is built upon the classical high order methods for conservation laws and
kinetic equations, yet it optimizes the computational efficiency by
dynamically building low rank global basis and updating the corresponding
coefficients via a SVD truncation procedure. While the SVD truncation
significantly reduces the computational storage and CPU time, it also destroys
the conservation property in the truncation step.

###  3.2 An orthogonal projection with mass, momentum and energy conservation

To preserve the local mass, momentum and energy in the truncation step,
following a similar idea in [8], we define a weighted inner product space and
the associated norm for functions in $v$ as follows,

$L^{2}_{w}(\mathbb{R})=\\{f:\mathbb{R}\rightarrow\mathbb{R}:\|f\|_{w}<\infty\\}.$

where

$\langle f,g\rangle_{w}=\int f(v)g(v)w(v)dv,\quad\|f\|_{w}=\sqrt{\langle
f,f\rangle_{w}}.$ (3.9)

The weight function $w(v)$ is chosen to have exponential decay so that
$N\doteq\text{span}(1,v,v^{2})\subset L^{2}_{w}(\mathbb{R})$. For example, we
can take $w(v)=\exp(-v^{2}/2)$. Now we are ready to define an orthogonal
decomposition of $f$ with respect to the subspace $N$.

In the context of Vlasov dynamics, we consider $f(x,v)$ with a $v$-direction
projection. We first scale function

$\tilde{f}=\frac{1}{w(v)}f.$ (3.10)

Then we find $P_{N}(\tilde{f})\in N$ s.t.

$\int
P_{N}(\tilde{f}(x,v))g(v)w(v)dv=\int\tilde{f}(x,v)g(v)w(v)dv=\int{f}(x,v)g(v)dv,\quad\forall
g\in N.$ (3.11)

By taking $g=1,v,v^{2}$, $P_{N}(\tilde{f}(x,v))w(v)$ preserves the mass
density $\rho(x)$, current density $J(x)$ and kinetic energy density
$\kappa(x)$ of $f$ defined in (2.3)-(2.5). In fact, we have the _conservative_
decomposition of $f$ as in the following Proposition.

######  Proposition 3.1.

Let $w(v)=\exp(-v^{2}/2)$ and the weighted inner product defined as in (3.9).
$f\in L^{2}_{w}(\mathbb{R})$ can be decomposed as

$f=w(v)(P_{{N}}(\tilde{f})+(I-P_{{N}})(\tilde{f}))\doteq
w(v)(\tilde{f}_{1}+\tilde{f}_{2})\doteq f_{1}+f_{2},$ (3.12)

where $f_{1}$ admits an explicit representation

$f_{1}=w(v)\tilde{f}_{1}=w(v)(c_{1}+c_{2}v+c_{3}(v^{2}-c)),$ (3.13)

with $c=\frac{\langle 1,v^{2}\rangle_{w}}{\langle 1,1\rangle_{w}}$,
$c_{1}=\frac{\rho(x)}{\|1\|_{w}^{2}}$, $c_{2}=\frac{J(x)}{\|v\|_{w}^{2}}$, and
$c_{3}=\frac{2\kappa(x)-\rho(x)c}{\|v^{2}-c\|_{w}^{2}}$, and $\rho(x)$, $J(x)$
and $\kappa(x)$ are as defined in (2.3)-(2.5). $f_{1}$ preserves the mass,
momentum and kinetic energy density of $f$, while $f_{2}$ has zero mass,
current and kinetic energy density.

###### Proof.

We first compute $c$, so that $\\{1,v,v^{2}-c\\}$ is a set of orthogonal
basis. By taking $g=1,v,v^{2}-c$ in (3.11), we can determine the constants
$c_{1}$, $c_{2}$, and $c_{3}$. By construction of the orthogonal projection,
$f_{1}$ preserves the mass, momentum and kinetic energy density of $f$, while
$f_{2}$ has zero of them. ∎

###  3.3 Low rank truncation with mass, momentum and energy conservation

Similar to the weighted inner product (3.9), the projection operator (3.11)
and the conservative decomposition of $f$ in (3.12), we define their
corresponding discrete analogue. Consider ${\bf f},{\bf
g}\in\mathbb{R}^{N_{v}}$, where ${\bf f}$ and ${\bf g}$ can be viewed as
function values of $f$ and $g$ at $v_{\text{grid}}$ (3.3), respectively. We
define the standard $l^{2}$ inner product as

$\langle{\bf f},{\bf g}\rangle=h_{v}\sum_{j}f_{j}g_{j},$ (3.14)

where $h_{v}$ is the mesh size in $v$-direction, serving as the quadrature
weights for the uniform $v_{\text{grid}}$. We also define weighted inner
product and the associated norm as

$\langle{\bf f},{\bf g}\rangle_{\bf
w}=\sum_{j}f_{j}g_{j}w_{j},\quad\quad\|{\bf f}\|_{\bf w}=\sqrt{\langle{\bf
f},{\bf f}\rangle_{\bf w}}$ (3.15)

${\bf w}\in\mathbb{R}^{N_{v}}$ with $w_{j}=w(v_{j})h_{v}$ is the quadrature
weights for $v$-integration with weight function $w(v)$. Correspondingly, we
let

$l^{2}_{\bf w}=\\{{\bf f}\in\mathbb{R}^{N_{v}}:\|{\bf f}\|_{\bf w}<\infty\\}.$

We denote ${\bf 1}_{v}\in\mathbb{R}^{N_{v}}$ the vector of all ones, ${\bf v}$
and ${\bf v}^{2}$ $\in\mathbb{R}^{N_{v}}$ the coordinates and coordinates
square of $v_{\text{grid}}$ (3.3), respectively. Consider the subspace
$\mathcal{N}\doteq\text{span}\\{{\bf 1}_{v},{\bf v},{\bf v}^{2}\\}\subset
l^{2}_{w}$, a conservative low rank truncation of numerical solution ${\bf
f}\in\mathbb{R}^{N_{x}\times N_{v}}$ written in the low rank form of (3.4) can
be obtained from steps below.

  1. 1.

Compute macroscopic quantities of ${\bf f}$. We compute the discrete
macroscopic charge, current and kinetic energy density ${\bm{\rho}}$, ${\bf
J}$ and ${\bm{\kappa}}\in\mathbb{R}^{N_{x}}$ by quadrature

$\displaystyle\left(\begin{array}[]{l}{\bm{\rho}}\\\ {\bf J}\\\
{\bm{\kappa}}\end{array}\right)=\sum_{l=1}^{r}C_{l}\left\langle{\bf
U}^{(2)}_{l},\left(\begin{array}[]{l}{\bf 1}_{v}\\\ {\bf v}\\\ \frac{1}{2}{\bf
v}^{2}\end{array}\right)\right\rangle\ {\bf U}^{(1)}_{l}.$ (3.22)

They are the discrete analog of (2.3), (2.4) and (2.6), i.e.

$\bm{\rho}(i)=\langle{\bf f}(i,:),{\bf 1}_{v}\rangle,\quad{\bf
J}(i)=\langle{\bf
f}(i,:),\bm{v}\rangle,\quad\bm{\kappa}(i)=\langle\bm{f}(i,:),\frac{1}{2}\bm{v}^{2}\rangle,\quad\quad
i=1,\cdots,{N_{x}}.$

  2. 2.

Scale. Similar to (3.10), we scale ${\bf f}$ as

$\tilde{\bf f}=\frac{1}{{\bf w}}\star{\bf f}=\sum_{l=1}^{r}\left(C_{l}\ \ {\bf
U}_{l}^{(1)}\otimes\left(\frac{1}{{\bf w}}\star{\bf
U}_{l}^{(2)}\right)\right).$ (3.23)

  3. 3.

Project. We perform an orthogonal projection of $\tilde{\bf f}$ with respect
to the inner product (3.15) onto subspace $\mathcal{N}$, i.e.

$\langle P_{\mathcal{N}}(\tilde{\bf f}),{\bf g}\rangle_{\bf
w}=\langle\tilde{\bf f},{\bf g}\rangle_{\bf w},\quad\forall{\bf
g}\in\mathcal{N}.$ (3.24)

Similar to Proposition 3.1, ${\bf w}\star P_{\mathcal{N}}(\tilde{\bf f})$
preserves the mass, momentum and kinetic energy densities in the discrete
sense. We show below in Proposition 3.2 a conservative decomposition of ${\bf
f}$.

######  Proposition 3.2.

${\bf f}$ can be decomposed as

${\bf f}={\bf w}\star(P_{\mathcal{N}}(\tilde{{\bf
f}})+(I-P_{\mathcal{N}})(\tilde{{\bf f}}))\doteq{\bf w}\star(\tilde{{\bf
f}}_{1}+\tilde{{\bf f}}_{2})\doteq{\bf f}_{1}+{\bf f}_{2}$ (3.25)

where $f_{1}$ can be represented as a rank three tensor

${\bf f}_{1}={\bf w}\star\tilde{{\bf f}}_{1}=\sum_{k=1}^{3}\
{\mathcal{M}}_{k}\otimes({\bf w}\star{\bf V}_{k}),$ (3.26)

where ${\bf V}_{1}={\bf 1}_{v}$, ${\bf V}_{2}={\bf v}$, and ${\bf V}_{3}={\bf
v}^{2}-c{\bf 1}_{v}$, and $c=\frac{\langle{\bf 1}_{v},{\bf v}^{2}\rangle_{\bf
w}}{\|{\bf 1}_{v}\|_{\bf w}^{2}}$,

${\mathcal{M}}_{1}=\frac{\bm{\rho}}{\|{\bf 1}_{v}\|_{\bf
w}^{2}},\quad{\mathcal{M}}_{2}=\frac{\bf J}{\|{\bf v}\|_{\bf
w}^{2}},\quad{\mathcal{M}}_{3}=\frac{2{\bm{\kappa}}-c{\bm{\rho}}}{\|{\bf
v}^{2}-c{\bf 1}_{v}\|_{\bf w}^{2}}.$ (3.27)

${\bm{\rho}}$, ${\bf J}$ and ${\bm{\kappa}}$ are the discrete mass, momentum
and kinetic energy density of ${\bf f}$ as in (3.22). ${\bf f}_{1}$ preserves
the discrete mass, momentum and kinetic energy density of ${\bf f}$, while
${\bf f}_{2}$ has zero of them, i.e. $\langle{\bf f}_{2}(i,:),{\bf
1}_{v}\rangle=\langle{\bf f}_{2}(i,:),{\bf v}\rangle=\langle{\bf
f}_{2}(i,:),{\bf v}^{2}\rangle=0,\quad i=1,\cdots,{N_{x}}.$

###### Proof.

The proof follows in the same manner as Proposition 3.1, but considering the
discrete weighted inner product and the associated norm (3.15). ∎

  4. 4.

Truncate in $l^{2}_{\bf w}$. To ensure the conservation in the truncation
process, we propose to perform an SVD truncation of $\tilde{\bf f}_{2}$ in the
decomposition (3.25) with respect to the weighted inner product (3.15).
Because of the weight function, we scale $\tilde{\bf f}_{2}$ by ${\sqrt{\bf
w}}\star\tilde{\bf f}_{2}$, where $\sqrt{\cdot}$ is in the element-wise sense
and $\star$ is for the element-wise multiplication. Then, a standard SVD
truncation procedure is applied, denoted as
$\mathcal{T}_{\varepsilon}({\sqrt{\bf w}}\star\tilde{\bf f}_{2})$. Here and
below, $\varepsilon$ is truncation threshold. Lastly, we rescale
$\mathcal{T}_{\varepsilon}({\sqrt{\bf w}}\star\tilde{\bf f}_{2})$ back and
obtain the truncated $\tilde{\bf f}_{2}$. Such a weighted truncation writes

$\mathcal{T}^{\bf w}_{\varepsilon}(\tilde{\bf f}_{2})=\frac{1}{\sqrt{\bf
w}}\star\mathcal{T}_{\varepsilon}(\sqrt{\bf w}\star\tilde{\bf f}_{2}).$ (3.28)

In summary, ${\bf f}_{2}$ is truncated to

${\bf w}\star\mathcal{T}^{\bf w}_{\varepsilon}(\tilde{\bf f}_{2})=\sqrt{\bf
w}\star\mathcal{T}_{\varepsilon}(\sqrt{\bf w}\star\tilde{\bf f}_{2})=\sqrt{\bf
w}\star\mathcal{T}_{\varepsilon}(\frac{1}{\sqrt{\bf w}}\star{\bf f}_{2}).$

  5. 5.

Update. We obtain the low rank truncation of ${\bf f}$ with local mass,
momentum and energy conservation, denoted as

$T_{c}({\bf f})={\bf f}_{1}+{\bf w}\star\mathcal{T}^{\bf
w}_{\varepsilon}(\tilde{\bf f}_{2}).$ (3.29)

We call the proposed truncation (3.29) the conservative truncation. The
following proposition guarantees that the local conservation of mass,
momentum, and kinetic energy density of ${\bf f}$ is preserved in the proposed
truncation procedure of $T_{c}({\bf f})$.

######  Proposition 3.3.

$T_{c}({\bf f})$ has the same discrete charge, current and kinetic energy
density as ${\bf f}$.

###### Proof.

With Proposition 3.2, it is sufficient to show that ${\bf
w}\star\mathcal{T}^{\bf w}_{\varepsilon}(\tilde{\bf f}_{2})$ has zero charge,
current and kinetic energy density. Since $\tilde{\bf f}_{2}$ is orthogonal to
$\mathcal{N}$ with respect to the weighted inner product (3.15), and its
truncation is performed in the same inner product space, $\mathcal{T}^{\bf
w}_{\varepsilon}(\tilde{\bf f}_{2})$ is also orthogonal to $\mathcal{N}$, i.e.

$\langle\mathcal{T}^{\bf w}_{\varepsilon}(\tilde{\bf f}_{2}),g\rangle_{\bf
w}=0,\quad g=1,\,{\bf v},\,{\bf v}^{2},$

i.e. ${\bf w}\star\mathcal{T}^{\bf w}_{\varepsilon}(\tilde{\bf f}_{2})$ has
zero discrete charge, current and kinetic energy density. ∎

Next we establish the local conservation of mass and momentum in the low rank
tensor approach with the conservative truncation (3.29). Since the full
algorithm (without truncation) does not have energy conservation property, the
low rank tensor scheme does not preserve energy conservation. If the
associated full-rank method is able to locally preserve the energy density, so
will the corresponding the low rank method with the proposed conservative
truncation.

######  Proposition 3.4.

(Local mass and momentum conservation of the low rank tensor approach with
conservative truncation.) If the discrete differential operators
$D_{x},\,D_{v}$ employed are conservative, i.e., can be written in a flux
difference form, and linear, i.e., can preserve linear relations, then the
proposed low rank method with an SSP multi-step time integrator preserves the
mass and momentum locally; that is the schemes for ${\bm{\rho}}^{n}$ and ${\bf
J}^{n}$, from integrating the scheme on ${\bf f}^{n}$, are consistent and
conservative discretization of the macroscopic moment equations (2.7)-(2.8).

###### Proof.

Without loss of generality, we prove the proposition for the low rank method
with the forward Euler time integrator for simplicity. By taking discrete
moments of (3.6) and from (3.22), we derive the discrete evolution equations
for the charge density ${\bm{\rho}}$ and current density ${\bf J}$ in
$\mathbb{R}^{N_{x}}$,

$\displaystyle{\bm{\rho}}^{n+1,*}$ $\displaystyle={\bm{\rho}}^{n}-\Delta
tD_{x}({\bf J}^{n}),$ (3.30) $\displaystyle{\bf J}^{n+1,*}$
$\displaystyle={\bf J}^{n}-\Delta
t(D_{x}({\bm{\sigma}}^{n})-{\bm{\rho}}^{n}\star{\bf E}^{n}),$ (3.31)

where we have used the fact that $D_{x}$ and $D_{v}$ are linear, and the
summation by parts

$\langle D_{v}{\bf U}^{(2)}_{l},{\bf v}\rangle=-\langle{\bf
U}^{(2)}_{l},\mathbf{1}_{v}\rangle$ (3.32)

holds by assuming the basis ${\bf U}^{(2)}_{l}$, $l=1,\ldots,r^{n},$ vanishes
at the boundary. Combining (3.30)-(3.31) with the property of the conservative
truncation from Proposition 3.3, we have

$\displaystyle{\bm{\rho}}^{n+1}$
$\displaystyle={\bm{\rho}}^{n+1,*}={\bm{\rho}}^{n}-\Delta tD_{x}({\bf
J}^{n}),$ (3.33) $\displaystyle{\bf J}^{n+1}$ $\displaystyle={\bf
J}^{n+1,*}={\bf J}^{n}-\Delta t(D_{x}(\sigma^{n})-{\bm{\rho}}^{n}\star{\bf
E}^{n}),$ (3.34)

which are consistent and conservative discretization of the macroscopic moment
equations (2.7)-(2.8). ∎

######  Remark 3.5.

If a WENO finite difference discretization is employed for $D_{v}$, then the
local conservation of momentum cannot be guaranteed. This is because the WENO
discretization does not satisfy the summation by parts property (3.32). A WENO
finite difference discretization can be employed for $D_{x}$ for local
conservation.

######  Remark 3.6.

If only the mass conservation is desired, then we can modify the subspace
$\mathcal{N}=\text{span}\\{{\bf 1}_{v}\\}$ for projection; similarly, if only
the mass and momentum conservation is desired, we let
$\mathcal{N}=\text{span}\\{{\bf 1}_{v},{\bf v}\\}$. We denote the
corresponding projections as $P_{1}$ and $P_{2}$, respectively, and denote the
projection to $\mathcal{N}=\text{span}\\{{\bf 1}_{v},{\bf v},{\bf v}^{2}\\}$
as $P_{3}$. Numerical performances of different projection operators will be
assessed in the Section 5 for numerical experiments.

######  Remark 3.7.

The choice of weight function $w(v)$ can affect the performance. In
particular, if the weight function does not have sufficient decay, then the
solution may not be close enough to zero at the boundary of the velocity
domain, and large error could be incurred. On the other hand, when rescaling
$f_{2}$ by dividing ${w}$, small weights at the boundary can introduce
numerical instability.

Finally, the proposed conservative truncation algorithm is summarized in the
Algorithm 1.

  * •

Input: the pre-compressed low rank solution at time $t^{n+1}$:

${\bf f}^{n+1,*}=\sum_{l=1}^{r^{*}}C_{l}^{*}\ {\bf U}_{l}^{(1),*}\otimes{\bf
U}_{l}^{(2),*}.$

  * •

Output: the compressed low rank solution $f^{n+1}$ with the same density,
current density, and kinetic energy functions as $f^{n+1,*}$.

  1. 1.

Compute ${\bm{\rho}}^{n+1}$, ${\bf J}^{n+1}$, ${\bm{\kappa}}^{n+1}$ of ${\bf
f}^{n+1,*}$ from (3.22).

  2. 2.

Compute ${\mathcal{M}}^{n+1}$ by (3.27) from ${\bm{\rho}}^{n+1}$, ${\bf
J}^{n+1}$, ${\bm{\kappa}}^{n+1}$.

  3. 3.

Compute ${\bf f}_{1}$ from (3.26).

  4. 4.

Perform the truncation of ${\bf f}_{2}$ by (3.28).

  5. 5.

Update the compressed low rank solution by (3.29).

Algorithm 1 The conservative truncation procedure for the 1D1V VP system.

##  4 2D2V Vlasov-Poisson system by the HT format

We extend the proposed conservative algorithm to the 2D2V case by the HT
format. Below, we briefly review the fundamentals of the HT format for
efficiently representing tensors in $d$ dimensions, and the low rank tensor
method with the HT format for solving the 2D2V VP system (2.1).

$f_{t}+v_{1}f_{x_{1}}+v_{2}f_{x_{2}}+E_{1}f_{v_{1}}+E_{2}f_{v_{2}}=0,$ (4.1)

where the electric field $(E_{1},E_{2})$ is solved from the coupled Poisson’s
equation. The macroscopic equations can be obtained from taking moments of
(4.1) in the form of (2.7)-(2.9).

###  4.1 HT format for high order tensors

We denote the dimension index $D=\\{1,2,\ldots,d\\}$ and define a _dimension
tree_ $\mathcal{T}$ which is a binary tree containing a subset $\alpha\subset
D$ at each node. Furthermore, $\mathcal{T}$ has $D$ as the root node and
$\\{1\\},\ \\{2\\},\ldots,\ \\{d\\}$ as the leaf nodes. Each non-leaf node
$\alpha$ has two children nodes denoted as $\alpha_{l}$ and $\alpha_{r}$ with
$\alpha=\alpha_{l}\bigcup\alpha_{r}$ and
$\alpha_{l}\bigcap\alpha_{r}=\emptyset$. For example, the dimension tree
$\mathcal{T}$ given in Figure 4.1 can be used to approximate
$f(x_{1},x_{2},v_{1},v_{2})$ in (4.1) in the HT format. The efficiency of the
HT format lies in the nestedness property [19]: for a non-leaf node $\alpha$
with two children nodes $\alpha_{l},\,\alpha_{r}$, then

$\text{range}(\mathcal{M}^{(\alpha)}({\bf
a}))\subset\text{range}(\mathcal{M}^{(\alpha_{l})}({\bf
a})\otimes\mathcal{M}^{(\alpha_{r})}({\bf a})),$ (4.2)

which implies that there exists a third order tensor ${\bf
B}^{(\alpha)}\in\mathbb{R}^{r_{\alpha_{l}}\times r_{\alpha_{r}}\times
r_{\alpha}}$, known as the transfer tensor, such that

${\bf
U}_{l_{\alpha}}^{(\alpha)}=\sum_{l_{\alpha_{l}}=1}^{r_{\alpha_{l}}}\sum_{l_{\alpha_{r}=1}}^{r_{\alpha_{l}}}{\bf
B}^{(\alpha)}_{l_{\alpha_{l}},l_{\alpha_{r}},l_{\alpha}}{\bf
U}_{l_{\alpha_{l}}}^{(\alpha_{l})}\otimes{\bf
U}_{l_{\alpha_{r}}}^{(\alpha_{r})},\quad l_{\alpha}=1,\ldots,r_{\alpha}.$
(4.3)

In other words, the frame vectors at the parent node can be recovered by those
at the two children nodes $\alpha_{l},\,\alpha_{r}$ with the transfer tensor.
By recursively making use of (4.3), a tensor in the HT format stores a frame
at each leaf node and a third order transfer tensor at each non-leaf node
based on a dimension tree. Denote ${\bf
r}=\\{r_{\alpha}\\}_{\alpha\in\mathcal{T}}$ as the hierarchical ranks. The
storage of the HT format scales as $\mathcal{O}(dr^{3}+drN)$, where
$r=\max{\bf r}$. If $r$ is reasonably low, then the HT format avoids the curse
of dimensionality. In summary, the HT format is fully characterized by the
three key components, including a dimension tree, frames at leaf nodes and
transfer tensors at non-leaf nodes, see Figure 4.1 for the data layout.

$\\{1,\,2,\,3,\,4\\}$$\\{1,\,2\\}$$\\{1\\}$$\\{2\\}$$\\{3,\,4\\}$$\\{3\\}$$\\{4\\}$

${\bf B}^{(1,2,3,4)}$${\bf B}^{(1,2)}$${\bf U}^{(1)}$${\bf U}^{(2)}$${\bf
B}^{(3,4)}$${\bf U}^{(3)}$${\bf U}^{(4)}$

Figure 4.1: Dimension tree $\mathcal{T}$ to express fourth-order tensors in
the HT format.

###  4.2 A conservative low rank tensor method in HT for the 2D2V VP system

We follow the low rank tensor method for updating the 2D2V VP solution in
[17], but propose to perform a conservative hierarchical HOSVD truncation with
preservation of charge, current and kinetic energy density. We assume at each
time step, the solution ${\bf f}$ is expressed as the fourth order tensor in
the HT format with dimension tree $\mathcal{T}$ together with frames ${\bf
U}^{(1)},{\bf U}^{(2)},{\bf U}^{(3)},{\bf U}^{(4)}$ at four leaf nodes,
corresponding to directions $x_{1},\,x_{2},\,v_{1},\,v_{2}$, respectively, and
transfer tensors ${\bf B}^{(1,2,3,4)},{\bf B}^{(1,2)},{\bf B}^{(3,4)}$, see
Figure 4.1. In particular,

${\bf f}=\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf U}_{l_{12}}^{(1,2)}\otimes{\bf
U}_{l_{34}}^{(3,4)},$ (4.4)

with

${\bf U}_{l_{12}}^{(1,2)}=\sum_{l_{1}=1}^{r_{1}}{\sum_{l_{2}=1}^{r_{2}}}{\bf
B}^{(1,2)}_{l_{1},l_{2},l_{12}}{\bf U}_{l_{1}}^{(1)}\otimes{\bf
U}_{l_{2}}^{(2)},\quad l_{12}=1,\ldots,r_{12},$ (4.5)

and

${\bf U}_{l_{34}}^{(3,4)}=\sum_{l_{3}=1}^{r_{3}}{\sum_{l_{4}=1}^{r_{4}}}{\bf
B}^{(3,4)}_{l_{3},l_{4},l_{34}}{\bf U}_{l_{3}}^{(3)}\otimes{\bf
U}_{l_{4}}^{(4)},\quad l_{34}=1,\ldots,r_{34}.$ (4.6)

Further, the electric field $E_{1}$ and $E_{2}$ are represented in the second
order HT format.

Now we are ready to formulate a conservative HT truncation of VP solutions in
a low rank tensor format. We denote pre-compressed solution ${\bf f}$ in the
HT format as in (4.4). The idea of the conservative truncation is similar to
the 1D1V case. We first introduce a weighted inner product space for ${\bf
f},{\bf g}\in\mathbb{R}^{N_{v_{1}}\times N_{v_{2}}}$

$\langle{\bf f},{\bf g}\rangle_{\bf
w}=h_{v_{1}}h_{v_{2}}\sum_{j_{1}=1}^{N_{v_{1}}}\sum_{j_{2}=1}^{N_{v_{2}}}f_{j_{1},j_{2}}g_{j_{1},j_{2}}w_{j_{1},j_{2}},$
(4.7)

where $h_{v_{1}}$ and $h_{v_{2}}$ are the mesh size for the $v_{1}$ and
$v_{2}$ grids, respectively. Here

${\bf w}\doteq{\bf w}^{(1)}\otimes{\bf w}^{(2)}\in\mathbb{R}^{N_{v_{1}}\times
N_{v_{2}}},\quad{\bf w}^{(1)}\in\mathbb{R}^{N_{v_{1}}},{\bf
w}^{(2)}\in\mathbb{R}^{N_{v_{2}}},$ (4.8)

where ${\bf w}^{(1)}$, ${\bf w}^{(2)}$ are vectors as point values of the
weight function $w(v)=\exp(-\frac{{v}^{2}}{2})$ at corresponding $v_{1}$ and
$v_{2}$ grids. Then we seek a decomposition ${\bf f}={\bf f}_{1}+{\bf f}_{2}$,
where ${\bf f}_{1}$ comes from the rescaling and orthogonal projection of
${\bf f}$ with respect to the weighted inner product (4.7) onto the subspace

$\mathcal{N}\doteq\text{span}\\{{\bf 1}_{v_{1}\otimes v_{2}},{\bf
v}_{1}\otimes{\bf 1}_{v_{2}},{\bf 1}_{v_{1}}\otimes{\bf v}_{2},{\bf
v}_{1}^{2}\otimes{\bf 1}_{v_{2}}+{\bf 1}_{v_{1}}\otimes{\bf v}_{2}^{2}\\}$
(4.9)

for the conservation of mass, momentum and kinetic energy density, see
Proposition 4.1 below for details. For further discussion, we also introduce a
standard inner product for ${\bf f},{\bf g}\in\mathbb{R}^{N_{v_{1}}\times
N_{v_{2}}}$

$\langle{\bf f},{\bf
g}\rangle=h_{v_{1}}h_{v_{2}}\sum_{j_{1}=1}^{N_{v_{1}}}\sum_{j_{2}=1}^{N_{v_{2}}}f_{j_{1},j_{2}}g_{j_{1},j_{2}}.$
(4.10)

######  Proposition 4.1.

Let ${\bf f}_{1}$ come from the scaling and orthogonal projection of rescaled
${\bf f}$ with respect to the weighted inner product (4.7) onto the subspace
$\eqref{eq: 2d2v_subspa}$, in a similar spirit to the 1D1V case (3.25). Assume
${\bf f}$ is written in the low rank HT format (4.39). ${\bf f}_{1}$ can be
represented in low rank HT format (consistently with the subscript $1$ in the
notations),

$\displaystyle{\bf f}_{1}=\tilde{P}({\bf f})$ $\displaystyle\doteq$
$\displaystyle\sum_{l_{12}=1}^{4}{\sum_{l_{34}=1}^{4}}({\bf
B}_{1}^{(1,2,3,4)})_{l_{12},l_{34}}({\bf U}_{1}^{(1,2)})_{l_{12}}\otimes({\bf
U}_{1}^{(3,4)})_{l_{34}},$ (4.11)
$\displaystyle\stackrel{{\scriptstyle=}}{{\eqref{eq: B1}}}$
$\displaystyle\sum_{k=1}^{4}{({\bf U}_{1}^{(1,2)})_{k}}\otimes{({\bf
U}_{1}^{(3,4)})_{k}},$ (4.12)

where we introduce the notation of $\tilde{P}$ as the scaled orthogonal
projection. The specifications of the frame vectors and transfer tensors are
outlined below.

  * •

For ${\bf f}_{1}$, the hierarchical ranks are $r_{12}=r_{34}=4$,
$r_{3}=r_{4}=3$. $r_{1}$ and $r_{2}$ are the same as those for ${\bf f}$.

  * •

The transfer tensor, ${\bf B}_{1}^{(1,2,3,4)}$, of size $4\times 4$, is an
identity matrix,

${\bf B}_{1}^{(1,2,3,4)}={\bf I}_{4\times 4}.$ (4.13)

$({\bf U}_{1}^{(3,4)})_{k}$, $k=1,\ldots,4,$ is constructed from an
orthonormal set of basis $\\{{\bf V}_{1},\cdots{\bf V}_{4}\\}$ in the
$v_{1}-v_{2}$ dimensions defined in (4.40). We have $r_{3}=r_{4}=3$ and

$({\bf U}_{1}^{(3,4)})_{k}=\sum_{l_{3}=1}^{3}{\sum_{l_{4}=1}^{3}}({\bf
B}_{1}^{(3,4)})_{l_{3},l_{4},k}({\bf w}^{(1)}\star({\bf
U}_{1}^{(3)})_{l_{3}})\otimes({\bf w}^{(2)}\star({\bf U}_{1}^{(4)})_{l_{4}})$
(4.14)

with ${\bf w}^{(1)}$ and ${\bf w}^{(2)}$ the same as in (4.8). The frame
vectors for node $3$ are

$\displaystyle({\bf U}_{1}^{(3)})_{1}=\frac{1}{c_{1}}{\bf
1}_{v_{1}},\quad({\bf U}_{1}^{(3)})_{2}=\frac{1}{c_{2}}{\bf v}_{1},\quad({\bf
U}_{1}^{(3)})_{3}=\frac{1}{c_{3}}({\bf v}_{1}^{2}-c{\bf 1}_{v_{1}}),$ (4.15)

where $c=\frac{\langle{\mathbf{1}_{v_{1}},{\bf v}_{1}^{2}}\rangle_{\bf
w}}{\langle{\mathbf{1}_{v_{1}},\mathbf{1}_{v_{1}}}\rangle_{\bf w}}$ is the
orthogonalization constant of the basis, and $c_{l}$, $l=1,2,3$ are
normalization constants for the corresponding basis of ${\bf 1}_{v_{1}}$,
${\bf v}_{1}$ and ${\bf v}_{1}^{2}-c{\bf 1}_{v_{1}}$. We have the same frame
vectors for the node $4$ but for $v_{2}$, assuming that the weight function
and domain in $v_{2}$ is the same as $v_{1}$,

$\displaystyle({\bf U}_{1}^{(4)})_{1}=\frac{1}{c_{1}}{\bf
1}_{v_{2}},\quad({\bf U}_{1}^{(4)})_{2}=\frac{1}{c_{2}}{\bf v}_{2},\quad({\bf
U}_{1}^{(4)})_{3}=\frac{1}{c_{3}}({\bf v}_{2}^{2}-c{\bf 1}_{v_{2}})$ (4.16)

The transfer tensor ${\bf B}_{1}^{(3,4)}$ is a tensor of size $3\times 3\times
4$. It has zero elements, except the following specification for $({\bf
B}_{1}^{(3,4)})_{l_{3},l_{4},l_{34}}$

$({\bf B}_{1}^{(3,4)})_{1,1,1}=({\bf B}_{1}^{(3,4)})_{2,1,2}=({\bf
B}_{1}^{(3,4)})_{1,2,3}=1,\quad({\bf B}_{1}^{(3,4)})_{3,1,4}=({\bf
B}_{1}^{(3,4)})_{1,3,4}=\frac{1}{\sqrt{2}}.$ (4.17)

  * •

$({\bf U}_{1}^{(1,2)})_{k}$, $k=1,\cdots 4$, is in the following form,

$({\bf U}_{1}^{(1,2)})_{k}=\sum_{l_{1}=1}^{r_{1}}{\sum_{l_{2}=1}^{r_{2}}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},k}({\bf U}_{1}^{(1)})_{l_{1}}\otimes({\bf
U}_{1}^{(2)})_{l_{2}},$ (4.18)

with the same frame vectors as ${\bf f}$ on the leaf nodes $1$ and $2$,
meaning that ${\bf U}_{1}^{(1)}$ and ${\bf U}_{1}^{(2)}$ are the same as ${\bf
U}^{(1)}$ and ${\bf U}^{(2)}$, respectively. ${\bf B}_{1}^{(1,2)}$, of size
$r_{1}\times r_{2}\times 4$, has its elements as,

$({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},k}=\sum_{l_{12}=1}^{r_{12}}\sum_{l_{34}=1}^{r_{34}}{\bf
B}^{(1,2)}_{l_{1},l_{2},l_{12}}{\bf S}_{k,l_{34}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}},\quad k=1,\cdots 4,$ (4.19)

with $r_{12}$, $r_{34}$ and transfer tensors ${\bf B}^{(1,2)}$ and ${\bf
B}^{(1,2,3,4)}$ from the HT representation of ${\bf f}$, and

${\bf S}_{k,l_{34}}=\langle{\bf U}_{l_{34}}^{(3,4)},{\bf V}_{k}\rangle,\quad
k=1,\cdots 4,$ (4.20)

where ${\bf U}_{l_{34}}^{(3,4)}$, $l_{34}=1,\ldots,r_{34}$, is the frame
tensor for node (3, 4) of ${\bf f}$. The inner product in the sense of (4.10)
can be evaluated in a dimension-by-dimension manner.

Finally ${\bf f}_{1}$ has the same discrete macroscopic charge, current and
kinetic energy density as with $\bf f$, which are denoted as $\bm{\rho}$,
${\bf J}_{1}$, ${\bf J}_{2}$, and $\bm{\kappa}$. The discrete macroscopic
charge, current and kinetic energy density of ${\bf f}_{1}$ are

${\bm{\rho}}=(c_{1})^{2}({\bf
U}_{1}^{(1,2)})_{1}=(c_{1})^{2}\sum_{l_{1},l_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},1}{\bf U}^{(1)}_{l_{1}}\otimes{\bf
U}^{(2)}_{l_{2}},$ (4.21) ${{\bf J}_{1}}=c_{1}c_{2}({\bf
U}_{1}^{(1,2)})_{2}=c_{1}c_{2}\sum_{l_{1},l_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},2}{\bf U}^{(1)}_{l_{1}}\otimes{\bf
U}^{(2)}_{l_{2}},$ (4.22) ${{\bf J}_{2}}=c_{1}c_{2}({\bf
U}_{1}^{(1,2)})_{3}=c_{1}c_{2}\sum_{l_{1},l_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},3}{\bf U}^{(1)}_{l_{1}}\otimes{\bf
U}^{(2)}_{l_{2}},$ (4.23) ${\bm{\kappa}}=\frac{1}{\sqrt{2}}c_{1}c_{3}({\bf
U}_{1}^{(1,2)})_{4}+c\bm{\rho}=\frac{1}{\sqrt{2}}c_{1}c_{3}\sum_{l_{1},l_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},4}{\bf U}^{(1)}_{l_{1}}\otimes{\bf
U}^{(2)}_{l_{2}}+c\bm{\rho}.$ (4.24)

They are the same as those of ${\bf f}$, which can be obtained by the HT
tensor contraction

$\displaystyle\left(\begin{array}[]{l}{\bm{\rho}}\\\ {\bf J}_{1}\\\ {\bf
J}_{2}\\\ {\bm{\kappa}}\end{array}\right)$
$\displaystyle=\sum_{l_{12}}\sum_{l_{34}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}\left\langle{\bf
U}_{l_{34}}^{(3,4)},\left(\begin{array}[]{l}{\bf 1}_{v_{1}\otimes v_{2}}\\\
{\bf v}_{1}\otimes{\bf 1}_{v_{2}}\\\ {\bf 1}_{v_{1}}\otimes{\bf v}_{2}\\\
\frac{1}{2}{\bf v}_{1}^{2}\otimes{\bf 1}_{v_{2}}+\frac{1}{2}{\bf
1}_{v_{1}}\otimes{\bf v}_{2}^{2}\end{array}\right)\right\rangle{\bf
U}_{l_{12}}^{(1,2)},$ (4.33) $\displaystyle=\sum_{l_{12}}\sum_{l_{34}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}\left(\begin{array}[]{l}S_{1,l_{34}}(c_{1})^{2}\\\
S_{2,l_{34}}c_{1}c_{2}\\\ S_{3,l_{34}}c_{1}c_{2}\\\
S_{4,l_{34}}c_{1}c_{3}\frac{1}{\sqrt{2}}+S_{1,l_{34}}(c_{1})^{2}c\end{array}\right){\bf
U}_{l_{12}}^{(1,2)}.$ (4.38)

###### Proof.

We first scale ${\bf f}$ by ${\bf w}$ in (4.8)

$\tilde{\bf f}=\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf U}_{l_{12}}^{(1,2)}\otimes\frac{1}{\bf
w}\star{\bf U}_{l_{34}}^{(3,4)}.$ (4.39)

Let $\tilde{{\bf U}}_{l_{34}}^{(3,4)}\doteq\frac{1}{\bf w}\star{\bf
U}_{l_{34}}^{(3,4)}$, where $\star$ is the element-wise multiplication in the
corresponding dimensions. Then we perform the orthogonal projection of
$\tilde{\bf f}$ onto the subspace of (4.9) to obtain $\tilde{\bf f}_{1}$.
Finally we rescale back to ${\bf f}_{1}={\bf w}\star\tilde{\bf f}_{1}$.

To perform the orthogonal projection, the orthogonalization of basis in (4.9)
gives the frame vectors for the node 3 and 4 as specified in (4.15) and the
transfer tensor ${\bf B}_{1}^{(3,4)}$ as specified in (4.17). This gives a set
of orthonormal basis of $({\bf U}_{1}^{(3,4)})_{l}$, $l=1,2,3,4$. We denote
them as ${\bf V}_{l}$, $l=1,2,3,4$,

$\displaystyle{\bf V}_{1}$ $\displaystyle=({\bf U}_{1}^{(3)})_{1}\otimes({\bf
U}_{1}^{(4)})_{1},\quad{\bf V}_{2}=({\bf U}_{1}^{(3)})_{2}\otimes({\bf
U}_{1}^{(4)})_{1}$ $\displaystyle{\bf V}_{3}$ $\displaystyle=({\bf
U}_{1}^{(3)})_{1}\otimes({\bf U}_{1}^{(4)})_{2},\quad{\bf
V}_{4}=\frac{1}{\sqrt{2}}({\bf U}_{1}^{(3)})_{3}\otimes({\bf
U}_{1}^{(4)})_{1}+\frac{1}{\sqrt{2}}({\bf U}_{1}^{(3)})_{1}\otimes({\bf
U}_{1}^{(4)})_{3}.$ (4.40)

To perform the orthogonal projection to (4.9), only the bases of node $(3,4)$
are affected. In particular, let ${\bf P}^{(3,4)}$ be the projection operator
at the node $(3,4)$, then

$\displaystyle{\bf P}^{(3,4)}\tilde{\bf f}$ $\displaystyle={\bf
P}^{(3,4)}\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf U}_{l_{12}}^{(1,2)}\otimes\frac{1}{\bf
w}\star{\bf U}_{l_{34}}^{(3,4)}$
$\displaystyle=\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf U}_{l_{12}}^{(1,2)}\otimes\left({\bf
P}^{(3,4)}(\frac{1}{\bf w}\star{\bf U}_{l_{34}}^{(3,4)})\right),$ (4.41)

where

${\bf P}^{(3,4)}(\frac{1}{\bf w}\star{\bf U}_{l_{34}}^{(3,4)})={\bf
P}^{(3,4)}\tilde{{\bf U}}^{(3,4)}_{l_{34}}=\sum_{k=1}^{4}{\langle\tilde{{\bf
U}}^{(3,4)}_{l_{34}},{\bf V}_{k}\rangle_{{\bf w}}}{\bf
V}_{k}\doteq\sum_{k=1}^{4}{\bf S}_{k,l_{34}}{\bf V}_{k},$ (4.42)

with

${\bf S}_{k,l_{34}}={\langle\tilde{{\bf U}}^{(3,4)}_{l_{34}},{\bf
V}_{k}\rangle_{{\bf w}}}=\langle\frac{1}{\bf w}\star{\bf
U}_{l_{34}}^{(3,4)},{\bf V}_{k}\rangle_{{\bf w}}=\langle{\bf
U}_{l_{34}}^{(3,4)},{\bf V}_{k}\rangle.$ (4.43)

Plug (4.42) into (4.41),

$\displaystyle{\bf P}^{(3,4)}\tilde{\bf f}$
$\displaystyle=\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf
U}_{l_{12}}^{(1,2)}\otimes(\sum_{k=1}^{4}{\bf S}_{k,l_{34}}{\bf V}_{k}),$
$\displaystyle=\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}\left(\sum_{l_{1}=1}^{r_{1}}\sum_{l_{2}=1}^{r_{2}}({\bf
U}_{l_{1}}^{(1)}\otimes{\bf U}_{l_{2}}^{(2)}){\bf
B}^{(1,2)}_{l_{1},l_{2},l_{12}}\right)\otimes(\sum_{k=1}^{4}{\bf
S}_{k,l_{34}}{\bf V}_{k}),$
$\displaystyle=\sum_{k=1}^{4}\left(\sum_{l_{1}=1}^{r_{1}}\sum_{l_{2}=1}^{r_{2}}\underbrace{\left(\sum_{l_{12}=1}^{r_{12}}{\sum_{l_{34}=1}^{r_{34}}}{\bf
B}^{(1,2,3,4)}_{l_{12},l_{34}}{\bf B}^{(1,2)}_{l_{1},l_{2},l_{12}}{\bf
S}_{k,l_{34}}\right)}_{\doteq({\bf B}_{1}^{(1,2)})_{l_{1},l_{2},k}}({\bf
U}_{l_{1}}^{(1)}\otimes{\bf U}_{l_{2}}^{(2)})\right)\otimes{\bf V}_{k},$
(4.44)
$\displaystyle=\sum_{k=1}^{4}\left(\sum_{l_{1}=1}^{r_{1}}\sum_{l_{2}=1}^{r_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},k}({\bf U}_{l_{1}}^{(1)}\otimes{\bf
U}_{l_{2}}^{(2)})\right)\otimes{\bf V}_{k},$ (4.45)

where we let $({\bf B}_{1}^{(1,2)})_{l_{1},l_{2},k}$ as specified in (4.19),
and

$({\bf
U}_{1}^{(1,2)})_{k}\doteq\left(\sum_{l_{1}=1}^{r_{1}}\sum_{l_{2}=1}^{r_{2}}({\bf
B}_{1}^{(1,2)})_{l_{1},l_{2},k}({\bf U}_{l_{1}}^{(1)}\otimes{\bf
U}_{l_{2}}^{(2)})\right)$ (4.46)

as in (4.18).

The macroscopic charge, momentum and energy density of ${\bf f}_{1}$,
(4.21)-(4.24) can be derived, from the form of ${\bf f}_{1}$ in (4.12), the
form of ${\bf U}_{1}^{(1,2)}$ in (4.46), and

${\bf 1}_{v_{1}\otimes v_{2}}=(c_{1})^{2}{\bf V}_{1},\quad{\bf
v}_{1}\otimes{\bf 1}_{v_{2}}=c_{1}c_{2}{\bf V}_{2},\quad{\bf
1}_{v_{1}}\otimes{\bf v}_{2}=c_{1}c_{2}{\bf V}_{3},$ $\frac{1}{2}{\bf
v}_{1}^{2}\otimes{\bf 1}_{v_{2}}+\frac{1}{2}{\bf 1}_{v_{1}}\otimes{\bf
v}_{2}^{2}=\frac{1}{\sqrt{2}}c_{1}c_{3}{\bf V}_{4}+c(c_{1})^{2}{\bf V}_{1},$

which is due to (4.40). The agreement of macroscopic charge, momentum and
energy density of ${\bf f}$ and ${\bf f}_{1}$ is a direct consequence of
manipulation of equalities in (4.18), (4.19), (4.38) and orthonormal property
of the basis in ${\bf U}_{1}^{(3,4)}$.

${\bf B}_{1}^{(1,2,3,4)}$ in (4.13)${\bf B}_{1}^{(1,2)}$ in (4.19)${\bf
U}^{(1)}$${\bf U}^{(2)}$${\bf B}_{1}^{(3,4)}$ in (4.17)${\bf U}_{1}^{(3)}$ in
(4.15)${\bf U}_{1}^{(4)}$ in (4.16) Figure 4.2: The data layout of ${\bf
f}_{1}$. Here ${\bf U}^{(1)}$ and ${\bf U}^{(2)}$ are the same as those vector
frames for ${\bf f}$.

∎

Once the orthogonal projection is performed, ${\bf f}_{1}$ stays untouched for
the truncation. As with the 1D1V case, we will perform HOSVD to truncate the
remainder ${\bf f}_{2}\doteq{\bf f}-{\bf f}_{1}$. If a standard HT truncation
is directly applied to ${\bf f}_{2}$, unlike the 1D1V case the conservation
property cannot be guaranteed. Below, we elaborate and investigate such an
issue, and further develop a conservative projection procedure after the HT
truncation to ensure charge, current and kinetic energy density conservation.

We start with a brief description of a naive weighted hierarchical HOSVD
truncation procedure, as a direct analog of the 1D1V case. First, we scale
${\bf f}_{2}$ according to the weights and define $\tilde{\bf
f}_{2}=\frac{1}{\sqrt{\bf w}}\star{\bf f}_{2}$. Note that the rescaling is
computed dimension-by-dimension. The standard hierarchical HOSVD root-to-leaf
truncation with threshold $\varepsilon$ is applied to $\tilde{\bf f}_{2}$ and
${\mathcal{T}}_{\varepsilon}(\tilde{\bf f}_{2})$ is obtained. Finally, ${\bf
f}_{2}$ is defined as

${\bf f}_{2}={\sqrt{\bf w}}\star{\mathcal{T}}_{\varepsilon}(\tilde{\bf
f}_{2}).$ (4.47)

The issue with such procedure is the loss of macroscopic conservation in the
root-to-leaf truncation. In particular, the HT truncation
${\mathcal{T}}_{\varepsilon}(\tilde{\bf f}_{2})$ can be represented as, see
[14],

${\mathcal{T}}_{\varepsilon}(\tilde{\bf
f}_{2})=(\pi_{1}\otimes\pi_{2}\otimes\pi_{3}\otimes\pi_{4})(\pi_{12}\otimes\pi_{34})\,{\bf
f}_{2},$ (4.48)

where $\pi_{\alpha}$ denotes the orthogonal projection on the subspace spanned
by $\tilde{r}_{\alpha}$ leading left singular vectors $\\{\tilde{{\bf
U}}^{(\alpha)}_{i}\\}_{i=1}^{\tilde{r}_{\alpha}}$ of the matricization
${\mathcal{M}}^{(\alpha)}({\bf f}_{2})$. However, the bases at node $(3,4)$
are no longer orthogonal to $\\{{\bf V}_{1},\,{\bf V}_{2},\,{\bf V}_{3},\,{\bf
V}_{4}\\}$ in (4.40) due to the truncation $\pi_{3}$ and $\pi_{4}$ at leaf
nodes. Hence ${\bf f}_{2}$ as in (4.47) may have nonzero charge, current and
kinetic energy density, breaking the conservation property. To fix the issue,
we propose to apply the operator $({\bf I}-\tilde{P})$, to ${\sqrt{\bf
w}}\star{\mathcal{T}}_{\varepsilon}(\tilde{\bf f}_{2})$ to ensure zero charge,
current and kinetic energy density of truncated ${\bf f}_{2}$. Here
$\tilde{P}$ is the same as that specified in (4.11) in Proposition 4.1. We
introduce the following notation,

$\widetilde{{\mathcal{T}}_{\varepsilon}}({\bf f}_{2})\doteq({\bf
I}-\tilde{P})\left({\sqrt{\bf w}}\star{\mathcal{T}}_{\varepsilon}(\tilde{\bf
f}_{2})\right).$ (4.49)

Finally, the conservative HT truncation of ${\bf f}$ is done as follows

$T_{c}({\bf f})\doteq{\bf f}_{1}+\widetilde{{\mathcal{T}}_{\varepsilon}}({\bf
f}_{2}).$ (4.50)

The following proposition is a straightforward consequence from the orthogonal
projection.

######  Proposition 4.2.

$\widetilde{{\mathcal{T}}_{\varepsilon}}({\bf f}_{2})$ has zero charge
density, zero current density, and zero kinetic energy density. Hence,
$T_{c}({\bf f})$ preserves the charge, current, and kinetic energy densities
($\bm{\rho}$, ${\bf J}_{1}$, ${\bf J}_{2}$, $\bm{\kappa}$) of the original
${\bf f}$.

###### Proof.

The zero charge, current and kinetic energy density is a direct consequence of
the ${\bf I}-\tilde{P}$ projection operator. From this fact, together with
from Proposition 4.1, $T_{c}({\bf f})$ preserves the charge, current, and
kinetic energy densities ($\bm{\rho}$, ${\bf J}_{1}$, ${\bf J}_{2}$,
$\bm{\kappa}$) of the original ${\bf f}$. ∎

######  Proposition 4.3.

(Local mass and momentum conservation for the 2D2V VP system.) If the discrete
differential operators $D_{x},\,D_{v}$ employed are conservative, i.e., can be
written in a flux difference form, and linear, i.e., can preserve linear
relations, then the proposed low rank method with an SSP multi-step time
integrator preserves the mass and momentum locally; that is the schemes for
$\bm{\rho}$, ${\bf J}_{1}$, ${\bf J}_{2}$, from integrating the scheme on
${\bf f}^{n}$, are consistent and conservative discretization of the
macroscopic moment equations (2.7)-(2.8).

###### Proof.

The proof is based on a conservative and linear discretization of the discrete
differential operators $D_{x},\,D_{v}$ employed and from the conservative
truncation procedure as proposed above. The details are similar to that of
Proposition 3.4 and hence omitted for brevity. ∎

######  Remark 4.4.

To ensure the local conservation, ${\bf f}_{1}$ should not be further
compressed. In the numerical simulation, we still truncate ${\bf f}_{1}$ at
nodes (1, 2) with threshold $10^{-15}$ to remove the redundancy from the add
basis procedure. Thus the local conservation property is preserved on the same
scale of machine precision, i.e. $10^{-15}$.

We summarized the conservative truncation procedure as Algorithm 2 below for
the 2D2V VP solution.

  * •

Input: the pre-compressed low rank solution ${\bf f}$ in the HT format with
dimension tree given in Figure 4.1 (a) and the associated data including the
frame tensors ${\bf U}^{(1)}$, ${\bf U}^{(2)}$, ${\bf U}^{(3)}$, and ${\bf
U}^{(4)}$ at nodes $(1)$, $(2)$, $(3)$, and $(4)$, respectively, and transfer
tensors ${\bf B}^{(1,2)}$, ${\bf B}^{(3,4)}$, and ${\bf B}^{(1,2,3,4)}$ at
nodes $(1,2)$, $(3,4)$, and $(1,2,3,4)$, respectively.

  * •

Output: the compressed low rank solution $T_{c}({\bf f})$ in the HT format
with the same charge, current, and kinetic energy density functions as ${\bf
f}$.

  1. 1.

Compute the rescaled orthogonal projection to obtain ${\bf
f}_{1}=\tilde{P}({\bf f})$, in the HT format with data layout in Figure 4.2.
The transfer tensor ${\bf B}_{1}^{(1,2,3,4)}$ is the identity matrix of size
$4\times 4$, $({\bf B}_{1}^{(1,2)})_{l_{12}}$ is a matrix of size $r_{1}\times
r_{2}\times 4$ from (4.19), and ${\bf B}_{1}^{(3,4)}$ is a matrix of size
$3\times 3\times 4$ from (4.17), with the frame tensors ${\bf
U}_{1}^{(1)}={\bf U}^{(1)}$, ${\bf U}_{1}^{(2)}={\bf U}^{(2)}$, ${\bf
U}_{1}^{(3)}$ and ${\bf U}_{1}^{(3)}$ from (4.15) and (4.16) respectively. At
the end of this step, we truncate node (1, 2) of ${\bf f}_{1}$ with threshold
$10^{-15}$ to remove redundant basis in ${\bf U}_{1}^{(1,2)}$.

  2. 2.

Perform the HOSVD truncation, together with an orthogonal projection operator,
to ${\bf f}_{2}\doteq{\bf f}-{\bf f}_{1}$ to ensure zero charge, current and
energy densities:

    1. (a)

Compute ${\bf f}_{2}={\bf f}-{\bf f}_{1}$, and scale it to obtain $\tilde{\bf
f}_{2}=\frac{1}{\sqrt{\bf w}}\star{\bf f}_{2}$.

    2. (b)

Apply the standard HOSVD truncation to $\tilde{\bf f}_{2}$, and apply
rescaling to obtain ${\sqrt{\bf w}}\star{\mathcal{T}}_{\varepsilon}(\tilde{\bf
f}_{2})$.

    3. (c)

Apply ${\bf I}-\tilde{P}$ to ${\sqrt{\bf
w}}\star{\mathcal{T}}_{\varepsilon}(\tilde{\bf f}_{2})$ to obtain
$\widetilde{{\mathcal{T}}_{\varepsilon}}({\bf f}_{2})$, i.e. (4.49), with the
same $\tilde{P}$ operator as in the previous step.

  3. 3.

Update the compressed low rank solution $T_{c}({\bf f})={\bf
f}_{1}+\widetilde{{\mathcal{T}}_{\varepsilon}}({\bf f}_{2})$ from (4.50).

Algorithm 2 The conservative truncation procedure for the 2D2V VP solution.

##  5 Numerical results

In this section we present a collection of numerical examples to demonstrate
the efficacy of the proposed conservative low rank tensor method for
simulating the VP system. In the simulations, fifth order upwind finite
difference methods are employed for spatial discretization, together with a
second order SSP multi-step method denoted by SSPML2 for temporal
discretization. The numerical solutions of high dimensions are represented in
the HT format [18]. We also compare the proposed conservative low rank methods
against the nonconservative version in terms of efficiency and ability to
conserve the physical invariants. Unless otherwise noted, we let the weight
function $w({\bf v})=\exp(-\frac{|{\bf v}|^{2}}{2})$.

###  5.1 1D1V Vlasov-Poisson system

######  Example 5.1.

(Weak Laudau damping.) We consider the weak Landau damping test with initial
condition

$f(x,v,t=0)=\frac{1}{\sqrt{2\pi}}\left(1+\alpha\cos\left(kx\right)\right)\exp\left(-\frac{v^{2}}{2}\right),$
(5.1)

where $\alpha=0.01$ and $k=0.5$.

The computational domain is set to be $[0,L_{x}]\times[-L_{v},L_{v}]$ with
$L_{x}=2\pi/k$ and $L_{v}=6$. We set $\varepsilon=10^{-5}$ for truncation. In
Figure 5.3, we report the time histories of the electric energy and numerical
ranks of the solutions computed by the conservative and non-conservative
methods for comparision. It is observed that both methods are able to predict
the correct damping rate of the electric energy, and meanwhile, the numerical
ranks of the conservative method are bigger than those of the non-conservative
counterpart. In Figure 5.4, the time histories of the relative deviation of
total mass, momentum and energy are plotted. The conservative method is found
to be able to conserve the total mass and momentum up to machine precision
regardless of the mesh size used. As mentioned above, the proposed
conservative method cannot conserve the total energy, as the time integrator
employed is not energy conserving. With the mesh refinement, the conservarion
error of total energy decreases. Meanwhile, the non-conservative method can
conserve the total mass and energy up to the magnitude of truncation threshold
$\varepsilon$, but the total momentum is conserved to the machine precision,
which is attributed to the symmetry of the solution in the velocity direction.
Noteworthy, though both methods cannot conserve the total energy, the
conservative one does a better job in energy conservation compared to the non-
conservative counterpart, and it is because that the kinetic energy is
preserved in the truncation.

Figure 5.3: Example 5.1. The time evolution of the electric energy (a, c) and
the rank of the numerical solutions (b, d). Conservative method (a, b) and
non-conservative method (c, d). $\varepsilon=10^{-5}$.

Figure 5.4: Example 5.1. The time evolution of relative deviation of total
mass (a, d), absolute total momentum (b, d), and relative deviation of total
energy (c, f). Conservative method (a, b, c) and non-conservative method (d,
e, f). $\varepsilon=10^{-5}$.

######  Example 5.2.

(Strong Laudau damping.) We consider the strong Landau damping test with the
initial condition (5.1) and a bigger perturbation parameter $\alpha=0.5$.

Two truncation thresholds $\varepsilon=10^{-3},\,10^{-4}$ are used to compare
the performance of the proposed conservative method with the non-conservative
one. In Figure 5.2, we report the time evolution of the electric energy
together with the ranks of the numerical solutions for $\varepsilon=10^{-3}$.
We observe that the conservative method is able to capture correctly the
nonlinear dynamics of the strong Landau damping as opposed to the non-
conservative method. This is because the truncation error due to the large
threshold used greatly pollutes the accuracy for the non-conservative method,
while by design the conservative method exactly conserves the mass and
momentum densities in the low rank setting. Such conservation help resolve the
nonlinear Vlasov dynamics with a relatively large truncation threshold. As
observed in Figure 5.6, the conservative method can conserve the total mass
and momentum up to the machine precision regardless of the mesh size used.
Again, the energy conservation is not observed, but the conservation error
decreases with mesh refinement, which is not the case for the non-conservative
method. Then we consider a smaller truncation threshold $\varepsilon=10^{-4}$,
and the truncation error is reduced accordingly. Both methods generate
consitent results as plotted in Figure 5.7. As with the weak case, the
conservative method has slightly larger ranks than the non-conservative
method. We have a similar observation of the methods in conserving the
invariants in Figure 5.8, as that in Figure 5.6.

Figure 5.5: Example 5.2. The time evolution of the electric energy (a, c) and
the rank of the numerical solutions (b, d). Conservative method (a, b) and
non-conservative method (c, d). $\varepsilon=10^{-3}$.

Figure 5.6: Example 5.2. The time evolution of relative deviation of total
mass (a, d), absolute total momentum (b, e), and relative deviation of total
energy (c, f). Conservative method (a, b, c) and non-conservative method (d,
e, f). $\varepsilon=10^{-3}$.

Figure 5.7: Example 5.2. The time evolution of the electric energy (a, c) and
the ranks of the numerical solutions (b, d). Conservative method (a, b) and
non-conservative method (c, d). $\varepsilon=10^{-4}$.

Figure 5.8: Example 5.2. The time evolution of relative deviation of total
mass (a, d), absolute total momentum (b, e), and relative deviation of total
energy (c, f). Conservative method (a, b, c) and non-conservative method (d,
e, f). $\varepsilon=10^{-4}$.

######  Example 5.3.

(Bump on tail.) In this example, we simulate the bump-on-tail test with the
initial condition

$f(x,v,t=0)=\left(1+\alpha\cos\left(kx\right)\right)\left(n_{p}\exp\left(-\frac{v^{2}}{2}\right)+n_{b}\exp\left(-\frac{(v-u)^{2}}{2v_{t}}\right)\right),$
(5.2)

where $\alpha=0.04$, $k=0.3$, $n_{p}=\frac{9}{10\sqrt{2\pi}}$,
$n_{b}=\frac{2}{10\sqrt{2\pi}}$, $u=4.5$, $v_{t}=0.5$.

In the simulation, we set $\varepsilon=10^{-4}$ for truncation. The weight
function $w(v)=\exp(-\frac{v^{2}}{3})$ is chosen for the conservative method.
Note that unlike the previous Landau damping examples, the solution of bump on
tail does not have the symmetry. In Figure 5.9, we report the time evolution
of the electric energy as well as the ranks of the solutions. Consistent
numerical results are observed. In Figure 5.10, we plot the time evolution of
the relative deviation of the total mass, momentum and energy for both
methods. The non-conservative method is able to conserve the invariants on the
scale of the truncation threshold. Note that the exact momentum conservation
is not observed as expected. Meanwhile, the conservative method is able to
conserve the total mass and momentum up to the machine precision and has
smaller conservation errors of total energy compared to the non-conservative
method. Note that the conservation error of the total momentum increases after
$t=15$ for the conservative method with the coarse mesh size $32\times 64$
which is ascribed to the boundary error. In Figure 5.11, we report the contour
plots of the solutions by the two methods. The numerical solutions are
observed to qualitatively match each other. Last, we test the performance of
projector $P_{\mathcal{N}}$ with different subspaces $\mathcal{N}$ as
discussed in Remark 3.6. In particular, we denote by $P_{1}$, $P_{2}$ and
$P_{3}$ the orthogonal projectors with the weighted inner product (3.15) onto
the subspaces $\text{span}\\{1_{v}\\}$, $\text{span}\\{1_{v},{\bf v}\\}$,
$\text{span}\\{1_{v},{\bf v},{\bf v}^{2}\\}$, respectively. In Figure 5.12, we
plot the time evolution of the electric energy, ranks, total mass, total
momentum, and total energy of the solutions from the conservative method with
$P_{1}$, $P_{2}$ and $P_{3}$. The three methods generate consistent results
for the time evolution of the electric energy, and the ranks are comparable.
By construction all three methods can conserve the total mass, and the methods
with $P_{2}$ and $P_{3}$ can conserve the total momentum. Further, the method
with $P_{3}$ does the best job in conserving the total energy as $P_{3}$
preserves the kinetic energy for truncation.

Figure 5.9: Example 5.3. The time evolution of the electric energy (a, c) and
the ranks of the numerical solutions (b, d). Conservative method (a, b) and
non-conservative method (c, d). $\varepsilon=10^{-4}$.

Figure 5.10: Example 5.3. The time evolution of relative deviation of total
mass (a, d), relative deviation of total momentum (b, e), and relative
deviation of total energy (c, f). Conservative method (a, b, c) and non-
conservative method (d, e, f). $\varepsilon=10^{-4}$.

Figure 5.11: Example 5.3. Contour plots of the solutions at $t=30$.
Conservative method with mesh $N_{x}\times N_{v}=64\times 128$ (a). Non-
conservative method with mesh $N_{x}\times N_{v}=64\times 128$ (b).
Conservative method with mesh $N_{x}\times N_{v}=128\times 256$ (c). Non-
conservative method with mesh $N_{x}\times N_{v}=128\times 256$ (d).
$\varepsilon=10^{-4}$.

Figure 5.12: Example 5.3. Comparison of projections $P_{1},\,P_{2},\,P_{3}$ of
the proposed conservative low rank method. The time evolution of electric
energy (a), ranks of the numerical solutions (b), relative deviation of total
mass (c), relative deviation of total momentum (e), and relative deviation of
total energy (f). $N_{x}\times N_{v}=64\times 128$. $\varepsilon=10^{-4}$.

####  5.1.1 2D2V Vlasov-Poisson system

######  Example 5.4.

(Weak Landau damping.) We consider the 2D2V weak Landau damping, the dynamics
of which is similar to the 1D1V case. The initial condition is

$f({\bf x},{\bf
v},t=0)=\frac{1}{(2\pi)^{d/2}}\left(1+\alpha\sum_{m=1}^{d}\cos\left(kx_{m}\right)\right)\exp\left(-\frac{|{\bf
v}|^{2}}{2}\right),$ (5.3)

where $d=2$, $\alpha=0.01$, and $k=0.5$.

We set the computation domain as $[0,L_{x}]^{2}\times[-L_{v},L_{v}]^{2}$,
where $L_{x}=\frac{2\pi}{k}$ and $L_{v}=6$, and the truncation threshold
$\varepsilon=10^{-5}$. We simulate the problem with both conservative and non-
conservative methods, and the solutions are represented in the fourth order HT
format. In Figures 5.13-5.14, we report the time evolution of the electric
energy, hierarchical ranks of the numerical solution, relative deviation of
total mass and energy together with absolute total momentum $J_{1}$ and
$J_{2}$. It is observed that both methods are able to predict the damping rate
of the electric energy. Furthermore, the conservative method is able to
conserve the total mass and momentum $J_{1}$ and $J_{2}$ up to the machine
precision and enjoys better total energy conservation compared to the non-
conservative counterpart. On the other hand, the hierarchical ranks of the
solution from the conservative method, especially $r_{1}$ and $r_{2}$, are
larger than that from the non-conservative method. This is because ${\bf
f}_{1}$ is constructed without compression to guarantee the local conservation
(in fact ${\bf f}_{1}$ is truncated with threshold $10^{-15}$ in the
simulation), and then the solution tensor is not compressed in the ${\bf x}$
direction.

Figure 5.13: Example 5.4. Conservative low rank method. The time evolution of
electric energy (a), hierarchical ranks of the xnumerical solution of mesh
size $N_{x}\times N_{v}=64\times 128$ (b), relative deviation of total mass
(c), absolute total momentum $J_{1}$ (d), absolute total momentum $J_{2}$ (e),
and relative deviation of total energy (f). $\varepsilon=10^{-5}$. In (b),
$r_{12}$ and $r_{34}$ are close, $r_{1}$ and $r_{2}$ are close, and $r_{3}$
and $r_{4}$ are close.

Figure 5.14: Example 5.4. Non-conservative low rank method. The time evolution
of electric energy (a), hierarchical ranks of the numerical solution of mesh
size $N_{x}\times N_{v}=64\times 128$ (b), relative deviation of total mass
(c), absolute total momentum $J_{1}$ (d), absolute total momentum $J_{2}$ (e),
and relative deviation of total energy (f). $\varepsilon=10^{-5}$. In (b),
$r_{12}$ and $r_{34}$ are close. $r_{1}$, $r_{2}$, $r_{3}$, and $r_{4}$ are
close.

######  Example 5.5.

We consider the 2D2V two-stream instability with initial condition

$f({\bf x},{\bf
v},t=0)=\frac{1}{2^{d}(2\pi)^{d/2}}\left(1+\alpha\sum_{m=1}^{d}\cos\left(kx_{m}\right)\right)\prod_{m=1}^{d}\left(\exp\left(-\frac{(v_{m}-v_{0})^{2}}{2}\right)+\exp\left(-\frac{(v_{m}+v_{0})^{2}}{2}\right)\right),$
(5.4)

where $d=2$, $\alpha=0.001$, $v_{0}=2.4$, and $k=0.2$.

The computation domain is set as $[0,L_{x}]^{2}\times[-L_{v},L_{v}]^{2}$,
where $L_{x}=\frac{2\pi}{k}$ and $L_{v}=8$, and the truncation threshold is
set as $\varepsilon=10^{-5}$. In Figures 5.15-5.16, respectively for
conservative and non-conservative methods, we report the time evolution of the
electric energy, hierarchical ranks of the numerical solution, relative
deviation of total mass and energy together with absolute total momentum
$J_{1}$ and $J_{2}$. The observation is similar to the previous example that
the proposed conservative method is able to conserve the total mass and
momentum, and meanwhile, the hierarchical ranks of the solution tensor from
the conservative method are larger than that from the non-conservative
counterpart.

Figure 5.15: Example 5.5. Conservative low rank method. The time evolution of
the electric energy (a), hierarchical ranks of the numerical solutions (b),
relative deviation of total mass (c), absolute total momentum $J_{1}$ (d),
absolute total momentum $J_{2}$ (e), and relative deviation of total energy
(f). $\varepsilon=10^{-5}$.

Figure 5.16: Example 5.5. Non-conservative low rank method. The time evolution
of electric energy (a), hierarchical ranks of the numerical solutions (b),
relative deviation of total mass (c), absolute total momentum $J_{1}$ (d),
absolute total momentum $J_{2}$ (e), and relative deviation of total energy
(f). $\varepsilon=10^{-5}$.

##  6 Conclusion

In this paper, we proposed a conservative truncation procedure for a low-rank
tensor approach for performing a grid-based Vlasov simulations. The basic idea
is initialized in the 1D1V setting, and is further developed to the 2D2V
setting with the HT tensor decompositions. The newly developed conservative
low rank tensor algorithm is theoretically proved to be a locally conservative
scheme to the macroscopic equations for charge and current densities, and is
numerically verified to globally conserve the total charge and current.
Further development of the low rank tensor algorithm with local energy
conservation is subject to our future work.

## References

  * [1] F. Allmann-Rahn, R. Grauer, and K. Kormann.  A parallel low-rank solver for the six-dimensional vlasov-maxwell equations.  arXiv preprint arXiv:2201.03471, 2022.
  * [2] C. K. Birdsall and A. B. Langdon.  Plasma physics via computer simulation.  CRC press, 2004.
  * [3] J. D. Carroll and J.-J. Chang.  Analysis of individual differences in multidimensional scaling via an N-way generalization of “Eckart-Young” decomposition.  Psychometrika, 35(3):283–319, 1970.
  * [4] J. Dawson.  Particle simulation of plasmas.  Rev. Mod. Phys., 55(2):403, 1983.
  * [5] L. De Lathauwer, B. De Moor, and J. Vandewalle.  A multilinear singular value decomposition.  SIAM J. Matrix Anal. Appl., 21(4):1253–1278, 2000.
  * [6] A. Dektor and D. Venturi.  Dynamically orthogonal tensor methods for high-dimensional nonlinear pdes.  Journal of Computational Physics, 404:109125, 2020.
  * [7] V. Ehrlacher and D. Lombardi.  A dynamical adaptive tensor method for the Vlasov–Poisson system.  J. Comput. Phys., 339:285–306, 2017.
  * [8] L. Einkemmer and I. Joseph.  A mass, momentum, and energy conservative dynamical low-rank scheme for the vlasov equation.  Journal of Computational Physics, page 110495, 2021.
  * [9] L. Einkemmer and C. Lubich.  A low-rank projector-splitting integrator for the Vlasov–Poisson equation.  SIAM J. Sci. Comput., 40(5):B1330–B1360, 2018.
  * [10] L. Einkemmer and C. Lubich.  A quasi-conservative dynamical low-rank algorithm for the vlasov equation.  SIAM Journal on Scientific Computing, 41(5):B1061–B1081, 2019.
  * [11] L. Einkemmer, A. Ostermann, and C. Piazzola.  A low-rank projector-splitting integrator for the vlasov–maxwell equations with divergence correction.  J. Comput. Phys., 403:109063, 2020.
  * [12] F. Filbet and E. Sonnendrucker.  Comparison of eulerian Vlasov solvers.  Computer Physics Communications, 150(3):247–266, 2003.
  * [13] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu.  Strong stability preserving Runge-Kutta and multistep time discretizations.  World Scientific, 2011.
  * [14] L. Grasedyck.  Hierarchical singular value decomposition of tensors.  SIAM J. Matrix Anal. Appl., 31(4):2029–2054, 2010.
  * [15] M. Griebel.  A parallelizable and vectorizable multi-level algorithm on sparse grids.  In W. Hackbusch, editor, Parallel algorithms for partial differential equations, volume 31 of Notes on numerical fluid mechanics, pages 94–100. 1991.
  * [16] W. Guo and Y. Cheng.  A sparse grid discontinuous Galerkin method for high-dimensional transport equations and its application to kinetic simulations.  SIAM J. Sci. Comput., 38(6):A3381–A3409, 2016.
  * [17] W. Guo and J.-M. Qiu.  A low rank tensor representation of linear transport and nonlinear vlasov solutions and their associated flow maps.  arXiv preprint arXiv:2106.08834, 2021.
  * [18] W. Hackbusch.  Tensor spaces and numerical tensor calculus, volume 42.  Springer Science & Business Media, 2012.
  * [19] W. Hackbusch and S. Kühn.  A new scheme for the tensor representation.  J. Fourier Anal. Appl., 15(5):706–722, 2009.
  * [20] R. A. Harshman et al.  Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multimodal factor analysis.  UCLA Working Papers in Phonetics, pages 1–84, 1970.
  * [21] F. L. Hitchcock.  The expression of a tensor or a polyadic as a sum of products.  J. Math. Phys., 6(1-4):164–189, 1927.
  * [22] T. G. Kolda and B. W. Bader.  Tensor decompositions and applications.  SIAM Rev., 51(3):455–500, 2009.
  * [23] K. Kormann.  A semi-Lagrangian Vlasov solver in tensor train format.  SIAM J. Sci. Comput., 37(4):B613–B632, 2015.
  * [24] K. Kormann and E. Sonnendrücker.  Sparse grids for the Vlasov–Poisson equation.  In Sparse Grids and Applications-Stuttgart 2014, pages 163–190. Springer, 2016.
  * [25] I. Oseledets.  Tensor-train decomposition.  SIAM J. Sci. Comput., 33(5):2295–2317, 2011.
  * [26] I. V. Oseledets and S. V. Dolgov.  Solution of linear systems and matrix inversion in the TT-format.  SIAM J. Sci. Comput., 34(5):A2718–A2739, 2012.
  * [27] I. V. Oseledets and E. E. Tyrtyshnikov.  Breaking the curse of dimensionality, or how to use SVD in many dimensions.  SIAM J. Sci. Comput., 31(5):3744–3759, 2009.
  * [28] Z. Peng, R. G. McClarren, and M. Frank.  A low-rank method for two-dimensional time-dependent radiation transport calculations.  Journal of Computational Physics, 421:109735, 2020.
  * [29] S. Smolyak.  Quadrature and interpolation formulas for tensor products of certain classes of functions.  In Dokl. Akad. Nauk SSSR, volume 4, pages 240–243, 1963.
  * [30] Z. Tao, W. Guo, and Y. Cheng.  Sparse grid discontinuous Galerkin methods for the Vlasov-Maxwell system.  J. Comput. Phys: X, 3:100022, 2019.
  * [31] L. R. Tucker.  Some mathematical notes on three-mode factor analysis.  Psychometrika, 31(3):279–311, 1966.
  * [32] C. Zenger.  Sparse grids.  In Parallel Algorithms for Partial Differential Equations, Proceedings of the Sixth GAMM-Seminar, volume 31, 1990.


# Yuanfudao at SemEval-2018 Task 11: Three-way Attention and Relational
Knowledge for Commonsense Machine Comprehension

Liang Wang Meng Sun Wei Zhao Kewei Shen Jingming Liu
Yuanfudao Research
Beijing, China
<EMAIL_ADDRESS>

###### Abstract

This paper describes our system for _SemEval-2018 Task 11: Machine
Comprehension using Commonsense Knowledge_ Ostermann et al. (2018b). We use
Three-way Attentive Networks (_TriAN_) to model interactions between the
passage, question and answers. To incorporate commonsense knowledge, we
augment the input with relation embedding from the graph of general knowledge
_ConceptNet_ Speer et al. (2017). As a result, our system achieves state-of-
the-art performance with $83.95\%$ accuracy on the official test data. Code is
publicly available at https://github.com/intfloat/commonsense-rc.

##  1 Introduction

It is well known that humans have a vast amount of commonsense knowledge
acquired from everyday life. For machine reading comprehension, natural
language inference and many other NLP tasks, commonsense reasoning is one of
the major obstacles to make machines as intelligent as humans.

A large portion of previous work focus on commonsense knowledge acquisition
with unsupervised learning Chambers and Jurafsky (2008); Tandon et al. (2017)
or crowdsourcing approach Singh et al. (2002); Wanzare et al. (2016).
_ConceptNet_ Speer et al. (2017), _WebChild_ Tandon et al. (2017) and
_DeScript_ Wanzare et al. (2016) etc are all publicly available knowledge
resources. However, resources based on unsupervised learning tend to be noisy,
while crowdsourcing approach has scalability issues. There is also some
research on incorporating knowledge into NLP tasks such as reading
comprehension Lin et al. (2017); Yang and Mitchell (2017) neural machine
translation Zhang et al. (2017a) and text classification Zhang et al. (2017b)
etc. Though experiments show performance gains over baselines, these gains are
often quite marginal over the state-of-the-art system without external
knowledge.

In this paper, we present Three-way Attentive Networks(_TriAN_) for multiple-
choice commonsense reading comprehension. The given task requires modeling
interactions between the passage, question and answers. Different questions
need to focus on different parts of the passage, attention mechanism is a
natural choice and turns out to be effective for reading comprehension. Due to
the relatively small size of training data, _TriAN_ use word-level attention
and consists of only one layer of LSTM Hochreiter and Schmidhuber (1997).
Deeper models result in serious overfitting and poor generalization
empirically.

To explicitly model commonsense knowledge, relation embeddings based on
ConceptNet Speer et al. (2017) are used as additional features. ConceptNet is
a large-scale graph of general knowledge from both crowdsourced resources and
expert-created resources. It consists of over $21$ million edges and $8$
million nodes. ConceptNet shows state-of-the-art performance on tasks like
word analogy and word relatedness.

Besides, we also find that pretraining our network on other datasets helps to
improve the overall performance. There are some existing multiple-choice
English reading comprehension datasets contributed by NLP community such as
_MCTest_ Richardson et al. (2013) and _RACE_ Lai et al. (2017). Although those
datasets don’t focus specifically on commonsense comprehension, they provide a
convenient way for data augmentation. Augmented data can be used to learn
shared regularities of reading comprehension tasks.

Combining all of the aforementioned techniques, our system achieves
competitive performance on the official test set.

Figure 1: _TriAN_ Model Architecture.

##  2 Model

The overall architecture of _TriAN_ is shown in Figure 1. It consists of an
input layer, an attention layer and an output layer.

Input Layer.  A training example consists of a passage
$\\{P_{i}\\}_{i=1}^{|P|}$, a question $\\{Q_{i}\\}_{i=1}^{|Q|}$, an answer
$\\{A_{i}\\}_{i=1}^{|A|}$ and a label $y^{*}\in\\{0,1\\}$ as input. $P$, $Q$
and $A$ are all sequences of word indices. For a word $P_{i}$ in the given
passage, the input representation of $P_{i}$ is the concatenation of several
vectors:

  * •

GloVe embeddings. Pretrained $300$-dimensional _GloVe_ vector
$\mathbf{E}_{P_{i}}^{glove}$.

  * •

Part-of-speech and named-entity embeddings. Randomly initialized
$12$-dimensional part-of-speech embedding $\mathbf{E}_{P_{i}}^{pos}$ and
$8$-dimensional named-entity embedding $\mathbf{E}_{P_{i}}^{ner}$.

  * •

Relation embeddings. Randomly initialized $10$-dimensional relation embedding
$\mathbf{E}_{P_{i}}^{rel}$. The relation is determined by querying ConceptNet
whether there is an edge between $P_{i}$ and any word in question
$\\{Q_{i}\\}_{i=1}^{|Q|}$ or answer $\\{A_{i}\\}_{i=1}^{|A|}$. If there exist
multiple different relations, just randomly choose one.

  * •

Handcrafted features. We also add logarithmic term frequency feature and co-
occurrence feature $\mathbf{f}_{P_{i}}$. Term frequency is calculated based on
English Wikipedia. Co-occurrence feature is a binary feature which is true if
$P_{i}$ appears in question $\\{Q_{i}\\}_{i=1}^{|Q|}$ or answer
$\\{A_{i}\\}_{i=1}^{|A|}$.

The input representation for $P_{i}$ is $\mathbf{w}_{P_{i}}$:

$\mathbf{w}_{P_{i}}=[\mathbf{E}_{P_{i}}^{glove};\mathbf{E}_{P_{i}}^{pos};\mathbf{E}_{P_{i}}^{ner};\mathbf{E}_{P_{i}}^{rel};\mathbf{f}_{P_{i}}]$
(1)

In a similar way, we can get input representation for question
$\mathbf{w}_{Q_{i}}$ and answer $\mathbf{w}_{A_{i}}$.

Attention Layer.  We use word-level attention to model interactions between
the given passage $\\{P_{i}\\}_{i=1}^{|P|}$, the question
$\\{Q_{i}\\}_{i=1}^{|Q|}$ and the answer $\\{A_{i}\\}_{i=1}^{|A|}$. First,
let’s define a sequence attention function Chen et al. (2017):

$\displaystyle Att_{seq}(\mathbf{u},\
\\{\mathbf{v}_{i}\\}_{i=1}^{n})=\sum_{i=1}^{n}$
$\displaystyle{\alpha_{i}\mathbf{v}_{i}}$ (2)
$\displaystyle\alpha_{i}=\mathrm{softmax}_{i}(f(\mathbf{W}_{1}\mathbf{u})^{T}$
$\displaystyle f(\mathbf{W}_{1}\mathbf{v}_{i}))$

$\mathbf{u}$ and $\mathbf{v}_{i}$ are vectors and $\mathbf{W}_{1}$ is a
matrix. $f$ is a non-linear activation function and is set to $ReLU$.

Question-aware passage representation
$\\{\mathbf{w}_{P_{i}}^{q}\\}_{i=1}^{|P|}$ can be calculated as:
$\mathbf{w}_{P_{i}}^{q}=Att_{seq}(\mathbf{E}_{P_{i}}^{glove},\\{\mathbf{E}_{Q_{i}}^{glove}\\}_{i=1}^{|Q|})$.
Similarly, we can get passage-aware answer representation
$\\{\mathbf{w}_{A_{i}}^{p}\\}_{i=1}^{|A|}$ and question-aware answer
representation $\\{\mathbf{w}_{A_{i}}^{q}\\}_{i=1}^{|A|}$. Then three BiLSTMs
are applied to the concatenation of those vectors to model the temporal
dependency:

$\displaystyle\mathbf{h}^{q}$
$\displaystyle=\mathrm{BiLSTM}(\\{\mathbf{w}_{Q_{i}}\\}_{i=1}^{|Q|})$ (3)
$\displaystyle\mathbf{h}^{p}$
$\displaystyle=\mathrm{BiLSTM}(\\{[\mathbf{w}_{P_{i}};\mathbf{w}_{P_{i}}^{q}]\\}_{i=1}^{|P|})$
$\displaystyle\mathbf{h}^{a}$
$\displaystyle=\mathrm{BiLSTM}(\\{[\mathbf{w}_{A_{i}};\mathbf{w}_{A_{i}}^{p};\mathbf{w}_{A_{i}}^{q}]\\}_{i=1}^{|A|})$

$\mathbf{h}^{p},\mathbf{h}^{q},\mathbf{h}^{a}$ are the new representation
vectors that incorporates more context information.

Output Layer.  Question sequence and answer sequence representation
$\mathbf{h}^{q},\mathbf{h}^{a}$ are summarized into fixed-length vectors with
self-attention Yang et al. (2016). Self-attention function is defined as
follows:

$\displaystyle Att_{self}(\\{\mathbf{u}_{i}\\}_{i=1}^{n})=\sum_{i=1}^{n}$
$\displaystyle{\alpha_{i}\mathbf{u}_{i}}$ (4)
$\displaystyle\alpha_{i}=\mathrm{softmax}_{i}(\mathbf{W}_{2}^{T}\mathbf{u}_{i})$

Then we have question representation
$\mathbf{q}=Att_{self}(\\{\mathbf{h}_{i}^{q}\\}_{i=1}^{|Q|})$, answer
representation $\mathbf{a}=Att_{self}(\\{\mathbf{h}_{i}^{a}\\}_{i=1}^{|A|})$
and passage representation
$\mathbf{p}=Att_{seq}(\mathbf{q},\\{\mathbf{h}_{i}^{p}\\}_{i=1}^{|P|})$. The
final output $y$ is based on their bilinear interactions:

$y=\sigma(\mathbf{p}^{T}\mathbf{W}_{3}\mathbf{a}+\mathbf{q}^{T}\mathbf{W}_{4}\mathbf{a})$
(5)

Model Learning.  We first pretrain _TriAN_ on _RACE_ dataset for $10$ epochs.
Then our model is fine-tuned on official training data. Standard cross entropy
function is used as the loss function to minimize.

##  3 Experiments

###  3.1 Setup

Data.  For data preprocessing, we use _spaCy_
111https://github.com/explosion/spaCy for tokenization, part-of-speech tagging
and named-entity recognition. Statistics for official dataset _MCScript_
Ostermann et al. (2018a) are shown in Table 1. _RACE_
222http://www.cs.cmu.edu/~glai1/data/race/ dataset is used for network
pretraining. English stop words are ignored when computing handcrafted
features. Input word embeddings are initialized with $300$-dimensional _GloVe_
Pennington et al. (2014) vectors
333http://nlp.stanford.edu/data/glove.840B.300d.zip.

| train | dev | test
---|---|---|---
# of examples | 9731 | 1411 | 2797
Table 1: Official dataset statistics.

Hyperparameters.  Our model _TriAN_ is implemented based on _PyTorch_
444http://pytorch.org/. Models are trained on a single GPU(Tesla P40) and each
epoch takes about $80$ seconds. Only the word embeddings of top $10$ frequent
words are fine-tuned during training. The dimension of both forward and
backward LSTM hidden state is set to $96$. Dropout rate is set to $0.4$ for
both input embeddings and BiLSTM outputs Srivastava et al. (2014). For
parameter optimization, we use _Adamax_ Kingma and Ba (2014) with an initial
learning rate $2\times 10^{-3}$. Learning rate is then halved after $10$ and
$15$ training epochs. The model converges after $50$ epochs. Gradients are
clipped to have a maximum L2 norm of $10$. Minibatch with batch size $32$ is
used. Hyperparameters are optimized by random search strategy Bergstra and
Bengio (2012). Our model is quite robust over a wide range of hyperparameter
configurations.

###  3.2 Main Results

The experimental results are shown in Table 2. Human performance is shared by
task organizers. For _TriAN-ensemble_ , we average the output probabilities of
9 models trained with the same datasets and network architecture but different
random seeds. _TriAN-ensemble_ is the model that we used for official
submission.

model | dev | test
---|---|---
Random | 50.00% | 50.00%
TriAN-RACE | 64.78% | 64.28%
TriAN-single | 83.84% | 81.94%
TriAN-ensemble | 85.27% | 83.95%
HFL | – | 84.13%
Human | – | 98.00%
Table 2: Main results. _TriAN-RACE_ only use _RACE_ dataset for training;
_HFL_ is the 1st place team for _SemEval-2018 Task 11_. The evaluation metric
is accuracy.

From Table 2, we can see that even though _RACE_ dataset contains nearly
$100k$ questions, _TriAN-RACE_ achieves quite poor results. The accuracy on
development set is only $64.78\%$, which is worse than most participants’
systems. However, pretraining acts as a way of implicit knowledge transfer and
is beneficial for overall performance, as will be seen in Section 3.3. The
accuracy of our system _TriAN-ensemble_ is very close to the 1st place team
_HFL_ with $0.18\%$ difference. Yet there is still a large gap between machine
learning models and human.

We also compared the performances of shallow and deep TriAN models. On
datasets such as _SQuAD_ Rajpurkar et al. (2016), deep models typically works
better than shallow ones. Notice that the attention layer in our proposed
_TriAN_ model can be stacked multiple times if we treat the output vectors of
BiLSTMs as new input representations.

model | dev | test
---|---|---
1-layer TriAN-single | 83.84% | 81.94%
2-layer TriAN-single | 82.71% | 80.55%
Table 3: Accuracy comparison of shallow and deep TriAN models.

Maybe a little bit surprising, Table 3 shows that _2-layer TriAN_ model
performs worse than _1-layer TriAN_. One possible explanation is that the
labeled dataset is relatively small and deeper models tend to easily overfit.

###  3.3 Ablation Study

The input representation consists of several components: part-of-speech
embedding, relation embedding and handcrafted features etc. We conduct an
ablation study to investigate the effects of each component. The results are
in Table 4.

model | dev | test
---|---|---
TriAN-single | 83.84% | 81.94%
w/o pretraining | 82.71% | 80.51%
w/o ConceptNet | 82.78% | 81.08%
w/o POS | 82.84% | 81.27%
w/o features | 82.92% | 81.35%
w/o NER | 83.60% | 81.87%
Table 4: Ablation study for input representation.

Pretraining on _RACE_ dataset turns out to be the most important factor.
Without pretraining, the accuracy drops by more than $1\%$ on both development
and test set. Relation embeddings based on ConceptNet make approximately $1\%$
difference. Part-of-speech and named-entity embeddings are also helpful. In
fact, combining input representations from multiple sources has been a
standard practice for reading comprehension tasks.

At attention layer, our proposed _TriAN_ involves applying several attention
functions to model interactions between different text sequences. It would be
interesting to examine the importance of each attention function, as shown in
Table 5.

model | dev | test
---|---|---
TriAN-single | 83.84% | 81.94%
w/o passage-question attention | 83.51% | 82.20%
w/o passage-answer attention | 83.07% | 81.39%
w/o question-answer attention | 83.23% | 81.84%
w/o attention | 81.93% | 80.65%

Table 5: Ablation study for attention. The last one “ _w/o attention_ ”
removes all word-level attentions.

Interestingly, removing any of the three word-level sequence attentions does
not seem to hurt the performance much. In fact, removing passage-question
attention even results in higher accuracy on test set than _TriAN-single_.
However, if we remove all word-level attentions, the performance drastically
drops by $1.9\%$ on development set and $1.3\%$ on test set.

###  3.4 Discussion

Even though our system is built for commonsense reading comprehension, it
doesn’t have any explicit knowledge reasoning component. Relation embeddings
based on ConceptNet only serve as additional input features. Methods like
event calculus Mueller (2014) are more rigorous mathematically and resemble
the way of how human brain works. The problem of event calculus is that it
requires large amounts of domain-specific axioms and therefore doesn’t scale
well.

Another limitation is that our system relies on hard-coded commonsense
knowledge bases, just like most systems for commonsense reasoning. For humans,
commonsense knowledge comes from constant interactions with the real-world
environment. From our point of view, it is quite hopeless to enumerate all of
them.

There are a lot of reading comprehension datasets available. When the size of
training data is relatively small like this _SemEval-2018 task_ , transfer
learning among different datasets is a useful technique. This paper shows that
pretraining is a simple and effective method. However, it still remains to be
seen whether there is a better alternative approach.

##  4 Conclusion

In this paper, we present the core ideas and design philosophy for our system
_TriAN_ at _SemEval-2018 Task 11: Machine Comprehension using Commonsense
Knowledge_. We build upon recent progress on neural models for reading
comprehension and incorporate commonsense knowledge from ConceptNet.
Pretraining and handcrafted features are also proved to be helpful. As a
result, our proposed model _TriAN_ achieves near state-of-the-art performance.

## Acknowledgements

We would like to thank SemEval 2018 task organizers and several anonymous
reviewers for their helpful comments.

## References

  * Bergstra and Bengio (2012) James Bergstra and Yoshua Bengio. 2012.  Random search for hyper-parameter optimization.  _Journal of Machine Learning Research_ , 13(Feb):281–305.
  * Chambers and Jurafsky (2008) Nathanael Chambers and Dan Jurafsky. 2008.  Unsupervised learning of narrative event chains.  _Proceedings of ACL-08: HLT_ , pages 789–797.
  * Chen et al. (2017) Danqi Chen, Adam Fisch, Jason Weston, and Antoine Bordes. 2017.  Reading wikipedia to answer open-domain questions.  _arXiv preprint arXiv:1704.00051_.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780.
  * Kingma and Ba (2014) Diederik Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_.
  * Lai et al. (2017) Guokun Lai, Qizhe Xie, Hanxiao Liu, Yiming Yang, and Eduard Hovy. 2017.  Race: Large-scale reading comprehension dataset from examinations.  _arXiv preprint arXiv:1704.04683_.
  * Lin et al. (2017) Hongyu Lin, Le Sun, and Xianpei Han. 2017.  Reasoning with heterogeneous knowledge for commonsense machine comprehension.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , pages 2032–2043.
  * Mueller (2014) Erik T Mueller. 2014.  _Commonsense reasoning: an event calculus based approach_.  Morgan Kaufmann.
  * Ostermann et al. (2018a) Simon Ostermann, Ashutosh Modi, Michael Roth, Stefan Thater, and Manfred Pinkal. 2018a.  MCScript: A Novel Dataset for Assessing Machine Comprehension Using Script Knowledge.  In _Proceedings of the 11th International Conference on Language Resources and Evaluation (LREC 2018)_ , Miyazaki, Japan.
  * Ostermann et al. (2018b) Simon Ostermann, Michael Roth, Ashutosh Modi, Stefan Thater, and Manfred Pinkal. 2018b.  SemEval-2018 Task 11: Machine Comprehension using Commonsense Knowledge.  In _Proceedings of International Workshop on Semantic Evaluation(SemEval-2018)_ , New Orleans, LA, USA.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning. 2014.  Glove: Global vectors for word representation.  In _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_ , pages 1532–1543.
  * Rajpurkar et al. (2016) Pranav Rajpurkar, Jian Zhang, Konstantin Lopyrev, and Percy Liang. 2016.  Squad: 100,000+ questions for machine comprehension of text.  _arXiv preprint arXiv:1606.05250_.
  * Richardson et al. (2013) Matthew Richardson, Christopher JC Burges, and Erin Renshaw. 2013.  Mctest: A challenge dataset for the open-domain machine comprehension of text.  In _Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing_ , pages 193–203.
  * Singh et al. (2002) Push Singh, Thomas Lin, Erik T Mueller, Grace Lim, Travell Perkins, and Wan Li Zhu. 2002.  Open mind common sense: Knowledge acquisition from the general public.  In _OTM Confederated International Conferences” On the Move to Meaningful Internet Systems”_ , pages 1223–1237. Springer.
  * Speer et al. (2017) Robert Speer, Joshua Chin, and Catherine Havasi. 2017.  Conceptnet 5.5: An open multilingual graph of general knowledge.  In _AAAI_ , pages 4444–4451.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: A simple way to prevent neural networks from overfitting.  _The Journal of Machine Learning Research_ , 15(1):1929–1958.
  * Tandon et al. (2017) Niket Tandon, Gerard de Melo, and Gerhard Weikum. 2017.  Webchild 2.0: fine-grained commonsense knowledge distillation.  _Proceedings of ACL 2017, System Demonstrations_ , pages 115–120.
  * Wanzare et al. (2016) Lilian DA Wanzare, Alessandra Zarcone, Stefan Thater, and Manfred Pinkal. 2016.  Descript: A crowdsourced corpus for the acquisition of high-quality script knowledge.  In _The International Conference on Language Resources and Evaluation._
  * Yang and Mitchell (2017) Bishan Yang and Tom Mitchell. 2017.  Leveraging knowledge bases in lstms for improving machine reading.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , volume 1, pages 1436–1446.
  * Yang et al. (2016) Zichao Yang, Diyi Yang, Chris Dyer, Xiaodong He, Alex Smola, and Eduard Hovy. 2016\.  Hierarchical attention networks for document classification.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 1480–1489.
  * Zhang et al. (2017a) Jiacheng Zhang, Yang Liu, Huanbo Luan, Jingfang Xu, and Maosong Sun. 2017a.  Prior knowledge integration for neural machine translation using posterior regularization.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , volume 1, pages 1514–1523.
  * Zhang et al. (2017b) Ye Zhang, Matthew Lease, and Byron C Wallace. 2017b.  Exploiting domain knowledge via grouped weight sharing with application to text categorization.  _arXiv preprint arXiv:1702.02535_.