---
abstract: 'In this paper, we extend the dynamical low-rank approximation method to the space of finite signed measures. Under this framework, we derive stochastic low-rank dynamics for stochastic differential equations (SDEs) coming from classical stochastic dynamics or unraveling of Lindblad quantum master equations. We justify the proposed method by error analysis and also numerical examples for applications in solving high-dimensional SDE, stochastic Burgers’ equation, and high-dimensional [Lindblad]{} equation.'
address:
- 'Department of Mathematics, Duke University, Box 90320, Durham, NC 27708 USA'
- 'Departments of Mathematics, Physics, and Chemistry, Duke University, Box 90320, Durham, NC 27708 USA'
author:
- Yu Cao
- Jianfeng Lu
bibliography:
- 'ref.bib'
title: 'Stochastic dynamical low-rank approximation method'
---

Dynamical low-rank approximation, stochastic differential equation, [Lindblad]{} equation, model reduction.

Introduction
============

Many problems in computational physics are challenging to solve due to *curse of dimensionality*, such as high dimensional master equations and many-body quantum dynamics. In attempt to resolve the difficulty, many ideas have been proposed: model reduction method [@antoulas_survey_2001; @Chatterjee_00_pod], Monte Carlo method [@Gisin92; @Dalibard92; @Tully], [*etc*.]{} In many situations, several methods of dimension reduction need to be combined together. For instance, after applying Monte Carlo method to some deterministic dynamics by simulating a stochastic differential equation (SDE) instead, the dimension of that SDE may still be very large. Then it is attractive to further apply model reduction method in order to capture the main dynamical flows. This is our motivation to study model reduction method for high-dimensional SDEs arising from high-dimensional PDEs or matrix ODEs. In particular, our main motivation comes from two important physical systems: Fokker-Planck equation [@Pavliotis] and [Lindblad]{} equation [@Lindblad76; @Gorini76; @Breuer].

[Fokker-Planck]{} equations and [Lindblad]{} equations are the governing master equations to describe system evolution for open classical and quantum systems respectively under Markovian approximation. Both are challenging to solve when the dimension becomes large. To resolve this problem, it is standard to consider Monte Carlo (particle) method based on stochastic differential equations, with statistical average of sample trajectories to obtain the quantity of interest. For the quantum case, such methods are known as “unraveling” and “stochastic wave-function method” for Lindblad equation [@Gisin92; @Dalibard92].

More specifically, suppose we would like to solve the Fokker-Planck equation $\partial_t \mu_t = {\mathcal{A}}_t \mu_t$ where $\mu_t$ is the probability distribution (or measure to be more general) and ${\mathcal{A}}_t$ is a time-parametrized operator mapping a probability distribution to its tangent space. In the particle based methods, one simulates a SDE $X_t$ with infinitesimal generator ${\mathcal{A}}_t^{*}$ (the adjoint operator of ${\mathcal{A}}_t$) and with initial condition $X_0$ drawn from $\mu_0$ and thus the distribution of $X_t$ is exactly $\mu_t$. Similarly, as the quantum analog, [Lindblad]{} equation has the form $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{\rho}}}= {\mathcal{L}}({{\boldsymbol{\rho}}})$ where ${{\boldsymbol{\rho}}}$ denotes the density matrix of a quantum system and ${\mathcal{L}}$ is the generator for a completely positive dynamical semigroup [@Lindblad76; @Gorini76]. [[One may solve it by sampling a SDE $X_t$ such that ${\mathbb{E}}\bigl[X_t X^{\dagger}_t\bigr]$ is exactly the solution of [Lindblad]{} equation (see Lemma \[lemma::lb\_unravel\] below for more details).]{}]{} [[There are various choices of the SDEs, [*e.g.*]{}, quantum state diffusion (QSD) [@Gisin92] and linear quantum state diffusion (LQSD) [@Brun00].]{}]{} While it is also possible to use other stochastic processes such as jumping process [@Dalibard92; @Brun00], we will limit the scope of our consideration to Monte Carlo methods based on diffusion processes.

The Monte Carlo method for both [Fokker-Planck]{} equation and [Lindblad]{} equation can be described under the same framework:

Given a ${{\mathbb{C}}^{{n}}}$-valued SDE $X_t$, one would like to approximate ${\mathbb{E}}\left[f(X_t)\right] \equiv \int f\ {\,\mathrm{d}}\mu_t$ for a collection of prescribed functions $f \in {\mathsf{F}}$, where $\mu_t$ is the distribution of $X_t$.

In the case of Fokker-Planck equation, ${\mathsf{F}}$ could be a collection of smooth functions; in the case of [Lindblad]{} equation, ${\mathsf{F}}$ could be a singleton set $\left\{ f(x) = x x^{\dagger} \right\}$ where $x\in {{\mathbb{C}}^{{n}}}$. In the sequel, we shall consider $f(x) = x x^{\dagger}$ only, which turns out to be an interesting and useful choice: in the case of Fokker-Planck equation, choosing such $f$ means one would like to calculate the second moment of measure $\mu_t$; in the case of [Lindblad]{} equation, choosing such $f$ means one would like to compute the density matrix.

To reduce the computational complexity, a popular approach is model reduction, that is, to retrieve the dynamics by only capturing the evolution of a lower-dimensional object. For our case, there are two directions:

(i) find a low-rank approximation for $X_t$, or,

(ii) find a low-rank approximation for $\mu_t$.

[[In the literature, there are several methods taking the first approach, in the flavor of [Karhunen-Lo[\` e]{}ve]{} expansion (KLE): for instance, proper orthogonal method (POD) [@Chatterjee_00_pod; @sapsis_dynamically_2011], dynamical orthogonal (DO) method [@sapsis_dynamically_2011; @Sapsis09; @Sapsis12; @Musharbash15] and dynamical bi-orthogonal method (DyBO) [@Cheng13_theory; @Cheng13_complexity]. It is clear that in the above framework for our cases, the realization of randomness in $X_t$ is not important, whereas the distribution $\mu_t$ is the key for accurate approximation.]{}]{} Hence, it is natural to consider the low-rank approximation for $\mu_t$, [*i.e.*]{}, on the space of probability measures. Then the problem is formulated as follows:

[[ Given a collection of prescribed test functions ${\mathsf{F}}$ and the time evolution equation of probability measures $\partial_t \mu_t = {\mathcal{A}}_t \mu_t$ on ${{\mathbb{C}}^{{n}}}$, one would like to find low-rank approximation $\mu_{{LR},t}\approx \mu_t$ such that $$\sup_{f\in {\mathsf{F}}} {\left\lVert\int f\ {\,\mathrm{d}}\mu_t - f\ {\,\mathrm{d}}\mu_{{LR},t} \right\rVert}$$ is small. ]{}]{}

As a remark, in [Section \[sec::sdlr\]]{}, we shall use the space of finite signed measures, instead of probability measures to avoid the technicality; please see the discussion in [Section \[sec::sdlr\]]{} for details.

Our work is motivated by extending the (deterministic) dynamical low-rank approximation, introduced by Koch and Lubich in [@Lubich07] for matrix ODEs, to the stochastic case. The main idea in the dynamical low-rank approximation has been illustrated in the context of matrix ODE [@Lubich07], summarized in the next paragraph.

Consider a matrix ODE system $M(t)\in {\mathbb{C}}^{n\times n}$, $$\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{M}(t) = F(t, M(t)).$$ The dynamical low-rank approximation method in [@Lubich07] consists of two steps. Firstly, identify a sub-manifold $\mathcal{M}_{r}\subset {\mathbb{C}}^{n\times n}$ and approximate the matrix ODE solution $M(t)$ by $M_{{LR}}(t)\in \mathcal{M}_r$ for all $t$; secondly, the time-evolution $M_{{LR}}(t)$ is given by $$\label{eqn::tsp_ode}
  \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} M_{{LR}}(t) = \operatorname*{arg\,min}_{v \in {\mathcal{T}}_{M_{{LR}}(t)} {\mathcal{M}}_r} d\bigl(v, F(t, M_{{LR}}(t))\bigr),$$ where ${\mathcal{T}}_{M_{{LR}}(t)} {\mathcal{M}}_r$ is the tangent space of $\mathcal{M}_r$ at the current location $M_{{LR}}(t)$ and $d(v_1, v_2) := {\left\lVertv_1 - v_2\right\rVert}$ is a metric on tangent space; thus the evolution is constructed as close as possible to the solution of matrix ODE by projecting $F(t, M_{{LR}}(t))$ onto the tangent space ${\mathcal{T}}_{M_{{LR}}(t)}{\mathcal{M}}_r$.

In our proposed method, we adopt this idea to the space of finite signed measures on ${{\mathbb{C}}^{{n}}}$ with bounded second moment, denoted by ${\mathcal{M}}$. The subspace ${\mathcal{M}}_r$ in this case is defined as the space of finite signed measures supported on a linear subspace of ${{\mathbb{C}}^{{n}}}$ with dimension at most $r$. Then, we hope to approximate $\mu_t$ by $\mu_{{LR},t}\in
{\mathcal{M}}_r$. The time-evolution equation of the low-rank approximation is given by $$\partial_t \mu_{{LR},t} \equiv {\mathcal{A}}_{{LR},t} \mu_{{LR},t} := \operatorname*{arg\,min}_{\nu \in {\mathcal{T}}_{\mu_{{LR},t}}{\mathcal{M}}_r}  d_{{\mathsf{F}}}\left(\nu,\ {\mathcal{A}}_{t} \mu_{{LR},t}\right),$$ where $d_{{\mathsf{F}}}$ is a [pseudometric]{} defined in [Eq. ]{} below and ${\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r$ is the tangent space of ${\mathcal{M}}_r$ at $\mu_{{LR},t}$. We will refer this method as the [stochastic dynamical [low-rank]{} approximation]{} method (or [SDLR]{} in abbreviation).

As a concrete example, let ${\mathsf{F}}= \left\{ f(x) = x x^{\dagger} \right\}$ be a singleton set, consisting only one test function (which maps to ${\mathbb{C}}^{n\times n}$). Assume that the time evolution equation $\partial_t \mu_t  = {\mathcal{A}}_{t} \mu_t$ is the Fokker-Planck equation of a SDE of the form $$\label{eqn::sde}
{\,\mathrm{d}}X_t = a(X_t, t){\,\mathrm{d}}t + \sum_{j=1}^{N} b_j(X_t, t) {\,\mathrm{d}}W_j,$$ where $X_t\in {{\mathbb{C}}^{{n}}}$, $a$ and $b_j$ are functions ${{\mathbb{C}}^{{n}}}\times {[0,T]}\rightarrow {{\mathbb{C}}^{{n}}}$ and $W_j$ are independent real-valued standard Brownian motions. With some additional assumptions [[and restrictions]{}]{}, one could obtain the low-rank dynamics given by [Eq. ]{} (or equivalently [Eq. ]{}). The details are given in [Theorem \[thm::fx=xx\]]{}, which is one of the main results in this paper.

As already mentioned above, the stochastic dynamical low-rank approximation is also motivated by developing efficient methods for the Lindblad quantum master equations. In that context, the deterministic low-rank approximation has been studied by Le Bris and Rouchon to find low-rank approximation of [Lindblad]{} equation [@LeBris13]. In the subsequent work [@LeBris15], Le Bris, Rouchon and Roussel also introduced an unraveling scheme for the low-rank quantum master equation obtained in [@LeBris13]. The unraveling of Lindblad equations and its connection with the low-rank approximation will be discussed in Section \[sec::discuss::unravel\]. In particular, as another main result of this paper, we establish a commuting diagram of unraveling and low-rank approximation, with the proposed [SDLR]{} method.

The rest of the paper is organized as follows. In [Section \[sec::sdlr\]]{}, we shall formulate the [stochastic dynamical [low-rank]{} approximation]{} method in the space of finite signed measures. Then in [Section \[sec::example\]]{}, we will provide a concrete example, in which a low-rank dynamics for [Fokker-Planck]{} equation is derived, as well as the low-rank dynamics of the SDE for that [Fokker-Planck]{} equation. The comparison of our method and DO method will also be given at the end of [Section \[sec::example\]]{}. [[ Consistency of our low-rank approximation and error analysis will be provided in [Section \[sec::error\]]{}. In [Section \[sec::discuss::unravel\]]{}, we shall establish the connection between the action of dynamical low-rank approximation and the action of unraveling. We will prove a commuting relation between them.]{}]{} Then numerical results will be presented in [Section \[sec::numerics\]]{} to demonstrate the performance. In [Section \[sec::conclusion\]]{}, we will give a brief summary and some potential follow-up work.

Throughout this paper, ${\left\lVert\cdot\right\rVert}_{{HS}}$ means Hilbert-Schmidt norm (or Frobenius norm as in linear algebra). Bold letters, like ${{\boldsymbol{U}}}$ and ${{\boldsymbol{G}}}$ [*etc*.]{}, represent matrices (linear operators). Also, ${\left[\cdot, \cdot\right]_{-}}$ is commutator and ${\left[\cdot, \cdot\right]_{+}}$ is anti-commutator.

Stochastic dynamical low-rank approximation method {#sec::sdlr}
==================================================

As we recalled in the introduction, the dynamical low-rank approximation method [@Lubich07], developed for deterministic ODE dynamics, involves the identification of an approximate sub-manifold and projection onto the tangent space by solving a minimization problem. In this section, we will adopt this idea to formulate the dynamical low-rank approximation in the space of finite signed measures on ${{\mathbb{C}}^{{n}}}$ with bounded second moment. This low-rank approximation method offers an abstract framework, for instance, to approximate both [Fokker-Planck]{} and [Lindblad]{} equations via low-rank dynamics combined with the particle methods. Thus the proposed method will be named *[stochastic dynamical [low-rank]{} approximation]{}* ([SDLR]{}) method. A concrete example and corresponding low-rank dynamics will be given in [Section \[sec::example\]]{} below.

Problem setup and low-rank approximation
----------------------------------------

Consider the measure space $({{\mathbb{C}}^{{n}}}, \mathcal{B})$, where $\mathcal{B}$ is the $\sigma$-algebra of Borel sets on ${{\mathbb{C}}^{{n}}}$. Denote ${\mathcal{M}}$ the collection of finite signed measures with bounded second moment on this measure space: $${\mathcal{M}}:= \biggl\{\, \mu \text{ is a finite signed measure}\ \bigg\vert \  \int {\lvertx\rvert}^2\ {\left\lvert\mu\right\rvert}({\,\mathrm{d}}x)  < \infty \, \biggr\}\,,$$ where the positive measure ${\left\lvert\mu\right\rvert}$ is variation of measure $\mu$.

Consider a given differentiable trajectory $\mu_t\in {\mathcal{M}}$ solving $$\partial_t \mu_t = {\mathcal{A}}_t \mu_t,$$ where ${\mathcal{A}}_t: {\mathcal{M}}\rightarrow {\mathcal{M}}$ is a given time-dependent (linear) operator. In the context of Fokker-Planck equation, ${\mathcal{A}}_t$ is the adjoint operator of the infinitesimal generator of the corresponding SDE. In the context of [Lindblad]{} equation, ${\mathcal{A}}_t$ is the adjoint operator of the infinitesimal generator of the SDE-type unraveling scheme of that [Lindblad]{} equation (see [Section \[sec::discuss::unravel\]]{}).

The low-rank approximation of ${\mathcal{M}}$, denoted by ${\mathcal{M}}_r$, is a subset of ${\mathcal{M}}$, which contains all measures in ${\mathcal{M}}$ with support on a $r$-dimensional linear subspace of ${{\mathbb{C}}^{{n}}}$. Such low-rankness is used to deal with the problem of high dimensionality of ${{\mathbb{C}}^{{n}}}$. As a remark, the low-rankness we explore here is not in the sense of taking an ansatz of the measure in the space ${\mathcal{M}}$ as a linear combination of a few prescribed measures as a basis (which would be a usual Galerkin approximation in the space of measures). Instead, the low-rankness here means that the measure $\mu$ is mostly concentrated on a $r$-dimensional linear subspace of ${{\mathbb{C}}^{{n}}}$, where $r \ll n$. Intuitively, this approximation would work well for some dissipative dynamics for which the measure is contracted to some low-dimensional space as time evolves (see [Section \[sec::numerics\]]{} for numerical demonstration).

Let us characterize the structure of ${\mathcal{M}}_r$. For any $\mu_{{LR}}\in {\mathcal{M}}_r$, by definition, it is supported on a $r$-dimensional linear subspace, whose orthonormal basis is denoted by $\left\{U_1, U_2, \cdots, U_r \right\}$. Then one could define a linear mapping ${{\boldsymbol{U}}}: {\mathbb{C}}^r\rightarrow {{\mathbb{C}}^{{n}}}$ by $${{\boldsymbol{U}}}: y\in {\mathbb{C}}^r \rightarrow \sum_{j=1}^{r} U_j y_j, \qquad \text{ or in matrix form } {{\boldsymbol{U}}}= \begin{bmatrix} U_1 & U_2 & \cdots & U_r \end{bmatrix}.$$ Let us denote the $r$-dimensional Stiefel manifold on ${{\mathbb{C}}^{{n}}}$ by $${V_r({{\mathbb{C}}^{{n}}})}:=  \left\{ {{\boldsymbol{U}}}: {\mathbb{C}}^r \rightarrow {{\mathbb{C}}^{{n}}}\text{ is linear, and }\ {{\boldsymbol{U}}}^{\dagger} {{\boldsymbol{U}}}= 
 {\mathrm{Id}}_{r\times r} \right\}.$$ Then, ${\mathcal{M}}_r$ could be viewed as a collection of measures in ${\mathcal{M}}$ with support on ${\text{Ran}}({{\boldsymbol{U}}})$ with bounded second moment, where ${\text{Ran}}({{\boldsymbol{U}}})$ is the range of some linear operator ${{\boldsymbol{U}}}\in {V_r({{\mathbb{C}}^{{n}}})}$. The restriction of $\mu_{{LR}}$ on $\text{Ran}({{\boldsymbol{U}}})$ can be represented as a finite signed measure on ${\mathbb{C}}^r$, given by the pullback $$\theta(E) := \mu_{{LR}}({{\boldsymbol{U}}}E) 
= \mu_{{LR}}(x \in {\text{Ran}}( {{\boldsymbol{U}}}): x\in {{\boldsymbol{U}}}E) = \mu_{{LR}}(x \in {\text{Ran}}( {{\boldsymbol{U}}}): {{\boldsymbol{U}}}^{\dagger} x\in E)$$ for any Borel set $E \subset {\mathbb{C}}^r$. Hence, for any measurable set $F \subset \text{Ran}({{\boldsymbol{U}}})\in \mathcal{B}$, $$\mu_{{LR}}(F) = \mu_{{LR}} ({{\boldsymbol{U}}}{{\boldsymbol{U}}}^{\dagger} F) = \theta( {{\boldsymbol{U}}}^{\dagger} F).$$ Thus, there exists a one-to-one correspondence between ${\mathcal{M}}_r$ and ${V_r({{\mathbb{C}}^{{n}}})}\oplus {{\mathcal{M}}_{{\mathbb{C}}^r}}$ where ${{\mathcal{M}}_{{\mathbb{C}}^r}}$ denotes the space of finite signed measures on ${\mathbb{C}}^r$, with bounded second moment. The low-rank dynamics on ${\mathcal{M}}_r$ that we shall consider is equivalent to the dynamics of a pair $\bigl({{\boldsymbol{U}}}(t), \theta_t\bigr)\in {V_r({{\mathbb{C}}^{{n}}})}\oplus {{\mathcal{M}}_{{\mathbb{C}}^r}}$. We remark that we use in the general framework finite signed measure instead of probability measure to avoid the subtleties arising from the geometry of probability measures (due to the positivity), see for example [@Ambrosio; @Lott17]. In practice, we will guarantee that the resulting dynamics yields probability measure by imposing more constraints on the low-rank approximation, see [Section \[sec::example\]]{}.

Tangent space projection
------------------------

It is well-known that the tangent space of ${V_r({{\mathbb{C}}^{{n}}})}$ at ${{\boldsymbol{U}}}\in {V_r({{\mathbb{C}}^{{n}}})}$ is given by (see for example [@Higham96 Theorem 1.2]), $$\label{eq:tangUV}
{\mathcal{T}}_{{{\boldsymbol{U}}}} {V_r({{\mathbb{C}}^{{n}}})}= \left\{ i {{\boldsymbol{G}}}{{\boldsymbol{U}}}: {{\boldsymbol{G}}}^{\dagger} = {{\boldsymbol{G}}}\text{ is a linear operator on }{{\mathbb{C}}^{{n}}}\right\}.$$ A self-contained proof is provided in \[sec::stiefel\_tangent\] for readers’ convenience. Thus for a differentiable trajectory ${{\boldsymbol{U}}}(t)\in {V_r({{\mathbb{C}}^{{n}}})}$, we have $\frac{{\mathrm{d}}}{{\mathrm{d}}t} {{\boldsymbol{U}}}(t) = i {{\boldsymbol{G}}}(t) {{\boldsymbol{U}}}(t)$ for some Hermitian matrices ${{\boldsymbol{G}}}(t)$.

[[Consider]{}]{} any differentiable trajectory $\theta_t$ that $$\label{eq:thetaeqn}
  \partial_t \theta_t = {\mathcal{A}}_{\theta,t} \theta_t,$$ where ${\mathcal{A}}_{\theta,t}: {{\mathcal{M}}_{{\mathbb{C}}^r}}\rightarrow {{\mathcal{M}}_{{\mathbb{C}}^r}}$ is some operator. The tangent space of ${\mathcal{M}}_{r}$ at $({{\boldsymbol{U}}}(t), \theta_t)$ is fully characterized by $${\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}\oplus {\mathcal{T}}_{\theta_t} {{\mathcal{M}}_{{\mathbb{C}}^r}}.$$

It is straightforward to adopt the idea of tangent space projection ([*i.e.*]{}, [Eq. ]{}) to our situation. Consider a natural [pseudometric]{} on ${\mathcal{M}}$ $$\label{eqn::metric}
d_{{\mathsf{F}}} \left(\nu_1, \nu_2\right):= \sup_{f\in {\mathsf{F}}} {\left\lVert\int f\ {\,\mathrm{d}}\nu_1 - f\ {\,\mathrm{d}}\nu_2 \right\rVert},\qquad \nu_{1,2} \in  {\mathcal{M}}.$$ Recall that ${\mathsf{F}}$ is a collection of test functions and ${\left\lVert\cdot\right\rVert}$ is some suitable norm associated with functions in ${\mathsf{F}}$. In [Section \[sec::example\]]{}, we will choose $${\mathsf{F}}= \left\{f = xx^{\dagger}\right\}$$ be a singleton and the norm ${\left\lVert\cdot\right\rVert}$ is chosen as Hilbert-Schmidt norm. An equivalent choice is that $${\mathsf{F}}= \Bigl\{\,f(x) = {\left\langlex, {{\boldsymbol{O}}}x\right\rangle} \;\big\vert \;
\text{Hermitian matrix } {{\boldsymbol{O}}}\text{ satisfies } {\lVert{{\boldsymbol{O}}}\rVert}_{{HS}} \le
1\Bigr\},$$ and the norm is simply the absolute value. The equivalence of these two choices is proved in Lemma \[lemma::choice\_equivalent\_F\] in \[sec:equivalence\]. In fact, from the perspective of quantum mechanics, this equivalence is natural, since finding a good approximation of density matrix $\int f\ {\,\mathrm{d}}\nu \equiv {\mathbb{E}}_{\nu}\bigl[xx^{\dagger}\bigr]$ is equivalent to finding a good approximation of all observations as ${\left\langle{{\boldsymbol{O}}}\right\rangle}_{\text{avg}} := {\mathbb{E}}_{\nu} \bigl[ {\left\langlex, {{\boldsymbol{O}}}x\right\rangle}\bigr] = \operatorname{tr}\left( {\mathbb{E}}_{\nu} \bigl[ {{\boldsymbol{O}}}xx^{\dagger}\bigr]\right)
= {\left\langle{{\boldsymbol{O}}}, {\mathbb{E}}_{\nu}\bigl[xx^{\dagger}\bigr]\right\rangle}_{{HS}}$ where observable ${{\boldsymbol{O}}}$ is a Hermitian matrix.

The tangent space projection of $\partial_t \mu_t = {\mathcal{A}}_{t} \mu_t$ to the tangent space ${\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r$ is then given by $$\label{eqn::sdlr_minimization}
\begin{split}
  \partial_t \mu_{{LR},t} \equiv {\mathcal{A}}_{{LR},t} \mu_{{LR},t} := & \operatorname*{arg\,min}_{\nu \in {\mathcal{T}}_{\mu_{{LR},t}}{\mathcal{M}}_r}  d_{{\mathsf{F}}}\left(\nu,\ {\mathcal{A}}_{t} \mu_{{LR},t}\right) \\
= & \operatorname*{arg\,min}_{{\widetilde{{\mathcal{A}}}}_{{LR},t}:\ {\widetilde{{\mathcal{A}}}}_{{LR},t} \mu_{{LR},t} \in {\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r} \sup_{f\in {\mathsf{F}}} {\left\lVert\int f\ {\,\mathrm{d}}\bigl({\widetilde{{\mathcal{A}}}}_{{LR},t}\mu_{{LR}, t}\bigr)  - f\ {\,\mathrm{d}}\bigl({\mathcal{A}}_t \mu_{{LR},t}\bigr) \right\rVert}. \\
\end{split}$$ In the second line, the minimization problem is reformulated from finding tangent vector $\nu$ to finding differential operators ${\widetilde{{\mathcal{A}}}}_{{LR},t}$. Though the notation is slightly abused, the variational problem above should still be clear.

Equivalently, using the adjoint operators ${\widetilde{{\mathcal{A}}}}_{{LR},t}^{*}$ and ${\mathcal{A}}_{t}^{*}$, it can be written as $$\label{eqn::sdlr_minimization_equiv}
\begin{split}
 \partial_t \mu_{{LR},t} = & \operatorname*{arg\,min}_{{\widetilde{{\mathcal{A}}}}_{{LR},t}^{*}:\ {\widetilde{{\mathcal{A}}}}_{{LR},t} \mu_{{LR},t} \in {\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r} \sup_{f\in {\mathsf{F}}} {\left\lVert\int \bigl({\widetilde{{\mathcal{A}}}}_{{LR},t}^{*} f\bigr)\ {\,\mathrm{d}}\mu_{{LR}, t} - \bigl({\mathcal{A}}_t^{*} f\bigr)\ {\,\mathrm{d}}\mu_{{LR},t}\right\rVert} \\
=& \operatorname*{arg\,min}_{{\widetilde{{\mathcal{A}}}}_{{LR},t}^{*}:\ {\widetilde{{\mathcal{A}}}}_{{LR},t} \mu_{{LR},t} \in {\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r} \sup_{f\in {\mathsf{F}}} {\left\lVert{\mathbb{E}}_{\mu_{{LR},t}} \left[{\widetilde{{\mathcal{A}}}}_{{LR},t}^{*} f - {\mathcal{A}}_t^{*} f \right] \right\rVert}. \\
\end{split}$$ In [Section \[sec::example\]]{}, we shall parametrize the infinitesimal generator ${\widetilde{{\mathcal{A}}}}_{{LR},t}^{*}$ by some functions to simplify the minimization problem and also to avoid the vagueness of generic infinitesimal generator for general stochastic processes.

SDLR method for SDEs with $N$ driving Brownian motions and test function $f(x) = x x^{\dagger}$  {#sec::example}
================================================================================================

Based on the framework of SDLR, we may explore various low-rank dynamics. In this section, $\mu_t$ is considered to be the probability measure of $X_t$ given in [Eq. ]{} and thus $\partial_t \mu_t = {\mathcal{A}}_t \mu_t $ corresponds to the Fokker-Planck equation of SDE $X_t$ in [Eq. ]{}. The test function space ${\mathsf{F}}$ is taken to be a singleton with the only element $f(x) = x x^{\dagger}$. We will assume these in the sequel without explicit mentioning.

To apply the method developed in last section, we impose some further restrictions on ${\mathcal{T}}_{\mu_{{LR},t}} {\mathcal{M}}_r$ and thus on the choice of low-rank dynamics. We will make some further comments on these restrictions after we derive the resulting low-rank dynamics by applying [SDLR]{} method.

1.  In the tangent space of ${V_r({{\mathbb{C}}^{{n}}})}$, for Hermitian matrix ${{\boldsymbol{G}}}(t)$ in [Eq. ]{}, we [[consider only those ${{\boldsymbol{G}}}(t)$ such that]{}]{} $$\label{eqn::assmp_g}
      {{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\text{h.c.}},$$ where projection operator ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} := {{\boldsymbol{U}}}(t) {{\boldsymbol{U}}}^{\dagger}(t)$ and orthogonal projection operator ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} := {\mathrm{Id}}- {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}$.

2.  In the tangent space of ${{\mathcal{M}}_{{\mathbb{C}}^r}}$, [[consider only those ${\mathcal{A}}_{\theta,t}$ corresponding]{}]{} to the Fokker-Planck equation of some SDE on ${\mathbb{C}}^r$ (cf. [Eq. ]{}). This basically means that we parametrize the infinitesimal generator ${\mathcal{A}}_{{LR},t}^{*}$ (and hence ${\mathcal{A}}_{{LR},t}$) by a collection of functions. The exact form will be given below in Lemma \[lemma::generator\]. Because the space of all possible infinitesimal generators is opaque and quite large, hence, choosing infinitesimal generators of a particular form is necessary in practice.

Following the second [[constraint]{}]{} above, let us denote $Y_t(\omega)$ the corresponding SDE on ${\mathbb{C}}^r$ whose infinitesimal generator is ${\mathcal{A}}_{\theta,t}$. Consider the following family of stochastic dynamics with ${{\boldsymbol{G}}}$, $A$ and $B$’s to be chosen $$\label{eqn::sdeY}
\begin{split}
  \dot{{{\boldsymbol{U}}}}(t) &=  i {{\boldsymbol{G}}}(t) {{\boldsymbol{U}}}(t),\qquad 
  \text{with } {{\boldsymbol{G}}}^{\dagger}(t) = {{\boldsymbol{G}}}(t),\  {{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\text{h.c.}}\\
  {\,\mathrm{d}}Y &= A(Y,t) {\,\mathrm{d}}t + \sum_{j=1}^{M} B_j(Y,t) {\,\mathrm{d}}W_j. \\
  \end{split}$$ Define $X_{{LR}, t}(\omega) := {{\boldsymbol{U}}}(t) Y_t(\omega)$ and denote $\mu_{{LR}, t}$ the probability measure induced by random variable $X_{{LR},t}(\omega)$. It is straightforward to check that $X_{{LR},t}(\omega)$ satisfies the SDE $${\,\mathrm{d}}X_{{LR}} = a_{{LR}}\left(X_{{LR}}, t\right){\,\mathrm{d}}t + \sum_{j=1}^{M}b_{{LR},j} \left(X_{{LR}}, t\right){\,\mathrm{d}}W_j,
  \label{eqn::sde_lowrank}$$ where, for $1\le j\le M$, $$\begin{cases}
    \begin{aligned}
      a_{{LR}}(X_{{LR},t}, t) & =  
      i {{\boldsymbol{G}}}(t) X_{{LR},t} + {{\boldsymbol{U}}}(t) A({{\boldsymbol{U}}}^{\dagger}(t) X_{{LR},t}, t),
    \end{aligned} \\
    b_{{LR}, j}(X_{{LR},t}, t)= {{\boldsymbol{U}}}(t) B_j({{\boldsymbol{U}}}^{\dagger}(t) X_{{LR},t}, t).
  \end{cases}$$ Note that $M$ does not have to be the same as $N$ and ${\mathrm{d}}W_j$ in [Eq. ]{} are not necessarily the same Brownian motions in [Eq. ]{}; we use ${\mathrm{d}}W_j$ for both to save notation. Since we are interested in the error in the weak sense, how the randomness is achieved does not matter; what is important is the infinitesimal generator which does not depend on the particular realization of the Brownian motion.

\[lemma::generator\] The adjoint infinitesimal generator, acting on $f(x) = x x^{\dagger}$, for SDE $X_t$ is $$\begin{split}
\bigl({\mathcal{A}}^{*}_t f\bigr)(x) &=
x a^{\dagger}(x,t) + a(x,t) x^{\dagger} +  \sum_{j=1}^{N} b_j(x,t) b_j^{\dagger}(x,t). \\
\end{split}$$ Hence for low-rank dynamics $\mu_{{LR},t}$, ${\mathcal{A}}_{{LR},t}^{*}$ is a family of generators parametrized by ${{\boldsymbol{G}}}(t)$, $A(y)$ and $B_j(y)$ ($1\le j\le M$), with the following form $$\begin{split}
 \bigl({\mathcal{A}}_{{LR},t}^{*} f\bigr)(x) &
= x \left( i {{\boldsymbol{G}}}(t) x + {{\boldsymbol{U}}}(t) A( {{\boldsymbol{U}}}^{\dagger}(t) x, t)\right)^{\dagger} + \left( i {{\boldsymbol{G}}}(t) x + {{\boldsymbol{U}}}(t) A( {{\boldsymbol{U}}}^{\dagger}(t) x, t)\right) x^{\dagger} \\
& \qquad + \sum_{j=1}^{M} {{\boldsymbol{U}}}(t) B_j( {{\boldsymbol{U}}}^{\dagger}(t) x, t) B_j^{\dagger}({{\boldsymbol{U}}}^{\dagger}(t) x, t) {{\boldsymbol{U}}}^{\dagger}(t). \\
\end{split}$$

Note that ${\mathcal{A}}_{{LR},t}^{*}$ depends on the current ${{\boldsymbol{U}}}(t)$, which is natural. The proof is straightforward by applying [It[ô]{}]{} formula to $f(x) = xx^{\dagger}$.

We are now ready to apply [Eq. ]{} and to find the time-evolution of the optimal low-rank dynamics $\mu_{{LR},t}$.

\[thm::fx=xx\] Assume that $\int f(y)\ \theta_t({\,\mathrm{d}}y) \equiv {\mathbb{E}}_{Y_t\sim \theta_t} \bigl[Y_t Y^{\dagger}_t\bigr]$ is an invertible $r\times r$ matrix for any $t\in {[0,T]}$, then

(i) The following choices of ${{\boldsymbol{U}}}(t)$ and $\theta_t$ give an optimal low-rank dynamics $$\label{eqn::low_rank_xx}
    \begin{cases} 
    M  = N\\
    A(y, t) =  {{\boldsymbol{U}}}^{\dagger}(t) a({{\boldsymbol{U}}}(t) y, t)  \\
    B_j(y, t) = {{\boldsymbol{U}}}^{\dagger}(t) b_j({{\boldsymbol{U}}}(t) y, t) \\
    {{\boldsymbol{G}}}(t) = (-i) {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \biggl( {\mathbb{E}}_{\theta_t} \bigl[a({{\boldsymbol{U}}}(t) y, t) ({{\boldsymbol{U}}}(t) y)^{\dagger}\bigr] \\ \hspace{10em} + \sum_{j=1}^{N} {\mathbb{E}}_{\theta_t}\bigl[b_j({{\boldsymbol{U}}}(t) y, t) b_j^{\dagger}({{\boldsymbol{U}}}(t) y, t) \bigr] \biggr) {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}[yy^{\dagger}]^{-1} {{\boldsymbol{U}}}^{\dagger}(t) + {\text{h.c.}}\end{cases}$$

(ii) The solution is optimal in the sense that, for any given ${{\boldsymbol{U}}}(t)$ and $\theta_t$, $$\begin{aligned}
          \min_{{\widetilde{{\mathcal{A}}}}_{{LR},t}}\; d_{{\mathsf{F}}}\bigl({\widetilde{{\mathcal{A}}}}_{{LR},t} \mu_{{LR},t},\
          {{\mathcal{A}}}_{t} \mu_{{LR},t} \bigr) & =
          d_{{\mathsf{F}}}\bigl( {{\mathcal{A}}}_{{LR},t}
          \mu_{{LR},t},\ {\mathcal{A}}_{t} \mu_{{LR},t} \bigr)\\
          & = {\left\lVert\sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t}
            \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr]
            {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\right\rVert}_{{HS}},
        \end{aligned}$$ where the adjoint of ${{\mathcal{A}}}_{{LR},t}$ ([*i.e.*]{}, ${{\mathcal{A}}}^{*}_{{LR},t}$) is given by [Eq. ]{}, ${\widetilde{{\mathcal{A}}}}_{{LR},t}$ is any operator acting on probability measures whose corresponding infinitesimal generator ${\widetilde{{\mathcal{A}}}}_{{LR},t}^{*}$ is parametrized by $A$, $B_j$ and ${{\boldsymbol{G}}}$ as in Lemma \[lemma::generator\].

(iii) Time-evolution equation of ${{\boldsymbol{U}}}(t)$ and $\theta_t$ is given by the following system: $$\left\{
          \label{eqn::low_rank_xx_v2}
          \begin{split}
            {\,\mathrm{d}}Y &= {{\boldsymbol{U}}}^{\dagger}(t) a({{\boldsymbol{U}}}(t) Y, t)\ {\,\mathrm{d}}t 
            + \sum_{j=1}^{N} {{\boldsymbol{U}}}^{\dagger}(t) b_j({{\boldsymbol{U}}}(t) Y, t)\ {\,\mathrm{d}}W_j\,,\\
            \frac{{\,\mathrm{d}}{{\boldsymbol{U}}}}{{\,\mathrm{d}}t} &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \left({\mathbb{E}}\bigl[a({{\boldsymbol{U}}}(t) Y, t) Y^{\dagger}\bigr] + \sum_{j=1}^{N} {\mathbb{E}}\bigl[b_j({{\boldsymbol{U}}}(t) Y, t) b_j^{\dagger}({{\boldsymbol{U}}}(t) Y, t)\bigr] {{\boldsymbol{U}}}(t) \right) {\mathbb{E}}\bigl[Y Y^{\dagger}\bigr]^{-1}\,. \\
          \end{split}\right.$$ If rank $r = n$, then ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} = 0$ and ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = {\mathrm{Id}}$. Consequently, $\frac{{\,\mathrm{d}}{{\boldsymbol{U}}}}{{\,\mathrm{d}}t} = 0$ and the infinitesimal generator for SDE of $X_{{LR},t}$ is the same as that for $X_t$. That means, when full rank is used, the original SDE is recovered; equivalently the Fokker-Planck equation for $\mu_t$ is recovered.

Since the time evolution of ${{\boldsymbol{U}}}(t)$ and $\theta_t$ can be fully recovered from solving ODE-SDE coupled system of $Y$ and ${{\boldsymbol{U}}}$ (with multiple replica of $Y$), we shall refer [Eq. ]{} as the resulting dynamical low-rank approximation as well when no confusion arises. After all, we will not solve equation $\partial_t \mu_{{LR},t} = {\mathcal{A}}_{{LR},t} \mu_{{LR},t}$ directly; instead, we shall use Monte Carlo method, [*i.e.*]{}, solving [Eq. ]{} to estimate $\mu_{{LR},t}$ by the empirical measure.

The main idea of the proof is to use stationary conditions with respect to ${{\boldsymbol{G}}}$, $A$ and $B_j$ to derive the low-rank dynamics. Recall that ${\mathcal{P}}_{{{\boldsymbol{U}}}} := {{\boldsymbol{U}}}{{\boldsymbol{U}}}^{\dagger}$ and ${\mathcal{Q}}_{{{\boldsymbol{U}}}} := {\mathrm{Id}}- {{\boldsymbol{U}}}{{\boldsymbol{U}}}^{\dagger}$. Fix time $t$, by [Eq. ]{} and , we know $$\begin{split}
 d_{{\mathsf{F}}}\bigl({\mathcal{A}}_{{LR},t} \mu_{{LR},t},\ {\mathcal{A}}_{t} \mu_{{LR},t} \bigr) 
= {\left\lVert{\mathbb{E}}_{\mu_{{LR},t}} \bigl[{\mathcal{A}}_{{LR},t}^{*} f - {\mathcal{A}}_{t}^{*} f \bigr] \right\rVert}_{{HS}} 
=: {\left\lVert{{\boldsymbol{C}}}\right\rVert}_{{HS}}, \\
\end{split}$$ where $$\begin{split}
 {{\boldsymbol{C}}}&:={\mathbb{E}}_{\mu_{{LR},t}}\left[x (i {{\boldsymbol{G}}}(t) x + {{\boldsymbol{U}}}(t) A({{\boldsymbol{U}}}^{\dagger}(t) x, t) - a(x,t))^{\dagger} + {\text{h.c.}}\right] \\
 & \qquad + {\mathbb{E}}_{\mu_{{LR},t}}\left[ \sum_{j=1}^{M} 
 {{\boldsymbol{U}}}(t) B_j( {{\boldsymbol{U}}}^{\dagger}(t) x, t)   
 B_j^{\dagger}({{\boldsymbol{U}}}^{\dagger}(t) x, t) 
 {{\boldsymbol{U}}}^{\dagger}(t)  
 - \sum_{j=1}^{N} b_j(x,t) b_j^{\dagger}(x,t)\right].\\
 \end{split}$$ To minimize $d_{{\mathsf{F}}}\bigl({\mathcal{A}}_{{LR},t} \mu_{{LR},t},\ {\mathcal{A}}_{t} \mu_{{LR},t} \bigr)$, it is equivalent to minimize ${\left\lVert{{\boldsymbol{C}}}\right\rVert}_{{HS}}^2 = {\left\langle{{\boldsymbol{C}}}, {{\boldsymbol{C}}}\right\rangle}_{{HS}}$. The first order stationary conditions of ${\left\langle{{\boldsymbol{C}}}, {{\boldsymbol{C}}}\right\rangle}_{{HS}}$ with respect to $A$, $B_j$ ($1\le j \le M$) and ${{\boldsymbol{G}}}(t)$ gives $$\label{eqn::statn_codn}
\left\{
\begin{split}
0 &= {\left\langle{{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{C}}}{{\boldsymbol{U}}}(t), \ {\mathbb{E}}_{\theta_t}\bigl[(\delta A) y^{\dagger} + y (\delta A)^{\dagger}\bigr]\right\rangle}_{{HS}} \\
0 &= {\left\langle{{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{C}}}{{\boldsymbol{U}}}(t), \ {\mathbb{E}}_{\theta_t}\bigl[(\delta B_j) B_j^{\dagger} + B_j (\delta B_j)^{\dagger}\bigr]\right\rangle}_{{HS}} \\
0 &= \operatorname{tr}\left( {\left[{{\boldsymbol{C}}},  {\mathbb{E}}_{\mu_{{LR},t}}[ x x^{\dagger}]\right]_{-}} \delta {{\boldsymbol{G}}}\right),
\end{split}\right.$$ where $\delta A$ and $\delta B_j$ are perturbations of functions $A$ and $B_j$ respectively and $\delta {{\boldsymbol{G}}}$ is a Hermitian matrix as perturbation of ${{\boldsymbol{G}}}$ and ${\left[\cdot, \cdot\right]_{-}}$ is the commutator.

Consider the third condition in [Eq. ]{}. Let us complete the basis of ${{\mathbb{C}}^{{n}}}$ by extending $\text{Ran}({{\boldsymbol{U}}}(t))$, denoted by $\{U_1, U_2, \cdots, U_r, U_{r+1}, \cdots \cdots\}$. Since $\delta {{\boldsymbol{G}}}$ is arbitrary among all possible perturbations, consider the special choice $\delta {{\boldsymbol{G}}}= \lambda U_j {\left\langleU_k, \cdot\right\rangle} + {\text{h.c.}}$ where $j \le r$ and $k > r$. Note that when both $j, k \le r$ or $j, k > r$, by [[our restriction in [Eq. ]{}]{}]{}, ${\left\langleU_j, {{\boldsymbol{G}}}(t) U_k\right\rangle} = 0$. Hence as the perturbation of ${{\boldsymbol{G}}}$, $\delta {{\boldsymbol{G}}}$ should preserve this property. That means, we cannot choose $\delta {{\boldsymbol{G}}}$ with nonzero entries for $j, k \le r$ nor $j, k > r$. Denote ${{\boldsymbol{D}}}\equiv \bigl[{{\boldsymbol{C}}}, {\mathbb{E}}_{\mu_{{LR},t}}[ x x^{\dagger}]\bigr]_{-}$, which is anti-Hermitian, that is, ${{\boldsymbol{D}}}^{\dagger} = - {{\boldsymbol{D}}}$. Plugging the expression of $\delta {{\boldsymbol{G}}}$ into the third condition in [Eq. ]{}, we could easily compute that $$0 = \lambda {\left\langleU_k, {{\boldsymbol{D}}}U_j\right\rangle} + \lambda^{*} {\left\langleU_j, {{\boldsymbol{D}}}U_k\right\rangle} = \lambda {\left\langleU_k, {{\boldsymbol{D}}}U_j\right\rangle} + \lambda^{*} {\left\langleU_k, {{\boldsymbol{D}}}^{\dagger} U_j\right\rangle}^{*} = \lambda {\left\langleU_k, {{\boldsymbol{D}}}U_{j}\right\rangle}  - c.c.$$ Hence $\text{Im}\left(\lambda {\left\langleU_k, {{\boldsymbol{D}}}U_j\right\rangle}\right) = 0$. Since $\lambda\in {\mathbb{C}}$ is arbitrary, ${\left\langleU_k,  {{\boldsymbol{D}}}U_j\right\rangle} = 0$ for all $j\le r$ and $k > r$. That is to say, ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{D}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = 0$. By plugging the expression of ${{\boldsymbol{D}}}\equiv \bigl[{{\boldsymbol{C}}}, {\mathbb{E}}_{\mu_{{LR},t}}[ x x^{\dagger}]\bigr]_{-}$, we could compute that $$\begin{split} 
0 
&= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t} \bigl[yy^{\dagger}\bigr] {{\boldsymbol{U}}}^{\dagger}(t)\\
&= {\mathbb{E}}_{\theta_t} \Biggl[ i {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {{\boldsymbol{U}}}(t) y y^{\dagger} 
- {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} a({{\boldsymbol{U}}}(t)y,t) y^{\dagger}\\
& \hspace{4em} -{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\sum_{j=1}^{N} b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t) {{\boldsymbol{U}}}(t) \Biggr] {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr] {{\boldsymbol{U}}}^{\dagger}(t) \\
&= i {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t} \bigl[y y^{\dagger}\bigr]^2 {{\boldsymbol{U}}}^{\dagger}(t) 
- {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[a({{\boldsymbol{U}}}(t)y,t) y^{\dagger}\bigr] {\mathbb{E}}_{\theta_t}\bigl[y y^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t) \\
&- {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \sum_{j=1}^{N} {\mathbb{E}}_{\theta_t}\bigl[ b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t) \bigr] {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}[ y y^{\dagger}] {{\boldsymbol{U}}}^{\dagger}(t).
 \end{split}$$ Multiply both sides by ${{\boldsymbol{U}}}(t){\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]^{-2}$ on the right, the last equation yields $$i {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {{\boldsymbol{U}}}(t)={\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \left( {\mathbb{E}}_{\theta_t} \bigl[ a({{\boldsymbol{U}}}(t)y,t)  y^{\dagger}\bigr] + \sum_{j=1}^{N} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t) \bigr] {{\boldsymbol{U}}}(t)\right) {\mathbb{E}}_{\theta_t}[y y^{\dagger}]^{-1}.$$ Then multiply ${{\boldsymbol{U}}}^{\dagger}(t)$ on the right, and then divide both side by $i$, one could obtain $$\begin{split}
 {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} 
 =&  (-i) {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \left( {\mathbb{E}}_{\theta_t} \bigl[a({{\boldsymbol{U}}}(t)y,t) ({{\boldsymbol{U}}}(t) y)^{\dagger}\bigr] + \sum_{j=1}^{N} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t) \bigr] \right) {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}[y y^{\dagger}]^{-1} {{\boldsymbol{U}}}^{\dagger}(t).\\
 \end{split}$$ Thus we have already obtained the expression of ${{\boldsymbol{G}}}(t) \equiv {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\text{h.c.}}$ (cf. [Eq. ]{}). [ As a remark, up to here, we have not yet used any information nor assumption about $A$ and $B_j$. ]{}

Then we plug the expression of ${{\boldsymbol{G}}}(t)$ into ${{\boldsymbol{C}}}$, $$\begin{split}
{{\boldsymbol{C}}}&= {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[ \bigl({{\boldsymbol{U}}}(t) A(y,t) - a({{\boldsymbol{U}}}(t)y,t)\bigr)y^{\dagger} {{\boldsymbol{U}}}^{\dagger}\bigr] + {\text{h.c.}}\\
& \qquad + \sum_{j=1}^{M} {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t} \bigl[B_j(y,t) B_j^{\dagger}(y,t)\bigr] {{\boldsymbol{U}}}^{\dagger}(t)  \\
& \qquad - \sum_{j=1}^{N} \left( {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} 
 + {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\right). \\
\end{split}$$ One could observe that ${{\boldsymbol{C}}}= {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}$ with $$\left\{
\begin{split}
{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} &= {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[ \bigl({{\boldsymbol{U}}}(t) A(y,t) - a({{\boldsymbol{U}}}(t)y,t)\bigr)y^{\dagger} {{\boldsymbol{U}}}^{\dagger}\bigr] + {\text{h.c.}}+ \sum_{j=1}^{M} {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t} \bigl[B_j(y,t) B_j^{\dagger}(y,t)\bigr] {{\boldsymbol{U}}}^{\dagger}(t) \\
&  - \sum_{j=1}^{N}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}, \\
{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} &= - \sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}.
\end{split}\right.$$ Notice that ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}$ does not depend on any parameter we choose, that is, it is independent of $A$, $B_j$ and ${{\boldsymbol{G}}}(t)$. Recall that we would like to minimize $${\left\langle{{\boldsymbol{C}}}, {{\boldsymbol{C}}}\right\rangle}_{{HS}} 
= {\left\langle{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}, {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}\right\rangle}_{{HS}} 
+ {\langle{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}, {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\rangle}_{{HS}}.$$ The second term is non-negative and we cannot minimize it further. As for the first term, by choosing $$\left\{
\begin{split}
M &= N, \\
A(y,t) &=  {{\boldsymbol{U}}}^{\dagger}(t) a({{\boldsymbol{U}}}(t) y, t), \\
B_j(y,t) &= {{\boldsymbol{U}}}^{\dagger}(t) b_j({{\boldsymbol{U}}}(t)y,t), \\
\end{split}\right.$$ one could easily verify that ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = 0$. With such choice, $${\left\lVert{{\boldsymbol{C}}}\right\rVert}_{{HS}} = {\left\lVert\sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\right\rVert}_{{HS}}.$$ Hence the above choice must be optimal (although it does not imply uniqueness). However, the choice of ${{\boldsymbol{G}}}$ is indeed unique under [[our restriction]{}]{}. One could straightforwardly verify that this solution satisfies the first two equations in the first order stationary condition ([*i.e.*]{}, [Eq. ]{}); that is to say, the above choice yields ${{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{C}}}{{\boldsymbol{U}}}(t) = 0$. Though we don’t have to use first two parts in [Eq. ]{} to derive low-rank dynamics, it is still nice to observe the consistency.

Thus we have proved the first and second part of this theorem. The third part follows easily from the first part.

Let us come back to the two restrictions we made at the beginning of this section.

-   The reason to impose the condition ${{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\text{h.c.}}$ is to remove the redundant degree of freedom: It is easy to observe that ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}$ is redundant and does not play any role in $\dot{{{\boldsymbol{U}}}}(t)$. Thus, we might as well let ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \equiv 0$. Further if we assume ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t){\mathcal{P}}_{{{\boldsymbol{U}}}(t)} \equiv 0$, then it directly implies that ${{\boldsymbol{U}}}^{\dagger}(t) \dot{{{\boldsymbol{U}}}}(t) = 0$, which is used as orthogonal constraint in [@Lubich07]. Conversely, if ${{\boldsymbol{U}}}^{\dagger}(t) \dot{{{\boldsymbol{U}}}}(t) = 0$, then ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = 0$, which shows that such constraint is similar to the constraint ${{\boldsymbol{U}}}^{\dagger}(t) \dot{{{\boldsymbol{U}}}}(t) = 0$ for matrix ODEs.

    Another reason comes from the above proof. With such constraint, the first order stationary condition with respect to $\delta {{\boldsymbol{G}}}$ yields a unique expression for ${{\boldsymbol{G}}}(t)$. This indicates that redundancy has been removed via the above constraint.

-   The main reason for imposing constraint in tangent space of ${{\mathcal{M}}_{{\mathbb{C}}^r}}$ is that Fokker-Planck type generator automatically helps to preserve the positivity of measure. In fact, the whole tangent space is rather big and too opaque to handle, since it might involve generators for other stochastic processes, [*e.g.*]{}, jump processes. The mixture of jump process and diffusion makes it more challenging to derive a simple low-rank dynamics; which could be an interesting future research direction.

It is a good place to compare our approach with the dynamical orthogonal (DO) method [@Sapsis09]. Using the current notations, the ansatz in DO method is taken to be $$X_{{LR},t}(\omega) = \bar{X}_t + {{\boldsymbol{U}}}(t) Y_t(\omega),
\label{eqn::do_assmp}$$ where $\bar{X}_t$ is deterministic and for all $t$, ${\mathbb{E}}[Y_t(\omega)]= 0$, ${{\boldsymbol{U}}}^{\dagger}(t) \dot{{{\boldsymbol{U}}}}(t) = 0$, ${{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{U}}}(t) = {\mathrm{Id}}$. By re-deriving the low-rank dynamics following the proof in that paper, one could obtain that $$\label{eqn::do}
\begin{split}
  \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}\bar{X} &= {\mathbb{E}}\bigl[a(X_{{LR},t},t)\bigr] \\
  \dot{{{\boldsymbol{U}}}} &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}\bigl[a(X_{{LR},t},t) Y^{\dagger}\bigr] {\mathbb{E}}\bigl[Y Y^{\dagger}\bigr]^{-1}\\
  {\,\mathrm{d}}Y &= {{\boldsymbol{U}}}^{\dagger}(t) \left(a(X_{{LR},t},t) - {\mathbb{E}}\bigl[a(X_{{LR},t},t)\bigr]\right) {\,\mathrm{d}}t + \sum_{j=1}^{N} {{\boldsymbol{U}}}^{\dagger}(t) b_j(X_{{LR},t},t){\,\mathrm{d}}W_j.\\
\end{split}$$ Compared with [Eq. ]{}, the expressions of time-evolution of $Y$ and ${{\boldsymbol{U}}}(t)$ are almost the same except: (1) the extra term ${{\boldsymbol{U}}}^{\dagger}(t) {\mathbb{E}}[a(X_{{LR},t},t)]$ in DO method, due to the zero-mean constraint in $Y$; and (2) the non-trivial difference that in [Eq. ]{}, we have a term $$\sum_{j=1}^{N} {\mathbb{E}}[b_j(X_{{LR},t},t) b_j^{\dagger}(X_{{LR},t},t)]
{{\boldsymbol{U}}}(t),$$ which could be understood as the [It[ô]{}]{} correction term due to the second moment. Some numerical experiments comparing [SDLR]{} method and DO method will be presented in [Section \[sec::numerics\]]{} for high-dimensional geometric Brownian motion and stochastic Burgers’ equation.

Error analysis of the stochastic dynamical low-rank approximation {#sec::error}
=================================================================

In this section, we provide some error analysis for the low-rank dynamics that we derive in [Section \[sec::example\]]{}. [Theorem \[thm::recover\]]{} indicates that the low-rank dynamics in [Eq. ]{} (or [Eq. ]{}) is optimal under our ansatz, in the sense that if $X_t$ is itself a SDE, whose range is supported on a rank $r$ subspace, then $X_{{LR},t} = X_t$ by choosing rank-$r$ low-rank dynamics [[with some additional assumptions.]{}]{} In other words, the low-rank dynamics we derive is consistent. We also prove an inequality to bound error propagation in [Theorem \[thm::bound\]]{} based on [Gr[ö]{}nwall’s inequality]{}.

Consistency: recovering low-rank dynamics {#subsec::recover}
-----------------------------------------

\[lemma::vanish\_diffusion\] Suppose $\mu_t$ is the measure induced by $X_t$ which solves a SDE of the form as in [Eq. ]{} $${\,\mathrm{d}}X_t = a(X_t, t) {\,\mathrm{d}}t + \sum_{j=1}^N b_j(X_t, t) {\,\mathrm{d}}W_j,$$ with initial condition $X_0 \sim \mu_0$. Assume that $\mu_t$ is supported on a $r$-dimensional linear subspace, whose basis forms a linear operator ${{\boldsymbol{V}}}(t):{\mathbb{C}}^r\rightarrow {{\mathbb{C}}^{{n}}}$ with ${{\boldsymbol{V}}}(t)^{\dagger} {{\boldsymbol{V}}}(t) = {\mathrm{Id}}_{r\times r}$. Then, for $x\in \operatorname{Ran}({{\boldsymbol{V}}}(t))$, the coefficients of the SDE satisfy $$a(x,t) = \left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} \right) x + {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} a(x,t), \qquad \text{and} \qquad b_j(x,t) = {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} b_j(x,t), \qquad \forall j,$$ and hence $$\sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{V}}}(t)}  {\mathbb{E}}_{\mu_t} \bigl[ b_j(x,t) b_j^{\dagger}(x,t) \bigr] {\mathcal{Q}}_{{{\boldsymbol{V}}}(t)} = 0.$$

Since the range of $X_t$ is $\operatorname{Ran}({{\boldsymbol{V}}}(t))$ by assumption, ${{\boldsymbol{V}}}(t) {{\boldsymbol{V}}}(t)^{\dagger} X_t = X_t$. Take derivative for both sides (note that ${{\boldsymbol{V}}}(t)$ is deterministic) $$\label{eqn::recover_lowrank}
    \begin{split}
    {\,\mathrm{d}}X_t &= {\,\mathrm{d}}\bigl({\mathcal{P}}_{{{\boldsymbol{V}}}(t)} X_t\bigr) 
    = \left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{\mathcal{P}}_{{{\boldsymbol{V}}}(t)}\right) X_t {\,\mathrm{d}}t + {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} {\,\mathrm{d}}X_t \\
    &=  \left(\left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{\mathcal{P}}_{{{\boldsymbol{V}}}(t)}\right) X_t + {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} a(X_t, t)\right) {\,\mathrm{d}}t + \sum_{j=1}^{N} {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} b_j(X_t, t) {\,\mathrm{d}}W_j. \\
    \end{split}$$ Matching it with the SDE of $X_t$, we obtain $a(X_t, t) = \left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{\mathcal{P}}_{{{\boldsymbol{V}}}(t)}\right) X_t + {\mathcal{P}}_{{{\boldsymbol{V}}}(t)}
    a(X_t, t)$ and $b_j(X_t, t) = {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} b_j(X_t, t)$. The conclusion of the Lemma follows.

\[thm::recover\] Besides the same assumptions of Lemma \[lemma::vanish\_diffusion\], we further assume that $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\boldsymbol{V}}}(t) = i {{\boldsymbol{F}}}(t) {{\boldsymbol{V}}}(t)$ with some Hermitian matrix ${{\boldsymbol{F}}}(t)$ which satisfies $${{\boldsymbol{F}}}(t) = {\mathcal{Q}}_{{{\boldsymbol{V}}}(t)} {{\boldsymbol{F}}}(t) {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} + {\text{h.c.}}$$ Then if $\mu_0 = \mu_{{LR},0}$ and ${{\boldsymbol{V}}}(0) = {{\boldsymbol{U}}}(0)$, we have $$\mu_t = \mu_{{LR},t},\quad \text{and} \quad {{\boldsymbol{V}}}(t) = {{\boldsymbol{U}}}(t)$$ for all $t\in{[0,T]}$ as long as the low-rank dynamics in [Eq. ]{} exists (in particular, ${\mathbb{E}}_{\theta_t}[yy^{\dagger}]$ remains invertible).

Fix time $t$ and assume $X_t = X_{{LR},t}$ in distribution and ${{\boldsymbol{V}}}(t) = {{\boldsymbol{U}}}(t)$. Next, we shall show that $X_{{LR},t}$ and $ X_t$ satisfy the same SDE locally and $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{V}}}(t) = \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{U}}}(t)$.

By [Theorem \[thm::fx=xx\]]{}, we could straightforwardly compute that $$\begin{split}
    {\,\mathrm{d}}X_{{LR},t} &= {\,\mathrm{d}}\bigl( {{\boldsymbol{U}}}(t) Y_t \bigr) = \left( \frac{{\,\mathrm{d}}{{\boldsymbol{U}}}(t)}{{\,\mathrm{d}}t}  Y_t \right) {\,\mathrm{d}}t + {{\boldsymbol{U}}}(t) {\,\mathrm{d}}Y_t  \\
    &= \frac{{\,\mathrm{d}}{{\boldsymbol{U}}}(t)}{{\,\mathrm{d}}t}  {{\boldsymbol{U}}}^{\dagger}(t) X_{{LR},t} {\,\mathrm{d}}t + {{\boldsymbol{U}}}(t) \left({{\boldsymbol{U}}}^{\dagger}(t) a(X_{{LR},t},t)\ {\,\mathrm{d}}t + \sum_{j=1}^{N} {{\boldsymbol{U}}}^{\dagger}(t) b_j(X_{{LR},t}, t) {\,\mathrm{d}}W_j \right) \\
    &= \left(\frac{{\,\mathrm{d}}{{\boldsymbol{U}}}(t)}{{\,\mathrm{d}}t}  {{\boldsymbol{U}}}^{\dagger}(t) X_{{LR},t} + {{\boldsymbol{U}}}(t) {{\boldsymbol{U}}}^{\dagger}(t) a(X_{{LR},t},t) \right) {\,\mathrm{d}}t + \sum_{j=1}^{N} {{\boldsymbol{U}}}(t) {{\boldsymbol{U}}}^{\dagger}(t) b_j(X_{{LR},t}, t) {\,\mathrm{d}}W_j \\
    &= \left( \left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}\right) X_{{LR},t} + {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} a(X_{{LR},t}, t) \right) {\,\mathrm{d}}t + \sum_{j=1}^{N} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} b_j(X_{{LR},t}, t) {\,\mathrm{d}}W_j .\\
    \end{split}$$ In the last step, we have used the fact that ${{\boldsymbol{U}}}(t) \frac{{\,\mathrm{d}}{{\boldsymbol{U}}}^{\dagger}(t)}{{\,\mathrm{d}}t} X_{{LR},t} = 0$ due to [[the restriction]{}]{} on Hermitian matrix ${{\boldsymbol{G}}}(t)$ (cf. [Eq. ]{}). Comparing with [Eq. ]{}, we obtain that the SDEs of $X_t$ and $X_{{LR}, t}$ coincide since ${{\boldsymbol{V}}}(t) = {{\boldsymbol{U}}}(t)$, as long as $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{\mathcal{P}}_{{{\boldsymbol{V}}}(t)} = \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}$. That means, we still need to verify $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{V}}}(t) = \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{U}}}(t)$.

Moreover, by [Eq. ]{}, $$\begin{split}
    \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{U}}}(t) &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} \left({\mathbb{E}}\bigl[a(X_{{LR}},t) Y^{\dagger}\bigr] + \sum_{j=1}^{N} {\mathbb{E}}\bigl[b_j(X_{{LR}},t) b_j^{\dagger}(X_{{LR}},t)\bigr] {{\boldsymbol{U}}}(t) \right) {\mathbb{E}}[Y Y^{\dagger}]^{-1}\\
    &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}\bigl[a(X_{{LR}},t) Y^{\dagger}\bigr] {\mathbb{E}}[Y Y^{\dagger}]^{-1}, \qquad  \bigl(\text{by Lemma }\ref{lemma::vanish_diffusion} \text{ and } {{\boldsymbol{U}}}(t) = {{\boldsymbol{V}}}(t) \bigr) \\
    &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}\left[ \left(\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {\mathcal{P}}_{{{\boldsymbol{V}}}(t)}\right) X_{{LR},t}  Y^{\dagger} + {\mathcal{P}}_{{{\boldsymbol{V}}}(t)} a(X_{{LR},t},t) Y^{\dagger} \right] {\mathbb{E}}[Y Y^{\dagger}]^{-1}, \qquad (\text{by {Lemma~\ref{lemma::vanish_diffusion}}} ) \\
    &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\left( \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {\mathcal{P}}_{{{\boldsymbol{V}}}(t)}\right) {{\boldsymbol{U}}}(t) = i {{\boldsymbol{F}}}(t) {{\boldsymbol{V}}}(t)  = \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{V}}}(t).
    \end{split}$$ Therefore $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{V}}}(t) = \frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{U}}}(t)$ and the conclusion holds.

Error bound
-----------

We shall quantify error defined by $$E(t) := {\left\lVert {\mathbb{E}}_{\mu_t} \bigl[ xx^{\dagger} \bigr] - {\mathbb{E}}_{\mu_{{LR},t}}
    \bigl[ xx^{\dagger}\bigr] \right\rVert}_{{HS}},$$ which measures the difference of second moment for two measures $\mu_t$ and $\mu_{{LR},t}$ in Hilbert-Schmidt norm. Recall that we have chosen the single test function $f(x) = x x^{\dagger}$.

\[thm::bound\] Assume that:

1.  Throughout the time-evolution of low-rank dynamics for $t\in {[0,T]}$, $$d_{{\mathsf{F}}}\bigl( {{\mathcal{A}}}_{{LR},t}
            \mu_{{LR},t},\ {\mathcal{A}}_{t} \mu_{{LR},t} \bigr) \equiv {\left\lVert\sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\right\rVert}_{{HS}} \le {\epsilon}^2.$$ Recall that $\theta_t$ is the pullback of the restriction of $\mu_{{LR},t}$ to $\operatorname{Ran}({{\boldsymbol{U}}}(t))$.

2.  There exists a function $\gamma(t)$ such that, for $f(x) = xx^{\dagger}$, $${\left\lVert {\mathbb{E}}_{\mu_t} \bigl[ {\mathcal{A}}_{t}^* f \bigr] -
                {\mathbb{E}}_{\mu_{{LR},t}}\bigl[{\mathcal{A}}_t^{*} f \bigr] \right\rVert}_{{HS}} \le \gamma(t)
            E(t) \equiv \gamma(t) {\left\lVert {\mathbb{E}}_{\mu_t} \bigl[ f \bigr] -
                {\mathbb{E}}_{\mu_{{LR},t}} \bigl[ f\bigr] \right\rVert}_{{HS}}.$$

Then, the error satisfies the integral inequality $$E(t) \le E(0) + {\epsilon}^2 t + \int_{0}^{t} \gamma(s) E(s) {\,\mathrm{d}}s,$$ and thus by [Gr[ö]{}nwall’s inequality]{}, $$E(t) \le \bigl(E(0) + {\epsilon}^2 t\bigr) \exp\left(\int_{0}^{t} \gamma(s){\,\mathrm{d}}s \right).$$

The proof follows from standard error analysis for time evolution equations: $$\begin{split}
    E(t) & \equiv {\left\lVert{\mathbb{E}}_{\mu_t} \bigl[ f \bigr] - {\mathbb{E}}_{\mu_{{LR},t}}\bigl[ f\bigr]\right\rVert}_{{HS}} \\
    & \le E(0) + {\left\lVert\int_{{{\mathbb{C}}^{{n}}}} f(x) \bigl(\mu_t - \mu_0\bigr)({\mathrm{d}}x) - f(x) \bigl(\mu_{{LR},t} - \mu_{{LR},0}\bigr)({\mathrm{d}}x)   \right\rVert}_{{HS}}\\
    & = E(0) + {\left\lVert\int_{{{\mathbb{C}}^{{n}}}} f(x) \left(\int_{0}^{t}{\mathcal{A}}_s \mu_s{\,\mathrm{d}}s \right)({\mathrm{d}}x) - f(x) \left( \int_{0}^{t} {\mathcal{A}}_{{LR},s} \mu_{{LR},s}{\,\mathrm{d}}s \right)({\mathrm{d}}x)   \right\rVert}_{{HS}}\\
    &= E(0) + {\left\lVert\int_{0}^{t} \left( \int_{{{\mathbb{C}}^{{n}}}} f(x) \left({\mathcal{A}}_s \mu_s \right)({\mathrm{d}}x) - f(x) \left({\mathcal{A}}_{{LR},s} \mu_{{LR},s} \right)({\mathrm{d}}x)\right)\   {\mathrm{d}}s \right\rVert}_{{HS}}  \\
    &= E(0) + {\left\lVert\int_{0}^{t} \left( \int_{{{\mathbb{C}}^{{n}}}} \left({\mathcal{A}}_{s}^{*} f\right) (x) \mu_s({\mathrm{d}}x) - \left({\mathcal{A}}_{{LR},s}^{*}  f\right)(x)  \mu_{{LR},s} ({\mathrm{d}}x)\right)\ {\mathrm{d}}s   \right\rVert}_{{HS}} \\
    &\le E(0) + \int_{0}^{t}  {\left\lVert{\mathbb{E}}_{\mu_s} \bigl[{\mathcal{A}}_s^{*} f \bigr] - {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_s^{*} f\bigr] + {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_s^{*} f\bigr] - {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_{{LR},s}^{*} f\bigr]  \right\rVert}_{{HS}}\ {\mathrm{d}}s \\
    &\le E(0) + \int_{0}^{t}  {\left\lVert{\mathbb{E}}_{\mu_s} \bigl[{\mathcal{A}}_s^{*} f \bigr] - {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_s^{*} f\bigr]\right\rVert}_{{HS}}\ {\mathrm{d}}s + \int_{0}^{t} {\left\lVert {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_s^{*} f\bigr] - {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_{{LR},s}^{*} f\bigr]  \right\rVert}_{{HS}}\ {\mathrm{d}}s \\
    & = E(0) + \int_{0}^{t}  {\left\lVert{\mathbb{E}}_{\mu_s} \bigl[{\mathcal{A}}_s^{*} f \bigr] - {\mathbb{E}}_{\mu_{{LR},s}}\bigl[{\mathcal{A}}_s^{*} f\bigr]\right\rVert}_{{HS}}\ {\mathrm{d}}s 
    + \int_{0}^{t} d_{{\mathsf{F}}}\left({\mathcal{A}}_s \mu_{{LR},s}, {\mathcal{A}}_{{LR},s} \mu_{{LR},s} \right)\ {\mathrm{d}}s \\ 
    &\le E(0) + \int_{0}^{t} \gamma(s) E(s)\ {\,\mathrm{d}}s + \int_{0}^{t} {\epsilon}^2\ {\mathrm{d}}s. \\
    \end{split}$$ Thus we have proved the integral inequality; the rest of the conclusion follows from [Gr[ö]{}nwall’s inequality]{}.

The estimate above holds for arbitrary ${\epsilon}$, however in practice, the above bound is most useful when ${\epsilon}$ is a small number. This corresponds to that the diffusivity function ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} b_j({{\boldsymbol{U}}}(t)y,t) =\mathcal{O}({\epsilon})$. For instance, SDE system with small noise falls into this type. When ${\epsilon}$ is large, of course, this indicates that low-rank approximation fails. One might need to use higher rank to get an accurate approximation of the dynamics, or even use full rank ([*i.e.*]{}, solving the original dynamics).

From the proof, the assumption (2) in [Theorem \[thm::bound\]]{} is natural, in order to use [Gr[ö]{}nwall’s inequality]{}. To better illustrate the assumption, [[we provide here a concrete example to give explicit form of $\gamma(t)$.]{}]{}

If SDE in [Eq. ]{} has linear drift and diffusion functions, that is, $$a(x,t) = {\boldsymbol{\Lambda}}(t) x, \qquad b_j(x,t) = {\boldsymbol{\Theta}}_j(t) x,$$ where bold Greek letters ${\boldsymbol{\Lambda}}$ and ${\boldsymbol{\Theta}}_j$ are time-dependent matrices on ${{\mathbb{C}}^{{n}}}$.

By Lemma \[lemma::generator\], one could straightforwardly find that $$\gamma(t) = 2{\left\lVert{\boldsymbol{\Lambda}}(t) \right\rVert}_{2} + \sum_{j=1}^{N} {\left\lVert{\boldsymbol{\Theta}}_j(t)\right\rVert}_2^2$$ satisfies the assumption.

Another example of $\gamma(t)$ for unraveling scheme will be given in the next section below.

Consider that if the rank $r = {n}$, then ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} = 0$ for all $t$, so that ${\epsilon}= 0$. If we further let $E(0) = 0$, then $E(t) = 0$ for all $t$. This result is consistent with intuition, since in this case the low-rank dynamic is exactly the original SDE. Though the rank $r$ does not appear explicitly in the error estimate, it is implicitly hidden inside ${\epsilon}$. Loosely speaking, the larger the rank $r$, the smaller the ${\epsilon}$. [[ Even though the relation between ${\epsilon}$ and $r$ is not analytically given, [Theorem \[thm::bound\]]{} can still be useful in practical simulation and in designing adaptive scheme. Numerically what we need is to set up error tolerance ${\epsilon}$, and then compute ${\left\lVert\sum_{j=1}^{N} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}\right\rVert}_{{HS}}$ on-the-fly in the dynamics; if this quantity is close to ${\epsilon}^2$, then it indicates that the rank chosen is not large enough anymore and we should adaptively increase the rank $r$ in order to control the error.]{}]{} This idea has been used in [@LeBris13] in numerically solving [Lindblad]{} equation by the deterministic low-rank approximation.

Connections to unraveling of Lindblad equations {#sec::discuss::unravel}
===============================================

In this section, we shall discuss the relationship between unraveling of [Lindblad]{} equation and dynamical low-rank approximation method; in particular, we establish a commuting diagram for the action of unraveling and the action of dynamical low-rank approximation under certain conditions.

[Lindblad]{} equation and stochastic unraveling
-----------------------------------------------

[Lindblad]{} equation, one of the most popular quantum master equations for open quantum systems, has the following form [@Lindblad76] $$\label{eqn::lb}
\dot{{{\boldsymbol{\rho}}}} = {\mathcal{L}}({{\boldsymbol{\rho}}}) =  -i {\left[{{\boldsymbol{H}}}, {{\boldsymbol{\rho}}}\right]_{-}} + 
\sum_{k} \left( {{\boldsymbol{L}}}_{k} {{\boldsymbol{\rho}}}{{\boldsymbol{L}}}_{k}^{\dagger} - \frac{1}{2} {\left[{{\boldsymbol{L}}}_{k}^{\dagger} {{\boldsymbol{L}}}_{k},  {{\boldsymbol{\rho}}}\right]_{+}}  \right),$$ where Hermitian operator ${{\boldsymbol{H}}}$ is Hamiltonian, ${{\boldsymbol{L}}}_k$ are Lindblad operators, [[ ${\left[\cdot, \cdot\right]_{-}}$ is the commutator]{}]{} and ${\left[\cdot, \cdot\right]_{+}}$ is the anti-commutator. The second term on the right hand side models the interaction of the system with the environment.

![Diagram for unraveling and dynamical low-rank approximation for [Lindblad]{} equations.[]{data-label="fig::diagram"}](diagram.pdf){width="\textwidth"}

The [Lindblad]{} equation is often challenging to solve numerically due to its high-dimensionality. There are two major dimension-reduction approaches from the literature, summarized in the [Figure \[fig::diagram\]]{}: (stochastic) unraveling method (II) and (deterministic) dynamical low-rank approximation method ([slowromancap1@]{}).

[**Step ([slowromancap1@]{}): Dynamical low-rank approximation method**]{}. This method for [Lindblad]{} equation has been studied in [@LeBris13]: The ansatz $ {{\boldsymbol{\rho}}}_{{LR}}(t) = {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t)$ was used, where ${{\widetilde{{\boldsymbol{U}}}}}(t)$ satisfies the orthonormality constraint ${{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {{\widetilde{{\boldsymbol{U}}}}}(t) = {\mathrm{Id}}$; and ${{\widetilde{{\boldsymbol{\sigma}}}}}(t)$ is a $r\times r$ strictly positive matrix with trace one. Tilde is used to distinguish the ${{\boldsymbol{U}}}(t)$ and ${{\boldsymbol{\sigma}}}(t)$ for deterministic dynamical low-rank approximation and those in the [SDLR]{} method. Dynamical low-rank approximation method would lead to a coupled ODE system for ${{\widetilde{{\boldsymbol{U}}}}}(t)$ and ${{\widetilde{{\boldsymbol{\sigma}}}}}(t)$, which approximates the [Lindblad]{} equation of ${{\boldsymbol{\rho}}}(t)$. We shall revisit the result from [@LeBris13] below in [Theorem \[thm::lebris\]]{}, while dropping the trace-preserving constraint to ${{\widetilde{{\boldsymbol{\sigma}}}}}$. We will discuss the trace-preserving constraint further in [Section \[sub:trace\]]{}.

\[thm::lebris\] Consider the subspace $${\widetilde{{\mathcal{M}}}}_r := \left\{ {{\widetilde{{\boldsymbol{U}}}}}{{\widetilde{{\boldsymbol{\sigma}}}}}{{\widetilde{{\boldsymbol{U}}}}}^{\dagger}: {{\widetilde{{\boldsymbol{U}}}}}^{\dagger} {{\widetilde{{\boldsymbol{U}}}}}= {\mathrm{Id}}_{r\times r},\ {{\widetilde{{\boldsymbol{\sigma}}}}}> 0 \right\},$$ to approximate the manifold of positive matrices. Then by dynamical low-rank approximation, [*i.e.*]{}, by solving $$\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\boldsymbol{\rho}}}_{{LR}}(t) = \operatorname*{arg\,min}_{v \in {\mathcal{T}}_{{{\boldsymbol{\rho}}}_{{LR}}(t)} {\mathcal{M}}_r} {\left\lVert{\mathcal{L}}\bigl({{\boldsymbol{\rho}}}_{{LR}}(t)\bigr) - v \right\rVert}_{{HS}},$$ with further restriction on the tangent space that $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\widetilde{{\boldsymbol{U}}}}}(t) = i {{\boldsymbol{G}}}(t){{\widetilde{{\boldsymbol{U}}}}}(t)$, where Hermitian matrix ${{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} + {\text{h.c.}}$, we have the unique low-rank dynamics, $$\label{eqn::lebris}
\left\{
\begin{split}
\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\widetilde{{\boldsymbol{\sigma}}}}}(t) &= {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {\mathcal{L}}\bigl({{\boldsymbol{\rho}}}_{{LR}}(t)\bigr) {{\widetilde{{\boldsymbol{U}}}}}(t)\\
\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\widetilde{{\boldsymbol{U}}}}}(t) &= {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {\mathcal{L}}\bigl({{\boldsymbol{\rho}}}_{{LR}}(t) \bigr) {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t)^{-1}. \\
\end{split}\right.$$

The proof is almost the same as the derivation in [@LeBris13]. It is provided in \[app::proof\_lebris\] for readers’ convenience. It should be noticed that if the condition ${{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\text{h.c.}}$ is not imposed, the low-rank dynamics is not unique as additional degrees of freedom exist, called $A$ in [@LeBris13]. The above dynamics could be viewed as a special case by choosing $A= 0$ in that reference. Just to be clear, note that this $A$ is not the drift function in this paper; it is a notation used in [@LeBris13]. However, the major difference comes from the trace-preserving condition in [@LeBris13], which takes the subspace as ${\widetilde{{\mathcal{M}}}}_r := \bigl\{ {{\widetilde{{\boldsymbol{U}}}}}{{\widetilde{{\boldsymbol{\sigma}}}}}{{\widetilde{{\boldsymbol{U}}}}}^{\dagger}: {{\widetilde{{\boldsymbol{U}}}}}^{\dagger} {{\widetilde{{\boldsymbol{U}}}}}=
  {\mathrm{Id}}_{r\times r},\ {{\widetilde{{\boldsymbol{\sigma}}}}}> 0,\ \operatorname{tr}({{\widetilde{{\boldsymbol{\sigma}}}}}) = 1 \bigr\}$. As a consequence, $\lambda {\mathrm{Id}}_{r\times r}$ term in equation (5) in [@LeBris13] does not appear in our expression.

**Step ([slowromancap2@]{}): Unraveling**. Recall that unraveling means a stochastic wave-equation that recovers the evolution of density matrices ${{\boldsymbol{\rho}}}(t)$ in expectation. More specifically, it looks for a stochastic process $X_t(\omega)$ on the Hilbert space ${{\mathbb{C}}^{{n}}}$, such that the expectation ${{\boldsymbol{\rho}}}(t):={\mathbb{E}}\bigl[X_t(\omega) X_t^{\dagger}(\omega)\bigr]$ solves the [Lindblad]{} equation. Multiple choices exist for the unraveling stochastic process: quantum state diffusion (QSD) [@Gisin92], linear quantum state diffusion (LQSD) [@Brun00] and quantum jump process [@Dalibard92]. We restrict the stochastic unraveling to SDE type in the discussion below, that is, we shall only consider unraveling of the form as in [Eq. ]{}. Still various choices exist, while the coefficients satisfy the relation stated in the following Lemma.

\[lemma::lb\_unravel\] If $X_t(\omega)$ in [Eq. ]{} is a stochastic unraveling of a [Lindblad]{} equation $\dot{{{\boldsymbol{\rho}}}} = {\mathcal{L}}({{\boldsymbol{\rho}}})$ if and only if $$\label{eqn::lb_equiv_cond}
      a(x,t) x^{\dagger} + x a^{\dagger}(x,t) + \sum_{j=1}^{N} b_j(x,t) b_j^{\dagger}(x,t) = {\mathcal{L}}\bigl(xx^{\dagger}\bigr),\qquad \forall x\in {{\mathbb{C}}^{{n}}},\ \forall t\in {[0,T]}.$$

Among various options of unraveling, two most popular choices are:

-   Linear quantum state diffusion (LQSD) refers to the choice $$\left\{\begin{split}
            a(x) &= \left(-i {{\boldsymbol{H}}}- \frac{1}{2} \sum_{k} {{\boldsymbol{L}}}_k^{\dagger} {{\boldsymbol{L}}}_k \right) x\\
            b_{k,1}(x) &= \frac{1}{\sqrt{2}} {{\boldsymbol{L}}}_k x \\
            b_{k,2}(x) &= \frac{i}{\sqrt{2}} {{\boldsymbol{L}}}_k x.\\
          \end{split}\right.$$ So that the autonomous SDE is given by $${\,\mathrm{d}}X_t = a(x) {\,\mathrm{d}}t +\sum_{k} \left( b_{k,1}(x){\,\mathrm{d}}W_{k,1} + b_{k,2}(x)  {\,\mathrm{d}}W_{k,2}\right),$$ where $W_{k,1}$ and $W_{k,2}$ are independent standard (real-valued) Brownian motions. If we have $N$ [Lindblad]{} operators, then we have $2N$ (real-valued) Brownian motions in SDE. If one combines $W_{k,1}$ and $W_{k,2}$ together to form a complex-valued Brownian motion $W_{k,1} + i W_{k,2}$, then there are only $N$ diffusion terms.

-   Quantum state diffusion (QSD) refers to the choice $$\left\{\begin{split}
            a(x) &= \left(-i {{\boldsymbol{H}}}+ \sum_{k} {\left\langlex, {{\boldsymbol{L}}}_k^{\dagger} x\right\rangle} {{\boldsymbol{L}}}_k - \frac{1}{2} {{\boldsymbol{L}}}_k^{\dagger} {{\boldsymbol{L}}}_k - \frac{1}{2} {\left\lvert{\left\langlex, {{\boldsymbol{L}}}_k^{\dagger} x\right\rangle}\right\rvert}^2 \right) x\\
            b_{k,1}(x) &= \frac{1}{\sqrt{2}} \bigl({{\boldsymbol{L}}}_k - {\left\langlex, {{\boldsymbol{L}}}_k x\right\rangle} \bigr) x \\
            b_{k,2}(x) &=  \frac{i}{\sqrt{2}} \bigl({{\boldsymbol{L}}}_k - {\left\langlex, {{\boldsymbol{L}}}_k x\right\rangle} \bigr) x. \\
          \end{split}\right.$$

[[One could easily verify that the above two choices satisfy [Eq. ]{}.]{}]{}

[**Step ([slowromancap3@]{}):**]{} In [@LeBris15], the authors also considered an unraveling scheme resulting from the (deterministic) low-rank approximation to [Lindblad]{} equations. As the low-rank approximation preserves the structure of the equation, the unraveling is similar to the above discussions. This low-rank unraveling could be used in control variate for Monte Carlo method for [Lindblad]{} equations.

A commuting diagram for low-rank approximation and unraveling
-------------------------------------------------------------

The route in [@LeBris15] is from [Lindblad]{} equation to low-rank quantum master equation and then to the unraveling of low-rank QME (step ([slowromancap1@]{}) to step ([3]{}) in [Figure \[fig::diagram\]]{}). It is thus natural to consider the alternative route: [*i.e.*]{}, finding the unraveling of [Lindblad]{} equation and then applying [stochastic dynamical [low-rank]{} approximation]{} method that we developed in [Section \[sec::sdlr\]]{} and \[sec::example\]); in [Figure \[fig::diagram\]]{}, this refers to the route from step ([slowromancap2@]{}) to step ([slowromancap4@]{}). One immediate question is whether these two routes commute, in the sense that they end up with the same equation for low-rank unraveling. Due to the non-uniqueness of unraveling scheme, in general, the answer is negative. Perhaps, a more specific and reasonable question to ask is that given the unraveling of [Lindblad]{} equation, after applying [SDLR]{} method, whether its statistical average recovers the low-rank QME (such as that derived in [@LeBris13] using deterministic low-rank approximation).

The answer is positive with slight modification in (deterministic) dynamical low-rank approximation method in [@LeBris13]. See [Theorem \[thm::lebris\]]{} for details of the [[modification in constraints]{}]{} as well as the resulting low-rank QME. Moreover, it turns out that such commuting diagram does not depend on the unraveling scheme chosen for [Lindblad]{} equation. The result is summarized in [Theorem \[thm::commuting\]]{}.

\[thm::commuting\] For any unraveling scheme in SDE form of [Lindblad]{} equation (see [Lemma \[lemma::lb\_unravel\]]{}), the low-rank unraveling obtained after applying the [stochastic dynamical [low-rank]{} approximation]{} method, is an unraveling scheme of low-rank quantum master equation obtained via (deterministic) dynamical low-rank approximation method given by [Eq. ]{} in Theorem \[thm::lebris\].

By applying [SDLR]{} (step ([slowromancap4@]{}) in [Figure \[fig::diagram\]]{}) to unraveling scheme of [Lindblad]{} equation and using [Eq. ]{}, one could obtain the following result. $$\left\{\begin{split}
\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\boldsymbol{U}}}(t) &= {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathcal{L}}\bigl({\mathbb{E}}\bigl[ {{\boldsymbol{U}}}(t) Y_t Y_t ^{\dagger} {{\boldsymbol{U}}}(t)^{\dagger} \bigr]\bigr) {{\boldsymbol{U}}}(t) {\mathbb{E}}\bigl[Y_t Y_t^{\dagger} \bigr]^{-1}\\
{\,\mathrm{d}}Y_t &=  {{\boldsymbol{U}}}^{\dagger}(t) a( {{\boldsymbol{U}}}(t) Y_t, t)\ {\,\mathrm{d}}t + \sum_{j=1}^{N} {{\boldsymbol{U}}}^{\dagger}(t) b_j({{\boldsymbol{U}}}(t) Y_t, t)\ {\,\mathrm{d}}W_j. \\
\end{split}\right.$$

Next we will verify that the above coupled ODE-SDE system does play the role of unraveling of low-rank QME in [Eq. ]{}, which is step ([slowromancap3@]{}) in [Figure \[fig::diagram\]]{}. By denoting ${{\boldsymbol{\sigma}}}(t) := {\mathbb{E}}\bigl[Y_t Y^{\dagger}_t\bigr]$, and ${{\boldsymbol{\rho}}}_{{LR}}(t) := {\mathbb{E}}\bigl[ X_{{LR},t} X_{{LR},t}^{\dagger} \bigr]\equiv {{\boldsymbol{U}}}(t) {{\boldsymbol{\sigma}}}(t) {{\boldsymbol{U}}}^{\dagger}(t)$, one could find that $$\begin{split}
\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\boldsymbol{\sigma}}}(t) &=  
{\mathbb{E}}\left[Y_t a^{\dagger}({{\boldsymbol{U}}}(t) Y_t,t) {{\boldsymbol{U}}}(t) + {{\boldsymbol{U}}}^{\dagger}(t) a({{\boldsymbol{U}}}(t) Y_t, t) Y_t^{\dagger} + \sum_{j=1}^{N} {{\boldsymbol{U}}}^{\dagger}(t) b_j({{\boldsymbol{U}}}(t) Y_t, t) b_j^{\dagger}({{\boldsymbol{U}}}(t)Y_t, t) {{\boldsymbol{U}}}(t) \right]\\
&= {{\boldsymbol{U}}}^{\dagger}(t){\mathbb{E}}\left[ \bigl({{\boldsymbol{U}}}(t) Y_t\bigr) a^{\dagger}({{\boldsymbol{U}}}(t) Y_t,t) 
+  a({{\boldsymbol{U}}}(t) Y_t, t) \bigl({{\boldsymbol{U}}}(t) Y_t\bigr)^{\dagger}
+ \sum_{j=1}^{N}  b_j({{\boldsymbol{U}}}(t) Y_t, t) b_j^{\dagger}({{\boldsymbol{U}}}(t)Y_t, t)
\right] {{\boldsymbol{U}}}(t) \\
&= {{\boldsymbol{U}}}^{\dagger}(t) {\mathbb{E}}\left[ {\mathcal{L}}\bigl({{\boldsymbol{U}}}(t) Y_t Y_t^{\dagger} {{\boldsymbol{U}}}^{\dagger}(t)\bigr) \right] {{\boldsymbol{U}}}(t)
\qquad 
\text{(use {Eq. \eqref{eqn::lb_equiv_cond}})}  \\
& = {{\boldsymbol{U}}}^{\dagger}(t) {\mathcal{L}}\bigl({{\boldsymbol{U}}}(t) {{\boldsymbol{\sigma}}}(t) {{\boldsymbol{U}}}^{\dagger}(t) \bigr) {{\boldsymbol{U}}}(t).
\end{split}$$ The time-evolution equation for ${{\boldsymbol{U}}}(t)$ can be rewritten, in terms of ${{\boldsymbol{\sigma}}}(t)$, as $$\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\boldsymbol{U}}}(t) ={\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\mathcal{L}}\bigl({{\boldsymbol{U}}}(t) {{\boldsymbol{\sigma}}}(t) {{\boldsymbol{U}}}^{\dagger}(t) \bigr) {{\boldsymbol{U}}}(t) {{\boldsymbol{\sigma}}}(t)^{-1}.$$ By comparing these two equations with [Eq. ]{}, one could conclude that the low-rank SDE, after applying [SDLR]{} method, exactly recovers the [Eq. ]{}, which means, the low-rank unraveling is the unraveling for low-rank QME given in [Eq. ]{} and the diagram in [Figure \[fig::diagram\]]{} indeed commutes in this sense. Also, note that in the above calculation, we haven’t used any specific choice of unraveling scheme, thus the conclusion is independent of unraveling scheme chosen for [Lindblad]{} equation.

It might not be surprising that the diagram commutes under the above conditions. If we consider ${{\boldsymbol{\rho}}}_{{LR}}(t) = {\mathbb{E}}\bigl[X_{{LR},t} X_{{LR},t}^{\dagger}\bigr]$, by our conditions, $X_{{LR},t} = {{\boldsymbol{U}}}(t) Y_t$, hence ${{\boldsymbol{\rho}}}_{{LR}}(t) = {{\boldsymbol{U}}}(t) {\mathbb{E}}\bigl[Y_t Y_t^{\dagger}\bigr]
  {{\boldsymbol{U}}}^{\dagger}(t)$, which is consistent with the ansatz used in [@LeBris13]. What is interesting is that the commuting diagram result is independent of any unraveling for [Lindblad]{} equation, which somewhat shows that the low-rank dynamics we derived from the perspective of dynamical low-rank approximation in space of signed measures, preserves the structure of [Lindblad]{} super-operator ${\mathcal{L}}$.

Recall that in Theorem \[thm::bound\], we have studied how the error between the original SDE and the low-rank approximation from [SDLR]{} method propagates with respect to time. For the case of unraveling of [Lindblad]{} equation, by using Lemma \[lemma::generator\] and Lemma \[lemma::lb\_unravel\], one could show that the growth rate for the error between the unraveling scheme and our corresponding low-rank approximation is bounded by $$\gamma(t)  = {\left\lVert{\mathcal{L}}\right\rVert}_{{HS}}:= \sup_{{{\boldsymbol{F}}}= {{\boldsymbol{F}}}^{\dagger},\ {\left\lVert{{\boldsymbol{F}}}\right\rVert}_{{HS}} = 1} {\left\lVert{\mathcal{L}}({{\boldsymbol{F}}})\right\rVert}_{{HS}}.$$

Discussion on trace-preserving restrictions {#sub:trace}
-------------------------------------------

In [@LeBris13; @LeBris15], the trace of low-rank approximated density matrix is required to be one, along the time-evolution. However, we did not consider such condition above in this section, nor in the derivation of low-rank dynamics in [Theorem \[thm::fx=xx\]]{} via [SDLR]{} method. One natural question is that what happens if trace-preserving condition is imposed in [SDLR]{} method. In our setting up and for the [Lindblad]{} equation case, the trace-preserving constraint in [SDLR]{} method should be $$\operatorname{tr}\bigl({{\boldsymbol{\rho}}}_{{LR}}(t)\bigr) \equiv \operatorname{tr}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[y y^{\dagger}\bigr] {{\boldsymbol{U}}}^{\dagger}(t)\right) = 1, \qquad \text{or equivalently} \qquad \operatorname{tr}\left({\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]\right) = 1,$$ because ${{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{U}}}(t) = {\mathrm{Id}}_{r\times r}$. From the constraint that $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} \operatorname{tr}\left({\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]\right) = 0$, one could derive that $$\label{eqn::trace_constraint}
 \operatorname{tr}\left( {\mathbb{E}}_{\theta_t}\bigl[ y A(y, t)^{\dagger} + A(y, t) y^{\dagger} + \sum_{j=1}^{M} B_j(y,t) B_j^{\dagger}(y,t) \bigr]   \right) = 0.$$

Recall that in the proof of [Theorem \[thm::fx=xx\]]{}, such constraint does not affect the first-order stationary condition with respect to ${{\boldsymbol{G}}}(t)$, hence ${{\boldsymbol{G}}}(t)$ still has the same form. Thus, the optimization problem in the derivation is still $$\min_{A,\ \bigl\{B_j\bigr\}_{j=1}^{M} } {\left\langle{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}, {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}\right\rangle},$$ where $$\begin{split}
{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} &= {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[ \bigl({{\boldsymbol{U}}}(t) A(y,t) - a({{\boldsymbol{U}}}(t)y,t)\bigr)y^{\dagger} {{\boldsymbol{U}}}^{\dagger}\bigr] + {\text{h.c.}}+ \sum_{j=1}^{M} {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t} \bigl[B_j(y,t) B_j^{\dagger}(y,t)\bigr] {{\boldsymbol{U}}}^{\dagger}(t) \\
& \qquad  - \sum_{j=1}^{N}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t} \bigl[b_j({{\boldsymbol{U}}}(t)y,t) b_j^{\dagger}({{\boldsymbol{U}}}(t)y,t)\bigr] {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} \\
&= {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\left[ A(y,t) y^{\dagger} + y A^{\dagger}(y,t) + \sum_{j=1}^{M} B_j(y,t) B_j^{\dagger}(y,t) \right] {{\boldsymbol{U}}}^{\dagger}(t) \\
& \qquad - {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathbb{E}}_{\theta_t}\left[ a({{\boldsymbol{U}}}(t) y, t) ({{\boldsymbol{U}}}(t) y)^{\dagger} + ({{\boldsymbol{U}}}(t) y) a^{\dagger}({{\boldsymbol{U}}}(t) y, t) + \sum_{j=1}^{N} b_j({{\boldsymbol{U}}}(t) y, t) b_j^{\dagger}({{\boldsymbol{U}}}(t) y, t) \right]{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} \\
&= {{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\left[ A(y,t) y^{\dagger} + y A^{\dagger}(y,t) + \sum_{j=1}^{M} B_j(y,t) B_j^{\dagger}(y,t) \right] {{\boldsymbol{U}}}^{\dagger}(t)  - {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\mathcal{L}}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t)\right) {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}.
\end{split}$$ In the last step, we have used [Eq. ]{}.

If we assume it is possible to achieve ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = 0$, then $${\mathbb{E}}_{\theta_t}\left[ A(y,t) y^{\dagger} + y A^{\dagger}(y,t) + \sum_{j=1}^{M} B_j(y,t) B_j^{\dagger}(y,t) \right]   = {{\boldsymbol{U}}}^{\dagger}(t) {\mathcal{L}}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t)\right) {{\boldsymbol{U}}}(t).$$ The constraint in [Eq. ]{} requires the the trace of left hand side is zero, while on the right hand side, in general, $$\operatorname{tr}\left\{ {{\boldsymbol{U}}}^{\dagger}(t) {\mathcal{L}}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t)\right) {{\boldsymbol{U}}}(t)\right\} = \operatorname{tr}\left\{ {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}  {\mathcal{L}}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t)\right) \right\} \neq 0,$$ even though $\operatorname{tr}\left\{ {\mathcal{L}}\left({{\boldsymbol{U}}}(t) {\mathbb{E}}_{\theta_t}\bigl[yy^{\dagger}\bigr]{{\boldsymbol{U}}}^{\dagger}(t)\right) \right\} \equiv 0$. One could conclude that in general, the minimization problem $$\min_{A,\ \bigl\{B_j\bigr\}_{j=1}^{M} } {\left\langle{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}, {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}\right\rangle} > 0,$$ under the constraint in [Eq. ]{}. This could be anticipated since one has to work on a smaller space during minimization. Hence, in general, the low-rank unraveling from [SDLR]{} with trace-preserving constraint is not optimal in the sense of [Theorem \[thm::fx=xx\]]{}.

In fact, from the optimization problem, it is not straightforward how to minimize the ${\left\langle{\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}, {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{C}}}{\mathcal{P}}_{{{\boldsymbol{U}}}(t)}\right\rangle}$ under the trace constraint for functions $A$ and $B_j$. After all, we need to work on the quadratic variational problem on functional space with trace constraint, not on matrices as in [@LeBris13].

While it is perhaps desirable to have a trace-preserving dynamical low-rank approximation of the density matrix, [[our choice of not considering trace-preserving constraint can be justified via better approximating measurement outcome.]{}]{} Since the expected measurement outcome for observable ${{\boldsymbol{O}}}$ is $\operatorname{tr}({{\boldsymbol{O}}}{{\boldsymbol{\rho}}}) = {\left\langle{{\boldsymbol{O}}}, {{\boldsymbol{\rho}}}\right\rangle}$, to get an accurate approximation of the expectation, it is sufficient that ${{\boldsymbol{\rho}}}_{{LR}}(t)$ is close to ${{\boldsymbol{\rho}}}(t)$ in Hilbert-Schmidt norm, without the requirement of trace-preserving constraint. One could also normalize the resulting density matrix from the [SDLR]{} as a postprocessing step.

Methods selection and control variate
-------------------------------------

From the commuting diagram, one might question the usefulness of low-rank unraveling in practice. Since in the low-rank unraveling, at each time step, one needs to store $\mathcal{O}({n}\times r  + N_s\times r)$ data for ${{\boldsymbol{U}}}(t)$ and random variable $Y$, where $N_s$ is the sample size. For the deterministic low-rank dynamics, it only needs to store $\mathcal{O}({n}\times r + r^2)$ for ${{\widetilde{{\boldsymbol{U}}}}}(t)$ and ${{\widetilde{{\boldsymbol{\sigma}}}}}(t)$.

When $r = \mathcal{O}(1)$ is useful to approximate the full dynamics, solving the deterministic low-rank approximation of [Lindblad]{} equation is a better choice as one does not need to simulate many sample paths to get statistical averages and the simulation of deterministic low-rank dynamics has smaller memory and computational cost.

On the other hand, when rank $r$ requires to be large in order to approximate the system accurately, then it becomes inefficient to solve the deterministic low-rank approximation for [Lindblad]{} equation. It is advantageous to turn to stochastic approximation and to use unraveling of [Lindblad]{} equation.

However, one could consider using control variate method [@LeBris15] to facilitate the simulation, that is, to use $$\overline{{{\boldsymbol{\rho}}}} = \overline{{{\boldsymbol{\rho}}}}_{MC} + \lambda({{{\boldsymbol{\rho}}}}_{{LR}}  - \overline{{{\boldsymbol{\rho}}}}_{\text{LRMC}}),$$ where $\overline{{{\boldsymbol{\rho}}}}_{MC}$ is obtained via unraveling scheme; ${{\boldsymbol{\rho}}}_{{LR}}$ is obtained by solving deterministic low-rank dynamics of [Lindblad]{} equation; $\overline{{{\boldsymbol{\rho}}}}_{\text{LRMC}}$ is obtained by solving low-rank unraveling scheme. The rank $r = \mathcal{O}(1)$ to ensure that one could simulate both deterministic and stochastic low-rank dynamics. The low-rank dynamics cannot provide accurate approximation, however, one could still simulate low-rank dynamics to achieve reduction of variance, by choosing correct parameter $\lambda$. Please refer to [@LeBris15] for details and numerical performance.

Generally speaking, for Fokker-Planck equation and [Lindblad]{} equation with continuum state space (e.g., probability on $\mathbb{R}^d$ or density matrix over $L^2(\mathbb{R}^d)$), one would not want to solve a deterministic low-rank dynamics for them directly, since that is still a PDE in potentially high dimension when $r$ is not so small (even when $r = 6$, solving such a PDE in $6$ dimension is already rather challenging with standard methods). When the low-rank approximation is accurate, one could choose to use low-rank SDE (derived by [SDLR]{} method) to achieve reduction of complexity in model; otherwise, one has to simulate the original SDE. Of course, this discussion only involves which model to solve; detailed numerical methods and algorithmic implementation would still make a significant difference to the overall performance.

Numerical experiments {#sec::numerics}
=====================

In this section, we will validate our method using numerical examples of some high-dimensional SDEs: high-dimensional geometric Brownian motion, stochastic Burgers’ equation and unraveling of quantum damped harmonic oscillator. For the first two examples, we will also compare [SDLR]{} method with DO method. It could be observed that [SDLR]{} method has comparable performance in approximating the mean when the rank is chosen correctly, compared with DO method; and it performs better in approximating the second moment. We will measure the relative error as the indicator of performance of low rank approximation:

-   the relative error for the mean (linear) is defined as ${\left\lVert{\mathbb{E}}_{\mu_t} \bigl[x\bigr] -
      {\mathbb{E}}_{\mu_{{LR},t}}\bigl[x\bigr]\right\rVert}/{\left\lVert{\mathbb{E}}_{\mu_t} \bigl[x\bigr]\right\rVert}$;

-   the relative error for the second moment (quadratic) is defined as ${\left\lVert {\mathbb{E}}_{\mu_{t}} \bigl[xx^{\dagger}\bigr] -
      {\mathbb{E}}_{\mu_{{LR},t}}\bigl[
      xx^{\dagger}\bigr]\right\rVert}_{{HS}}/{\left\lVert{\mathbb{E}}_{\mu_t}\bigl[xx^{\dagger}\bigr]\right\rVert}_{{HS}}$.

Let us comment on some details of the numerical implementation. For simplicity, Euler-Maruyama method is used as stochastic integrator. An order-one deterministic numerical scheme in [@Higham96] is employed to preserve orthogonality of ${{\boldsymbol{U}}}(t)$. The inverse of ${\mathbb{E}}_{\theta_t} \bigl[ yy^{\dagger}\bigr]$ in [Eq. ]{} will cause numerical instability when its condition number is large. This problem also appears in dynamical orthogonal method [@Sapsis09]. Paper [@Babaee17] suggested to use pseudo-inverse to maintain the algorithmic stability, which is also adopted here in the numerical simulation. If the chosen rank is representative (not over-estimating the rank), then the pseudo-inverse should be simply the matrix inverse.

Geometric Brownian motion
-------------------------

Consider the geometric Brownian motion of the form $${\,\mathrm{d}}X_t = {\boldsymbol{\Lambda}} X_t{\,\mathrm{d}}t + {\boldsymbol{\Theta}} X_t{\,\mathrm{d}}W_t,$$ where $X_t \in {\mathbb{C}}^{20}$, ${\boldsymbol{\Lambda}}$, ${\boldsymbol{\Theta}} \in {\mathbb{C}}^{20\times 20}$.

A rank-5 initial condition is used with $X_{0} = x_k$ with probability $p_k$ where $\{x_k\}_{k=1}^{5}$ are randomly generated orthogonal vectors in ${\mathbb{C}}^{20}$ and $p_k \propto  \text{Poisson}(k-1,0.5)$ for $1\le k \le 5$. $\text{Poisson}(k,\lambda)$ represents the probability density function of Poisson distribution with rate $\lambda$ at value $k$. ${\boldsymbol{\Lambda}}$ is of the form ${\boldsymbol{Q}} {{\boldsymbol{D}}}{\boldsymbol{Q}}^{\dagger}$ where ${{\boldsymbol{D}}}$ is a randomly generated diagonal matrix with diagonal elements uniformly distributed in the interval $[-4.5, -0.5]$ and ${\boldsymbol{Q}}$ is a randomly generated orthogonal matrix; ${\boldsymbol{\Theta}} = \sqrt{0.05}\, {\mathrm{Id}}_{20 \times 20}$. It is not difficult to prove that for such geometric Brownian motion, the second moment $ {\mathbb{E}}_{\mu_t} \bigl[xx^{\dagger}\bigr]$ decays to $0$ as $t\rightarrow\infty$ (known as mean-square stability). The result is visualized in Figures \[fig::gbm\_spectrum\], \[fig::gbm\_SDLR\] for both [SDLR]{} method and DO method. It should be remarked that to faithfully compare [SDLR]{} with DO under the same rank $r$, [*e.g.*]{}, $r = 5$, the matrix ${{\boldsymbol{U}}}(t)$ in [SDLR]{} has dimension $20\times 5$, whereas it has dimension $20\times 4$ for DO method since $\bar{X}$ for DO method should account for one rank and contributes $n=20$ degrees of freedom.

In numerical experiment, the sample size is $10^5$ and time step is $\frac{1}{300}$. From [Figure \[fig::gbm\_spectrum\]]{}, the mean-square stability is clearly observed. In [Figure \[fig::gbm\_SDLR\]]{}, for [SDLR]{} method, when rank increases, the relative error decreases for both mean and second moment. Since the fifth eigenvalue is extremely small compared with others, it is reasonable that choosing rank $5$ does not significantly improve the accuracy further. For DO method, even for small rank, the relative error of mean is small; the relative error for the second moment decreases as rank increases, which is expected. As can be seen, DO method captures the mean better and [SDLR]{} method captures the second moment better. This finding is consistent with the theoretical derivation.

![Spectrum of ${\mathbb{E}}_{\mu_t} \bigl[xx^{\dagger}\bigr]$ in $\log 10$ scale (the five largest eigenvalues) for high-dimensional geometric Brownian motion.[]{data-label="fig::gbm_spectrum"}](GBM_spectrum.pdf){width="48.00000%"}

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving high-dimensional geometric Brownian motion.[]{data-label="fig::gbm_SDLR"}](GBM_relative_error_of_mean_for_SDLR.pdf "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving high-dimensional geometric Brownian motion.[]{data-label="fig::gbm_SDLR"}](GBM_relative_error_of_second_moment_for_SDLR.pdf "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving high-dimensional geometric Brownian motion.[]{data-label="fig::gbm_SDLR"}](GBM_relative_error_of_mean_for_DO.pdf "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving high-dimensional geometric Brownian motion.[]{data-label="fig::gbm_SDLR"}](GBM_relative_error_of_second_moment_for_DO.pdf "fig:"){width="\textwidth"}

Stochastic Burgers’ equation
----------------------------

We will use the stochastic Burgers’ equation $h(z, t,\omega)$ of the form $$\begin{split}
& {\,\mathrm{d}}h = (\nu \partial_{z}^2 h - h \partial_z h) {\,\mathrm{d}}t + g(z) {\,\mathrm{d}}W \\
& h(z, 0, \omega) = h_0(z, \omega) \qquad h(0, t, \omega) = h(1, t, \omega) ,
\end{split}$$ for $(z, t) \in [0,1]\times [0,T]$. This example is adapted from [@Ozen16_SPDE]. Notice that ${\,\mathrm{d}}W$ is chosen as a scalar Brownian motion, independent of spatial coordinate $z$.

Due to the periodic boundary condition, by separating $(t, \omega)$ from variable $z$, $$h(z,t,\omega) = \sum_{k\in {\mathbb{Z}}} X_k(t,\omega) e^{2\pi i k z}.$$ Then stochastic Burgers’ equation could be viewed as a SDE on Hilbert space $L^2[0,1]$ with basis functions $\{e^{2\pi i k z} \}_{k\in {\mathbb{Z}}}$. The SDE has the form $$\begin{split}
{\,\mathrm{d}}X_{k}(t,\omega) &= \left(- (2\pi k)^2 \nu X_{k}(t,\omega) - \sum_{k'} X_{k-k'}(t,\omega) X_{k'}(t,\omega) (2\pi i k') \right)\ {\,\mathrm{d}}t + g_k\ {\,\mathrm{d}}W \\
X_{k}(0,\omega) &= {\left\langlee^{2\pi i k z}, h_0(z,\omega)\right\rangle} \qquad g_k = {\left\langlee^{2\pi i k z }, g(z)\right\rangle},
\end{split}$$ where the inner product ${\left\langlef_1, f_2\right\rangle} = \int_{0}^{1} f_1^{*}(z) f_2(z)\ {\,\mathrm{d}}z$ and ${\,\mathrm{d}}W$ is independent of the mode $k$.

In the direct numerical simulation, we truncate $k$ by letting $X_{k}(t,\omega) = 0$ for $|k|> \frac{n-1}{2}$ where $n$ is an odd positive integer. Then $X_{k}(t,\omega)$ can be stored in a ${\mathbb{C}}^{n}$ vector. With careful choice of initial condition and let $n\rightarrow \infty$, we would expect to have the solution of such truncated SDE converge to the true solution. Adapted from Example 4.1 in Ref. [@Ozen16_SPDE], let $\nu = 0.01$, and $$g(z) = \gamma \cos(2\pi z) = \frac{\gamma}{2} \left(e^{2\pi i z} + e^{2\pi i (-1) z}\right), \qquad \gamma = \frac{1}{10}.$$

As for initial condition, a rank-5 case is considered $$h_0(z, \cdot) = \left\{ \begin{split}
1, &\qquad  \text{with probability } p_1 \\
\sqrt{2} \sin(2\pi \floor{k/2} z), &\qquad \text{with probability } p_k,\ k = 2, 4\\
\sqrt{2} \cos(2\pi \floor{k/2} z), &\qquad \text{with probability } p_k,\ k = 3, 5 \\ 
\end{split}  \right.$$ where $p_k\propto \text{Poisson}(k-1,0.5)$ just like last example and $\floor{ \cdot }$ is the floor function.

In [Figure \[fig::sbe\_dimension\]]{}, it could be observed that the numerical result is stable with respect to dimension for truncation, [*i.e.*]{}, $n$. Therefore, it is justifiable to simply solve the truncated system with $n = 21$ by [SDLR]{} and DO method.

For $n = 21$, sample size $10^4$, time step $\frac{1}{200}$, the spectrum and relative error are visualized in Figures \[fig::sbe\_spec\] and \[fig::sbe\_error\] respectively. From the spectrum, the second moment tends to behave like a rank one matrix, since the largest eigenvalue almost keeps a constant while other eigenvalues roughly exponentially decay. In [Figure \[fig::sbe\_error\]]{}, for both methods, when rank increases, the relative error decreases, which is consistent with expectation. It could be seen that the performance of [SDLR]{} method and DO method is similar and comparable, in calculating both ${\mathbb{E}}_{\mu_t} \bigl[x\bigr]$ and ${\mathbb{E}}_{\mu_t} \bigl[xx^{\dagger}\bigr]$. As for more detailed performance comparison between these two methods, [Figure \[fig::sbe\_error\]]{} is not very informative, especially for second moment. Therefore, we additionally provide Table \[table::error\_SBE\], which gives the ratio between relative error for [SDLR]{} method and relative error for DO method at time $t = 1$, for various rank $r$. Then, it could be observed that DO method performs better in approximating the mean and [SDLR]{} method performs better in approximating the second moment. As a reminder, the result in Table \[table::error\_SBE\] can only be interpreted qualitatively, due to random fluctuation in simulation.

![Spectrum of ${\mathbb{E}}_{\mu_t} \bigl[xx^{\dagger}\bigr]$ in $\log 10$ scale (the five largest eigenvalues) for stochastic Burgers’ equation.[]{data-label="fig::sbe_spec"}](SBE_spectrum.pdf){width="48.00000%"}

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving stochastic Burgers’ equation.[]{data-label="fig::sbe_error"}](SBE_relative_error_of_mean_for_SDLR "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving stochastic Burgers’ equation.[]{data-label="fig::sbe_error"}](SBE_relative_error_of_second_moment_for_SDLR "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving stochastic Burgers’ equation.[]{data-label="fig::sbe_error"}](SBE_relative_error_of_mean_for_DO "fig:"){width="\textwidth"}

 

[0.48]{} ![Relative error for [SDLR]{} method and DO method in solving stochastic Burgers’ equation.[]{data-label="fig::sbe_error"}](SBE_relative_error_of_second_moment_for_DO "fig:"){width="\textwidth"}

   $r$   ratio for relative error of ${\mathbb{E}}_{\mu_t}\big[x\bigr]$   ratio for relative error of ${\mathbb{E}}_{\mu_t} \big[xx^{\dagger}\bigr]$
  ----- ---------------------------------------------------------------- ----------------------------------------------------------------------------
    3                                1.291                                                                  0.572
    4                                1.264                                                                  0.593
    5                                1.389                                                                  0.816

  : This table summarizes the ratio between relative error for [SDLR]{} method and the relative error for DO method, at time $t = 1$ for stochastic Burgers’ equation.[]{data-label="table::error_SBE"}

Quantum damped harmonic oscillator
----------------------------------

This example is to solve a simple quantum damped harmonic oscillator, see [*e.g.*]{}, [@Fujii12]. The [Lindblad]{} equation is given by

$$\dot{{{\boldsymbol{\rho}}}} = -i {\left[ \omega {\boldsymbol{d}}^{\dagger} {\boldsymbol{d}},  {{\boldsymbol{\rho}}}\right]_{-}} +
 \gamma_1 \left( {\boldsymbol{d}} {{\boldsymbol{\rho}}}{\boldsymbol{d}}^{\dagger} - \frac{1}{2} {\left[{\boldsymbol{d}}^{\dagger} {\boldsymbol{d}},  {{\boldsymbol{\rho}}}\right]_{+}}  \right) + 
 \gamma_2 \left( {\boldsymbol{d}}^{\dagger}{{\boldsymbol{\rho}}}{\boldsymbol{d}} - \frac{1}{2} {\left[{\boldsymbol{d}} {\boldsymbol{d}}^{\dagger}, {{\boldsymbol{\rho}}}\right]_{+}} \right),$$

where $\omega$ is angular frequency, ${\boldsymbol{d}}^{(\dagger)}$ are annihilation (creation) operator for harmonic oscillator. Adopting the bra-ket notation, let $\left\{\ {\left\lvertk\right\rangle}\ \right\}_{k=0, 1, \cdots, n-1}$ be the orthonormal basis of quantum states. Then the effect of operator ${\boldsymbol{d}}^{(\dagger)}$ is ${\boldsymbol{d}}^{\dagger}{\lvertk\rangle} = \sqrt{k+1}{\lvertk+1\rangle}$ and ${\boldsymbol{d}}{\lvertk\rangle} = \sqrt{k} {\lvertk-1\rangle}$. The special state ${\left\lvert0\right\rangle}$ usually refers to ground state and ${\left\lvertk\right\rangle}$ ($k\ge 1$) are known as excited states, in physical literatures.

We truncate the system by $n=21$ states, [*i.e.*]{}, from ground state ${\lvert0\rangle}$ to excited state ${\lvert20\rangle}$. Set $\omega = 1.0$ and initial condition is chosen as $\psi(0) = {\lvertk\rangle}$ with probability $p_k\propto \text{Poissson}(k, 0.5)$ for $k=0,1,\cdots, 4$; hence the density matrix ${{\boldsymbol{\rho}}}(0)$ is again a rank-5 matrix. This model has been tested for two parameter sets $\gamma_1 = 0.2$, $\gamma_2 = 0$ and $\gamma_1 = 0$, $\gamma_2 = 0.2$. In the first case, the environment acts as an annihilation operator to the system so that the system is moving to a lower energy state, whereas in the second case, the environment acts as a creation operator so that the system is moving up to a higher energy state. The above [Lindblad]{} equation is solved by both QSD and LQSD unraveling schemes, as well as the [SDLR]{} method for QSD and LQSD correspondingly (see [Section \[sec::discuss::unravel\]]{} for expressions of QSD and LQSD). For this example, the solution of [Lindblad]{} equation via numerical ODE integrator is used as the reference and treated as the exact solution.

The sample size is $3\times 10^5$ and time step is $\frac{1}{500}$. The results are visualized in [Figure \[fig::qdho\_test1\]]{} and \[fig::qdho\_test2\]. The spectrum is consistent with the physical intuition. In [Figure \[fig::qdho\_test1\]]{}, the highest eigenvalue roughly indicates the probability at ground state ${\lvert0\rangle}$, which steadily increases. This is consistent with the functioning of annihilation operator, [*i.e.*]{}, moving the quantum state to lower energy state. The spectrum in [Figure \[fig::qdho\_test2\]]{} can be explained in a similar way. When the rank increases from 3 to 5, relative error decreases in all figures. The solution at rank $5$ is comparable with the original unraveling SDE system ([*i.e.*]{}, $r=21$). In the first case $\gamma_1 = 0.2$ and $\gamma_2 = 0.0$, since the system is lowering down to ground state, it is expected that low-rank approximation should work better.

Conclusions and outlook {#sec::conclusion}
=======================

In this paper, we have proposed a tangent space projection method in the space of finite signed measures; it is termed as [stochastic dynamical [low-rank]{} approximation]{} method (or [SDLR]{} in abbreviation). Then by applying [SDLR]{} method, we have derived the low-rank dynamics of SDE in [Theorem \[thm::fx=xx\]]{}, obtaining an ODE-SDE coupled system as a low-rank approximation of the original high-dimensional SDE in the weak sense. We have also established a commuting diagram for the action of stochastic unraveling and dynamical low-rank approximation. This method has further been validated by error analysis. Three numerical examples have been provided in [Section \[sec::numerics\]]{} to demonstrate the good performance of this low-rank approximation method for dissipative physical systems.

There are some continuing interesting questions to explore. For instance, whether [SDLR]{} method can be extended to infinite-dimensional Hilbert space instead of ${{\mathbb{C}}^{{n}}}$. As is well-known, many interesting quantum master equations (as well as their unraveling schemes) evolve on an infinite-dimensional Hilbert space. Also, many SPDEs are essentially SDEs on infinite dimensional Hilbert space, [*e.g.*]{}, stochastic Burgers’ equation in the form given in last section. Though finite truncation is a must in numerical simulation in practice, it is still desirable to see whether our scheme could be applied to such SPDEs directly in theory. Another interesting question is to develop adaptive schemes for [SDLR]{} method that automatically adjust the rank on-the-fly. We shall leave these questions to future works.

Acknowledgment {#acknowledgment .unnumbered}
==============

This work is partially supported by the National Science Foundation under award DMS-1454939.

Additional proofs
=================

Tangent space of Stiefel manifold ${V_r({{\mathbb{C}}^{{n}}})}$ {#sec::stiefel_tangent}
---------------------------------------------------------------

Stiefel manifold on ${{\mathbb{C}}^{{n}}}$ is a collection of ${n}\times r$ complex-valued matrices with orthonormal columns; it is denoted by ${V_r({{\mathbb{C}}^{{n}}})}$. For a differentiable trajectory ${{\boldsymbol{U}}}(t)$ in the Stiefel manifold, one has ${{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{U}}}(t) = {\mathrm{Id}}_{r\times r}$ for all time $t$, thus $\dot{{{\boldsymbol{U}}}}^{\dagger}(t) {{\boldsymbol{U}}}(t) + {{\boldsymbol{U}}}^{\dagger}(t) \dot{{{\boldsymbol{U}}}}(t) = 0$ and consequently, tangent space ${\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}= \left\{{{\boldsymbol{V}}}\in {\mathbb{C}}^{{n}\times r}: {{\boldsymbol{V}}}^{\dagger} {{\boldsymbol{U}}}(t) + {{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{V}}}= 0 \right\}.$ Then we shall prove the following lemma.

The tangent space of Stiefel manifold is given by $${\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}= \left\{ i {{\boldsymbol{G}}}{{\boldsymbol{U}}}(t):\ {{\boldsymbol{G}}}^{\dagger} = {{\boldsymbol{G}}}\right\}.$$

Let us first prove ${\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}\subset \left\{i {{\boldsymbol{G}}}{{\boldsymbol{U}}}(t):\ {{\boldsymbol{G}}}\text{ is a Hermitian matrix} \right\}$. The other direction is trivial.

Suppose the columns of ${{\boldsymbol{U}}}(t)$ are denoted by $U_1, U_2, \cdots U_r$ respectively. One could complete the basis and find an orthonormal set $\left\{U_{r+1}, \cdots U_{{n}-1}, U_{{n}}\right\}$. Let ${{\boldsymbol{U}}}_{\perp} = \begin{bmatrix} U_{r+1} & \cdots &  U_{{n}-1} & U_{{n}}\end{bmatrix}$. Then $\begin{bmatrix} {{\boldsymbol{U}}}(t) & {{\boldsymbol{U}}}_{\perp}\end{bmatrix}$ is a unitary matrix. Let us extend the matrix ${{\boldsymbol{V}}}\in {\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}$ to be $\begin{bmatrix} {{\boldsymbol{V}}}& {\widetilde{{{\boldsymbol{V}}}}} \end{bmatrix}$ and define matrix ${\widetilde{{{\boldsymbol{G}}}}}$ by $\begin{bmatrix}{{\boldsymbol{V}}}& {\widetilde{{{\boldsymbol{V}}}}} \end{bmatrix} \begin{bmatrix} {{\boldsymbol{U}}}(t) & {{\boldsymbol{U}}}_{\perp}\end{bmatrix}^{\dagger} =: i {\widetilde{{{\boldsymbol{G}}}}}$. The choice of ${\widetilde{{{\boldsymbol{V}}}}}_{{n}\times ({n}-r)}$ does not play any role in the proof; it is only introduced to conveniently define ${\widetilde{{{\boldsymbol{G}}}}}$. Then $$i {\widetilde{{{\boldsymbol{G}}}}} = {{\boldsymbol{V}}}{{\boldsymbol{U}}}^{\dagger}(t) + {\widetilde{{{\boldsymbol{V}}}}} {{\boldsymbol{U}}}_{\perp}^{\dagger}.$$ Then one could calculate that $$\left\{\begin{split}
& {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\widetilde{{{\boldsymbol{G}}}}} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = (-i) {{\boldsymbol{U}}}(t) \left({{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{V}}}\right) {{\boldsymbol{U}}}^{\dagger}(t) \\
& {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)}  {\widetilde{{{\boldsymbol{G}}}}} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} = (-i) {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{V}}}{{\boldsymbol{U}}}^{\dagger}(t), 
\end{split}\right.$$ where ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} := {{\boldsymbol{U}}}(t) {{\boldsymbol{U}}}^{\dagger}(t)$ and ${\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} := {\mathrm{Id}}- {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}$. It could be straightforwardly verified that ${\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\widetilde{{{\boldsymbol{G}}}}} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)}$ is Hermitian, due to the assumption that ${{\boldsymbol{V}}}\in {\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}$. Let us define Hermitian matrix ${{\boldsymbol{G}}}$ by $${{\boldsymbol{G}}}:= {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\widetilde{{{\boldsymbol{G}}}}} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {\widetilde{{{\boldsymbol{G}}}}} {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} + {\mathcal{P}}_{{{\boldsymbol{U}}}(t)} {\widetilde{{{\boldsymbol{G}}}}}^{\dagger} {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} .$$ Then $$i {{\boldsymbol{G}}}{{\boldsymbol{U}}}(t) = \left( {{\boldsymbol{U}}}(t) \bigl({{\boldsymbol{U}}}^{\dagger}(t) {{\boldsymbol{V}}}\bigr) {{\boldsymbol{U}}}^{\dagger}(t) + {\mathcal{Q}}_{{{\boldsymbol{U}}}(t)} {{\boldsymbol{V}}}{{\boldsymbol{U}}}^{\dagger}(t) \right) {{\boldsymbol{U}}}(t) = {{\boldsymbol{V}}}.$$ Thus, ${\mathcal{T}}_{{{\boldsymbol{U}}}(t)} {V_r({{\mathbb{C}}^{{n}}})}\subset \left\{i {{\boldsymbol{G}}}{{\boldsymbol{U}}}(t):\ {{\boldsymbol{G}}}\text{ is a Hermitian matrix} \right\}$ and the proof is completed.

Equivalent choices of test functions {#sec:equivalence}
------------------------------------

\[lemma::choice\_equivalent\_F\] In defining [pseudometric]{} in [Eq. ]{}, the following two choices are equivalent,

1.  ${\mathsf{F}}= \left\{f(x) = xx^{\dagger}\right\}$ with Hilbert-Schmidt norm.

2.  $${\mathsf{F}}= \left\{f(x)={\left\langlex, {{\boldsymbol{O}}}x\right\rangle} \;\big\vert \;
      \text{Hermitian matrix } {{\boldsymbol{O}}}\text{ satisfies } {\lVert{{\boldsymbol{O}}}\rVert}_{{HS}} \le
      1\right\}$$ and the associated norm is simply absolute value.

Denote the [pseudometric]{} defined by choice (1) as $d_{{\mathsf{F}}_1}$ and the one defined by choice (2) as $d_{{\mathsf{F}}_2}$.

Firstly, for $\nu = \nu_1 - \nu_2 \in {\mathcal{M}}$, $$\begin{split}
 {\left\lvert\int {\left\langlex, {{\boldsymbol{O}}}x\right\rangle}\ {\,\mathrm{d}}\nu\right\rvert} &= {\left\lvert\int \operatorname{tr}(xx^{\dagger} {{\boldsymbol{O}}})\ {\,\mathrm{d}}\nu\right\rvert} = {\left\lvert\operatorname{tr}\left(\int xx^{\dagger}\ {\,\mathrm{d}}\nu {{\boldsymbol{O}}}\right)\right\rvert} \\
&\le {\left\lVert\int xx^{\dagger}\ {\,\mathrm{d}}\nu\right\rVert}_{{HS}} {\left\lVert{{\boldsymbol{O}}}\right\rVert}_{{HS}} \le {\left\lVert\int xx^{\dagger}\ {\,\mathrm{d}}\nu\right\rVert}_{{HS}}  = d_{{\mathsf{F}}_1}(\nu_1, \nu_2).
\end{split}$$ Then, it implies that $d_{{\mathsf{F}}_2}(\nu_1, \nu_2) \le d_{{\mathsf{F}}_1}(\nu_1, \nu_2)$. Secondly, we shall prove that the equality could be reached. The equality is reached when ${{\boldsymbol{O}}}= \frac{\int xx^{\dagger}\ {\,\mathrm{d}}\nu}{{\left\lVert\int xx^{\dagger}\ {\,\mathrm{d}}\nu\right\rVert}_{{HS}}}$, which is also a Hermitian matrix. Thus these two choices are equivalent.

Proof of [Theorem \[thm::lebris\]]{} {#app::proof_lebris}
------------------------------------

At fixed time $t$, the tangent space of ${{\boldsymbol{\rho}}}_{{LR}}(t) \equiv {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) \in {\mathcal{M}}_r$ is parametrized by Hermitian matrices ${{\boldsymbol{G}}}(t)$ and ${{\boldsymbol{\eta}}}(t)$ with the form $$\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{{\boldsymbol{\rho}}}}_{{LR}}(t) = i {{\boldsymbol{G}}}(t) {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) + {{\widetilde{{\boldsymbol{U}}}}}(t) {{\boldsymbol{\eta}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) -i  {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t)  {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {{\boldsymbol{G}}}(t),$$ where ${{\boldsymbol{G}}}(t) = {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} + {\text{h.c.}}$, while ${{\boldsymbol{\eta}}}(t)$ is any Hermitian matrix.

Denote ${{\boldsymbol{C}}}\equiv {\mathcal{L}}\bigl({{\boldsymbol{\rho}}}_{{LR}}(t)\bigr) - \left(i {{\boldsymbol{G}}}(t) {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) + {{\widetilde{{\boldsymbol{U}}}}}(t) {{\boldsymbol{\eta}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t)  - i  {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t)  {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {{\boldsymbol{G}}}(t)\right)$. Then we need to minimize ${\left\lVert{{\boldsymbol{C}}}\right\rVert}_{{HS}}^2 = {\left\langle{{\boldsymbol{C}}}, {{\boldsymbol{C}}}\right\rangle}_{{HS}}$ for varying ${{\boldsymbol{G}}}(t)$ and ${{\boldsymbol{\eta}}}(t)$. The first order stationary conditions with respect to ${{\boldsymbol{G}}}(t)$ and ${{\boldsymbol{\eta}}}(t)$ yield $$\left\{\begin{split}
0 &=\operatorname{tr}\left( {\left[ {{\boldsymbol{\rho}}}_{{LR}}(t),  {{\boldsymbol{C}}}\right]_{-}} \delta {{\boldsymbol{G}}}\right), \\
0 &= \operatorname{tr}\left({{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {{\boldsymbol{C}}}{{\widetilde{{\boldsymbol{U}}}}}(t) \delta {{\boldsymbol{\eta}}}\right).\\
\end{split}\right.$$ Similar to the argument in the proof of [Theorem \[thm::fx=xx\]]{}, since $\delta {{\boldsymbol{G}}}$ is any Hermitian matrix satisfying $\delta {{\boldsymbol{G}}}= {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} \delta {{\boldsymbol{G}}}{\mathcal{P}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)}+ {\text{h.c.}}$, from the first part in the last equation, one could obtain $$\begin{split}
0 &= {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {\left[ {{\boldsymbol{\rho}}}_{{LR}}(t),  {{\boldsymbol{C}}}\right]_{-}}  {\mathcal{P}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} \\
&= - {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {{\boldsymbol{C}}}{{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t) {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t). \\
\end{split}$$ After replacing ${{\boldsymbol{C}}}$ by its definition and some simplification, one could obtain $${\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {{\boldsymbol{G}}}(t) {\mathcal{P}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} = -i {\mathcal{Q}}_{{{\widetilde{{\boldsymbol{U}}}}}(t)} {\mathcal{L}}\bigl( {{\boldsymbol{\rho}}}_{{LR}}(t) \bigr) {{\widetilde{{\boldsymbol{U}}}}}(t) {{\widetilde{{\boldsymbol{\sigma}}}}}(t)^{-1} {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t).$$ Then by the fact that $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t} {{\widetilde{{\boldsymbol{U}}}}}(t) = i {{\boldsymbol{G}}}(t) {{\widetilde{{\boldsymbol{U}}}}}(t)$, the time-evolution equation for $\frac{{\,\mathrm{d}}}{{\,\mathrm{d}}t}{{\widetilde{{\boldsymbol{U}}}}}(t)$ could be derived easily.

From the second part of first order stationary condition ([*i.e.*]{}, with respect to ${{\boldsymbol{\eta}}}(t)$), one could deduce that ${{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {{\boldsymbol{C}}}{{\widetilde{{\boldsymbol{U}}}}}(t) = 0$. After plugging in the expression of ${{\boldsymbol{C}}}$, it follows immediately that $${{\boldsymbol{\eta}}}(t) = {{\widetilde{{\boldsymbol{U}}}}}^{\dagger}(t) {\mathcal{L}}\big({{\boldsymbol{\rho}}}_{{LR}}(t) \bigr){{\widetilde{{\boldsymbol{U}}}}}(t).$$

[[ Since there is only one solution satisfying stationary conditions, the optimal solution for pair $({{\boldsymbol{G}}}(t), {{\boldsymbol{\eta}}}(t))$ is unique.]{}]{}
