
==== Front
J Sci Comput
J Sci Comput
Journal of Scientific Computing
0885-7474
1573-7691
Springer US New York

39148670
2629
10.1007/s10915-024-02629-8
Article
Implicit Low-Rank Riemannian Schemes for the Time Integration of Stiff Partial Differential Equations
http://orcid.org/0000-0002-8410-1372
Sutti Marco 1
Vandereycken Bart Bart.Vandereycken@unige.ch

2
1 grid.19188.39 0000 0004 0546 0241 Mathematics Division, National Center for Theoretical Sciences, National Taiwan University, Taipei, Taiwan, ROC
2 https://ror.org/01swzsf04 grid.8591.5 0000 0001 2175 2154 Section of Mathematics, University of Geneva, Geneva, Switzerland
13 8 2024
13 8 2024
2024
101 1 319 5 2023
17 6 2024
17 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
We propose two implicit numerical schemes for the low-rank time integration of stiff nonlinear partial differential equations. Our approach uses the preconditioned Riemannian trust-region method of Absil, Baker, and Gallivan, 2007. We demonstrate the efficiency of our method for solving the Allen–Cahn and the Fisher–KPP equations on the manifold of fixed-rank matrices. Our approach allows us to avoid the restriction on the time step typical of methods that use the fixed-point iteration to solve the inner nonlinear equations. Finally, we demonstrate the efficiency of the preconditioner on the same variational problems presented in Sutti and Vandereycken, 2021.

Keywords

Implicit methods
Numerical time integration
Riemannian optimization
Stiff PDEs
Manifold of fixed-rank matrices
Variational problems
Preconditioning
Trust-region method
Allen–Cahn equation
Fisher–KPP equation
Mathematics Subject Classification

65F08
65F55
65L04
65F45
65N22
65K10
58C05
http://dx.doi.org/10.13039/501100004663 Ministry of Science and Technology, Taiwan 111-2124-M-002-014 112-2124-M-002-009 Swiss National Science Foundation192129 Vandereycken Bart issue-copyright-statement© Springer Science+Business Media, LLC, part of Springer Nature 2024
==== Body
pmcIntroduction

The topic of this paper is the efficient solution of large-scale variational problems arising from the discretization of partial differential equations (PDEs), both time-independent and time-dependent. In the first part of the paper, we use the preconditioned Riemannian trust-region (RTR) method of Absil, Baker, and Gallivan [1] to solve the nonlinear equation derived from an implicit scheme for numerical time integration. All the calculations are performed on a low-rank manifold, which allows us to approximate the solution with significantly fewer degrees of freedom. In the second part of the paper, we solve variational problems derived from the discretization of elliptic PDEs. These are large-scale finite-dimensional optimization problems arising from the discretization of infinite-dimensional problems. Variational problems of this type have been considered as benchmarks in several nonlinear multilevel algorithms [26, 29, 69, 76].

A common way to speed up numerical computations is by approximating large matrices using low-rank methods. This is particularly useful for high-dimensional problems, which can be solved using low-rank matrix and tensor methods. The earliest examples are low-rank solvers for the Lyapunov equation, AX+XA⊤=C, and other matrix equations; see, e.g., [38, 40, 64]. The low-rank approximation properties for these problems are also reasonably well-understood from a theoretical point of view. Grasedyck [25, Remark 1] showed that the solution X to a Sylvester equation AX-XB+C=0 could be approximated up to a relative accuracy of ε using a rank r=O(log2(cond2(A))log2(1/ε)), under some hypotheses on A, B, and C. Typically, to obtain a low-rank approximation of the unknown solution X, an iterative method that directly constructs the low-rank approximation is used. This work uses techniques that achieve a low-rank approximation through Riemannian optimization [2, 10, 21]. To ensure critical points have a low-rank representation, the optimization problem (which may be reformulated from the original) is limited to the manifold Mr of matrices with fixed rank r. Some early references on this manifold include [27, 34, 74]. Retraction-based optimization on Mr was studied in [62, 63]. Optimization on Mr has gained a lot of momentum during the last decade, and examples of such methods are [50, 66, 74] for matrix and tensor completion, [63] for metric learning, [37, 52, 75] for matrix and tensor equations, and [57, 58] for eigenvalue problems. These optimization problems are ill-conditioned in discretized PDEs, making simple first-order methods like gradient descent excessively slow.

Riemannian Preconditioning

This work employs preconditioning techniques on Riemannian manifolds, which are similar to preconditioning techniques in the unconstrained case; see, e.g., [56]. Several authors have tackled preconditioning in the Riemannian optimization framework; the following overview is not meant to be exhaustive. In [37, 58, 75], for example, the gradient is preconditioned with the inverse of the local Hessian. Solving these Hessian equations is done by a preconditioned iterative scheme, mimicking the class of quasi or truncated Newton methods. We also refer to [72] for a recent overview of geometric methods for obtaining low-rank approximations. The work most closely related to the present paper is [75], which proposed a preconditioner for the manifold of symmetric positive semidefinite matrices of fixed rank. Boumal and Absil [11] developed a preconditioner for Riemannian optimization on the Grassmann manifold. Mishra and Sepulchre [51] investigated the connection between quadratic programming and Riemannian gradient optimization, particularly on quotient manifolds. The method of [51] proved efficient, especially in quadratic optimization with orthogonality and rank constraints. Related to this preconditioned metric approach are those of [52, 55], and more recently [14], who extend the preconditioned metric from the matrix case to the tensor case using the tensor train (TT) format for the tensor completion problem. On tensor manifolds, [37] developed a preconditioned version for Riemannian gradient descent and the Richardson method, using the Tucker and TT formats.

Trust-Region Methods

The Riemannian trust-region (RTR) method of [1] embeds an inner truncated conjugate gradient (tCG) method to solve the so-called trust-region minimization subproblem. The tCG solver naturally lends itself to preconditioning, and the preconditioner is typically a symmetric positive definite operator that approximates the inverse of the Hessian matrix. Ideally, it has to be cheap to compute. Preconditioning with the projected Euclidean Hessian was done for symmetric positive semidefinite matrices with fixed rank [75]. In contrast, we develop it here for any, i.e., typically non-symmetric, fixed-rank matrix.

We follow the steps outlined in [75] to find the preconditioner, namely: find the Euclidean Hessian, find the Riemannian Hessian operator, vectorize it to get the Hessian matrix, a linear and symmetric matrix; the inverse of the Hessian matrix should make a good candidate for a preconditioner; apply the preconditioner.

Low-Rank Approximations for Time-Dependent PDEs

Various approaches have been employed to address the low-rank approximation of time-dependent partial differential equations (PDEs). One such method is the dynamical low-rank approximation (DLRA) [34, 48], which optimally evolves a system’s low-rank approximation for common time-dependent PDEs. For example, suppose we are given a discretized dynamical system as a first-order differential equation (ODE), namely,

The DLRA idea is to replace the derivative of W with respect to time, , in the ODE with the tangent vector in the tangent space t Mr at W, TWMr, that is closer to the right-hand side G(W) . Recent developments of the DLRA include, but are not limited to, [15, 16, 32, 33, 54], and [9].

Another approach is the dynamically orthogonal Runge–Kutta of [45, 46], and its more recent developments [17, 18, 22, 23, 61, 71].

The step-truncation methods of Rodgers, Dektor, and Venturi [59, 60] form another class of methods for the low-rank approximation of time-dependent problems. In [60], they study implicit rank-adaptive algorithms involving rank truncation onto a tensor or matrix manifold where the discrete implicit integrator is obtained by fixed-point iteration. To accelerate convergence, this iteration is warm started with one time step of a conventional time-stepping scheme.

Recently, Massei et al. [49] also investigated the low-rank numerical integration of the Allen–Cahn equation. However, their approach is very different from ours since they use hierarchical low-rank matrices.

Contributions and Outline

The most significant contributions of this paper are two implicit numerical time integration schemes that can be used to solve stiff nonlinear time-dependent partial differential equations (PDEs). While also employing an implicit time-stepping scheme for the time evolution, as in [60], instead of using a fixed-point iteration method for solving the nonlinear equations derived by the time integration scheme, we use a preconditioned RTR method (named PrecRTR) on the manifold of fixed-rank matrices. Our preconditioner for the RTR subproblem on the manifold of fixed-rank matrices can be regarded as an extension of the preconditioner of [75] for the manifold of symmetric positive semidefinite matrices of fixed rank. We apply our low-rank implicit numerical time integration schemes for solving two time-dependent, stiff nonlinear PDEs: the Allen–Cahn and Fisher–KPP equations. Additionally, we consider the two variational problems already studied in [26, 29, 69, 76]. The numerical experiments demonstrate the efficiency of the preconditioned algorithm in contrast to the non-preconditioned algorithm.

The remaining part of this paper is organized as follows. Section 2 introduces the problem settings and the objective functions object of study of this work. In Sect. 3, we recall some preliminaries on the Riemannian optimization framework and the RTR method and give an overview of the geometry of the manifold of fixed-rank matrices. In Sect. 4, we recall more algorithmic details of the RTR method. Sections 5 and 6 present the core contribution of this paper: an implicit Riemannian low-rank scheme for the numerical integration of stiff nonlinear time-dependent PDEs, the Allen–Cahn equation and the Fisher–KPP equation. Other numerical experiments on the two variational problems from [69] are presented and discussed in Sect. 7. Finally, we draw some conclusions and propose a research outlook in Sect. 8. The discretization details for the Allen–Cahn and the Fisher–KPP equations are provided in Appendices A and B, respectively. More details about the derivation of the preconditioner for the RTR method on the manifold of fixed-rank matrices are given in Appendix C.

Notation

The space of n×r matrices is denoted by Rn×r. By X⊥∈Rn×(n-r) we denote an orthonormal matrix whose columns span the orthogonal complement of span(X), i.e., X⊥⊤X=0 and X⊥⊤X⊥=Ir. In the formulas throughout the paper, we typically use the Roman capital script for operators and the italic capital script for matrices. For instance, PX indicates a projection operator, while PX is the corresponding projection matrix.

The directional derivative of a function f at x in the direction of ξ is denoted by Df(x)[ξ]. With ‖·‖F, we indicate the Frobenius norm of a matrix.

Even though we did not use a multilevel algorithm in this work, we want to maintain consistency with the notation used in [69]. Consequently, we denote by ℓ the discretization level. Hence, the total number of grid points on a two-dimensional square domain is given by 22ℓ. This notation was adopted in [69] due to the multilevel nature of the Riemannian multigrid line-search (RMGLS) algorithm. In contrast, here we omit the subscripts ·h and ·H because they were due to the multilevel nature of RMGLS. We use Δ to denote the Laplacian operator, and the spatial discretization parameter is denoted by hx. The time step is represented by h.

The Problem Settings and Cost Functions

In this section, we present the optimization problems studied in this paper. The first two problems are time-dependent, stiff PDEs: the Allen–Cahn and the Fisher–KPP equations. The last two problems are the same considered in [69].

The Allen–Cahn equation in its simpler form reads1 ∂w∂t=εΔw+w-w3,

where w≡w(x,t), x∈Ω=[-π,π)2, with periodic boundary conditions, and t≥0. We reformulate it as a variational problem, which leads us to considerminwF(w):=∫Ωεh2‖∇w‖2+(1-h)2w2+h4w4-w~·wdxdy.

The second problem considered is the Fisher–KPP equation with homogeneous Neumann boundary conditions, for which we construct the cost functionF(W)=12Tr(W⊤Mm⊤MmW)-Tr((W(n-1))⊤Mp⊤MmW)+2hTr(W(n))∘2-W(n)⊤MmWRω.

We emphasize that this is the only example we do not formulate as a variational problem. We refer the reader to Sect. 6 for the details about this cost function.

Thirdly, we study the following variational problem, studied in [26, 29, 76], and called LYAP in [69, Sect. 5.1],2 minwF(w(x,y))=∫Ω12‖∇w(x,y)‖2-γ(x,y)w(x,y)dxdysuch thatw=0on∂Ω,

where ∇=(∂∂x,∂∂y), Ω=[0,1]2 and γ is the source term. The variational derivative (Euclidean gradient) of F is3 δFδw=-Δw-γ.

A critical point of (2) is thus also a solution of the elliptic PDE -Δw=γ. We refer the reader to [69, Sect. 5.1.1] or [67, Sect. 7.4.1.1] for the details about the discretization.

Finally, we consider the variational problem from [69, Sect. 5.2]:4 minwF(w)=∫Ω12‖∇w‖2+λw2(13w+12)-γwdxdysuch thatw=0on∂Ω.

For γ, we chooseγ(x,y)=ex-2y∑j=152j-1sin(jπx)sin(jπy).

which is the same right-hand side adopted in [69]. The variational derivative of F isδFδw=-Δw+λw(w+1)-γ=0.

Regardless of the specific form of the functional F, all the problems studied in this paper have the general formulationminWF(W)s.t.W∈{X∈Rn×n:rank(X)=r},

where F denotes the discretization of the functional F.

More details about each problem are provided later in the dedicated sections.

Riemannian Optimization Framework and Geometry

As anticipated above, in this paper we use the Riemannian optimization framework [2, 21]. This approach exploits the underlying geometric structure of the low-rank constrained problems, thereby allowing the constraints to be explicitly taken into account. In practice, the optimization variables in our discretized problems are constrained to a smooth manifold, and we perform the optimization on the manifold.

Specifically, in this paper, we use the RTR method of [1]. A more recent presentation of the RTR method can be found in [10]. In the next section, we introduce some fundamental geometry concepts used in Riemannian optimization, which are needed to formulate the RTR method, whose pseudocode is provided in Sect. 4.

Geometry of the Manifold of Fixed-Rank Matrices

The manifold of fixed-rank matrices is defined asMr={X∈Rm×n:rank(X)=r}.

Using the singular value decomposition (SVD), one has the equivalent characterizationMr={UΣV⊤:U∈Strm,V∈Strn,Σ=diag(σ1,σ2,…,σr)∈Rr×r,σ1≥⋯≥σr>0},

where Strm is the Stiefel manifold of m×r real matrices with orthonormal columns, and diag(σ1,σ2,…,σr) is a square matrix with σ1,σ2,…,σr on its main diagonal.

Tangent Space and Metric

The following proposition shows that Mr is a smooth manifold with a compact representation for its tangent space.

Proposition 1

([74, Prop. 2.1]) The set Mr is a smooth submanifold of dimension (m+n-r)r embedded in Rm×n. Its tangent space TXMr at X=UΣV⊤∈Mr is given by5 TXMr=UU⊥Rr×rRr×(n-r)R(m-r)×r0(m-r)×(n-r)VV⊥⊤.

In addition, every tangent vector ξ∈TXMr can be written as6 ξ=UMV⊤+UpV⊤+UVp⊤,

with M∈Rr×r, Up∈Rm×r, Vp∈Rn×r such that Up⊤U=Vp⊤V=0.

The orthogonality conditions Up⊤U=Vp⊤V=0 are also known as gauging conditions [72, §9.2.3]. Since Mr⊂Rm×n, we represent tangent vectors in (5) and (6) as matrices of the same dimensions.

The Riemannian metric is the restriction of the Euclidean metric on Rm×n to the submanifold Mr, i.e.,gX(ξ,η)=⟨ξ,η⟩=Tr(ξ⊤η),withX∈Mrandξ,η∈TXMr.

Projectors

Defining PU=UU⊤ and PU⊥=I-PU for any U∈Strm, where Strm is the Stiefel manifold of m-by-r orthonormal matrices, the orthogonal projection onto the tangent space at X is [74, (2.5)]7 PX:Rm×n→TXMr,Z↦PUZPV+PU⊥ZPV+PUZPV⊥.

Since this projector is a linear operator, we can represent it as a matrix. The projection matrix PX∈Rn2×n2 representing the operator PX can be written asPX:=PV⊗PU+PV⊗PU⊥+PV⊥⊗PU.

Riemannian Gradient

The Riemannian gradient of a smooth function f:Mr→R at X∈Mr is defined as the unique tangent vector gradf(X) in TXMr such that∀ξ∈TXMr,⟨gradf(X),ξ⟩=Df(X)[ξ],

where Df denotes the directional derivatives of f. More concretely, for embedded submanifolds, the Riemannian gradient is given by the orthogonal projection onto the tangent space of the Euclidean gradient of f seen as a function on the embedding space Rm×n; see, e.g., [2, (3.37)]. Then, denoting ∇f(X) the Euclidean gradient of f at X, the Riemannian gradient is given by8 gradf(X)=PX(∇f(X)).

Riemannian Hessian

The Riemannian Hessian is defined by (see, e.g., [2, def. 5.5.1], [10, def. 5.14])Hessf(x)[ξx]=∇ξxgradf(x),

where ∇ξx is the Levi-Civita connection. If M is a Riemannian submanifold of the Euclidean space Rn, as it is the case for the manifold of fixed-rank matrices, it follows that [10, cor. 5.16]∀ξ∈TxM,Hessf(x)[ξ]=PX(Dgradf(x)[ξ]).

In practice, this is what we use in the calculations.

Retraction

To map the updates in the tangent space onto the manifold, we make use of so-called retractions. A retraction RX is a smooth map from the tangent space to the manifold, RX:TXMr→Mr, used to map tangent vectors to points on the manifold. It is, essentially, any smooth first-order approximation of the exponential map of the manifold; see, e.g., [3]. To establish convergence of the Riemannian algorithms, it is sufficient for the retraction to be defined only locally. An excellent survey on low-rank retractions is given in [4]. In our setting, we have chosen the metric projection, which is provided by a truncated SVD.

The RTR Method

As we anticipated above, to solve the implicit equation resulting from the time-integration scheme, we employ the RTR method of [1]. For reference, we provide the pseudocode for RTR in Algorithm 1. Step 4 in Algorithm 1 uses the truncated conjugate gradient (tCG) of [65, 70]. This method lends itself very well to being preconditioned.

Algorithm 1 RTR method of [1]

Riemannian Gradient and Riemannian Hessian

In general, in the case of Riemannian submanifolds, the full Riemannian Hessian of an objective function f at x∈M is given by the projected Euclidean Hessian plus the curvature part9 Hessf(x)[ξ]=Px∇2f(x)Px+Pxp(``curvature terms'')Pxp.

This suggests using Px∇2f(x)Px as a preconditioner in the RTR scheme; see [75, §6.1] for further details.

For the LYAP problem, the Riemannian gradient is given bygradF(X)=PX(hx2(AX+XA-Γ)).

The directional derivative of the gradient, i.e., the Euclidean Hessian applied to ξ∈TXMr, isHessF(X)[ξ]=DgradF(X)[ξ]=hx2(Aξ+ξA).

The orthogonal projection of the Euclidean Hessian followed by vectorization yieldsvec(PX(DgradF(x)[ξ]))=hx2PXvec(Aξ+ξA)=hx2PX(A⊗I+I⊗A)vec(ξ)=hx2PX(A⊗I+I⊗A)PXvec(ξ),

where the second PX is inserted for symmetrization. From here we can read the symmetric n2-by-n2 matrix10 HX=hx2PX(A⊗I+I⊗A)PX.

The inverse of this matrix (10) should be a good candidate for a preconditioner. In Appendix C, we present the derivation of the preconditioner for the tCG subsolver on the manifold of fixed-rank matrices.

In general, the preconditioner from above cannot be efficiently inverted because of the coupling with the nonlinear terms. Nonetheless, numerical experiments in Sect. 7 show that it remains an efficient preconditioner even for problems with a (mild) nonlinearity.

The Allen–Cahn Equation

The Allen–Cahn equation is a reaction-diffusion equation originally studied for modeling the phase separation process in multi-component alloy systems [5, 6]. It later turned out that the Allen–Cahn equation has a much wider range of applications. Recently, [78] provided a good review. Applications include mean curvature flows [41], two-phase incompressible fluids [77], complex dynamics of dendritic growth [44], image inpainting [20, 47], and image segmentation [7, 43].

The Allen–Cahn equation in its simplest form is given by (1). It is a stiff PDE with a low-order polynomial nonlinearity and a diffusion term εΔw. As in [60], we set ε=0.1, and we solve (1) on a square domain [-π,π)2 with periodic boundary conditions, and we also use the same initial condition as in [60, (77)–(78)], namely,11 w0(x,y)=u(x,y)-u(x,2y)+u(3x+π,3y+π)-2u(4x,4y)+2u(5x,5y).

whereu(x,y)=e-tan2(x)+e-tan2(y)sin(x)sin(y)1+e|csc(-x/2)|+e|csc(-y/2)|.

We emphasize that with this choice, the matrix W0 which discretizes the initial condition (11) has no low-rank structure and will be treated as a dense matrix. Nonetheless, thanks to the Laplacian’s smoothing effect as time evolution progresses, the solution W can be well approximated by low-rank matrices [49, §4.2]. In particular, for large simulation times, the solution converges to either -1 or 1 in most of the domain, giving rise to four flat regions that can be well approximated by low rank; see panels (e) and (f) of Fig. 1.

Spatial Discretization

We discretize (1) in space on a uniform grid with 256×256 points. In particular, we use the central finite differences to discretize the Laplacian with periodic boundary conditions. This results in the matrix ODE12

where W:[0,T]→R256×256 is a matrix that depends on t, ∘3 denotes the elementwise power of a matrix (so-called Hadamard power, defined by W∘α=[wijα]), and A is the second-order periodic finite difference differentiation matrix13 A=1hx2-2111-21⋱⋱⋱1-2111-2.

This matrix ODE is an initial value problem (IVP) in the form of [72, (48)]14

where G:=εAW+WA+W-W∘3 is the right-hand side of (12).

Reference Solution

To get a reference solution Wref, we solve the (full-rank) IVP problem (14) with the classical explicit fourth-order Runge–Kutta method (ERK4), with a time step h=10-4. Figure 1 illustrates the time evolution of the solution to the Allen–Cahn equation at six different simulation times. It is apparent that the solution evolves from an initial condition with many peaks and valleys to a solution with four flat regions occupying most of the domain.Fig. 1 Time evolution of the solution w to the Allen–Cahn equation, computed with ERK4, h=10-4

As a preliminary study, we monitor the discrete L2-norm of the right-hand side of (1) for this reference solution and the numerical rank history of Wref. From panel (a) of Fig. 2, it appears that after t≈13, ‖∂w/∂t‖L2(Ω)≤10-3, which means that the solution w enters a stationary phase; see also last two panels of Fig. 1. Panel (b) of Fig. 2 plots the numerical rank of Wref versus time, with relative singular value tolerance of 10-10. The numerical rank exhibits a rapid decay during the first ≈2 seconds, then varies between 13 and 17 during the rest of the simulation. The rank decreases as the diffusion term comes to dominate the system.Fig. 2 Preliminary numerical study of Wref

Low-Rank Implicit Time-Stepping Scheme

As mentioned above, we employ the implicit Euler method for the time integration of (12), which gives15 Wk+1=Wk+h·G(Wk+1),

and, additionally, we want Wk to be of low rank. This is achieved by using our preconditioned RTR method (PrecRTR) on the manifold of fixed-rank matrices to solve for Wk+1 the nonlinear equation (15).

Since our strategy is optimization, and since we wish to maintain some coherence with the LYAP and NPDE problems presented in Sect. 2 and [69], we build a variational problem whose first-order optimality condition will be exactly (15). This leads us to consider the problem16 minwF(w):=∫Ωεh2‖∇w‖2+(1-h)2w2+h4w4-w~·wdxdy.

It is interesting to note that this cost function is very similar to the NPDE functional [69, (5.11)]. Here, w~ is the solution at the previous time step and plays a similar role as γ in the NPDE functional (it is constant w.r.t. w). The only additional term w.r.t. the NPDE problem is the term h/4w4. Moreover, in contrast to LYAP and NPDE, we need to solve this optimization problem many times, i.e., at every time step, to describe the time evolution of w.

We aim to obtain good low-rank approximations on the whole interval [0, T] . However, it is clear from our preliminary study (panel (b) of Fig. 2) that at the beginning of the time evolution, the numerical solution is not really low rank due to the initial condition chosen. For this reason, in our numerical experiments, we consider the dense matrix until t0=0.5, and only then do we start our rank-adaptive method. Indeed, according to the rank history of the reference solution (panel (b) of Fig. 2), at t=0.5 the numerical rank has already dropped to 23. We deemed this to be a good starting time for our low-rank algorithm. Our procedure is summarized in Algorithm 2.

Algorithm 2 Low-rank Riemannian implicit Euler for the Allen–Cahn equation

The discretizations of the objective function F(w) and its gradient are detailed in Appendix A.

Numerical Experiments

The algorithm was implemented in MATLAB and is publicly available at https://github.com/MarcoSutti/PrecRTR. The RTR method of [1] was executed using solvers from the Manopt package [12] with the Riemannian embedded submanifold geometry from [74]. We conducted our experiments on a desktop machine with Ubuntu 22.04.1 LTS and MATLAB R2022a installed, with Intel Core i7-8700 CPU, 16GB RAM, and Mesa Intel UHD Graphics 630.

For the time integration, we use the time steps h={0.05,0.1,0.2,0.5,1}, and we monitor the error ‖w-wref‖L2(Ω). Figure 3 reports on the results. Panel (a) shows the time evolution of the error ‖w-wref‖L2(Ω), while panel (b) shows that the error decays linearly in h, as expected.Fig. 3 Panel (a): error versus time for the preconditioned low-rank evolution of the Allen–Cahn equation (1) with initial condition (11). Panel (b): error at T=15 versus time step h

Figure 4 reports on the numerical rank history for the preconditioned low-rank evolution of the Allen–Cahn equation, with h=0.05.Fig. 4 Rank versus time for the preconditioned low-rank evolution of the Allen–Cahn equation (1) with initial condition (11), with h=0.05

Discussion/Comparison with Other Solvers

From the results reported on Fig. 3, it is evident that even with very large time steps, we can still obtain relatively good low-rank approximations of the solution, especially at the final time T=15. For example, compare with Fig. 4 in [60], where the biggest time step is h=0.01 — i.e., one hundred times smaller than our largest time step. Moreover, factorized formats are not mentioned in [60]. In contrast, we always work with the factors to reduce computational costs.

In [60], the authors study implicit rank-adaptive algorithms based on performing one time step with a conventional time-stepping scheme, followed by an implicit fixed-point iteration step involving a rank truncation operation onto a tensor or matrix manifold. Here, we also employ an implicit time-stepping scheme for the time evolution. Still, instead of using a fixed-point iteration method for solving the nonlinear equations, we use our preconditioned RTR (PrecRTR) on the manifold of fixed-rank matrices. This way, we obtain a preconditioned dynamical low-rank approximation of the Allen–Cahn equation.

In general, implicit methods are much more effective for stiff problems, but they are also more expensive than their explicit counterparts since solutions of nonlinear systems replace function evaluations. Nonetheless, the additional computational overhead of the implicit method is compensated by the fact that we can afford a larger time step, as demonstrated by Fig. 3. Moreover, the cost of solving the inner nonlinear equations remains moderately low thanks to our preconditioner.

On the one hand, another known issue with explicit methods is that the time step needs to be in proportion to the smallest singular value of the solution [72, §9.5.3]. On the other hand, when using fixed point iterations, one still obtains a condition on the time step size, which depends on the Lipschitz constant of the right-hand side term, namelyh<1‖A‖∞‖G‖∞,

where A is the matrix of the coefficients defining the stages (the Butcher tableau); see, e.g., [36, (3.13)] and [19]. This condition appears to be a restriction on the time step, not better than the restrictions for explicit methods to be stable. This shows that (quoting from [36, §3.2]) “fixed point iterations are unsuitable for solving the nonlinear system defining the stages. For solving the nonlinear system, other methods like the Newton method should be used”. A similar condition also holds for the method of Rodgers and Venturi, see [60, (31)]:h<1LG.

Their paper states: “Equation (31) can be seen as a stability condition restricting the maximum allowable time step h for the implicit Euler method with fixed point iterations.” This makes a case for using the Newton method instead of fixed point iteration to find a solution to the nonlinear equation.

As we observed from the MATLAB profiler,1 as ℓ increases, the calculation of the preconditioner becomes dominant in the running time. We are in the best possible situation since the preconditioner dominates the cost.

Finally, we emphasize that low-rank Lyapunov solvers (see [64] for a review) cannot be used to solve this kind of problem due to a nonlinear term in the Hessian, and PrecRTR proves much more effective than the RMGLS method of [69]. However, the latter may remain useful in all those problems for which an effective preconditioner is unavailable.

Of course, our method proves efficient when the rank is low and the time step is not too small. Otherwise, if these conditions are not met, there is no advantage over using full-rank matrices.

The Fisher–KPP Equation

The Fisher–KPP equation is a nonlinear reaction-diffusion PDE, which in its simplest form reads [53, (13.4)]17 ∂w∂t=∂2w∂x2+r(ω)w(1-w),

where w≡w(x,t;ω), r(ω) is a species’s reaction rate or growth rate. It is called “stochastic”2 Fisher–KPP equation in the recent work of [18].

It was originally studied around the same time in 1937 in two independent, pioneering works. Fisher [24] studied a deterministic version of a stochastic model for the spread of a favored gene in a population in a one-dimensional habitat, with a “logistic” reaction term. Kolmogorov, Petrowsky, and Piskunov provided a rigorous study of the two-dimensional equation and obtained some fundamental analytical results, with a general reaction term. We refer the reader to [35] for an English translation of their original work.

The Fisher–KPP equation can be used to model several phenomena in physics, chemistry, and biology. For instance, it can be used to describe biological population or chemical reaction dynamics with diffusion. It has also been used in the theory of combustion to study flame propagation and nuclear reactors; see [53, §13.2] for a comprehensive review.

Boundary and Initial Conditions

Here, we adopt the same boundary and initial conditions as in [18]. The reaction rate is modeled as a random variable that follows a uniform law r∼U1/4,1/2. We consider the spatial domain x∈[0,40] and the time domain t∈[0,10]. We impose homogeneous Neumann boundary conditions, i.e.,∀t∈[0,10],∂w∂x(0,t)=0,∂w∂x(40,t)=0.

These boundary conditions represent the physical condition of zero diffusive fluxes at the two boundaries. The initial condition is “stochastic”, of the formw(x,0;ω)=a(ω)e-b(ω)x2,

where a∼U1/5,2/5 and b∼U1/10,11/10. The random variables a, b, and r are all independent, and we consider Nr=1000 realizations.

Reference Solution with the IMEX-CNLF Method

To obtain a reference solution, we use the implicit-explicit Crank–Nicolson leapfrog scheme (IMEX-CNLF) for time integration [30, Example IV.4.3]. This scheme treats the linear diffusion term with Crank–Nicolson, an implicit method. In contrast, the nonlinear reaction term is treated explicitly with leapfrog, a numerical scheme based on the implicit midpoint method.

For the space discretization, we consider 1000 grid points in x, while for the time discretization, we use 1601 points in time, so that the time step is h=10/(1601-1)=0.00625.

Let w(i) denote the spatial discretization of the ith realization. At a given time t, each realization is stored as a column of our solution matrix, i.e.,

Moreover, let Rω be a diagonal matrix whose diagonal entries are the r(ω)(i) coefficients for every realization indexed by i, i=1,2,…,Nr. Indeed,

The IMEX-CNLF scheme applied to (17) gives the algebraic equation18 (I-hA)W(n+1)=(I+hA)W(n-1)+2hW(n)Rω-2h(W(n))∘2Rω,

where A is the matrix that discretizes the Laplacian with a second-order centered finite difference stencil and homogeneous Neumann boundary conditions, i.e.,19 A=1hx2-221-21⋱⋱⋱1-212-2.

For ease of notation, we call Mm=I-hA and Mp=I+hA, so that (18) becomes20 MmW(n+1)=MpW(n-1)+2hW(n)Rω-2h(W(n))∘2Rω,

Panels (a) and (b) of Fig. 5 show the 1000 realizations at t=0 and at t=10, respectively. Panel (c) reports on the numerical rank history. To compute the numerical rank, we use MATLAB’s default tolerance, about 10-11.Fig. 5 Fisher–KPP reference solution computed with an IMEX-CNLF scheme. Panel (a): all the 1000 realizations at t=0. Panel (b): all the 1000 realizations at t=10. Panel (c): numerical rank history

Low-Rank Crank–Nicolson Leapfrog (LR-CNLF) Scheme

To obtain a low-rank solver for the Fisher–KPP PDE, we proceed similarly as for the Allen–Cahn equation low-rank solution. We build a cost function F(W) , so that its minimization gives the solution to (20), i.e.,minWF(W):=12MmW-MpW(n-1)+2h(W(n))∘2-W(n)RωF2.

Developing and keeping only the terms that depend on W, we get the cost function21 F(W)=12Tr(W⊤Mm⊤MmW)-Tr((W(n-1))⊤Mp⊤MmW)+2hTr(W(n))∘2-W(n)⊤MmWRω.

Appendix B provides further details about the low-rank formats of (21), its gradient, and Hessian.

Numerical Experiments

We monitor the following quantities:the numerical rank of the solution WLR-CNLF given by the low-rank solver;

the discrete L2-norm of the error ‖wLR-CNLF-wCNLF‖L2(Ω)=hx·‖WLR-CNLF-WCNLF‖F.

As was done in the previous section for the reference solution, here we also consider 1000 realizations. We apply our technique with rank adaption, with tolerance for rank truncation of 10-8. The inner PrecRTR is halted once the gradient norm is less than 10-8. Figure 6 reports on the numerical experiments. Panel (a) reports on the numerical rank history of WLR-CNLF; the numerical rank of the reference solution WCNLF from Sect. 6.2 is also plotted as reference. Panel (b) reports the history of the discrete L2-norm of the error versus time for several h. It is clear that the low-rank approximation improves for smaller values of h, but this improvement is counterbalanced by an error that accumulates as the simulation progresses.Fig. 6 Panel (a): rank history for the LR-CNLF method compared to the reference solution, for h=0.00625. Panel (b): discrete L2-norm of the error versus time, for several h

Numerical Experiments for LYAP and NPDE

This section focuses on the numerical properties of PrecRTR, our preconditioned RTR on the manifold of fixed-rank matrices on the variational problems from [69], recalled in Sect. 2. These are large-scale finite-dimensional optimization problems arising from the discretization of infinite-dimensional problems. These problems have been used as benchmarks in several nonlinear multilevel algorithms [26, 29, 76]. For further information on the theoretical aspects of variational problems, we recommend consulting [13, 42].

We consider two scenarios: in the first one, we let Manopt automatically take care of the trust-region radius Δ¯, while we fix Δ¯=0.5 in the second one. The tolerance on the norm of the gradient in the trust-region method is set to 10-12, and we set the maximum number of outer iterations nmax outer=300.

Tables

Tables 1, 2, 3 report the numerical results of the LYAP problem, while Tables 4, 5, 6 correspond to the NPDE problem. In all the tables, ℓ indicates the level of discretization, i.e., the total number of grid points on a two-dimensional square domain equals 22ℓ. The quantities ‖ξ(end)‖F and r(W(end)) are the Frobenius norm of the gradient and the residual, respectively, both evaluated at the final iteration of the simulation, W(end).

Tables 1 and 4 report the monitored quantities ‖ξ(end)‖F and r(W(end)) for our PrecRTR, for the LYAP and NPDE problems, respectively. CPU times are in seconds and were obtained as an average over 10 runs.

When the maximum number of outer PrecRTR iterations nmax outer=300 is reached, we indicate this in bold text. We also set a limit on the cumulative number of inner iterations ∑ninner: the inner tCG solver is stopped when ∑ninner first exceeds 30000. This is also highlighted using the bold text in the following tables.Table 1 Preconditioned RTR for the LYAP problem

ℓ	Size	Rank 5	Rank 10	
Time	‖ξ(end)‖F	r(W(end))	Time	‖ξ(end)‖F	r(W(end))	
10	1 048 576	0.21	1.0150×10-13	9.7480×10-8	0.85	6.2481×10-14	4.2204×10-11	
11	4 194 304	0.49	2.9645×10-14	4.8741×10-8	1.53	5.7690×10-13	2.0374×10-11	
12	16 777 216	1.01	3.8413×10-14	2.4371×10-8	2.93	1.0921×10-13	1.0478×10-11	
13	67 108 864	1.56	7.3017×10-14	1.2185×10-8	5.74	1.3556×10-13	5.2396×10-12	
14	268 435 456	3.80	1.5082×10-13	6.0927×10-9	10.87	9.3753×10-14	2.6045×10-12	
15	1 073 741 824	7.48	2.7525×10-13	3.0464×10-9	25.02	2.4835×10-13	1.3177×10-12	

Table 2 Effect of preconditioning: dependence on ℓ for LYAP

Prec	ℓ	Rank 5	Rank 10	
10	11	12	13	14	15	10	11	12	13	14	15	
No	nouter	51	54	61	59	162	92	300	103	61	63	62	59	
∑ninner	4 561	9 431	21 066	36 556	30 069	30 096	27 867	30 025	33 818	45 760	44 467	38 392	
maxninner	1 801	3 191	7 055	9 404	1 194	1 851	2 974	3 385	8 894	24 367	24 537	25 013	
Yes	nouter	41	45	50	52	56	60	44	64	62	53	56	56	
∑ninner	44	45	50	52	56	60	69	104	82	60	69	56	
maxninner	4	1	1	1	1	1	9	9	8	8	8	1	

Tables 2 and 5 report on the effect of preconditioning as the problem size ℓ increases for LYAP and NPDE, respectively. The reductions in the number of iterations of the inner tCG between the non-preconditioned (rows 3–5) and the preconditioned (last three rows) versions are impressive. Moreover, for the preconditioned method (last three rows in the tables), both tables demonstrate that nouter and ∑ninner depend (quite mildly) on the problem size ℓ, while maxninner is basically constant.

For NPDE, in both the non-preconditioned and preconditioned methods, the numbers of iterations are typically higher than those for the LYAP problem, which is plausibly due to the nonlinearity of the problem.

Finally, Tables 3 and 6 report the results for varying rank and fixed problem size ℓ=12. The stopping criteria are the same as above. It is remarkable that, for PrecRTR for the LYAP problem, all three monitored quantities basically do not depend on the rank. For NPDE, there is some more, but still moderate, dependence on the rank.Table 3 Effect of preconditioning: dependence on the rank with fixed size ℓ=12, for LYAP

Prec	Iterations	Rank	
1	2	5	10	15	20	
No	nouter	53	53	61	61	300	62	
∑ninner	17 650	18 775	21 066	33 818	12 816	33 292	
maxninner	6 276	7 225	7 055	8 894	3 794	6 928	
Yes	nouter	51	51	50	49	49	48	
∑ninner	51	51	50	49	49	48	
maxninner	1	1	1	1	1	1	

Table 4 Preconditioned RTR for the NPDE problem

ℓ	Size	Rank 5	Rank 10	
Time	‖ξ(end)‖F	r(W(end))	Time	‖ξ(end)‖F	r(W(end))	
Rank 5	
10	1 048 576	0.45	2.0719×10-14	1.5614×10-5	1.17	1.7303×10-14	1.8660×10-7	
11	4 194 304	0.89	2.7106×10-14	7.8072×10-6	2.10	6.0181×10-14	9.3301×10-8	
12	16 777 216	1.65	5.2974×10-14	3.9036×10-6	4.73	5.9537×10-14	4.6650×10-8	
13	67 108 864	2.84	1.2492×10-13	1.9518×10-6	8.91	1.1536×10-13	2.3325×10-8	
14	268 435 456	5.89	2.4349×10-13	9.7591×10-7	19.67	2.6992×10-13	1.1663×10-8	
15	1 073 741 824	12.96	6.4490×10-13	4.8796×10-7	45.71	5.8336×10-13	5.8313×10-9	

Table 5 Effect of preconditioning: dependence on ℓ for NPDE

Prec	ℓ	Rank 5	Rank 10	
10	11	12	13	14	15	10	11	12	13	14	15	
No	nouter	53	57	61	79	68	68	63	87	76	68	62	65	
∑ninner	4 603	9 505	13 817	41 144	47 186	38 079	6 610	38 858	30 567	31 028	39 803	39 337	
maxninner	2 022	3 595	7 735	14 195	28 410	32 433	1 487	11 550	6 035	10 598	22 468	30 118	
Yes	nouter	50	56	61	63	65	66	53	58	63	69	69	71	
∑ninner	57	64	69	72	74	75	78	84	90	98	97	100	
maxninner	6	7	7	7	7	7	10	10	11	11	10	10	

Table 6 Effect of preconditioning: dependence on the rank for NPDE with fixed level ℓ=12

Prec	Iterations	Rank	
1	2	5	10	15	20	
No	nouter	59	57	61	76	62	60	
∑ninner	9 183	16 044	13 826	30 567	61 339	31 192	
maxninner	3 569	4 642	7 744	6 035	31 627	8 540	
Yes	nouter	59	61	61	63	60	61	
∑ninner	78	90	69	90	90	104	
maxninner	11	10	7	11	11	13	

Conclusions and Outlook

In this paper, we have shown how to combine an efficient preconditioner with optimization on low-rank manifolds. Unlike classical Lyapunov solvers, our optimization strategy can treat nonlinearities. Moreover, compared to iterative methods that perform rank-truncation at every step, our approach allows for much larger time steps as it does not need to satisfy a fixed-point Lipschitz restriction. We illustrated this technique by applying it to two time-dependent nonlinear PDEs — the Allen–Cahn and the Fisher–KPP equations. In addition, the numerical experiments for two time-independent variational problems demonstrate the efficiency in computing good low-rank approximations with a number of tCG iterations in the trust region subsolver which is almost independent of the problem size.

Future research may focus on higher-order methods, such as more accurate implicit methods. Additionally, we may explore higher-dimensional problems, problems in biology, and stochastic PDEs.

Low-Rank Formats for the Allen–Cahn Equation (ACE)

Objective Functional

Discretizing (16) similarly as in [69, §5.2.1], we obtain22 F=hx2∑i,j=02ℓ-1εh2(∂wxij2+∂wyij2)+1-h2wij2+h4wij4-w~ijwij.

To obtain the factored format of the discretized objective functional, we consider the factorizations W=UΣV⊤, and W~=U~Σ~V~⊤.

The first term and the fourth term in (22) have the same factorized form as those seen in [69, §5.2.1]. The only slight change is due to the periodic boundary conditions adopted here. As a consequence, the matrix L that discretizes the first-order derivatives.3 with periodic boundary conditions becomesL=1hx-11-11⋱⋱1-11.

Note the presence of the unitary coefficient in the lower-left corner. The reader can easily verify that A=L⊤L, where A is the matrix (13) obtained by discretizing the Laplacian with central finite differences and periodic boundary conditions. We recall from [69] that, given this matrix, the first-order derivatives of W can be computed as∂Wx=LWand∂Wy=WL⊤.

For the second term in (22), we have∑i,j1-h2wij2=1-h2Tr(W⊤W)=1-h2‖Σ‖F2.

For the third term, it is easier to consider the full-rank format23 h4∑i,jwij4.

We call G~=U~Σ~ and G=UΣ. Finally, the discretized objective functional in factorized matrix form isF=hx2εh2‖(LU)Σ‖F2+‖(LV)Σ‖F2+1-h2‖Σ‖F2+h4∑i,jwij4-Tr((G~⊤G)(V⊤V~)).

Table 7 summarizes the asymptotic complexities for the ACE cost function. In this and the following similar tables, we indicate the sizes of the matrices in the order in which they appear in the product. If all the matrices in a term are the same size, we only indicate that once. Matrices without any specific structure are stored as dense, unless otherwise specified.Table 7 Asymptotic complexities for ACE cost function

Product	Factor sizes	Notes on structure and storage	Cost	
‖(LU)Σ‖F2	n×n, n×r, r×r	L sparse banded, Σ sparse diagonal	O(nr)	
‖(LV)Σ‖F2	n×n, n×r, r×r	L sparse banded, Σ sparse diagonal	O(nr)	
‖Σ‖F2	r×r	Σ sparse diagonal	O(r)	
∑i,jwij4	n×n		O(n2)	
Tr((G~⊤G)(V⊤V~))	O(nr2+r3)	
G~⊤G	n×r, n×r		O(nr2)	
V⊤V~	n×r, n×r		O(nr2)	
(G~⊤G)(V⊤V~)	r×r, r×r		O(r3)	
Tr(·)	r×r		O(r)	

Gradient

The gradient of F (16) is the variational derivativeδFδw=-εhΔw+(1-h)w+hw3-w~.

The discretized Euclidean gradient in matrix form is given byG=hx2-εh(AW+WA)+(1-h)W+hW∘3-W~,

with A as in (13).

For the term W∘2=W⊙W, we perform the element-wise multiplication in factorized form as explained in [39, §7] and store the result in the format U∘2Σ∘2V∘2⊤, i.e.,W∘2=W⊙W=(U∗⊤U)(Σ⊗Σ)(V∗⊤V)⊤=U∘2Σ∘2V∘2⊤,

where ∗⊤ denotes a transposed variant of the Khatri–Rao product [31]. Then for W∘3=W⊙W⊙W we consider the factorized formatW∘3=W∘2⊙W=(U∘2∗⊤U)(Σ∘2⊗Σ)(V∘2∗⊤V)⊤=U∘3Σ∘3V∘3⊤,

Substituting the formats W=UΣV⊤, W∘3=U∘3Σ∘3V∘3⊤, and W~=U~Σ~V~⊤, we get the factorized form of the Euclidean gradient G=UGΣGVG⊤, whereUG=-εhA+(1-h)IUUU∘3U~,ΣG=hx2blkdiagΣ,(-εh)Σ,hΣ∘3,-Σ~,

andVG=VAVV∘3V~.

The gradient G is just an augmented matrix, analogously to the discretized gradient in factored format for the NPDE problem; see [69, §5.2.2]. The operations needed to form G are summarized in Table 8.Table 8 Asymptotic complexities for ACE gradient

Product	Factor sizes	Notes on structure and storage	Cost	
-εhA+(1-h)IU	n×n, n×r	A, I sparse banded	O(nr)	
AV	n×n, n×r	A sparse banded	O(nr)	
U∘2=U∗⊤U	n×r		O(nr2)	
Σ∘2=Σ⊗Σ	r×r	Σ sparse diagonal	O(r2)	
V∘2=V∗⊤V	n×r		O(nr2)	
U∘3=U∘2∗⊤U	n×r2, n×r		O(nr3)	
Σ∘3=Σ∘2⊗Σ	r2×r2, r×r	Σ∘2, Σ sparse diagonal	O(r3)	
V∘3=V∘2∗⊤V	n×r2, n×r		O(nr3)	

Hessian

The discretized Euclidean Hessian is (compare [68, §7.4.2.3])HW[η]=hx2-εh(Aη+ηA)+3hW∘2⊙η+(1-h)η.

The factored form of the discretized Euclidean Hessian is HW[η]=UHW[η]SHW[η]VHW[η]⊤, withUHW[η]=-εhA+(1-h)IUηUηU⊙,SHW[η]=hx2blkdiagSη,(-εh)Sη,3hΣ⊙,VHW[η]=VηAVηV⊙.

where η=UηSηVη⊤ is a tangent vector in TWMr, and W∘2⊙η=U⊙Σ⊙V⊙⊤ (Table 9).Table 9 Asymptotic complexities for ACE Hessian

Product	Factor sizes	Notes on structure and storage	Cost	
-εhA+(1-h)IUη	n×n, n×r	A sparse banded, I sparse diagonal	O(nr)	
AVη	n×n, n×r	A sparse banded	O(nr)	
U∘2=U∗⊤U	n×r		O(nr2)	
Σ∘2=Σ⊗Σ	r×r	Σ sparse diagonal	O(r2)	
V∘2=V∗⊤V	n×r		O(nr2)	
U⊙=U∘2∗⊤Uη	n×r2, n×r		O(nr3)	
Σ⊙=Σ∘2⊗Ση	r2×r2, r×r	Σ∘2, Ση sparse diagonal	O(r3)	
V⊙=V∘2∗⊤Vη	n×r2, n×r		O(nr3)	

Remark 1

We have the following relationships between the cost function, the gradient, and their discretized counterparts: i.e., we can first discretize F to obtain F, and then the Euclidean gradient of F is G. This is equivalent to computing the variational derivative δFδw first, and then discretizing it to obtain G.

Low-Rank Formats for the Fisher–KPP Equation (FKPPE)

Cost Function

The FKPPE cost function (21) involves the quantities W, W(n), and W(n-1). In low-rank matrix format, these are factorized asW=UΣV⊤,W(n)=U(n)Σ(n)(V(n))⊤,W(n-1)=U(n-1)Σ(n-1)(V(n-1))⊤,

where all the U and V factors have size n-by-r, while the Σ factors are stored as sparse diagonal r-by-r matrices. As in A.1 for the Allen–Cahn equation, the square Hadamard power of W(n) is factorized as (W(n))∘2=U∘2Σ∘2V∘2⊤, where U∘2, V∘2∈Rn×r2, and Σ∘2 is a sparse diagonal r2-by-r2 matrix.

We call the operations GW=UΣ, G(n)=U(n)Σ(n), G(n-1)=U(n-1)Σ(n-1), and G⊙=U⊙Σ⊙. We point out that Mm⊤Mm is a symmetric sparse banded matrix with bandwidth 2. This implies that the number of nonzero elements is 2(n-2)+2(n-1)+n=5n-6≪n2, which allows for efficient matrix-matrix products.

Refer to Table 10 for details on the computational costs for evaluating (21) in low-rank format.Table 10 Asymptotic complexities for FKPPE cost function

Product	Factor sizes	Notes on structure and storage	Cost	
12TrGW⊤(Mm⊤Mm)GW	O(nr2)	
GW⊤(Mm⊤Mm)	n×r, n×n	Mm⊤Mm sparse banded	O(nr)	
5n-6 nonzero coefficients	
-Tr(V⊤V(n-1))((G(n-1))⊤(Mp⊤Mm)GW)	O(nr2+r3)	
V⊤V(n-1)	n×r, n×r		O(nr2)	
(G(n-1))⊤(Mp⊤Mm)GW	n×r, n×n	Mp⊤Mm sparse banded	O(nr2)	
(V⊤V(n-1))((G(n-1))⊤(Mp⊤Mm)GW	r×r, r×r	5n-6 nonzero coefficients	O(r3)	
-2hTr((G(n))⊤MmGW)(V⊤RωV(n))	O(nr2+r3)	
(G(n))⊤MmGW	n×r, n×n, n×r	Mm sparse banded	O(nr2)	
V⊤Rω	n×r, n×n	3n-2 nonzero coefficients	O(nr)	
(V⊤Rω)V(n)	r×n, n×r	Rω sparse diagonal	O(nr2)	
((G(n))⊤MmGW)(V⊤RωV(n))	r×r	n nonzero coefficients	O(r3)	
2hTr(G⊙⊤MmGW)(V⊤RωV⊙)	O(nr3+r5)	
G⊙⊤MmGW	n×r2, n×n, n×r	Mm sparse banded	O(nr2)	
(V⊤Rω)V⊙	r×n, n×r2	3n-2 nonzero coefficients	O(nr3)	
(G⊙⊤MmGW)(V⊤RωV⊙)	r2×r, r×r2	Rω sparse diagonal	O(r5)	

Gradient

The Euclidean gradient of the FKPPE cost function F(W) isG=MmW-MpW(n-1)+2h(W(n))∘2-W(n)Rω⊤Mm.

In low-rank format we have G=UGΣGVG⊤, whose factors areUG=(Mm⊤Mm)U(Mm⊤Mp)U(n-1)Mm⊤U∘2Mm⊤U(n),ΣG=blkdiagΣ,-Σ(n-1),2hΣ∘2,-2hΣ(n),VG=VV(n-1)RωV∘2RωV(n)⊤.

Table 11 summarizes the asymptotic complexities for the FKPPE factorized gradient.Table 11 Asymptotic complexities for the FKPPE gradient

Product	Factor sizes	Notes on structure and storage	Cost	
(Mm⊤Mm)U	n×n, n×r	Mm⊤Mm and Mm⊤Mp sparse banded	O(nr)	
(Mm⊤Mp)U(n-1)	n×n, n×r	5n-6 nonzero coefficients	O(nr)	
Mm⊤U∘2	n×n, n×r2	Mm sparse banded	O(nr2)	
Mm⊤U(n)	n×n, n×r	3n-2 nonzero coefficients	O(nr)	
RωV∘2	n×n, n×r2	Rω sparse diagonal	O(nr2)	
RωV(n)	n×n, n×r	n nonzero coefficients	O(nr)	

Hessian

The discretized Euclidean Hessian of (21) isHessF(W)[η]=(Mm⊤Mm)η.

The factored form of the discretized Euclidean Hessian isHW[η]=(Mm⊤Mm)UηSηVη⊤.

where η=UηSηVη⊤ is a tangent vector in TWMr, in a SVD-like format. The only operation needed is the product (Mm⊤Mm)Uη, whose cost is O(nr).

Derivation of the Preconditioner

As mentioned in Sect. 4, the tCG trust-region subsolver can be preconditioned with the inverse of (10). However, inverting the matrix HX directly would be too computationally expensive, taking O(n6) in this case. A suitable preconditioner can be used to solve this problem, thereby reducing the number of iterations required by the tCG solver. This appendix provides the derivation of such a preconditioner.

Applying the Preconditioner

In practice, applying the preconditioner in X∈Mr means solving (without explicitly inverting the matrix) for ξ∈TXM the system24 HXvec(ξ)=vec(η),

where HX is defined in (10) and η∈TXM is a known tangent vector. This equation is equivalent toPX(hx2(Aξ+ξA))=η.

Using definition (7) of the orthogonal projector onto TXMr, we obtainPU(Aξ+ξA)PV+PU⊥(Aξ+ξA)PV+PU(Aξ+ξA)PV⊥=η,

which is equivalent to the system25 PU(Aξ+ξA)PV=PUηPV,PU⊥(Aξ+ξA)PV=PU⊥ηPV,PU(Aξ+ξA)PV⊥=PUηPV⊥.

The main difference w.r.t. [75] is that here, in general, the tangent vectors are not symmetric. Using the matrix representations (6) of the tangent vectors ξ and η at X=UΣV⊤ξ=UMξV⊤+UpξV⊤+U(Vpξ)⊤,η=UMηV⊤+UpηV⊤+U(Vpη)⊤,

with Mξ∈Rr×r, Upξ∈Rm×r, Vpξ∈Rn×r such that (Upξ)⊤U=(Vpξ)⊤V=0. Analogously, for the tangent vector η we have the constraints Mη∈Rr×r, Upη∈Rm×r, Vpη∈Rn×r such that (Upη)⊤U=(Vpη)⊤V=0, also called gauging conditions [72, §9.2.3].

After some manipulations (see Appendix C.3), system (25) can be written as26 U⊤AUMξ+U⊤AUpξ+MξV⊤AV+(Vpξ)⊤AV=Mη,PU⊥AUMξ+PU⊥AUpξ+UpξV⊤AV=Upη,MξV⊤APV⊥+U⊤AU(Vpξ)⊤+(Vpξ)⊤APV⊥=(Vpη)⊤,

where Mξ, Upξ, and Vpξ are the unknown matrices.

The solution flow of system (26) is as follows. From the second and the third equations of (26), we get Upξ and Vpξ depending on Mξ, then we insert the expressions obtained in the first equation to get Mξ.

We introduce orthogonal matrices Q and Q~ to diagonalize U⊤AU and V⊤AV, respectively,27 D=QU⊤AUQ⊤,D~=Q~V⊤AVQ~⊤,

and use them to define the following matricesU^=UQ⊤,V^=VQ~⊤,M^ξ=QMξQ~⊤,M^η=QMηQ~⊤,U^pξ=UpξQ~⊤,V^pξ=VpξQ⊤.

With these transformations, the first equation in (26) becomes (see Appendix C.3.1 for the details)28 DM^ξ+M^ξD~+U^⊤AU^pξ+(V^pξ)⊤AV^=M^η.

By using the same transformations, we can also rewrite the second equation in (26) as (see Appendix C.3.2 for the details)29 PU⊥(AU^M^ξ+AU^pξ+U^pξD~)=U^pη,

with the condition U^⊤U^pξ=0. The ith column of this equation is4PU⊥(A+d~iI)U^pξ(:,i)=U^pη(:,i)-PU⊥AU^M^ξ(:,i),U^⊤U^pξ(:,i)=0,

where d~i, for i=1,…,n, are the diagonal entries of D~. We rewrite this equation as a saddle-point system30 A+d~iIU^U^⊤0U^pξ(:,i)y=U^pη(:,i)-PU⊥AU^M^ξ(:,i)0,

for all y∈Rr. This saddle-point system can be efficiently solved with the techniques described in Sect. C.2.

Let us defineTi:=A+d~iIU^U^⊤0,b1i:=U^pη(:,i)0,b2i:=-PU⊥AU^M^ξ(:,i)0.

The solution of (30) is given byU^pξ(:,i)=(Ti-1b1i+Ti-1b2i)(1:n),

where the notation (1 : n) means that we only keep the first n entries of the vector. In other terms, we have31 U^pξ(:,i)=Ti-1(U^pη(:,i))-Ti-1(PU⊥AU^)M^ξ(:,i).

Here, Ti-1 denotes solving for U^pξ(:,i) the ith saddle-point system, corresponding to (30).

For the third equation in (26), we proceed analogously as above; see Appendix C.3.3 for the details. After some manipulations, we obtain the saddle-point system32 A+diIV^V^⊤0V^pξ(:,i)z=V^pη(:,i)-PV⊥AV^M^ξ⊤(:,i)0,

for all z∈Rr. The solution is33 V^pξ(:,i)=T~i-1(V^pη(:,i))-T~i-1(PV⊥AV^)M^ξ⊤(:,i).

Here, T~i-1 denotes solving for V^pξ(:,i) the ith saddle-point system, corresponding to (32).

We now go back to the first equation in (26), in its form given in (28). To treat the term U^⊤AU^pξ appearing in (28), let us define the vectorsvi=U^⊤ATi-1(U^pη(:,i)),wi=U^⊤ATi-1(PU⊥AU^)M^ξ(:,i).

We emphasize that the vector vi is known, while wi is not because M^ξ is unknown. With these definitions, and (31), one can easily verify that34 U^⊤AU^pξ=v1-w1,…,vr-wr.

Similarly, for treating the term V^⊤AV^pξ, we define the vectorsv~i=V^⊤AT~i-1(V^pη(:,i)),w~i=V^⊤AT~i-1(PV⊥AV^)M^ξ⊤(:,i),

then35 V^⊤AV^pξ=v~1-w~1,…,v~r-w~r.

Inserting (34) and (35) into (28), we obtainDM^ξ+M^ξD~+v1-w1,…,vr-wr+v~1-w~1,…,v~r-w~r⊤=M^η.

The vectors wi and w~i contain the unknown matrix M^ξ, so we leave them on the left-hand side, while since vi and v~i are known, we move them to the right-hand side. By letting Ki:=d~iIr-U^⊤ATi-1(PU⊥AU^) and K~i:=diIr-V^⊤AT~i-1(PV⊥AV^) for i=1,…,r, we get36 K1M^ξ(:,1),…,KrM^ξ(:,r)+M^ξ(1,:)K~1⊤⋮M^ξ(r,:)K~r⊤=R,

where the matrix on the right-hand side is defined byR:=M^η-v1,…,vr-v~1,…,v~r⊤.

Now, we need to isolate M^ξ in (36). We vectorize the first term on the left-hand side of (36) and getvecK1M^ξ(:,1)⋯KrM^ξ(:,r)=Kvec(M^ξ),

where K=blkdiag(K1,…,Kr), a block-diagonal matrix with the Ki on the main diagonal. Vectorizing the second term on the left-hand side of (36), we obtainvecM^ξ(1,:)K~1⊤⋮M^ξ(r,:)K~r⊤=ΠvecK~1M^ξ⊤(1,:)⋯K~rM^ξ⊤(r,:)=ΠK~vec(M^ξ⊤),

where K~=blkdiag(K~1,…,K~r)∈Rr2×r2 is a block-diagonal matrix, and Π is the perfect shuffle matrix defined by vec(X⊤)=Πvec(X) [28, 73]. Wrapping it up, from (36) we obtain the vectorized equation37 K+ΠK~Πvec(M^ξ)=vec(R).

The matrix K+ΠK~Π is of size r2-by-r2, so it can be efficiently inverted if the rank r is really low. We solve this equation for M^ξ, and then we use (31) and (33) to find U^pξ and V^pξ, respectively. Finally, undoing the transformations done by Q and Q~, we find the components of ξMξ=Q⊤M^ξQ~,Upξ=U^pξQ~,Vpξ=V^pξQ,

and thus the tangent vector ξ such that (24) is satisfied.

Efficient Solution of the Saddle-Point System

We use a Schur complement idea to efficiently solve the saddle-point problem (30) and invert the Ti. See the techniques described in [8]. LetBi=U^pη(:,i)PU⊥AU^∈Rn×r+1.

The system Ti(X)=Bi can be solved for X by eliminating the (negative) Schur complement Si=U^⊤(A+d~iI)-1U^. This gives38 Ni=Si-1(U^⊤(A+d~iI)-1Bi),

39 Xi=(A+d~iI)-1Bi-(A+d~iI)-1U^Ni.

We use Cholesky factorization to solve the system (38). For solving (39), we use a sparse solver for (A+d~iI)-1Bi and (A+d~iI)-1U^, for example, MATLAB backslash (A+d~iI)\Bi and (A+d~iI)\U^. Equation (37) is a linear system of size r2.

Algebraic Manipulations

First Equation in (25)

The first equation in system (25) becomesPU(A(UMξV⊤+UpξV⊤+U(Vpξ)⊤)+(UMξV⊤+UpξV⊤+U(Vpξ)⊤)A)PV=PU(UMηV⊤+UpηV⊤+U(Vpη)⊤)PV.

Using the gauging conditions (Upξ)⊤U=(Vpξ)⊤V=0 and (Upη)⊤U=(Vpη)⊤V=0, we obtain

Left-multiplying by U⊤ and right-multiplying by V we get the first equation in (26), i.e.,U⊤AUMξ+U⊤AUpξ+MξV⊤AV+(Vpξ)⊤AV=Mη.

With the transformations introduced in (27), it becomesQU⊤AUQ⊤QMξQ~⊤+QU⊤AUpξQ~⊤+QMξQ~⊤Q~V⊤AVQ~⊤+Q(Vpξ)⊤AVQ~⊤=QMηQ~⊤,DM^ξ+M^ξD~+U^⊤AU^pξ+(V^pξ)⊤AV^=M^η.

which is the first equation in system (26) in the form (28).

Second Equation in (25)

Analogously for the second equation in (25), i.e.,PU⊥(Aξ+ξA)PV=PU⊥ηPV,PU⊥(A(UMξV⊤+UpξV⊤+U(Vpξ)⊤)+(UMξV⊤+UpξV⊤+U(Vpξ)⊤)A)PV=PU⊥(UMηV⊤+UpηV⊤+U(Vpη)⊤)PV

we obtain

then, using PU⊥=I-UU⊤,

Right-multiplying by V we get the second equation in system (26), i.e.,PU⊥AUMξ+PU⊥AUpξ+UpξV⊤AV=Upη.

With the transformations introduced in (27), this becomesPU⊥AUQ⊤QMξQ~⊤+PU⊥AUpξQ~⊤+UpξQ~⊤Q~V⊤AVQ~⊤=UpηQ~⊤,PU⊥AU^M^ξ+PU⊥AU^pξ+U^pξD~=U^pη,PU⊥(AU^M^ξ+AU^pξ+U^pξD~)=U^pη,

which is the second equation in (25) in the form (29).

Third Equation in (25)

The third equation in (25)PU(Aξ+ξA)PV⊥=PUηPV⊥,

becomesPU(A(UMξV⊤+UpξV⊤+U(Vpξ)⊤)+(UMξV⊤+UpξV⊤+U(Vpξ)⊤)A)PV⊥=PU(UMηV⊤+UpηV⊤+U(Vpη)⊤)PV⊥.PUA(UMξV⊤+UpξV⊤+U(Vpξ)⊤)PV⊥+PU(UMξV⊤+UpξV⊤+U(Vpξ)⊤)APV⊥=PU(UMηV⊤+UpηV⊤+U(Vpη)⊤)PV⊥.

Using the gauging conditions (Upξ)⊤U=(Vpξ)⊤V=0, we obtain

Left-multiplying by U⊤ and noting that (Vpη)⊤PV⊥=(Vpη)⊤, we obtain the third equation in system (26), namely,MξV⊤APV⊥+U⊤AU(Vpξ)⊤+(Vpξ)⊤APV⊥=(Vpη)⊤.

Acknowledgements

The authors would like to thank the anonymous referee whose valuable comments contributed to improving an earlier version of this paper. M.S. would also like to thank Jhih-Huang Li for the occasional discussions, which significantly contributed to overcoming obstacles encountered while working on this project.

Funding

Open access funding provided by University of Geneva The work of M.S. was supported by the National Center for Theoretical Sciences and the National Science and Technology Council of Taiwan (R.O.C.), under the contracts 111-2124-M-002-014- and 112-2124-M-002-009-. The work of B.V. was supported by the Swiss National Science Foundation (grant number 192129).

Data Availability

The code and datasets generated and analyzed during the current study are available in the PrecRTR repository, https://github.com/MarcoSutti/PrecRTR.

Declarations

Conflict of interest

The authors confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. The authors further confirm that the order of authors listed in the manuscript has been approved by all of us. The authors also hereby certify that there is not any actual or potential conflict of interest.

1 Not reported here, but reproducible with the distributed code.

2 “Stochastic” might be too big of a term since no Brownian motion is involved. It is just a PDE with random coefficients for the initial condition and the reaction rate.

3 Sometimes known as forward difference matrix.

4 We use the MATLAB notation (:,i) to denote the ith column extraction from a matrix.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Absil PA Baker CG Gallivan KA Trust-region methods on riemannian manifolds Found. of Comput. Math. 2007 7 303 330 10.1007/s10208-005-0179-9
Absil, P.A., Baker, C.G., Gallivan, K.A.: Trust-region methods on riemannian manifolds. Found. of Comput. Math. 7, 303–330 (2007). 10.1007/s10208-005-0179-9
2. Absil PA Mahony R Sepulchre R Optimization Algorithms on Matrix Manifolds 2008 Princeton, NJ Princeton University Press
Absil, P.A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ (2008)
3. Absil PA Malick J Projection-like retractions on matrix manifolds SIAM J. Optim. 2012 22 1 135 158 10.1137/100802529
Absil, P.A., Malick, J.: Projection-like retractions on matrix manifolds. SIAM J. Optim. 22(1), 135–158 (2012). 10.1137/100802529
4. Absil PA Oseledets IV Low-rank retractions: a survey and new results Comput. Optim. Appl. 2015 62 1 5 29 10.1007/s10589-014-9714-4
Absil, P.A., Oseledets, I.V.: Low-rank retractions: a survey and new results. Comput. Optim. Appl. 62(1), 5–29 (2015). 10.1007/s10589-014-9714-4
5. Allen SM Cahn JW Ground state structures in ordered binary alloys with second neighbor interactions Acta Metall. 1972 20 3 423 433 10.1016/0001-6160(72)90037-5
Allen, S.M., Cahn, J.W.: Ground state structures in ordered binary alloys with second neighbor interactions. Acta Metall. 20(3), 423–433 (1972). 10.1016/0001-6160(72)90037-5
6. Allen SM Cahn JW A correction to the ground state of FCC binary ordered alloys with first and second neighbor pairwise interactions Scr. Metall. 1973 7 12 1261 1264 10.1016/0036-9748(73)90073-2
Allen, S.M., Cahn, J.W.: A correction to the ground state of FCC binary ordered alloys with first and second neighbor pairwise interactions. Scr. Metall. 7(12), 1261–1264 (1973). 10.1016/0036-9748(73)90073-2
7. Beneš M Chalupecký V Mikula K Geometrical image segmentation by the Allen-Cahn equation Appl. Numer. Math. 2004 51 2 187 205 10.1016/j.apnum.2004.05.001
Beneš, M., Chalupecký, V., Mikula, K.: Geometrical image segmentation by the Allen-Cahn equation. Appl. Numer. Math. 51(2), 187–205 (2004). 10.1016/j.apnum.2004.05.001
8. Benzi M Golub GH Liesen J Numerical solution of saddle point problems Acta. Numer. 2005 14 1 137 10.1017/S0962492904000212
Benzi, M., Golub, G.H., Liesen, J.: Numerical solution of saddle point problems. Acta. Numer. 14, 1–137 (2005). 10.1017/S0962492904000212
9. Billaud-Friess M Falcó A Nouy A A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds BIT Numer. Math. 2022 62 2 387 408 10.1007/s10543-021-00884-x
Billaud-Friess, M., Falcó, A., Nouy, A.: A new splitting algorithm for dynamical low-rank approximation motivated by the fibre bundle structure of matrix manifolds. BIT Numer. Math. 62(2), 387–408 (2022). 10.1007/s10543-021-00884-x
10. Boumal N An Introduction to Optimization on Smooth Manifolds 2023 Cambridge Cambridge University Press
Boumal, N.: An Introduction to Optimization on Smooth Manifolds. Cambridge University Press, Cambridge (2023). 10.1017/9781009166164
11. Boumal N Absil PA Low-rank matrix completion via preconditioned optimization on the Grassmann manifold Linear Algebra Appl. 2015 475 200 239 10.1016/j.laa.2015.02.027
Boumal, N., Absil, P.A.: Low-rank matrix completion via preconditioned optimization on the Grassmann manifold. Linear Algebra Appl. 475, 200–239 (2015). 10.1016/j.laa.2015.02.027
12. Boumal N Mishra B Absil PA Sepulchre R Manopt, a Matlab toolbox for optimization on manifolds J. Mach. Lear. Res. 2014 15 1455 1459
Boumal, N., Mishra, B., Absil, P.A., Sepulchre, R.: Manopt, a Matlab toolbox for optimization on manifolds. J. Mach. Lear. Res. 15, 1455–1459 (2014)
13. Brenner S Scott R The Mathematical Theory of Finite Element Methods. Texts in Applied Mathematics 2007 New York Springer
Brenner, S., Scott, R.: The Mathematical Theory of Finite Element Methods. Texts in Applied Mathematics. Springer, New York (2007)
14. Cai, J.F., Huang, W., Wang, H., Wei, K.: Tensor completion via tensor train based low-rank quotient geometry under a preconditioned metric (2022)
15. Ceruti G Kusch J Lubich C A rank-adaptive robust integrator for dynamical low-rank approximation BIT Numer. Math. 2022 62 4 1149 1174 10.1007/s10543-021-00907-7
Ceruti, G., Kusch, J., Lubich, C.: A rank-adaptive robust integrator for dynamical low-rank approximation. BIT Numer. Math. 62(4), 1149–1174 (2022). 10.1007/s10543-021-00907-7
16. Ceruti G Lubich C An unconventional robust integrator for dynamical low-rank approximation BIT Numer. Math. 2022 62 1 23 44 10.1007/s10543-021-00873-0
Ceruti, G., Lubich, C.: An unconventional robust integrator for dynamical low-rank approximation. BIT Numer. Math. 62(1), 23–44 (2022). 10.1007/s10543-021-00873-0
17. Charous, A., Lermusiaux, P.: Dynamically orthogonal differential equations for stochastic and deterministic reduced-order modeling of ocean acoustic wave propagation. In: OCEANS 2021: San Diego – Porto, pp. 1–7 (2021). 10.23919/OCEANS44145.2021.9705914
18. Charous, A., Lermusiaux, P.F.J.: Stable rank-adaptive dynamically orthogonal runge-kutta schemes (2022)
19. Dieudonné J Foundations of Modern Analysis 1960 New York Academic Press
Dieudonné, J.: Foundations of Modern Analysis. Academic Press, New York (1960)
20. Dobrosotskaya JA Bertozzi AL A wavelet-laplace variational technique for image deconvolution and inpainting IEEE Trans. Image Process. 2008 17 5 657 663 10.1109/TIP.2008.919367 18390372
Dobrosotskaya, J.A., Bertozzi, A.L.: A wavelet-laplace variational technique for image deconvolution and inpainting. IEEE Trans. Image Process. 17(5), 657–663 (2008). 10.1109/TIP.2008.91936718390372
21. Edelman A Arias TA Smith ST The geometry of algorithms with orthogonality constraints SIAM J. Matrix Anal. Appl. 1998 20 2 303 353 10.1137/S0895479895290954
Edelman, A., Arias, T.A., Smith, S.T.: The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl. 20(2), 303–353 (1998). 10.1137/S0895479895290954
22. Feppon F Lermusiaux PFJ Dynamically orthogonal numerical schemes for efficient stochastic advection and lagrangian transport SIAM Rev. 2018 60 3 595 625 10.1137/16M1109394
Feppon, F., Lermusiaux, P.F.J.: Dynamically orthogonal numerical schemes for efficient stochastic advection and lagrangian transport. SIAM Rev. 60(3), 595–625 (2018). 10.1137/16M1109394
23. Feppon F Lermusiaux PFJ A geometric approach to dynamical model order reduction SIAM J. Matrix Anal. Appl. 2018 39 1 510 538 10.1137/16M1095202
Feppon, F., Lermusiaux, P.F.J.: A geometric approach to dynamical model order reduction. SIAM J. Matrix Anal. Appl. 39(1), 510–538 (2018). 10.1137/16M1095202
24. Fisher RA The wave of advance of advantageous genes Ann. Eug. 1937 7 4 355 369 10.1111/j.1469-1809.1937.tb02153.x
Fisher, R.A.: The wave of advance of advantageous genes. Ann. Eug. 7(4), 355–369 (1937). 10.1111/j.1469-1809.1937.tb02153.x
25. Grasedyck L Existence of a low rank or H-matrix approximant to the solution of a Sylvester equation Numer. Linear Algebra with Appl. 2004 11 4 371 389 10.1002/nla.366
Grasedyck, L.: Existence of a low rank or -matrix approximant to the solution of a Sylvester equation. Numer. Linear Algebra with Appl. 11(4), 371–389 (2004). 10.1002/nla.366
26. Gratton S Sartenaer A Toint PL Recursive trust-region methods for multiscale nonlinear optimization SIAM J. Optim. 2008 19 1 414 444 10.1137/050623012
Gratton, S., Sartenaer, A., Toint, P.L.: Recursive trust-region methods for multiscale nonlinear optimization. SIAM J. Optim. 19(1), 414–444 (2008). 10.1137/050623012
27. Helmke U Moore JB Optimization and Dynamical Systems 1994 London Springer-Verlag
Helmke, U., Moore, J.B.: Optimization and Dynamical Systems. Springer-Verlag, London (1994). 10.1007/978-1-4471-3467-1
28. Henderson HV Searle SR The vec-permutation matrix, the vec operator and Kronecker products: a review Linear Multilinear Algebra 1981 9 4 271 288 10.1080/03081088108817379
Henderson, H.V., Searle, S.R.: The vec-permutation matrix, the vec operator and Kronecker products: a review. Linear Multilinear Algebra 9(4), 271–288 (1981). 10.1080/03081088108817379
29. Henson VE Bouman CA Stevenson RL Multigrid methods nonlinear problems: an overview Computational Imaging 2003 US International Society for Optics and Photonics, SPIE 36 48
Henson, V.E.: Multigrid methods nonlinear problems: an overview. In: Bouman, C.A., Stevenson, R.L. (eds.) Computational Imaging, vol. 5016, pp. 36–48. International Society for Optics and Photonics, SPIE, US (2003). 10.1117/12.499473
30. Hundsdorfer W Verwer J Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations 2003 1 Berlin, Heidelberg Springer
Hundsdorfer, W., Verwer, J.: Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, 1st edn. Springer, Berlin, Heidelberg (2003). 10.1007/978-3-662-09017-6
31. Khatri CG Rao CR Solutions to some functional equations and their applications to characterization of probability distributions Sankhyā Ser. A 1968 30 2 167 180
Khatri, C.G., Rao, C.R.: Solutions to some functional equations and their applications to characterization of probability distributions. Sankhyā Ser. A 30(2), 167–180 (1968)
32. Kieri E Lubich C Walach H Discretized dynamical low-rank approximation in the presence of small singular values SIAM J. Numer. Anal. 2016 54 2 1020 1038 10.1137/15M1026791
Kieri, E., Lubich, C., Walach, H.: Discretized dynamical low-rank approximation in the presence of small singular values. SIAM J. Numer. Anal. 54(2), 1020–1038 (2016). 10.1137/15M1026791
33. Kieri E Vandereycken B Projection methods for dynamical low-rank approximation of high-dimensional problems Comput. Methods Appl. Math. 2019 19 1 73 92 10.1515/cmam-2018-0029
Kieri, E., Vandereycken, B.: Projection methods for dynamical low-rank approximation of high-dimensional problems. Comput. Methods Appl. Math. 19(1), 73–92 (2019). 10.1515/cmam-2018-0029
34. Koch O Lubich C Dynamical low-rank approximation SIAM J. Matrix Anal. Appl. 2007 29 2 434 454 10.1137/050639703
Koch, O., Lubich, C.: Dynamical low-rank approximation. SIAM J. Matrix Anal. Appl. 29(2), 434–454 (2007). 10.1137/050639703
35. Kolmogorov, A.N., Petrowsky, I.G., Piskunov, N.S.: Studies of the diffusion with the increasing quantity of the substance; its application to a biological problem. In: O.A. Oleinik (ed.) I.G. Petrowsky Selected Works. Part II: Differential Equations and Probability Theory, Classics of Soviet Mathematics, vol. 5, first edn., chap. 6, pp. 106–132. CRC Press, London (1996). 10.1201/9780367810504
36. Kressner, D.: Advanced numerical analysis (2015). https://www.epfl.ch/labs/anchp/wp-content/uploads/2018/05/AdvancedNA2015.pdf
37. Kressner D Steinlechner M Vandereycken B Preconditioned low-rank riemannian optimization for linear systems with tensor product structure SIAM J. Sci. Comput. 2016 38 4 A2018 A2044 10.1137/15M1032909
Kressner, D., Steinlechner, M., Vandereycken, B.: Preconditioned low-rank riemannian optimization for linear systems with tensor product structure. SIAM J. Sci. Comput. 38(4), A2018–A2044 (2016). 10.1137/15M1032909
38. Kressner D Tobler C Preconditioned low-rank methods for high-dimensional elliptic PDE eigenvalue problems Comput. Methods Appl. Math. 2011 11 3 363 381 10.2478/cmam-2011-0020
Kressner, D., Tobler, C.: Preconditioned low-rank methods for high-dimensional elliptic PDE eigenvalue problems. Comput. Methods Appl. Math. 11(3), 363–381 (2011). 10.2478/cmam-2011-0020
39. Kressner D Tobler C Algorithm 941: Htucker–a matlab toolbox for tensors in hierarchical tucker format ACM Trans. Math. Softw. 2014 40 3 22:1 22:22 10.1145/2538688
Kressner, D., Tobler, C.: Algorithm 941: Htucker–a matlab toolbox for tensors in hierarchical tucker format. ACM Trans. Math. Softw. 40(3), 22:1-22:22 (2014). 10.1145/2538688
40. Kürschner, P.: Efficient low-rank solution of large-scale matrix equations. Ph.D. thesis, Aachen (2016)
41. Laux T Simon TM Convergence of the Allen-Cahn equation to multiphase mean curvature flow Commun. Pure Appl. Anal. 2018 71 8 1597 1647 10.1002/cpa.21747
Laux, T., Simon, T.M.: Convergence of the Allen-Cahn equation to multiphase mean curvature flow. Commun. Pure Appl. Anal. 71(8), 1597–1647 (2018). 10.1002/cpa.21747
42. Le Dret H Lucquin B Partial Differential Equations: Modeling, Analysis and Numerical Approximation 2016 Basel Birkhäuser
Le Dret, H., Lucquin, B.: Partial Differential Equations: Modeling, Analysis and Numerical Approximation. Birkhäuser, Basel (2016)
43. Lee D Lee S Image segmentation based on modified fractional Allen-Cahn equation Math. Probl. Eng. 2019 2019 3980181 10.1155/2019/3980181
Lee, D., Lee, S.: Image segmentation based on modified fractional Allen-Cahn equation. Math. Probl. Eng. 2019, 3980181 (2019). 10.1155/2019/3980181
44. Lee HG Kim J An efficient and accurate numerical algorithm for the vector-valued Allen-Cahn equations Comput. Phys. Commun. 2012 183 10 2107 2115 10.1016/j.cpc.2012.05.013
Lee, H.G., Kim, J.: An efficient and accurate numerical algorithm for the vector-valued Allen-Cahn equations. Comput. Phys. Commun. 183(10), 2107–2115 (2012). 10.1016/j.cpc.2012.05.013
45. Lermusiaux P Evolving the subspace of the three-dimensional multiscale ocean variability: massachusetts bay J. Mar. Syst. 2001 29 1 385 422 10.1016/S0924-7963(01)00025-2
Lermusiaux, P.: Evolving the subspace of the three-dimensional multiscale ocean variability: massachusetts bay. J. Mar. Syst. 29(1), 385–422 (2001). 10.1016/S0924-7963(01)00025-2. (Three-Dimensional Ocean Circulation: Lagrangian measurements and diagnostic analyses)
46. Lermusiaux, P.F.J., Robinson, A.R.: Data assimilation via error subspace statistical estimation. Part i: theory and schemes. Mon. Weather Rev. 127(7), 1385–1407 (1999). https://doi.org/10.1175/1520-0493(1999)1271385:DAVESS2.0.CO;2
47. Li Y Jeong D Choi J Lee S Kim J Fast local image inpainting based on the Allen-Cahn model Digit. Signal Process. 2015 37 65 74 10.1016/j.dsp.2014.11.006
Li, Y., Jeong, D., Choi, J., Lee, S., Kim, J.: Fast local image inpainting based on the Allen-Cahn model. Digit. Signal Process. 37, 65–74 (2015). 10.1016/j.dsp.2014.11.006
48. Lubich C Oseledets IV A projector-splitting integrator for dynamical low-rank approximation BIT Numer. Math. 2014 54 1 171 188 10.1007/s10543-013-0454-0
Lubich, C., Oseledets, I.V.: A projector-splitting integrator for dynamical low-rank approximation. BIT Numer. Math. 54(1), 171–188 (2014). 10.1007/s10543-013-0454-0
49. Massei S Robol L Kressner D Hierarchical adaptive low-rank format with applications to discretized partial differential equations Numer. Linear Algebra Appl. 2022 29 6 e2448 10.1002/nla.2448
Massei, S., Robol, L., Kressner, D.: Hierarchical adaptive low-rank format with applications to discretized partial differential equations. Numer. Linear Algebra Appl. 29(6), e2448 (2022). 10.1002/nla.2448
50. Mishra B Meyer G Bach F Sepulchre R Low-rank optimization with trace norm penalty SIAM J. Optim. 2013 23 4 2124 2149 10.1137/110859646
Mishra, B., Meyer, G., Bach, F., Sepulchre, R.: Low-rank optimization with trace norm penalty. SIAM J. Optim. 23(4), 2124–2149 (2013)
51. Mishra B Sepulchre R Riemannian preconditioning SIAM J. Optim. 2016 26 1 635 660 10.1137/140970860
Mishra, B., Sepulchre, R.: Riemannian preconditioning. SIAM J. Optim. 26(1), 635–660 (2016). 10.1137/140970860
52. Mishra, B., Vandereycken, B.: A Riemannian approach to low-rank algebraic Riccati equations. In: 21st International Symposium on Mathematical Theory of Networks and Systems, pp. 965–968. Groningen, The Netherlands (2014)
53. Murray JD Mathematical Biology I. An Introduction 2002 3 New York, NY Springer
Murray, J.D.: Mathematical Biology I. An Introduction, 3rd edn. Springer, New York, NY (2002). 10.1007/b98868
54. Musharbash E Nobile F Vidličková E Symplectic dynamical low rank approximation of wave equations with random parameters BIT Numer. Math. 2020 60 4 1153 1201 10.1007/s10543-020-00811-6
Musharbash, E., Nobile, F., Vidličková, E.: Symplectic dynamical low rank approximation of wave equations with random parameters. BIT Numer. Math. 60(4), 1153–1201 (2020). 10.1007/s10543-020-00811-6
55. Ngo, T., Saad, Y.: Scaled Gradients on Grassmann Manifolds for Matrix Completion. In: Pereira, F., Burges, C., Bottou, L., Weinberger, K. (eds.) Adv. Neural Inf. Process. Syst., vol. 25. Curran Associates Inc., USA (2012)
56. Nocedal J Wright SJ Numerical Optimization 2006 2 New York, NY Springer
Nocedal, J., Wright, S.J.: Numerical Optimization, 2nd edn. Springer, New York, NY (2006). 10.1007/978-0-387-40065-5
57. Rakhuba M Novikov A Oseledets I Low-rank Riemannian eigensolver for high-dimensional Hamiltonians J. Comput. Phys. 2019 396 718 737 10.1016/j.jcp.2019.07.003
Rakhuba, M., Novikov, A., Oseledets, I.: Low-rank Riemannian eigensolver for high-dimensional Hamiltonians. J. Comput. Phys. 396, 718–737 (2019)
58. Rakhuba M Oseledets I Jacobi-Davidson method on low-rank matrix manifolds SIAM J. Sci. Comput. 2018 40 2 A1149 A1170 10.1137/17M1123080
Rakhuba, M., Oseledets, I.: Jacobi-Davidson method on low-rank matrix manifolds. SIAM J. Sci. Comput. 40(2), A1149–A1170 (2018)
59. Rodgers A Dektor A Venturi D Adaptive integration of nonlinear evolution equations on tensor manifolds J. Sci. Comput. 2022 92 2 39 10.1007/s10915-022-01868-x
Rodgers, A., Dektor, A., Venturi, D.: Adaptive integration of nonlinear evolution equations on tensor manifolds. J. Sci. Comput. 92(2), 39 (2022). 10.1007/s10915-022-01868-x
60. Rodgers A Venturi D Implicit step-truncation integration of nonlinear PDEs on low-rank tensor manifolds J. Sci. Comput. 2022 97 2 33 10.48550/ARXIV.2207.01962
Rodgers, A., Venturi, D.: Implicit step-truncation integration of nonlinear PDEs on low-rank tensor manifolds. J. Sci. Comput. 97(2), 33 (2022). 10.48550/ARXIV.2207.01962. arXiv:2207.01962
61. Sapsis TP Lermusiaux PF Dynamically orthogonal field equations for continuous stochastic dynamical systems Phys. D: Nonlinear Phenom. 2009 238 23 2347 2360 10.1016/j.physd.2009.09.017
Sapsis, T.P., Lermusiaux, P.F.: Dynamically orthogonal field equations for continuous stochastic dynamical systems. Phys. D: Nonlinear Phenom. 238(23), 2347–2360 (2009). 10.1016/j.physd.2009.09.017
62. Shalit, U., Weinshall, D., Chechik, G.: Online Learning in The Manifold of Low-Rank Matrices. In: Lafferty, J., Williams, C., Shawe-Taylor, J., Zemel, R., Culotta, A. (eds.) Adv. Neural Inf. Process. Syst., vol. 23. Curran Associates Inc., USA (2010)
63. Shalit U Weinshall D Chechik G Online learning in the embedded manifold of low-rank matrices J. Mach. Lear. Res. 2012 13 429 458
Shalit, U., Weinshall, D., Chechik, G.: Online learning in the embedded manifold of low-rank matrices. J. Mach. Lear. Res. 13, 429–458 (2012)
64. Simoncini V Computational methods for linear matrix equations SIAM Rev. 2016 58 3 377 441 10.1137/130912839
Simoncini, V.: Computational methods for linear matrix equations. SIAM Rev. 58(3), 377–441 (2016). 10.1137/130912839
65. Steihaug T The conjugate gradient method and trust regions in large scale optimization SIAM J. Numer. Anal. 1983 20 3 626 637 10.1137/0720042
Steihaug, T.: The conjugate gradient method and trust regions in large scale optimization. SIAM J. Numer. Anal. 20(3), 626–637 (1983). 10.1137/0720042
66. Steinlechner M Riemannian optimization for high-dimensional tensor completion SIAM J. Sci. Comput. 2016 38 5 S461 S484 10.1137/15M1010506
Steinlechner, M.: Riemannian optimization for high-dimensional tensor completion. SIAM J. Sci. Comput. 38(5), S461–S484 (2016)
67. Sutti, M.: Riemannian algorithms on the stiefel and the fixed-rank manifold. Ph.D. thesis, University of Geneva (2020). ID: unige:146438
68. Sutti, M., Vandereycken, B.: RMGLS: A MATLAB algorithm for Riemannian multilevel optimization. Available online (2020). 10.26037/yareta:zara3a5aivcsfk6uhq4oovjxhe
69. Sutti M Vandereycken B Riemannian multigrid line search for low-rank problems SIAM J. Sci. Comput. 2021 43 3 A1803 A1831 10.1137/20M1337430
Sutti, M., Vandereycken, B.: Riemannian multigrid line search for low-rank problems. SIAM J. Sci. Comput. 43(3), A1803–A1831 (2021). 10.1137/20M1337430
70. Toint, P.L.: Towards an efficient sparsity exploiting Newton method for minimization. In: Sparse Matrices and Their Uses, pp. 57–88. Academic Press, London, England (1981)
71. Ueckermann M Lermusiaux P Sapsis T Numerical schemes for dynamically orthogonal equations of stochastic fluid and ocean flows J. Comput. Phys. 2013 233 272 294 10.1016/j.jcp.2012.08.041
Ueckermann, M., Lermusiaux, P., Sapsis, T.: Numerical schemes for dynamically orthogonal equations of stochastic fluid and ocean flows. J. Comput. Phys. 233, 272–294 (2013). 10.1016/j.jcp.2012.08.041
72. Uschmajew, A., Vandereycken, B.: Geometric methods on low-rank matrix and tensor manifolds, chap. 9, pp. 261–313. Springer International Publishing, Cham (2020). 10.1007/978-3-030-31351-7_9
73. Van Loan CF The ubiquitous Kronecker product J. Comput. and Appl. Math. 2000 123 1–2 85 100 10.1016/S0377-0427(00)00393-9
Van Loan, C.F.: The ubiquitous Kronecker product. J. Comput. and Appl. Math. 123(1–2), 85–100 (2000)
74. Vandereycken B Low-rank matrix completion by riemannian optimization SIAM J. Optim. 2013 23 2 1214 1236 10.1137/110845768
Vandereycken, B.: Low-rank matrix completion by riemannian optimization. SIAM J. Optim. 23(2), 1214–1236 (2013). 10.1137/110845768
75. Vandereycken B Vandewalle S A riemannian optimization approach for computing low-rank solutions of Lyapunov equations SIAM J. Matrix Anal. Appl. 2010 31 5 2553 2579 10.1137/090764566
Vandereycken, B., Vandewalle, S.: A riemannian optimization approach for computing low-rank solutions of Lyapunov equations. SIAM J. Matrix Anal. Appl. 31(5), 2553–2579 (2010). 10.1137/090764566
76. Wen Z Goldfarb D A line search multigrid method for large-scale nonlinear optimization SIAM J. Optim. 2009 20 3 1478 1503 10.1137/08071524X
Wen, Z., Goldfarb, D.: A line search multigrid method for large-scale nonlinear optimization. SIAM J. Optim. 20(3), 1478–1503 (2009). 10.1137/08071524X
77. Yang X Feng JJ Liu C Shen J Numerical simulations of jet pinching-off and drop formation using an energetic variational phase-field method J. Comput. Phys. 2006 218 1 417 428 10.1016/j.jcp.2006.02.021
Yang, X., Feng, J.J., Liu, C., Shen, J.: Numerical simulations of jet pinching-off and drop formation using an energetic variational phase-field method. J. Comput. Phys. 218(1), 417–428 (2006). 10.1016/j.jcp.2006.02.021
78. Yoon S Jeong D Lee C Kim H Kim S Lee HG Kim J Fourier-spectral method for the phase-field equations Mathematics 2020 8 8 1385 10.3390/math8081385
Yoon, S., Jeong, D., Lee, C., Kim, H., Kim, S., Lee, H.G., Kim, J.: Fourier-spectral method for the phase-field equations. Mathematics 8(8), 1385 (2020). 10.3390/math8081385
