\section{Introduction}

Many efforts in structural biology aim to predict, or derive insights from, the structure of a macromolecule (such as a protein, RNA, or DNA), represented as a set of positions associated with atoms or groups of atoms in 3D Euclidean space. These problems can often be framed as functions mapping the input domain of structures to some property of interest---for example, predicting the quality of a structural model or determining whether two molecules will bind. Thanks to their importance and difficulty, such problems, which we broadly refer to as \textit{learning from structure}, have recently developed into an exciting and promising application area for deep learning \citep{gainza2020deciphering, graves2020review, ingraham2019generative, pereira2016boosting, townshend2019end, won2019assessment}.

Successful applications of deep learning are often driven by techniques that leverage the problem structure of the domain---for example, convolutions in computer vision \citep{cohen2016inductive} and attention in natural language processing \cite{vaswani2017attention}. What are the relevant considerations in the domain of learning from structure?
Using proteins as the most common example, the spatial \textit{arrangement} and \textit{orientation} of the amino acids govern the topology, dynamics and function of the molecule \citep{berg2002biochemistry}. These properties are in turn mediated by key pairwise and higher-order residue-residue \textit{interactions} \citep{del2006residue, hammes2006relating}. We refer to these as the \textit{geometric} and \textit{relational} aspects of the problem domain, respectively.

Current state-of-the-art methods for learning from structure are successful by leveraging one of these two aspects. Commonly, such methods either employ graph neural networks (GNNs), which are expressive in terms of relational reasoning \citep{battaglia2018relational} or convolutional neural networks (CNNs), which operate directly on spatial features. Here, we present a unifying architecture that bridges spatial and graph-based methods to leverage \textit{both} aspects of the problem domain.

We do so by introducing \textit{geometric vector perceptrons} (GVPs), a drop-in replacement for standard multi-layer perceptrons (MLPs) in aggregation and feed-forward layers of GNNs. GVPs operate directly on both scalar and \textit{geometric} features---features that transform as a vector under rotation of spatial coordinates. GVPs therefore allow for the embedding of geometric information at nodes and edges without reducing such information to hand-picked scalars that may not fully capture complex geometric relationships. We postulate that our approach makes it easier for a network to learn functions whose significant features are both geometric and relational.

Our method (GVP-GNN) can be applied to any problem where the input domain is a structure of a single macromolecule or of molecules bound to one another. In this work, we demonstrate our approach on two important problems connected to protein structure: model quality assessment and computational protein design. Model quality assessment (MQA) aims to select the best structural model of a protein from a large pool of candidate structures and is a crucial step in protein structure prediction \citep{cheng2019estimation}. Computational protein design (CPD) is the conceptual inverse of structure prediction, aiming to infer an amino acid sequence that will fold into a given structure. Our method outperforms existing methods on both tasks.

\section{Related Work}
Current methods for learning from protein structure generally use one of three classes of structural representations, which we outline below. We also describe representative and state-of-the-art examples for MQA and CPD to set the stage for our experiments later.

\paragraph{Sequential representation}
One way of representing a protein structure is as a sequence of feature vectors, one for each amino acid, that can be input into a 1D convolutional network. These representations intrinsically encode neither geometric nor relational information and instead rely on hand-selected
features to represent the 3D structural neighborhood of each amino acid. Such features can include contact-based features \citep{olechnovivc2017voromqa}, orientations or positions collectively projected to local coordinates \cite{karasikov2019smooth, wang2018computational}, and physics-inspired energy terms \citep{o2018spin2, uziela2017proq3d}. These methods were among the earliest developed for learning from structure, and a number of them (ProQ3D \citep{uziela2017proq3d}, VoroMQA \citep{olechnovivc2017voromqa}, SBROD \citep{karasikov2019smooth}) remain competitive in assessments of MQA.

\paragraph{Graph-based representations}
Graph-based methods represent amino acids or individual atoms as nodes and draw edges based on chemical bonds or spatial proximity, and are well-suited for complex relational reasoning \citep{battaglia2018relational}. These methods reduce the challenging task of representing the collective structural neighborhood of an amino acid to that of representing individual edges. The CPD method ProteinSolver \citep{strokach2020fast}, for example, uses pairwise Euclidean and sequence distance to represent edges. However, to encode more complex geometric information, a choice of scalar representation scheme is again necessary
\textit{i.e.}, vector-valued features whose components are intrinsically linked need to be decomposed into individual scalars.
For example, the MQA method ProteinGCN \citep{sanyal2020proteingcn} projects pairwise relative positions into local coordinates. The CPD method Structured Transformer \citep{ingraham2019generative} additionally uses a quaternion representation of pairwise relative orientations.

\paragraph{Voxel-based representations}
Many structural biology methods voxelize proteins or neighborhoods of amino acids into occupancy maps and then perform hierarchical 3D convolutions. These architectures circumvent the need for hand-picked representations of geometric information altogether and instead directly operate on approximate positions in 3D space. However, CNNs do not directly encode relational information within the structure. In addition, standard CNNs do not intrinsically account for the symmetries of the problem space. To address the latter, early methods, such as the MQA method 3DCNN \citep{derevyanko2018deep}, learned rotation invariance via data augmentation. More recent MQA methods such as Ornate \citep{pages2019protein} and many end-to-end CPD models \citep{anand2020protein, zhang2019prodconn} instead define voxels based on a local coordinate system.

Voxel-based and graph-based architectures, catering to geometric and relational aspects of the problem domain, respectively, have complementary strengths and weaknesses. To address this, some recent work on small molecules aims to leverage geometric features in GNNs through voxelized convolutions at each graph node \citep{spurek2019geometric}, message transformation based on encoded angle information \citep{klicpera2020directional} or products with position vectors to account for vector features \citep{cho2019enhanced}. Our contribution is similar in spirit but more general and expressive in methodology. With the introduction of GVPs, we provide a conceptually simple augmentation of GNNs that allows for arbitrary scalar and vector features at all steps of graph propagation, and has the desired equivariance properties with respect to rotations of atomic structures in space.

\section{Methods}
In this section, we formally introduce the geometric vector perceptron, discuss datasets and evaluation metrics for quality assessment and protein design, and describe the protein representations and model architectures used for those tasks.

\subsection{Geometric vector perceptrons}
Given a tuple $(\mathbf{s}, \mathbf{V})$ of scalar features $\mathbf{s} \in \mathbb{R}^n$ and vector features $\mathbf{V}\in\mathbb{R}^{\nu\times 3}$,
the geometric vector perceptron computes new features $(\mathbf{s}',\mathbf{V}') \in \mathbb{R}^m \times \mathbb{R}^{\mu \times 3}$. As outlined in pseudocode in Algorithm~\ref{alg:gvp} and illustrated in Figure~\ref{fig:gvp}A, the computation is analogous to the linear combination of inputs in normal dense layers, with the addition of a norm operation and an intermediary hidden layer, which allow information in the vector channels to be extracted as scalars.

\begin{algorithm}
\caption{Geometric vector perceptron}\label{alg:gvp}
\begin{algorithmic}[H]
\State \textbf{Input}: Scalar and vector features $(\textbf{s},\textbf{V})$ $\in$ $\mathbb{R}^n \times \mathbb{R}^{\nu \times 3}$ .
\State \textbf{Output}: Scalar and vector features $(\textbf{s}',\textbf{V}')$ $\in$ $\mathbb{R}^m \times \mathbb{R}^{\mu \times 3}$ .
\State $h \gets \max \left(\nu, \mu\right)$

\BState \textbf{GVP}:
\State \ \ \ \ $\text{V}_h \gets \textbf{W}_h \textbf{V} \quad \in \mathbb{R}^{h \times 3} $
\State \ \ \ \ $\text{V}_\mu \gets \textbf{W}_\mu \text{V}_h \quad \in \mathbb{R}^{\mu \times 3} $
\State \ \ \ \ $\text{s}_h \gets \norm{\text{V}_h} \ (\text{row-wise}) \quad \in \mathbb{R}^{h}$
\State \ \ \ \ $\text{v}_\mu \gets \norm{\text{V}_\mu} \ (\text{row-wise}) \quad \in \mathbb{R}^{\mu}$
\State \ \ \ \ $\text{s}_{h+n} \gets \text{concat}\left(\text{s}_h, \textbf{s}\right) \quad \in \mathbb{R}^{h+n}$
\State \ \ \ \ $\text{s}_m \gets \textbf{W}_m \text{s}_{h+n} + \textbf{b}  \quad \in \mathbb{R}^{m}$
\State \ \ \ \ $\textbf{s}' \gets \sigma\left(\text{s}_m \right) \quad \in \mathbb{R}^{m}$
\State \ \ \ \ $\textbf{V}' \gets  \sigma^+\left(\text{v}_\mu\right) \odot \text{V}_\mu, (\text{row-wise multiplication}) \quad \in \mathbb{R}^{\mu \times 3}$
\BState  \textbf{return} $(\textbf{s}',\textbf{V}')$
\end{algorithmic}
\end{algorithm}

$\mathbf{W}_h \in \mathbb{R}^{h \times \nu}$, $\mathbf{W}_\mu \in \mathbb{R}^{\mu \times h}$, and $\mathbf{W}_m \in \mathbb{R}^{m \times (n+h)}$ are weight matrices and $\mathbf{b}\in\mathbb{R}^m$ is a bias vector.
$\mathbb{\sigma}$ and $\mathbb{\sigma^+}$ denote activation functions. In our experiments, we choose ReLU and sigmoid for $\mathbb{\sigma}$ and $\mathbb{\sigma^+}$, respectively, although other nonlinearities may be used.
Additionally, $h$ may also be defined independently, but we use the maximum of $\nu, \mu$ for convenience.

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{schematic_v6.png}
    \caption{
     \textbf{(A)} Schematic of the geometric vector perceptron illustrating the computation shown in Algorithm~\ref{alg:gvp}. Given a tuple of scalar and vector input features $(\mathbf{s},\mathbf{V})$, the perceptron computes an updated tuple $(\mathbf{s}',\mathbf{V}')$. $\mathbf{s}'$ is a function of both $\mathbf{s}$ and $\mathbf{V}$.
    \textbf{(B)} Illustration of the structure-based prediction tasks. In model quality assessment (top), the goal is to predict a quality score given the 3D structure of a candidate model. Individual atoms are represented as colored spheres. The quality score measures the accuracy of a candidate structure with respect to an experimentally determined structure (shown in gray). In computational protein design (bottom), the goal is to predict an amino acid sequence that would fold into a given protein backbone structure.
   }
    \label{fig:gvp}
\end{figure}

Despite its conceptual simplicity, the GVP offers some desirable properties. First, the vector and scalar outputs of the GVP are equivariant and invariant, respectively, with respect to an arbitrary composition of rotations and reflections in 3D Euclidean space described by $R$ \textit{i.e.},

\begin{equation}
\label{eq:eq}
\text{GVP}(\left(\mathbf{s}, R(\mathbf{V}))\right) = (\textbf{s'}, R(\textbf{V}'))
\end{equation}

This is due to the fact that the only operations on vector-valued inputs are scalar multiplication, linear combination, and the $L_2$ norm. We include a formal proof in Supplementary Information.

In addition, a GVP $G_s$ defined with $n, \mu=0$---that is, the part of a GVP that transforms vector features to scalar features---is able to $\epsilon$-approximate a function $F: \mathbb{R}^{\nu \times 3} \rightarrow \mathbb{R}$ that is invariant with respect to rotations and reflections in 3D under mild assumptions.

\newtheorem*{theorem}{Theorem}
\label{th:eps}
\begin{theorem}
    Let R describe an arbitrary rotation or reflection in $\mathbb{R}^3$.
    For $\nu \ge 3$ let $\Omega^\nu \subset \mathbb{R}^{\nu \times 3}$ be the set of all $\mathbf{V} = \begin{bmatrix}\mathbf{v}_1, & \ldots, & \mathbf{v}_\nu\end{bmatrix}^T\in \mathbb{R}^{\nu \times 3}$ such that $\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3$ are linearly independent and $0<||\mathbf{v}_i||_2 \le b$ for all $i$ and some finite $b>0$.

    Then for any continuous $F: \Omega^\nu \rightarrow \mathbb{R}$ such that $F(R(\mathbf{V})) = F(\mathbf{V})$ and for any $\epsilon>0$, there exists a form $f(\mathbf{V})=\mathbf{w}^TG_s(\mathbf{V})$ such that $|F(\mathbf{V})-f(\mathbf{V})|< \epsilon$ for all $\mathbf{V}\in\Omega^\nu$.
\end{theorem}
\begin{proof} In Supplementary Information. \end{proof}

As a corollary, a GVP with nonzero $n,\mu$ is also able to approximate similarly-defined functions over the full input domain $\mathbb{R}^n \times \mathbb{R}^{\nu \times 3}$.

Finally, the GVP layer is nearly as fast as normal dense layers, incurring additional overhead only in reshaping and concatenation operations.

In addition to the GVP layer itself, we use a version of dropout that drops entire vector channels at random (as opposed to components within vector channels). We also introduce layer normalization for the vector features as
\begin{equation}
\mathbf{v}_i \leftarrow \mathbf{v}_i / \sqrt{\frac{1}{|\{j\}|}\sum_j||\mathbf{v}_j||_2^2} \quad \forall \ i \quad ,
\end{equation}
where $\mathbf{v}_i \in \mathbb{R}^3$ are the individual row vectors of the vector feature matrix $\mathbb{R}^{\nu \times 3}$.
This vector layer norm has no trainable parameters, but we continue to use normal layer normalization on scalar channels with trainable parameters $\gamma, \beta$.

\subsection{Dataset and evaluation}
We benchmark GVP-augmented graph neural networks on two distinct challenges in structural biology: model quality assessment and protein design. Figure~\ref{fig:gvp}B illustrates these tasks.
\subsubsection{Model quality assessment}
Model quality assessment aims to select the best structural model of a protein from a large pool of candidate structures.\footnote{We refer to models of protein structure as "structural models" or "candidate structures" to avoid confusion with the term "model" as used in the ML community.} The performance of different MQA methods is regularly assessed in the biennial blind prediction competition CASP, of which the 13th was the most recent \citep{cheng2019estimation}. For a number of recently solved but unreleased structures, called \textit{targets}, structure generation programs produce a large number of candidate structures. MQA methods are ranked by how well they are able to predict the GDT-TS score of a candidate structure. The GDT-TS is a  measure of how similar two protein backbones are after global alignment; specifically, it is the mean of
\begin{equation}
    \text{GDT-N} = \frac{1}{\text{\# residues}}\sum_\text{all residues $i$} \mathbbm{1}\left[d(\text{target}_i, \text{candidate}_i)<N \text{ Angstroms}\right]
\end{equation}
for $N=1,2,4,8$, where $d$ is the Euclidean distance. In addition to accurately predicting the absolute quality of a candidate structure, a good MQA method should also be able to accurately assess the \textit{relative} model qualities among a pool of candidates for a given target so that the best ones can be selected, perhaps for further refinement. Therefore, MQA methods are commonly evaluated on two metrics: a global correlation between the predicted and ground truth scores, pooled across all targets, and the average \textit{per-target} correlation among only the candidate structures for a specific target \citep{cao2016protein, derevyanko2018deep, pages2019protein}.

To fulfill this aim, we train our MQA model on an absolute loss and a pairwise loss. That is, for each training step we intake pairs $i,j$ where $i,j$ are candidate structures for the same target and compute
\begin{equation}
    \mathscr{L} = H(y^{(i)}-\hat{y}^{(i)}) + H(y^{(j)}-\hat{y}^{(j)}) +  H\left((y^{(i)}-y^{(j)})-(\hat{y}^{(i)}-\hat{y}^{(j)})\right)
\end{equation}
where $H$ is the Huber loss. When reshuffling at the beginning of each epoch, we also randomly pair up the candidate structures for each target. Interestingly, adding the pairwise term also improves global correlation, likely because the much larger number of possible pairs makes it more difficult to overfit.

We train and validate on the set of candidate structures generated in the CASP 5-10 assessments, which collectively contain 528 targets and 79200 candidate structures. For testing, we predict model quality for a total of 20880 stage 1 and stage 2 candidate structures from CASP 11 (84 targets) and 12 (40 targets).

\subsubsection{Protein design}
\label{sec:dataset_cpd}
Computational protein design is the conceptual inverse of structure prediction, aiming to infer an amino acid sequence that will fold into a given structure. In comparison to model quality assessment, computational protein design is more difficult to unambiguously benchmark, as some structures may correspond to a large space of sequences and others may correspond to none at all.
Therefore, the proxy metric of \textit{native sequence recovery}---splitting the set of all known native structures in the PDB and attempting to design the sequences corresponding to held-out structures---is typically used to benchmark CPD models \citep{li2014direct, o2018spin2, wang2018computational}. Drawing an analogy between sequence design and language modelling, \citet{ingraham2019generative} also evaluate the model \textit{perplexity} on held-out native sequences. Both metrics rest on the implicit assumption that native sequences are optimized for their structures \citep{kuhlman2000native} and should be assigned high probabilities.

To best approximate real-world applications that may require design of novel structures, the held-out evaluation set should bear minimal similarity to the training structures. We use the CATH 4.2 dataset curated by \citet{ingraham2019generative} in which all available structures with 40\% nonredudancy are partitioned by their CATH (class, architecture, topology/fold, homologous superfamily) classification. The canonical training, validation, and test splits consist of 18204, 608, and 1120 structures, respectively.

\subsection{Model architecture}

\paragraph{GVP-GNNs}
In a computation graph, a GVP can be placed at any node typically occupied by a dense layer. To most directly augment relational representations with geometric information, we modify GNNs by permitting all edge and node embeddings to have geometric vector and scalar features, and then use GVPs at all graph propogation and point-wise feed-forward steps.

\paragraph{Representations of proteins}
A protein structure is a sequence of amino acids, where each amino acid consists of four \textit{backbone} atoms\footnote{C$\alpha$, C, N, and O.}
and a set of \textit{sidechain} atoms located in 3D Euclidean space. Here we represent only the backbone  because our MQA benchmark corresponds to the assessment of backbone structure. In CPD, the sidechains are by definition unknown.

Let X$_i$ be the position of atom X in the $i$th amino acid (where X is a C$\alpha$, C, N, or O atom).  We represent backbone structure as a graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ where each node $\mathfrak{v}_i \in \mathcal{V}$ corresponds to an amino acid and has embedding $\mathbf{h}^{(i)}_\mathfrak{v}$ with the following features:
\begin{itemize}
    \item Scalar features $\{\sin, \cos\} \times {\phi, \psi, \omega}$, where $\phi, \psi, \omega$ are the dihedral angles computed from C$_{i-1}$, N$_{i}$, C$\alpha_{i}$, C$_{i}$, and N$_{i+1}$.
    \item The \textit{forward} and \textit{reverse} unit vectors in the directions of C$\alpha_{i+1} - \text{C}\alpha_i$ and C$\alpha_{i-1} - \text{C}\alpha_i$, respectively.
    \item The unit vector in the imputed direction of C$\beta_i-\text{C}\alpha_i$.\footnote{C$\beta$ is the second carbon from the carboxyl carbon C.} This is computed by assuming tetrahedral geometry and normalizing
    \begin{equation*}
        \sqrt{\frac{1}{3}}(\mathbf{n} \times \mathbf{c})/||\mathbf{n} \times \mathbf{c}||_2-\sqrt{\frac{2}{3}}(\mathbf{n} + \mathbf{c})/||\mathbf{n}+\mathbf{c}||_2
    \end{equation*}
    where $\mathbf{n}=\text{N}_i-\text{C}\alpha_i$ and $\mathbf{c}=\text{C}_i-\text{C}\alpha_i$. The three unit vectors unambiguously define the orientation of each amino acid residue.

    \item A one-hot representation of amino acid identity, when available.
\end{itemize}
The set of edges is $\mathcal{E}=\{\mathbf{e}_{i\rightarrow j}\}_{i\neq j}$ for all $i, j$ where $\mathfrak{v}_i$ is among the $k=30$ nearest neighbors of $\mathfrak{v}_j$ as measured by the distance between their C$\alpha$ atoms. This is motivated by the fact that if two entities have a relational dependency they are likely to be close together in space.
Each edge has an embedding $\mathbf{h}^{(i\rightarrow j)}_e$ with the following features:
\begin{itemize}
    \item The unit vector in the direction of C$\alpha_i-\text{C}\alpha_j$.
    \item The radial basis encoding of the distance $\text{RBF}(||\text{C}\alpha_i-\text{C}\alpha_j||_2)$.
    \item A sinusoidal encoding of $i-j$ as described in \citet{vaswani2017attention}, representing distance along the backbone.
\end{itemize}
In our notation, each feature vector $\mathbf{h}$ is a concatenation of scalar and vector features as described above. Collectively, these features are sufficient for a complete description of the protein backbone. In particular, whereas previous graph-based representations depended on scalar edge embeddings to represent relative orientations and positions, we are able to directly encode the absolute orientations of each amino acid and relative positions in the equivalent of fewer real-valued channels.

\paragraph{Network architecture}
Our GVP-GNN takes as input the protein graph defined above and performs graph propagation steps that update the node embeddings according to:
\begin{eqnarray}
    \mathbf{h}^{(i\rightarrow j)}_m &:=& g\left(\text{concat}\left(\mathbf{h}_\mathfrak{v}^{(i)}, \mathbf{h}^{(i\rightarrow j)}_e\right)\right) \\
    \mathbf{h}^{(j)}_\mathfrak{v} &\leftarrow& \text{LayerNorm}\left(\mathbf{h}^{(j)}_\mathfrak{v} + \frac{1}{||\{i:\mathbf{e}_{i\rightarrow j}\in E\}||}\text{Dropout}\left(\sum_{i:\mathbf{e}_{i\rightarrow j}\in E} \mathbf{h}_m^{(i\rightarrow j)}\right)\right)
\end{eqnarray}
where $g$ is a composition of GVPs. We do not update edge embeddings and do not use a global graph embedding. Between node update steps, we also use a residual feed-forward point-wise update layer:
\begin{equation}
    \mathbf{h}^{(i)}_\mathfrak{v} \leftarrow \text{LayerNorm}\left(\mathbf{h}^{(i)}_\mathfrak{v} + \text{Dropout}\left(g\left(\mathbf{h}_\mathfrak{v}^{(i)}\right)\right)\right)
\end{equation}
In model quality assessment, the network performs regression against the true quality score of a candidate structure, a global scalar property. To obtain a single global representation, we apply a node-wise GVP to reduce all features to scalars after all graph propagation steps, and then average the representations across all nodes.\footnote{We also tried learning a weighted average of nodes, but this led to increased overfitting.} A final dense network with dropout and layer normalization then outputs the network's prediction.

In computational protein design, the network learns a generative model over the space of protein sequences conditioned on the given backbone structure. Following \citet{ingraham2019generative}, we frame this as an \textit{autoregressive} task and use a masked encoder-decoder architecture to capture the joint distribution over all positions: for each position $i$, the network models the distribution over amino acids at $i$ based on the complete structure graph, as well as the sequence information at positions $j < i$. The encoder first performs graph propagation on the structural information only. Then, sequence information is added to the graph, and the decoder performs further graph propagation where incoming messages $\mathbf{h}_m^{(i\rightarrow j)}$ for $i\ge j$ are computed only with the encoder embeddings. Finally, we use one last GVP with 20-way scalar output and softmax activation to output the probability of the amino acids.

The MQA network, structure encoder, and masked decoder each use 3 graph propagation layers. The feed-forward modules consist of $2$ GVPs and the message-gather function uses 3.

\section{Experiments}

\subsection{Model quality assessment}
We compare our method against state-of-the-art methods on the CASP 11-12 test set in Table~\ref{tab:casp11-12}.

 These include representatives of voxel-based methods (3DCNN and Ornate), a graph-based method (GraphQA), and three methods that use sequential representations. All of these methods learn solely from protein structure,\footnote{There are two versions of GraphQA; we compare against the one using only structure information.} with the exception of ProQ3D, which in addition uses sequence information on related proteins that is not always available. We include ProQ3D because it is an improved version of the best single-model method in CASP 11 and CASP 12 \citep{uziela2017proq3d}. Our method outperforms all other structural methods in both global and per-target correlation, and even performs better than ProQ3D on all but one benchmark.
\begin{table}
  \caption{Comparison with state-of-the-art methods on CASP 11 and 12 in terms of global (Glob) and mean per-target (Per) Pearson correlation coefficients (higher is better). ProQ3D is set aside as the sole non-structure-only method. The top performing structure-only method for each metric is in bold, as is the top performing-method overall (if different). Our method generally improves over all other methods.}
  \label{tab:casp11-12}
  \centering

  \begin{tabular}{lccccccccccc}
    \toprule
    & \multicolumn{5}{c}{CASP 11} & & \multicolumn{5}{c}{CASP 12} \\
    \cmidrule(r){2-6} \cmidrule(r){8-12}
     & \multicolumn{2}{c}{Stage 1} & & \multicolumn{2}{c}{Stage 2}      & & \multicolumn{2}{c}{Stage 1} & & \multicolumn{2}{c}{Stage 2}                  \\
    \cmidrule(r){2-3} \cmidrule(r){5-6} \cmidrule(r){8-9} \cmidrule(r){11-12}
    Method     & Glob     & Per & & Glob & Per & & Glob & Per & & Glob & Per \\
    \midrule
    Ours       & \textbf{0.84} & \textbf{0.66} & & \textbf{0.87} & \textbf{0.45} & & \textbf{0.79} & \textbf{0.73} & & \textbf{0.82} & \textbf{0.62} \\
    3DCNN \citep{derevyanko2018deep}   & 0.59 & 0.52 & & 0.64 & 0.40 & & 0.49 & 0.44 & & 0.61 & 0.51 \\
    Ornate \citep{pages2019protein}    & 0.64 & 0.47 & & 0.63 & 0.39 & & 0.55 & 0.57 & & 0.67 & 0.49 \\
    GraphQA \citep{baldassarre2020graphqa}             & 0.83 & 0.63 & & 0.82 & 0.38 & & 0.72 & 0.68 & & 0.81 & 0.61 \\
    VoroMQA \citep{olechnovivc2017voromqa} & 0.69 & 0.62 & & 0.65 & 0.42 & & 0.46 & 0.61 & & 0.61 & 0.56 \\
    SBROD \citep{karasikov2019smooth}      & 0.58 & 0.65 & & 0.55 & 0.43 & & 0.37 & 0.64 & & 0.47 & 0.61 \\
    \midrule
    ProQ3D \citep{uziela2017proq3d}        & 0.80 & \textbf{0.69} & & 0.77 & 0.44 & & 0.67 & 0.71 & & 0.81 & 0.60 \\
    \bottomrule
  \end{tabular}
\end{table}

The CASP 11-12 datasets have been the most well-benchmarked in the recent MQA literature. However, for completeness, we also evaluate our method on CASP 13 (Table \ref{tab:casp13}). Because of its recency, many target structures remain publicly unavailable. We use the stage 2 candidate structures of a subset of 20 targets previously used for benchmarking \citep{baldassarre2020graphqa}. Our method achieves improved results over all other methods, including ProteinGCN \citep{sanyal2020proteingcn}, a more recent graph-based method. Because of the small sample size, we emphasize that these results, although promising, should be considered preliminary until further structures for CASP 13 are available.

\begin{table}
  \caption{Performance on the 20 publicly available CASP 13 targets, stage 2 in terms of global and mean per-target Pearson correlation coefficient. As before, the top structure-only method is in bold, as is the top method overall (if different). Our method outperforms all other methods.}
  \label{tab:casp13}
  \centering
  \begin{tabular}{lccccc}
    \toprule
    Method     & Global     & Mean per-target \\
    \midrule
    Ours   &  \textbf{0.881} & \textbf{0.685} \\
    ProteinGCN \citep{sanyal2020proteingcn}    & 0.723 & 0.603 \\
    VoroMQA \citep{olechnovivc2017voromqa}     & 0.769 & 0.665 \\
    \midrule
    ProQ3D \citep{uziela2017proq3d}           & 0.849 & 0.671  \\
    \bottomrule
  \end{tabular}
\end{table}

\subsection{Protein design}
\label{sec:results_cpd}

Our method achieves state-of-the-art performance on CATH 4.2, representing a substantial improvement both in terms of perplexity and sequence recovery over Structured Transformer \citep{ingraham2019generative}, which was trained on the same training set (Table \ref{tab:ingraham}). Following \citet{ingraham2019generative}, we report evaluation on short (100 or fewer amino acid residues) and single-chain subsets of the CATH 4.2 test set, containing 94 and 103 proteins, respectively, in addition to the full test set. Although Structured Transformer leverages an attention mechanism on top of a graph-based representation of proteins, the authors note in ablation studies that removing attention appeared to increase performance. We therefore retrain and compare against a version of Structured Transformer with the attention layers replaced with standard graph propagation operations that we call Structured GNN. Our method also improves upon this model.

\begin{table}
  \caption{Performance on the CATH 4.2 test set and its short and single-chain subsets in terms of per-residue perplexity (lower is better) and recovery (higher is better). Recovery is reported as the median over all structures of the mean recovery of 100 sequences per structure. Our method performs better than Structured Transformer and a variant of it, Structured GNN, in which we replaced the attention mechanisms with standard graph propagation operations (see \ref{sec:results_cpd}).}
  \label{tab:ingraham}
  \centering
  \begin{tabular}{lccccccc}
    \toprule
     & \multicolumn{3}{c}{Perplexity} & & \multicolumn{3}{c}{Recovery \%}  \\
    \cmidrule(r){2-4} \cmidrule(r){6-8}
    Method     & Short & Single-chain & All & & Short & Single-chain & All \\
    \midrule
    Ours & \textbf{7.10} & \textbf{7.44} & \textbf{5.29} & & \textbf{32.1} & \textbf{32.0} & \textbf{40.2} \\
    Structured GNN       & 8.31 & 8.88 & 6.55 & & 28.4 & 28.1 & 37.3 \\
    Structured Transformer \citep{ingraham2019generative} & 8.54 & 9.03 & 6.85 & & 28.3 & 27.6 & 36.4 \\
    \bottomrule
  \end{tabular}
\end{table}

We emphasize that the underlying architecture of Structured GNN is comparable to ours, with the exception of GVPs and geometric vector channels. In particular, the underlying protein graph is built in a similar manner, except with structural information encoded in solely \textit{scalar} channels. For example, \textit{relative} orientations are encoded in terms of quaternion coefficients (4 scalar channels per \textit{edge}), whereas we represent \textit{absolute} orientations with 3 vector channels per \textit{node}. Therefore, our improvement over Structured GNN is the most direct indication of the benefit of our approach that combines geometric and relational information. Additionally, due to the efficiency of our representation, we are able to achieve this performance gain with a modest but notable decrease in parameter count---1.01M in our model versus 1.38M in Structured GNN and 1.53M in Structured Transformer.

\section{Conclusion}
In this work, we developed the first architecture designed specifically for learning on dual relational and geometric representations of 3D macromolecular structure. At its core, our method, GVP-GNN, augments graph neural networks with computationally fast layers that perform expressive geometric reasoning over Euclidean vector features. Our method possesses desirable theoretical properties and empirically outperforms existing architectures on learning quality scores and sequence designs, respectively, from protein structure.

In further work, we hope to apply our architecture to other important structural biology problem areas, including protein complexes, RNA structure, and protein-ligand interactions. Our results more generally highlight the promise of domain-aware architectures in specialized applications of deep learning, and we hope to continue developing and refining such architectures in the domain of learning from structure.

\section*{Acknowledgements}
We thank Tri Dao and all members of the Dror group for feedback and discussions.

\section*{Funding}
We acknowledge support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program, and Intel Corporation. SE is supported by a Stanford Bio-X Bowes fellowship. RJLT is supported by the U.S. Department of Energy, Office of Science Graduate Student Research (SCGSR) program.

\medskip
\begin{small}
\bibliographystyle{abbrvnat}

\section*{Supplementary Information}

\subsection*{Equivariance and invariance}
The vector and scalar outputs of the GVP are equivariant and invariant, respectively, with respect to an arbitrary composition of rotations and reflections in 3D Euclidean space described by $R$ \textit{i.e.},
\begin{equation*}
\text{GVP}(\left(\mathbf{s}, R(\mathbf{V}))\right) = (\mathbf{s'}, R(\mathbf{V}'))
\end{equation*}
\begin{proof}

We can write the transformation described by $R$ as multiplying $\mathbf{V}$ with a unitary matrix $\mathbf{U} \in \mathbb{R}^{3 \times 3}$ from the right. The L$_2$-norm, scalar multiplications, and nonlinearities are defined row-wise as in Algorithm 1. We consider scalar and vector outputs separately. The scalar output, as a function of the inputs, is
\begin{equation*}
    \mathbf{s'} = \sigma\left( \mathbf{W}_m\begin{bmatrix}\norm{\mathbf{W}_h\mathbf{V}}\\ \mathbf{s}\end{bmatrix} + \mathbf{b}\right)
\end{equation*}
Since $\norm{\textbf{W}_h \mathbf{V} \mathbf{U}} = \norm{\textbf{W}_h \mathbf{V}}$, we conclude $\mathbf{s'}$ is invariant.
Similarly the vector output is
\begin{equation*}
    \mathbf{V'} = \sigma^+\left(\norm{\textbf{W}_\mu \textbf{W}_h \mathbf{V}}\right)\odot \textbf{W}_\mu \textbf{W}_h \mathbf{V}
\end{equation*}
The row-wise scaling can also be viewed as left-multiplication by a diagonal matrix $\mathbf{D}$. Since $\norm{\textbf{W}_\mu \textbf{W}_h \mathbf{V}} = \norm{\textbf{W}_\mu \textbf{W}_h \mathbf{VU}}$, $\mathbf{D}$ is invariant. Since
\begin{equation*}
    \mathbf{DW}_\mu\mathbf{W}_h\mathbf{(VU)} = \left(\mathbf{DW}_\mu\mathbf{W}_h\mathbf{V}\right)\mathbf{U}
\end{equation*}
we conclude that $\mathbf{V'}$ is equivariant.
\end{proof}

\subsection*{The ability to approximate rotation-invariant functions}

The GVP inherits an analogue of the Universal Approximation property \cite{cybenko1989approximation} of standard dense layers. If $R$ describes an arbitrary rotation or reflection in 3D Euclidean space,\footnote{More precisely, if $R$ describes a unitary transformation.} we show that a GVP can approximate arbitrary scalar-valued functions invariant under $R$ and defined over $\Omega^\nu \subset \mathbb{R}^{\nu\times 3}$, the bounded subset of $\mathbb{R}^{\nu \times 3}$ whose elements can be canonically oriented based on three linearly independent vector entries. Without loss of generality, we assume the first three vector entries can be used.

The machinery corresponding to such approximations corresponds to a GVP $G_s$ with only vector inputs, only scalar outputs, and a sigmoidal nonlinearity $\sigma$, followed by a dense layer. This can also be viewed as the sequence of matrix multiplication with $\mathbf{W}_h$, taking the L$_2$-norm, and a dense network with one hidden layer. Such machinery can be extracted from any two consecutive GVPs (assuming a sigmoidal $\sigma$). The theorem is restated from the main text:

\begin{theorem}
    Let R describe an arbitrary rotation or reflection in $\mathbb{R}^3$.
    For $\nu \ge 3$ let $\Omega^\nu \subset \mathbb{R}^{\nu \times 3}$ be the set of all $\mathbf{V} = \begin{bmatrix}\mathbf{v}_1, & \ldots, & \mathbf{v}_\nu\end{bmatrix}^T\in \mathbb{R}^{\nu \times 3}$ such that $\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3$ are linearly independent and $0<\norm{\mathbf{v}_i} \le b$ for all $i$ and some finite $b>0$. Then for any continuous $F: \Omega^\nu \rightarrow \mathbb{R}$ such that $F(R( \mathbf{V})) = F(\mathbf{V})$ and for any $\epsilon>0$, there exists a form $f(\mathbf{V})=\mathbf{w}^TG_s(\mathbf{V})$ such that $|F(\mathbf{V})-f(\mathbf{V})|< \epsilon$ for all $\mathbf{V}\in\Omega$.
\end{theorem}
\begin{proof}
    The idea is to write $F$ as a composition $F=\Tilde{F} \circ \omega$ and $\omega = h \circ y$.  We show that multiplication with $\mathbf{W}_h$ and and taking the L$_2$-norm can compute $y$, and that the dense network with one hidden layer can approximate $\Tilde{F} \circ h$.

    Call an element $\mathbf{V} \in \Omega^\nu$ \textit{oriented} if $\mathbf{v}_1 = x_1 \mathbf{e}_x$, $\mathbf{v}_2 = x_2\mathbf{e}_x + y_2\mathbf{e}_y$, and $\mathbf{v}_3 = x_3\mathbf{e}_x + y_3\mathbf{e}_y+z_3\mathbf{e}_z$, with $x_1, y_2, z_3>0$.
    Define $\omega: \Omega^\nu \rightarrow \mathbb{R}^{3\nu-3}$ to be the \textit{orientation function} that orients its input and then extracts the vector of $3\nu-3$ coefficients, $[x_1, x_2, y_2, x_3, y_3, z_3, \ldots, x_i, y_i, z_i, \ldots]^T$. These elements can be written as
    \begin{eqnarray*}
        x_1 &=& \norm{\mathbf{v}_1}  \\
        x_i &=& \mathbf{v}_i \cdot \mathbf{v}_1/x_1, \quad i \ge 2\\
        y_2 &=& \sqrt{\norm{\mathbf{v}_2}^2 -x_2^2} \\
        y_i &=& (\mathbf{v}_i \cdot \mathbf{v}_2 - x_ix_2)/y_2, \quad i \ge 3\\
        z_3 &=& \sqrt{\norm{\mathbf{v}_3}^2 - x_3^2 - y_3^2}\\
        z_i &=& (\mathbf{v}_i \cdot \mathbf{v}_3 - x_ix_3 - y_iy_3)/z_3, \quad i \ge 4
    \end{eqnarray*}
    and are invariant under rotation and reflection, because they are defined using only the norms and inner products of the $\mathbf{v}_i$. Then $F = \tilde{F}\circ \omega$, where $\tilde{F}: [-b, b]^{3\nu-3}\rightarrow \mathbb{R}$.

    The key insight is that if we construct $\mathbf{W}_h$ such that the rows of $\mathbf{W}_h\mathbf{V}$ are the original vectors $\mathbf{v}_i, \forall i$ and all differences $\mathbf{v}_i - \mathbf{v}_j, \forall i, j\le \min(i,3)$, then we can compute $\omega(\mathbf{V})$ from the row-wise norms of $\mathbf{W}_h\mathbf{V}$. That is, $\omega = h \circ y$ where $\mathbf{y} = y(\mathbf{V}) = \norm{^.} \odot (\mathbf{W}_h\mathbf{V})\in\mathbb{R}^{4\nu-6}$ and $h$ is an application of the cosine law. The GVP precisely computes $\mathbf{y}$ as an intermediate step: $G_s(\mathbf{V}) = \sigma \odot (\mathbf{W}_m\mathbf{y} +\mathbf{b})$.
    It remains to show that there exists a form $\tilde{f}(\mathbf{y}) = \mathbf{w}^T[\sigma \odot(\mathbf{W}_h\mathbf{y} + \mathbf{b})]$ that $\epsilon$-approximates $\tilde{F} \circ h : [-2b, 2b]^{4\nu-6}\rightarrow\mathbb{R}$. Up to a translation and uniform scaling of the hypercube, this is the result of the Universal Approximation Theorem \citep{cybenko1989approximation}.
\end{proof}


\section{Introduction}
American Sign Language (ASL) is the fourth most commonly used language in the United States and is the language most commonly used by Deaf people in the United States and the English-speaking regions of Canada  \cite{ASL-Info-2001}. Unfortunately, until recently, ASL received little research. This is due, in part, to its delayed recognition as a language until William C. Stokoe's publication in 1960 \cite{Stokoe-1960}. Limited data has been a long-standing obstacle to ASL research and computational modeling. The lack of large-scale datasets has prohibited many modern machine-learning techniques, such as Neural Machine Translation, from being applied to ASL. In addition, the modality required to capture sign language (i.e. video) is complex in natural settings (as one must deal with background noise, motion blur, and the curse of dimensionality). Finally, when compared with spoken languages, such as English, there has been limited research conducted into the linguistics of ASL.

We realize a simplified version of Liddell and Johnson's Movement-Hold (MH) Model \cite{Movement-Hold-1989} using a Probabilistic Graphical Model (PGM). We trained our model on ASLing, a dataset collected from three fluent ASL signers. We evaluate our PGM against other models to determine its ability to model ASL. Finally, we interpret various aspects of the PGM and draw conclusions about ASL phonetics. The main contributions of this paper are

\begin{enumerate} [topsep=0pt, itemsep=0pt,partopsep=1ex,parsep=1ex]
    \item A PGM which models the the body configuration in ASL
    \item A video dataset containing over 800 phrases comprising of over 4,000 signs.
\end{enumerate}

\subsection{The Movement-Hold Model}
The Stokoe Model was the first attempt to analyze the internal structure of signs and challenge the notion that signs are holistic iconic gestures. Stokoe analyzed signs using three parameters: location, handshape, and movement. He assumed each parameter was produced simultaneously \cite{txtbk-Stokoe-Model}.
The foundational Stokoe's model is unable to recognize sequential contrasts occurring between signs because it does not require the parameters of a sign to be produced in a prescribed order \cite{txtbk-Stokoe-Problems}.

In response to the limitations of the Stokoe Model, Liddell and Johnson published their Movement-Hold (MH) Model.
The MH Model partitions each sign into sequential hold and movement segments. A bundle of articulatory features is associated with each segment and contains information regarding the segment's handshape, location, palm orientation, and non-manual markers (e.g. facial expressions).
A segment with no changes to the articulation bundle is a hold segment and can be represented with an H for a full-length hold or an X for a shorter hold. A segment in which some of the articulatory features are transitioning is a movement segment and can be represented with an M. The MH Model recognizes at least nine sign structures, including Hold (H), X M H, and H M X M H, but not H M. The MH Model was novel in its approach to modeling the sequential aspect of signs \cite{txtbk-Movement-Hold-Model}.

The MH Model is a significant departure from previous models because its basic units of movements and holds are produced sequentially, allowing for precise representation of sequences within signs.

\begin{figure}[H]
\small
    \centering

    \frame{\begin{tabular}{cccccccc}
    $H_1$ & $H_1$ & $H_1$ & $H_2$ & $H_2$ & $E$ & $E$ & $E$\\
    $H_1$ & $H_1$ & $H_1$ & $E$ & $E$ & $E$ & $E$ & $E$\\
    $H_1$ & $H_2$ & $H_3$ & $H_4$ & $H_4$ & $H_4$ & $E$ & $E$
    \end{tabular}}
    \caption{Some example sign structures realized by our model. Each sign is represented in eight frames instead of twenty-five for easier viewing. Each $H_i$ is a distinct body configuration with $E$ being an ending token that is used for padding. The first example illustrates an $HMH$. The second demonstrates an $H$. Finally, the third illustrates an $XMXMXMH$}
    \label{fig:MH_markov}
\end{figure}

\subsection{Prior Works}
While many computational approaches to modeling sign language, including the previously published Interspeech sign language recognition paper \cite{Dreuw-2007} focus on representing signs holistically. A prior approach from speech recognition that gained popularity within sign language modeling is the use of the Hidden Markov Models (HMM) and its variants to model subunits of signs.  Influenced by the MH Model, Theodorakis et al. \cite{Theodorakis-2014} used HMMs to represents signs which they segmented into dynamic and static portions.
Fang et al. \cite{Fang-2004} modeled Chinese Sign Language using HMMs and temporal clustering.
Vogler et al. \cite{Vogler-1999} trained HMMs to recognize movement and hold phonemes in an attempt to increase the scalability of ASL recognition systems.
Bauer et al. \cite{Bauer-2001} take a different approach in their attempt to create an automatic German Sign Language recognition system, using a K-means algorithm to define subunits.
With all four of these papers focusing on creating a sign language recognition system, there is more focus on recognition success over the phonological model. However, those algorithms do not always identify recognizable phonemes.
For this reason, building from the ground up, we focus specifically on creating a phonological model and interpreting the model. Additionally, our model is entirely unsupervised, probabilistic, and capable of being trained end-to-end.

\begin{figure*}[ht!]
    \centering
    \subfloat{{\includegraphics[width=5cm]{figures/Edith.png} }}%
    \subfloat{{\includegraphics[width=5cm]{figures/caroline.png} }}%
    \subfloat{{\includegraphics[width=5cm]{figures/michael.png} }}%
    \caption{Example frames from our database, ASLing. The green points indicate the features extracted by our Convolutional Pose Machine: the algorithm used to extract body keypoint locations.}%
    \label{fig:example}%
\end{figure*}

\section{Methodology}

We aim to computationally model the phonetic structure within ASL. We restrict ourselves to a simplified form of the MH Model. Our primary deviations from the MH Model are
\begin{enumerate}[topsep=0pt, itemsep=0pt,partopsep=1ex,parsep=1ex]
    \item we are modelling the global positions (specifically the wrist, elbow, shoulder, and head positions), and
    \item we are explicitly modeling holds and implicitly modelling movement between the holds and compressing it into a single frame.
\end{enumerate}

\subsection{Dataset}
ASLing is an American Sign Language corpus produced by three fluent signers of American Sign Language. Two signers identify as Deaf, and one signer identifies as Hard of Hearing.
At the time of writing, ASLing contained over 800 phrases. The dataset had over 1,200 signs, with each sign produced an average of about 3 times. The phrases were chosen because they each contained at least one word with a high frequency in English. Each sign was glossed and each phrase was translated into English. Each phrase was tagged with a level of noise (none, low, medium, high, or broken) that described the performance of our pose predictor. The signers were instructed to record with as much of their upper body in frame as possible, in a location with good lighting, and while wearing contrastive clothing.
The data was collected in phrases instead of individual signs to allow for more natural signing.

\subsection{Feature Extraction and Data Preparation}

The data collected was stored in a video format. Each sign was extracted from the videos (resulting in $M=1515$ signing portions). From each frame of the signing portions, 2D estimates of key points were extracted using a Convolutional Pose Machine \cite{CPM} (specifically the head, shoulders, elbows, and wrist locations). These poses were then normalized (setting the head to the origin and setting the unit to the average distance between the head and shoulders) resulting in $D=14$ features. Finally, the sign were padded to a length of $P=25$ frames with a vector of zeros (as this was the length of the longest sign). Note that the vector of zeros was used as the end token.

\subsection{Phonetic Model}

Through the MH Model, the production of signs can be treated as a sequential process (see Figure \ref{fig:MH_markov} for an example). The entire process can be posed as a Probabilistic Graphical Model (PGM).
PGMs have several benefits, including their ease of training and interpretability. Additionally, PGMs can be used to generate synthetic data by sampling from the model \cite{pgm}.
We chose specifically to use a Dynamic Bayesian Network to model the sequential aspect of signing. Specifically, let
{\footnotesize
\begin{align*}

    N = &  \text{The number of body configurations prototypes}\\
    M = &  \text{The number of signs}\\
    P = &  \text{The number of frames in each sign}\\
    D = &  \text{The dimensionality of frame features}\\

    \mu_i = &  \text{The $i^{th}$ idealized body configuration prototypes}\\
    \sigma = &  \text{The amount of dispersion around each dimension from the}\\  &\text{idealized body configuration prototypes during sign production}\\
    \pi = &  \text{The probability of each body configuration prototype being at}\\ & \text{the start of a sign}\\
    T_{i, j} = &  \text{The transition probability of moving from body configuration}\\ & \text{prototype $i$ to body configuration prototype $j$}\\

    c^w_i = &  \text{The body configuration prototype chosen for frame $i$ in sign $w$}\\
    x^w_i = &  \text{The observed production of the chosen body configuration}\\ & \text{prototype for frame $i$ in sign $w$}
\end{align*}

}%

In addition to these, we have the hyperparameters $\alpha, \mu_\mu, \mu_\sigma, \sigma_\mu, \sigma_\sigma$. We specifically chose $\alpha = 1, \mu_\sigma = 1, \sigma_\sigma = 10, \mu_\mu = 0, \sigma_\mu = 10$.
Our full model can then be stated as
\begin{align}
\sigma | \mu_\sigma, \sigma_\sigma  & \sim LogNormal(\mu_\sigma, \sigma_\sigma)\\
T_n  | \alpha & \sim Dir(\alpha \mathbf{1}_N)\\
\mu_0  & = \mathbf{0}_D\\
\mu_{n>0} | \mu_\mu, \sigma_\mu  & \sim \mathcal{N}(\mu_\mu, Diag(\sigma_\mu))\\
\pi | \alpha & \sim Dir(\alpha \mathbf{1}_N)\\
c^w_{f=0}  | \pi & \sim Categorical(\pi)\\
c^w_{f>0} | T, c^w_{f-1} & \sim Categorical(T_{c^w_{f-1}})\\
x^w_f | \mu, c^w_f, \sigma & \sim \mathcal{N}(\mu_{c^w_f}, Diag(\sigma))
\end{align}
This process is shown in the plate diagram in Figure \ref{fig:plate}. Note that the 0th body configuration is defined as a special end phoneme (which is indicated by the 0 feature vector).
\begin{figure}[H]
    \centering

    \begin{tikzpicture}

        \node[latent] (pi) {$\pi$};
        \node[latent, right=of pi] (T) {$T$};

        \node[left=0.5cm of pi] (alphapi) {$\alpha$};
        \node[right=0.5cm of T] (alphaT) {$\alpha$};

        \node[latent, below=0.5cm of pi] (c0) {$c^w_0$};
        \node[obs, below=.75cm of c0] (x0) {$x^w_0$};
        \node[latent, right=of c0] (c1) {$c^w_1$};
        \node[obs, below=.75cm of c1] (x1) {$x^w_1$};
        \node[latent, right=of c1] (c2) {$c^w_2$};
        \node[obs, below=.75cm of c2] (x2) {$x^w_2$};

        \node[right=0.2 of c2] (ellipsec) {\ldots};
        \node[right=0.2 of x2] (ellipsex) {\ldots};

        \node[above right = 0.2 and 0.2 of ellipsec] (nulla) {};
        \node[below left  = 0.2 and 0.2 of x0] (nullb) {};

        \node[latent, below=1.25cm of x0] (mu) {$\mu$};
        \node[latent, below=1.25cm of x2] (sigma) {$\sigma$};
        \node[left=0.5cm of sigma] (musigma) {$\mu_\sigma$};
        \node[right=0.5cm of sigma] (sigmasigma) {$\sigma_\sigma$};
        \node[left=0.5cm of mu] (mumu) {$\mu_\mu$};
        \node[right=0.5cm of mu] (sigmamu) {$\sigma_\mu$};

        \draw[every loop]

            (musigma) edge node {} (sigma)
            (sigmasigma) edge node {} (sigma)
            (alphapi) edge node {} (pi)
            (mumu) edge node {} (mu)
            (sigmamu) edge node {} (mu)
            (alphaT) edge node {} (T)

            (pi) edge node {} (c0)
            (T) edge node {} (c1)
            (T) edge node {} (c2)
            (c0) edge node {} (c1)
            (c1) edge node {} (c2)

            (c0) edge node {} (x0)
            (c1) edge node {} (x1)
            (c2) edge node {} (x2)
            (mu) edge node {} (x0)
            (mu) edge node {} (x1)
            (mu) edge node {} (x2)
            (sigma) edge node {} (x0)
            (sigma) edge node {} (x1)
            (sigma) edge node {} (x2);

    \plate [] {prototypes} {(mu)} {$N$}; %
    \plate [] {dimensions} {(sigma)} {$D$}; %
    \plate [] {signs} {(c0)(c1)(c2)(x0)(c1)(c2)(ellipsec)(ellipsex)(nulla)(nullb)} {$M$};
    \plate [] {frames} {(c0)(c1)(c2)(x0)(c1)(c2)(ellipsec)(ellipsex)} {$P$};
    \end{tikzpicture}
    \caption{The plate diagram of our PGM. The non-shaded nodes are latent variables. Shaded nodes are observed variables. $\alpha, \mu_\mu, \sigma_\mu, \mu_\sigma, \sigma_\sigma$ are hyperparameters.}
    \label{fig:plate}
\end{figure}

\begin{figure}
    \centering
    \subfloat{{\includegraphics[width=3cm, trim={0.7cm 0.7cm 0 0},clip]{figures/all.png} }}%
    \subfloat{{\includegraphics[width=3cm, trim={0.7cm 0.7cm 0 0},clip]{figures/poses.png} }}%
    \caption{A sampling of body poses exhibited in the data (left) and the learned body configuration prototypes from the 5 body configuration prototype model (right).}
    \label{fig:my_label}
\end{figure}

\subsubsection{Training}

We wish to obtain estimates of $\theta = \pi, T, \mu, \sigma$. With these parameters, one can run the generative process in full and derive maximum a posteriori (MAP) estimates of the model conditioned on data. As we wish only to obtain MAP estimates of these, we will use Expectation Maximization (EM).

To do this, we follow the standard EM procedure:
\begin{enumerate} [topsep=0pt, itemsep=0pt,partopsep=1ex,parsep=1ex]
    \item Calculate the optimal cluster assignments given the data and our current estimation of $\theta$.
    \item Calculate an updated (MAP) estimate of $\theta$ using our cluster assignments and prior on $\theta$.
\end{enumerate}

To compute the optimal cluster assignments, we calculate $P(Z | X)$:
\begin{align}
    P(c_0^w = i | x_0^w) &\propto e^{\frac{-1}{2} (x_0^w-\mu_i)^T diag(\sigma)^{-1} (x_0^w-\mu_i)} \pi_i\\
    P(c_f^w = i | c_{f-1}^w = j, x_f^w) &\propto e^{\frac{-1}{2} (x_f^w-\mu_i)^T diag(\sigma)^{-1} (x_f^w-\mu_i)} T_{j, i}
\end{align}
From this, one can trivially compute the optimal cluster assignments ($Z = c$) given $\theta$ in $\mathcal{O}(MPN)$.

To compute the optimal estimate of $\theta$ given our cluster assignments and priors we use calculate $P(\theta | X, Z)$
\begin{align}
    P(\pi | c_0^w = i) \propto& \frac{\pi_i}{B(\alpha)} \prod_j \pi_j^{\alpha_j - 1}\\
    P(T_j | c_f^w = i, c_{f-1}^w = j) \propto& \frac{T_{j, i}}{B(\alpha)} \prod_k \left(T_{j, k}\right)^{\alpha_k - 1} \\
    P(\mu_i | x_f^w, c_f^w = i) \propto& e^{\frac{-1}{2} (x_f^w-\mu_i)^T diag(\sigma)^{-1} (x_f^w-\mu_i)} * \nonumber\\
                                       & e^{\frac{-1}{2} (\mu_i-\mu_\mu)^T diag(\sigma_\sigma)^{-1} (\mu_i^w-\mu_\mu)}\\
    P(\sigma |\mu_i, x, c) \propto& P(x, c, \mu | \sigma) * \nonumber\\
                                  &\prod_i \frac{1}{\sigma_i \sigma_{\sigma_i}} exp\left(-\frac{\ln(\sigma_i) - \mu_{\sigma_i}}{2 \sigma_{\sigma_i}}\right)
\end{align}

where $P(x, c, \mu | \sigma)$ was left unexpanded in the interest of space. These can each be optimized numerically (we chose to optimize these in the order of $\pi, T, \mu, \sigma$). These two steps were repeated until convergence.

\begin{figure}[h]
    \centering
    \includegraphics[height=4cm]{figures/viz.png}%
    \caption{A T-SNE visualization of the various models presented in this paper}
    \label{fig:tsne}
\end{figure}

\section{Results and Evaluation}

After this model was trained to convergence, we performed the following analysis
\subsection{Evaluation}

\begin{table*}[h]
    \centering
    \begin{tabular}{|c|cccc|}
    \hline
        Model & Ours & GMM & GMM-LDA & VAE \\
    \hline
        Binary Cross Entropy of Descriminator & $\mathbf{0.65 \pm 0.01}$ & $0.05 \pm 0.01$ & $0.03 \pm 0.01$ & $0.04 \pm 0.02$\\
    \hline
    \end{tabular}
    \caption{Binary Cross Entropy (BCE) of a Discriminator trained on distinguishing Synthetic data from Real data. A low BCE means the synthetic generator was unable to adequately fool a discriminator. Note the ability to fool a discriminator corresponds with the ability to produce life-like signs.}
    \label{tab:bce}
	\vspace{-0.6cm}
\end{table*}

We chose to evaluate our model by determining how capable it was of producing signs. We compared our model against several other models:

\begin{enumerate} [topsep=0pt, itemsep=0pt,partopsep=1ex,parsep=1ex]
    \item A GMM-LDA model of sign language
    \item A GMM model of sign language
    \item A (Seq2Seq) Variational Auto Encoder (VAE) \cite{vae}
\end{enumerate}

It should be noted that there was not enough data to successfully train a generative adversarial network\cite{gan} (GAN). The former two models can be seen as ablations of our model.  The VAE was chosen to determine how our model compared with similar deep models. The probabilistic models are described in further detail at the end of this section.

Each of these algorithms were evaluated on its ability to fool a descriminator. We chose to use a simple neural network as a descriminator (specifically a GRU network). This network was trained on the binary problem of distinguishing between generated data (from the generative processes listed above) and real data. The results are shown in Table \ref{tab:bce}.
In this evaluation, our model outperformed all other models by a wide margin. This was likely because the ablated versions were incapable of capturing the dynamic nature of ASL.

In addition to these quantitative experiments, we also used T-SNE\cite{tsne} to project simulated data from each model to two dimensions (shown in Figure \ref{fig:tsne}). From this, we can observe that the LDA-GMM and GMM models underfit (likely due to the high bias of the model). It is also evident that the VAE did not fit the data well. Finally, one can observe that our model produced several signs which did not coincide with real data. This was likely due to our model implicitly modeling movement within the MH Model. Even so, all signs present in the corpus had similar signs generated by our model. This implies that our model was capable of modeling all body configurations present in the corpus.

\subsubsection{GMM-LDA and GMM}
The GMM-LDA model we used is a topic model\cite{lda} over the body configurations used within a sign. In this model, the topics are distributions of body configuration prototypes and the words are the body configuration prototypes themselves. The GMM model is similar to the GMM-LDA. Their model specifications are given in Figure \ref{fig:gmmodels}.
\begin{figure}
	\centering
	\begin{tabular}{cc}
		\frame{\includegraphics[width=0.4\linewidth]{figures/LDAGMM.png}} &
		\frame{\includegraphics[width=0.4\linewidth]{figures/GMM.png}}
	\end{tabular}
	\caption{The model specifications for the LDA-GMM (left) and GMM (right). $f$ is the frame index, $t$ is the topic index, and $T$ is the number of topics. $\alpha, \beta$ are concentration hyperparameters, $\pi$ is the relative frequency of the different body configuration prototypes, and $\psi, \theta$ are the word distributions of the topics and the topic choice for the signs respectively}
	\label{fig:gmmodels}
	\vspace{-0.5cm}
\end{figure}

\subsection{Interpretation}

To interpret the model, we chose to fit a simplified model containing only five body configuration prototypes (including the end state). Note that without loss of generality, this section assumes the signer is right-hand dominant.
The body configuration prototypes discovered appear to be the following body configurations (where the listed location is the location of articulation for the right hand): end state, head, right chest, left hand, and neutral signing space (i.e. the location used for fingerspelling). We examined $\pi$ to determine that the most common starting position was the right chest area. The least common starting position was the head (we are not counting the end state in the least common starting positions).
In addition, one can determine the average length of each hold. This can be done by examining the self transition probabilities $T[j, j]$ and running Geometric Trials (where the first transition away from the given state is a success):

\begin{align}
\mathbb{E}(\text{length of hold }i) = \frac{1}{1 - T[i, i]}
\end{align}

From this, one can determine most holds last for approximately 8.16 frames (with the standard deviation between the lengths of the different body configuration prototypes being only 0.98 frames). This means that each hold is approximately 800 milliseconds. Note that the holds derived by the non-simplified model are substantially shorter (in the order of 400 milliseconds).
We can also determine the popularity of each body configuration prototype by evaluating the expected number of times each will arise within a sign:
\begin{align}
\mathbb{E}(\text{count of }j) = \sum_{i=0}^{19} P(c_i = j) = \sum_{i=0}^{19} (\pi T^i)[j]
\end{align}

This leads to the conclusion that the ``ending state" will be be in an average of 11 frames. Additionally, the most common pose is the right chest (showing up an average of 2.91 frames per sign) and the least common is the head (showing up an average of 1.86 times per sign).
We can also examine the probability of any given body configuration prototype ending a sign by examining the probability of transitioning to the 0th prototype. Therefore, the right chest is the most common ending body configuration prototype and the left hand is the least common ending body configuration prototype.
Finally, we can examine the amount of variation during production between each joint compared to the body configuration prototype. From examining this, one can observe that the head and shoulder exhibit the least variation. The elbows and left hand exhibit slightly more variation. The right hand exhibits the most variation. This is to be expected as the right hand performs the most motion when signing.

\begin{figure}
    \centering
    \subfloat{{\includegraphics[height=3.5cm]{figures/initials.png} }}%

    \subfloat{{\includegraphics[height=3.5cm]{figures/transition_log.png} }}%
    \caption{The initial probabilities and (natural) log of transition probabilities learned by our simplified phonetic model. Note the high probability of self transitions. Self transitions mean that one stayed in the body configuration prototype and thus naturally maps to the concept of a Hold (see Figure \ref{fig:MH_markov})}
    \label{fig:probs}
\end{figure}

\section{Conclusions}

In this paper we successfully modeled the phonetic structure of ASL using a simplified Movement-Hold Model. We evaluated our model against multiple ablated versions and a deep counterpart. Finally, we interpreted various aspects of the model to demonstrate its utility and interpretability. Nevertheless, our model did have several limitations. First, we did not explicitly model movement. This likely caused the artifacts seen in Figure \ref{fig:tsne}. Additionally, our model did not model the local movement of the hand (specifically handshape and palm orientation). We plan on not only addressing these issues in future works, but also modeling non-manuals and incorporating 3D to allow for easier interpretation.

\bibliographystyle{IEEEtran}


\section{Introduction}\label{sec:introduction}
Recently, the state $X(4140)$ was confirmed in the measurement of the $B^+ \to J/\psi \phi K^+$ reaction at LHCb~\cite{Aaij:2016nsc,Aaij:2016iza}. However, one surprise is that the $X(4140)$ deduced from the analysis, with quantum numbers $J^{PC}=1^{++}$, has a width $\Gamma \approx 83\pm 21^{+21}_{-14}$~MeV, substantially larger than that claimed in the former experiments~\cite{PDG2016}.

In this work, we show that the low invariant $J/\psi \phi$ mass region,  requires the contributions of a narrow $X(4140)$, and a wide $X(4160)$ resonance which couples to $J/\psi \phi$ but is mostly made by a $D_s^* \bar{D}_s^*$ molecule as discussed in Ref.~\citenum{Molina:2009ct}. As a consequence of analyticity and driven by the molecular nature of $X(4160)$, the $J/\psi\phi$ mass spectrum develops the strong cusp around the $D_s^* \bar{D}_s^*$ threshold~\cite{Aaij:2016nsc,Aaij:2016iza}.

\section{The formalism}
\label{sec: formalism}
\begin{figure}
\includegraphics[width=0.48\textwidth]{BtoDsDs.eps}
\includegraphics[width=0.48\textwidth]{BtoJpsiphi_quarklevel.eps}
\caption{Microscopic quark level production of $K^- D_s^*\bar{D}_s^*$ (left) and $B^- \to K^- J/\psi \phi$ (right) in $B^-$ decay.}  \label{fig:weakdecay}
\end{figure}
\begin{figure}
\includegraphics[width=0.6\textwidth]{BtoDsDs_FSI.eps}
\includegraphics[width=0.35\textwidth]{BtoJpsiphi_FSI.eps}
\caption{(a) The tree diagram, and (b) the $D_s^*\bar{D}_s^*$ final state interaction for $B^- \to K^- D_s^*\bar{D}_s^*$ in the presence of the $X(4160)$ resonance. (c) The mechanism for $B^- \to K^- J/\psi \phi$ driven by the $X(4160)$ resonance.}  \label{fig:DsDs}
\end{figure}

As the first step, the dominant process of $B^- \to K^- D_s^*\bar{D}_s^*$ (we take the complex conjugate reaction to deal with $b$ quark rather than $\bar{b}$ quark) at the quark level proceeds as shown in the left panel of Fig.~\ref{fig:weakdecay}, involving external emission. Obviously in the neighborhood of the resonance the tree level term of Fig.~\ref{fig:DsDs}(a) is small compared to the resonant term of Fig.~\ref{fig:DsDs}(b), but we keep it in the calculations. For the production of $J/\psi \phi$ with this mechanism, the tree level of Fig.~\ref{fig:DsDs}(a) does not contribute and then we are led to the diagram of Fig.~\ref{fig:DsDs}(c)

The $B^- \to K^- J/\psi \phi$ reaction can also proceed through the mechanism as shown in the right panel of Fig.~\ref{fig:weakdecay}
involving internal conversion. Yet, the internal conversion is penalized by color factors with respect to the external emission of the left diagram of Fig.~\ref{fig:weakdecay}, and hence this term, or rescattering of this term like that in Fig.~\ref{fig:DsDs}(c), but with $J/\psi \phi$ intermediate state instead of
 $D_s^*\bar{D}_s^*$, which would involve the extra factor $g_{J/\psi\phi}/g_{D_s^*\bar{D}_s^*}$ versus the amplitude of \ref{fig:DsDs}(c), can be safely neglected.

We can then write the amplitude for the $B^-\to K^-D_s^*\bar{D}_s^*$,
\begin{eqnarray}
\mathcal{M}_{B\to K D_s^*\bar{D}_s^*} &=& A\left[ 1+G_{D_s^*\bar{D}_s^*}(M_{\rm inv}(D_s^*\bar{D}_s^*)) \right. \nonumber \\
&&\times \left. t_{D_s^*\bar{D}_s^*\to D_s^*\bar{D}_s^*}(M_{\rm inv}(D_s^*\bar{D}_s^*))\right] \left( \vec{\epsilon}\cdot \vec{k}\, \vec{\epsilon}^{\,\prime}\cdot \vec{k}-\frac{1}{3}\vec{k}^2\vec{\epsilon}\cdot \vec{\epsilon}^{\,\prime}\right), \label{eq:amp1}
\end{eqnarray}
and  the amplitude for the $B^-\to K^-J/\psi\phi$,
\begin{eqnarray}
\mathcal{M}_{B\to K J/\psi\phi}
 &=& A\times G_{D_s^*\bar{D}_s^*}(M_{\rm inv}(J/\psi\phi))\nonumber \\
&& \times  t_{D_s^*\bar{D}_s^*\to J/\psi\phi}(M_{\rm inv}(J/\psi\phi)) \left( \vec{\epsilon}\cdot \vec{k}\, \vec{\epsilon}^{\,\prime}\cdot \vec{k}-\frac{1}{3}\vec{k}^2\vec{\epsilon}\cdot \vec{\epsilon}^{\,\prime}\right) \nonumber \\
&& + \frac{B \, M^3_{X(4140)}(\vec{\epsilon}_{J/\psi} \times \vec{\epsilon}_\phi) \cdot \vec{k}}{M^2_{\rm inv}(J/\psi\phi)-M^2_{X(4140)}+iM_{X(4140)}\Gamma_{X(4140)}} \nonumber  \\
 &=& \mathcal{M}^{X(4160)} +\mathcal{M}^{X(4140)}\label{eq:amp2}
\end{eqnarray}
where $\vec{\epsilon}$, $\vec{\epsilon}^{\,\prime}$ are the polarization vectors of $D_s^*$ and $\vec{D}_s^*$ (or $J/\psi$ and $\phi$), and we evaluate it in the frame of reference where the $D_s^*\bar{D}_s^*$ (or $J/\psi\phi$) system is at rest. For the $J/\psi\phi$ production, we have included the term $\mathcal{M}^{X(4140)}$ to account for the production of $J/\psi\phi$ via the $1^{++}$ $X(4140)$ resonance. Here we take $M_{X(4140)}=4132$~MeV and $\Gamma_{X(4140)}=19$~MeV. The $A$ and $B$ are free parameters to be fitted to the data. The details of the loop functions $G$ and the transition amplitude $t$ in Eqs.~(\ref{eq:amp1}, \ref{eq:amp2}) are given in Ref.~\citenum{Wang:2017mrt}. Finally, the mass distribution can be written as,
\begin{eqnarray}
\frac{d\Gamma}{dM_{\rm inv}} &=& \frac{1}{(2\pi)^3}\frac{1}{4M_{B^-}^2}|\vec{k}'|\, \tilde{p}|\mathcal{M}|^2, \label{eq:dwidth}
\end{eqnarray}
where $\vec{k}'$ is the $K^-$ momentum in the $B^-$ rest frame, and $\tilde{p}$ the $D_s^*$ (or $J/\psi$) momentum in the $D_s^*\bar{D}_s^*$ (or $J/\psi\phi$) rest frame.

\section{Results}
A suitable fit to the data is obtained as shown in the left panel of Fig.~\ref{fig:result_Jpsiphi}, Our results (red solid curve, labeled as `Full') are in good agreement with the experimental data at low $J/\psi\phi$ invariant masses, and we also present the uncertainties of our results in the figure.
As we can see in the figure, we obtain a contribution from the $X(4140)$ (blue dotted curve) that is dominant at low invariant masses, and is responsible for the peak observed in the experiment around 4135~MeV. The $X(4160)$  (green dashed curve) is responsible for most of the strength and produces a broader peak around 4170~MeV. And finally a cusp appears at the $D_s^*\bar{D}_s^*$ threshold as it shows up in the experiment. This cusp comes from the factor $G_{D_s^*\bar{D}_s^*}(M_{\rm inv})$ and reflects the analytical structure of this function with a discontinuity of the derivative at threshold. A similar feature is also found in Refs.~ \citenum{Dai:2018nmw,Wang:2018djr}. One must stress that this factor appears here as a consequence of analyticity, with a large strength due to the $D_s^*\bar{D}_s^*$ molecular structure of the $X(4160)$. In an analysis like the one of Refs.~\citenum{Aaij:2016nsc,Aaij:2016iza}, where a sum of amplitudes for resonance excitation and some background are fitted to the data, this factor is not considered, and as a consequence the cusp around $D_s^*\bar{D}_s^*$ in the data is missed in the fit.

\begin{figure}
\includegraphics[width=0.48\textwidth]{Dwidth_jpsiphi.eps}
\includegraphics[width=0.48\textwidth]{Dwidth_DssDss.eps}
\caption{Left: the $J/\psi\phi$ mass distribution of the  $B^- \to K^- J/\psi \phi$ decay. Right: The $D_s^*\bar{D}_s^*$ mass distribution of the $B^- \to K^- D_s^*\bar{D}_s^*$ decay. The bands come considering the errors in the fitted parameters, and represents the 68\% confidence-level.}
\label{fig:result_Jpsiphi}
\end{figure}

To finish the work, and as a test of the explanation given here, we present in the right panel of Fig.~\ref{fig:result_Jpsiphi} the $D_s^*\bar{D}_s^*$ mass distribution above threshold obtained with the same parameters as in Fig.~\ref{fig:result_Jpsiphi}.
As we can see in this figure, there is a peak close to threshold, which should not be mis-identified with a new resonance, but it is the reflection of the $X(4160)$ which in our fit has the mass at 4169~MeV.

\section{Summary}
The present work and the prediction, tied to the interpretation given for the $B^-\to K^-J/\psi\phi$ spectrum, should act as an incentive to measure this reaction and learn about properties of the $X(4140)$ and $X(4160)$. In addition, the $J/\psi \phi$ mass distribution of  the $e^+e^- \to \gamma J/\psi\phi$ reaction reported by BESIII~\cite{Ablikim:2017cbv} is also compatible with the existence of the $X(4140)$ state, appearing as a peak, and a strong cusp structure at the  $D_s^*\bar{D}_s^*$ threshold, resulting from the molecular nature of the $X(4160)$ state\cite{Wang:2018djr}.

\section*{Acknowledgements}
This work is partly supported by the National Natural
Science Foundation of China under Grant Nos.  11375024, 11522539, 11735003, 11505158, 11475015, and  11647601.  It is also supported by
the Youth Innovation Promotion Association CAS (No.
2016367),
and Academic Improvement Project of Zhengzhou University.
This work is also partly supported by the Spanish Ministerio de Economia
y Competitividad
and European FEDER funds under the contract number FIS2014-57026-REDT,
FIS2014-51948-C2-1-P, and FIS2014-51948-C2-2-P, and the Generalitat
Valenciana in the program Prometeo II-2014/068.

\bibliographystyle{ws-procs9x6}


\section{Introduction}\label{sec:intro}

When transmitting information across a network, latency can be caused either by \emph{propagation delays}, which are due to the transmission medium and are proportional to the distance the signal must travel, or by \emph{processing delays}, which are due to encoding, decoding, buffering, filtering, or other signal processing that must happen on either end of the transmission.

In scenarios where distances are relatively short, such as swarms of unmanned aerial vehicles (UAVs) communicating over a wireless network, propagation delays are negligible and it is reasonable to assume that latency is due entirely to processing delays. With UAVs, processing delays are often fixed and known, since they are dictated by the hardware capabilities of the UAVs and the bandwidth of the communication channel.

We consider the problem of controlling a set of $N$ \emph{dynamically decoupled} plants, which we refer to as \emph{agents}. A four-agent example is depicted in Fig.~\ref{fig:cartoon}.

\begin{figure}[ht]
	\centering
	\includegraphics{blk_cartoon2}
    \caption{Example of a dynamically decoupled system. The plants $\mathcal{G}_i$ have corresponding controllers $\mathcal{K}_i$ that share information along edges of a directed graph. Implied links are not shown, e.g., $\mathcal{K}_2\to\mathcal{K}_1$ and $\mathcal{K}_3\to\mathcal{K}_4$.}
    \vspace{-1mm}
	\label{fig:cartoon}
\end{figure}

We make standard linear-quadratic-Gaussian (LQG) assumptions. That is, (1) the $\mathcal{G}_i$ and $\mathcal{K}_i$ are assumed to be continuous-time linear time-invariant (LTI) systems, (2) the exogenous noise is assumed to be Gaussian and uncorrelated between agents, and (3) the objective function is quadratic in the state and inputs.
Let $x_i$ denote the state of the agent $i$, and aggregate the agents' states into a global state vector $x$. Proceed similarly for the inputs $u$, measurements $y$, and exogenous noise $w$. We can represent the global dynamics compactly as in Fig.~\ref{fig:centralized_LQG}.

\begin{figure}[ht]
	\centering
	\includegraphics{global_dynamics2}
	\caption{Classical four-block LQG problem with state-space representation of the plant. The goal is to design $\mathcal{K}$ to minimize the $\Htwo$ norm of the map from $w$ to $z$. Blue matrices a have block-diagonal structure due to the dynamically decoupled assumptions.}
	\label{fig:centralized_LQG}
\end{figure}

Due to our dynamically decoupled assumption, the aggregated plant matrices $A$, $B_1$, $B_2$, $C_2$, $D_{21}$ are block-diagonal. We allow the quadratic cost to couple the states and inputs of various agents, so the cost matrices $C_1$ and $D_{12}$ have no special structure.

We assume inter-controller communication happens instantaneously across the directed graph that represents the communication network (dashed lines in Fig.~\ref{fig:cartoon}), but there is a \emph{processing delay} of $\tau$ seconds for all incoming transmissions. We define the connectivity matrix $S \in \{0,1\}^{N\times N}$ as $S_{ij}=1$ if there is a directed path from controller $j$ to controller $i$ and $S_{ij}=0$ otherwise. The processing delays impose block structural constraints on $\mathcal{K}$. Specifically, $\mathcal{K}_{ij}$ is: zero if $S_{ij}=0$ (no communication), delayed by $\tau$ if $S_{ij}=1$ and $i\ne j$ (processing delay), and unconstrained if $S_{ij}=1$ and $i=j$ (no delay). For example, the connectivity and controller structures for the example of Fig.~\ref{fig:cartoon} are
\[
S\!=\!\squeezemat{1mm}\bmat{1&1&1&0\\ 1&1&1&0\\ 1&1&1&0\\ 1&1&1&1}\!,\,\,
\mathcal{K} \!\sim\! \left[\begin{array}{cccc}\mathcal{K}_{11}&e^{-s\tau}\mathcal{K}_{12}&e^{-s\tau}\mathcal{K}_{13}&0\\ e^{-s\tau}\mathcal{K}_{21}&\mathcal{K}_{22}&e^{-s\tau}\mathcal{K}_{23}&0\\ e^{-s\tau}\mathcal{K}_{31}&e^{-s\tau}\mathcal{K}_{32}&\mathcal{K}_{33}&0\\ e^{-s\tau}\mathcal{K}_{41}&e^{-s\tau}\mathcal{K}_{42}&e^{-s\tau}\mathcal{K}_{43}&\mathcal{K}_{44}\end{array}\right]\!.
\]
In this example, the controller $\mathcal{K}_3$ receives (delayed) measurements from $\mathcal{G}_1$ and $\mathcal{G}_2$, instantaneous measurement from $\mathcal{G}_3$, and no information from $\mathcal{G}_4$. In other words,
\[
u_3(t) = f_3\bigl( y_1(0:t-\tau), y_2(0:t-\tau), y_3(0:t) \bigr)
\]
for some function $f_3$. In the remainder of the paper, we derive the optimal structured $\mathcal{K}$ as a function of the connectivity $S$, delay $\tau$, and state-space matrices of the four-block plant (Fig.~\ref{fig:centralized_LQG}). We also provide a detailed and intuitive \emph{agent-level} description of the optimal $\mathcal{K}_i$.

\paragraph{Summary of relevant literature.}

The sparsity and delay constraints on our plant and controller are an example of a \textit{quadratically invariant} (QI) system~\cite{rotkowitz2005characterization,rotkowitz2010convexity}. Therefore, the problem of finding the optimal structured linear controller amounts to solving a convex (but infinite-dimensional) optimization problem.

When the structure does not contain delays, the exact optimal controller may be computed using vectorization~\cite{rotkowitz_vectorization,vamsi_elia}. In general, vectorization does not produce a minimal realization, nor does it elucidate structure, such as a controller-estimator separation or an interpretation of the signals communicated between subcontrollers.
Explicit solutions to special cases of output-feedback have been reported, such as: triangular~\cite{lessard2015optimal,tanaka2014optimal}, broadcast~\cite{lessard2012decentralized}, or dynamically decoupled~\cite{kim2015explicit,kashyap2019explicit} cases.

When delays are present, we distinguish the cases of discrete vs.\ continuous-time delay. In discrete time, the delay operator $z^{-1}$ is rational, so delays may be absorbed into the plant and the problem reduces to the non-delayed case~\cite{spdel}. In continuous time, this reduction is not possible because the delay operator $e^{-s\tau}$ is now irrational. One approach is to use a Pad\'{e} approximation for the delay~\cite{yan1996teleoperation} followed by vectorization. Alternatively, a Galerkin-style finite-dimensional approximation~\cite{scherer02,voulgaris_stabilization} can be used.

The present work seeks an explicit solution for the delayed case that provides structure and intuition without resorting to approximations or vectorization.
Some special cases have been addressed in the literature. For centralized problems with a fixed loop delay (\emph{dead time}), a loop-shifting technique involving FIR blocks can transform the problem into an equivalent LQG problem with a finite-dimensional LTI plant~\cite{mirkin2003every,mirkin2003extraction}. A similar idea was used in the discrete-time case to decompose the structure into dead time and FIR components, which can be optimized separately~\cite{lampHtwoDelay}.

The aforementioned loop-shifting approach was also extended to the \emph{adobe delay} case, where the feedback loop contains both a delayed and a non-delayed path~\cite{mirkin2011dead,mirkin2012h2}. This approach was used to obtain explicit optimal controllers for bilateral teleoperation problems, where the controllers communicate across a delayed channel~\cite{kristalny2012decentralized,cho2012h2}. The idea was also generalized to haptic interfaces that have two-way communication with a shared virtual environment~\cite{kristalny2013decentralized}. These are special cases of the general problem we will treat, described in Section~\ref{sec:intro}, where $S$ is the all-ones matrix.

\paragraph{Overview.}

Our work generalizes the agent-level solution for the non-delayed case reported in~\cite{kashyap2019explicit} by leveraging the loop-shifting decomposition for adobe delays from~\cite{mirkin2012h2}. The paper is organized as follows.
Section~\ref{sec:notation} covers our notation and assumptions.
Section~\ref{sec:result} presents our solution to the general problem as a function of the connectivity matrix $S$ and processing delay $\tau$. Our results are presented in continuous time with an infinite-horizon cost but can easily be generalized to discrete time and/or a finite-horizon cost. We also characterize the signals that should be transmitted between subcontrollers and show that each subcontroller has an intuitive observer-regulator structure.
Finally, Section~\ref{sec:diss} shows that we recover the expected limiting cases in the limits of zero or infinite processing delay.

\section{Preliminaries}\label{sec:notation}

\paragraph{Notation.}
State-space notation for transfer functions:
\begin{align*}
	\mathcal{G}(s) = \left[\begin{array}{c|c}%
		A&B\\ \hlinet C&D\end{array}\right]= D+C(sI-A)^{-1}B.
\end{align*}
We let $N$ denote the total number of agents and $[N] \defeq \{1,\dots,N\}$. The $i^\text{th}$ subsystem has state dimension $n_i$, input dimension $m_i$, and measurement dimension $p_i$. The global state dimension is $n \defeq n_1+\cdots+n_N$ and similarly for $m$ and $p$. The identity matrix of size $k$ is denoted $I_k$. We write $\blkdiag(\{X_i\})$ to denote the block-diagonal matrix formed by the blocks $\{X_1,\dots,X_n\}$ and $\diag(X)$ to denote the block-diagonal matrix formed by the diagonal blocks of $X$.
The zeros used throughout are matrix or vector zeros and their sizes are determined from context.
The symbol $\otimes$ denotes the Kronecker product.

For $i\in [N]$, we write $\underline{i} \subseteq [N]$ to denote the \textit{descendants} of node $i$, i.e., the set of nodes $j$ such that there is a directed path from $i$ to $j$. Likewise, $\bar{i} \subseteq [N]$ denotes the \textit{ancestors} of node $i$. Similarly, $\bar{\bar{i}}$ and $\underline{\underline{i}}$ denote the \textit{strict ancestors} and \textit{strict descendants}, respectively. We also write $s^c = [N]\setminus s$ to denote the complement.
For example, in the graph of Fig.~\ref{fig:cartoon}, we have $\underline{2} = \{1,2,3,4\}$, $\underline{\underline{2}} = \{1,3,4\}$, and $\bar{\bar{2}} = \{1,3\}$.
We also use this notation to index matrices.
For example, if $X$ is a $4\times 4$ block matrix associated with Fig.~\ref{fig:cartoon}, then
$X_{1\bar{\bar{2}}} = \bmat{X_{11} & X_{13}}$.

We will require the use of specific partitions of the identity matrix. We define $I_{n} \defeq \blkdiag(\{I_{n_i}\})$ and for each agent $i \in [N]$, we define $E_{n_i} \defeq (I_{n})_{:i}$ (the $i^\textup{th}$ block column of $I_{n}$). Similar to the descendant and ancestor definitions, $n_{\underline{i}}\defeq\sum_{k\in \underline{i}}n_k$ and $n_{\bar{i}}\defeq\sum_{k\in \bar{i}}n_k$. The dimensions of $E_{{n_{\bar{i}}}}$ and $E_{{n_{\underline{i}}}}$ are determined by the context of use. Finally, $1_{n}$ is the $n\times 1$ matrix of $1$'s.

\paragraph{Processing delay notation~\cite{mirkin2012h2}.}
The \textit{adobe delay} matrix $\Lambda_m^i \defeq \blkdiag(I_{m_i},e^{-s\tau}I_{m_{\underline{\underline{i}}}})$ leaves block $i$ unchanged and imposes a delay of $\tau$ on all strict descendants of $i$.
The completion operator $\pi_\tau\{\cdot\}$ acts on a state-space system delayed by $\tau$ and returns the FIR system with support on $[0,\tau]$ that completes it:
\[
\squeezemat{2pt}
\pi_\tau\Biggl\{\left[\begin{array}{c|c}A&B\\ \hlinet C&D\end{array}\right]\!e^{-s\tau}\Biggr\}
\defeq\left[\begin{array}{c|c}A&e^{-A\tau}\!B\\ \hlinet C&D\end{array}\right]-\left[\begin{array}{c|c}A&B\\ \hlinet C&D\end{array}\right]e^{-s\tau}\!.
\]
We also define the function $\Gamma : (\Omega,\Lambda_m^i) \mapsto (\tilde{\Omega},\Pi_u,\Pi_b)$, which maps a four-block plant $\Omega$ defined as
\begin{equation}\label{eq:fourblock}
	\Omega \defeq
	\bmat{ \mathcal{U} & \mathcal{V} \\ \W & \mathcal{G} } \defeq
	\left[\begin{array}{c|cc}
		A & B_1 & B_2\\ \hlinet
		C_1 & 0 & D_{12} \\
		C_2 & D_{21} & 0 \end{array}\right]
\end{equation}
and adobe delay matrix $\Lambda_m^i$ to a modified plant $\tilde{\Omega}$ and FIR systems $\Pi_u$ and $\Pi_b$. See Appendix~\ref{sec:appendix_gamma} for details.

\paragraph{Problem statement.}
Given a four-block plant~\eqref{eq:fourblock} and associated connectivity matrix $S$ as described in Section~\ref{sec:intro} (refer to Fig.~\ref{fig:centralized_LQG}), the plants are dynamically decoupled but the cost function may couple the states and inputs of the different agents. Iit follows that $\W$ and $\mathcal{G}$ are block diagonal but $\mathcal{U}$ and $\mathcal{V}$ need not. Let $\mathbb{S}_\tau$ be the set of causal structured LTI controllers with processing delay $\tau$, as described in Section~\ref{sec:intro}.
The problem is to find $\mathcal{K}$ to
\begin{equation}\label{opt}
\begin{aligned}
\underset{\mathcal{K}}{\minimize}
\qquad & \normm{\,\mathcal{U} + \mathcal{V}\mathcal{K}(I-\ \mathcal{G}\mathcal{K})^{-1}\W\,}_2^2 \\
\subject \qquad & \mathcal{K} \in \mathbb{S}_{\tau}
\text{ and $\mathcal{K}$ stabilizes $\Omega$.}
\end{aligned}
\end{equation}
We now describe and explain the technical assumptions we make on the four-block plant $\Omega$ and connectivity $S$.

\paragraph{Riccati assumptions.} Four matrices $(A,B,C,D)$ are said to satisfy the \emph{Riccati assumptions}~\cite{kim2015explicit,mirkin2012h2} if:
\begin{enumerate}[label*=R\arabic*., ref=R\arabic*]
	\item $C^\tp D=0$ and $D^\tp D>0$. \label{harf}
	\item $(A,B)$ is stabilizable.
	\item $\begin{bmatrix}A-j\omega I&B\\C&D\end{bmatrix}$ is full column rank for all $\omega \in \R$.
\end{enumerate}
If the Riccati assumptions hold, there is a unique stabilizing solution to the associated algebraic Riccati equation, which we write as $(X,F)=\ric(A,B,C,D)$. Thus $X \succeq 0$ satisfies
$
A^\tp X+XA+C^\tp C -XB(D^\tp D)^{-1}B^\tp X=0
$ with $A+BF$ Hurwitz, where $F\defeq -(D^\tp D)^{-1}B^\tp X$.

\begin{assumption}[\textbf{System assumptions}]
	\label{Ass:System}
	For the $N$ interacting agents, we will assume the following.
	\begin{enumerate}[label*=\arabic{thm}.\arabic*., ref=\arabic{thm}.\arabic*]
		\item $A_{ii}$ is Hurwitz for all $i\in[N]$, i.e., $A$ is Hurwitz. \label{Ass:System_hurwitz}
		\item The Riccati assumptions hold for $(A,B_2,C_1,D_{12})$ and for $(A_{ii}^\tp ,C_{2_{ii}}^\tp ,B_{1_{ii}}^\tp ,D_{{21}_{ii}}^\tp )$ for all $i\in[N]$. \label{Ass:System_riccati}
		\item $D_{12_{:\underline{i}}}^\tp D_{12_{:\underline{i}}}=I$ and $D_{21_{ii}}D_{21_{ii}}^\tp=I$ for all $i\in[N]$. \label{Ass:System_D1221}
	\end{enumerate}
\end{assumption}
\noindent Assumption~\ref{Ass:System_hurwitz} is an assumption of nominal stability, carried over from Kim et~al.~\cite{kim2015explicit}. Assumption~\ref{Ass:System_riccati} ensures the necessary condition that the centralized LQR problem and the individual agents' estimation problems are non-degenerate. Assumption~\ref{Ass:System_D1221} simplifies the exposition of the delayed problem~\cite{mirkin2011dead,mirkin2012h2}, though the results still hold for the general case~\cite[Rem.~3.2]{mirkin2011dead}.

\section{Main Result}\label{sec:result}

Our main result is an agent-level description of the optimal controller that solves~\eqref{opt}. That is, we provide explicit state-space formulas for each $\mathcal{K}_i$ and describe which signals should be transmitted and received between agents.

\begin{thm}
	\label{thm:3}
	Consider a general instance of the structured optimal control problem described in Section~\ref{sec:intro} and suppose Assumption~\ref{Ass:System} holds. Consider the four-block sub-plant for agent $i$ and its descendants:
	\begin{align*}
			{\Omega}_i \defeq
			\left[\begin{array}{c c}
				{\mathcal{U}}_{:i} & {\mathcal{V}}_{:\underline{i}}\\
				{\W}_{ii} & {\mathcal{G}}_{i\underline{i}}\end{array}\right] \defeq
			\left[\begin{array}{c|cc}
				A_{\underline{ii}} & B_{1_{\underline{i}i}} & B_{2_{\underline{ii}}} \\[2pt] \hlinet
				{C}_{1_{:\underline{i}}} & 0 & D_{12_{:\underline{i}}} \\
				C_{2_{i\underline{i}}} & D_{21_{ii}} & 0\end{array}\right].
	\end{align*}
	Now apply
	$(\tilde{\Omega}_i,\Pi_{u_i},\Pi_{b_i})=\Gamma(\Omega_i,\Lambda_m^i)$. Define the estimation gain $\tilde{F}^i$ and control gain $L^i$ by
	\begin{subequations}\label{eq:Riccati}
		\begin{align}
		(\tilde{X}^i,\tilde{F}^i) &\defeq \ric\bigl(A_{\underline{ii}},\tilde{B}_{2_{\underline{ii}}},\tilde{C}_{1_{:\underline{i}}},D_{12_{:\underline{i}}}\bigr),\\
		(Y^i,{L^i}^\tp) &\defeq
		\ric(A_{ii}^\tp,C_{2_{ii}}^\tp,B_{1_{ii}}^\tp,D_{{21}_{ii}}^\tp).
		\end{align}
    \end{subequations}
	A realization of the optimal controller $\mathcal{K}_i$ for agent $i$ that solves~\eqref{opt} is given by the state-space equations
\begin{subequations}\label{eq:state_agent_Ktau_reduced}
    \begin{align}
        \dot \eta_{i,\underline{i}} &= (A_{\underline{ii}}+(LC)^i)\eta_{i,\underline{i}} + (LC)^i \sum_{k\in \bar{\bar i}}\eta_{k,\underline{i}}(t-\tau) \notag\\
        &  \qquad\qquad + \tilde{B}_{2_{\underline{ii}}} \tilde{\nu}_{i,\underline{i}}- E_{n_{\underline{i}}}^{\tp} E_{n_i}L^i (y_i+\Pi_{b_i} \nu_{i,\underline{i}}),\\
        \tilde{\nu}_{i,\underline{i}} &= \tilde{F}^i \eta_{i,\underline{i}}, \\
        u_i &= E_{m_i}^\tp E_{m_{\underline{i}}}\biggl( \nu_{i,\underline{i}}+ \sum_{k\in\bar{\bar{i}}}\tilde{\nu}_{k,\underline{i}}(t-\tau) \biggr),
    \end{align}
\end{subequations}
where $(LC)^i\defeq E_{n_{\underline{i}}}^{\tp} E_{n_i} L^i C_{2_{ii}} E_{n_{i}}^{\tp} E_{n_{\underline{i}}}$ and $\nu_{i,\underline{i}} \defeq \Pi_{u_i} \tilde{\nu}_{i,\underline{i}}$.
\end{thm}

\begin{proof}
	See Appendix \ref{App:A}.
\end{proof}

We use a slight abuse of notation in~\eqref{eq:state_agent_Ktau_reduced}. These are time-domain equations but $\Pi_{b_i}$ and $\Pi_{u_i}$ are FIR transfer matrices. Products such as $\Pi_{b_i} \nu_{i,\underline{i}}$ are to be interpreted as the time-domain output of $\Pi_{b_i}(s)$ with input $\nu_{i,\underline{i}}(t)$.

Fig.~\ref{fig:Ki} shows the signals received and transmitted by each subcontroller. Fig.~\ref{fig:local} shows a more detailed block diagram. The optimal controller is an interconnection of state-space systems, FIR blocks, and delays.

\section{Limiting cases}\label{sec:diss}

\begin{figure}[ht]
	\centering
	\includegraphics{blkdiag_mini2}
	\caption{Block diagram representation of a local plant $\mathcal{G}_i$ and associated controller $\mathcal{K}_i$ (dashed red box). Controller $i$ receives local measurements $y_i$ and information $(\eta,\tilde\nu)$ from its strict ancestors $\bar{\bar{i}}$ (which incurs a processing delay $e^{-s\tau}$). The controller then computes  and transmits information to its strict descendants $\underline{\underline{i}}$.}
	\label{fig:Ki}
\end{figure}

\begin{figure*}[htb]
	\centering
	\includegraphics{blkdiag_local1}
	\caption{Expanded diagram of the optimal controller shown in Fig.~\ref{fig:Ki} that reveals the observer-regulator structure of $\mathcal{K}_i$ from Theorem~\ref{thm:3}. The green blocks represent blocks that depend on the processing delay $\tau$.
		%
		Agent $i$ updates its local innovation $\eta_{i,\underline i}$ and computes its local input $u_i$ using the local measurement $y_i$ and the delayed innovations $\eta_{k,\underline k}$ and partial inputs $\tilde v_{k,\underline k}$ from  strict ancestors $k\in \bar{\bar{i}}$.
		%
		For simplicity, we used the following notation for the observer dynamics: $\mathcal{T}_{\underline{ii}} \defeq  (sI-A_{\underline{ii}}-(LC)^i)^{-1}$, where $(LC)^i\defeq E_{n_{\underline{i}}}^{\tp} E_{n_i} L^i C_{2_{ii}} E_{n_{i}}^{\tp} E_{n_{\underline{i}}}$.}
	\label{fig:local}
\end{figure*}

We now study the behavior of our optimal controller as the processing delay $\tau$ varies. Setting $\tau=0$, nodes belonging to the same connected component can be treated as a single node due to the instantaneous communication assumption. For example, if $\tau=0$ in Fig.~\ref{fig:cartoon}, nodes $1$, $2$, and $3$ become a single node, and we are left with the equivalent two-node graph $S = \left(\begin{smallmatrix}1 & 0 \\ 1 & 1 \end{smallmatrix}\right)$. This special case was solved in~\cite{kashyap2019explicit} and we recover this result by setting $\tau=0$ in Theorem~\ref{thm:3}. At the other extreme, the case $\tau\to\infty$ leads to a controllers that do not communicate at all. These results are formalized in Corollary~\ref{corr:1} below.

\begin{cor}\label{corr:1}
Consider the setting of Theorem~\ref{thm:3}.
\begin{enumerate}[label*=\arabic*), ref=\arabic*]
    \item If there is no processing delay ($\tau=0$), The controller $\mathcal{K}_{0}$ that solves~\eqref{opt} is
    \begin{subequations}
		\begin{align*}
            \hspace{-2mm}\dot \eta_{i,\underline{i}} &= A_{\underline{ii}} \eta_{i,\underline{i}} +
            (LC)^i\sum_{k \in \bar{i}}\eta_{k,\underline{i}} + {B}_{2_{\underline{ii}}} \nu_{i,\underline{i}} - E_{n_{\underline{i}}}^{\tp} E_{n_i}L^iy_i\\
            \hspace{-1mm}u_i &= E_{m_i}^\tp\sum_{j\in\bar{i}} F_j \eta_j,
        \end{align*}
    \end{subequations}
    where $\nu_{i,\underline{i}} \defeq {F}^i \eta_{i,\underline{i}}$ and $F_j\defeq E_{m_{\underline{j}}}F^{j}E_{n_{\underline{j}}}^\tp$ and the remaining parameters are defined in Theorem~\ref{thm:3}.

    \item If there is infinite processing delay ($\tau\to\infty$), the controller $\mathcal{K}_{\infty}$ that  solves~\eqref{opt} is
	\begin{subequations}\label{eq:Kinfty}
		\begin{align}
			\dot \eta_{i,i} &= A_{ii} \eta_{i,i} + L^iC_{2_{ii}}\eta_{i,i}+ B_{2_{ii}} u_i-L^iy_i\\
			u_i &= F_\infty^i\eta_{i,i},
		\end{align}
	\end{subequations}
	where
    $(X_\infty^i,F_\infty^i) \defeq		\ric(A_{ii},B_{2_{ii}},C_{1_{:{i}}},D_{12_{:{i}}})$.
\end{enumerate}
\end{cor}

\begin{proof}
    When $\tau \to 0$, we have $\Pi_{u_i} = I$ and $\Pi_{b_i} =0$. Moreover, $\tilde\Omega_i = \Omega_i$ and $\tilde F^i = F^i$.
    When $\tau \to \infty$, we also have $\Pi_{u_i} \to I$ and $\Pi_{b_i} \to 0$ but this time, $\tilde B_{2_{\underline{ii}}} \to \diag(B_{2_{\underline{ii}}})$. The corresponding control gain $\tilde{F}^i$ reduces to an augmented nominal gain $E_{m_i}F_{\infty}^{i}E_{n_i}^\tp$ after elimination of the uncontrollable ($\underline{\underline{i}}$) and unobservable ($\underline{i}^c$) modes to obtain~(\ref{eq:Kinfty}). Letting $\tau\to\infty$ is equivalent to solving $N$ separate LQG problems (the global cost matrices are block-diagonal).
\end{proof}

The effect of changing the processing delay $\tau$ is illustrated in the block diagram of Fig.~\ref{fig:local}. When $\tau$ varies, only the green blocks are affected. The remaining parts of the controller can be designed without knowing $\tau$.

\section{Conclusion}
\label{sec:conclufuture}
In this work, we considered an optimal control problem where dynamically decoupled agents communicate over a network and incur processing delays to receive transmissions from neighboring controllers. We described the structure of the optimal control strategies at the level of individual agents (Fig.~\ref{fig:local}), which shows which signals should be transmitted and which parts of the controller depend on the processing delay.