\section{Introduction}\label{sec:intro}
Discovering hidden rules from sequential observations has been an essential topic in machine learning, with a large variety of applications such as physics simulation \citep{sanchez2020learning}, autonomous driving \citep{diehl2019graph}, ECG analysis \citep{golany2021ecg} and event analysis \citep{chen2020learning}, to name a few. A standard scheme is to consider sequential data at each timestamp to be holistic and homogeneous under some ideal assumptions (i.e., only the temporal behavior of one entity is involved in a sequence), under which data/observation is treated as a collection of slices at a different time from a unified system.
A series of sequential learning models fall into this category, including variants of recurrent neural networks (RNNs) \citep{cho2014properties,hochreiter1997long}, neural differential equations (DEs) \citep{chen2018neural,kidger2020neural,rusch2021unicornn,zhu2021neural} and spatial/temporal attention-based approaches \citep{vaswani2017attention,fan2019multi,song2017end}.
These variants fit well into the scenarios agreeing with the aforementioned assumptions and are proved effective in learning or modeling for relatively simple applications with clean data sources. 

In the real world, a system may not only describe a single and holistic entity but also consist of several \emph{distinguishable} interacting but simple subsystems, where each subsystem corresponds to a physical entity. For example, we can think of the movement of a solar system as the mixture of distinguishable subsystems of the sun and surrounding planets, while interactions between these celestial bodies over time are governed by the laws of gravity. Back centuries ago, physicists and astronomers made enormous efforts to discover the rule of celestial movements from the records of every single body and eventually delivered the neat yet elegant differential equations (DEs) depicting principles of moving bodies and interactions therein. Likewise, nowadays, researchers also developed a series of machine learning models for sequential data with distinguishable partitions \citep{qin2017dual}. Two widely adopted strategies for learning the interactions between subsystems are graph neural networks \citep{iakovlev2020learning,ha2021unraveling,kipf2018neural,yildiz2022learning,xhonneux2020continuous} and attention mechanism \citep{vaswani2017attention,lu2019understanding,goyal2019recurrent}, while the interactions are typically encoded with ``messages'' between nodes and pair-wise ``attention scores'', respectively.

It is worth noting an even more difficult scenario: 
\begin{itemize}
    \item \emph{The data/observation is so cluttered that cannot be readily distinguished into separate parts}. 
\end{itemize}
This can be either due to the way of data collection (e.g., videos consisting of multiple objects) or because there are no explicit physical entities originally (e.g., weather time series). 
To tackle this, a fair assumption can be introduced that complex observations can be decoupled into several relatively independent modules in the feature space, where each module corresponds to a \emph{latent entity}. Latent entities may not have exact physical meanings, but learning procedures can greatly benefit from such decoupling, as this assumption can be viewed as strong regularization to the system. This assumption has been successfully incorporated in several models for learning from \emph{regularly} sampled sequential data by emphasizing ``independence'' to some extent between channels or groups in the feature space \citep{li2018independently,yu2020rhyrnn,goyal2019recurrent,madan2021fast}. 
Another successful counterpart in parallel benefiting from this assumption is transformer \citep{vaswani2017attention} which stacks multiple layers of self-attention and point-wise feedforward networks. In transformers, each attention head can be viewed as a relatively independent module, and interaction happens throughout the head re-weighting procedure following the attention scores. \citet{lu2019understanding} presented an interpretation from a dynamic point of view by regarding a basic layer in the transformer as one step of integration governed by differential equations derived from interacting particles. \citet{vuckovic2020mathematical} extended this interpretation with more solid mathematical support by viewing the forward pass of the transformer as applying successive Markov kernels in a particle-based dynamic system. 

We note, however, despite the ubiquity of this setting, there is barely any previous investigation focusing on learning for \emph{irregularly sampled} and \emph{cluttered} sequential data. The aforementioned works either fail to handle the irregularity \citep{goyal2019recurrent,li2018independently} or neglect the independence/modularity assumption in the latent space \citep{chen2018neural,kidger2020neural}. 
In this paper, inspired by recent advances of neural controlled dynamics \citep{kidger2020neural} and novel interpretation of attention mechanism \citep{vuckovic2020mathematical}, we take a step to propose an effective approach addressing this problem under the dynamic setting. To this end, our approach explicitly learned to decouple a complex system into several latent sub-systems and utilizes an additional meta-system capturing the evolution of interactions over time.
Specifically, taking into account the meta-system capturing interactions evolving in a constrained set (e.g., simplex), we further characterized such interactions using projected differential equations (ProjDEs) with neural-friendly projection operators. We argued our \textbf{contributions} as follows:
\begin{itemize}
    \item We provide a novel modeling strategy for sequential data from a system decoupling perspective;
    \item We propose a novel and natural interpretation of evolving interactions as a ProjDE-based meta-system, with insights into projection operators in the sense of Bregman divergence;
    \item Our approach is parameter-insensitive and more compatible with other modules and data, thus being flexible to be integrated into various tasks.
\end{itemize}
Extensive experiments were conducted on either regularly or irregularly sampled sequential data, including both synthetic and real-world settings. It was observed that our approach achieved prominent performance compared to the state-of-the-art on a wide spectrum of tasks.


\section{Related Work}

\paragraph{Sequential learning.} Traditionally, learning with sequential data can be performed using variants of recurrent neural networks (RNNs) \citep{hochreiter1997long,cho2014properties,li2018independently} under the Markov setting. While such RNNs are generally designed for regular sampling frequency, a more natural line of counterparts lies in the continuous time domain allowing irregularly sampled time series as input. As such, a variety of RNN-based methods are developed by introducing exponential decay on observations \citep{che2018recurrent,mei2017neural}, incorporating an underlying Gaussian process \citep{li2016scalable,futoma2017learning}, or integrating some latent evolution under ODEs \citep{rubanova2019latent,de2019gru}. A seminal work interpreting forward passing in neural networks as an integration of ODEs was proposed in \citet{chen2018neural}, followed by a series of relevant works \citep{liu2019neural,li2020scalable,dupont2019augmented}. As integration over ODEs allows for arbitrary step length, it is natural modeling of irregular time series and proved powerful in many machine learning tasks (e.g., bioinformatics \citep{golany2021ecg}, physics \citep{nardini2021learning} and computer vision \citep{park2021vid}). \cite{kidger2020neural} studied a more effective way of injecting observations into the system via a mathematical tool called Controlled differential Equation, achieving state-of-the-art performance on several benchmarks. Some variants of neural ODEs have also been extended to discrete structure \citep{chamberlain2021grand,xhonneux2020continuous} and non-Euclidean setting \citep{chamberlain2021beltrami}.

\paragraph{Learning with independence.} Independence or modular property serves as strong regularization or prior in some learning tasks under static setting \citep{wang2020learning,liu2020independence}. In the sequential case, some early attempts over RNNs emphasized implicit ``independence'' in the feature space between dimensions or channels \citep{li2018independently,yu2020rhyrnn}. As independence assumption commonly holds in vision tasks (with distinguishable objects), \citet{pang2020complex,li2020hoi} proposed video understanding schemes by decoupling the spatiotemporal patterns. For a more generic case where the observations are collected without any prior, \citet{goyal2019recurrent} devised a sequential learning scheme called recurrent independence mechanism (RIM), 
and its generalization ability was extensively studied in \citet{madan2021fast}. \citet{lu2019understanding} investigated self-attention mechanism \citep{vaswani2017attention} and interpreted it as a nearly independent multi-particle system with interactions therein. \citet{vuckovic2020mathematical} further provided more solid mathematical analysis with the tool of Markov kernel. The study of such a mechanism in the dynamical setting was barely observed.

\textbf{Learning dynamics under constraints.} It is practically significant as a series of real-world systems evolve within some manifolds, such as fluid \cite{vinuesa2022enhancing}, coarse-grained dynamics \cite{kaltenbach2020incorporating}, and molecule modeling \cite{chmiela2020accurate}. While some previous research incorporates constraints from a physical perspective \cite{kaltenbach2020incorporating,linot2020deep}, an emerging line is empowered by machine learning to integrate or even discover the constraints \cite{kolter2019learning,lou2020neural,goldt2020modeling}. To ensure a system evolves in constraints, efficient projections or pseudo-projections are required, about which Bregman divergence provides rich insights \cite{martins2016softmax,krichene2015efficient,lim2016efficient}. Despite these results, to our best knowledge, there is barely any related investigation about neural-friendly projections.


\section{Methodology}

\subsection{Background}\label{sec:background}
In this section, we briefly review three aspects related to our approach. Our approach is built upon the basic sub-system derived from \emph{Neural Controlled Dynamics} \citep{kidger2020neural}, while the interactions are modeled at an additional meta-system analogous to \emph{Self-attention} \citep{lu2019understanding,vuckovic2020mathematical}, and further interpreted and generalized using the tool of \emph{Projected Differential Equations} \citep{dupuis1993dynamical}.

\paragraph{Neural Controlled Dynamics.} Continuous-time dynamics can be expressed using differential equations $\mathbf{z}'(t)=d\mathbf{z}/dt =f(\mathbf{z}(t),t)$, where $\mathbf{z}\in\mathbb{R}^d$ and $t$ are a $d$-dimension state and the time, respectively. Function $f:\mathbb{R}^d\times\mathbb{R}_+\rightarrow \mathbb{R}^d$ governs the evolution of the dynamics. Given the initial state $\mathbf{z}(t_0)$, the state at any time $t_1$ can be evaluated with:
\begin{equation}\label{eq: ode}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1}f(\mathbf{z}(s),s)\mathrm{d}s
\end{equation}
In practice, we aim at learning the dynamics from a series of observations or controls $\{\mathbf{x}(t_k)\in\mathbb{R}^b | k=0,1,...\}$ by parameterizing the dynamics with $f_\theta(\cdot)$ where $\theta$ is the unknown parameter to be learned. Thus, a generic dynamics incorporating outer signals $\mathbf{x}$ can be written as:
\begin{equation}\label{eq: node}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1}f_\theta(\mathbf{z}(s),\mathbf{x}(s),s)\mathrm{d}s
\end{equation} 
Rather than directly injecting $\mathbf{x}$ as in Eq.~(\ref{eq: node}), Neural Controlled Differential Equation (CDE) proposed to deal with outer signals with a Riemann–Stieltjes integral \citep{kidger2020neural}:
\begin{equation}\label{eq: ncde}
    \mathbf{z}(t_1)=\mathbf{z}(t_0)+\int_{t_0}^{t_1} \mathbf{F}_\theta(\mathbf{z}(s))\mathbf{x}'(s)\mathrm{d}s
\end{equation}
where $\mathbf{F}_\theta:\mathbb{R}^d \rightarrow \mathbb{R}^{d\times b}$ is a learnable vector field and $\mathbf{x}'(s)=\mathrm{d}\mathbf{x}/\mathrm{d}s$ is the derivative of signal $\mathbf{x}$ w.r.t. time $s$, thus ``$\mathbf{F}_\theta(\mathbf{z}(s))\mathbf{x}'(s)$'' is a matrix-vector multiplication. During implementation, \citet{kidger2020neural} argued that a simple cubic spline interpolation on $\mathbf{x}$ allows dense calculation of $\mathbf{x}'(t)$ at any time $t$ and exhibits promising performance. In \cite{kidger2020neural}, it is also mathematically shown that incorporating observations/controls following Eq. (\ref{eq: ncde}) is with greater representation ability compared to Eq. (\ref{eq: node}), hence achieving state-of-the-art performance on several public tasks.

\paragraph{Self-attention.} It is argued in \citet{lu2019understanding,vuckovic2020mathematical} that a basic unit in Transformer \citep{vaswani2017attention} with self-link consisting of one self-attention layer and point-wise feedforward layer amounts to simulating a multi-particle dynamical system. Considering such a layer with $n$ attention-heads (corresponding to $n$ particles), given an attention head index $i\in\{1,2,...,n\}$, the update rule of the $i$th unit at depth $l$ reads:
\begin{subequations}\label{eq: att}
    \begin{alignat}{2}                  &\tilde{\mathbf{z}}_{l,i}=\mathbf{z}_{l,i}+\mathrm{MHAtt}_{W^l_{\mathrm{att}}}\left(\mathbf{z}_{l,i},\left[\mathbf{z}_{l,1},...,\mathbf{z}_{l,n}\right]\right) \label{eq: att_dif} \\
    &\mathbf{z}_{l+1,i}=\tilde{\mathbf{z}}_{l,i}+\mathrm{FFN}_{W^l_{\mathrm{ffn}}}\left(\tilde{\mathbf{z}}_{l,i}\right) \label{eq: att_conv}
    \end{alignat}
\end{subequations}
where $\mathrm{MHAtt}_{W^l_{\mathrm{att}}}$ and $\mathrm{FFN}_{W^l_{\mathrm{ffn}}}$ are multi-head attention layer and feedforward layer with parameters $W^l_{\mathrm{att}}$ and $W^l_{\mathrm{ffn}}$, respectively. Eq.~(\ref{eq: att}) can then be interpreted as an interacting multi-particle system:
\begin{equation}\label{eq: att_overall}
    \frac{\mathrm{d}\mathbf{z}_i(t)}{\mathrm{d}t}=F(\mathbf{z}_i(t),[\mathbf{z}_1(t),...,\mathbf{z}_n(t)],t)+G(\mathbf{z}_i(t))
\end{equation}
where function $F$ corresponding to Eq.~(\ref{eq: att_dif}) represents the diffusion term and $G$ corresponding to Eq.~(\ref{eq: att_conv}) stands for the convection term. Notably, the attention score obtained via $\mathrm{softmax}$ in Eq.~(\ref{eq: att_dif}) is regarded as a Markov kernel. 
Readers are referred to \citet{lu2019understanding,vuckovic2020mathematical} for more details.

\paragraph{Projected DEs.} It is a tool depicting the behavior of dynamics where solutions are constrained within a (convex) set. Concretely, given a closed polyhedral $\mathcal{K}\subset \mathbb{R}^n$ and a mapping $H:\mathcal{K}\rightarrow \mathbb{R}^n$, we can introduce an operator $\Pi_{\mathcal{K}}: \mathbb{R}^n\times\mathcal{K} \rightarrow \mathbb{R}^n$ which is defined by means of directional derivatives as:
\begin{equation}\label{eq: derivative}
    \Pi_{\mathcal{K}}(\mathbf{a},H(\mathbf{a}))=\lim_{\alpha\rightarrow 0_+}\frac{P_{\mathcal{K}}(\mathbf{a}+\alpha H(\mathbf{a}))-\mathbf{a}}{\alpha}
\end{equation}
where $P_{\mathcal{K}}(\cdot)$ is a projection onto $\mathcal{K}$ in terms of Euclidean distance:
\begin{equation}\label{eq:proj-l2}
    \lVert P_{\mathcal{K}}(\mathbf{a})-\mathbf{a}\rVert _2=\inf_{\mathbf{y}\in\mathcal{K}}\lVert\mathbf{y}-\mathbf{a}\rVert_2
\end{equation}
Intuitively, Eq.~(\ref{eq: derivative}) pictures the dynamics of $\mathbf{a}$ driven by function $H$, but constrained within $\mathcal{K}$. Whenever $\mathbf{a}$ reaches beyond $\mathcal{K}$, it would be projected back using Eq.~(\ref{eq:proj-l2}). By extending Eq.~(\ref{eq: derivative}), \cite{dupuis1993dynamical,zhang1995stability} considered the projected differential equations as follows:
\begin{equation}
    \frac{\mathrm{d}\mathbf{a}(t)}{\mathrm{d}t}=\Pi_{\mathcal{K}}(\mathbf{a},H(\mathbf{a}))
\end{equation}
which allows for discontinuous dynamics on $\mathbf{a}$.

\subsection{Learning to decouple}\label{sec:method}
Our method is built upon the assumption that cluttered sequential observations are composed of several relatively independent sub-systems and, therefore, explicitly learns to decouple them as well as to capture the mutual interactions with a meta-system in parallel.
Let the cluttered observations/controlls be $\mathbf{c}(t)\in\mathbb{R}^k$ at time $t$ for $t=1,...,T$, where $T$ is the time horizon. We employ $k$ distinct mappings with learnable parameters (e.g., MLP) to obtain respective controls to each sub-system: $\mathbf{x}_i(t)=p_i(\mathbf{c}(t))\in\mathbb{R}^m$ for $i=1,...,n$. A generic dynamics of the proposed method can be written as:
\label{sec:method}
\begin{subequations}\label{eq:dec}
    \begin{alignat}{2}
        \frac{\mathrm{d}\mathbf{z}_i(t)}{\mathrm{d}t}&=f_{i}\left(\mathbf{z}_i(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t)\right],\mathbf{x}_i(t),\mathbf{a}(t)\right) \label{eq:dec_sub}\\
        \frac{\mathrm{d}\textbf{a}(t)}{\mathrm{d}t}&=\Pi_\mathcal{S}\left(\mathbf{a}(t),g(\mathbf{a}(t),\left[\mathbf{z}_1
        (t),...,\mathbf{z}_n(t)\right])\right) \label{eq:dec_att}
    \end{alignat}
\end{subequations}
where Eq.~(\ref{eq:dec_sub}) and Eq.~(\ref{eq:dec_att}) refer to the $i$th sub-system describing the evolution of a single latent entity and meta-system depicting the interactions, respectively. $\mathbf{z}_i(t)\in\mathbb{R}^q$ is the hidden state for the $i$th subsystem, and $\mathbf{a}$ is a tensor governs the dynamics of the interactions. Here $\Pi_\mathcal{S}(\cdot)$ is a projection operator, which projects the evolving trajectory into set $\mathcal{S}$. We introduce such an operator as it is assumed that interactions among latent entities should be constrained following some latent manifold structure. $f_i(\cdot)$ and $g(\cdot)$ are both learnable functions and also the essential roles for capturing the underlying complex dynamics.

\begin{remark}
It is seen the projection operator $\Pi_\mathcal{S}(\cdot)$ and the set $\mathcal{S}$ play important roles in Eq.~(\ref{eq:dec_att}). For $\Pi_\mathcal{S}(\cdot)$, while previous works of ProjDEs only consider L2-induced projection, we propose novel interpretation and extension under Bregman divergence. For $\mathcal{S}$, we consider a probabilistic simplex following the setting in \citet{lu2019understanding,vuckovic2020mathematical}, though it can be any polyhedral.
\end{remark}

According to Eq.~(\ref{eq:dec}), we fully decouple a complex system into several components. Although we found some decoupling counterparts in the context of RNNs \citep{li2018independently,yu2020rhyrnn} and attention-like mechanism \citep{lu2019understanding,goyal2019recurrent}, their decoupling could not be applied to our problem.
We elaborate on the details of implementing Eq.~(\ref{eq:dec}) in the following.

\paragraph{Learning sub-systems.} Sub-systems corresponding to the latent entities seek to model relatively independent dynamics separately. Specifically, we employ the way of integrating $\mathbf{x}_i$s into Eq.~(\ref{eq:dec_sub}) in a controlled dynamical fashion as in the state-of-the-art method \citep{kidger2020neural}:
\begin{equation}\label{eq: sub_generic}
    \mathrm{d}\mathbf{z}_i(t)= \mathbf{F}_i \left( \mathbf{z}_i(t),\mathbf{a}(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t) \right] \right) \mathrm{d}\mathbf{x}_i(t)
\end{equation}
where $\mathbf{F}_i(\cdot)\in\mathbb{R}^{q\times m}$ is a learnable vector field. Concretely, if we let $\mathbf{z}(t)=\left[\mathbf{z}_i(t),...,\mathbf{z}_n(t) \right]$ be the tensor collecting all sub-systems, the $i$th sub-system in a self-attention fashion reads:
\begin{equation}\label{eq: sub}
    \mathrm{d}\mathbf{z}_i(t)= \mathbf{F}( \left[ \mathbf{A}(t) \cdot \mathbf{z}(t) \right] _i) \mathrm{d}\mathbf{x}_i(t)
\end{equation}
where $[\cdot]_i$ takes the $i$th slice from a tensor. Note timestamp $t$ can be arbitrary, resulting in irregularly sampled sequential data. To address this, we follow the strategy in \citet{kidger2020neural} by performing cubic spline interpolation on $\mathbf{x}_i$ over observed timestamp $t$, resulting in $\mathbf{x}_i(t)$ at dense time $t$. Note that for all sub-systems, different from Eq.~(\ref{eq: sub_generic}) we utilize an identical function/network $\mathbf{F}(\cdot)$ as in Eq.~(\ref{eq: sub}), but with different control sequence $\mathbf{x}_i(t)=p_i(\mathbf{c}(t))$. Since in our implementation, $p_i(\cdot)$ is a lightweight network such as MLP, this can significantly reduce the parameter size.

\paragraph{Learning interactions.} In our approach, interactions between latent entities are modeled separately as another meta-system. This is quite different from some related methods \citep{lu2019understanding, vuckovic2020mathematical} where sub-systems and interactions are treated as one holistic step of forward integration. For the meta-system describing the interactions in Eq.~(\ref{eq:dec_att}), two essential components are involved: domain $\mathcal{S}$ and the projection operator $\Pi$. In the context of ProjDEs, a system is constrained as $\mathbf{a}(t)\in\mathcal{S}$ for any $t$. In terms of interactions, a common choice of $\mathcal{S}$ is the stochastic simplex which can be interpreted as a transition kernel \citep{vuckovic2020mathematical}. We allow follow this setting by defining $\mathcal{S}$ be a row-wise stochastic $(n-1)$-simplices:
\begin{equation}\label{eq: simplex}
    \mathcal{S}\triangleq\{\mathbf{A}\in\mathbb{R}^{n\times n}|\mathbf{A1=1},\mathbf{A}_{ij}\geq 0\}
\end{equation}
where $\mathbf{1}$ is a vector with all $1$ entries. $\mathbf{A}=\mathrm{mat}(\mathbf{a})$ is a $n\times n$ matrix. In the sequel, we will use the notation $\mathbf{A}$ throughout.
Thus the meta-system capturing the interactions can be implemented as follows:
\begin{equation}\label{eq: att_implement}
    \frac{d\mathbf{A}(t)}{dt}=\Pi_\mathcal{S}\left(\mathbf{A}(t),g(\mathbf{A}(t),\left[\mathbf{z}_1(t),...,\mathbf{z}_n(t)\right])\right)
\end{equation}
For the projection operator, we consider two versions shown in Eq.~(\ref{eq:projs}). In Eq.~(\ref{eq:pj_soft}), we give a row-wise projection onto the $(n-1)$-simplex with entropic regularization \citep{amos2019differentiable}, which has a well-known closed-form solution $\mathrm{softmax}(\cdot)$ appearing in attention mechanism. In Eq.~(\ref{eq:pj_sparse}), we adopt a standard L2-induced projection identical to Eq.~(\ref{eq:proj-l2}), which leads to sparse solutions \cite{wainwright2008graphical}. Intuitively, the projection of a point onto a simplex in terms of L2 distance tends to lie on a facet or a vertex of a simplex, thus being sparse.
\begin{subequations}\label{eq:projs}
    \begin{equation}\label{eq:pj_soft}
    P_\mathcal{S}^{\text{soft}}(\mathbf{A}_{j,:})=\argmin_{\mathbf{B}\in\mathcal{S}}\mathbf{A}_{j,:}^\top \mathbf{B}_{:,j}-\mathbb{H}^{\text{entr}
    }(\mathbf{B}_{:,j}) 
    \end{equation}
    \begin{equation}\label{eq:pj_sparse}
    \begin{split}
        P_\mathcal{S}^{\text{sparse}}(\mathbf{A}_{j,:})&=\argmin_{\mathbf{B}\in\mathcal{S}}\mathbf{A}_{j,:}^\top \mathbf{B}_{:,j}-\mathbb{H}^{\text{gini}}(\mathbf{B}_{:,j}) \\
        &=\argmin_{\mathbf{B}\in\mathcal{S}}|\mathbf{A}_{j,:}-\mathbf{B}_{:,j}|^2
    \end{split}
    \end{equation}
\end{subequations}
where $\mathbb{H}^{\text{entr}}(\cdot)$ and $\mathbb{H}^{\text{gini}}(\mathbf{y})= \frac{1}{2} \sum_i\mathbf{y}_i (\mathbf{y}_i - 1)$ are the standard entropy and the gini-entropy, respectively. $\mathbf{A}_{j,:}$ and $\mathbf{B}_{:,j}$ are the $i$th row and column of $\mathbf{A}$ and $\mathbf{B}$, respectively. 
While the solution to Eq.~(\ref{eq:pj_soft}) is $\mathrm{softmax}(\mathbf{A})$, Eq.~(\ref{eq:pj_sparse}) also has closed-form solution shown in Appendix \ref{app:soft_sparse}.

Comparing Eq.~(\ref{eq:pj_soft}) to the standard Euclidean projection in Eq.~(\ref{eq:pj_sparse}), we note the entropic regularization $\mathbb{H}(\cdot)$ in Eq.~(\ref{eq:pj_soft}) allows for a smoother trajectory by projecting any $\mathbf{A}$ into the interior of $(n-1)$-simplex.
We visualize the two versions of projections in Eq.~(\ref{eq:projs}) onto $1$-simplex from some random points in Fig.~\ref{fig:softVSsparse}. One can readily see that Eq.~(\ref{eq:pj_sparse}) is an exact projection such that points far from the simplex are projected onto the boundary. However, $\mathrm{softmax}$ is smoother by projecting all points onto a relative interior of 1-simplex without sudden change. In the context of Bregman divergence, different distances can facilitate efficient convergence under different ``L-relative smoothness'' \cite{dragomir2021fast}, which can potentially accelerate the learning of dynamics. We leave this to our future work.

We further discuss some neural-friendly features of Eq.~(\ref{eq:pj_soft}) and (\ref{eq:pj_sparse}) facilitating the neural computation:

\textbf{(1)} First, the neural computational graph can be simplified using projection Eq.~(\ref{eq:pj_soft}). Though Eq.~(\ref{eq: att_implement}) using projection Eq.~(\ref{eq:pj_soft}) defines a projected dynamical system directly on $\mathbf{A}$, we switch to update the system using $\mathbf{L}$ as follows, which is considered to further ease the forward integration. This is achieved by instead modeling the dynamics of the feature before fed into $\mathrm{softmax}(\cdot)$:
\begin{subequations}\label{eq: transform}
\begin{small}
    \begin{alignat}{2}
        \mathbf{A}(t) &= \text{Softmax}(\mathbf{L}(t)) \\
        \mathbf{L}(t) &= \mathbf{L}(0) + \int_0^t \frac{\mathrm{d}}{\mathrm{d}s} \frac{\mathbf{Q}(\mathbf{z}(s)) \cdot \mathbf{K}^\top(\mathbf{z}(s))}{\sqrt{d_k}} \mathrm{d} s, \\
        \mathbf{L}(t + \Delta t) &= \mathbf{L}(t) + \Delta t \cdot \frac{\mathrm{d}}{\mathrm{d}s} \frac{\mathbf{Q}(\mathbf{z}(s)) \cdot \mathbf{K}^\top(\mathbf{z}(s))}{\sqrt{d_k}} \bigg|_{s = t}
    \end{alignat}
\end{small}
\end{subequations}
where $\mathbf{Q}(\cdot)$ and $\mathbf{K}(\cdot)$ correspond to the query and key in the attention mechanism, respectively. $\mathbf{L}(0) = \mathbf{Q}(\mathbf{z}(0)) \cdot \mathbf{K}^\top(\mathbf{z}(0)) / \sqrt{d_k}$. We show that updating the dynamic of $\mathbf{L}$ following Eq.~(\ref{eq: transform}) is equivalent to directly updating $\mathbf{A}$ in Appendix \ref{app:equi}.

\textbf{(2)} Second, both the solution to projection Eq.~(\ref{eq:pj_sparse}) and its gradient w.r.t. $\mathbf{A}$ are in closed form. See Proposition \ref{prop:sol_sparse} and Proposition \ref{prop:gradient_sparse} in Appendix \ref{app:soft_sparse} for more details. This, in turn, eases the computational flow in the neural architecture with high efficiency and stability.

Though only two versions of projections are discussed under Bregman divergence, we believe they are sufficiently distinguishable for analyzing the behavior of ProjDEs. For generic neural-friendly projections, we leave them to our future work.

\begin{figure}[tb]
  \centering
    \includegraphics[width=.30\textwidth]{figs/soft_vs_sparse.pdf}
  \caption{Comparsion of softmax and L2 projection onto a simplex. We see that the softmax projection trends to project onto the ``center'' of the simplex while the L2 projection trends to project onto the corner.}
  \label{fig:softVSsparse}
\end{figure}

\textbf{Integration.} We employ the standard Euler's discretization for performing the forward integration by updating $\mathbf{z}$ and $\mathbf{A}$ simultaneously with a sufficiently small time step. We term our approach a \textbf{d}ecoupling-based \textbf{n}eural \textbf{s}ystem (\textbf{DNS}) using projection Eq.~(\ref{eq:pj_soft}) and \textbf{DNS\textsubscript{G}} using projection Eq.~(\ref{eq:pj_sparse}), respectively.

\begin{figure*}[tb]
\captionsetup[subfigure]{labelformat=empty}
    \centering
    \foreach \t in {5,...,12}{
        \begin{subfigure}{.11\textwidth}
            \centering
            \includegraphics[width=.98\textwidth, trim=100 100 100 100,clip]{figs/softmax_att/trajectory_\t.pdf}\\
            \includegraphics[width=.98\textwidth, trim=100 100 100 100,clip]{figs/softmax_att/attention_\t.pdf}
            \caption{$t=$\t}
        \end{subfigure}
    }
  \vspace{-8pt}  
  \caption{A figure showing the corresponding three-body trajectory (on the top), as well as the evolution over time on interactions (at the bottom) between three \textbf{latent sub-systems} in a Three-Body environment. Timestamp from 5 to 12.}
  \label{fig: att}
\end{figure*}
\section{Experiments}
We evaluate the performance of DNS on multiple synthetic and real-world datasets. More details about the dataset and implementation details can be found in Appendix \ref{app: experiment details} and \ref{app: dataset}. Throughout all the tables consisting of the results, ``-'' indicates ``not applicable'' since RIM cannot handle irregular cases. 

\begin{remark}
    In all the experiments, the input feature is treated holistically without any distinguishable parts. For example, in the Three Body dataset, the input is a 9-dimensional vector, with every 3 dimensions (coordinates) from a single object. However, this prior is not fed into any models in comparison. Thus, we do not compare to models integrated with strong prior such as \cite{kipf2018neural}.
\end{remark}

\textbf{Baselines.} We compare DNS with several selected models capturing interactions or modeling irregular time series, including \textbf{CT-GRU} \citep{mozer2017discrete} using state-decay decay mechanisms, \textbf{RIM} \citep{goyal2019recurrent} updating almost independent modules discretely, and \textbf{NeuralCDE} \citep{kidger2020neural} which reports state-of-the-art performance on several benchmarks. 
\begin{table}[t]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Trajectory prediction}. MSE loss of the three body dataset ($\times 10^{-2}$).}
\label{tab: three_body}
\begin{tabular}{lccr}
\toprule
    Model  & Regular & Irregular \\
\midrule
CT-GRU    & 1.8272   & 2.4811 \\
NeuralCDE & 3.3297   & 5.0077 \\
RIM       & 2.4510   & -      \\
\midrule
DNS       & \textbf{1.7573}  & \textbf{2.2164} \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\textbf{Adapting DNS to the noisy case.}  
To allow DNS fitting to noisy and uncertain circumstances, we create a variant by slightly modifying it.
This variant is obtained by replacing cubic spline interpolation over $\mathbf{x}_i(t)$ with natural smoothing spline \citep{green1993nonparametric}, in consideration of incorporating smoother controls and alleviating data noise. This version is termed as \textbf{DNS\textsubscript{S}}. 

\subsection{Three body}
The three-body problem is characterized by a chaotic dynamical system for most randomly initial conditions. A small perturbation may cause drastic changes in the movement. Taking into account the problem's complexity, it is particularly suitable for testing our approach. 
In this experiment, we consider a trajectory predicting problem given the noisy historical motion of three masses, where gravity causes interactions between them. 
Therefore, models need to (implicitly) learn both Newton's laws of motion for modeling sub-system dynamics and Newton's law of universal gravitation to decouple the latent interaction. This dataset consists of 50k training samples and 5k test samples. For each sample, 
8 historical locations for the regular setting and 6 historical locations (randomly sampled from 8) for the irregular setting in the 3-dimensional space of three bodies are given to predict 3 subsequent locations. To equip with the cluttered setting, the correspondence between dimensions and bodies will not be fed into the learning models, hence a 9-dimensional observation at each time stamp. Models' performance is summarized in Table~\ref{tab: three_body}. 
\begin{figure}[t]
    \centering
    \includegraphics[width=.4\textwidth, trim = {100 40 100 0, clip}]{figs/three-body-focus.pdf}
    \caption{A figure showing the focus of 3 sub-systems on 9-dimensional input of Three Body. The strength of focus is reflected by the thickness of the lines.}
    \label{fig:body_focus}
\end{figure}
We can conclude that DNS outperformed all the selected counterparts in both regular and irregular settings. Notably, although our method is built on NeuralCDE, with the decoupling, the performance can be significantly improved. See Table~\ref{tab: appendix three body} in Appendix \ref{app:three body} for more detailed results.
\begin{table*}[t]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Link prediction}. Accuracy on Spring ($\%$). \textsc{Clean}, \textsc{Noisy}, and \textsc{Short} correspond to settings with clean, noisy, and short portion data, respectively. Detailed results for \textsc{Clean} and \textsc{Noisy} are separately summarized in Tab.~\ref{tab:spring app} and Tab.~\ref{tab: appendix noisy spring} in the appendix.}
\label{tab: spring}
\begin{tabular}{lcc|cc|cc}
\toprule
    \multirow{2}{*}{Model} & \multicolumn{2}{c}{Clean} & \multicolumn{2}{c}{Noisy} & \multicolumn{2}{c}{Short}\\
    \cline{2-7}
    & Regular & Irregular & Train\&Test & Test & 50\% & 25\% \\
    \midrule
    CT-GRU & 92.89$\pm$0.52   &  88.47$\pm$0.34 & 92.71$\pm$0.55 & 92.80$\pm$0.53 & -& -\\
    NeuralCDE & 92.47$\pm$0.06   &  89.74$\pm$0.18 & 90.76$\pm$0.08 & 89.61$\pm$0.09 & -& -\\
    RIM & 89.73$\pm$0.07   &  -  &89.65$\pm$0.14 & 89.64$\pm$0.10  & -& - \\
    \midrule
    DNS\textsubscript{G} & 94.31$\pm$0.48 & \textbf{94.25$\pm$0.29} & - & - & 92.58 & 92.31 \\
    DNS\textsubscript{S} & - & - & \textbf{93.67}$\pm$\textbf{0.57} & \textbf{92.99}$\pm$\textbf{1.30} & -& -\\ 
    DNS & \textbf{94.44$\pm$0.69}  & 93.60$\pm$1.21 & 93.42$\pm$1.05 & 89.56$\pm$0.42 & 91.11 & 92.13\\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table*}

\textbf{Visualization and analysis.} We visualize dynamics $\mathbf{A}$ of DNS along the movements of three body system. See Fig. \ref{fig: att} for results. We set the time stamps starting from 5 to 12 to make visualization more informative. It is seen in the beginning ($t=5,6$ or even earlier), $\mathbf{A}$ remains stable as the three bodies are apart from each other without intensive interactions. At $t=7$, $\mathbf{A}$ demonstrates obvious change when two bodies start to the coil. Another body joins in this party at $t=8$, yielding another moderate change of $\mathbf{A}$. When flying apart, one body seems more independent, while another two keep entangled together. These are well reflected via the meta-system $\mathbf{A}$. To further see how the holistic 9-dimensional input is decoupled into sub-systems $\mathbf{z}_i$, we visualize the sub-system focus in Fig.~\ref{fig:body_focus} (also see Appendix \ref{app:cam}). Interestingly, latent entities (sub-systems) do not correspond to physical entities (three bodies). Instead, the first sub-system puts more focus on the whole input, but the remaining two sub-systems concentrate on the x-axis and y/z-axis, respectively. Though counterintuitive, this unexpected decoupling exhibits good performance. We will investigate how to decouple out physical entities from cluttered observations in our future work.
\begin{table}[tb]
\centering
\begin{small}
\begin{sc}
\caption{\textbf{Link prediction}. Ablation study. ($\%$). 
}
\label{tab:abl}
\begin{tabular}{lc}
\toprule
Control & Accuracy ($\%$) \\
\midrule
No encoding       & 91.57                \\ 
MLP(2$\times$input)     & 91.51                \\ 
MLP(16$\times$input)      & 91.17                \\  
DNS (8$\times$MLP(2$\times$input))      & \textbf{95.38}                \\ 
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{table}

\begin{figure*}[tb]
\captionsetup[subfigure]{labelformat=empty}
\centering
    \foreach \t in {1,2,3,4,8,13,19,28}{
        \begin{subfigure}{.11\textwidth}
        \centering
        \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/traj_\t.pdf}\\ \vspace{-5pt}
         \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/soft_\t.pdf}\\ \vspace{-5pt}
        \includegraphics[width=.98\textwidth, trim={25 25 25 0}]{figs/spring_att/sparse_\t.pdf}
        \vspace{-15pt}
        \caption{$t=$\t}
        \end{subfigure}
    }
    \vspace{-10pt}
\caption{Visualization of the evolution of the meta-systems of DNS and DNS\textsubscript{G} on Spring dataset. On each time stamp $t$, from top to bottom, we show the trajectory of the 5 balls, the meta-system state of DNS, and the meta-system state of DNS\textsubscript{G}, respectively.}
\label{fig:vis_spring}
\end{figure*}



\subsection{Spring}
We experiment with the capability of DNS in decoupling the independence in complex dynamics controlled by simple physics rules. We use a simulated system in which particles are connected by (invisible) springs \citep{kuramoto1975self,kipf2018neural}. Each pair of particles has an equal probability of having an interaction or not. Our task is to use observed trajectory to predict whether there are springs between any pair of two particles, which is analogous to the task of link prediction under a dynamical setting. This can be inferred from whether two trajectories change coherently. The spring dataset consists of 50k training examples and 10k test examples. Each sample has a length of 49. We test a variety of combinations of the number of sub-systems and dimensions of the hidden state. Experimental results are in Table~\ref{tab: spring}. To test the models' noise resistance, we add Gaussian noise to the spring dataset and obtain the noisy spring dataset. We set two scenarios, ``Train\&Test'' and ``Test'', corresponding to injecting noise at both training and test phases and only at testing phases, respectively. Experimental results are in Table~\ref{tab: spring}.

\textbf{Clean Spring.} From \textsc{Clean} part of Table \ref{tab: spring}, we see variants of DNS stably outperform all the selected counterparts by a large margin. Especially, under the irregularly sampled data, DNS and DNS\textsubscript{G} have a remarkable performance gap with all other methods and maintain reliability as in the regular setting. We believe this is significant since learning from irregularly sampled data is typically much more difficult than learning from normal data.

\textbf{Noisy Spring.} According to \textsc{Noisy} part of Table \ref{tab: spring}, DNS\textsubscript{S} is quite reliable in noisy cases. It seems a smoothing procedure on the controls can be helpful under massive uncertainty. Also, we see that adding noise tends to damage the performance of all methods. This also raises one of our future research directions to investigate how to handle different controls. Without applying a smooth cubic spline, DNS can still have a good performance, which indicates that by decoupling, the model focuses on learning latent interaction patterns, and patterns are less susceptible to noise. 

\textbf{Visualization and analysis.}  We also visualize state $\mathbf{A}$ of meta-systems over time in Fig.~\ref{fig:vis_spring} for Spring. From top to bottom, the first, second and third rows correspond to the trajectory of particles, meta-system state of DNS, and meta-system state of DNS\textsubscript{G}. One interesting thing we note is that the interactions in DNS\textsubscript{G} almost concentrate on the starting portion of all the time stamps. At $t=8$ and after, there is no interaction at all. Though not obvious, this also happens to DNS in the sense that $\mathbf{A}$ tends to be diagonal. We suppose this is because DNS and DNS\textsubscript{G} only need a portion of data from the start to determine the existence of a link rather than looking into all the redundant time stamps. 

\textbf{Short Spring.} We thus verify this by training and testing both variants with 50\% and 25\% of data cropped from the starting time stamp and summarize results in \textsc{Short} part of Table \ref{tab: spring}. It is seen that incomplete data in this task only slightly impact the performance. And this can be surprisingly reflected in the evolution of meta-systems. This also aligns with the intuition that \emph{Link prediction} needs fewer data than \emph{Trajectory prediction} as in Three Body.
\begin{table}[tb]
\begin{center}
\begin{small}
\begin{sc}
\caption{\textbf{Video classification}. Accuracy of the human actions dataset ($\%$). \textsc{Norm} and \textsc{Unnorm} refer to normalized and unnormalized inputs, respectively. Detailed results with superscript $^\dagger$ and $^\ddagger$ are in Tab.~\ref{tab:action app} and Tab.~\ref{tab:v_unnorm}, respectively.}
\label{tab: video}
\begin{tabular}{lccc}
\toprule
    \multirow{2}{*}{Model} &  \multicolumn{1}{c}{Norm} & \multicolumn{2}{c}{Unnorm}\\
    \cline{2-4}
     & Irreg & Reg & Irreg \\
\midrule
CT-GRU    & 67.30$\pm$6.19$^\dagger$   & 60.33$^\ddagger$ & 66.67$^\ddagger$\\
NeuralCDE & 89.73$\pm$3.38$^\dagger$   & 70.33$^\ddagger$ & 59.17$^\ddagger$\\
RIM       & -      & 55.50$^\ddagger$ & - \\
\midrule
DNS       & \textbf{91.35}$\pm$\textbf{3.48}$^\dagger$  & \textbf{97.00}$^\ddagger$ & \textbf{95.33}$^\ddagger$\\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table}

\textbf{Ablation study.} Since our method merely incorporates an extra meta-system and a control encoder for modeling the interaction compared to standard NeuralCDE, we conduct experiments under different settings to see how different encoders and hidden state dimensions can contribute to improving NeuralCDE. To ensure fairness, we cast a 2-layer MLP with different output sizes (2 and 16 times of input size) as in DNS to obtain varying sizes of controls. Results are summarized in Table \ref{tab:abl} (detailed in Tab.~\ref{tab: appendix spring abl}). We see that with an extra control encoder, there is no obvious performance difference among these settings. However, once the interaction meta-system is imposed, DNS can achieve quite significant performance gain. This, in turn, shows the necessity of the proposed meta-system for explicitly modeling the evolving interactions.

\subsection{Human actions}
The recognition of human actions dataset contains three types of human actions, which are hand clapping, hand waving, and jogging \citep{schuldt2004recognizing}. For this dataset, we consider the limbs of the character as subsystems. When the character does one kind of action, subsystems interact in a specific pattern. We test the performance of all the selected models with the learnable backbone Resnet18 \citep{he2016deep}. We also test the compatibility of all methods with different dynamical ranges: \textsc{Norm} and \textsc{Unnorm} indicate pixel value in $[0,1]$ and $[0,255]$, respectively. Experimental results are summarized in Table \ref{tab: video}. DNS consistently outperforms all other methods, and exhibits strong compatibility to drastically changed ranges under \textsc{Unnorm} setting. Thus it is potentially more flexible to be integrated into various tasks with a large dynamical range (e.g., earthquake).

To view how the decoupling works for video recognition tasks, we visualize the strength of the learned parameters by mapping the 128-D feature into 6 latent sub-systems in Figure \ref{fig: cnn top feature} with re-ordered indices for better view. 
It can be seen that there are some obvious latent structures in the grouping of the parameters 128-D control to the system. Each sub-system mainly focuses on a small portion of the control, based on which we can infer that each sub-system models different components in inputted images.

\begin{figure}[tb]
  \centering
    \includegraphics[width=.45\textwidth, trim=0 20 0 20, clip]{figs/top40feature.pdf}
    \vspace{-10pt}
  \caption{A figure showing the importance of each feature vector entry for subsystems}
  \label{fig: cnn top feature}
\end{figure}

\section{Conclusion}
In this paper, we propose a method for modeling cluttered and irregularly sampled sequential data. Our method is built upon the assumption that complex observation may be derived from relatively simple and independent latent sub-systems, wherein the interactions also evolve over time. We devise a strategy to explicitly decouple such latent sub-systems and a meta-system governing the interaction. Inspired by recent findings of projected differential equations and the tool of Bregman divergence, we present a novel interpretation of our model and pose some potential future directions. Experiments on various tasks demonstrate the prominent performance of our method over previous state-of-the-art methods.


