Source-linked AI summary
Dynamical low-rank approximations for the optimal control of matrix ODEs
Benjamin Carrel, Tommaso Vanzan
TL;DR
This paper addresses efficient open-loop optimal control for systems governed by matrix ODEs. It develops two dynamical low-rank approaches, derives gradient error estimates, and reports a lower-order memory requirement than a classical full-rank iteration.
Problem
The paper addresses efficient numerical solution of open-loop optimal control problems constrained by matrix ODEs.
Source: §1 · p. 1Method
OtP projects the state and adjoint equations separately, whereas PtO first approximates the state equation and derives optimality conditions with an adjoint in the state’s tangent space.
Source: §1.1 · p. 2Results
Gradient error estimates depend on initial-state approximation, normal vector-field components, dynamics stability, and approach-specific terminal initialization or curvature terms.
Source: §1.1 · p. 2Takeaways & Limitations
Each gradient-descent iteration requires O(Nt(m + n + r)r) memory, potentially substantially less than the O(Ntmn) required by a classical full-rank iteration.
Source: §5.3 · p. 19In OtP, the zero terminal adjoint is not rank-r, so initialization requires a rank-r perturbation or a full-rank warm start followed by truncation.
Source: §4.1 · p. 8
Abstract
from arXiv · showhide
We consider the open-loop optimal control of a dynamical system whose state is described by a time-dependent matrix $X(t)\in\mathbb{R}^{m\times n}$. Gradient-based optimization methods may incur substantial computational and memory costs, particularly for large-scale systems and fine temporal discretizations. To alleviate this computational burden, we develop a mathematical framework to use Dynamical Low Rank Approximations (DLRA) in optimal control problems. We study two distinct paradigms: the first one, called Optimize-then-Project (OtP), first derives the optimality system and then approximates both the state and adjoint equations using DLRA. The second one, called Project-then-Optimize (PtO), first replaces the state equation with its DLR approximation, and then derives the optimality conditions for the resulting optimization problem. The two approaches differ fundamentally concerning the role of the adjoint variable. In OtP, the adjoint variable evolves on the manifold of fixed rank matrices, whereas in PtO, the adjoint variable belongs to the tangent space to the fixed-rank manifold at the evolving state variable. This latter formulation leads to a covariant ODE describing the dynamics of the adjoint variable as the underlying tangent space evolves in time. For both approaches we establish rigorous error estimates for the approximated gradients, identifying how their accuracy depends on the low-rank approximations of initial and terminal conditions, the components of the state and adjoint vector fields normal to the relevant tangent spaces, the stability properties of the dynamics, and the manifold curvature. PtO further requires tailored time discretization strategies for the adjoint covariant ODE, which we show to be robust against small singular values of the state variable. To illustrate the framework, we apply these methods to the optimal control of a random parabolic PDE.
1. Introduction
The paper develops a mathematical framework for using dynamical low-rank approximations in open-loop matrix-ODE optimal control, where dense state and adjoint storage can be costly. It compares Optimize-then-Project and Project-then-Optimize, establishes gradient-error estimates, and discusses time integration and a random parabolic PDE application.
1. Introduction
- Dense storage of state and adjoint matrices over Nt time steps can cost O(mnNt), motivating low-rank representations for large-scale optimal control.The introduction cites matrix and tensor ODE examples including semidiscrete PDEs, quantum density matrices, kinetic equations, and random PDEs.Source: §1 · pp. 1–2
1.1. Contributions and outline.
- OtP projects the full state and adjoint optimality equations independently, whereas PtO first projects the state equation and yields an adjoint in the evolving state’s tangent space.The PtO adjoint therefore satisfies a covariant ODE, and the two approaches are not equivalent.Source: §1.1 · p. 2
- Gradient-error estimates for both approaches depend on low-rank initialization, normal vector-field components, and dynamical stability, with additional terminal-initialization error for OtP and curvature-related terms for PtO.Both estimates recover the full-rank gradient when r = min{m, n}.Source: §1.1 · p. 2
- For PtO’s backward covariant ODE, the authors propose a projected Euler method and prove robustness with respect to small singular values of the state.The factor equations are generally not stiff even as the smallest singular value tends to zero, but require the state and its time derivative at every instant.Source: §1.1 · pp. 2–3
- The framework is applied to optimal control of a random parabolic PDE with random coefficients and initial data, with experiments assessing error bounds, minimization effects, and computational speed-ups.Source: §1.1 · p. 3
1.2. Layout.
- The paper’s stated organization proceeds from the model problem and fixed-rank DLRA background to the two approaches, time integration, and numerical application.Source: §1.2 · p. 3
2. Problem setting
The problem setting is an open-loop control problem for a matrix ODE with a tracking objective and control regularization. The authors establish existence of a minimizer under their assumptions and motivate low-rank gradient methods by the cost of forward and backward dense-matrix solves.
2. Problem setting
- The objective combines tracking error against a desired trajectory with Tikhonov regularization weighted by ν > 0, which penalizes large control inputs.Source: §2 · p. 3
- Under the stated assumptions on the vector field, the problem has a minimizer, and the optimality conditions comprise forward state and backward adjoint equations.The vector field is written as f(t, X) + B(t, X)z, with f and B measurable in time and globally Lipschitz continuous in X, subject to growth conditions.Source: §2 · pp. 4–5
- A gradient-descent iteration requires forward and backward matrix-equation solves, which can be computationally intensive when m and n are large.The manuscript therefore analyzes dynamical low-rank approximations to accelerate gradient computation.Source: §2 · p. 6
3. The manifold of fixed-rank matrices and dynamical low-rank approximations
The section introduces the fixed-rank matrix manifold and its tangent-space projector, then recalls DLRA as the projected evolution of a matrix ODE. The projected state equation is treated under an assumption of a unique continuous solution.
3. The manifold of fixed-rank matrices and dynamical low-rank approximations
- A rank-r matrix is represented on the fixed-rank manifold, and its tangent vectors admit a unique parametrization used to describe the tangent space.Source: §3 · p. 6
- The tangent-space projector has a closed form for Y = USV⊤, with weighted inner products ensuring its orthogonality in the chosen matrix metric.Source: §3 · p. 6
- DLRA chooses the approximate trajectory’s derivative by projecting the full vector field onto the tangent space at the current low-rank state.The variational minimization is equivalent to a Galerkin condition and the projected matrix ODE.Source: §3 · p. 7
- The projected state equation is assumed to admit a unique continuous solution on the maximal interval for every square-integrable control.Source: §3 · p. 7
4. Dynamical low-rank approximation and optimality systems
This section develops DLRA strategies for approximating optimal-control gradients, covering OtP and PtO formulations. It formulates the PtO adjoint as a covariant evolution in the state’s moving tangent space and establishes gradient-error estimates that account for manifold curvature, with consistency in the full-rank limit.
4.1. Optimize-then-Project.
- OtP independently projects the state and adjoint equations onto their respective fixed-rank tangent spaces, approximating the true gradient rather than necessarily optimizing an identifiable auxiliary functional.For the stated objective, the adjoint approximation requires a small rank-r terminal perturbation because the zero matrix is not on the fixed-rank manifold; computations can initialize by briefly integrating the full adjoint and truncating it.Source: §4 · p. 7, §4.1 · p. 8, §4.1.1 · p. 8
4.1.1. Error estimate for the gradient.
- The OtP gradient-error estimate assumes bounded normal components of the state and adjoint dynamics and one-sided Lipschitz conditions; in the full-rank limit, only the adjoint terminal-condition error remains.The terminal contribution can be made arbitrarily small by tuning ε0.Source: §4.1.1 · pp. 8–9
4.2. Project-then-Optimize.
- PtO replaces the state equation with its DLRA approximation and the resulting reduced optimal-control problem admits a minimizer, but need not have a unique one.Even when the vector field is linear in state and control, the state-dependent projector makes the solution map nonlinear, so the reduced objective is generally nonconvex.Source: §4.2 · p. 10
4.2.1. Derivation of optimality conditions.
- In PtO, perturbations of the state trajectory remain in its tangent space, and the adjoint Q is constrained to that space as it evolves with the state.Its equation describes a covariant derivative: it determines the tangential evolution and is complemented by the explicit tangent-space constraint.Source: §4.2.1 · pp. 11–12
- The PtO adjoint Q(t) stays in T_Y(t)M_r, but its derivative generally leaves that tangent space as it rotates, requiring a normal correction in the covariant equation.The tangent-space formulation is suitable for DLRA-type integration, whereas the equivalent ambient formulation can be integrated with a standard matrix ODE solver.Source: §4.2.1 · pp. 12–13
4.2.2. Error estimate for the gradient.
- Lemma 4.13 bounds the normal component F⊥=(I−P_Y)F under Assumption 4.3, with σr(Y) denoting the smallest singular value.Source: §4.2.2 · p. 13
- Theorem 4.15 bounds PtO gradient error, with additional geometric terms tracking tangent-space rotation κeff(t) and εY relative to σr(Y(t)).The estimate assumes the local bound in Assumption 4.14; κeff(t) measures the rate of tangent-space rotation along the state trajectory.Source: §4.2.2 · pp. 13–14
- In the full-rank limit, εY and bεQ vanish and the projection is the identity, so the estimate implies ∇bJ(z)=∇eJ(z).Source: §4.2.2 · p. 14
5. Time Discretization
The section specifies low-rank forward and backward time discretizations, analyzes projected Euler for PtO’s covariant adjoint equation, and estimates computational costs.
5.1. Forward sweep.
- The forward sweep uses the BUG integrator to update column and row bases before computing the core by Galerkin projection; it is first order regardless of inner-solver accuracy.The paper uses BUG in its numerical experiments, while noting that higher-order and rank-adaptive alternatives are possible.Source: §5.1 · pp. 14–15
5.2. Backward sweep.
- OtP integrates its adjoint backward with BUG, whereas PtO uses a covariant adjoint ODE advanced by projected Euler from a zero terminal value in the state’s tangent space.For the application-specific PtO right-hand side, the curvature term involves inverse singular values and may induce stiffness; the experiments report no significant time-step restriction.Source: §5.2 · p. 15, §5.2.2 · p. 18
5.2.1. Time integration of a covariant ODE.
- Projected Euler advances tangent factors by projection, avoiding the need for state-factor time derivatives that are generally unavailable at discrete time instants.Its local truncation error is O(∆t^2), and its updates approximate the factor ODEs to second order.Source: §5.2.1 · pp. 15–16
- With Lipschitz G and bounded effective curvature, projected Euler converges globally at first order, with constants depending on Q’s acceleration, curvature, and projector derivatives.The stated convergence analysis assumes the exact state trajectory; transferring it rigorously to a discrete trajectory is left for future work.Source: §5.2.1 · p. 17
- Both covariant-ODE integrators avoid time-step restrictions from small singular values, while projected Euler’s error depends on effective curvature, which remains moderate as rank increases.The final-time error and curvature behavior are shown for ranks 4, 8, 16, and 32.Source: Figure 1 · p. 18, §5.2.1 · p. 18
5.3. Computational cost.
- Storing the state and either adjoint representation costs O(Nt(m+n+r)r), versus O(Nt mn) for a full-rank iteration; computational savings also require direct low-rank evaluation of F and its adjoint.For the stated affine vector-field form with sparse operators, the state term costs O(Np(m+n)r^2), compared with O(Np mn) for the ambient matrix ODE.Source: §5.3 · p. 19
6. Application to parabolic optimal control problems with random coefficients
This section formulates a random parabolic optimal-control problem as a matrix ODE and tests OtP and PtO low-rank approximations. The experiments report lower computational and storage complexity for both approaches, alongside differences in gradient and adjoint accuracy and optimization behavior.
6. Application to parabolic optimal control problems with random coefficients.
- The random parabolic model seeks a single deterministic open-loop control minimizing expected state-tracking error plus a penalty on control energy.The reduced functional is convex and has a unique minimizer under the stated assumptions.Source: §6 · pp. 19–20
6.1. Semidiscretization and optimality conditions.
- Finite-element spatial discretization and positive-weight quadrature yield sample-wise state equations and a backward adjoint system with zero terminal conditions.The optimality system also includes the control optimality condition; the stiffness matrices are symmetric.Source: §6.1 · p. 20
6.2. Affine random expansion and matrix reformulation.
- An affine expansion of the random diffusion coefficient lets the ensemble of discrete states be represented by a matrix ODE.Its columns correspond to individual state realizations.Source: §6.2 · pp. 20–21
6.3. Numerical experiments.
- Both DLRA approaches converge at rate Δt in gradient and adjoint errors before truncation error dominates, while PtO has larger errors for fixed Δt and rank.The adjoint-control operator averages adjoint states, mitigating approximation errors so gradients are typically more accurate than the corresponding adjoints; the PtO adjoint's larger normal component is associated with larger error constants.Source: §6.3 · pp. 22–23, Figure 2 · p. 22
- Even small ranks achieve small objective gaps, but PtO with r = 8 appears to stall at a local minimum although the original example is convex.The experiments also report that effective curvature does not explode as the smallest singular value approaches zero and remains well behaved as rank increases.Source: §6.3 · pp. 22–23, Figure 3 · p. 23, Figure 4 · p. 24
- Both OtP and PtO have complexity of order O((m + n)r), compared with O(m × n) for the full-rank approach, and require equal storage for state and adjoint variables.Figure 5 compares computational times for 30 gradient-descent iterations and memory footprints across problem sizes and ranks.Source: §6.3 · p. 23, Figure 5 · p. 24
7. Conclusion
The paper presents a framework for using dynamical low-rank approximation to accelerate optimal control of time-dependent matrix systems. Its analysis relates gradient-approximation accuracy to normal vector-field components, system stability, and effective curvature; extending the estimates to time discretization remains future work.
7. Conclusion
- Gradient-approximation accuracy in both OtP and PtO depends jointly on normal vector-field components, system stability, and geometric factors such as effective curvature.The authors identify extending the error estimates to account for time discretization as ongoing work.Source: §7 · p. 24
Appendix A. Adjoint of the map H 7→DPY [H]Z
This appendix derives the adjoint of the map H ↦ DP_Y[H]Z and records simplifications when its arguments lie in the tangent or normal spaces.
Appendix A. Adjoint of the map H 7→DPY [H]Z
- The derivation treats the two contributions to the adjoint separately, expressing the tangent-space components U_p and V_p in terms of H.It then gives a simplified expression when W lies in the tangent space at Y and notes a further simplification when Z is normal.Source: Appendix A · pp. 26–27
- For the first contribution, the appendix uses the cyclic properties of the trace operator to derive the adjoint expression.Source: Appendix A · p. 26
- When W lies in the tangent space at Y, the adjoint maps into that tangent space and has a corresponding tangent-space parametrization.The appendix also identifies an additional simplification when Z lies in the normal space.Source: Appendix A · p. 27