Melvin Leok

Wu Lin, Valentin Duruisseaux, Melvin Leok et al.

Riemannian submanifold optimization with momentum is computationally challenging because, to ensure that the iterates remain on the submanifold, we often need to solve difficult differential equations. Here, we simplify such difficulties for a class of sparse or structured symmetric positive-definite matrices with the affine-invariant metric. We do so by proposing a generalized version of the Riemannian normal coordinates that dynamically orthonormalizes the metric and locally converts the problem into an unconstrained problem in the Euclidean space. We use our approach to simplify existing approaches for structured covariances and develop matrix-inverse-free $2^\text{nd}$-order optimizers for deep learning with low precision by using only matrix multiplications. Code: https://github.com/yorkerlin/StructuredNGD-DL

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

This paper provides new results for control of complex flight maneuvers for a quadrotor unmanned aerial vehicle (UAV). The flight maneuvers are defined by a concatenation of flight modes or primitives, each of which is achieved by a nonlinear controller that solves an output tracking problem. A mathematical model of the quadrotor UAV rigid body dynamics, defined on the configuration space $\SE$, is introduced as a basis for the analysis. The quadrotor UAV has four input degrees of freedom, namely the magnitudes of the four rotor thrusts; each flight mode is defined by solving an asymptotic optimal tracking problem. Although many flight modes can be studied, we focus on three output tracking problems, namely (1) outputs given by the vehicle attitude, (2) outputs given by the three position variables for the vehicle center of mass, and (3) output given by the three velocity variables for the vehicle center of mass. A nonlinear tracking controller is developed on the special Euclidean group $\SE$ for each flight mode, and the closed loop is shown to have desirable closed loop properties that are almost global in each case. Several numerical examples, including one example in which the quadrotor recovers from being initially upside down and another example that includes switching and transitions between different flight modes, illustrate the versatility and generality of the proposed approach.

Melvin Leok, Tatiana Shingel

The numerical analysis of variational integrators relies on variational error analysis, which relates the order of accuracy of a variational integrator with the order of approximation of the exact discrete Lagrangian by a computable discrete Lagrangian. The exact discrete Lagrangian can either be characterized variationally, or in terms of Jacobi's solution of the Hamilton-Jacobi equation. These two characterizations lead to the Galerkin and shooting-based constructions for discrete Lagrangians, which depend on a choice of a numerical quadrature formula, together with either a finite-dimensional function space or a one-step method. We prove that the properties of the quadrature formula, finite-dimensional function space, and underlying one-step method determine the order of accuracy and momentum-conservation properties of the associated variational integrators. We also illustrate these systematic methods for constructing variational integrators with numerical examples.

Dmitry V. Zenkov, Melvin Leok, Anthony M. Bloch

This paper discusses Hamel's formalism and its applications to structure-preserving integration of mechanical systems. It utilizes redundant coordinates in order to eliminate multiple charts on the configuration space as well as nonphysical artificial singularities induced by local coordinates, while keeping the minimal possible degree of redundancy and avoiding integration of differential-algebraic equations.

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

This paper presents an analytical model and a geometric numerical integrator for a system of rigid bodies connected by ball joints, immersed in an irrotational and incompressible fluid. The rigid bodies can translate and rotate in three-dimensional space, and each joint has three rotational degrees of freedom. This model characterizes the qualitative behavior of three-dimensional fish locomotion. A geometric numerical integrator, refereed to as a Lie group variational integrator, preserves Hamiltonian structures of the presented model and its Lie group configuration manifold. These properties are illustrated by a numerical simulation for a system of three connected rigid bodies.

Melvin Leok, Tatiana Shingel

We introduce a novel technique for constructing higher-order variational integrators for Hamiltonian systems of ODEs. In particular, we are concerned with generating globally smooth approximations to solutions of a Hamiltonian system. Our construction of the discrete Lagrangian adopts Hermite interpolation polynomials and the Euler-Maclaurin quadrature formula, and involves applying collocation to the Euler-Lagrange equation and its prolongation. Considerable attention is devoted to the order analysis of the resulting variational integrators in terms of approximation properties of the Hermite polynomials and quadrature errors. A performance comparison is presented on a selection of these integrators.

Xuefeng Shen, Melvin Leok

In this paper, we consider exponential integrators for semilinear Poisson systems. Two types of exponential integrators are constructed, one preserves the Poisson structure, and the other preserves energy. Numerical experiments for semilinear Possion systems obtained by semi-discretizing Hamiltonian PDEs are presented. These geometric exponential integrators exhibit better long time stability properties as compared to non-geometric integrators, and are computationally more efficient than traditional symplectic integrators and energy-preserving methods based on the discrete gradient method.

Valentin Duruisseaux, Thai Duong, Melvin Leok et al.

Incorporating prior knowledge of physics laws and structural properties of dynamical systems into the design of deep learning architectures has proven to be a powerful technique for improving their computational efficiency and generalization capacity. Learning accurate models of robot dynamics is critical for safe and stable control. Autonomous mobile robots, including wheeled, aerial, and underwater vehicles, can be modeled as controlled Lagrangian or Hamiltonian rigid-body systems evolving on matrix Lie groups. In this paper, we introduce a new structure-preserving deep learning architecture, the Lie group Forced Variational Integrator Network (LieFVIN), capable of learning controlled Lagrangian or Hamiltonian dynamics on Lie groups, either from position-velocity or position-only data. By design, LieFVINs preserve both the Lie group structure on which the dynamics evolve and the symplectic structure underlying the Hamiltonian or Lagrangian systems of interest. The proposed architecture learns surrogate discrete-time flow maps allowing accurate and fast prediction without numerical-integrator, neural-ODE, or adjoint techniques, which are needed for vector fields. Furthermore, the learnt discrete-time dynamics can be utilized with computationally scalable discrete-time (optimal) control strategies.

Evan S. Gawlik, Melvin Leok

We derive a family of high-order, structure-preserving approximations of the Riemannian exponential map on several matrix manifolds, including the group of unitary matrices, the Grassmannian manifold, and the Stiefel manifold. Our derivation is inspired by the observation that if $Ω$ is a skew-Hermitian matrix and $t$ is a sufficiently small scalar, then there exists a polynomial of degree $n$ in $tΩ$ (namely, a Bessel polynomial) whose polar decomposition delivers an approximation of $e^{tΩ}$ with error $O(t^{2n+1})$. We prove this fact and then leverage it to derive high-order approximations of the Riemannian exponential map on the Grassmannian and Stiefel manifolds. Along the way, we derive related results concerning the supercloseness of the geometric and arithmetic means of unitary matrices.

Evan Gawlik, Melvin Leok

We construct interpolation operators for functions taking values in a symmetric space -- a smooth manifold with an inversion symmetry about every point. Key to our construction is the observation that every symmetric space can be realized as a homogeneous space whose cosets have canonical representatives by virtue of the generalized polar decomposition -- a generalization of the well-known factorization of a real nonsingular matrix into the product of a symmetric positive-definite matrix times an orthogonal matrix. By interpolating these canonical coset representatives, we derive a family of structure-preserving interpolation operators for symmetric space-valued functions. As applications, we construct interpolation operators for the space of Lorentzian metrics, the space of symmetric positive-definite matrices, and the Grassmannian. In the case of Lorentzian metrics, our interpolation operators provide a family of finite elements for numerical relativity that are frame-invariant and have signature which is guaranteed to be Lorentzian pointwise. We illustrate their potential utility by interpolating the Schwarzschild metric numerically.

Jeremy Schmitt, Tatiana Shingel, Melvin Leok

In this paper, we present a variational integrator that is based on an approximation of the Euler--Lagrange boundary-value problem via Taylor's method. This can viewed as a special case of the shooting-based variational integrator. The Taylor variational integrator exploits the structure of the Taylor method, which results in a shooting method that is one order higher compared to other shooting methods based on a one-step method of the same order. In addition, this method can generate quadrature nodal evaluations at the cost of a polynomial evaluation, which may increase its efficiency relative to other shooting-based variational integrators. A symmetric version of the method is proposed, and numerical experiments are conducted to exhibit the efficacy and efficiency of the method.

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

Euler-Lagrange equations and variational integrators are developed for Lagrangian mechanical systems evolving on a product of two-spheres. The geometric structure of a product of two-spheres is carefully considered in order to obtain global equations of motion. Both continuous equations of motion and variational integrators completely avoid the singularities and complexities introduced by local parameterizations or explicit constraints. We derive global expressions for the Euler-Lagrange equations on two-spheres which are more compact than existing equations written in terms of angles. Since the variational integrators are derived from Hamilton's principle, they preserve the geometric features of the dynamics such as symplecticity, momentum maps, or total energy, as well as the structure of the configuration manifold. Computational properties of the variational integrators are illustrated for several mechanical systems.

Xuefeng Shen, Melvin Leok

Rigid body dynamics on the rotation group have typically been represented in terms of rotation matrices, unit quaternions, or local coordinates, such as Euler angles. Due to the coordinate singularities associated with local coordinate charts, it is common in engineering applications to adopt the unit quaternion representation, and the numerical simulations typically impose the unit length condition using constraints or by normalization at each step. From the perspective of geometric structure-preserving, such approaches are undesirable as they are either computationally less efficient, or interfere with the preservation of other geometric properties of the dynamics. In this paper, we adopt the approach used in constructing Lie group variational integrators for rigid body dynamics on the rotation group to the representation in terms of unit quaternions. In particular, the rigid body dynamics is lifted to unit quaternions, and the Lie group structure of unit quaternions is used to represent tangent vector intrinsically, thereby avoiding the use of a Lagrange multiplier. A Lie group variational integrator in the unit quaternion representation is derived, and numerical results are presented.

Helen Parks, Melvin Leok

The dynamical motion of mechanical systems possesses underlying geometric structures, and preserving these structures in numerical integration improves the qualitative accuracy and reduces the long-time error of the simulation. For a single mechanical system, structure preservation can be achieved by adopting the variational integrator construction. This construction has been generalized to more complex systems involving forces or constraints as well as to the setting of Dirac mechanics. Variational integrators have recently been applied to interconnected systems in Parks and Leok (2017), which are an important class of practically useful mechanical systems whose description in terms of Dirac structures and Dirac mechanical systems was elucidated in Jacobs and Yoshimura (2014). Since these interconnected systems are modeled as a collection of subsystems with forces of interconnection, we revisit some of the properties of forced variational integrators. In particular, we derive a class of Dirac variational integrators with forces that exhibit preservation properties that are critical when applying variational integrators to the discretization of interconnected Dirac systems. We close with a discussion of ongoing and future research based on these findings.

Melvin Leok, Gautam Wilkins

Uncertainty quantification seeks to provide a quantitative means to understand complex systems that are impacted by parametric uncertainty. The polynomial chaos method is a computational approach to solve stochastic partial differential equations (SPDE) by projecting the solution onto a space of orthogonal polynomials of the stochastic variables and solving for the deterministic coefficients. Polynomial chaos can be more efficient than Monte Carlo methods when the number of stochastic variables is low, and the integration time is not too large. When performing long-term integration, however, achieving accurate solutions often requires the space of polynomial functions to become unacceptably large. This paper presents an alternative approach, where sets of empirical basis functions are constructed by examining the behavior of the solution for fixed values of the random variables. The empirical basis functions are evolved over time, which means that the total number can be kept small, even when performing long-term integration. We introduce this method of empirical chaos expansion, and apply it to a number of model equations, demonstrating that the computational time scales linearly with the final integration time. That is not the case for polynomial chaos in general, since achieving accuracy for long-term integration usually requires larger polynomial bases, causing a nonlinear scaling with the final integration time. We also present an analytical method that uses the dynamics of the SPDE to evolve the empirical basis functions and demonstrate how it can be applied to extend the validity of empirical basis functions in time without the need to sample additional realizations of the original SPDE.

Yiming Zhang, Sitong Liu, Ke Li et al.

Diffusion models are a leading paradigm for data generation, but training-free editing typically re-runs the full denoising trajectory for every edit strength, making iterative refinement expensive. To address this issue, we instead edit near the data manifold, where small local updates can replace repeated re-synthesis. To enable this, we estimate a local manifold tangent space directly from perturbed samples and prove that this sample-based estimator closely approximates the true tangent. Building on this guarantee, we devise a Jacobian-free algorithm that constructs a tangent frame via small perturbations to the initial noise and alternates small tangent moves with diffusion-based projections. Updates within this frame follow principled on-manifold directions while suppressing off-manifold drift, enabling fine-grained edits without full re-diffusion or additional training. Edit strength is controlled by the number of steps for rapid, continuous adjustments that preserve fidelity and plug into existing samplers. Empirically, the resulting tangent directions yield smooth, semantic unsupervised traversals and effective CLIP-guided optimization, demonstrating practical interactive continuous editing.

James Hall, Melvin Leok

In this paper, we present a new variational integrator for problems in Lagrangian mechanics. Using techniques from Galerkin variational integrators, we construct a scheme for numerical integration that converges geometrically, and is symplectic and momentum preserving. Furthermore, we prove that under appropriate assumptions, variational integrators constructed using Galerkin techniques will yield numerical methods that are in a certain sense optimal, converging at the same rate as the best possible approximation in a certain function space. We further prove that certain geometric invariants also converge at an optimal rate, and that the error associated with these geometric invariants is independent of the number of steps taken. We close with several numerical examples that demonstrate the predicted rates of convergence.

Harsh Choudhary, Chandan Gupta, Vyacheslav Kungurtsev et al.

Many important physical systems can be described as the evolution of a Hamiltonian system, which has the important property of being conservative, that is, energy is conserved throughout the evolution. Physics Informed Neural Networks and in particular Hamiltonian Neural Networks have emerged as a mechanism to incorporate structural inductive bias into the NN model. By ensuring physical invariances are conserved, the models exhibit significantly better sample complexity and out-of-distribution accuracy than standard NNs. Learning the Hamiltonian as a function of its canonical variables, typically position and velocity, from sample observations of the system thus becomes a critical task in system identification and long-term prediction of system behavior. However, to truly preserve the long-run physical conservation properties of Hamiltonian systems, one must use symplectic integrators for a forward pass of the system's simulation. While symplectic schemes have been used in the literature, they are thus far limited to situations when they reduce to explicit algorithms, which include the case of separable Hamiltonians or augmented non-separable Hamiltonians. We extend it to generalized non-separable Hamiltonians, and noting the self-adjoint property of symplectic integrators, we bypass computationally intensive backpropagation through an ODE solver. We show that the method is robust to noise and provides a good approximation of the system Hamiltonian when the state variables are sampled from a noisy observation. In the numerical results, we show the performance of the method concerning Hamiltonian reconstruction and conservation, indicating its particular advantage for non-separable systems.

Harsh Choudhary, Vyacheslav Kungurtsev, Chandan Gupta et al.

Hamiltonian Neural Networks (HNNs) integrate physical priors into neural models by learning a system's Hamiltonian, improving generalization and sample efficiency. Identifying the system Hamiltonian from noisy observations of state variables is a challenging task. For simulations to faithfully reflect the long-term behavior of Hamiltonian systems, especially energy conservation, it is essential to use symplectic integrators, which preserve the system's geometric structure. This fidelity comes at a cost: implicit symplectic integrators are more computationally intensive and make backpropagation through the ODE solver non-trivial. However, by leveraging the fact that symplectic discretizations of the adjoint system yield the same sensitivities associated by backpropagation, we obtain an efficient method of training the Neural Network parameters. In our work, we explore this alternate method of HNN training under noisy observation of trajectories with our HNN model based on an implicit symplectic integrator. Computationally, a predictor-corrector based ODE solver and fixed point iteration help to mitigate the computational cost of the implicit timestepping, resulting in more efficient generation of gradient updates. We showcase the numerical advantage, in experiments, in system identification and energy preservation on a range of non-separable, chaotic systems and the efficient computation and memory complexity of our method. We also observe that the post-processing of the learned Hamiltonian using backward error analysis yields a modified Hamiltonian that is a more accurate approximation of the true Hamiltonian without the need to use more accurate discretizations of the flow map.

Kehan Long, Vikas Dhiman, Melvin Leok et al.

This paper considers safe control synthesis for dynamical systems with either probabilistic or worst-case uncertainty in both the dynamics model and the safety constraints. We formulate novel probabilistic and robust (worst-case) control Lyapunov function (CLF) and control barrier function (CBF) constraints that take into account the effect of uncertainty in either case. We show that either the probabilistic or the robust (worst-case) formulation leads to a second-order cone program (SOCP), which enables efficient safe and stable control synthesis. We evaluate our approach in PyBullet simulations of an autonomous robot navigating in unknown environments and compare the performance with a baseline CLF-CBF quadratic programming approach.

Jeremy M. Schmitt, Melvin Leok

Discrete Hamiltonian variational integrators are derived from Type II and Type III generating functions for symplectic maps, and in this paper we establish a variational error analysis result that relates the order of accuracy of the associated numerical methods with the extent to which these generating functions approximate the exact discrete Hamiltonians. We also introduce the notion of an adjoint discrete Hamiltonian, and relate it to the adjoint of the associated symplectic integrator. We show that when constructing discrete Lagrangians and discrete Hamiltonians using the Taylor variational integrator approach, the same underlying one-step method and quadrature rule does not necessarily lead to the same symplectic integrator, and the same observation holds when developing variational integrators based on averaging techniques. Numerical experiments also indicate that the resonance behavior of variational integrators also depend on the type of generating functions used, and we relate this resonance behavior to the ill-posedness of the boundary-value problems used to define the exact discrete Lagrangian and exact discrete Hamiltonian.

Evan S. Gawlik, Melvin Leok

We study schemes for interpolating functions that take values in the special orthogonal group $SO(n)$. Our focus is on interpolation schemes obtained by embedding $SO(n)$ in a linear space, interpolating in the linear space, and mapping the result onto $SO(n)$ via the closest point projection. The resulting interpolants inherit both the order of accuracy and the regularity of the underlying interpolants on the linear space. The values and derivatives of the interpolants admit efficient evaluation via either explicit formulas or iterative algorithms, which we detail for two choices of embeddings: the embedding of $SO(n)$ in the space of $n \times n$ matrices and, when $n=3$, the identification of $SO(3)$ with the set of unit quaternions. Along the way, we point out a connection between these interpolation schemes and geodesic finite elements. We illustrate the utility of these interpolation schemes by numerically computing minimum acceleration curves on $SO(n)$, a task which is handled naturally with $SO(n)$-valued finite elements having $C^1$-continuity.

Evan S. Gawlik, Melvin Leok

We derive iterative methods for computing the Fréchet derivative of the map which sends a full-rank matrix $A$ to the factor $U$ in its polar decomposition $A=UH$, where $U$ has orthonormal columns and $H$ is Hermitian positive definite. The methods apply to square matrices as well as rectangular matrices having more rows than columns. Our derivation relies on a novel identity that relates the Fréchet derivative of the polar decomposition to the matrix sign function $\mathrm{sign}(X) = X (X^2)^{-1/2}$ applied to a certain block matrix $X$.

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

This paper provides nonlinear tracking control systems for a quadrotor unmanned aerial vehicle (UAV) that are robust to bounded uncertainties. A mathematical model of a quadrotor UAV is defined on the special Euclidean group, and nonlinear output-tracking controllers are developed to follow (1) an attitude command, and (2) a position command for the vehicle center of mass. The controlled system has the desirable properties that the tracking errors are uniformly ultimately bounded, and the size of the ultimate bound can be arbitrarily reduced by control system parameters. Numerical examples illustrating complex maneuvers are provided.

Melvin Leok, Jingjing Zhang

We consider the continuous and discrete-time Hamilton's variational principle on phase space, and characterize the exact discrete Hamiltonian which provides an exact correspondence between discrete and continuous Hamiltonian mechanics. The variational characterization of the exact discrete Hamiltonian naturally leads to a class of generalized Galerkin Hamiltonian variational integrators, which include the symplectic partitioned Runge-Kutta methods. We also characterize the group invariance properties of discrete Hamiltonians which lead to a discrete Noether's theorem.

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

This paper presents an analytical model and a geometric numerical integrator for a rigid body connected to an elastic string, acting under a gravitational potential. Since the point where the string is attached to the rigid body is displaced from the center of mass of the rigid body, there exist nonlinear coupling effects between the string deformation and the rigid body dynamics. A geometric numerical integrator, refereed to as a Lie group variational integrator, is developed to numerically preserve the Hamiltonian structure of the presented model and its Lie group configuration manifold. These properties are illustrated by a numerical simulation.

Eugene G. Fahnestock, Taeyoung Lee, Melvin Leok et al.

We present a combination of tools which allows for investigation of the coupled orbital and rotational dynamics of two rigid bodies with nearly arbitrary shape and mass distribution, under the influence of their mutual gravitational potential. Methods for calculating that mutual potential and resulting forces and moments for a polyhedral body representation are simple and efficient. Discrete equations of motion, referred to as the Lie Group Variational Integrator (LGVI), preserve the structure of the configuration space, SE(3), as well as the geometric features represented by the total energy and the total angular momentum. The synthesis of these approaches allows us to simulate the full two body problem accurately and efficiently. Simulation results are given for two octahedral rigid bodies for comparison with other integration methods and to show the qualities of the results thus obtained. A significant improvement is seen over other integration methods while correctly capturing the interesting effects of strong orbit and attitude dynamics coupling, in multiple scenarios.

Taeyoung Lee, Melvin Leok, N. Harris McClamroch

We develop the equations of motion for full body models that describe the dynamics of rigid bodies, acting under their mutual gravity. The equations are derived using a variational approach where variations are defined on the Lie group of rigid body configurations. Both continuous equations of motion and variational integrators are developed in Lagrangian and Hamiltonian forms, and the reduction from the inertial frame to a relative coordinate system is also carried out. The Lie group variational integrators are shown to be symplectic, to preserve conserved quantities, and to guarantee exact evolution on the configuration space. One of these variational integrators is used to simulate the dynamics of two rigid dumbbell bodies.

Mathieu Desbrun, Anil N. Hirani, Melvin Leok et al.

We present a theory and applications of discrete exterior calculus on simplicial complexes of arbitrary finite dimension. This can be thought of as calculus on a discrete space. Our theory includes not only discrete differential forms but also discrete vector fields and the operators acting on these objects. This allows us to address the various interactions between forms and vector fields (such as Lie derivatives) which are important in applications. Previous attempts at discrete exterior calculus have addressed only differential forms. We also introduce the notion of a circumcentric dual of a simplicial complex. The importance of dual complexes in this field has been well understood, but previous researchers have used barycentric subdivision or barycentric duals. We show that the use of circumcentric duals is crucial in arriving at a theory of discrete exterior calculus that admits both vector fields and forms.

Melvin Leok

We introduce generalized Galerkin variational integrators, which are a natural generalization of discrete variational mechanics, whereby the discrete action, as opposed to the discrete Lagrangian, is the fundamental object. This is achieved by approximating the action integral with appropriate choices of a finite-dimensional function space that approximate sections of the configuration bundle and numerical quadrature to approximate the integral. We discuss how this general framework allows us to recover higher-order Galerkin variational integrators, asynchronous variational integrators, and symplectic-energy-momentum integrators. In addition, we will consider function spaces that are not parameterized by field values evaluated at nodal points, which allows the construction of Lie group, multiscale, and pseudospectral variational integrators. The construction of pseudospectral variational integrators is illustrated by applying it to the (linear) Schrodinger equation. G-invariant discrete Lagrangians are constructed in the context of Lie group methods through the use of natural charts and interpolation at the level of the Lie algebra. The reduction of these G-invariant Lagrangians yield a higher-order analogue of discrete Euler-Poincare reduction. By considering nonlinear approximation spaces, spatio-temporally adaptive variational integrators can be introduced as well.