# Dynamical low-rank approximations of solutions to the Hamilton-Jacobi-
Bellman equation

Martin Eigel
Weierstrass Institute for
Applied Analysis and Stochastics
Berlin, Germany
<EMAIL_ADDRESS>&Reinhold Schneider
Department of Mathematics
Technical University Berlin
Berlin, Germany
<EMAIL_ADDRESS>&David Sommer
Weierstrass Institute for
Applied Analysis and Stochastics
Berlin, Germany
<EMAIL_ADDRESS>
###### Abstract

We present a novel method to approximate optimal feedback laws for nonlinear
optimal control based on low-rank tensor train (TT) decompositions. The
approach is based on the Dirac-Frenkel variational principle with the
modification that the optimisation uses an empirical risk. Compared to current
state-of-the-art TT methods, our approach exhibits a greatly reduced
computational burden while achieving comparable results. A rigorous
description of the numerical scheme and demonstrations of its performance are
provided.

_Keywords_ dynamical low-rank approximation $\cdot$ feedback control $\cdot$
Hamilton-Jacobi-Bellman $\cdot$ Variational Monte Carlo $\cdot$ tensor product
approximation

##  1 Introduction

Feedback control is ubiquitous and indispensable in real dynamical systems.
Since the controlled system can in general not be expected to follow model
predictions exactly, system trajectories will eventually leave the forecasted
path, meaning that any preplanned series of controls (albeit an optimal one)
is based on wrong assumptions and therefore not only suboptimal, but
potentially dangerous. As an illustration, one might think of an astronaut who
calculated an optimal course to land on the moon but then does not modify the
forecasted actuation values of their rocket-drive when atmospheric effects
steer them off said course, which will leave them drifting to outer space. It
is therefore vital to deploy controls based on current state feedback, where
current means as frequently as possible in practical applications.

However, The problem of computing an optimal feedback control law for
nonlinear optimal control problems is notoriously difficult. This is because
the synthesis of such a feedback law requires solving the Hamilton-Jacobi-
Bellman (HJB) equation, which is a nonlinear parabolic partial differential
equation (PDE) of generally high dimension $d\gg 1$ [BC97]. Classical schemes
to solve the HJB equation such as Galerkin-schemes in linear ansatz spaces
suffer from the curse of dimensionality [KK18], i.e. an exponential complexity
growth. In practice, this means that the computation of a solution is often
infeasibly slow if it can be discretised and stored at all. Another severe
obstacle can be the low regularity of viscosity solutions, cf [BD+97]. In this
paper, our focus lies on the alleviation of the curse of dimensionality in
order to enable the numerical treatment of high-dimensional control problems.
We hence only consider problems where the lack of regularity is not present or
not pronounced enough to prevent a sufficiently accurate approximation.

The relevance of efficient numerical methods can be seen by the fact that true
feedback control methods - that is: methods solving the HJB equation - are
rarely used in practice due to the necessary computational effort. Control
problems arising e.g. in mechanical engineering often require new planning of
controls to be computed within seconds. There hence is a tight upper limit on
the time budget available for generating new controls. Therefore, most
engineers deploy variations of Model Predictive Control (MPC) where open-loop
controls are computed in such rapid succession that they effectively “close
the loop” [CA13]. This is a conservative approach since the feedback property
of the resulting controller is obtained solely by means of the measurements at
the discrete planning steps. In between two state measurements, the controller
is not in feedback form.

In this work, we present a novel method to tackle nonlinear optimal control
problems that $(1)$ yields a true feedback controller and $(2)$ has greatly
reduced computational cost compared to current state-of-the-art methods. Our
method is based on policy iteration, linearising the HJB equation (which is
then sometimes called the generalised Hamilton-Jacobi-Bellman or GHJB
equation) and a modification of the Dirac-Frenkel variational principle. This
then allows the computation of approximate solutions on a specified function
manifold, for which we choose the set of multivariate polynomials with a fixed
tensor train (TT) rank.

Tree based tensor networks and tensor trains in particular have already been
used for successful approximations of the value function in various works, see
e.g. [OSS21, KKD19, Fac+20]. These recent results are summarized in the PhD
thesis of Leon Sallandt [Sal21], which is still being finalised as this paper
is written. There, the approach is based on a Lagrangian (or dynamic
programming) perspective by computing the value functions for several initial
states and learning the global function from these values by regression using
a multi-polynomial TT model. With appropriate modifications, this approach can
already be combined with regression techniques performed e.g. by machine
learning methods, in particular artificial neural networks (NN). In the
present paper, we follow a different approach, exploiting the Riemannian
structure of the TT manifold [HRS12, Ste16] by an empirical version of the
Dirac-Frenkel principle.

The solution obtained by the abstract Dirac-Frenkel principle can be shown to
be quasi-optimal in some time interval $[0,T_{\text{DF}}]$ but deteriorates
from the best low-rank approximation after a certain time [Lub+13]. We expect
a similar behaviour in our case which may restrict the time interval in
practice. Combining both approaches - abstract and empirical - is an open
research question, which we aim to address in future work. Similarly, we defer
the stochastic control case to a forthcoming paper, confining ourselves to
deterministic control in this paper. We conjecture that the present approach
is even more advantageous in the stochastic case.

The rest of the paper is organised as follows: In Section 2 we provide a short
overview of the related literature, specifically the current state-of-the-art
of tensor based methods to solve the HJB equation. Section 3 introduces the
finite horizon optimal control problem in feedback form, which the rest of
this work revolves around. In Section 4 the tensor train format, the
corresponding manifold and the representation of the tangent space are
introduced. These are needed to formulate the Dirac-Frenkel variational
principle, which is introduced in its abstract form in Section 5. In Section
6, we combine the concepts of Sections 4 and 5 to develop our proposed DLRA
method for approximately solving the HJB equation. Numerical results that
illustrate the practical performance are presented and discussed in Section 7.
Finally, we close in Section 8 with an outlook on future work.

##  2 Related work

The Bellman equation governing the value function of an optimal control (OC)
problem was introduced as early as 1957 by Richard Bellman [Bel57]. Since
then, numerous sophisticated methods have been introduced to approximate
solutions, mostly based on the principle of dynamical programming, see e.g.
[Ber05] for a broad introduction to the subject. The alternative approach,
which we follow in this work, is to consider the infinitesimal version of the
Bellman equation, namely the Hamilton-Jacobi-Bellman equation [BC97], which is
a nonlinear parabolic PDE. In both cases, many methods rely on a fixed point
iteration of the equation, which in the OC and Reinforcement Learning (RL)
literature is called policy iteration [How60]. Alternatives are domain
splitting algorithms [FLS94], semi-Lagrangian methods [Fal87, FK14, TAK17],
data-based methods using Neural Networks [Luo+14], variational iterative
methods [KDK13], actor-critic methods [ZH21], tree-based methods [AS19] and
tropical methods [AGL08, AF18].

For a fixed starting value, an optimal control can be obtained by open-loop
approaches such as Pontryagin’s maximum principle [BGP61, Pon87]. In this way,
the value function can be evaluated pointwise by simply adding up the cost of
that optimal control. Controls of this type have been used to find the value
function e.g. in [KW17, NGK19, AKK21, OSS21]. In this work, we use optimal
open-loop controls as benchmarks to which we compare the feedback controller
computed by our method.

Since any solution method for the HJB equation has to deal with the curse of
dimensionality, some form of model order reduction has to take place in
practical applications. Possible function approximations can be obtained by
using neural networks [DLM19, NR21, IRZ21] or sparse polynomials [AKK21]. In
this work we use the TT format introduced to the mathematical community by
Oseledets [OT09, Ose11] for multivariate polynomials. A striking recent
example of the power of the low rank TT structure for function approximation
can be found in [RSN21], in which the authors use TTs with polynomial basis
functions to outperform state-of-the-art NNs on the solution of parabolic PDEs
by orders of magnitude, while requiring lower computational time. For further
details on TTs and more general hierarchical tensor networks, we refer the
reader to the survey articles [HS14, BSU16] and the standard textbooks [Hac12,
Hac14]. For recent applications of TTs as value function approximators, see
e.g. [OSS21, KKD19]. Solution methods based on high-dimensional polynomials
and tensor spaces have also been considered in [KK18, DKK21]. As a conjecture
for future work, block sparsity of the TTs appearing in optimal control
methods could be exploited to further reduce the sample complexity [TSG21].

In contrast to the aforementioned methods, our new approach is a dynamical low
rank approximation (DLRA) [KL07, KL10] of the value function. The main idea is
to approximate solutions to matrix- or tensor-valued ordinary differential
equations (ODEs) by projecting the right-hand side onto the tangent space of
the manifold of matrices/tensors of fixed (TT-)rank at the current
approximation. In this abstract setting, the projection is usually decomposed
into orthogonal parts of the tangent space after which a splitting scheme is
applied, resulting in so called projector-splitting schemes [Lub+13, LOV15,
KLW16, CL20, CKL21]. The obtained approximation is quasi-optimal on a finite
time domain, a property known as the Dirac-Frenkel variational principle, or
Dirac-Frenkel/McLachlan variational principle [Mur35, McL64]. DLR
approximations to parabolic PDEs have been studied in [Con20, Bac+21], but –
to the best of our knowledge – this work is the first application of DLR
methods to a finite horizon optimal control problem and in particular to the
nonlinear HJB equation. In order to derive an abstract DLR problem on the TT
manifold, we use a Variational Monte Carlo (VMC) approach [EST20, Bay+21]. In
our setting it can be understood as an empirical least squares tensor
regression based on random samples.

##  3 The optimal control problem

Throughout this work, we consider a deterministic dynamical system

$\displaystyle\dot{x}(t)$ $\displaystyle=f(t,x(t))+g(t,x(t))u(t),\quad
t\in[t_{0},T],$ (1) $\displaystyle x(t_{0})$ $\displaystyle=x_{0},$ (2)

with initial time $t_{0}\in[0,T]$, initial condition
$x_{0}\in\Omega\subset\mathbb{R}^{d}$, control $u\in
L^{2}(0,T;\mathbb{R}^{m})$, free dynamics
$f:[0,T]\times\Omega\rightarrow\Omega$ and control interface
$g:[0,T]\times\Omega\rightarrow\mathbb{R}^{d\times m}$. To ensure existence
and uniqueness of solutions (for admissible controls $u$), we assume $f$ and
$g$ to be smooth (possibly nonlinear) functions. A total cost is associated
with the triple $(t_{0},x_{0},u)\in[0,T]\times\Omega\times
L^{2}(0,T;\mathbb{R}^{m})$ in terms of the cost functional

$\displaystyle\mathcal{J}(t_{0},x_{0},u)=\int_{t_{0}}^{T}c(t,x(t))+u(t)^{\intercal}R(t)u(t)\mathrm{d}t+c_{T}(x(T)),$
(3)

where the running cost $c:[0,T]\times\Omega\rightarrow\mathbb{R}_{+}$ and the
terminal cost $c_{T}:\Omega\rightarrow\mathbb{R}_{+}$ are non-negative,
coercive and smooth functionals. Moreover,
$R:[0,T]\rightarrow\mathbb{R}^{m\times m}$ is continuous, $R(t)$ is positive
definite for all $t$, and the trajectory $x(\cdot)$ is subject to (1)+(2) with
the given control $u$. The function mapping time-state pairs to optimal future
costs is called the value function. It is canonically defined as

$\displaystyle
V^{*}:[0,T]\times\Omega\rightarrow\mathbb{R},\quad(t_{0},x_{0})\mapsto\inf_{u\in
L^{2}(0,T;\mathbb{R}^{m})}\mathcal{J}(t_{0},x_{0},u).$

If the dynamics and cost terms satisfy sufficient regularity conditions, the
value function is given as the viscosity solution of the well known Hamilton-
Jacobi-Bellman equation.

######  Theorem 1 (see e.g. [BC97, BD+97]).

Let $\ell(t,x,u)=c(t,x)+u^{\intercal}R(t)u$ and assume there are
$\sigma,\delta\geq 1$ with $\sigma<\delta$, $\ell_{0}>0$. Moreover, for every
compact $K\subset\mathbb{R}^{d}$ there exists some $f_{K}>0$ such that

$\displaystyle\|f(x,u)\|$ $\displaystyle\leq
f_{K}(1+\|x\|^{\sigma}),\quad\textnormal{ for all }~{}(x,u)\in
K\times\mathbb{R}^{m},$ $\displaystyle|\ell(x,u)|$
$\displaystyle\geq\ell_{0}\|a\|^{\delta}\quad\textnormal{ for all
}~{}(x,u)\in\mathbb{R}^{d}\times\mathbb{R}^{m}.$

Then, the value function $V^{*}$ is the unique viscosity solution of the HJB
equation

$\displaystyle\dfrac{\partial}{\partial
t}V^{*}(t,x)+\min_{u\in\mathbb{R}^{m}}\left[\nabla_{x}V^{*}(t,x)^{\intercal}(f(t,x)+g(t,x)u)+c(t,x)+u^{\intercal}R(t)u\right]$
$\displaystyle=0$ (4) $\displaystyle V^{*}(T,\cdot)$
$\displaystyle=c_{T}(\cdot).$ (5)

Note that the HJB equation is an infinitesimal version of the Bellman
equation, which we state for the sake of completeness.

######  Theorem 2 ([BC97, BD+97]).

For all $x_{0}\in\Omega$ and $0\leq t_{0}\leq t_{1}\leq T$, we have

$\displaystyle V^{*}(t_{0},x_{0})=\inf_{u\in
L^{2}(t_{0},t_{1};\mathbb{R}^{m})}\left[\int_{t_{0}}^{t_{1}}\ell(t,x(t),u(t))\mathrm{d}t+V^{*}(t_{1},x(t_{1}))\right],$
(6)

where $x(t)$ satisfies (1) with initial condition $x(t_{0})=x_{0}$ and control
$u$.

Now consider feedback controls of the form $u(t)=\alpha(t,x(t))$, where
$\alpha:[0,T]\times\Omega\rightarrow\mathbb{R}^{m}$ is continuous on $[0,T]$
and Lipschitz in $\Omega$. We call such functions $\alpha$ admissible feedback
laws (or equivalently admissible policies) and we denote the set of admissible
policies by $\mathcal{A}$. Next, we define the policy evaluation function
$\mathcal{J}_{\alpha}$ via the associated cost

$\displaystyle\mathcal{J}^{\alpha}(t_{0},x_{0})=\int_{t_{0}}^{T}c(t,x(t))+\alpha(t,x(t))^{\intercal}R(t)\alpha(t,x(t))\mathrm{d}t+c_{T}(x(T)).$
(7)

An optimal policy $\alpha^{*}$ is a policy which achieves minimal costs for
any starting values, i.e.

$\displaystyle\mathcal{J}^{\alpha^{*}}(t_{0},x_{0})=\min_{\alpha\in\mathcal{A}}\mathcal{J}^{\alpha}(t_{0},x_{0}),\quad\textnormal{
for all }~{}(t_{0},x_{0})\in[0,T]\times\Omega.$

The goal of optimal (feedback) control is to approximate such an optimal
policy. If the value function is known and partially differentiable, an
optimal policy can be obtained immediately, as the following theorem states.

######  Theorem 3 ([BC97]).

An optimal policy is given by

$\displaystyle\alpha^{*}(t,x)=-\dfrac{1}{2}R(t)^{-1}g(t,x)^{\intercal}\nabla_{x}V^{*}(t,x),$
(8)

if the gradient of $V^{*}$ exists.

Hence, the problem of synthesizing an optimal policy is the problem of finding
the value function, which involves solving the HJB equation (4) or the Bellman
equation (6).

In the following, we always assume that the conditions of Theorem 1 are
satisfied so that that an optimal policy is given by the value function via
(8). Hence, we can identify the value function with the policy evaluation
function of that optimal policy denoted by $V^{*}=\mathcal{J}^{\alpha^{*}}$.
Our goal is to approximate the value function successively on small
subintervals, moving backwards in time from $t=T$ to $t=0$. This approach is
based on Bellman’s principle. However, in contrast to comparable recent work
[OSS21], we use the HJB equation (4) on each subinterval instead of the
Bellman equation. In particular, we define suitable approximate solutions to
the HJB equation by means of the Dirac-Frenkel variational principle. While
theoretical simplicity is lost to some extend, computational simplicity is
gained in return. This is mainly because DLR approximations of (4) can be
computed very efficiently since samples do not have to be propagated through
the dynamics to evaluate the integral in (6).

Assume now that $0=t_{0}<t_{1}<\ldots t_{m}=T$ is an equidistant
discretisation and consider a partitioning $\\{[t_{i},t_{i+1}]\\}_{i=0}^{m}$
of the time interval $[0,T]$. An immediate consequence of Bellman’s principle
is that an optimal policy for the whole time domain $[0,T]$ must also be
optimal on any subinterval $[t_{i},t_{i+1}]$. Conversely, a policy that is
optimal on all subintervals is also optimal on the whole interval. This
enables to learn the value function by moving backwards in time and
(approximately) computing the restrictions
$V^{*}(t,x)\big{|}_{t\in[t_{i},t_{i+1}]}$ for $i=m-1,\ldots,0$. In the
following, we denote by $V^{*}_{i}=V^{*}\big{|}_{[t_{i},t_{i+1}]}$ for
$i=0,\ldots,m-1$ the restrictions of the value function to a particular
subinterval and set $V^{*}_{m}=c_{T}$. Approximations of $V^{*}_{i}$ are
denoted by $\hat{V}_{i}$ and the approximation $\hat{V}$ of $V^{*}$ on the
whole time domain is defined by $\hat{V}=\hat{V}_{i}$ on $[t_{i},t_{i+1}]$.
Algorithm 1 summarizes the idea of successive backward approximation, which we
deploy to approximate the value function.

Data: Time discretisation points $0=t_{0}<\ldots<t_{m}=T$, approximation
$\hat{V}_{m}$ of the terminal cost.

Result: Approximation $\hat{V}$ of the value function.

for  _$i=m-1,m-2,\ldots,0$_ do

Compute approximate solution $\hat{V}_{i}$ of the HJB-eq. (4) on
$[t_{i},t_{i+1}]$ with terminal condition $\hat{V}_{i+1}$.

Set $\hat{V}=\hat{V}_{i}$ on $[t_{i},t_{i+1}]$.

end for

Algorithm 1 Bellman-based backwards scheme to approximate the value function

TT approximations of the value function by means of such a backwards scheme
were already presented e.g. in [OSS21]. In that work however, the integral
formulation (6) is used exclusively, sampling trajectories $x(t)$ for given
controls and adding up the costs. In contrast, the DLR approximation method
used here allows to directly work with the HJB equation (4).

##  4 Tensor trains as function approximators

For practical computations, the approximations $\hat{V}_{i}$ from Algorithm 1
have to be confined to a finite-dimensional functions space. To this end,
consider a set of one-dimensional basis functions
$\phi_{1},\ldots,\phi_{n}:\mathbb{R}\rightarrow\mathbb{R}$ and functions
$v:\mathbb{R}^{d}\rightarrow\mathbb{R}$ of the form

$\displaystyle
v(x)=A\phi(x)=\sum_{i_{1},\ldots,i_{d}=1}^{n}A_{i_{1},\ldots,i_{d}}\phi_{i_{1}}(x_{1})\cdot\ldots\cdot\phi_{i_{d}}(x_{d}),$
(9)

with coefficient tensor $A\in\mathbb{R}^{n\times n\times\ldots\times n}$ of
order $d$. Usually, the basis functions $\phi_{i}$ are (orthonormal)
polynomials. Consequently, $v$ is a multivariate polynomial with a storage
complexity of $\mathcal{O}(n^{d})$ for its coefficient tensor. The TT format
provides a possibility to alleviate this exponential complexity by assuming
some low-rank structure. A TT representation of $A$ is any decomposition of
the form

$\displaystyle A_{i_{1},\ldots,i_{d}}=U^{1}_{i_{1}}\cdot\ldots\cdot
U^{d}_{i_{d}},$ (10)

where

$\displaystyle U^{1}\in\mathbb{R}^{n\times
r_{1}},~{}U^{\mu}\in\mathbb{R}^{r_{\mu-1}\times n\times
r_{\mu}}~{}\textnormal{for}~{}i=2,\ldots,d-1,~{}U^{d}\in\mathbb{R}^{r_{d-1}\times
n}$

are called the components of the representation and $i_{\mu}$ denotes the
middle index of the component, i.e.
$U^{\mu}_{i_{\mu}}\in\mathbb{R}^{r_{\mu-1}\times r_{\mu}}$. The rank of the
specific representation is given by the tuple $(r_{1},\ldots,r_{d-1})$. The
TT-rank $\mathbf{r}=(r_{1},\ldots,r_{d-1})$ of $A$ is defined as the (entry-
wise) minimal rank tuple such that a TT representation (10) with the
corresponding ranks exists. Such a minimal TT representation exists for any
tensor. In fact, the minimal rank entry $r_{\mu}$ is equal to the matrix rank
of the $\mu$-th unfolding of $A$ (for details we refer to [HRS12a]). The TT
representation exhibits a storage complexity of
$\mathcal{O}(nd\max(r_{1},\dots,r_{d-1})^{2})$, scaling only linearly in the
dimension $d$, and hence avoiding the curse of dimensionality, provided that
the ranks stay bounded. It is important to note that even for fixed rank
$\mathbf{r}$, a decomposition of the form (10) is not unique. For any
$\mu=1,\ldots,d-1$ we can set $U_{\mu}\rightarrow U_{\mu}S$ and
$U_{\mu+1}\rightarrow S^{-1}U_{\mu+1}$ for invertible
$S\in\mathbb{R}^{r_{\mu}\times r_{\mu}}$ without changing the tensor. A unique
representation is then given by requiring left- and right-orthogonality of the
components in the sense of the following definition.

######  Definition 1.

For a component $U_{\mu}\in\mathbb{R}^{r_{\mu-1}\times n\times r_{\mu}}$,
define the left and right unfolding

$\displaystyle L(U^{\mu})\in\mathbb{R}^{r_{\mu-1}n\times r_{\mu}},\qquad
R(U^{\mu})\in\mathbb{R}^{r_{\mu-1}\times r_{\mu}n},$

by suitable matrix reshaping (for details regarding the order, see e.g.
[Ste16]). A component $U^{\mu}$ is called left- or right-orthogonal if

$\displaystyle
L(U^{\mu})^{\intercal}L(U^{\mu})=I_{d}\in\mathbb{R}^{r_{\mu}\times
r_{\mu}},\quad\textnormal{or}\quad
R(U^{\mu})R(U^{\mu})^{\intercal}=I_{d}\in\mathbb{R}^{r_{\mu-1}\times
r_{\mu-1}},$

respectively. A TT representation
$A_{i_{1},\ldots,i_{d}}=U^{1}_{i_{1}}\cdot\ldots\cdot U^{d}_{i_{d}}$ of a
tensor $A$ is called $\mu$-orthogonal if $U^{1},\ldots,U^{\mu-1}$ are left
orthogonal and $U^{\mu+1},\ldots,U^{d}$ are right orthogonal. In that case,
$U^{\mu}$ is called the core of the representation.

Left and right orthogonality of all but one component imposes
$\sum_{\mu=1}^{d-1}r_{\mu}^{2}$ additional conditions on the representation.
Hence, the $\mu$-orthogonal TT representation of $A$ is unique for any $\mu$.

For a given TT rank $\mathbf{r}$, we define the set

$\displaystyle\mathcal{M}_{\mathbf{r}}=\\{A\in\mathbb{R}^{n\times\ldots\times
n}:A~{}\textnormal{ has TT rank }\mathbf{r}\\}.$

It is noteworthy that $\mathcal{M}_{\mathbf{r}}$ is a smooth manifold in
$\mathbb{R}^{n\times\ldots\times n}$ [HRS12a]. With a chosen suitable basis
$\\{\phi_{1},\ldots,\phi_{n}\\}$, we define a set of function approximations

$\displaystyle\mathfrak{F}_{\mathbf{r}}=\left\\{v:\mathbb{R}^{d}\rightarrow\mathbb{R}:v~{}\textnormal{admits
a
representation}~{}\eqref{eq:tensor_func}~{}\textnormal{with}~{}A\in\mathcal{M}_{\mathbf{r}}\right\\}.$

Note that by identification of a function with its coefficient tensor,
$\mathfrak{F}_{\mathbf{r}}$ forms a smooth manifold in the $n^{d}$-dimensional
linear space $\mathfrak{F}=\bigcup_{\mathbf{r}}\mathfrak{F}_{\mathbf{r}}$ in
the same way that $\mathcal{M}_{\mathbf{r}}$ forms a smooth manifold in
$\mathbb{R}^{n\times\ldots\times n}$. In order to do perform an optimisation
on $\mathfrak{F}_{\mathbf{r}}$, or $\mathcal{M}_{\mathbf{r}}$, respectively,
we require a representation of the tangent space
$\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})$ of $\mathcal{M}_{\mathbf{r}}$ in
$U$. Throughout this work, we use the following representation.

######  Theorem 4 ([HRS12a] or [Ste16]).

Let $U\in\mathcal{M}_{\mathbf{r}}$ be $d$-orthogonal. The tangent space
$\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})$ of $\mathcal{M}_{\mathbf{r}}$ in
the point $U$ is given by $\tau(X)$, where

$\displaystyle X$ $\displaystyle=U^{\ell}_{1}\times\ldots\times
U^{\ell}_{d-1}\times\mathbb{R}^{r_{d-1}\times n\times r_{d}},$ $\displaystyle
U_{\mu}^{\ell}$ $\displaystyle=\\{W^{\mu}\in\mathbb{R}^{r_{\mu-1}\times
n\times
r_{\mu}}:L(U^{\mu})^{\intercal}L(W^{\mu})=0\in\mathbb{R}^{r_{\mu}\times
r_{\mu}}\\},$

and

$\displaystyle\tau:X\longrightarrow\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}}),\quad\tau(W^{1},\ldots,W^{d})=\delta
U$ $\displaystyle\delta
U_{i_{1},\ldots,i_{d}}=\sum_{\mu=1}^{d}U^{1}_{i_{1}}\cdot\ldots\cdot
U^{\mu-1}_{i_{\mu-1}}W^{\mu}_{i_{\mu}}U^{\mu+1}_{i_{\mu+1}}\cdot\ldots\cdot
U^{d}_{i_{d}}.$ (11)

The tangent space has the same dimension as the underlying manifold. The
previously mentioned ambiguity in the representation is now eliminated due to
the gauging conditions $L(U^{\mu})^{\intercal}L(W^{\mu})=0$ in
$U^{\ell}_{\mu}$.

######  Corollary 1.

Each of the spaces $U_{\mu}^{\ell}$ has dimension
$r_{\mu-1}nr_{\mu}-r_{\mu}^{2}$ and hence the tangent space has dimension

$\displaystyle
n_{X}\coloneqq\textnormal{dim}(X)=\sum_{\mu=1}^{d}r_{\mu-1}nr_{\mu}-\sum_{\mu=1}^{d-1}r_{\mu}^{2}.$

Using the representation (11) for elements $\delta U$ of the tangent space of
$U$, a simple form for the sum $U+\delta U$ can be obtained.

######  Lemma 1 (see [Ste16]).

Let $U\in\mathcal{M}_{\mathbf{r}}$ be $d$-orthogonal and denote its component
tensors by $U_{1},\ldots,U_{d}$. Let $\delta
U\in\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})$ be given by $(\delta
U_{1},\ldots,\delta U_{d})\in X$. Then,

$\displaystyle U(i_{1},\ldots,i_{d})+\delta U(i_{1},\ldots,i_{d})=$
$\displaystyle\begin{bmatrix}\delta
U_{1}(i_{1})&U_{1}(i_{1})\end{bmatrix}\begin{bmatrix}U_{2}(i_{2})&0\\\ \delta
U_{2}(i_{2})&U_{2}(i_{2})\end{bmatrix}\ldots$
$\displaystyle\ldots\begin{bmatrix}U_{d-1}(i_{d-1})&0\\\ \delta
U_{d-1}(i_{d-1})&U_{d-1}(i_{d-1})\end{bmatrix}\begin{bmatrix}U_{d}(i_{d})\\\
U_{d}(i_{d})+\delta U_{d}(i_{d})\end{bmatrix}.$

This can easily be verified by multiplying out the matrix products. In
particular, the sum $U+\delta U$ has at most TT-rank $2\mathbf{r}$.

##  5 The Dirac-Frenkel variational principle

The Dirac-Frenkel variational principle [Mur35] provides a principled way to
approximate tensor valued ODEs of the form

$\displaystyle\dot{A}(t)$ $\displaystyle=F(t,A(t)),$ (12) $\displaystyle A(0)$
$\displaystyle=A_{0},$ (13)

where $A(t)\in\mathbb{R}^{n\times\ldots\times n}$, on the manifold
$\mathcal{M}_{\mathbf{r}}$. More precisely, given an approximation
$Y_{0}\in\mathcal{M}_{\mathbf{r}}$ of the initial condition $A_{0}$, an
approximation $Y(t)\in\mathcal{M}_{\mathbf{r}}$ of $A(t)$ is defined as the
solution of the TT-valued ODE

$\displaystyle\dot{Y}(t)$
$\displaystyle=\operatorname*{arg\,min}_{\vartheta\in\mathcal{T}_{\mathcal{M}_{\mathbf{r}}}(Y(t))}\|\vartheta-F(t,Y(t))\|,$
(14) $\displaystyle Y(0)$ $\displaystyle=Y_{0}.$ (15)

The minimum in (14) is attained by the orthogonal projection of the right-hand
side onto the tangent space, leading to

$\displaystyle\dot{Y}(t)=P_{\mathcal{T}_{\mathcal{M}_{r}}(Y(t))}F(t,Y(t)).$
(16)

In this abstract setting, error bounds can be derived, which we quote for the
sake of completeness.

######  Theorem 5.

[Lub+13] Suppose that $\dot{A}(t)\leq\mu$ and that a continuously
differentiable best approximation $X(t)\in\mathcal{M}_{\mathbf{r}}$ to $A(t)$
exists for $t\in[0,T]$. Let $\delta>0$ be such that the smallest nonzero
singular value of every matrix unfolding of $X(t)$ is greater or equal to
$\rho$, and assume that the best-approximation error is bounded by
$\|X(t)-A(t)\|\leq c\rho$ for $t\in[0,T]$ with a constant $c$ depending only
on the dimension $d$. Then, the approximation error of the dynamical low-rank
approximation defined by (20) with initial value $Y(0)=X(0)$ is bounded by

$\displaystyle\|Y(t)-X(t)\|\leq 2\beta e^{\beta
t}\int_{0}^{t}\|X(s)-A(s)\|\mathrm{d}s,$

with $\beta=C\mu\rho-1$ for $t\in[0,T]$, as long as the right-hand side
remains bounded by $c\rho$. The constant $C$ is only dependent on $d$ and is
given in [Lub+13].

In recent years there have been numerous works on the numerical treatment of
ODEs of this type, see [KL07, LO13, KLW16] for an introduction in the matrix
case and [LOV15, CL20, CKL21] for more recent tensor-based research
directions. Generally, these methods rely on a splitting of the projector
$P_{\mathcal{M}_{\mathbf{r}}(Y(t))}$ into orthogonal parts of the tangent
space, so-called projector splitting algorithms. The norm $\|.\|$ governing
(14) and hence the projector is usually the Frobenius norm. This is in
contrast to our work, where $\|.\|$ is an empirical norm111the details of
which are provided in the next chapter. Carrying over results from the
treatment of the abstract Dirac-Frenkel principle to the empirical case
(specifically the projector splitting schemes) is an important direction of
future work, that we do not yet address in this paper.

##  6 Dynamical low-rank approximation of the HJB equation

Based on the preceding review of tools that we require, we now return to the
HJB equation (4) on $[t_{i},t_{i+1}]$ with terminal condition
$\hat{V}_{i+1}(t_{i+1},\cdot)$. The goal is to obtain an approximation
$\hat{V}_{i}$ of the value function on the current interval. Inserting (8)
into the HJB (4) leads to a coupled problem:
Find $V$ such that

$\displaystyle\dfrac{\partial}{\partial
t}V(t,x)+\nabla_{x}V(t,x)^{\intercal}(f(t,x)+g(t,x)\alpha(t,x))+c(t,x)+\alpha(t,x)^{\intercal}R(t)\alpha(t,x)$
$\displaystyle=0,$ (17) $\displaystyle V(t_{i+1},\cdot)$
$\displaystyle=\hat{V}_{i+1}(t_{i+1},\cdot),$ (18)

where $\alpha$ satisfies

$\displaystyle\alpha(t,x)=-\dfrac{1}{2}R(t)^{-1}g(t,x)^{\intercal}\nabla_{x}V(t,x).$
(19)

To compute $\hat{V}_{i}$, we use a fixed point iteration of the coupled
problem, iteratively solving (17)+(18) for fixed $\alpha$ and then updating
$\alpha$ via (19). This procedure is known as policy iteration in the optimal
control literature. We depict a conceptual summary in algorithm 2. If the
solutions to (17) are exact, it converges under mild assumptions on dynamics
and cost terms [SL79]. In order to track the convergence of the scheme under
approximations, we introduce on $L^{2}((t_{i},t_{i+1});L^{2}(\Omega,\rho))$
and $L^{2}((0,T);L^{2}(\Omega,\rho))$ the norms

$\displaystyle\|v\|_{i}^{2}=\int_{t_{i}}^{t_{i+1}}\|v(t,\cdot)\|^{2}_{L^{2}(\Omega,\rho)}\mathrm{d}t,\qquad\|v\|^{2}=\sum_{i=0}^{m-1}\|v\|^{2}_{i},$

and stop the iteration once the $\|.\|_{i}$-difference of two consecutive
approximations becomes smaller than a specified threshold.

Data: Interval $[t_{i},t_{i+1}]$, terminal condition $\hat{V}_{i+1}$,
admissible policy $\alpha$, error tolerance $\delta$.

Result: Approximation $\hat{V}_{i}$.

while  _norm change of $\hat{V}_{i}>\delta$_ do

Compute approximate solution $\hat{V}_{i}$ of (17)+(18) on $[t_{i},t_{i+1}]$.

Update $\alpha\propto\nabla_{x}\hat{V}_{i}$ according to (19).

end while

Algorithm 2 Policy iteration on subinterval

It remains to be shown how to compute the approximations $\hat{V}_{i}$. To
ease notation, without loss of generality we consider the interval
$[t_{0},t_{1}]$ instead of $[t_{i},t_{i+1}]$ for the remainder of this
chapter. We construct $\hat{V}_{0}$ as a dynamical low-rank approximation of
(17) in the tensor train format.

Let the terminal condition $\hat{V}_{1}\in\mathfrak{F}_{\mathbf{r}}$ and
consider for given $\alpha$ the following problem:
Find $V$ such that

$\displaystyle\dfrac{\partial}{\partial t}V(t,\cdot)$
$\displaystyle=\operatorname*{arg\,min}_{\vartheta\in\mathcal{T}_{V(t,\cdot)}(\mathfrak{F}_{\mathbf{r}})}\|\vartheta-\tilde{F}(t,V(t,\cdot))\|^{2}_{L^{2}(\Omega,\rho)},$
(20) $\displaystyle V(t_{1},\cdot)$ $\displaystyle=\hat{V}_{1}(t_{1},\cdot),$
(21)

where
$\tilde{F}(t,V(t,\cdot))(\cdot)=-\nabla_{x}V(t,\cdot)^{\intercal}(f(t,\cdot)+g(t,\cdot)\alpha(t,\cdot))-c(t,\cdot)+\alpha(t,\cdot)^{\intercal}R(t)\alpha(t\cdot)$.
Note that this essentially means that the time derivative of $V_{0}$ is
approximated in the tangent space of the current solution. By a simple time
inversion $t\rightarrow t_{0}+(t_{1}-t)$, the terminal condition can be turned
into an initial condition. Crucially, any solution to (20) stays on the
manifold $\mathfrak{F}_{\mathbf{r}}$ and can therefore be identified with a
time-dependent coefficient tensor $A(t)\in\mathcal{M}_{\mathbf{r}}$ via
$V(t,x)=A(t)\phi(x)$. Denoting the coefficient tensor of $\hat{V}_{1}$ by
$\hat{A}_{1}$, we see that the abstract problem (20)+(21) is equivalent to the
TT-valued ODE

$\displaystyle\dot{A}(t)$
$\displaystyle=\operatorname*{arg\,min}_{B\in\mathcal{T}_{A(t)}(\mathcal{M}_{\mathbf{r}})}\|B\phi-F(t,A(t)\phi)\|^{2}_{L^{2}(\Omega,\rho)},$
(22) $\displaystyle A(t_{0})$ $\displaystyle=A_{1},$ (23)

where $F$ arises from time inversion of $\tilde{F}$.

In general, the $L^{2}$-integral on the right-hand side of (22) is difficult
to compute. Nevertheless, we can easily carry out a pointwise evaluation of
the basis functions $\phi$ as well as the other terms in $F(t,A(t)\phi)$. In
practice, we hence replace the exact $L^{2}$-norm with a Monte Carlo
approximation

$\displaystyle\|v\|^{2}_{L^{2}(\Omega,\rho,M)}=\dfrac{1}{M}\sum_{k=1}^{M}|v(x_{k})|^{2},\qquad
x_{k}\sim\rho,$

for $v\in L^{2}(\Omega,\rho)$. This turns the right-hand side of the ODE into
an empirical risk minimisation. We eventually arrive at

$\displaystyle\dot{A}(t)$
$\displaystyle=\operatorname*{arg\,min}_{B\in\mathcal{T}_{A(t)}(\mathcal{M}_{\mathbf{r}})}\dfrac{1}{M}\sum_{k=1}^{M}|B\phi(x_{k})-F(t,A(t)\phi(x_{k}))(x_{k})|^{2},\qquad
t\in(t_{0},t_{1})$ (24) $\displaystyle A(t_{0})$ $\displaystyle=\hat{A}_{1}.$
(25)

Statistical bounds for the error of the empirical minimiser in (24) compared
to the best $L^{2}$-approximation
$\Phi^{*}(t)=\operatorname*{arg\,min}_{\Phi\in
L^{2}(\Omega,\rho)}\|\Phi-F(t,A(t)\phi(\cdot))(\cdot)\|_{L^{2}(\Omega,\rho)}^{2}$
are given in [Eig+19].

A crucial observation is that the minimisation on the right-hand side is a
linear problem since the optimisation is over the linear tangent space.
Implementation details on how the minimum in (24) for a given $t$ can be
computed are given in Appendix A. Since the fit is linear, issues of local
minima are avoided which for instance occur in the alternating linear scheme
(ALS) [HRS12a] and other nonlinear optimisation methods. Alternating methods
can still be applied here to divide the problem into smaller sub-problems and
reduce the computational burden, leading (in their simplest form) to an
effective Lie-Trotter type splitting of the right hand side. A more detailed
examination of this topic is however beyond the scope of this paper and might
be addressed in future work.

The numerical realisation of (24) poses an additional hurdle. While the true
solution $A(t)$ always stays on the manifold $\mathcal{M}_{\mathbf{r}}$, it is
straightforward to see that any one step with a numerical integrator, e.g. a
Runge-Kutta method, leads to leaving it. This is due to the fact that by Lemma
1 any sum $U+\delta U$ where $U\in\mathcal{M}_{\mathbf{r}}$ and $\delta
U\in\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})$ has rank $2\mathbf{r}$ in
general. We therefore need to retract back onto the manifold after each step
of the integrator by truncating the ranks appropriately. To make this precise,
let $t_{0}=t^{(0)}<t^{(1)}<t^{(2)}<\ldots<t^{(L)}=t_{1}$ be a micro-
discretisation of the macro-interval $[t_{0},t_{1}]$ with equidistant step
size $\tau$ and define a numerical approximation $A_{\ell}$ of $A(t^{(\ell)})$
by the explicit Euler scheme

$\displaystyle A_{0}$ $\displaystyle=\hat{A}_{1},$ $\displaystyle A_{\ell+1}$
$\displaystyle=\mathcal{R}(A_{\ell}+\tau\Delta
A_{\ell}),\qquad\ell=0\ldots,L-1.$

Here, $\Delta A_{\ell}$ is the solution to the minimisation problem on the
right-hand side of (24) if $A_{\ell}$ is substituted for $A(t)$, the addition
$A_{\ell}+\tau\Delta A_{\ell}$ is performed like in Lemma 1, and $\mathcal{R}$
denotes the rank-truncation of a TT with rank $2\mathbf{r}$ back to a tensor
of rank $\mathbf{r}$. This truncation is performed by a TT-SVD with fixed rank
[OT09]. Once all $A_{\ell}$ are obtained in this way, we define
$\hat{V}_{0}(t,x)$ by linear interpolation, i.e.

$\displaystyle\hat{V}_{0}(t,x)=A_{\ell}\phi(x)+\dfrac{t-t^{(\ell)}}{\tau}(A_{\ell+1}-A_{\ell})\phi(x)\quad\textnormal{for}\quad
t\in[t^{(\ell)},t^{(\ell+1)}],$

or by simply always setting it to

$\displaystyle\hat{V}_{0}(t,x)=A_{\ell}\phi(x)\quad\textnormal{for}\quad
t\in[t^{(\ell)},t^{(\ell+1)}].$

Now, if $\hat{V}^{\textnormal{old}}_{0}$ is the approximation from the
previous policy iteration step, one could compute the empirical approximation
to the $\|.\|_{0}$-norm

$\displaystyle\|\hat{V}_{0}-\hat{V}^{\textnormal{old}}_{0}\|_{0,L,M}^{2}$
$\displaystyle=\dfrac{1}{L-1}\sum_{\ell=0}^{L-1}\|\hat{V}_{0}(t^{(\ell)},\cdot)-\hat{V}^{\textnormal{old}}_{0}(t^{(\ell)},\cdot)\|^{2}_{L^{2}(\Omega,\rho,M)}$
$\displaystyle=\dfrac{1}{(L-1)M}\sum_{\ell=0}^{L-1}\sum_{k=1}^{M}|\hat{V}_{0}(t^{(\ell)},x_{k})-\hat{V}^{\textnormal{old}}_{0}(t^{(\ell)},x_{k})|^{2}$

and stop the iteration once this norm difference becomes smaller than the
threshold $\delta$. However, since we are first and foremost interested in
obtaining a nearly optimal control $\alpha$, we instead add the change in the
controls $\alpha\propto\nabla_{x}\hat{V}_{0}$ and
$\alpha^{\textnormal{old}}\propto\nabla_{x}\hat{V}^{\textnormal{old}}_{0}$ and
stop the iteration once

$\displaystyle\dfrac{1}{(L-1)M}\sum_{\ell=0}^{L-1}\sum_{k=1}^{M}|\hat{V}_{0}(t^{(\ell)},x_{k})-\hat{V}^{\textnormal{old}}_{0}(t^{(\ell)},x_{k})|^{2}+\|\alpha(t^{(\ell)},x_{k})-\alpha^{\textnormal{old}}(t^{(\ell)},x_{k})\|^{2}<\delta.$
(26)

The reason for this is that the $L^{2}$-norm is agnostic to errors in the
gradients, which may arise due to overfitting. By requiring (26), we demand
that not only $\hat{V}_{0}$ but also the relevant part of the gradient
$\nabla_{x}\hat{V}_{0}$ converges. In that sense, the left-hand side of (26)
can be seen as an empirical approximation of an $H^{1}$-norm of
$\hat{V}_{0}-\hat{V}_{0}^{\textnormal{old}}$, where the norms for the
gradients are now weighted by $R$ and $g$ to represent only the gradient parts
relevant for the control.

##  7 Numerical tests

This chapter is concerned with numerical experiments that illustrate the
performance of the proposed DLR approximation222All computations are carried
out on an Intel Xeon Gold 6154 CPU 3.00GHz, openSUSE Leap 15.2 distribution..
We consider a problem of the form

$\displaystyle\dot{x}=Ax+\mathrm{nl}(x)+gu,$

where $x\in\mathbb{R}^{d}$, $A\in\mathbb{R}^{d\times d}$,
$g\in\mathbb{R}^{d}$, $u$ is scalar and $\mathrm{nl}$ is a smooth nonlinear
function with $\mathrm{nl}(0)=0$. In particular, the optimal control problem
is derived from a modified one dimensional heat equation

$\displaystyle\dfrac{\partial}{\partial t}{x}(s,t)$
$\displaystyle=\sigma\dfrac{\partial^{2}}{\partial
s^{2}}x(s,t)+x(s,t)^{3}+g(s)u(t),\quad$ $\displaystyle\text{ for
}(s,t)\in[-1,1]\times(0,T),$ $\displaystyle x(s,0)$
$\displaystyle=\tilde{x}_{0}(s),\quad$ $\displaystyle\text{ for }s\in[-1,1],$
$\displaystyle\dfrac{\partial}{\partial s}x(-1,t)$
$\displaystyle=\dfrac{\partial}{\partial s}x(1,t)=0\quad$ $\displaystyle\text{
for }t\in(0,T),$

with unstable reaction term $x(s,t)^{3}$, diffusion $\sigma>0$, scalar control
$u$ and initial state $\tilde{x}_{0}$. Note that due to the instability
introduced by the reaction term, this problem is generally more difficult to
control than most other canonically treated examples like viscous Burgers’
type equations, Allen-Kahn or degenerate Zeldovich equations [KK18, OSS21]
since the quadratic regulator usually provides a strong and mostly stable
controller for these types of problems. This however is not the case for the
nonlinear reaction problem defined above. We hence omit the mentioned
alternative examples and just note that our method can be applied with them as
well, although the difference to the linear quadratic regulator would turn out
to be small.

Our goal is to find a control $u$ such that the quadratic cost functional

$\displaystyle\tilde{J}(0,x_{0},u)=\int_{0}^{T}\|x(\cdot,t)\|_{L^{2}(\Omega)}^{2}+\gamma
u(t)^{2}~{}\mathrm{d}t+c_{T}\|x(\cdot,T)\|_{L^{2}(\Omega)},$

is minimal with $\gamma,c_{T}>0$. A semi-discretisation of the PDE with finite
differences at $d$ equidistant points $-1=s_{1}<\ldots<s_{d}=1$ leads to a an
ODE of the form

$\displaystyle\dot{x}$ $\displaystyle=Ax+x^{3}+gu,$ (27) $\displaystyle x(0)$
$\displaystyle=x_{0},$ (28)

with $x_{0}=(\tilde{x}_{0}(s_{1}),\ldots,\tilde{x}_{0}(s_{n}))^{\intercal}$,
$x(t)\in\mathbb{R}^{d}$,
$g=(g(s_{1}),\ldots,g(s_{d}))^{\intercal}\in\mathbb{R}^{d}$ and
$A\in\mathbb{R}^{d\times d}$ is given by

$\displaystyle A=\dfrac{\sigma}{h^{2}}\begin{pmatrix}-2&2&&\\\ 1&-2&1&\\\
&\ddots&\ddots&\ddots&\\\ &&1&\ddots&1\\\ &&&2&-2\end{pmatrix},\quad
h=s_{1}-s_{0}=\dfrac{2}{d-1}.$

The $x$-dependent term in the cost functional can be approximated using a
simple quadrature rule with nodes $s_{1},\ldots,s_{d}$ (here, we use the
rectangle rule with an additional node at the last grid point $s_{d}$). This
yields the new cost functional

$\displaystyle J(0,x_{0},u)=\int_{0}^{T}x(t)^{\intercal}Qx(t)+\gamma
u(t)^{2}~{}\mathrm{d}t+c_{T}x(T)^{\intercal}Qx(T),$ (29)

where

$\displaystyle Q=h\begin{pmatrix}1&&\\\ &\ddots&\\\ &&1\end{pmatrix},$

and $x(t)$ is understood to be the solution of $\dot{x}=Ax+x^{3}+gu$ with
starting value $x_{0}$. The control problem is now to find a control $u$ for
the nonlinear system (27) such that (29) is minimal for every starting value
$x_{0}$.

To specify the control problem, we choose the parameters $\sigma=1$,
$\gamma=0.1$, $c_{T}=1$ and $g=\chi_{[-0.4,0.4]}$ and discretise with $n=12$
equidistant grid points. The time horizon is $T=0.3$ and the time step size
$\tau=t_{i+1}-t_{i}$ is $0.001$, which is used for both the macro-intervals as
well as the micro-intervals of the policy iteration (see Section 6). The same
step size is also used to discretise the integral in (29) when computing the
costs. As a threshold for the policy iteration, $\delta=10^{-6}$ is set. We
choose $\Omega=(-2,2)^{d}$ and let $\rho$ be the uniform distribution on
$\Omega$. For the TT approximations, we use the first $n$
$H^{2}_{\textnormal{mix}}(\Omega)$-orthonormal polynomials
$\phi_{1},\ldots,\phi_{n}$ as basis functions (up to degree $n-1$). Here,
$H^{2}_{\textnormal{mix}}(\Omega)$ denotes the tensorised space
$\bigotimes_{\mu=1}^{d}H^{2}((-2,2))$, $H^{2}$ is the Sobolev space of twice
weakly differentiable functions. We set $n=9$, yielding a maximal polynomial
degree of $8$ in the basis. The rank of the TT manifold is chosen to be

$\displaystyle\mathbf{r}=(3,5,5,5,5,5,5,5,5,5,3).$

Note that by this the dimension of the approximation space is reduced from
$n^{d}=9^{12}>282$ trillion to a manageable number of degrees of freedom $\leq
nd\max_{\mu=1,\ldots,d}(r_{\mu})^{2}=2700$. The number of sample points used
to approximate the $L^{2}$-norm in (24) is chosen as

$\displaystyle M=6\cdot nd\max_{\mu=1,\ldots,d}(r_{\mu})^{2}=6\cdot 8\cdot
12\cdot 5^{2}=16200,$

which is a generous upper bound for the number of degrees of freedom of the
fit.

As a benchmark for assessing the performance of our method, we use the TT-
based approach from [OSS21] with the same hyper-parameters. To make this
precise, instead of solving (24) by means of our dynamical low-rank scheme,
$\hat{V}_{i}$ is approximated in each policy iteration step by sampling the
trajectories $x_{k}(t)$, $t\in[t_{i},t_{i+1}]$ of all sample points. With
this, the integrals

$\displaystyle\hat{V}_{i}(t_{i},x_{k})=\int_{t_{i}}^{t_{i+1}}\ell(t,x_{k}(t),\alpha(t,x_{k}(t))\mathrm{d}t+\hat{V}_{i+1}(t_{{i+1}},x_{k}(t_{i+1}))$

are evaluated subsequently. An approximation of $V_{i}(t_{i},\cdot)$ is then
obtained via a nonlinear fit of a rank-$\mathbf{r}$ TT to the resulting data-
target pairs $(x_{k},y_{k}=\hat{V}_{i}(t_{i},x_{k}))_{k=1}^{M}$, which is
performed by the ALS. Note that the authors in [OSS21] suggest replacing the
upper integral bound $t_{i+1}$ with $t_{i+l}$, $l>1$, where the trajectory on
$[t_{i+1},t_{i+l}]$ is controlled by the already computed (nearly optimal)
controls from previous steps, to remove the error associated with
$\hat{V}_{i+1}$ from the computation of $\hat{V}_{i}$. Since this greatly
increases the computational complexity, we stick with the above mentioned
“one-step scheme” and refer to this benchmark method as the Bellman method,
since it explicitly utilises Bellman’s principle in the form of the terminal
cost $\hat{V}_{i+1}$. Our method, utilising Dynamical Low Rank Approximation,
will be called the DLRA method. Even for the DLRA method we have found it
beneficial for stable convergence to compute some $\hat{V}_{i}$ with the
Bellman method before starting the dynamical low rank solver. In this example
the first $10$ of the $300$ approximations are computed in this way.

######  Remark 1.

In both the nonlinear fit required for the Bellman method and the linear fit
of our DLRA method, we add a regularisation term
$\delta\|\hat{V}_{i}(t_{i},\cdot)\|^{2}_{H^{2}_{\textnormal{mix}}(\Omega)}$ to
the minimisation functional. Due to the multilinear structure of the TT and
our choice of the basis functions as
$H^{2}_{\textnormal{mix}}(\Omega)$-orthonormal, this leads to local
minimisation problems of the form

$\displaystyle\min_{\mathbf{c}}\|M\mathbf{c}-\mathbf{y}\|_{2}^{2}+\delta\|\mathbf{c}\|^{2}_{F}$

in ALS (compare to [OSS21a]). Here, $\mathbf{c}\in\mathbb{R}^{r_{\mu-1}\times
n_{\mu}\times r_{\mu}}$ denotes the core that is currently optimised and
$\|.\|_{F}$ denotes the Frobenius norm. In both methods, we use
$\delta=10^{-10}$ but in ALS we successively lower $\delta$ via

$\displaystyle\delta\rightarrow\max(0.9,\|M\mathbf{c}-\mathbf{y}\|_{2}^{2}/\|\mathbf{y}\|^{2}_{2})\cdot\delta$

after every sweep. This is a purely heuristical rule to make sure the
regularisation is relaxed once the attractor of the global minimum is found.

######  Remark 2.

For the DLRA method we add an additional regularisation term

$\displaystyle\delta_{0}|B\phi(0)|^{2}$

to the minimisation in (24) since we know that the right-hand side satisfies
$F(t,A(t)\phi(0))=0$. Note that this can be realised by simply adding the
point $x_{M+1}=0$ to the set of samples $\\{x_{k}\\}_{k=1}^{M}$. Since this is
a hard constraint on the true solution, we set $\delta_{0}=10^{10}$.

As a second, classical benchmark, we consider the linear quadratic regulator
(LQR), resulting from linearising the problem around $x=0$. Since this
controller does not see the unstable reaction term, we expect poor performance
compared to both the Bellman and the DLRA method.

To compare the practical performance of the methods, two different sets of
initial conditions $\tilde{x}_{0}$ are generated. For the first set, we sample
a polynomial degree between $2$ and $20$ and then again randomly sample the
coefficients of a univariate polynomial of that degree. Denoting this
polynomial $p$, we then set $\tilde{x}_{0}(s)=(s-1)^{2}(s+1)^{2}p(s)$ to make
sure $\tilde{x}_{0}$ satisfies the Neumann boundary conditions. Finally, in
order to have interesting trajectories (27) for which the $x^{3}$-term
requires strong control beyond LQR, we normalise such that
$\max_{s\in[-1,1]}|\tilde{x}_{0}(s)|=1.9$. The second set of initial
conditions is generated by simply setting $\tilde{x}_{0}(s)\equiv c$ for
constants $c\in[1,2)$.

Figures 1 and 3 show the control values $u(t)$ along one trajectory of each
type of initial conditions. Figures 2 and 4 depict the mean costs over 500
randomly sampled initial conditions in each of the two cases, where we have
omitted those initial conditions for which the open-loop solver used to
compute the optimal control did not converge. Examining the graphs, we note
that the Bellman method and the DLRA method achieve similar, almost optimal
performance over the chosen test sets. Interestingly, the DLRA method actually
slightly outperforms the full Bellman method and is often closer to the
optimal control trajectories, which for instance can be seen in Figure 3. We
attribute this to the generalisation error of the Bellman method: even if the
value function approximation should be more accurate – due to a projection
directly onto the manifold – the associated optimisation is nonlinear and may
get stuck in local optima. In the DLRA method, we avoid this problem by coping
only with linear minimisation problems.

Figure 1: Control values $u(t)$ of the different controllers along the
trajectory of a fixed polynomial initial condition $x_{0}$. Figure 2: Sample-
mean costs
$\overline{\mathcal{J}_{u}(0,\cdot)}\approx\dfrac{1}{N}\sum_{k}^{N}J_{u}(0,x^{(k)}_{0})$
of polynomial initial conditions $x_{0}^{(k)}$ with the different controllers.
Figure 3: Control values $u(t)$ of the different controllers along the
trajectory of a fixed uniform initial condition
$x_{0}=c\cdot(1,\ldots,1)^{\intercal}$. In this example $c=1.28$ is used.
Figure 4: Sample-mean costs
$\overline{\mathcal{J}_{u}(0,\cdot)}\approx\dfrac{1}{N}\sum_{k}^{N}J_{u}(0,x^{(k)}_{0})$
of uniform initial conditions $x_{0}^{(k)}$ with the different controllers.

###  7.1 Computational cost and a hybrid approach

The distinct advantage of the DLRA method is its greatly reduced computational
cost. Table 1 contains the computation times for the two methods (Bellman and
DLRA), as well as their mean costs on the set of polynomial initial
conditions, with the same hyper-parameters and maximal polynomial degrees of
$4$, $6$ and $8$, respectively. We observe that the two methods achieve
comparable performance for degrees $6$ and $8$. However, the DLRA method
achieves this performance in roughly one tenth of the time that the Bellman
method requires. We stress again that the version we used is the fastest
version of the Bellman method available, since we employ the one-step scheme.
As discussed in [OSS21a, OSS21], this method also suffers from error
propagation due to a large number of time steps. Moreover our proposed method
projects onto the tangent space, whereas Bellman always tries to project onto
the tensor manifold.

The DLRA method performs significantly worse for a lower polynomial degree of
$4$. We attribute this to an effect that can be seen already for degree $8$ in
Figure 3. The DLRA controller drifts away from the true optimal control the
further it moves away from the terminal time $t=T$. This error seems to
originate from two main factors: for one, the true value function
$V^{*}(t,\cdot)$ successively moves further away from the manifold
$\mathcal{M}_{\mathbf{r}}$ even if the terminal condition satisfies
$c_{T}\in\mathcal{M}_{\mathbf{r}}$. To visualise that the true solution does
not stay on the manifold, the relative norm error of the last tangent fit in
each policy iteration is plotted over time in Figure 5. Note that these errors
should be close to $0$ if the solution to the GHJB equation is an element of
the manifold. Instead, the errors increase monotonically over time. The second
major source of error is the retraction after every Euler step. Both sources
of errors get worse for lower degrees because of the restricted manifold.
Hence, a degree of $4$, which is perfectly feasible for the Bellman method,
produces bad results with the DLRA method. Note that the observed behaviour is
expected.

| Bellman | DLRA | Hybrid
---|---|---|---
| comp. time | mean cost | comp. time | mean cost | comp. time | mean cost
pol. deg. $4$ | 3078.44 | 1.8822 | 333.29 | 2.6147 | 909.65 | 1.8804
pol. deg. $6$ | 4270.33 | 1.8801 | 421.52 | 1.8802 | 1851.93 | 1.8798
pol. deg. $8$ | 5967.91 | 1.8800 | 499.96 | 1.8799 | – | –
Table 1: Computation time of the methods in seconds as well as mean costs of
polynomial initial conditions for different maximal polynomial degrees of the
basis functions. The mean optimal cost is 1.8793.

This observation leads to a natural formulation of a hybrid method, possibly
alleviating the main weaknesses of both methods. These are the high
computational complexity for the Bellman method and error accumulation for the
DLRA method. The hybrid method uses DLRA updates but after each $m$ steps,
instead of computing $\hat{V}_{i}$ with the regular DLRA update, it performs a
full Bellman update [OSS21] with an $m$-step scheme

$\displaystyle\hat{V}_{i}(t_{i},x_{k})=\int_{t_{i}}^{t_{i+m}}\ell(t,x_{k}(t),\alpha(t,x_{k}(t))dt+\hat{V}_{i+m}(t_{{i+m}},x_{k}(t_{i+m})).$
(30)

For $m=1$ this method is equivalent to the Bellman method, for $m$ greater
than the number of total time steps it is equivalent to the DLRA method. For
any intermediate $m$ it periodically performs one costly but accurate Bellman
update in between fast DLRA updates. Since the maximal number of consecutive
DLRA steps is now $m$, the DLRA solver is prevented from drifting too far away
from the real solution, before being corrected again by the Bellman update,
yielding a new (more accurate) initial condition. Note in particular that the
evaluation of (30) does not include any $\hat{V}_{j}$ computed with the DLRA
method. Hence, after every $m$ steps, the accumulated error of the DLRA steps
is reset to $0$. Globally, only the error of the $m$-step Bellman updates (30)
accumulates.

The results for the hybrid method with $m=10$ are depicted in Table 1 for
degrees 4 and 6. We remark that for polynomials of degree 4, the hybrid scheme
provides an essential improvement with respect to accuracy when compared to
both Bellman and DLRA. There is an improvement for degree 6 but compared to
DLRA this effect is not pronounced. Surprisingly, for a sufficiently accurate
model, DLRA alone was sufficiently accurate. The case of degree 8 is omitted
since the DLRA controller is already nearly optimal in that case. The periodic
$10$-step Bellman updates with intermediate DLRA steps are sufficient to
outperform the full 1-step Bellman method, but at much lower computational
costs. From the perspective of the DLRA method, the periodic Bellman updates
enable the use of more restricted manifolds.

Figure 5: Relative norm error of the last tangent fit (see Appendix A) in each
policy iteration step of the DLRA solver over time and for different
polynomial degrees.

##  8 Concluding remarks

In this paper we present a novel method to approximate optimal feedback laws
for optimal control problems. The proposed method utilizes a tensor train
compression to break the curse of dimensionality of a multivariate polynomial
ansatz space. Moreover, it employs an empirical version of the Dirac-Frenkel
variational principle to solve the HJB equation. The method was tested
numerically on a canonical benchmark example which is difficult to control
with standard methods, and demonstrated to achieve near optimal performance
with greatly reduced computation time compared to state-of-the-art methods.

In the experiments it comes as no surprise that the proposed method works
quite well for short time intervals. However, it is striking that we can also
observe that with a sufficiently good model – meaning an adequate polynomial
degree in our case – the method even performs well on a large time horizon.
Although we have not considered infinite horizon problems yet, as long as we
know stabilizing controls, the present approach probably is applicable as
well. Moreover, for large time horizons we have presented a robust hybrid
method.

We would like to point out that the present successful approach strongly
exploits the explicit knowledge about the geometry of the considered model
class, i.e. (multi-)polynomial tensor trains in our setting. This advantage is
something which cannot be easily transferred to a neural network setting.

We expect the method to also perform favourably with higher dimensional
problems, which might be a future research topic. We predict that this will
require some form of rank adaptivity to retain the computational advantage
over state-of-the-art methods while achieving similar levels of accuracy. Rank
adaptivity can be incorporated very naturally in the proposed DLRA method:
instead of the full retraction onto the manifold $\mathcal{M}_{\mathbf{r}}$
after every step of the solver, one could round the TT based on an adaptive
threshold. Analysing the effect of a changing manifold on the Dirac-Frenkel
variational principle might be an interesting topic for future work.

As a second direction, the method could be applied to stochastic optimal
control problems. There, the GHJB equation (17) gets an additional Laplacian
term $\propto\Delta_{x}V(t,x)$, turning it into a Kolmogorov-Backward type
equation. Equations of this type for instance govern the time development of
observables of Itô diffusion processes. The application of our method to such
problems is currently being investigated.

## Acknowledgements

Martin Eigel acknowledges the partial support of the DFG SPP 1886 “Polymorphic
Uncertainty Modelling for the Numerical Design of Structures”. David Sommer
acknowledges support by the ProFIT project “ReLkat – Reinforcement Learning
for complex automation engineering”.

##  Appendix A Details of the empirical risk minimisation

We detail how to reduce the minimisation in (24) to a standard system of
linear equations. To achieve this, we use the characterization of the tangent
space given by Theorem 4 and represent an element $\delta U\in\tau(X)$ of the
tangent space as a vector $\mathbf{x}\in\mathbb{R}^{n_{X}}$. The first step
towards this representation is the parametrisation of the spaces
$U_{\mu}^{\ell}$.

###  A.1 A parametrisation of the tangent space

By Theorem 4, $U_{\mu}^{\ell}$ is precisely the set of all
$r_{\mu-1}n_{\mu}\times r_{\mu}$-matrices whose columns are orthogonal to the
columns of $L(U_{\mu})$. Let $QR=L(U_{\mu})$ be the QR decomposition and
denote the orthonormal columns of $Q$ by $q_{1},\ldots,q_{r_{\mu}}$. By the
Gram-Schmidt procedure, we can expand the columns to an orthonormal basis of
$\mathbb{R}^{r_{\mu-1}n_{\mu}}$ and denote the additional vectors by
$\hat{q}_{1}=q_{r_{\mu}+1},\ldots,\hat{q}_{r_{\mu-1}n_{\mu}-r_{\mu}}=q_{r_{\mu-1}n_{\mu}}$.
Now, let $W_{\mu}\in U_{\mu}^{\ell}$ and denote its $j$-th column by $w_{j}$.
Then there are coefficients $c_{j,k}$ such that

$\displaystyle w_{j}=\sum_{k=1}^{r_{\mu-1}n_{\mu}-r_{\mu}}c_{j,k}\hat{q}_{k}.$

In total we get $(r_{\mu-1}n_{\mu}-r_{\mu})r_{\mu}$ coefficients $c_{j,k}$,
which are stored in a vector

$\displaystyle\mathbf{x}_{\mu}$
$\displaystyle=(c_{1,1},\ldots,c_{1,r_{\mu-1}n_{\mu}-r_{\mu}},c_{2,1},\ldots,c_{2,r_{\mu-1}n_{\mu}-r_{\mu}},\ldots,\ldots,c_{r_{\mu},1},\ldots,c_{r_{\mu},r_{\mu-1}n_{\mu}-r_{\mu}})^{\intercal}$
$\displaystyle\in\mathbb{R}^{(r_{\mu-1}n_{\mu}-r_{\mu})r_{\mu}}=\mathbb{R}^{r_{\mu-1}n_{\mu}r_{\mu}-r^{2}_{\mu}}.$

From now on we always identify an element of $U_{\mu}^{\ell}$ with its
coefficient vector $\mathbf{x}_{\mu}$. Elements of
$C_{d}=\mathbb{R}^{r_{d-1}\times n_{d}\times r_{d}}$ are represented in the
same manner with the only difference that the sum in each column
representation goes from $k=1$ to $r_{d-1}n_{d}$ and the $\hat{q}_{i}$ can be
chosen as the canonical basis in $\mathbb{R}^{r_{d-1}n_{d}}$.

We eventually can represent an element of the tangent space $\delta
U\in\tau(X)$ by the concatenation of its coefficients vectors,

$\displaystyle\delta
U\cong\mathbf{x}=(\mathbf{x}^{\intercal}_{1},\ldots,\mathbf{x}^{\intercal}_{d})^{\intercal}\in\mathbb{R}^{n_{X}}.$

Since this becomes important when solving the regression problem (24) on the
tangent space later on, we define a “lift”

$\displaystyle\mathcal{L}_{\mu}:\mathbb{R}^{r_{\mu-1}n_{\mu}r_{\mu}-r_{\mu}^{2}}\longrightarrow\mathbb{R}^{r_{\mu-1}n_{\mu}r_{\mu}},$

which maps the coefficient vector of the gauged representation to the
vectorised entries of the corresponding tensor in $U_{\mu}^{\ell}$. This is
achieved by means of a
$r_{\mu-1}n_{\mu}r_{\mu}\times(r_{\mu-1}n_{\mu}r_{\mu}-r_{\mu}^{2})$-lifting
matrix

$\displaystyle
Z_{\mu}=\begin{pmatrix}Q_{\mu}^{\perp}&\mathbf{0}&\ldots&\ldots&\mathbf{0}\\\
\mathbf{0}&Q_{\mu}^{\perp}&\mathbf{0}&&\\\ &&\ddots&&\\\ &&&&\mathbf{0}\\\
\mathbf{0}&&&\mathbf{0}&Q_{\mu}^{\perp}\end{pmatrix},\qquad
Q_{\mu}^{\perp}=[\hat{q}_{1},\ldots,\hat{q}_{r_{\mu-1}n_{\mu}-r_{\mu}}]\in\mathbb{R}^{r_{\mu-1}n_{\mu}\times(r_{\mu-1}n_{\mu}-r_{\mu})}.$

By construction, $Z_{\mu}\mathbf{x}_{\mu}$ is the concatenation of the columns
of $L(W_{\mu})$. Hence, $\mathcal{L}_{\mu}(\mathbf{x}_{\mu})$ can be obtained
by $\mathcal{L}_{\mu}(\mathbf{x}_{\mu})=Z_{\mu}\mathbf{x}_{\mu}$.

###  A.2 Solving the system of linear equations

We examine problem (24) in a more general setting. Let
$U\in\mathcal{M}_{\mathbf{r}}$ be $d$-orthogonal and consider the minimisation
problem

$\displaystyle\min_{T\in\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})}\sum_{k=1}^{M}\left|T[x^{(k)}]-y^{(k)}\right|^{2},$
(31)

where $(x^{(k)},y^{(k)})_{k=1}^{M}\subset\mathbb{R}^{d}\times\mathbb{R}$ is a
set of data-target pairs. This leads to

$\displaystyle~{}\quad\min_{T\in\mathcal{T}_{U}(\mathcal{M}_{\mathbf{r}})}\sum_{k=1}^{M}|T[x^{(k)}]-y^{(k)}|^{2}$
$\displaystyle=\min_{W\in X}\sum_{k=1}^{M}|\tau(W)[x^{(k)}]-y^{(k)}|^{2}$
$\displaystyle=\min_{W\in
X}\sum_{k=1}^{M}\left|\sum_{\mu=1}^{d}\left[\sum_{j_{1},\ldots,j_{d}}^{n_{1},\ldots,n_{d}}U_{1}(j_{1})\cdot\ldots\cdot
W_{\mu}(j_{\mu})\cdot\ldots\cdot
U_{d}(j_{d})\phi_{j_{1}}(x_{1}^{(k)})\ldots\phi_{j_{d}}(x_{d}^{(k)})\right]-y^{(k)}\right|^{2}$
$\displaystyle=\min_{W\in
X}\left\|\mathfrak{O}(W)-\textbf{y}\right\|_{2}^{2},$

where $\mathbf{y}=(y^{(1)},\ldots,y^{(M)})^{\intercal}$ and the operator
$\mathfrak{O}:X\longrightarrow\mathbb{R}^{M}$ is defined by

$\displaystyle\mathfrak{O}(W)$
$\displaystyle=\sum_{\mu=1}^{d}\mathfrak{C}_{\mu}(W_{\mu}),\quad\textnormal{for}\quad
W=(W_{1},\ldots,W_{d}),$ (32)

with

$\displaystyle\mathfrak{C}_{\mu}:\mathbb{R}^{r_{\mu-1}\times n_{\mu}\times
r_{\mu}}\longrightarrow\mathbb{R}^{M},$
$\displaystyle(\mathfrak{C}_{\mu}(W_{\mu}))_{k}=\left[\sum_{j_{1},\ldots,j_{d}}^{n_{1},\ldots,n_{d}}U_{1}(j_{1})\cdot\ldots\cdot
W_{\mu}(j_{\mu})\cdot\ldots\cdot
U_{d}(j_{d})\phi_{j_{1}}(x_{1}^{(k)})\ldots\phi_{j_{d}}(x_{d}^{(k)})\right].$

Note that $\mathfrak{C}_{\mu}$ is a linear tensor operator in
$\mathbb{R}^{M\times r_{\mu-1}\times n_{\mu}\times r_{\mu}}$, which we can
transfer into a matrix $C_{\mu}\in\mathbb{R}^{M\times
r_{\mu-1}n_{\mu}r_{\mu}}$ by successive unfolding

$\displaystyle\mathbb{R}^{M\times r_{\mu-1}\times n_{\mu}\times
r_{\mu}}\longrightarrow\mathbb{R}^{M\times r_{\mu-1}n_{\mu}\times
r_{\mu}}\longrightarrow\mathbb{R}^{M\times r_{\mu-1}n_{\mu}r_{\mu}}.$

At the first stage, the operator $\mathfrak{C}_{\mu}$ acts on a tensor
$W_{\mu}\in U_{\mu}^{\ell}$. At the second stage, it acts on the left
unfolding $L(W_{\mu})$. And at the third stage, the matrix $C_{\mu}$ acts on
the concatenation of the columns of $L(W_{\mu})$. By the previous section,
this concatenation is given by $Z_{\mu}\mathbf{x}_{\mu}$. We hence have
$\mathfrak{C}_{\mu}(W_{\mu})=C_{\mu}Z_{\mu}\mathbf{x}_{\mu}$, leading to

$\displaystyle\mathfrak{O}(W)=\sum_{\mu=1}^{d}C_{\mu}Z_{\mu}x_{\mu}=A\mathbf{x},$

where $A=\left[C_{1}Z_{1},\ldots,C_{d-1}Z_{d-1},C_{d}\right]$ (note that
$Z_{d}\equiv I_{d}$). We have thus transformed (31) to a standard system of
linear equations

$\displaystyle\hat{\mathbf{x}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathbb{R}^{n_{X}}}\|A\mathbf{x}-\mathbf{y}\|^{2}\Longleftrightarrow
A^{\intercal}A\hat{\mathbf{x}}=A^{\intercal}\mathbf{y}.$

Once a solution is found by standard methods, we recover $W$ from
$\hat{\mathbf{x}}$ by reshaping of the component vectors $\mathbf{x}_{\mu}$.

######  Remark 3.

We would like to make two remarks about the implementation. First, note that
the matrices $Z_{\mu}$ do not have to be stored in order to compute the
product $C_{\mu}Z_{\mu}$ since we can compute

$\displaystyle C_{\mu}Z_{\mu}$
$\displaystyle=[C_{\mu}[0:r_{\mu-1}n_{\mu}]|\ldots|C_{\mu}[(r_{\mu}-1)(r_{\mu-1}n_{\mu}):r_{\mu}r_{\mu-1}n_{\mu}]]\cdot\textnormal{diag}(Q_{\mu}^{\perp},\ldots,Q_{\mu}^{\perp})$
$\displaystyle=[C_{\mu}[0:r_{\mu-1}n_{\mu}]Q_{\mu}^{\perp}|\ldots|C_{\mu}[(r_{\mu}-1)(r_{\mu-1}n_{\mu}):r_{\mu}r_{\mu-1}n_{\mu}]Q_{\mu}^{\perp}].$

Second, note that $U_{\mu}^{\ell}$ is $0$-dimensional if
$r_{\mu-1}n_{\mu}-r_{\mu}=0$. In this case, the space consists only of the
tensor of constant zeros $\mathbf{0}\in\mathbb{R}^{r_{\mu-1}\times
n_{\mu}\times r_{\mu}}$ and hence $W_{\mu}=\mathbf{0}$. No basis coefficients
$\mathbf{x}_{\mu}$ need to be computed. Consequently, the index $\mu$ can be
skipped entirely during optimisation. By this, $A$ and $\mathbf{x}$ become

$\displaystyle A$
$\displaystyle=[C_{1}Z_{1},\ldots,C_{\mu-1}Z_{\mu-1},C_{\mu+1}Z_{\mu+1},\ldots,C_{d-1}Z_{d-1},C_{d}],$
$\displaystyle\mathbf{x}$
$\displaystyle=(\mathbf{x}^{\intercal}_{1},\ldots,\mathbf{x}^{\intercal}_{\mu-1},\mathbf{x}^{\intercal}_{\mu+1},\ldots,\mathbf{x}^{\intercal}_{d})^{\intercal}.$

## References

  * [AF18] Marianne Akian and Eric Fodjo  “Probabilistic Max-Plus Schemes for Solving Hamilton-Jacobi-Bellman Equations”  In _Springer INdAM Series_ , 2018, pp. 183–209  DOI: 10.1007/978-3-030-01959-4_9
  * [AGL08] Marianne Akian, Stéphane Gaubert and Asma Lakhoua  “The max-plus finite element method for solving deterministic optimal control problems: basic properties and convergence analysis”  In _SIAM Journal on Control and Optimization_ 47.2 SIAM, 2008, pp. 817–848
  * [AS19] Alessandro Alla and Luca Saluzzi  “A HJB-POD approach for the control of nonlinear PDEs on a tree structure”, 2019
  * [AKK21] Behzad Azmi, Karl Kunisch and Dante Kalise  “Optimal Feedback Law Recovery by Gradient-Augmented Sparse Polynomial Regression”, 2021
  * [Bac+21] Markus Bachmayr, Henrik Eisenmann, Emil Kieri and André Uschmajew  “Existence of dynamical low-rank approximations to parabolic problems”  In _Mathematics of Computation_ , 2021
  * [BSU16] Markus Bachmayr, Reinhold Schneider and André Uschmajew  “Tensor Networks and Hierarchical Tensors for the Solution of High-Dimensional Partial Differential Equations”  In _Found. Comput. Math._ 16.6 Berlin, Heidelberg: Springer-Verlag, 2016, pp. 1423–1472  DOI: 10.1007/s10208-016-9317-9
  * [BC97] Martino Bardi and Italo Capuzzo-Dolcetta  “Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations”, 1997
  * [BD+97] Martino Bardi and Italo Capuzzo Dolcetta  “Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations”  Springer, 1997
  * [Bay+21] Christian Bayer, Martin Eigel, Leon Sallandt and Philipp Trunschke  “Pricing high-dimensional Bermudan options with hierarchical tensor formats”, 2021
  * [Bel57] Richard Bellman  “Dynamic Programming”  Dover Publications, 1957
  * [Ber05] Dimitri P. Bertsekas  “Dynamic Programming and Optimal Control”  Belmont, MA, USA: Athena Scientific, 2005
  * [BGP61] VG Boltyanskiy, Revaz Valer’yanovich Gamkrelidze and Lev Semenovich Pontryagin  “Theory of optimal processes”, 1961
  * [CA13] Eduardo F Camacho and Carlos Bordons Alba  “Model predictive control”  Springer science & business media, 2013
  * [CKL21] Gianluca Ceruti, Jonas Kusch and Christian Lubich  “A rank-adaptive robust integrator for dynamical low-rank approximation”, 2021  arXiv:2104.05247 [math.NA]
  * [CL20] Gianluca Ceruti and Christian Lubich  “An unconventional robust integrator for dynamical low-rank approximation”, 2020  arXiv:2010.02022 [math.NA]
  * [Con20] Dajana Conte  “Dynamical low-rank approximation to the solution of parabolic differential equations”  In _Applied Numerical Mathematics_ 156 Elsevier, 2020, pp. 377–384
  * [DLM19] Jérôme Darbon, Gabriel Provencher Langlois and Tingwei Meng  “Overcoming the curse of dimensionality for some Hamilton–Jacobi partial differential equations via neural network architectures”  In _Research in the Mathematical Sciences_ 7, 2019, pp. 1–50
  * [DKK21] Sergey Dolgov, Dante Kalise and Karl Kunisch  “Tensor Decomposition Methods for High-dimensional Hamilton-Jacobi-Bellman Equations”, 2021  arXiv:1908.01533 [math.OC]
  * [EST20] Martin Eigel, Reinhold Schneider and Philipp Trunschke  “Convergence bounds for empirical nonlinear least-squares”  In _arXiv preprint arXiv:2001.00639_ , 2020
  * [Eig+19] Martin Eigel, Reinhold Schneider, Philipp Trunschke and Sebastian Wolf  “Variational Monte Carlo—bridging concepts of machine learning and high-dimensional partial differential equations”  In _Advances in Computational Mathematics_ 45.5 Springer, 2019, pp. 2503–2532
  * [Fac+20] Konstantin Fackeldey, Mathias Oster, Leon Sallandt and Reinhold Schneider  “Approximative Policy Iteration for Exit Time Feedback Control Problems driven by Stochastic Differential Equations using Tensor Train format”, 2020  arXiv:2010.04465 [math.OC]
  * [Fal87] Maurizio Falcone  “A numerical approach to the infinite horizon problem of deterministic control theory”  In _Applied Mathematics and Optimization_ 15, 1987, pp. 1–13  DOI: 10.1007/BF01442644
  * [FK14] Maurizio Falcone and Dante Kalise  “A high-order semi-Lagrangian/finite volume scheme for Hamilton-Jacobi-Bellman-Isaacs equations”  In _IFIP Advances in Information and Communication Technology_ 443, 2014  DOI: 10.1007/978-3-662-45504-3_10
  * [FLS94] Maurizio Falcone, Piero Lanucara and Alessandra Seghini  “A splitting algorithm for Hamilton-Jacobi-Bellman equations”  In _Applied Numerical Mathematics_ 15.2 Elsevier, 1994, pp. 207–218
  * [Hac12] Wolfgang Hackbusch  “Tensor Spaces and Numerical Tensor Calculus”, 2012  DOI: 10.1007/978-3-642-28027-6
  * [Hac14] Wolfgang Hackbusch  “Numerical tensor calculus”  In _Acta numerica_ 23, 2014, pp. 651–742  DOI: 10.1017/S0962492914000087
  * [HS14] Wolfgang Hackbusch and Reinhold Schneider  “Tensor spaces and hierarchical tensor representations”  In _Extraction of quantifiable information from complex systems_ Springer, 2014, pp. 237–261
  * [HRS12] Sebastian Holtz, Thorsten Rohwedder and Reinhold Schneider  “On Manifolds of Tensors of Fixed TT-Rank”  In _Numer. Math._ 120.4 Berlin, Heidelberg: Springer-Verlag, 2012, pp. 701–731  DOI: 10.1007/s00211-011-0419-7
  * [HRS12a] Sebastian Holtz, Thorsten Rohwedder and Reinhold Schneider  “On manifolds of tensors of fixed TT-rank”  In _Numerische Mathematik_ 120, 2012, pp. 701–731
  * [How60] R.. Howard  “Dynamic Programming and Markov Processes”  Cambridge, MA: MIT Press, 1960
  * [IRZ21] Kazufumi Ito, Christoph Reisinger and Yufei Zhang  “A neural network based policy iteration algorithm with global H2-superlinear convergence for stochastic games on domains”  In _Found. Comput. Math._ 21, 2021, pp. 331–374
  * [KDK13] B Kafash, A Delavarkhalafi and SM Karbassi  “Application of variational iteration method for Hamilton–Jacobi–Bellman equations”  In _Applied Mathematical Modelling_ 37.6 Elsevier, 2013, pp. 3917–3928
  * [KK18] Dante Kalise and Karl Kunisch  “Polynomial Approximation of High-Dimensional Hamilton–Jacobi–Bellman Equations and Applications to Feedback Control of Semilinear Parabolic PDEs”  In _SIAM Journal on Scientific Computing_ 40.2 Society for Industrial & Applied Mathematics (SIAM), 2018, pp. A629–A652  DOI: 10.1137/17m1116635
  * [KKD19] Dante Kalise, Karl Kunisch and Sergey Dolgov  “Tensor Decomposition Methods for High-dimensional Hamilton-Jacobi-Bellman Equations”, 2019
  * [KW17] Wei Kang and Lucas C. Wilcox  “Mitigating the curse of dimensionality: sparse grid characteristics method for optimal feedback control and HJB equations”  In _Computational Optimization and Applications_ 68, 2017, pp. 289–315
  * [KLW16] Emil Kieri, Christian Lubich and Hanna Walach  “Discretized Dynamical Low-Rank Approximation in the Presence of Small Singular Values”  In _SIAM J. Numer. Anal._ 54, 2016, pp. 1020–1038
  * [KL07] Othmar Koch and Christian Lubich  “Dynamical Low-Rank Approximation”  In _SIAM J. Matrix Anal. Appl._ 29, 2007, pp. 434–454
  * [KL10] Othmar Koch and Christian Lubich  “Dynamical tensor approximation”  In _SIAM Journal on Matrix Analysis and Applications_ 31.5 SIAM, 2010, pp. 2360–2375
  * [LO13] Christian Lubich and Ivan Oseledets  “A projector-splitting integrator for dynamical low-rank approximation”  In _BIT_ 54, 2013  DOI: 10.1007/s10543-013-0454-0
  * [LOV15] Christian Lubich, Ivan Oseledets and Bart Vandereycken  “Time Integration of Tensor Trains”  In _SIAM Journal on Numerical Analysis_ 53, 2015, pp. 917–941  DOI: 10.1137/140976546
  * [Lub+13] Christian Lubich, Thorsten Rohwedder, Reinhold Schneider and Bart Vandereycken  “Dynamical Approximation By Hierarchical Tucker And Tensor-Train Tensors”  In _SIAM Journal on Matrix Analysis and Applications_ 34, 2013, pp. 470–494  DOI: 10.1137/120885723
  * [Luo+14] Biao Luo, Huai-Ning Wu, Tingwen Huang and Derong Liu  “Data-based approximate policy iteration for affine nonlinear continuous-time optimal control design”  In _Automatica_ 50.12 Elsevier, 2014, pp. 3281–3290
  * [McL64] A.D. McLachlan  “A variational solution of the time-dependent Schrodinger equation”  In _Molecular Physics_ 8.1 Taylor & Francis, 1964, pp. 39–44  DOI: 10.1080/00268976400100041
  * [Mur35] FD Murnaghan  “J. frenkel, wave mechanics; advanced general theory”  In _Bulletin of the American Mathematical Society_ 41.11 American Mathematical Society, 1935, pp. 776–776
  * [NGK19] Tenavi Nakamura-Zimmerer, Qi Gong and Wei Kang  “Adaptive Deep Learning for High-Dimensional Hamilton-Jacobi-Bellman Equations”, 2019
  * [NR21] Nikolas Nüsken and Lorenz Richter  “Solving high-dimensional Hamilton–Jacobi–Bellman PDEs using neural networks: perspectives from the theory of controlled diffusions and measures on path space”  In _Partial Differential Equations and Applications_ 2, 2021  DOI: 10.1007/s42985-021-00102-x
  * [Ose11] Ivan Oseledets  “Tensor-Train Decomposition”  In _SIAM J. Scientific Computing_ 33, 2011, pp. 2295–2317  DOI: 10.1137/090752286
  * [OT09] Ivan Oseledets and E. Tyrtyshnikov  “Breaking the Curse of Dimensionality, Or How to Use SVD in Many Dimensions”  In _SIAM J. Sci. Comput._ 31, 2009, pp. 3744–3759  DOI: 10.1137/090748330
  * [OSS21] Mathias Oster, Leon Sallandt and Reinhold Schneider  “Approximating optimal feedback controllers of finite horizon control problems using hierarchical tensor formats”, 2021  arXiv:2104.06108 [math.OC]
  * [OSS21a] Mathias Oster, Leon Sallandt and Reinhold Schneider  “Approximating the Stationary Bellman Equation by Hierarchical Tensor Products”, 2021  arXiv:1911.00279 [math.OC]
  * [Pon87] Lev Semenovich Pontryagin  “Mathematical theory of optimal processes”  CRC press, 1987
  * [RSN21] Lorenz Richter, Leon Sallandt and Nikolas Nüsken  “Solving high-dimensional parabolic PDEs using the tensor train format”, 2021
  * [Sal21] Leon Jasper Sallandt  “Computing High-Dimensional Value Functions of Optimal Feedback Control Problems using the Tensor-Train Format”, 2021
  * [SL79] George N. Saridis and C.. Lee  “An Approximation Theory of Optimal Control for Trainable Manipulators”  In _IEEE Transactions on Systems, Man, and Cybernetics_ 9, 1979, pp. 152–159
  * [Ste16] Michael Maximilian Steinlechner  “Riemannian Optimization for Solving High-Dimensional Problems with Low-Rank Tensor Structure”  Lausanne: EPFL, 2016, pp. 165  DOI: 10.5075/epfl-thesis-6958
  * [TAK17] Daniela Tonon, Maria Aronna and Dante Kalise  “Optimal Control: Novel Directions and Applications”, 2017  DOI: 10.1007/978-3-319-60771-9
  * [TSG21] Philipp Trunschke, Reinhold Schneider and Michael Götte  “A block-sparse Tensor Train Format for sample-efficient high-dimensional Polynomial Regression”  In _Frontiers in Applied Mathematics and Statistics_ Frontiers, 2021, pp. 57
  * [ZH21] Mo Zhou and Jiequn Han  “Actor-Critic Method for High Dimensional Static Hamilton–Jacobi–Bellman Partial Differential Equations based on Neural Networks”, 2021