# Robust Implicit Adaptive Low Rank Time-Stepping Methods for Matrix
Differential Equations

Daniel Appelö Department of Mathematics, Virginia Tech, Blacksburg, VA 24061
U.S.A<EMAIL_ADDRESS>Research supported by DOE Office of Advanced Scientific
Computing Research under the Advanced Research in Quantum Computing program,
subcontracted from award 2019-LLNL-SCW-1683, NSF DMS-2208164, and Virginia
Tech.  Yingda Cheng Department of Mathematics, Virginia Tech, Blacksburg, VA
24061 U.S.A<EMAIL_ADDRESS>Research supported by NSF DMS-2011838, DOE grant
DE-SC0023164 and Virginia Tech.

###### Abstract

In this work, we develop implicit rank-adaptive schemes for time-dependent
matrix differential equations. The dynamic low rank approximation (DLRA) is a
well-known technique to capture the dynamic low rank structure based on
Dirac–Frenkel time-dependent variational principle. In recent years, it has
attracted a lot of attention due to its wide applicability. Our schemes are
inspired by the three-step procedure used in the rank adaptive version of the
unconventional robust integrator (the so called BUG integrator) [2] for DLRA.
First, a prediction (basis update) step is made computing the approximate
column and row spaces at the next time level. Second, a Galerkin evolution
step is invoked using a base implicit solve for the small core matrix.
Finally, a truncation is made according to a prescribed error threshold. Since
the DLRA is evolving the differential equation projected on to the tangent
space of the low rank manifold, the error estimate of the BUG integrator
contains the tangent projection (modeling) error which cannot be easily
controlled by mesh refinement. This can cause convergence issue for equations
with cross terms.

To address this issue, we propose a simple modification, consisting of merging
the row and column spaces from the explicit step truncation method together
with the BUG spaces in the prediction step. In addition, we propose an
adaptive strategy where the BUG spaces are only computed if the residual for
the solution obtained from the prediction space by explicit step truncation
method, is too large. We prove stability and estimate the local truncation
error of the schemes under assumptions. We benchmark the schemes in several
tests, such as anisotropic diffusion, solid body rotation and the combination
of the two, to show robust convergence properties.

##  1 Introduction

In this work, we are interested in solving linear matrix differential
equations in the following form

$\frac{d}{dt}X(t)=F(X(t),t),\quad X(t)\in\mathbb{R}^{m_{1}\times m_{2}},\quad
X(0)=X_{0},$ (1)

where

$F(X(t),t)=\sum_{j=1}^{s}A_{j}X(t)B_{j}^{T}+G(t),$

and $A_{j}\in\mathbb{R}^{m_{1}\times m_{1}},B_{j}\in\mathbb{R}^{m_{2}\times
m_{2}}$ are sparse or structured matrices for which a fast matrix-vector
product is assumed to be known. Further, $G(t)\in\mathbb{R}^{m_{1}\times
m_{2}}$ is a given function that is assumed to have a known low rank
decomposition. The integer $s,$ which denotes the separation rank of
$F(\cdot),$ is assumed not to be too large. We focus on implicit methods,
which are particularly needed for computing solutions to stiff equations.

Equation (1) arises in many applications governed by partial differential
equations (PDE). For example, using the method of lines approach, we can cast
numerical discretizations of two-dimensional (or even higher dimensional)
linear time-dependent convection-diffusion-reaction equations in the form of
(1). If we represent the unknowns on a Cartesian grid, the elements of $X(t)$
will represent the unknowns (e.g. point values of the solution) at a given
time $t.$ The matrices $A_{j},B_{j}$ will, e.g., correspond to difference
operators or variable coefficients, and $G(t)$ encodes the given source terms
and boundary conditions. Compared to the traditional approach of numerical PDE
solvers, the matrix approach does not vectorize the unknowns into a single
vector at each time step, but rather we represent the unknowns in a matrix
format, retaining the relative relations in each spatial dimension.

In recent years, the matrix (or tensor in higher dimensions) approach received
increased attention in numerical analysis because it offers a way to tame _the
curse of dimensionality_ , enabling approximations to solutions of high
dimensional PDEs [10, 13]. The overarching idea is to expose low rank
structure in the solution manifold in order to drastically reduce storage and
computation. This is the idea of low rank tensor methods for PDEs [7, 1]. In
two dimensions, this approach reduces to the low rank matrix method, where it
is assumed that the matrix $X(t)$ has a small rank $r.$ As is well known, a
rank $r$ matrix of size $m_{1}\times m_{2}$ allows for a SVD decomposition
with only $(m_{1}+m_{2}-r)r$ degrees of freedom. This is much smaller than
$m_{1}m_{2}$, the degrees of freedom of the full matrix, when $r$ is small
compared to $m_{1}$ and $m_{2}$.

If the rank of $X(t)$ is fixed, i.e. we constrain the solution $X(t)$ to live
on the rank $r$ matrix manifold
$\mathcal{M}_{r}=\\{X\in\mathbb{R}^{m_{1}\times m_{2}},{\sf rank}(X)=r\\},$
then an effective numerical solution can be computed by the dynamic low rank
approximation (DLRA) [15, 16]. The DLRA solves
$\frac{d}{dt}X(t)=\Pi_{X(t)}F(X(t),t),$ where $\Pi_{X(t)}$ is the orthogonal
projection onto the tangent space $T_{X(t)}\mathcal{M}_{r}$ of the manifold
$\mathcal{M}_{r}$ at $X(t).$ The DLRA has found success in many applications,
including extensions to various tensor formats [18, 17]. However, the
constraint that $X(t)\in\mathcal{M}_{r}$ requires that the rank $r$ is
estimated a-priori and this is in contrast with the rank adaptive schemes.

A rank adaptive method adaptively choose the rank at each time step according
to a prescribed error tolerance. Several approaches to rank adaptivity can be
found in the literature. For example, a straightforward idea is the so-called
step truncation method [4] which evolves the low rank solution for one time
step by a traditional time stepping method in an ambient space of higher rank,
then performs a truncation (by SVD with given tolerance). This method is
intuitive and effective for explicit time stepping methods. For implicit
schemes [21], implementation of step truncation methods is nontrivial and
requires an effective iterative solver in adaptive low rank format [7, 1].
Another approach that has been suggested by multiple groups, and that is
similar to what we propose, is to modify the DLRA method by using rank
increase and decrease indicators and augment or truncate the spaces at each
time step [12, 11, 24]. In particular the ideas in [24], while being
specialized to time dependent Schrödinger, are similar to the approach we take
here. A promising idea in [2] uses the so called unconventional (BUG)
integrator [3] to achieve rank adaptivity. This approach was recently extended
to higher order in time in [19].

In the BUG integrator [3], two subproblems that come from the projector
splitting of DLRA [16] are solved. Then, the solutions to the two subproblems
are used to update the column and row spaces of the matrix; and a Galerkin
evolution step is performed in the resulting space, searching for the
coefficients of a small matrix. For the rank adaptive version of the
unconventional integrator [2], this translates to solving two matrix
differential equation of sizes $m_{1}\times r,m_{2}\times r$ in the first step
and one matrix differential equation of size $2r\times 2r$ in the Galerkin
step, followed by a truncated SVD step. The ideas in [3, 2] can be implemented
for implicit methods because the matrix size is fixed in each subproblems, so
a standard linear solver will suffice. However, the errors of both methods in
[3, 2] are subject to the tangent projection error or the so called _modeling
error_ [14].

It is important to understand that this modeling error is present in many
problems and can, as we will see below, result in non convergent numerical
methods. We now describe a trivial but prototypical example that will make the
DLRA/BUG method fail. Assume that $X(t)=U\Sigma V^{T}$ and that the column and
row space of the matrix $F(X(t),t)$ are both (respectively) orthogonal to $U$
and $V$, then $\Pi_{X(t)}F(X(t),t)=0,$ which means the solution to the DLRA
description will remain stationary. For example, this will happen for the
numerical discretization when $X(t)$ is a rank 1 even function in both
variables (say $X(t,x_{1},x_{2})=\exp(-x_{1}^{2})\exp(-x_{2}^{2})$) on a
square domain with center $(0,0)$, and $F(X(t),t)=AXB^{T}$, with $A={\sf
diag}(x_{1})$ and $B={\sf diag}(x_{2})$.

The main contribution of this paper is to propose a simple improvement to the
BUG integrator that enhances the robustness with respect to convergence. In
our method, we merge the row and column spaces from the explicit step
truncation method with the BUG spaces in the prediction step to alleviate the
tangent projection error associated with DLRA. As the BUG spaces can be
expensive to compute when implicit schemes are used, we propose an adaptive
acceleration strategy based on the residual of the underlying classic scheme.
The strategy is straightforward, if the residual obtained when only the step
truncation spaces are used, is larger than the truncation threshold then we
add in the BUG spaces. We observe that this strategy is robust in convergence
and can reduce the computational time for moderately stiff problems.

The rest of the paper is organized as follows. Section 2 reviews the
background on low rank time integrators. In Section 3, we describe the
proposed schemes and perform numerical analysis of their properties. Section 4
provides numerical experiments and in Section 5 we draw conclusions.

##  2 Background review

In this section, we will review basic concepts of low rank time stepping
methods for (1). In Section 2.1, we gather the notations used in the paper.
Section 2.2 reviews an explicit scheme truncation method based on the forward
Euler method, while Section 2.3 reviews the rank adaptive BUG scheme.

###  2.1 Notations and preliminaries

In this paper, we use $\|\cdot\|$ to denote the matrix Frobenius norm. The
Frobenius norm is a natural choice for matrix functions as $\|A\|$ coincides
with the $L^{2}$ vector norm of the vectorized matrix ${\sf vec}(A),$ which is
commonly used in the analysis of the standard ODE/PDE solvers. We also use
matrix inner product defined as follows, for $A,B\in\mathbb{R}^{m_{1}\times
m_{2}},$ we define $\langle
A,B\rangle=\sum_{i,j=1}^{m_{1},m_{2}}A_{ij}B_{i,j}.$ Then $\langle
A,A\rangle=\|A\|^{2}.$ The following property,

$\langle A,USV^{T}\rangle=\langle U^{T}AV,S\rangle,$ (2)

holds for any $U\in\mathbb{R}^{m_{1}\times r},S\in\mathbb{R}^{r\times
s},V\in\mathbb{R}^{m_{2}\times s}$. Here we stop to caution the reader that
throughout this paper a matrix $\Sigma,S,$ or $\mathcal{S}$ does not always
represent the diagonal matrix holding the singular values of another matrix.

The numerical solution at $t^{n}$ is denoted by
$\hat{X}^{n}\in\mathbb{R}^{m_{1}\times m_{2}}.$ In particular, this is a rank
$r_{n}$ matrix with singular value decomposition (SVD)
$\hat{X}^{n}=U^{n}\Sigma^{n}(V^{n})^{T},$ where
$U^{n}\in\mathbb{R}^{m_{1}\times r_{n}}$, $V^{n}\in\mathbb{R}^{m_{2}\times
r_{n}}$ have orthogonal columns, and $\Sigma^{n}={\sf
diag}(\sigma_{1},\cdots,\sigma_{r_{n}})$ is a $r_{n}\times r_{n}$ diagonal
matrix with diagonal entries $\sigma_{1}\geq\cdots\geq\sigma_{r_{n}}>0.$ Here,
the rank $r_{n}$ will be chosen adaptively by the numerical scheme so that
$\hat{X}^{n}$ approximate the true solution $X(t^{n})$ with a prescribed
accuracy.

We use $\mathcal{T}_{\epsilon}$ to denote a generic matrix approximation
operator with accuracy $\epsilon$ in the Frobenius norm, i.e.
$\|A-\mathcal{T}_{\epsilon}(A)\|\leq\epsilon.$ A prominent example is the
truncated SVD of a matrix. Namely, given a generic rank $r$ matrix
$A\in\mathbb{R}^{m_{1}\times m_{2}}$ with reduced SVD: $A=U\Sigma V^{T},$ then
$\mathcal{T}^{svd}_{\epsilon}(A)=U[:,1:s]\,{\sf
diag}(\sigma_{1},\cdots,\sigma_{r})\,V[:,1:s]^{T},$ where we have used the
standard MATLAB notation for submatrices and $r$ is chosen to be the smallest
integer so that $(\sum_{j=s+1}^{r}\sigma_{j}^{2})^{1/2}\leq\epsilon.$

In this paper, we will also encounter low rank matrix sum operations. This
frequently appearing operation in low rank methods can be efficiently computed
by Algorithm 1. The inputs to this algorithm are the SVD representations of
the matrices $X_{j}=U_{j}\Sigma_{j}(V_{j})^{T},j=1\ldots m,$ and the output is
the truncated SVD of their sum
$\mathcal{T}^{sum}_{\epsilon}(\sum_{j=1}^{m}X_{j})=\mathcal{U}\mathcal{S}\mathcal{V}^{T}.$
We first write $U=[U_{1},\ldots,U_{m}],$ $\Sigma={\sf
diag}(\Sigma_{1},\ldots,\Sigma_{m}),$ $V=[V_{1},\ldots,V_{m}]$. Then use the
column-pivoted QR decomposition, [6], denoted by $[Q,R,\Pi]={\sf qr}(A)$,
which computes $QR\Pi=A.$ In this factorization the columns of $Q$ are
orthogonal as usual, but the introduction of the permutation matrix $\Pi$,
makes it possible to guarantee that the diagonal entries of the upper
triangular matrix $R$ are strictly decreasing in magnitude. Assume that the
column pivoted QR procedure applied to to $U$ and $V$ yields the
factorizations $Q_{1}R_{1}\Pi_{1}$, and $Q_{2}R_{2}\Pi_{2}$. Then the
truncated matrix sum
$\mathcal{T}^{sum}_{\epsilon}(\sum_{j=1}^{m}X_{j})=Q_{1}\mathcal{U}\mathcal{S}(Q_{2}\mathcal{V})^{T}.$
is obtained by a truncated SVD on the small matrix
$\mathcal{T}^{svd}_{\epsilon}(R_{1}\Pi_{1}\,\Sigma\,\Pi_{2}^{T}R_{2}^{T})=\mathcal{U}\mathcal{S}\mathcal{V}^{T}$.

1

Input : low rank matrices $X_{j},j=1\ldots m$ or their SVD
$U_{j}\Sigma_{j}(V_{j})^{T},j=1\ldots m$

Output : truncated SVD of their sum
$\mathcal{T}^{sum}_{\epsilon}(\sum_{j=1}^{m}X_{j})=\mathcal{U}\mathcal{S}\mathcal{V}^{T}$

Parameter : tolerance $\epsilon$

2 Form $U=[U_{1},\ldots,U_{m}],\,\Sigma={\sf
diag}(\Sigma_{1},\ldots,\Sigma_{m}),\,V=[V_{1},\ldots,V_{m}]$.

3Perform column pivoted QR: $[Q_{1},R_{1},\Pi_{1}]={\sf qr}(U),$
$[Q_{2},R_{2},\Pi_{2}]={\sf qr}(V).$

4Compute the truncated SVD:
$\mathcal{T}^{svd}_{\epsilon}(R_{1}\Pi_{1}\,\Sigma\,\Pi_{2}^{T}R_{2}^{T})=\mathcal{U}\mathcal{S}\mathcal{V}^{T}$.

5Form $\mathcal{U}\leftarrow Q_{1}\mathcal{U},\mathcal{V}\leftarrow
Q_{2}\mathcal{V}$.

Return
$[\mathcal{U},\mathcal{S},\mathcal{V}]=\mathcal{T}^{sum}_{\epsilon}(\sum_{j=1}^{m}X_{j})$

Algorithm 1 Sum of low rank matrices

###  2.2 Explicit step truncation schemes

The primary assumption to guarantee the efficiency of any low rank solver is
that $r_{n}\ll m_{1},r_{n}\ll m_{2},$ i.e. the solution can be well
approximated by a low rank matrix. When such assumptions hold, the main steps
of the numerical algorithm should specify the evolution of
$U^{n},\Sigma^{n},V^{n}$ based on the matrix differential equation (1). For
explicit schemes, a simple approach is to numerically integrate $\hat{X}^{n}$
according to a standard time integrator to $t^{n+1}$ and then perform a
truncated SVD according to the error threshold to obtain $\hat{X}^{n+1}.$ This
is the so-called step truncation method. Methods of this type and their
variations have been discussed in [14, 4]. The main steps of the rank adaptive
forward Euler time scheme are highlighted in Algorithm 2. The first step of
Algorithm 2 is a standard forward Euler step for (1) with the truncated right
hand side in Line 2. This will result in a numerical solution
$\hat{X}^{n+1,pre}$ with higher rank than needed. Then, the second step in
Line 2 is to truncate the solution to a lower rank matrix. In both steps, we
use Algorithm 1 to perform the sum of the low rank matrices involved.

1

Input : numerical solution at $t^{n}:$ rank $r_{n}$ matrix $\hat{X}^{n}$ in
its SVD form $U^{n}\Sigma^{n}(V^{n})^{T}.$

Output : numerical solution at $t^{n+1}:$ rank $r_{n+1}$ matrix
$\hat{X}^{n+1}$ in its SVD form $U^{n+1}\Sigma^{n+1}(V^{n+1})^{T}.$

Parameter : time step $\Delta t$, error tolerance $\epsilon_{1},\epsilon_{2}$

2 (Evolution). $\hat{X}^{n+1,pre}=\hat{X}^{n}+\Delta
t\mathcal{T}^{sum}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})).$

(Truncation).
$\hat{X}^{n+1}=\mathcal{T}^{sum}_{\epsilon_{2}}(\hat{X}^{n+1,pre}).$

Algorithm 2 Forward Euler scheme $t^{n}\rightarrow t^{n+1}$

The convergence of this algorithm is well understood [21]. The local
truncation error is on the order of $O(\Delta t^{2}+\epsilon_{1}\Delta
t+\epsilon_{2}).$ If one choose $\epsilon_{1}=O(\Delta
t),\epsilon_{2}=O(\Delta t^{2}),$ we will obtain a first order accurate
solution. As for computational cost, we can see that the procedure involves QR
factorization of tall matrices and SVD for a small matrix. If we assume
$r_{n}=O(r),$ $s$ is small (i.e. $s=O(1)$) and $G$ has low rank (i.e. ${\sf
rank}(G)=O(r)$), the computational complexity is on the order of
$O(m_{1}r^{2}+m_{2}r^{2}+r^{3}).$ Algorithm 2 can be readily extended to
higher order by embedding in Runge-Kutta, multistep methods, and was also
considered with tangent projection in the projected Runge-Kutta schemes [14,
4, 8].

###  2.3 Dynamic low rank approximation

The DLRA modifies the equation, and solves

$\frac{d}{dt}X(t)=\Pi_{X(t)}F(X(t),t),$ (3)

where $\Pi_{X(t)}$ is the orthogonal projection onto the tangent space
$T_{X(t)}\mathcal{M}_{r}$ of the the rank $r$ matrix manifold
$\mathcal{M}_{r}=\\{X\in\mathbb{R}^{m_{1}\times m_{2}},{\sf rank}(X)=r\\}$ at
$X(t).$ The DLRA is particularly suited for a fixed rank calculation, if the
computation is constrained on $\mathcal{M}_{r}.$ It has found success in many
applications, including extensions to various tensor formats [18, 17] and
particularly in quantum mechanics [20].

We reminder the reader that if the SVD of $X\in\mathcal{M}_{r}$ is given by
$U\Sigma V^{T},$ then the tangent space of $\mathcal{M}_{r}$ at $X$ is given
by

$T_{X}\mathcal{M}_{r}=\left\\{[U\quad
U_{\perp}]\begin{bmatrix}\mathbb{R}^{r\times
r}&\mathbb{R}^{r\times(m_{2}-r)}\\\ \mathbb{R}^{(m_{1}-r)\times
r}&0^{(m_{1}-r)\times(m_{2}-r)}\end{bmatrix}[V\quad V_{\perp}]^{T}\right\\},$

where $U_{\perp},V_{\perp}$ are orthogonal complements of $U,V$ in
$\mathbb{R}^{m_{1}},\mathbb{R}^{m_{2}},$ respectively. As we can see later in
the PDE examples, typically if there are cross terms like $F(X)=AXB$ (e.g.
cross derivatives or variable coefficient problem with operators acting on
both left and right side of the matrix), then in general $AXB\notin
T_{X}\mathcal{M}_{r}.$ In this case, the _modeling error_ [14] (i.e. the error
from the tangent projection, which is $\|(I-\Pi_{X(t)})F(X,t)\|$) will be
evident in the numerical approximation, see for example error estimates in [3,
2].

Below we will review the rank adaptive BUG integrator in [2] for the
discretization of (3). First, two subproblems that come from the projector
splitting of DLRA [16] are solved in the $K-$ and $L-$step. Then, the
solutions to the two subproblems are used to update the column and row spaces
of the matrix; and a Galerkin evolution step is performed in the resulting
space, searching for the optimal solution with the Galerkin condition. The
idea of finding the appropriate subspace and computing the solution in that
space is similar to a widely used technique in numerical linear algebra called
projection methods [22].

For completeness, we describe the rank adaptive BUG integrator in Algorithm 3.
The original algorithm is written in a time continuous format, but here we
present it using an implicit Euler discretization of the $K-,L-,S-$steps to
facilitate the discussion in the remaining part of the paper.

Input : numerical solution at $t^{n}:$ rank $r_{n}$ matrix $\hat{X}^{n}$ in
its SVD form $U^{n}\Sigma^{n}(V^{n})^{T}.$

Output : numerical solution at $t^{n+1}:$ rank $r_{n+1}$ matrix
$\hat{X}^{n+1}$ in its SVD form $U^{n+1}\Sigma^{n+1}(V^{n+1})^{T}.$

Parameter : time step $\Delta t$, error tolerance $\epsilon_{2}$

1 (Prediction). $K-$step and $L-$step integrating from $t^{n}$ to $t^{n+1}$.

2Solve  $K^{n+1}-K^{n}=\Delta tF(K^{n+1}(V^{n})^{T},t)V^{n},\quad
K^{n}=U^{n}\Sigma^{n},$ to obtain $K^{n+1},$ and $[\tilde{U},\sim,\sim]={\sf
qr}([U^{n},K^{n+1}]).$

3Solve  $L^{n+1}-L^{n}=\Delta tF(U^{n}(L^{n+1})^{T},t)^{T}U^{n},\quad
L^{n}=V^{n}\Sigma^{n},$ to obtain $L^{n+1},$ and $[\tilde{V},\sim,\sim]={\sf
qr}([V^{n},L^{n+1}]).$

4(Galerkin Evolution). $S-$step: solve for $S^{n+1}$ from
$S^{n+1}-S^{n}=\Delta
t\tilde{U}^{T}F(\tilde{U}S^{n+1}\tilde{V}^{T},t)\tilde{V},\quad
S^{n}=\tilde{U}^{T}U^{n}\Sigma^{n}(V^{n})^{T}\tilde{V},$ to obtain $S^{n+1}.$

5(Truncation).
$\hat{X}^{n+1}=\tilde{U}\mathcal{T}_{\epsilon_{2}}^{svd}(S^{n+1})\tilde{V}^{T}=U^{n+1}\Sigma^{n+1}(V^{n+1})^{T}.$

Algorithm 3 Rank adaptive BUG integrator using implicit Euler [2].

##  3 Numerical method

We now present the proposed numerical methods, the Merge method and its
adapted version. We first describe the algorithms, and then discuss the
properties and the rationale behind the design of the schemes.

###  3.1 The Merge method

In the Merge method, we merge the column and row spaces from the explicit step
truncation method and the spaces from the $K-$ and $L-$ steps in the BUG
method. As mentioned above the method consists of three stages.

Prediction: Our first ingredient in the method is to predict the column and
row space using the explicit Euler scheme in combination with the BUG
prediction spaces from Algorithm 3. Precisely, the predicted column and row
spaces are defined as the column and row spaces of the collection of matrices

$[\hat{X}^{n},\,\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})),\,K^{n+1}(L^{n+1})^{T}].$

Here the truncation level $\epsilon_{1}$ can be set to zero (no truncation
used) or chosen as $\epsilon_{1}=C_{1}\Delta t$. Both choices gives first
order accurate solutions but depending on the form of (1) the added cost from
the $\epsilon_{1}=C_{1}\Delta t$ truncation may be offset by the reduction in
the dimension of the predicted row and column spaces. To complete the next
step we first orthogonalize the predicted column and row spaces. That is, we
use column pivoted QR to find orthogonal matrices
$\tilde{U}\in\mathbb{R}^{m_{1}\times
s_{1}},\tilde{V}\in\mathbb{R}^{m_{2}\times s_{2}}$ spanning these spaces. Here
$s_{1},s_{2}$ is bounded from above by $2r_{n}(s+1)+r_{G},$ where $r_{G}$ is
the rank of the source $G(t^{n}).$

Galerkin evolution: The Galerkin evolution step can be understood as to find
an approximation $\hat{X}^{n+1,pre}\in W_{\tilde{U},\tilde{V}},$ such that

$\langle\hat{X}^{n+1,pre},A\rangle=\langle\hat{X}^{n}+\Delta
tF(\hat{X}^{n+1,pre},t^{n+1}),A\rangle,\qquad\forall A\in
W_{\tilde{U},\tilde{V}},$ (4)

where $W_{\tilde{U},\tilde{V}}=\\{A\in\mathbb{R}^{m_{1}\times
m_{2}}:A=\tilde{U}\Sigma\tilde{V}^{T},\textrm{with
}\Sigma\in\mathbb{R}^{s_{1}\times s_{2}}\\}$ denote all size $m_{1}\times
m_{2}$ spaces with column and row spaces as $\tilde{U}$ and $\tilde{V}.$

To solve this problem we let
$\hat{X}^{n+1,pre}=\tilde{U}\tilde{\Sigma}\tilde{V}^{T},A=\tilde{U}\Sigma^{*}\tilde{V}^{T}.$
Then by (2), we get

$\langle\tilde{\Sigma},\Sigma^{*}\rangle=\langle\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1})))\tilde{V},\Sigma^{*}\rangle,\qquad\forall\Sigma^{*}\in\mathbb{R}^{s_{1}\times
s_{2}}.$

This means

$\displaystyle\tilde{\Sigma}$ $\displaystyle=$
$\displaystyle\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1})))\tilde{V}$ (5)
$\displaystyle=$
$\displaystyle\tilde{U}^{T}U^{n}\Sigma^{n}(V^{n})^{T}\tilde{V}^{T}+\Delta
t\sum_{j=1}^{s}\tilde{U}^{T}A_{j}\tilde{U}\tilde{\Sigma}\tilde{V}^{T}B_{j}^{T}\tilde{V}+\Delta
t\tilde{U}^{T}U^{n}_{G}\Sigma^{n}_{G}(V^{n}_{G})^{T}\tilde{V}.$

The equation (5) is a linear matrix equation (generalized Sylvester equation)
for the unknown $\tilde{\Sigma}.$ Since the problems we are considering have
low rank solutions, we expect the dimensions of the matrix $\tilde{\Sigma}$ to
be small. To solve this equation one can either use a direct solves or an
iterative solver (e.g. GMRES) [23].

Truncation: The truncation step uses the truncated SVD of
$\hat{X}^{n+1}=\mathcal{T}_{\epsilon_{2}}^{svd}(\tilde{U}\tilde{\Sigma}\tilde{V}^{T})=\tilde{U}\mathcal{T}_{\epsilon_{2}}^{svd}(\tilde{\Sigma})\tilde{V}^{T}.$
Here the truncation level is chosen as $\epsilon_{2}=C_{2}\Delta t^{2}$ to
ensure accuracy.

Algorithm 4 summarizes the Merge method.

1

Input : numerical solution at $t^{n}:$ rank $r_{n}$ matrix $\hat{X}^{n}$ in
its SVD form $\hat{U}^{n}\hat{\Sigma}^{n}(\hat{V}^{n})^{T}$.

Output : numerical solution at $t^{n+1}:$ rank $r_{n+1}$ matrix
$\hat{X}^{n+1}$ in its SVD form $U^{n+1}\Sigma^{n+1}(V^{n+1})^{T}.$

Parameters : truncation tolerances $\epsilon_{1},\epsilon_{2}$

2 (Merge Prediction). Compute the truncated SVD of the right hand side:
$[\mathcal{U},\mathcal{S},\mathcal{V}]=\mathcal{T}^{sum}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})),$.

3Compute the BUG prediction spaces $K^{n+1},L^{n+1}$ according to Algorithm 3.

4Orthogonalize the merged spaces by column pivoted QR to get the prediction
spaces $[\tilde{U},\sim,\sim]={\sf qr}([[\hat{U}^{n},\mathcal{U},K^{n+1}]),$
$[\tilde{V},\sim,\sim]={\sf qr}([\hat{V}^{n},\mathcal{V}^{n},L^{n+1}]).$

5(Galerkin Evolution). Find $\tilde{\Sigma}$ by solving
$\tilde{\Sigma}=\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1}))\tilde{V}.$

(Truncation).
$\hat{X}^{n+1}=\hat{U}^{n+1}\hat{\Sigma}^{n+1}(\hat{V}^{n+1})^{T}=\tilde{U}\mathcal{T}_{\epsilon_{2}}^{svd}(\tilde{\Sigma})\tilde{V}^{T}.$

Algorithm 4 Merge method $t^{n}\rightarrow t^{n+1}$

We now provide estimates of the computational cost of the different parts of
the computation. Here we assume $m_{1}=m_{2}=m$, $s$ terms in differential
equation and rank $r$.

  * •

An evaluation of $AU,(BV)^{T}$ costs $\mathcal{O}(mr)$ for each of $AU$, $BV$,
assuming $A$ and $B$ are sparse.

  * •

The truncation of the right hand side $\mathcal{T}^{sum}_{\epsilon_{1}}(F)$
has two components. Two QR orthogonalizations, each costs
$\mathcal{O}(2(m(rs)^{2}-\frac{(rs)^{2}}{3})$ and one SVD of the core dense
matrix, this costs $\mathcal{O}((rs)^{3})$.

  * •

The $K$ and $L$ solves used to compute the BUG spaces also require the
solution of matrix equations but the dimension of these matrix solves are for
a $m\times r$ matrix rather than a $sr\times sr$ matrix. For these equations
we have found that applying GMRES in vectorized form works well. It should be
noted that this is typically the most expensive part of the solvers.

  * •

The cost of solving for the Galerkin core matrix $\tilde{\Sigma}$ in

$\tilde{\Sigma}-\Delta
t\tilde{U}^{T}F(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1})\tilde{V}=\tilde{U}^{T}\hat{X}^{n}\tilde{V},$

depends on the algorithm used. For example, if GMRES is applied on the
vectorized version of this equation each evaluation of the left hand side
costs $\mathcal{O}(m(rs)^{2})$ if $F$ is a general function. If $F$ is in the
form described in equation (1) it is possible to pre-compute dense $rs\times
rs$ matrices (at a cost of $\mathcal{O}(ms(rs)^{2})$) and bring down the cost
of the evaluation of the left hand side to $\mathcal{O}(s(rs)^{3})$ per
evaluation.

Neglecting the cost associated with the growth of the GMRES Krylov subspace
this cost would then be multiplied with the number of iterations needed to
converge (which is difficult to estimate). We have found that when the
governing equation has elliptic terms it also works well to reformulate the
Galerkin equation into a fixed point iteration

$A_{2}^{-1}(I-\tilde{U}^{T}\Delta
tA_{1}\tilde{U})\tilde{\Sigma}^{k}-\tilde{\Sigma}^{k}(\Delta
tB_{2}\tilde{V})^{T}\tilde{V}B^{-1}_{1}=P(\tilde{\Sigma}^{k-1},\hat{X}^{n},\tilde{U},\tilde{V}).$

In the iteration we use a dense Sylvester solve with cost
$\mathcal{O}((rs)^{3})$ for each iteration. When the matrices $A_{1}$ and
$B_{2}$ correspond to approximations to second derivatives in the 1 and 2
direction and $A_{2}$ and $B_{1}$ are identity or diagonal positive definite
matrices we find that this iteration converges to machine precision in a
handful of iterations.

Finally, we would like to comment on the rank evolution of the numerical
solution. We find in our numerical experiments that the low rank method
generally track the rank growth of the implicit Euler scheme well, i.e. the
rank of $\hat{X}^{n}$ is on par with the rank of the implicit Euler solution
at $t^{n}.$ However, due to the merge in the prediction step, the predicted
rank is larger than the rank of the BUG space (which is equal to $2r$). This
do translate to larger computational cost in the Galerkin step compared to the
BUG solver. However we note that the upper bound in the BUG space $2r$ can be
insufficient to capture rapidly changing initial layers in the PDEs as pointed
out in [19] (where it was observed a larger than actual rank needs to be
imposed for the numerical initial solution). For the numerical experiments we
performed in this paper, we find that by merging with the spaces generated by
the explicit schemes, this issue seems to be addressed and we do not need to
impose large artificial rank for the numerical initial condition.

###  3.2 The Merge-adapt method

As will be shown in Lemma 3.3, the column and row spaces generated from the
current solution $\hat{U}^{n},\hat{V}^{n}$ combined with $\mathcal{U}$ and
$\mathcal{V}$ in
$[\mathcal{U},\mathcal{S},\mathcal{V}]=\mathcal{T}^{sum}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})),$
will be sufficiently accurate to yield a first order accurate solution as long
as the time step is sufficiently small compared to the Lipschitz constant $L$
of $\|F\|$. This time step restriction can be overcome by merging with the BUG
spaces $K^{n+1},L^{n+1}$, which are computed useing an implicit solver and can
handle the stiffness as shown in the Merge method in Algorithm 4. However,
since the computation of $K^{n+1}$ and $L^{n+1}$ requires a linear solve, it
is preferable to only add those when necessary. This motivates the design of
the Merge-adapt method, which is described in Algorithm 5. In Steps 1-4 in
Algorithm 5, we perform a calculation based purely on column and row spaces
generated by the explicit scheme for the term
$\mathcal{T}^{sum}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})).$ In Step 5, we do a
residual check, and if this fails, the method will fall back to the Merge
method as shown in Steps 6-8. Otherwise, we proceed to the next time step.

We, heuristically, argue that for this residual check can be used to robustly
maintain first order accuracy. For simplicity, assume that the differential
equation we would like to solve is

$\frac{d}{dt}X(t)=AX(t)B^{T}+G(t).$

Let $X^{n+1}$ be the solution obtained by the classic implicit Euler method,
i.e. it solves

${\sf vec}(X^{n+1})=\underbrace{(I-\Delta t(B\otimes A))^{-1}}_{C}{\sf
vec}(X^{n})+\Delta t\,{\sf vec}(G(t_{n+1})).$

The low rank solution is not, in general, expected to solve this equation but
we may introduce the residual $\hat{R}^{n+1}$ so that

${\sf vec}(\hat{X}^{n+1})=\underbrace{(I-\Delta t(B\otimes A))^{-1}}_{C}{\sf
vec}(\hat{X}^{n})+\Delta t\,{\sf vec}(G(t_{n+1}))+{\sf vec}(\hat{R}^{n+1}).$

Combining the equations we find

${\sf vec}(\hat{X}^{n+1})-{\sf vec}(X^{n+1})=C[{\sf vec}(\hat{X}^{n})-{\sf
vec}(X^{n})]+{\sf vec}(\hat{R}^{n+1}).$

Now if we assume $\|C\|\leq 1$ and the residual check pass, then the error
between the implicit Euler solution and the implicit adaptive low rank
solution $e^{n}={\sf vec}(\hat{X}^{n})-{\sf vec}(X^{n})$ satisfies
$\|e^{n+1}\|\leq\|e^{n}\|+C_{2}\Delta t^{2}.$ This is enough to guarantee that
the low rank solution is convergent because the implicit Euler solution is
first order in time. It is reasonable to assume that the condition $\|C\|\leq
1$ is satisfied for most discretizations of diffusion or advection-diffusion
equations [9].

Input : numerical solution at $t^{n}:$ rank $r_{n}$ matrix $\hat{X}^{n}$ in
its SVD form $\hat{U}^{n}\hat{\Sigma}^{n}(\hat{V}^{n})^{T}$.

Output : numerical solution at $t^{n+1}:$ rank $r_{n+1}$ matrix
$\hat{X}^{n+1}$ in its SVD form $U^{n+1}\Sigma^{n+1}(V^{n+1})^{T}.$

Parameters : truncation tolerances $\epsilon_{1},\epsilon_{2}$

1 (Cheap Prediction). Compute a first order prediction of the column and row
spaces of $\hat{X}^{n+1}$. Compute the truncated SVD of the right hand side:
$[\mathcal{U},\mathcal{S},\mathcal{V}]=\mathcal{T}^{sum}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n}))$.

2Orthogonalize by column pivoted QR to get the prediction spaces
$[\tilde{U},\sim,\sim]={\sf qr}([[\hat{U}^{n},\mathcal{U}]),$
$[\tilde{V},\sim,\sim]={\sf qr}([\hat{V}^{n},\mathcal{V}^{n}]).$

3(Galerkin Evolution). Find $\tilde{\Sigma}$ by solving
$\tilde{\Sigma}=\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1}))\tilde{V}.$

4(Truncation).
$\hat{X}^{n+1}=\hat{U}^{n+1}\hat{\Sigma}^{n+1}(\hat{V}^{n+1})^{T}=\tilde{U}\mathcal{T}_{\epsilon_{2}}^{svd}(\tilde{\Sigma})\tilde{V}^{T}.$

5(Residual Check). Compute $\hat{R}^{n+1}=\hat{X}^{n+1}-\hat{X}^{n}-\Delta
tF(\hat{X}^{n+1},t^{n+1})$. If $\|\hat{R}^{n+1}\|<\epsilon_{2}$ return the
solution $\hat{X}^{n+1}$. If not

6(Merge Predicion). Compute the BUG prediction spaces $K^{n+1},L^{n+1}$
according to Algorithm 3. Orthogonalize the merged spaces by column pivoted
QR: $[\tilde{U},\sim,\sim]={\sf qr}([[\hat{U}^{n},\mathcal{U},K^{n+1}]),$
$[\tilde{V},\sim,\sim]={\sf qr}([\hat{V}^{n},\mathcal{V}^{n},L^{n+1}]).$

7(Galerkin Evolution). Find $\tilde{\Sigma}$ by solving
$\tilde{\Sigma}=\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1}))\tilde{V}.$

8(Truncation).
$\hat{X}^{n+1}=\hat{U}^{n+1}\hat{\Sigma}^{n+1}(\hat{V}^{n+1})^{T}=\tilde{U}\mathcal{T}_{\epsilon_{2}}^{svd}(\tilde{\Sigma})\tilde{V}^{T}.$

Algorithm 5 Merge-adapt method $t^{n}\rightarrow t^{n+1}$

We note that the acceleration of Algorithm 5 over Algorithm 4 is based on the
assumption that residual check passes in Step 5 and the BUG spaces are not
computed. If the residual check fails, extra computation is performed in Steps
3-5 which may incur more computational cost. We observe that, for moderately
stiff problems with small $s$, the Merge-adapt method has computational
advantages. More discussion is provided in the numerical experiment section.

###  3.3 Analysis

In this section, we analyze properties of the proposed methods.

####  3.3.1 Stability

First, we consider conservative or dissipative systems for which $\langle
F(X,t),X\rangle\leq 0,\forall t,X.$ Then it follows that

$\frac{d}{dt}\|X\|^{2}=\langle F(X,t),X\rangle\leq 0,$

that is, the energy ($L^{2}$ norm) is monotonically decreasing.

######  Theorem 3.1.

If we have $\langle F(X,t),X\rangle\leq 0,\forall t,X,$ then the numerical
solutions from Algorithm 4 or 5 satisfy

$\|\hat{X}^{n+1}\|\leq\|\hat{X}^{n}\|.$

###### Proof.

For Algorithm 4, first note that
$\|X^{n+1,pre}\|=\|\tilde{U}\tilde{\Sigma}\tilde{V}^{T}\|=\|\tilde{\Sigma}\|.$
By (5) and (2),

$\displaystyle\|\tilde{\Sigma}\|^{2}=\langle\tilde{\Sigma},\tilde{\Sigma}\rangle=\langle\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1}))\tilde{V},\tilde{\Sigma}\rangle$
$\displaystyle=\langle\hat{X}^{n}+\Delta
tF(\tilde{U}\tilde{\Sigma}\tilde{V}^{T},t^{n+1}),\tilde{U}\tilde{\Sigma}\tilde{V}^{T}\rangle$
$\displaystyle\leq\langle\hat{X}^{n},\tilde{U}\tilde{\Sigma}\tilde{V}^{T}\rangle=\langle\hat{X}^{n},X^{n+1,pre}\rangle.$

Therefore, by Cauchy-Schwarz inequality,
$\|X^{n+1,pre}\|^{2}\leq\|\hat{X}^{n}\|\|X^{n+1,pre}\|,$ which gives
$\|X^{n+1,pre}\|\leq\|\hat{X}^{n}\|.$ Finally, because of the property of the
truncated SVD,
$\|\hat{X}^{n+1}\|=\|\mathcal{T}_{\epsilon_{2}}^{svd}(X^{n+1,pre})\|\leq\|X^{n+1,pre}\|\leq\|\hat{X}^{n}\|.$
The proof for Algorithm 5 is the same and is omitted. ∎

We can easily generalize this result to semi-bounded operator [9].
Specifically, if $\langle F(X,t),X\rangle\leq\alpha\|X\|^{2},\forall\,t,X,$
then

$\frac{d}{dt}\|X\|^{2}=\langle F(X,t),X\rangle\leq\alpha\|X\|^{2},$

which implies $\|X(t)\|\leq e^{\alpha t}\|X(0)\|.$

######  Theorem 3.2 (Stability for semi-bounded operator).

If we have $\langle F(X,t),X\rangle\leq\alpha\|X\|^{2},\forall\,t,X$ then the
numerical solutions from Algorithm 4 or 5 satisfy

$\|\hat{X}^{n+1}\|\leq e^{\alpha\Delta
t}\|\hat{X}^{n}\|,\quad\|\hat{X}^{n}\|\leq
e^{\alpha\,t^{n}}\|\hat{X}^{0}\|,\quad\textrm{if}\,\,\alpha\Delta t\leq 1.$

###### Proof.

The proof is similar to the proof of Theorem 3.1, so we only highlight the
difference. We have

$\langle\hat{X}^{n+1},\hat{X}^{n+1}\rangle\leq\langle\hat{X}^{n},\hat{X}^{n+1}\rangle+\Delta
t\alpha\|\hat{X}^{n+1}\|^{2},$

which implies

$\|\hat{X}^{n+1}\|\leq\frac{1}{1-\alpha\Delta t}\|\hat{X}^{n}\|,$

if $\alpha\Delta t\leq 1.$ The theorem follows using a similar arguments as in
Theorem 3.1. ∎

####  3.3.2 Convergence

A convergence estimate has been shown in Theorem 2 of [2] for the rank
adaptive BUG scheme. Assuming the Lipschitz continuity and boundedness of the
operator $F,$ the authors show that the BUG scheme has numerical error bounded
by the sum of the initial numerical error, the tangent projection error,
$\epsilon_{2}$ and first order in time error. The proof is based on the time
continuous version of the BUG schemes. We want to point out that for PDE
applications, in general the (differential) operators are not Lipschitz
bounded. Nevertheless, it is still of theoretical interests to investigate the
convergence properties of the schemes under such assumptions.

For simplicity, below we will estimate the local truncation error of the first
order implicit schemes with column and row space spanned by the cheap
prediction space, i.e. we are investigating the Steps 1-4 in Algorithm 5.
Since the Merge method uses a space that is union of the cheap prediction (as
shown below) and the BUG space, the local truncation error of Algorithm 4 will
be also be upper bounded by the estimate below.

######  Lemma 3.3 (Local truncation error with cheap prediction space).

Suppose $F$ is Lipschitz-continuous in both variables and bounded, i.e. there
exists a constant $L$ such that $\|F(Y,s)-F(Z,t)\|\leq L\|Y-Z\|+L|s-t|$ and
$\|F(X,t)\|\leq B$. If we denote $\mathcal{X}$ as the exact solution to (1)
with initial condition $\hat{X}^{n}=U^{n}\Sigma^{n}(V^{n})^{T}$ from $t^{n}$
to $t^{n+1}$, $\hat{X}^{n+1}$ the numerical solution obtained from from
Algorithm 5, Steps 1-4 with the same initial condition $\hat{X}^{n},$ then for
sufficiently small $\Delta t,$ we have the following error bound

$\|\mathcal{X}-\hat{X}^{n+1}\|\leq C\Delta t^{2}+2\Delta
t\epsilon_{1}+\epsilon_{2}.$ (6)

In (6), $C$ is a constant that only depends $L$ and $B$.

###### Proof.

First we prove (6) for $\epsilon_{1}=0$. Then
$\mathcal{T}^{sum}_{\epsilon_{1}=0}(F(\hat{X}^{n},t^{n}))=F(\hat{X}^{n},t^{n})$.
From (5), and left multiplying by $\tilde{U}$ and right multiplying by
$\tilde{V}^{T},$ we obtain

$\hat{X}^{n+1,pre}=\tilde{U}\tilde{\Sigma}\tilde{V}^{T}=\tilde{U}\tilde{U}^{T}(\hat{X}^{n}+\Delta
tF(\hat{X}^{n+1,pre},t^{n+1})))\tilde{V}\tilde{V}^{T}.$

Since the column space of $\hat{X}^{n}$ is a subset of the column space of
$\tilde{U},$ i.e. $R(\hat{X}^{n})\subseteq R(\tilde{U}),$ we have
$\tilde{U}\tilde{U}^{T}\hat{X}^{n}=\hat{X}^{n}.$ Similarly, because
$R((\hat{X}^{n})^{T})\subseteq R(\tilde{V}),$
$\hat{X}^{n}\tilde{V}\tilde{V^{T}}=\hat{X}^{n}.$ Therefore,
$\tilde{U}\tilde{U}^{T}\hat{X}^{n}\tilde{V}\tilde{V^{T}}=\hat{X}^{n}.$ By a
similar argument, $R(F(\hat{X}^{n},t^{n}))\subseteq R(\tilde{U})$ and
$R(F(\hat{X}^{n},t^{n})^{T})\subseteq R(\tilde{V})$ implies
$\tilde{U}\tilde{U}^{T}F(\hat{X}^{n},t^{n})\tilde{V}\tilde{V}^{T}=F(\hat{X}^{n},t^{n}).$
This gives

$\displaystyle\hat{X}^{n+1,pre}=\hat{X}^{n}+\Delta
t\tilde{U}\tilde{U}^{T}F(\hat{X}^{n+1,pre},t^{n+1})\tilde{V}\tilde{V}^{T}$
$\displaystyle=\hat{X}^{n}+\Delta
t\tilde{U}\tilde{U}^{T}F(\hat{X}^{n},t^{n})\tilde{V}\tilde{V}^{T}+\Delta
t\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}$
$\displaystyle=\hat{X}^{n}+\Delta tF(\hat{X}^{n},t^{n})+\Delta
t\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}.$
(7)

Therefore,

$\|\hat{X}^{n+1,pre}-\mathcal{X}\|\\\ =\|\hat{X}^{n}-\mathcal{X}+\Delta
tF(\hat{X}^{n},t^{n})+\Delta
t\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}\|.$
(8)

We have $\hat{X}^{n}-\mathcal{X}=-\Delta tF(X(t^{*}),t^{*}),$ where $t^{*}$ is
a point on the interval $[t^{n},t^{n+1}].$ Therefore,

$\displaystyle\|\hat{X}^{n}-\mathcal{X}+\Delta tF(\hat{X}^{n},t^{n})\|=\Delta
t\|-F(X(t^{*}),t^{*})+F(\hat{X}^{n},t^{n})\|$ $\displaystyle\leq\Delta
tL(\Delta t+\|X(t^{*})-\hat{X}^{n}\|)\leq\Delta t^{2}L(B+1).$ (9)

On the other hand,

$\displaystyle\|\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}\|$
$\displaystyle\leq\|\tilde{U}\tilde{U}^{T}\|_{2}\|F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n})\|\|\tilde{V}\tilde{V}^{T}\|_{2},$

where $\|\cdot\|_{2}$ is the matrix 2-norm. Since
$\|\tilde{U}\tilde{U}^{T}\|_{2}=\|\tilde{V}\tilde{V}^{T}\|_{2}=1$ and
$\|F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n})\|\leq L(\Delta
t+\|\hat{X}^{n}-\mathcal{X}\|+\|\hat{X}^{n+1,pre}-\mathcal{X}\|)\leq L(\Delta
t+B\Delta t+\|\hat{X}^{n+1,pre}-\mathcal{X}\|),$ we get

$\displaystyle\|\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n+1,pre},t^{n+1})-F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}\|\leq
L(\Delta t+B\Delta t+\|\hat{X}^{n+1,pre}-\mathcal{X}\|).$ (10)

Combining (8),(3.3.2),(10), we get

$\|\hat{X}^{n+1,pre}-\mathcal{X}\|\leq\Delta t^{2}L(2B+2)+L\Delta
t\|\hat{X}^{n+1,pre}-\mathcal{X}\|,$

Choosing $\Delta t$ small enough (i.e. $\Delta t\leq\frac{1}{2L}$), we proved

$\displaystyle\|\hat{X}^{n+1,pre}-\mathcal{X}\|\leq 4\Delta t^{2}L(B+1).$ (11)

Therefore, (6) follows by using
$\|\hat{X}^{n+1,pre}-\hat{X}^{n+1}\|=\|\hat{X}^{n+1,pre}-\mathcal{T}_{\epsilon_{2}}^{svd}(\hat{X}^{n+1,pre})\|\leq\epsilon_{2}.$

The proof for (6) with $\epsilon_{1}\neq 0$ largely follows the same process.
The only difference is that we no longer have
$\tilde{U}\tilde{U}^{T}F(\hat{X}^{n},t^{n})\tilde{V}\tilde{V}^{T}=F(\hat{X}^{n},t^{n})$
due to the truncation in the prediction step. Instead, we have
$\tilde{U}\tilde{U}^{T}\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n}))\tilde{V}\tilde{V}^{T}=\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})).$
In this case,

$\tilde{U}\tilde{U}^{T}F(\hat{X}^{n},t^{n})\tilde{V}\tilde{V}^{T}=\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n}))+\tilde{U}\tilde{U}^{T}(F(\hat{X}^{n},t^{n})-\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n})))\tilde{V}\tilde{V}^{T},$

and because $\|\tilde{U}\tilde{U}^{T}\|_{2}=\|\tilde{V}\tilde{V}^{T}\|_{2}=1,$
this leads to

$\|\tilde{U}\tilde{U}^{T}F(\hat{X}^{n},t^{n})\tilde{V}\tilde{V}^{T}-F(\hat{X}^{n},t^{n})\|\leq
2\|F(\hat{X}^{n},t^{n})-\mathcal{T}_{\epsilon_{1}}(F(\hat{X}^{n},t^{n}))\|=2\epsilon_{1}.$

By a similar argument, (6) follows. The details are omitted for brevity. ∎

The numerical convergence follows from Lemma 3.3 by a standard result for the
convergence of one-step method for ODE and is skipped for brevity. Now we
would like to comment on the results in this Lemma to understand property of
the scheme. As we can see from the results in Section 3.3.1, the schemes are
unconditionally stable due to the Galerkin step. Also, in Lemma 3.3, the
tangent projection error is absent because the DLRA approximation is not used.
However, the time step is limited by the Lipschitz constant and also the
constant in the estimate $C$ depends on the Lipschitz constant. For standard
error estimates for implicit schemes for stiff problems, the concept of
B-convergence [5] can remove such restrictions with the help of one-sided
Lipschitz continuity. Unfortunately, such results are not available in our
case.

##  4 Numerical results

We now present numerical results that illustrate the features of the methods
described above. In this section we will refer to the different methods in
figures and tables according to the following naming convention. We will
denote the Merge method by M, and the Merge-adapt method by MA. For many
problems we will also compare to the classic implicit Euler discretization as
IE or Implicit Euler. When computing errors we will compare to a solution
computed using Matlab’s built in ODE15s or ODE45 solvers with relative and
absolute tolerance set to $10^{-12}$. We always report relative errors in the
Frobenius norm. We note that in all of the examples below, the BUG method does
not converge since all or part of the right hand side of (1) is outside the
tangent projection space. We also note that in all the examples below we use
$\epsilon_{1}=0$ but that the results are very similar for
$\epsilon_{1}=\Delta t$.

The main PDE examples used in this section are simulations of solutions to the
equation

$\frac{\partial\rho}{\partial t}+r_{1}(x_{1})\frac{\partial\rho}{\partial
x_{2}}+r_{2}(x_{2})\frac{\partial\rho}{\partial x_{1}}=\\\
b_{1}(x_{2})\frac{\partial}{\partial
x_{1}}\left[a_{1}(x_{1})\frac{\partial\rho}{\partial
x_{1}}\right]+b_{2}(x_{2})\frac{\partial^{2}\left[a_{2}(x_{1})\rho\right]}{\partial
x_{1}\partial
x_{2}}+a_{3}(x_{1})\frac{\partial^{2}\left[b_{3}(x_{2})\rho\right]}{\partial
x_{1}\partial x_{2}}+a_{4}(x)\frac{\partial}{\partial
y}\left[b_{4}(y)\frac{\partial\rho}{\partial y}\right].$ (12)

Here $\rho=\rho(t,x_{1},x_{2})$ and we always use homogenous Dirichlet
boundary conditions. Equation (12) is discretized on the domain
$(x_{1},x_{2})\in[-1,1]^{2},$ using second order accurate finite differences.
Consider the grid $(x_{1,i},x_{2,j})=(-1+ih_{1},-1+jh_{2})$, with
$h_{1}=2/(m_{1}+1),$ $h_{2}=2/(m_{2}+1).$ The solution
$\rho(t_{n},x_{1,i},x_{2,j})$ is then approximated by the grid function
(matrix) $\hat{X}^{n}_{i,j}$ and derivatives are approximated as follow:

$\frac{\partial}{\partial
x_{1}}\left[a_{1}(x_{1})\frac{\partial\rho(t_{n},x_{1},x_{2})}{\partial
x_{1}}\right]\Bigg{|}_{x_{1,i},x_{2,j}}\approx
D_{+}^{1}\left[\frac{a_{1}(x_{1,i})+a_{1}(x_{1,i-1})}{2}\right]D_{-}^{1}\,X_{i,j}^{n},$
$\frac{\partial}{\partial
x_{1}}\left[a_{2}(x_{1})\frac{\partial\rho(t_{n},x_{1},x_{1})}{\partial
x_{2}}\right]\Bigg{|}_{x_{1,i},x_{2,j}}\approx
D_{0}^{1}a_{1}(x_{1,i})D_{0}^{2}\,X_{i,j}^{n}.$

with the remaining two terms discretized in the same way. Here the difference
operators in the 1-direction are defined

$2h_{1}D_{0}^{1}w_{i,j}\equiv w_{i+1,j}-w_{i-1,j},\ \
h_{1}D_{+}^{1}w_{i,j}\equiv w_{i+1,j}-w_{i,j},\ \ h_{1}D_{-}^{1}w_{i,j}\equiv
w_{i,j}-w_{i-1,j},$

with the operators in the 2-direction defined analogously.

###  4.1 Rotation with anisotropic diffusion

In this experiment we consider solid body rotation, $r_{1}(x_{1})=x_{1}$
$r_{2}(x_{2})=-x_{2}$, with anisotropic diffusion. The diffusion coefficients
are:

$\displaystyle a_{1}(x_{1})=a_{4}(x_{1})=\sqrt{\mu}(1+0.1\sin(\pi x_{1})),$
$\displaystyle a_{2}(x_{1})=\sqrt{\mu}(0.15+0.1\sin(\pi x_{1})),$
$\displaystyle a_{3}(x_{1})=\sqrt{\mu}(0.15+0.1\cos(\pi x_{1})),$
$\displaystyle b_{1}(x_{2})=b_{4}(x_{2})=\sqrt{\mu}(1+0.1\cos(\pi x_{2})),$
$\displaystyle b_{2}(x_{2})=\sqrt{\mu}(0.15+0.1\cos(\pi x_{2})),$
$\displaystyle b_{3}(x_{2})=\sqrt{\mu}(0.15+0.1\sin(\pi x_{2})),$

and $\mu=10^{-3}$. We start from the rank-1 initial data

$\rho(0,x_{1},x_{2})=e^{-\left(\frac{x_{1}}{0.3}\right)^{2}}e^{-\left(\frac{x_{2}}{0.1}\right)^{2}}$

and evolve the solution until time $t=\pi$ using a timestep $\Delta
t=\pi/n_{T}$, with $40,80,160,320$. We carry out the computation on three
different grids with $m_{1}=m_{2}=99,199,799$.

In Figure 1 we display 10 equidistant contours between 0.1 and 0.9 of the
solution. Results using the Merge method (left), the classic implicit Euler
(middle) and the BUG method (right) are presented. The solutions are displayed
at time $0.5\pi$ when they have rotated from being horizontal to vertical. The
BUG method does not see the rotation (it is outside the tangent space) and
remains stationary. The solution computed using the Merge method rotates as
expected and is almost identical to the classical implicit Euler solution.

Figure 1: Displayed are 10 equidistant contours between 0.1 and 0.9 of the solution computed using the Merge method (left), the classic implicit Euler (middle) and the BUG method (right). The solutions are displayed at time $0.5\pi$ when they have rotated from being horizontal to vertical. The BUG method does not see the rotation (it is outside the tangent space) and remains stationary. The solution computed using the Merge method rotates as expected and is almost identical to the classical implicit Euler solution.  $n_{T}$ |  M |  MA |  IE |  $\Delta t/h$ : $\mu\Delta t/h^{2}$
---|---|---|---|---
40 |  1.60(-1) |  1.60(-1) F = 5 |  1.60(-1) |  3.9 : 0.2
80 |  1.01(-1) [0.65] |  1.01(-1) [0.65] F = 3 |  1.01(-1) [0.65] |  2 : 0.098
160 |  6.01(-2) [0.75] |  6.01(-2) [0.75] F = 4 |  6.01(-2) [0.75] |  0.98 : 0.049
320 |  3.36(-2) [0.83] |  3.36(-2) [0.83] F = 0 |  3.36(-2) [0.83] |  0.49 : 0.025
40 |  1.60(-1) |  1.60(-1) F = 8 |  1.60(-1) |  7.9 : 0.79
80 |  1.02(-1) [0.65] |  1.02(-1) [0.65] F = 6 |  1.02(-1) [0.65] |  3.9 : 0.39
160 |  6.07(-2) [0.75] |  6.07(-2) [0.75] F = 4 |  6.07(-2) [0.75] |  2 : 0.2
320 |  3.40(-2) [0.83] |  3.40(-2) [0.83] F = 1 |  3.40(-2) [0.83] |  0.98 : 0.098
40 |  1.61(-1) |  1.61(-1) F = 13 |  1.61(-1) |  31 : 13
80 |  1.02(-1) [0.65] |  1.02(-1) [0.65] F = 11 |  1.02(-1) [0.65] |  16 : 6.3
160 |  6.09(-2) [0.75] |  6.09(-2) [0.75] F = 11 |  6.09(-2) [0.75] |  7.9 : 3.1
320 |  3.41(-2) [0.83] |  3.41(-2) [0.83] F = 9 |  3.41(-2) [0.83] |  3.9 : 1.6
Table 1: Solid body rotation with anisotropic diffusion. Displayed are the
errors (here 1.3(-1) means $1.3\cdot 10^{-1}$) for different timesteps along
with estimated rates of convergence (in brackets). The numbers on the far
right are the “hyperbolic and parabolic CFL numbers”. The top box is for
$m_{1}=m_{2}=99$, the middle for $m_{1}=m_{2}=199$ and the bottom for
$m_{1}=m_{2}=799$. The F indicates how many times the Merge-adapt needed to
add the BUG spaces.

In Table 1 we display the errors and rates of convergence for the Merge,
Merge-adapt and the implicit Euler method. All three methods behave almost
identically (the errors and rates of convergence only differ in the third or
fourth digit). However, the time it takes to compute the solution is quite
different. In Table 2 we display the CPU times for the different methods. The
computations were performed in Matlab on a MacBook Pro with an M2 chip and
16GB RAM. Reasonable effort went into optimizing the codes but the timings
should be seen as indicative rather than decisive. As can be seen the growth
in CPU time as a function of number of degrees of freedom is very mild for the
Merge and Merge-adapt methods in comparison to the classic implicit Euler
method. We further note that the Merge-adapt method has comparable cost with
the Merge method in this example because the equation involves relatively
complex dynamics, which means the explicit prediction is quite costly because
of the large number $s.$

$m_{1}=m_{2}$ |  M time$[s]$ |  MA time$[s]$ |  IE time$[s]$
---|---|---|---
99 |  0.9 |  0.9 |  5.0
199 |  1.4 |  1.5 |  27.1
799 |  12.8 |  11.4 |  926.2
Table 2: Displayed are the CPU times in seconds for the rotation with
anisotropic diffusion example. The results are all for $n_{T}=320$. As can be
seen the growth in CPU time as a function of number of degrees of freedom is
very mild for the Merge and Merge-adapt methods in comparison to the classic
implicit Euler method.

We also track the rank of the solution as a function of time. In Figure 2 we
display the rank for the different methods as a function of time. Once the
solution becomes more diagonal at time $\pi/4$ the rank is decreased and then
increased again as the rotation continues. For the Merge and Merge-adapt
method we use the truncation tolerance $\epsilon_{2}=\Delta t^{2}$ and in
order to compare with the classic implicit Euler we compute the SVD and count
the rank of that solution using the same threshold. We observe the ranks from
the low rank methods are slightly smaller than the implicit Euler method.

Figure 2: Displayed is the rank for the Merge, Merge-adapt and Implicit Euler
method for the problem with rotation and anisotropic diffusion. Here “Fail for
MA” is the indicator where the BUG space is needed for the Merge-adapt method.
For the Implicit Euler method we constantly use the threshold $\Delta t^{2}$
when computing the rank via the truncated SVD. All the computations are done
with $m_{1}=m_{2}=799$ and $n_{T}=320.$

###  4.2 Anisotropic diffusion

In this experiment we consider anisotropic diffusion without rotation
($r_{1}=r_{2}=0$). We take the diffusion coefficients to be constants with
$a_{1}=a_{4}=b_{1}=b_{4}=1$ and $a_{2}=a_{3}=b_{2}=b_{3}=0.3$. Starting from
the rank-1 initial data

$\rho(0,x_{1},x_{2})=\sin(\pi x_{1})\sin(\pi x_{2}),$

we evolve the solution until time $t=0.5$ using a base timestep of $\Delta
t=0.5/n_{T}$, with $40,80,160,320,640,1280$. We carry out the computation on
three different grids with $m_{1}=m_{2}=99,399,799$.

$n_{T}$ |  M |  MA |  IE |  $\Delta t/h$ : $\Delta t/h^{2}$
---|---|---|---|---
40 |  8.65(-2) |  8.65(-2) F = 37 |  9.31(-2) |  0.62 : 31
80 |  2.84(-2) [1.60] |  2.84(-2) [1.60] F = 70 |  4.39(-2) [1.08] |  0.31 : 16
160 |  9.94(-3) [1.51] |  9.94(-3) [1.51] F = 160 |  2.13(-2) [1.03] |  0.16 : 7.8
320 |  5.89(-3) [0.75] |  5.89(-3) [0.75] F = 296 |  1.05(-2) [1.01] |  0.078 : 3.9
640 |  3.78(-3) [0.63] |  3.78(-3) [0.63] F = 595 |  5.22(-3) [1.01] |  0.039 : 2
1280 |  2.19(-3) [0.78] |  2.19(-3) [0.78] F = 989 |  2.58(-3) [1.01] |  0.02 : 0.98
40 |  1.13(-1) |  1.13(-1) F = 40 |  9.30(-2) |  1.2 : 120
80 |  4.23(-2) [1.42] |  4.23(-2) [1.42] F = 80 |  4.38(-2) [1.08] |  0.62 : 62
160 |  1.33(-2) [1.66] |  1.33(-2) [1.66] F = 160 |  2.13(-2) [1.03] |  0.31 : 31
320 |  4.64(-3) [1.52] |  4.64(-3) [1.52] F = 320 |  1.05(-2) [1.01] |  0.16 : 16
640 |  2.88(-3) [0.69] |  2.88(-3) [0.69] F = 640 |  5.22(-3) [1.01] |  0.078 : 7.8
1280 |  1.87(-3) [0.61] |  1.87(-3) [0.61] F = 1280 |  2.58(-3) [1.01] |  0.039 : 3.9
40 |  1.03(-1) |  1.03(-1) F = 40 |  9.30(-2) |  5 : 2000
80 |  6.52(-2) [0.66] |  6.52(-2) [0.66] F = 80 |  4.38(-2) [1.08] |  2.5 : 1000
160 |  2.70(-2) [1.27] |  2.70(-2) [1.27] F = 160 |  2.13(-2) [1.03] |  1.2 : 500
320 |  9.69(-3) [1.48] |  9.69(-3) [1.48] F = 320 |  1.05(-2) [1.02] |  0.62 : 250
640 |  2.97(-3) [1.70] |  2.97(-3) [1.70] F = 640 |  5.21(-3) [1.01] |  0.31 : 1.25
1280 |  1.01(-3) [1.55] |  1.01(-3) [1.55] F = 1280 |  2.57(-3) [1.01] |  0.16 : 62
Table 3: Anisotropic diffusion. Displayed are the errors for different
timesteps along with estimated rates of convergence (in brackets). The numbers
on the far right are the “hyperbolic and parabolic CFL numbers”. The top box
is for $m_{1}=m_{2}=99$, the middle for $m_{1}=m_{2}=199$ and the bottom for
$m_{1}=m_{2}=799$. The F indicates how many times the Merge-adapt needed to
add the BUG spaces.

Figure 3: Displayed is the rank for the Merge, Merge-adapt and Implicit Euler
method for the problem with anisotropic diffusion initial data $\sin(\pi
x_{1})\sin(\pi x_{2}),$ (left) and $\sin(2\pi x_{1})\sin(2\pi x_{2})$ (right)
and the solid body rotation problem (right). Here “Fail for MA” is the
indicator where the BUG space is needed for the Merge-adapt method. For the
Implicit Euler method we constantly use the threshold $\Delta t^{2}$ when
computing the rank via the truncated SVD. All the computations are done with
$m_{1}=m_{2}=799$ and $n_{T}=320.$

In Table 3 we display the errors and rates of convergence. Again, the methods
are comparable although it should be noted that the rates of convergence for
the classic implicit Euler method is more uniform than for the low rank
methods. Curiously the error levels are smaller for the low rank method. It
should also be noted that for this example the Merge-adapt method requires the
addition of the BUG spaces almost every time step due to the stiffness of the
problem. We also track the rank of the solution as a function of time, the
results can be found in Figure 3. The rank evolution is very similar for all
methods.

We also consider the more oscillatory initial data

$\rho(0,x_{1},x_{2})=\sin(2\pi x_{1})\sin(2\pi x_{2}),$

which we again evolve until time $t=0.5$ using a base timestep of $\Delta
t=0.5/n_{T}$, with $40,80,160,320,640,1280$. Again, we carry out the
computation on three different grids with $m_{1}=m_{2}=99,399,799$.

This initial data decays more rapidly in time and as a consequence the
residual norm becomes smaller and the Merge-adapt method does not have to add
the BUG spaces as often. The errors and rates of convergence can be found in
Table 4, and the rank of the solution as a function of time, the results can
be found in Figure 2. It can be observed that for this example the classic
implicit Euler method at first is a bit more accurate than the low rank
methods but they catch up as the timestep decreases.

$n_{T}$ |  M |  MA |  IE |  $\Delta t/h$ : $\Delta t/h^{2}$
---|---|---|---|---
40 |  4.26(-1) |  4.30(-1) F = 18 |  7.58(-2) |  0.62 : 31
80 |  2.07(-1) [1.03] |  2.08(-1) [1.04] F = 35 |  3.80(-2) [0.99] |  0.31 : 16
160 |  7.90(-2) [1.39] |  7.96(-2) [1.38] F = 61 |  1.91(-2) [0.99] |  0.16 : 7.8
320 |  2.36(-2) [1.74] |  2.42(-2) [1.71] F = 104 |  9.67(-3) [0.98] |  0.078 : 3.9
640 |  5.14(-3) [2.19] |  5.15(-3) [2.23] F = 188 |  4.97(-3) [0.95] |  0.039 : 2
1280 |  1.00(-3) [2.35] |  1.01(-3) [2.34] F = 334 |  2.68(-3) [0.89] |  0.02 : 0.98
40 |  4.87(-1) |  4.89(-1) F = 20 |  7.57(-2) |  1.2 : 125
80 |  2.62(-1) [0.89] |  2.63(-1) [0.89] F = 40 |  3.80(-2) [0.99] |  0.62 : 62
160 |  1.13(-1) [1.20] |  1.13(-1) [1.21] F = 81 |  1.90(-2) [0.99] |  0.31 : 31
320 |  3.89(-2) [1.54] |  3.89(-2) [1.54] F = 134 |  9.60(-3) [0.98] |  0.16 : 16
640 |  1.05(-2) [1.87] |  1.06(-2) [1.87] F = 277 |  4.87(-3) [0.97] |  0.078 : 7.8
1280 |  2.01(-3) [2.39] |  2.01(-3) [2.39] F = 521 |  2.52(-3) [0.95] |  0.039 : 3.9
40 |  5.23(-1) |  5.23(-1) F = 39 |  7.58(-2) |  5 : 2000
80 |  3.49(-1) [0.58] |  3.49(-1) [0.58] F = 75 |  3.80(-2) [0.99] |  2.5 : 1000
160 |  1.80(-1) [0.95] |  1.80(-1) [0.95] F = 132 |  1.91(-2) [0.99] |  1.2 : 500
320 |  7.71(-2) [1.22] |  7.71(-2) [1.22] F = 242 |  9.71(-3) [0.97] |  0.62 : 250
640 |  2.89(-2) [1.41] |  2.89(-2) [1.41] F = 479 |  5.04(-3) [0.94] |  0.31 : 125
1280 |  9.17(-3) [1.65] |  9.17(-3) [1.65] F = 1162 |  2.78(-3) [0.85] |  0.16 : 62
Table 4: Anisotropic diffusion for a higher frequency initial data. Displayed
are the errors for different timesteps along with estimated rates of
convergence (in brackets). The numbers on the far right are the “hyperbolic
and parabolic CFL numbers”. The top box is for $m_{1}=m_{2}=99$, the middle
for $m_{1}=m_{2}=199$ and the bottom for $m_{1}=m_{2}=799$. The F indicates
how many times the Merge-adapt needed to add the BUG spaces.

###  4.3 Solid body rotation with isotropic diffusion

$m_{1}=m_{2}$ |  time MA [s] |  time M [s] |  $\Delta t/h$ : $\mu\Delta t/h^{2}$
---|---|---|---
999 |  1.9 [0.1%] |  3.1 |  0.39 : 0.020
1999 |  3.5 [0.1%] |  6.0 |  0.79 : 0.79
3999 |  7.1 [0.4%] |  15.5 |  1.57 : 0.31
7999 |  52.8 [72%] |  46.2 |  3.14 : 1.26
Table 5: CPU-times for the problem with solid body rotation with isotropic
diffusion. The numbers in brackets are the percentages for how often MA has a
residual that necessitates the addition of the BUG space.

In this experiment we time a problem with solid body rotation and isotropic
diffusion, corresponding to, $r_{1}(x_{1})=x_{1}$ $r_{2}(x_{2})=-x_{2}$,
$a_{1}=b_{1}=a_{4}=b_{4}=10^{-4}$, and $a_{j}=b_{j}=0,\,j=2,\ldots,3$. We
focus on the comparison of the computational cost of Merge and Merge-adapt
schemes for large mesh size. We start from the rank-1 initial data

$\rho(0,x_{1},x_{2})=e^{-\left(\frac{x_{1}}{0.3}\right)^{2}}e^{-\left(\frac{x_{2}}{0.1}\right)^{2}}$

and evolve the solution until time $t=0.25\pi$ using a base timestep of
$\Delta t=0.25\pi/n_{T}$, with $n_{T}=1000$. We carry out the computation on
very fine grids with $m_{1}=m_{2}=999,1999,3999,7999$. Here the size of the
diffusion is chosen so that the ratios $\Delta t/h$ : $\mu\Delta t/h^{2}$
change from smaller than one to larger than one. We report the CPU times for
the M and MA methods in Table 5. As can be seen the MA is faster than the M
method when the number of timesteps where the BUG space is used is small. The
more intensive use of the BUG space appears to happen when the ratio
$\mu\Delta t/h^{2}$ increases above one.

###  4.4 Solid body rotation

In this experiment we consider a hyperbolic problem with pure solid body
rotation, corresponding to, $r_{1}(x_{1})=x_{1}$ $r_{2}(x_{2})=-x_{2}$ and
$a_{j}=b_{j}=0,\,j=1,\ldots,4$. We start from the rank-1 initial data

$\rho(0,x_{1},x_{2})=e^{-\left(\frac{x_{1}}{0.3}\right)^{2}}e^{-\left(\frac{x_{2}}{0.1}\right)^{2}}$

and evolve the solution until time $t=\pi$ using a base timestep of $\Delta
t=\pi/n_{T}$, with $n_{T}=40,80,160,320$. We carry out the computation on
three different grids with $m_{1}=m_{2}=99,199,799$. The errors are listed in
Table 6. Again similar conclusions can be drawn, namely all three methods give
comparable errors. The rank evolution is reported in Figure 4. In this
problem, the rank evolution exhibits periodic pattern in time due to the
solution structure. The rank of the low rank methods are lower than the
implicit Euler method.

$n_{T}$ |  M |  MA |  IE |  $\Delta t/h$
---|---|---|---|---
40 |  2.47(-1) |  2.47(-1) F = 11 |  2.51(-1) |  3.9
80 |  1.71(-1) [0.53] |  1.71(-1) [0.53] F = 19 |  1.73(-1) [0.53] |  2
160 |  1.10(-1) [0.63] |  1.10(-1) [0.63] F = 4 |  1.10(-1) [0.64] |  0.98
320 |  6.57(-2) [0.74] |  6.57(-2) [0.74] F = 4 |  6.60(-2) [0.74] |  0.49
40 |  2.50(-1) |  2.50(-1) F = 14 |  2.54(-1) |  7.9
80 |  1.74(-1) [0.52] |  1.74(-1) [0.52] F = 27 |  1.76(-1) [0.53] |  3.9
160 |  1.12(-1) [0.62] |  1.12(-1) [0.62] F = 32 |  1.13(-1) [0.63] |  2
320 |  6.79(-2) [0.73] |  6.79(-2) [0.73] F = 12 |  6.81(-2) [0.73] |  0.98
40 |  2.51(-1) |  2.51(-1) F = 24 |  2.55(-1) |  31
80 |  1.75(-1) [0.52] |  1.75(-1) [0.52] F = 41 |  1.77(-1) [0.52] |  16
160 |  1.13(-1) [0.62] |  1.13(-1) [0.62] F = 61 |  1.14(-1) [0.63] |  7.9
320 |  6.85(-2) [0.72] |  6.85(-2) [0.72] F = 117 |  6.88(-2) [0.73] |  3.9
Table 6: Solid body rotation. Displayed are the errors for different timesteps
along with estimated rates of convergence (in brackets). The top box is for
$m_{1}=m_{2}=99$, the middle for $m_{1}=m_{2}=199$ and the bottom for
$m_{1}=m_{2}=799$. The F indicates how many times the Merge-adapt needed to
add the BUG spaces.  Figure 4: Displayed is the rank for the Merge, Merge-
adapt and Implicit Euler method for the problem with solid body rotation
problem. Here “Fail for MA” is the indicator where the BUG space is needed for
the Merge-adapt method. For the Implicit Euler method we constantly use the
threshold $\Delta t^{2}$ when computing the rank via the truncated SVD. All
the computations are done with $m_{1}=m_{2}=799$ and $n_{T}=320.$

##  5 Conclusions and future work

This work prototypes a class of implicit adaptive low rank time-stepping
schemes. By simply merging the explicit step truncation and the BUG spaces,
the modeling error is removed and the numerical schemes achieve robust
convergence upon mesh refinement. An adaptive strategy is proposed for the
prediction of row and column spaces, which is computationally advantageous for
moderately stiff problems. The immediate future work is the generalization to
higher order in time. More importantly, the ideas can be applied to tensor
differential equations, which will be investigated in the future.

##  6 Acknowledgement

We would like to thank Steven R. White for pointing out the reference [24].

## References

  * [1] M. Bachmayr.  Low-rank tensor methods for partial differential equations.  Acta Numerica, 32:1–121, 2023.
  * [2] G. Ceruti, J. Kusch, and C. Lubich.  A rank-adaptive robust integrator for dynamical low-rank approximation.  BIT Numerical Mathematics, 62(4):1149–1174, 2022.
  * [3] G. Ceruti and C. Lubich.  An unconventional robust integrator for dynamical low-rank approximation.  BIT Numerical Mathematics, 62(1):23–44, 2022.
  * [4] A. Dektor, A. Rodgers, and D. Venturi.  Rank-adaptive tensor methods for high-dimensional nonlinear PDEs.  Journal of Scientific Computing, 88(2):36, 2021.
  * [5] R. Frank, J. Schneid, and C. W. Ueberhuber.  The concept of B-convergence.  SIAM Journal on Numerical Analysis, 18(5):753–780, 1981.
  * [6] G. H. Golub and C. F. Van Loan.  Matrix computations.  JHU press, 2013.
  * [7] L. Grasedyck, D. Kressner, and C. Tobler.  A literature survey of low-rank tensor approximation techniques.  GAMM-Mitteilungen, 36(1):53–78, 2013.
  * [8] W. Guo and J.-M. Qiu.  A low rank tensor representation of linear transport and nonlinear vlasov solutions and their associated flow maps.  Journal of Computational Physics, 458:111089, 2022.
  * [9] B. Gustafsson, H.-O. Kreiss, and J. Oliger.  Time dependent problems and difference methods, volume 2.  Wiley Online Library, 2013.
  * [10] W. Hackbusch.  Tensor spaces and numerical tensor calculus, volume 42.  Springer, 2012.
  * [11] C. D. Hauck and S. Schnake.  A predictor-corrector strategy for adaptivity in dynamical low-rank approximations.  SIAM Journal on Matrix Analysis and Applications, 44(3):971–1005, 2023.
  * [12] M. Hochbruck, M. Neher, and S. Schrammer.  Rank-adaptive dynamical low-rank integrators for first-order and second-order matrix differential equations.  BIT Numerical Mathematics, 63(1):9, 2023.
  * [13] B. N. Khoromskij.  Tensor numerical methods in scientific computing, volume 19.  Walter de Gruyter GmbH & Co KG, 2018.
  * [14] E. Kieri and B. Vandereycken.  Projection methods for dynamical low-rank approximation of high-dimensional problems.  Computational Methods in Applied Mathematics, 19(1):73–92, 2019\.
  * [15] O. Koch and C. Lubich.  Dynamical low-rank approximation.  SIAM Journal on Matrix Analysis and Applications, 29(2):434–454, 2007.
  * [16] C. Lubich and I. V. Oseledets.  A projector-splitting integrator for dynamical low-rank approximation.  BIT Numerical Mathematics, 54(1):171–188, 2014.
  * [17] C. Lubich, I. V. Oseledets, and B. Vandereycken.  Time integration of tensor trains.  SIAM Journal on Numerical Analysis, 53(2):917–941, 2015.
  * [18] C. Lubich, B. Vandereycken, and H. Walach.  Time integration of rank-constrained tucker tensors.  SIAM Journal on Numerical Analysis, 56(3):1273–1290, 2018.
  * [19] J. Nakao, J.-M. Qiu, and L. Einkemmer.  Reduced augmentation implicit low-rank (RAIL) integrators for advection-diffusion and Fokker-Planck models.  11 2023.
  * [20] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig.  Time-evolution methods for matrix-product states.  Annals of Physics, 411:167998, 2019.
  * [21] A. Rodgers and D. Venturi.  Implicit integration of nonlinear evolution equations on tensor manifolds.  Journal of Scientific Computing, 97(2):33, 2023.
  * [22] Y. Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [23] V. Simoncini.  Computational methods for linear matrix equations.  SIAM Review, 58(3):377–441, 2016.
  * [24] M. Yang and S. R. White.  Time-dependent variational principle with ancillary krylov subspace.  Phys. Rev. B, 102:094315, Sep 2020.


University at Albany, State University of New York, USA
jmcurry@albany.eduhttps://orcid.org/ 0000-0003-2504-8388Supported by NSF
CCF-1850052 and NASA 80GRC020C0016 Florida State University, Tallahassee,
<EMAIL_ADDRESS>by NSF grant DMS-1722995 Florida State
University, Tallahassee,
Floridatneedham@fsu.eduhttps://orcid.org/0000-0001-6165-3433Supported by NSF
DMS-2107808 Max Planck Institute for Mathematics in the Sciences, Leipzig,
<EMAIL_ADDRESS>Graz University of Technology,
<EMAIL_ADDRESS>by the Austrian Science Fund (FWF): W1230
Justin Curry, Washington Mio, Tom Needham, Osman Okutan, Florian
Russold[300]Mathematics of computing Algebraic topology [300]Theory of
computation Computational geometry

# Convergence of Leray Cosheaves for Decorated Mapper Graphs

Justin Curry  Washington Mio  Tom Needham  Osman Berat Okutan  Florian
Russold111Corresponding Author

###### Abstract

We introduce decorated mapper graphs as a generalization of mapper graphs
capable of capturing more topological information of a data set. A decorated
mapper graph can be viewed as a discrete approximation of the cellular Leray
cosheaf over the Reeb graph. We establish a theoretical foundation for this
construction by showing that the cellular Leray cosheaf with respect to a
sequence of covers converges to the actual Leray cosheaf as the resolution of
the covers goes to zero.

###### keywords:

Leray cosheaves, Reeb graphs, Mapper, convergence

##  1 Introduction

Reeb11footnotetext: This is an abstract of a presentation given at CG:YRF
2023. It has been made public for the benefit of the community and should be
considered a preprint rather than a formally reviewed paper. Thus, this work
is expected to appear in a conference with formal proceedings and/or in a
journal. graphs and their discrete analogs—mapper graphs—are important tools
in computational topology [3, 14, 8, 2]. They are used for data visualization
[12, 10], for comparing scalar fields via distances between Reeb graphs [7, 1,
6], and data skeletonization [9], among other things. Given a continuous map
$f\colon X\rightarrow\mathbb{R}$, the Reeb graph $R_{f}$ summarizes the zero-
dimensional connectedness of $X$. In practice, we deal with maps $f\colon
P\rightarrow\mathbb{R}$ on discrete data sets $P$, where Reeb graphs are
replaced by mapper graphs. A mapper graph is constructed by choosing a cover
$\mathcal{U}=(U_{i})_{i\in I}$ of $f(X)$ or $f(P)$, taking the components or
clustering the data points of $f^{\minus 1}(U_{i})$, collapsing each of the
obtained components to a single vertex and connecting two vertices if their
corresponding components have common points. In [6, 11] it is shown that Reeb
graphs can be viewed as cosheaves and that the cellular Reeb cosheaf w.r.t. a
cover converges to the Reeb cosheaf if the resolution of the cover goes to
zero, establishing a theoretical justification for working with finite covers.
A limitation of this approach is that Reeb/mapper graphs fail to capture
higher-dimensional topological features. To overcome these limitations we
introduce decorated mapper graphs.

Figure 1: Pulling back the cover $\mathcal{U}$ of $f(X)$ along the induced map
$\hat{f}$ on the Reeb graph and refining it into components yields a cover
$\hat{\mathcal{U}}$ of $R_{f}$. The decorated mapper graph is the nerve
complex of this cover decorated by the homology
$H_{\bullet}\big{(}\pi_{f}^{\minus 1}(-)\big{)}$ of the preimages of
intersections in $\hat{\mathcal{U}}$ under the quotient map $\pi_{f}$. Since
$f=\hat{f}\circ\pi_{f}$, this is analogous to using components of preimages
under $f$.

A decorated mapper graph collects the homology (possibly inferred from a point
sample) in all degrees $H_{\bullet}(-)\coloneqq\bigoplus_{n\geq 0}H_{n}(-)$ of
the components of $f^{\minus 1}(U_{i})$ and $f^{\minus 1}(U_{i}\cap U_{j})$ on
the corresponding vertices and edges of the mapper graph; see Figure 1. It can
be viewed as a discrete approximation of the Leray cosheaf over the Reeb
graph. In this abstract we show that the cellular Leray cosheaf w.r.t. a cover
converges to the actual Leray cosheaf for any finite 1D cover $\mathcal{U}$ as
the resolution of the cover goes to zero.

##  2 Leray Cosheaves

The Leray (pre)cosheaf [4, 5, 15] parameterizes the homology of a space $X$
when viewed along a continuous map $f\colon X\rightarrow Y$. It does this by
recording for each pair of open subsets $V\subseteq W\subseteq Y$ their
homology and the induced map $H_{\bullet}\big{(}f^{\minus
1}(V)\big{)}\xrightarrow{H_{\bullet}(f^{\minus 1}(V)\subseteq f^{\minus
1}(W))}H_{\bullet}\big{(}f^{\minus 1}(W)\big{)}$. We further assume that $X$
is compact and $H_{\bullet}\big{(}f^{\minus 1}(V)\big{)}$ is finite-
dimensional for every $V\subseteq Y$.

######  Definition 2.1 (Leray (pre)cosheaf).

We define the graded Leray (pre)cosheaf $\mathcal{L}^{f}$ of a continuous map
$f\colon X\rightarrow Y$ by the following assignments: For all $V\subseteq
W\subseteq Y$ open

$\displaystyle\mathcal{L}^{f}(V)=\bigoplus_{n\in\mathbb{N}_{0}}\mathcal{L}^{f}_{n}(V)$
$\displaystyle\coloneqq\bigoplus_{n\in\mathbb{N}_{0}}H_{n}\big{(}f^{\minus
1}(V)\big{)}$ $\displaystyle\mathcal{L}^{f}(V\subseteq
W)=\bigoplus_{n\in\mathbb{N}_{0}}\mathcal{L}^{f}_{n}(V\subseteq W)$
$\displaystyle\coloneqq\bigoplus_{n\in\mathbb{N}_{0}}H_{n}\big{(}f^{\minus
1}(V)\subseteq f^{\minus 1}(W)\big{)}\hskip 2.0pt.$

In practice, we can not access the whole Leray cosheaf. We can only get a
cellular version [4, 4.1.6] given by its values on a finite cover of $f(X)$.
An open cover $\mathcal{U}=(U_{i})_{i\in I}$ of $f(X)$ defines a simplicial
complex
$\mathcal{N}_{\mathcal{U}}\coloneqq\\{\sigma=(i_{0},\ldots,i_{k})\text{}|\text{
}\mathbf{U}_{\sigma}\coloneqq U_{i_{0}}\cap\ldots\cap
U_{i_{k}}\neq\emptyset\\}$, called the nerve complex of $\mathcal{U}$. We call
$\mathcal{U}$ a finite 1D cover if it is finite and has a one-dimensional
nerve complex and denote by $\leq$ the face-relation of
$\mathcal{N}_{\mathcal{U}}$, i.e. $\sigma\leq\tau\iff\sigma\text{ is a face of
}\tau$ .

######  Definition 2.2 (Cellular Leray cosheaf).

We define the cellular Leray cosheaf $D_{\mathcal{U}}\mathcal{L}^{f}$ w.r.t. a
cover $\mathcal{U}$ as the cosheaf on $\mathcal{N}_{\mathcal{U}}$ given by the
following assignments: For all $\sigma\leq\tau\in\mathcal{N}_{\mathcal{U}}$

$\displaystyle D_{\mathcal{U}}\mathcal{L}^{f}(\sigma)$
$\displaystyle\coloneqq\mathcal{L}^{f}(U_{\sigma})$ $\displaystyle
D_{\mathcal{U}}\mathcal{L}^{f}(\sigma\leq\tau)$
$\displaystyle\coloneqq\mathcal{L}^{f}(U_{\tau}\subseteq U_{\sigma})\hskip
2.0pt.$

If $\pi_{f}\colon X\rightarrow R_{f}$ is the quotient map from $X$ to the Reeb
graph and $\hat{\mathcal{U}}$ a cover of $R_{f}$, then we define a decorated
mapper graph as $D_{\hat{\mathcal{U}}}\mathcal{L}^{\pi_{f}}$; see Figure 1.

##  3 Convergence

It is obvious that this process of discretization can lose information. This
raises the question: How well is $\mathcal{L}^{f}$ represented by
$D_{\mathcal{U}}\mathcal{L}^{f}$? To compare $\mathcal{L}^{f}$ and
$D_{\mathcal{U}}\mathcal{L}^{f}$, we define a process of transforming
$D_{\mathcal{U}}\mathcal{L}^{f}$ into a (pre)cosheaf on $Y$. We want to
approximate the value of $\mathcal{L}^{f}$ on any open set $V\subseteq Y$
given only the information in $D_{\mathcal{U}}\mathcal{L}^{f}$. To this end,
we define the subcomplex
$K_{V}\coloneqq\\{\sigma\in\mathcal{N}_{\mathcal{U}}\text{ }|\text{
}U_{\sigma}\cap V\neq\emptyset\\}\leq\mathcal{N}_{\mathcal{U}}$, which can be
viewed as a simplicial approximation of $V\subseteq\underset{\sigma\in
K_{V}}{\bigcup}U_{\sigma}$; see Figure 2. Moreover, if $V\subseteq W$, we
obtain an inclusion $\iota\colon K_{V}\xhookrightarrow{}K_{W}$.

######  Definition 3.1 (Continuous extension).

Let $\mathcal{U}$ be a finite 1D cover of $Y$. We define the continuous
extension $C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}$ of
$D_{\mathcal{U}}\mathcal{L}^{f}$ as a (pre)cosheaf on $Y$ via the following
assignments: For all $V\subseteq W\subseteq Y$ open

$\displaystyle C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V)$
$\displaystyle\coloneqq\bigoplus_{n\in\mathbb{N}_{0}}\Big{(}H_{0}\big{(}K_{V};D_{\mathcal{U}}\mathcal{L}^{f}_{n}|_{K_{V}}\big{)}\oplus
H_{1}\big{(}K_{V};D_{\mathcal{U}}\mathcal{L}^{f}_{n\minus
1}|_{K_{V}}\big{)}\Big{)}$ $\displaystyle
C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V\subseteq W)$
$\displaystyle\coloneqq\bigoplus_{n\in\mathbb{N}_{0}}\Big{(}H_{0}\big{(}\iota\big{)}_{n}\oplus
H_{1}\big{(}\iota\big{)}_{n\minus 1}\Big{)}$

where $H_{i}\big{(}K_{V};D_{\mathcal{U}}\mathcal{L}^{f}_{n}|_{K_{V}}\big{)}$
is the degree $i$ homology [4, 6.2.2] of the cosheaf
$D_{\mathcal{U}}\mathcal{L}^{f}_{n}$ restricted to $K_{V}$ and
$H_{i}(\iota)_{n}$ is the map induced on cosheaf homology by $\iota\colon
K_{V}\xhookrightarrow{}K_{W}$ cf. [13, A.17].

As shown in [4, 5, 15], if $\mathcal{U}$ is a finite 1D cover of $f(X)$ the
continuous extension yields
$C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V)\cong\underset{n\in\mathbb{N}_{0}}{\bigoplus}H_{n}\Big{(}f^{\minus
1}\big{(}\underset{\sigma\in K_{V}}{\bigcup}U_{\sigma}\big{)}\Big{)}$
approximating $H_{\bullet}\big{(}f^{\minus 1}(V)\big{)}$; see Figure 2.

Figure 2: The figure shows $D_{\mathcal{U}}\mathcal{L}^{f}_{0}\cong
D_{\mathcal{U}}\mathcal{L}^{f}$ on $\mathcal{N}_{\mathcal{U}}$ w.r.t. two
covers with different resolution as well as the restriction
$D_{\mathcal{U}}\mathcal{L}^{f}_{0}|_{K_{V}}$ to $K_{V}$ (in orange). Homology
is taken with coefficients in a field $k$ and unlabeled arrows represent
identity maps. For the coarse cover we get
$C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V)\cong k\ncong
H_{0}\big{(}f^{\minus 1}(V)\big{)}$ but for the finer one we get
$C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V)\cong k^{2}\cong
H_{0}\big{(}f^{\minus 1}(V)\big{)}$.

The following proposition shows that
$C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}(V\subseteq W)$ also gives us
the correct induced map.

######  Proposition 3.2.

Let $\mathcal{U}$ be a finite 1D cover of $f(X)$. Then, for all $V\subseteq
W\subseteq Y$ open, the following diagram commutes and the horizontal arrows
are isomorphisms:

$\leavevmode\hbox to319.22pt{\vbox
to114.94pt{\pgfpicture\makeatletter\hbox{\hskip
159.60762pt\lower-61.62509pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{\offinterlineskip{}{}{{{}}{{}}{{}}{{}}}{{{}}}{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-159.60762pt}{-53.31956pt}\pgfsys@invoke{
}\hbox{\vbox{\halign{\pgf@matrix@init@row\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding&&\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding\cr\hfil\hskip
83.00818pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-78.70264pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{${H_{0}\big{(}K_{V};D_{\mathcal{U}}\mathcal{L}_{n}^{f}|_{K_{V}}\big{)}\oplus
H_{1}\big{(}K_{V};D_{\mathcal{U}}\mathcal{L}_{n\minus
1}^{f}|_{K_{V}}\big{)}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip
83.00818pt\hfil&\hfil\hskip 95.37718pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-47.07167pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${H_{n}\Big{(}f^{\minus
1}\big{(}\underset{\sigma\in K_{V}}{\bigcup}U_{\sigma}\big{)}\Big{)}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
51.37721pt\hfil\cr\vskip 48.00005pt\cr\hfil\hskip
85.45264pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-81.1471pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{${H_{0}\big{(}K_{W};D_{\mathcal{U}}\mathcal{L}_{n}^{f}|_{K_{W}}\big{)}\oplus
H_{1}\big{(}K_{W};D_{\mathcal{U}}\mathcal{L}_{n\minus
1}^{f}|_{K_{W}}\big{)}}$} }}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
85.45264pt\hfil&\hfil\hskip 96.15497pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-47.84946pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{${H_{n}\Big{(}f^{\minus
1}\big{(}\underset{\sigma\in K_{W}}{\bigcup}U_{\sigma}\big{)}\Big{)}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip
52.155pt\hfil\cr}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}}{{{{}}}{{}}{{}}{{}}{{}}}} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}} {}{ {}{}{}}{}{ {}{}{}} {{{{{}}{
{}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{ }{}{{ {\pgfsys@beginscope
\pgfsys@setdash{}{0.0pt}\pgfsys@roundcap\pgfsys@roundjoin{} {}{}{} {}{}{}
\pgfsys@moveto{-2.07988pt}{2.39986pt}\pgfsys@curveto{-1.69989pt}{0.95992pt}{-0.85313pt}{0.27998pt}{0.0pt}{0.0pt}\pgfsys@curveto{-0.85313pt}{-0.27998pt}{-1.69989pt}{-0.95992pt}{-2.07988pt}{-2.39986pt}\pgfsys@stroke\pgfsys@endscope}}
}{}{}{{}}{}{}{{}}\pgfsys@moveto{9.05319pt}{26.50003pt}\pgfsys@lineto{55.47545pt}{26.50003pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{55.67543pt}{26.50003pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{29.8393pt}{30.25279pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$\scriptstyle{\cong}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{-74.15498pt}{17.13367pt}\pgfsys@lineto{-74.15498pt}{-40.05989pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{-1.0}{1.0}{0.0}{-74.15498pt}{-40.25987pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{{}{}}}{{}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-114.5712pt}{-13.41309pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{H_{0}(\iota)\oplus H_{1}(\iota)}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{107.45262pt}{10.1402pt}\pgfsys@lineto{107.45262pt}{-37.05989pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{-1.0}{1.0}{0.0}{107.45262pt}{-37.25987pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{{}{}}}{{}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-17.7089pt}{-16.15977pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{$\scriptstyle{H_{n}\Big{(}f^{\minus 1}\big{(}\underset{\sigma\in
K_{V}}{\bigcup}U_{\sigma}\big{)}\subseteq f^{\minus
1}\big{(}\underset{\sigma\in K_{W}}{\bigcup}U_{\sigma}\big{)}\Big{)}}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{ {}{}{}}{}{ {}{}{}}
{{{{{}}{ {}{}}{}{}{{}{}}}}}{}{{{{{}}{
{}{}}{}{}{{}{}}}}}{{}}{}{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.39998pt}\pgfsys@invoke{
}{}{}{}{}{{}}{}{}{{}}\pgfsys@moveto{11.49765pt}{-50.81956pt}\pgfsys@lineto{54.69766pt}{-50.81956pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}}}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{54.89764pt}{-50.81956pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} {
}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{30.67264pt}{-58.07233pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$\scriptstyle{\cong}$}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\hskip 2.0pt.$

To talk about convergence, we have to define a distance on (pre)cosheaves.
Assume $Y$ is a metric space and that all the involved (pre)cosheaves are
constructible [4, 11.0.10]. For every $\epsilon>0$ define
$V^{\epsilon}\coloneqq\\{y\in Y\text{ }|\text{ }d(y,V)<\epsilon\\}$. To a
(pre)cosheaf $F$ on $Y$ we now associate
$F^{\bullet}\coloneqq\\{(F^{\epsilon})_{\epsilon\geq 0}\hskip 2.0pt,\hskip
2.0pt(F_{\epsilon}^{\delta}\colon F^{\epsilon}\rightarrow
F^{\delta})_{\epsilon\leq\delta}\\}$, a parameterized family of (pre)cosheaves
on $Y$, by setting $F^{\epsilon}(V)\coloneqq F(V^{\epsilon})$. This allows us
to define the distance of two (pre)cosheaves $F$ and $G$ as the interleaving
distance of $F^{\bullet}$ and $G^{\bullet}$, i.e. $d(F,G)\coloneqq
d_{I}(F^{\bullet},G^{\bullet})$. Denote by
$\text{res}(\mathcal{U})\coloneqq\underset{U\in\mathcal{U}}{\text{sup}}\hskip
2.0pt\underset{x,y\in U}{\text{sup}}d(x,y)$ the resolution of an open cover
$\mathcal{U}$. We are now able to show that if the resolution of the cover
$\mathcal{U}$ goes to zero, $C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f}$
converges to $\mathcal{L}^{f}$.

######  Theorem 3.3.

If $f\colon X\rightarrow Y$ is continuous and $\mathcal{U}$ is a finite 1D
cover of $f(X)$, then

$d\big{(}C_{\mathcal{U}}D_{\mathcal{U}}\mathcal{L}^{f},\mathcal{L}^{f}\big{)}\leq\text{res}(\mathcal{U})\hskip
2.0pt.$

If we only consider $\mathcal{L}^{f}_{0}$, the degree-zero part of
$\mathcal{L}^{f}$, the construction in Definition 3.1 and Theorem 3.3
specializes to an abelianization of the convergence result in [11].

##  4 Future Work

Theorem 3.3 guarantees that, if we choose a cover of resolution $\leq\delta$,
the cellular Leray cosheaf is a $\delta$-approximation of the continuous one.
This result establishes a theoretical justification for working with finite
covers. Since we have to deal with maps from finite point sets in practice, in
future work we plan to investigate under what conditions we can infer the
cellular Leray cosheaf of a map $f\colon X\rightarrow Y$ from a finite sample
of $X$; cf. [2].

## References

  * [1] Ulrich Bauer, Xiaoyin Ge, and Yusu Wang.  Measuring distance between reeb graphs.  In Proceedings of the Thirtieth Annual Symposium on Computational Geometry, SOCG’14, page 464–473, New York, NY, USA, 2014. Association for Computing Machinery.  doi:10.1145/2582112.2582169.
  * [2] Adam Brown, Omer Bobrowski, Elizabeth Munch, and Bei Wang.  Probabilistic convergence and stability of random mapper graphs.  Journal of Applied and Computational Topology, 5(1):99–140, 2021\.  doi:10.1007/s41468-020-00063-x.
  * [3] Hamish Carr, Jack Snoeyink, and Ulrike Axen.  Computing contour trees in all dimensions.  Computational Geometry, 24(2):75–94, 2003.  Special Issue on the Fourth CGC Workshop on Computational Geometry.  doi:https://doi.org/10.1016/S0925-7721(02)00093-7.
  * [4] Justin Curry.  Sheaves, Cosheaves and Applications.  PhD thesis, University of Pennsylvania, 2014.  arXiv:1303.3255.
  * [5] Justin Curry, Robert Ghrist, and Vidit Nanda.  Discrete morse theory for computing cellular sheaf cohomology.  Foundations of Computational Mathematics, 16, 2013.  doi:10.1007/s10208-015-9266-8.
  * [6] Vin de Silva, Elizabeth Munch, and Amit Patel.  Categorified reeb graphs.  Discrete & Computational Geometry, 55(4):854–906, 2016.  doi:10.1007/s00454-016-9763-9.
  * [7] Barbara Di Fabio and Claudia Landi.  The edit distance for reeb graphs of surfaces.  Discrete & Computational Geometry, 55(2):423–461, 2016.  doi:10.1007/s00454-016-9758-6.
  * [8] Herbert Edelsbrunner, John Harer, and Amit K. Patel.  Reeb spaces of piecewise linear mappings.  In Proceedings of the Twenty-Fourth Annual Symposium on Computational Geometry, SCG ’08, page 242–250, New York, NY, USA, 2008. Association for Computing Machinery.  doi:10.1145/1377676.1377720.
  * [9] Xiaoyin Ge, Issam Safa, Mikhail Belkin, and Yusu Wang.  Data skeletonization via reeb graphs.  Advances in Neural Information Processing Systems, 24, 2011.
  * [10] P. Y. Lum, G. Singh, A. Lehman, T. Ishkanov, M. Vejdemo-Johansson, M. Alagappan, J. Carlsson, and G. Carlsson.  Extracting insights from the shape of complex data using topology.  Scientific Reports, 3(1):1236, 2013.  doi:10.1038/srep01236.
  * [11] Elizabeth Munch and Bei Wang.  Convergence between Categorical Representations of Reeb Space and Mapper.  In 32nd International Symposium on Computational Geometry (SoCG 2016), volume 51 of Leibniz International Proceedings in Informatics (LIPIcs), pages 53:1–53:16. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2016.  doi:10.4230/LIPIcs.SoCG.2016.53.
  * [12] Monica Nicolau, Arnold Levine, and Gunnar Carlsson.  Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival.  Proceedings of the National Academy of Sciences of the United States of America, 108:7265–70, 04 2011.  doi:10.1073/pnas.1102826108.
  * [13] Florian Russold.  Persistent sheaf cohomology, 2022.  arXiv:2204.13446.
  * [14] Gurjeet Singh, Facundo Memoli, and Gunnar Carlsson.  Topological Methods for the Analysis of High Dimensional Data Sets and 3D Object Recognition.  In Eurographics Symposium on Point-Based Graphics. The Eurographics Association, 2007.  doi:10.2312/SPBG/SPBG07/091-100.
  * [15] Hee Rhang Yoon and Robert Ghrist.  Persistence by parts: Multiscale feature detection via distributed persistent homology, 2020.  arXiv:2001.01623.