XiangLi VivekGovindan RohitPaturi SundararajanSrinivasan

# Speakers Unembedded: Embedding-free Approach to Long-form Neural Diarization

###### Abstract

End-to-end neural diarization (EEND) models offer significant improvements
over traditional embedding-based Speaker Diarization (SD) approaches but falls
short on generalizing to long-form audio with large number of speakers. EEND-
vector-clustering method mitigates this by combining local EEND with global
clustering of speaker embeddings from local windows, but this requires an
additional speaker embedding framework alongside the EEND module. In this
paper, we propose a novel framework applying EEND both locally and globally
for long-form audio without separate speaker embeddings. This approach
achieves significant relative DER reduction of 13% and 10% over the
conventional 1-pass EEND on Callhome American English and RT03-CTS datasets
respectively and marginal improvements over EEND-vector-clustering without the
need for additional speaker embeddings. Furthermore, we discuss the
computational complexity of our proposed framework and explore strategies for
reducing processing times.

###### keywords:

speaker diarization, end-to-end diarization, spectral clustering

##  1 Introduction

Speaker diarization addresses the ``who spoken when" problem by partitioning
an audio stream containing multiple speakers into homogeneous segments
associated with each speaker. Conventional diarization systems [1, 2, 3, 4, 5,
6] typically consist of a cascade of several separate modules: voice activity
detection to detect the speech frames, speaker embedding extraction to
transform the speech segments into discriminative representations, and
clustering to group speech regions by speaker identity. While effective for
long-form audio with an arbitrary number of speakers, these cascaded multi-
module approaches face challenges in handling overlapping speech and can
suffer from error propagation across the modules.

To overcome the limitations of cascaded approaches, end-to-end neural
diarization (EEND) was proposed in [7] which formulates speaker diarization as
a frame-wise multi-label classification task with permutation invariant
training [8]. EEND can naturally handle overlapping speech by allowing
multiple speakers to be active simultaneously and is also fully supervised
compared to the unsupervised clustering component of the cascaded approach.
However, despite its theoretical promise, EEND and its variants like EEND-SA
[9], EEND-EDA [10], etc have struggled to generalize to larger numbers of
speakers and arbitrarily long conversations.

In order to apply EEND models to longer audios and larger number of speakers,
recent works [11, 12] have proposed hybrid frameworks that integrate EEND with
conventional clustering-based approaches. These methods leverage the strong
diarization capability of EEND for speaker labeling over short local windows
while performing global clustering on speaker embeddings computed across the
local windows. This hybrid approach can handle both overlapping speech locally
and long conversations with an arbitrary number of speakers globally. Most of
the recent EEND improvements have focused on integrating additional embedding
[11, 12, 13, 14, 15, 16, 17] or attractor modules [18, 19, 20, 21, 22, 23,
24], which requires specialized model architectures, loss functions and data
requirements. Moreover, in some real-world scenarios, creating and storing
speaker embeddings may need to be avoided where possible due to privacy
considerations [25].

In this paper, we propose a novel embedding-free approach that doesn't require
any speaker embeddings and can still leverage the benefits of EEND and scale
it to long-form audios with arbitrary number of speakers. We achieve this by
utilizing a vanilla EEND model for both local diarization within the short
local windows as well as global diarization across local windows, hence named
local-global EEND. The proposed method consists of three steps: local EEND,
global EEND, and clustering. In the local step, long audio is split into
fixed-size windows, and EEND performs diarization within each window. The
global step solves the inter-window label permutation by re-applying EEND to
chunks formed by pairing speaker chunks across local windows. This generates
pairwise speaker scores which are used to build an affinity matrix for the
final clustering and global speaker labeling, without requiring any speaker
embeddings.

The rest of the paper details the local-global EEND approach, experimental
setup, results compared to the baselines, and a discussion on potential
computation improvements for the global step.

##  2 Local-global EEND

Figure 1: Local-global EEND framework. This assumes 3 local windows with
2-speaker local EEND, i.e $W$=3, $S_{local}$=2 resulting in $C$=12 pairwise-
speaker chunks for global EEND.

###  2.1 Local EEND

Figure 1 shows the schematic diagram of the proposed embedding-free approach
which can be divided into local EEND, global EEND and clustering steps.

The input audio is first split into $W$ windows with a fixed window length. In
each window $i$, frame-level acoustic features are extracted, denoted as
$\textbf{X}_{i}=\\{\textbf{x}_{i,t}\\}_{t=1}^{T},\textbf{x}_{i,t}\in\mathbb{R}^{F}$
where $t$ is the frame index, $T$ is the total number of frames in a window
and $F$ is the feature dimension of Mel-filterbank features in this work.
Speaker label $\textbf{y}_{i,t}=\\{y_{i,t,s}\\}_{s=1}^{S_{local}}$ denotes
speech activities for $S_{local}$ speakers at frame $t$ within window $i$ and
is defined as

$\displaystyle y_{i,t,s}=\begin{cases}0&\text{(Speaker $s$ is inactive at
$t$)}\\\ 1&\text{(Speaker $s$ is active at $t$)}\\\ \end{cases}$ (1)

The local EEND estimates frame-wise posteriors $P(y_{i,t,s}|{\textbf{X}_{i}})$
in each window using a vanilla EEND model. These posteriors are binarized
using a threshold $Th_{local}$ and median filtered [9] to obtain the local
speaker labels $\textbf{y}_{i,t}$.

###  2.2 Global EEND

In order to perform global SD, the global EEND step computes the speaker
similarities across the local windows using the same EEND model. In order to
compute these, the overlapping speaker frames within each local window are
first filtered out and the remaining frames of each speaker in a window are
paired with the frames of speakers in subsequent windows, resulting in new
chunks $\\{\hat{\textbf{X}}_{i}\\}_{i=1}^{C}$

$\hat{\textbf{X}}_{i}=concat(\textbf{x}_{j,t}|t=\\{m_{j,s}\\}_{s=1}^{M},\textbf{x}_{k,t}|t=\\{n_{k,s}\\}_{s=1}^{N})$
(1)

where $\textbf{x}_{j,t}$ and $\textbf{x}_{k,t}$ represent frame-level acoustic
features of window $j$ and window $k$ $(j\not=k)$, respectively.
$\\{m_{j,s}\\}_{s=1}^{M}$ represents the $M$ frame indices corresponding to
speaker $m$ in window $j$ and $\\{n_{k,s}\\}_{s=1}^{N}$ represents the $N$
frame indices corresponding to speaker $n$ in window $k$. $C$ is the total
number of pairwise-speaker chunks processed by global EEND, where

$C\leq W\times(W-1)/2\times{S_{local}}^{2}$ (2)

In the case where a speaker has limited or no non-overlapping frames, we
leverage the overlapping frames similarly to the EEND-vector clustering
approach.

EEND is applied to $\hat{\textbf{X}}_{i}$ to generate inter-window frame-level
speaker posteriors

$[\textbf{z}_{1},...,\textbf{z}_{M},\textbf{z}_{M+1},...,\textbf{z}_{M+N}]=EEND(\hat{\textbf{X}}_{i})$
(3)

where
$\\{\textbf{z}_{t}\\}_{t=1}^{M+N},\textbf{z}_{t}\in\mathbb{R}^{S_{local}}$ are
the inter-window frame-level posteriors. $[\textbf{z}_{1},...,\textbf{z}_{M}]$
and $[\textbf{z}_{M+1},...,\textbf{z}_{M+N}]$ are the posteriors corresponding
to the $M$ frames of speaker $m$ and $N$ frames of speaker $n$ respectively.
This process is repeated on every speaker pair across local windows as shown
in Figure 1.

###  2.3 Embedding-free clustering

The frame-wise posteriors $\\{\textbf{z}_{t}\\}_{t=1}^{M+N}$ are aggregated on
frames belonging to the same speaker, resulting in speaker-level posteriors
$\overline{\textbf{z}}_{m}$ and $\overline{\textbf{z}}_{n}$. Pairwise-speaker
similarity $S_{mn}$ is then calculated as

$\overline{\textbf{z}}_{m}=mean([\textbf{z}_{1},...,\textbf{z}_{M}])$ (4)
$\overline{\textbf{z}}_{n}=mean([\textbf{z}_{M+1},...,\textbf{z}_{N}])$ (5)
$S_{mn}=cosine\\_simalirity(\overline{\textbf{z}}_{m},\overline{\textbf{z}}_{n})$
(6)

Each $S_{mn}$ is an entry of the affinity matrix
$\textbf{S}\in\mathbb{R}^{S_{Global}\times S_{Global}}$ which will be used for
the final clustering, where $S_{Global}$ is the sum of number of speakers
detected in each local window with the local EEND, upper bound of which will
be $W\times S_{local}$.

In order to enhance the clustering performance as well as to save on
additional computations, we incorporate cannot-link constraints among
different speakers identified within the same local window obtained in the
local EEND step. This constraint is enforced by assigning a speaker similarity
of 0 between local speaker pairs. Spectral clustering is then employed to
group the speaker frames into $D$ speaker sets using the max eigengap
heuristic similar to [12, 6].

Table 1: Effect of sequence concatenation on CALLHOME2. $EEND_{vanilla}$ is the EEND model adapted on unmodified utterances. $EEND_{concat}$ is the EEND model adapted on speaker concatenated sequences described in 3.3. $EEND_{vanilla}+EEND_{concat}$ uses $EEND_{vanilla}$ for the local step and $EEND_{concat}$ for the global step of our local-global approach. System | Model | # of speakers in a session
---|---|---
|  | 2 | 3 | 4 | 5 | 6 | all
1-pass EEND | $EEND_{vanilla}$ | 7.53 | 14.91 | - | - | - | -
$EEND_{concat}$ | 7.36 | 17.74 | - | - | - | -
Local-global EEND | $EEND_{vanilla}$ | 7.99 | 12.21 | 16.39 | 17.10 | 26.12 | 12.48
$EEND_{concat}$ | 7.29 | 11.85 | 17.83 | 15.76 | 22.38 | 12.16
$EEND_{vanilla}+EEND_{concat}$ | 7.66 | 11.67 | 16.03 | 17.56 | 23.71 | 12.45

##  3 Experiments

In this section, we go over the datasets used, model architecture, settings
and techniques followed for efficient inference.

###  3.1 Data and metrics

For training the EEND model, we generate simulated mixtures by mixing
Switchboard-2 (Phase I & II & III), Switchboard Cellular (Part 1 & 2), and the
NIST Speaker Recognition Evaluation (2004 & 2005 & 2006 & 2008) with MUSAN
corpus [26], following the data generation procedure in [7]. Mixtures with up
to 3 speakers were created, with $\beta=[2,2,9]$ for mixture with 1, 2 and 3
speakers, respectively.

For model adaptation and evaluation, real telephone conversation dataset
CALLHOME [27], i.e., NIST SRE2000 (LDC2001S97, Disk-8) is used. It is widely
used as the benchmark for existing EEND-based approaches. The CALLHOME dataset
contains 500 sessions, each with 2 to 6 speakers. There are mostly two
dominant speakers in each conversation. We split the data into two subsets
according to [10] for adaptation (CALLHOME1) and evaluation (CALLHOME2).

As local-global EEND framework is designed for dealing with long
conversations, to showcase the effectiveness of this framework, evaluations on
other benchmarks with longer audios are reported as well, such as CALLHOME
American English (CHAE) [28] and RT03-CTS [29] which have an average duration
of 30 and 10 minutes respectively. We use the official eval splits for
evaluation on these datasets.

For evaluation metrics, we use the standard Diarization Error Rate (DER) [30]
with a collar tolerance of 250ms and included the overlapping speech segments
while scoring.

###  3.2 EEND model settings

To ensure a fair comparison with the existing hybrid baseline, we adopt the
front-end configuration from EEND-vector-clustering [12]. This involves the
extraction of 23-dimensional log-Mel-filterbank features, utilizing a frame
length of 25ms and a frame shift of 10ms. The window size $T$ is set at 300
(=30s) for both training and adaptation. The EEND architecture consists of 6
stacked self-attention-based Transformer layers, featuring eight attention
heads and a hidden size of 256. This aligns with the configuration employed in
[12]. In each window, the EEND model estimates the posteriors for up to 3
speakers.

During both training and adaptation, we employ the Adam optimizer [31]
alongside the Noam scheduler [32], incorporating 150,000 warm-up steps for
training. For adaptation, a fixed learning rate of $1\times 10^{-5}$ is
utilized. Both training and adaptation phases span 100 epochs.

###  3.3 Sequence concatenation during adaptation

In the global EEND step, we generate the chunk-level input by concatenating
frame-level acoustic features between every pair of speakers across local
windows. This process results in a new pairwise-speaker sequence that has not
been encountered in either the training or adaptation data. To enhance EEND's
generalization to this new input format, we incorporate this data generation
procedure during adaptation. First, the frame-level acoustic features from the
same speaker in each utterance are aggregated into several blocks after
discarding the overlapping speech frames. Subsequently, every two blocks are
concatenated to generate a new input utterance. We reformat the data using
this technique for half of the samples in each batch.

###  3.4 Efficiency improvement for inference

Real Time Factor (RTF) is a criteria used to measure the efficiency of SD
systems. It is calculated by dividing the time taken by the SD system by the
total duration of the spoken audio.

In the global step during inference, each speaker within a local window is
paired with every speaker in subsequent local windows, resulting in a
computational complexity of $O(n^{2})$. As illustrated in Figure 1, if there
are 3 local windows, each containing 2 speakers, there will be 12 inference
calls in the global EEND step. To enhance GPU efficiency and reduce RTF, we
propose batching multiple inference requests together. Additionally, we
explore different number of random frames ($N=128,64,32,16$) to minimize the
number of frames required for each speaker during global EEND inference,
thereby reducing computational load.

##  4 Results

In this section, we present the outcomes of our experiments, beginning with an
evaluation of the impact of sequence concatenation in the adaptation process.
Subsequently, we compare the proposed approach with existing baselines,
utilizing both ora- cle and estimated speaker counts during clustering. Our
analysis extends to additional benchmarks, such as CALLHOME American English
(CHAE) and RTCTS. Finally, we delve into an examination of the efficiency
improvements in global EEND inference.

###  4.1 Effect of sequence concatenation

We explored two types of input data formats during adaptation, resulting in
$EEND_{vanilla}$ and $EEND_{concat}$. The former denotes the model adapted
with original input sequences typically used in EEND model adaptation, while
the latter is adapted with a combination of original sequences and
concatenated sequences, as described in Section 3.3. Table 1 presents a
comparison between the local-global EEND and 1-pass EEND, utilizing the
different adaptation techniques. We used oracle speaker information during
global clustering and selected a binarization threshold between 0.3 to 0.7
that produced the best DER on the validation set for both 1-pass and local
EEND. We only evaluated 1-pass EEND on 2,3 speaker sessions since it was
trained to only detect a maximum of 3 speakers.

Local-global EEND outperforms 1-pass EEND on $EEND_{vanilla}$ by 18% in the
3-speaker session but shows a marginal performance degradation in the
2-speaker session whereas $EEND_{concat}$ outperforms the best 1-pass EEND by
3% and 21% in 2-speaker and 3-speaker sessions, respectively.

For 1-pass EEND, $EEND_{concat}$ only marginally benefits the 2-speaker
session but not the 3-speaker one, as expected due to the data mismatch
between adaptation and evaluation, where concatenation during adaptation only
occurs on pairwise speakers. When evaluating more speakers for local-global
EEND, $EEND_{concat}$ consistently performs better then $EEND_{vanilla}$,
except for the 4-speaker session.

In order to exactly match the local and global input conditions, we also
attempted to apply $EEND_{vanilla}$ in local EEND and $EEND_{concat}$ in
global EEND. This further improved performance in 3 and 4 speaker sessions but
$EEND_{concat}$ for both steps achieves the best overall DER across all
speaker sessions.

Table 2: Comparison with other baselines with oracle number of speakers on
CALLHOME2. The best scores are bolded.

System | # of speakers in a session |
---|---|---
| 2 | 3 | 4 | 5 | 6 | all
x-vector-clustering [10] | 8.93 | 19.01 | 24.48 | 32.14 | 34.95 | 18.98
EDA-EEND [10] | 8.35 | 13.20 | 21.71 | 33.00 | 41.07 | 15.43
EEND-vector-clust. (T=30s) [12] | 8.08 | 11.27 | 15.01 | 23.14 | 26.56 | 12.22
Local-global EEND | 7.29 | 11.85 | 17.83 | 15.76 | 22.38 | 12.16

Table 3: Comparison with other baselines with estimated number of speakers on
CALLHOME2. The best scores are bolded and the second best are underlined.

System | # of speakers in a session |
---|---|---
| 2 | 3 | 4 | 5 | 6 | all
x-vector-clustering [10] | 15.54 | 18.01 | 22.68 | 31.40 | 34.27 | 19.43
EDA-EEND [10] | 8.50 | 13.24 | 21.46 | 33.16 | 40.29 | 15.29
EEND-vector-clust. (T=30s) [12] | 7.96 | 11.93 | 16.38 | 21.21 | 23.10 | 12.49
EEND-EDA-local-global [18] | 7.11 | 11.88 | 14.37 | 25.95 | 21.95 | 11.84
Local-global EEND | 7.51 | 12.20 | 17.88 | 16.01 | 22.35 | 12.20

###  4.2 Comparison with other baselines

The performance of local-global EEND is compared with other baselines in Table
2 and Table 3, with oracle and estimated numbers of speakers respectively.
Results for EEND-vector-clustering with a window size of 30s are extracted
from their paper, and the outcomes for local-global EEND with $EEND_{concat}$
are reported.

In the case of oracle numbers of speakers, local-global EEND outperforms EEND-
vector-clustering [12] in all sessions except for 3 and 4-speaker sessions,
with only a marginal degradation in the 3-speaker session. Particularly
noteworthy is the substantial improvement of 32% and 15.7% in 5- and 6-speaker
sessions, respectively. The decline in the 4-speaker session may be attributed
to the sub-optimal window size, a phenomenon observed in the EEND-vector-
clustering paper as well, suggesting that different window sizes may result in
significant performance variations.

Table 3 presents results with estimated numbers of speakers, introducing
another strong baseline, EEND-EDA-local-global [18]. Compared to EEND-vector-
clustering, local-global EEND achieves superior performance across nearly all
sessions, with similar exceptions in the 4-speaker session. In comparison to
EEND-EDA-local-global, local-global EEND exhibits slightly inferior
performance in the $\\{2,3,6\\}$-speaker sessions but significantly
outperforms in the 5-speaker session. This discrepancy could stem from
differences in training data volume where the local-global EEND is trained
with only up to 3-speaker mixtures, whereas EEND-EDA-local-global is trained
on a larger training data including 100,000 additional 4-speaker mixtures.
Consequently, EEND-EDA-local-global achieves the best DER in the 4-speaker
session, aligning with this matched training scenario.

Table 4: DER (%) on long-form audio datasets ($Th_{local}$=0.5).

System | Dataset
---|---
| CHAE test | CHAE 109 | RTCTS test
1-pass EEND | 5.95 | 7.25 | 5.69
Local-global EEND | 5.20 | 7.02 | 5.14

###  4.3 Performance on other benchmarks

To provide a comprehensive evaluation of the local-global EEND system, we
extend our analysis to other well-established diarization benchmarks, namely
CALLHOME American English (CHAE) and RTCTS. These datasets feature longer-
duration audios, offering insights into the efficacy of local-global
diarization systems. Binarization is carried out using a fixed threshold
($Th_{local}$=0.5) to obtain diarization results, and clustering is performed
using the estimated number of speakers.

As depicted in Table 4, local-global EEND demonstrates superior performance
compared to 1-pass EEND across various datasets. Specifically, it outperforms
1-pass EEND by 12.7%, 9.7% and 3%, on the CHAE test set, RTCTS test set and
CHAE 109, respectively.

###  4.4 Inference efficiency improvements

Figure 2 shows the results on CHAE test set with different strategies to
improve inference efficiency. All the experiments are performed on NVIDIA A10G
GPU on AWS cloud (G5-2xLarge). Moving from sequential inference to batching
(batching 500 chunks), the RTF is reduced by 50%. Further RTF reduction is
gained from reducing computation by selection a subset of frames (N=128, 64,
32, 16) randomly for each speaker in global EEND. A subset of 64 frames can
produce a desirable RTF reduction by nearly 70% with no impact on DER.
Regarding the computational cost versus the input audio length, the local-
global EEND produces the RTF of
$\\{7.3e^{-3},1.5e^{-2},2.2e^{-2},5.0e^{-2}\\}$ for the audio length of
$\\{5,10,15,30\\}mins$, respectively.

Figure 2: RTF vs DER with different strategies on efficiency improvement,
including batching the inferences and minimizing the number of frames required
for each speaker. N indicates a subset of N random frames.

##  5 Conclusion

This paper introduces a novel embedding-free diarization methodology that
employs EEND in both local and global steps. The global clustering is
accomplished without the need for speaker embeddings, utilizing EEND on
concatenated pairwise speaker features across local windows to derive the
pairwise speaker similarities. This approach achieves significant relative DER
reduction of 13% and 10% over the conventional 1-pass EEND on CHAE and
RT03-CTS datasets respectively and even offers a marginal 3% relative DER
reduction over EEND-vector-clustering without the need for additional speaker
embeddings or loss functions. The paper also includes a discussion on the
computational complexity of the global EEND step and explores strategies for
reducing the processing times. By batching multiple chunk-level inferences and
minimizing the number of frames required for each speaker, the RTF can be
reduced by nearly 70% without the impact on diarization performance.

## References

  * [1] X. Anguera, S. Bozonnet, N. Evans, C. Fredouille, G. Friedland, and O. Vinyals, ``Speaker diarization: A review of recent research,'' _IEEE Transactions on audio, speech, and language processing_ , vol. 20, no. 2, pp. 356–370, 2012.
  * [2] G. Sell and D. Garcia-Romero, ``Speaker diarization with plda i-vector scoring and unsupervised calibration,'' in _2014 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2014, pp. 413–417.
  * [3] D. Snyder, P. Ghahremani, D. Povey, D. Garcia-Romero, Y. Carmiel, and S. Khudanpur, ``Deep neural network-based speaker embeddings for end-to-end speaker verification,'' in _2016 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2016, pp. 165–170.
  * [4] D. Garcia-Romero, D. Snyder, G. Sell, D. Povey, and A. McCree, ``Speaker diarization using deep neural network embeddings,'' in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 4930–4934.
  * [5] G. Sell, D. Snyder, A. McCree, D. Garcia-Romero, J. Villalba, M. Maciejewski, V. Manohar, N. Dehak, D. Povey, S. Watanabe _et al._ , ``Diarization is hard: Some experiences and lessons learned for the jhu team in the inaugural dihard challenge.'' in _Interspeech_ , 2018, pp. 2808–2812.
  * [6] Q. Wang, C. Downey, L. Wan, P. A. Mansfield, and I. L. Moreno, ``Speaker diarization with lstm,'' in _2018 IEEE International conference on acoustics, speech and signal processing (ICASSP)_. IEEE, 2018, pp. 5239–5243.
  * [7] Y. Fujita, N. Kanda, S. Horiguchi, K. Nagamatsu, and S. Watanabe, ``End-to-end neural speaker diarization with permutation-free objectives,'' in _Interspeech_ , 2019, pp. 4300–4304.
  * [8] D. Yu, M. Kolbæk, Z.-H. Tan, and J. Jensen, ``Permutation invariant training of deep models for speaker-independent multi-talker speech separation,'' in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 241–245.
  * [9] Y. Fujita, N. Kanda, S. Horiguchi, Y. Xue, K. Nagamatsu, and S. Watanabe, ``End-to-end neural speaker diarization with self-attention,'' in _2019 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2019, pp. 296–303.
  * [10] S. Horiguchi, Y. Fujita, S. Watanabe, Y. Xue, and K. Nagamatsu, ``End-to-end speaker diarization for an unknown number of speakers with encoder-decoder based attractors,'' in _Interspeech_ , 2020, pp. 269–273.
  * [11] K. Kinoshita, M. Delcroix, and N. Tawara, ``Integrating end-to-end neural and clustering-based diarization: Getting the best of both worlds,'' in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 7198–7202.
  * [12] ——, ``Advances in integration of end-to-end neural and clustering-based diarization for real conversational speech,'' in _Interspeech_ , 2021, pp. 3565–3569.
  * [13] N. Zeghidour, O. Teboul, and D. Grangier, ``Dive: End-to-end speech diarization via iterative speaker embedding,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2021, pp. 702–709.
  * [14] Y. Yu, D. Park, and H. K. Kim, ``Auxiliary loss of transformer with residual connection for end-to-end speaker diarization,'' in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 8377–8381.
  * [15] I. Fung, L. Samarakoon, and S. J. Broughton, ``Robust end-to-end diarization with domain adaptive training and multi-task learning,'' in _2023 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2023, pp. 1–7.
  * [16] C. Wang, J. Li, X. Fang, J. Kang, and Y. Li, ``End-to-end neural speaker diarization with absolute speaker loss,'' in _Interspeech_ , 2023, pp. 3577–3581.
  * [17] A. Plaquet and H. Bredin, ``Powerset multi-class cross entropy loss for neural speaker diarization,'' in _Interspeech_ , 2023, pp. 3222–3226.
  * [18] S. Horiguchi, S. Watanabe, P. Garcia, Y. Xue, Y. Takashima, and Y. Kawaguchi, ``Towards neural diarization for unlimited numbers of speakers using global and local attractors,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2021, pp. 98–105.
  * [19] M. Rybicka, J. Villalba, N. Dehak, and K. Kowalczyk, ``End-to-end neural speaker diarization with an iterative refinement of non-autoregressive attention-based attractors.'' in _Interspeech_ , 2022, pp. 5090–5094.
  * [20] Y. Fujita, T. Komatsu, R. Scheibler, Y. Kida, and T. Ogawa, ``Neural diarization with non-autoregressive intermediate attractors,'' in _ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2023, pp. 1–5.
  * [21] F. Hao, X. Li, and C. Zheng, ``End-to-end neural speaker diarization with an iterative adaptive attractor estimation,'' _Neural Networks_ , vol. 166, pp. 566–578, 2023.
  * [22] Z. Chen, B. Han, S. Wang, and Y. Qian, ``Attention-based encoder-decoder network for end-to-end neural speaker diarization with target speaker attractor,'' in _Interspeech_ , 2023, pp. 3552–3556.
  * [23] L. Samarakoon, S. J. Broughton, M. Härkönen, and I. Fung, ``Transformer attractors for robust and efficient end-to-end neural diarization,'' in _2023 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2023, pp. 1–8.
  * [24] F. Landini, M. Diez, T. Stafylakis, and L. Burget, ``Diaper: End-to-end neural diarization with perceiver-based attractors,'' _arXiv preprint arXiv:2312.04324_ , 2023.
  * [25] F. Teixeira, A. Abad, B. Raj, and I. Trancoso, ``Privacy-oriented manipulation of speaker representations,'' _arXiv preprint arXiv:2310.06652_ , 2023.
  * [26] D. Snyder, G. Chen, and D. Povey, ``Musan: A music, speech, and noise corpus,'' _arXiv preprint arXiv:1510.08484_ , 2015.
  * [27] M. Przybocki and A. Martin, ``2000 nist speaker recognition evaluation (ldc2001s97),'' _Philadelphia, New Jersey: Linguistic Data Consortium_ , 2001.
  * [28] A. D. G. Canavan and G. Zipperlen, ``Callhome american english speech ldc97s42,'' _Web Download. Philadelphia: Linguistic Data Consortium_ , 1997.
  * [29] J. G. Fiscus, G. Doddington, A. Le, G. Sanders, M. Przybocki, and D. Pallett, ``nist rich transcription evaluation data ldc2007s10,'' _Web Download. Philadelphia: Linguistic Data Consortium_ , 2007.
  * [30] NIST, ``The 2009 (rt-09) rich transcription meeting recognition evaluation plan,'' _http://www.itl.nist.gov/iad/ mig/tests/rt/2009/docs/rt09-meeting-eval-plan-v2.pdf_ , 2009.
  * [31] D. P. Kingma and J. Ba, ``Adam: A method for stochastic optimization,'' in _3rd International Conference on Learning Representations (ICLR)_ , 2015.
  * [32] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, ``Attention is all you need,'' _Advances in neural information processing systems_ , vol. 30, 2017.


# A note on representation-finite symmetric algebras

Karin Erdmann  Mathematical Institute, University of Oxford, ROQ, OX2 6GG,
Oxford UK<EMAIL_ADDRESS>
##  1\. Introduction

Tame symmetric algebras of period four occur as 2-blocks of group algebras
with quaternion defect groups, and algebras of quaternion type [4], and in
much more generality, as weighted surface algebras and their virtual mutations
and some other related algebras ( [6] (and [8], [9], and [10]). In [7] the
main result is the classification of tame symmetric algebras of period four
whose Gabriel quiver is 2-regular, and there is work in progress towards a
more general classification.

When the Gabriel quiver has few vertices, representation-finite algebras can
occur and it is necessary to clarify when this happens. This note gives a
characterisation, which will be used in [5]. Furthermore we describe the basic
algebras explicitly, this is based on existing literature on representation-
finite selfinjective algebras.

Recall that the syzygy $\Omega(M)$ of a module $M$ is the kernel of a minimal
projective cover of $M$. The module $M$ is periodic if $\Omega^{n}(M)\cong M$
for some $n\geq 1$, and if so, the minimal such $n$ is the period of $M$. The
algebra $\Lambda$ is periodic if $\Lambda$ is periodic as a bimodule, or
equivalently as a module for the enveloping algebra. It was proved in [3] that
all selfinjective representation-finite algebras are periodic. For general
background we refer to [2].

##  2\. A characterisation

We fix an algebraically closed field $K$. Since the properties we will study
are invariant under Morita equivalence, we can work with basic algebras, of
the form $KQ/I$ where $Q$ is a finite quiver which we can take to be
connected, and $I$ is an admissible ideal in the path algebra $KQ$.

For a vertex $i$ of $Q$, we write $i^{-}$ for the set of arrows ending at $i$,
and $i^{+}$ for the set of arrows starting at $i$. We work with right modules.
We denote by $P_{i}$ the indecomposable projective module corresponding to the
vertex $i$, and by $S_{i}$ its simple quotient.

We assume each simple module $S_{i}$ has $\Omega$-period diving $4$. Then
there is an exact sequence

$0\to S_{i}\to P_{i}\stackrel{{\scriptstyle
d_{3}}}{{\to}}P_{i}^{-}\stackrel{{\scriptstyle
d_{2}}}{{\to}}P_{i}^{+}\stackrel{{\scriptstyle d_{1}}}{{\to}}P_{i}\to S_{i}\to
0$

with ${\rm Im}(d_{k})\cong\Omega^{k}(S_{i})$ for $1\leq k\leq 3$. Here
$P_{i}^{+}=\oplus_{\alpha:i\to j}P_{j}$ and $P_{i}^{-}=\oplus_{\gamma:k\to
i}P_{k}$.

######  Lemma 2.1.

Let $\Lambda$ be an indecomposable basic symmetric algebra such that all
indecomposable non-projective modules have period dividing $4$. Then the
following are equivalent:
(a) There is a simple module $S$ such that $\Omega^{2}(S)$ is simple.
(b) The quiver of $\Lambda$ has an arrow $\alpha:i\to j$ such that
$i^{+}=\\{\alpha\\}=j^{-}$.
(c) There is a vertex $i$ of the quiver with $|P_{i}|=|P_{i}^{+}|=|P_{i}^{-}|$
where $|-|$ denotes dimension.
If these hold then $\Lambda$ has finite type.

Proof

(a) $\Rightarrow$ (b) Assume $\Omega^{2}(S_{i})\cong S_{j}.$ Then in the exact
sequence for $S_{i}$ we have $P_{j}\cong P_{i}^{+}$, so that there is a unique
arrow $\alpha$ say, starting at $i$ and it ends at $j$. Since $S_{i}$ has
period four, $\Omega^{2}(S_{j})\cong\Omega^{4}(S_{i})\cong S_{i}$. Consider
the exact sequence for $S_{j}$, the module $P_{j}^{-}$ of this sequence is
isomorphic to $P_{i}$. That is there is a unique arrow in $Q$ ending at $j$,
and this must be $\alpha$.

(b) $\Rightarrow$ (a) If (b) holds then
$\Omega(S_{i})\cong\alpha\Lambda\cong\Omega^{-1}(S_{j})$ and
$\Omega^{2}(S_{i})\cong S_{j}$.

(c) $\Leftrightarrow$ (a) Taking dimensions in the exact sequence of $S_{i}$
we have

$|\Omega^{2}(S_{i})|+|P_{i}|=|P_{i}^{+}|+|S_{i}|=|P_{i}^{-}|+|S_{i}|$

Hence $|\Omega^{2}(S_{i})|=1$ if and only if
$|P_{i}^{+}|=|P_{i}^{-}|=|P_{i}|$.

Any of the conditions above implies $\Lambda$ has finite type: Suppose
$\Omega^{2}(S_{i})\cong S_{j}$ for $S_{i}$ and $S_{j}$ simple, then also
$\Omega^{2}(S_{j})\cong S_{i}$. By part (b), the quiver has a connected
component either a vertex with a loop, or two vertices with a pair of arrows.
This must be all of $Q$ since the algebra is indecomposable. Any finite-
dimensional algebra with such a quiver is a Nakayama algebra and is of finite
type (see for example [11] Vol I, Ch I.10.3, or see in the proof of the
following Lemma). $\Box$

##  3\. The basic algebras

We describe now which symmetric algebras of finite type have only
indecomposable modules of $\Omega$-period dividing four. The algebras
appearing in the previous Lemma appear in (i) and (ii), and there is only one
other family of algebras.

######  Lemma 3.1.

Assume $\Lambda$ is symmetric indecomposable of finite type such that all
indecomposable non-projective modules have $\Omega$-period dividing $4$. Then
$\Lambda$ is one of the following.
(i) $\Lambda\cong K[T]/(T^{n})$ for some $n\geq 1$.
(ii) $\Lambda\cong KQ/I$ where $Q$ is the quiver

$\textstyle{1\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\alpha}$$\textstyle{2\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\beta}$

and $I=\langle(\alpha\beta)^{n}\alpha,(\beta\alpha)^{n}\beta\rangle$ for some
$n\geq 1$.
(iii) $\Lambda=KQ/I$ where $Q$

$\textstyle{1\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\rho}$$\scriptstyle{\alpha}$$\textstyle{2\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\beta}$

and $I=\langle\rho\alpha,\beta\rho,\alpha\beta-\rho^{t}\rangle$ for some
$t\geq 2$.

Proof If $\Lambda$ is a local algebra of finite type over an algebraically
closed field then $\Lambda$ is as in part (i), by elementary algebra. We
assume now that $\Lambda$ is not local.

Let $\Gamma_{s}$ be the stable Auslander-Reiten quiver of $\Lambda$. By work
of Riedtmann, this is isomorphic to a quotient of $\mathbb{Z}\Delta$ where
$\Delta$ is a Dynkin diagram of type $\mathbb{A},\mathbb{D},\mathbb{E}$. We
start by excluding types $\mathbb{D}$ and $\mathbb{E}$.

By [3], all selfinjective algebras of finite type are periodic as algebras;
the table 5.2 of [3] lists the (possible) periods, in terms of the types as
introduced in [1]:

The algebra has type $(\Delta,f,t)$ if the stable Auslander-Reiten quiver
$\Gamma_{s}$ is isomorphic to $\mathbb{Z}\Delta/\langle g\rangle$ where $g$ is
an automorphism of $\mathbb{Z}\Delta$ of the form $g=\zeta\tau^{-r}$, with
$\zeta$ an automorphism of order $t$, and $r\in\mathbb{N}$. Furthermore,
$f=\frac{r}{m_{\Delta}}$ where $m_{\Delta}$ is $n,2n-3,11,17$ or $29$ when
$\Delta$ is $\mathbb{A}_{n},\mathbb{D}_{n}$, $\mathbb{E}_{6},\mathbb{E}_{7}$
or $\mathbb{E}_{8}$ respectively. The algebra has $s$ simple modules where
$s=nf$.

Consider the types when $\Delta=\mathbb{D}_{n}$. By inspecting the list in 5.2
of [3] one sees that the period of the algebra always has a factor $2n-3$. To
have period $4$ this must be equal to $1$. But $2n-3=1$ means $n=2$ which is
excluded since $n\geq 4$ for $\mathbb{D}_{n}$. The same arguments shows that
algebras for types $\mathbb{D}_{3m},\mathbb{E}_{r}$ with $r=6,7,8$ and also
algebras of type $(\mathbb{A}_{2m+1},s,2)$ cannot have period four.

This leaves only the type $(\mathbb{A}_{n},s/n,1)$ where the algebra can have
period $4=\frac{2s}{(s,n+1)}.$. In this case, $f=r/m_{\Delta}=r/n$ but also
$f=s/n$ and therefore $r=s$. We assume that the algebra is symmetric,
therefore $\tau\cong\Omega^{2}$. We have period $4$, therefore $r=2$ and hence
$s=2$. Hence the quiver of the algebra has only two vertices and is connected.

To finish the proof, we will use the following properties for a finite-
dimensional symmetric algebra $\Lambda$.
(a) The quiver does not have a sink or a source.
(b) For vertices $i\neq j$ of $Q$, we have ${\rm soc}(e_{i}\Lambda)\cap
e_{i}\Lambda e_{j}=0$.
(c) If $0\neq\omega$ is a monomial in $KQ$ and belongs to the socle of
$\Lambda$ then it corresponds to a cyclic path, and all rotations of $\omega$
are non-zero and belong to the socle of $\Lambda$.

By part (a) the quiver $Q$ must contain the quiver given in (ii). The only
other possibility is the quiver given in (iii).Namely $Q$ cannot have double
arrows (this would lead to the a quotient which is a Kronecker algebra of
infinite type). There cannot be another loop since then the separated quiver
of the algebra would have infinite type.

We will now determine the ideal $I$ in the two cases.

In case (ii), using property (b) we see that the socle of $e_{1}\Lambda$ is
spanned by some $(\alpha\beta)^{n}$ and then by property (c), the monomial
$(\beta\alpha)^{n}$ spans the socle of $e_{2}\Lambda$. It follows that the
ideal $I$ is as stated (see also for example [11] vol I Ch I.10.3).

Consider now $Q$ as in (iii). First observe that in this case the simple
modules have $\Omega$-period $4$, to see this we can apply the Lemma (which
shows that the algebra is as in (ii)).

Let $H:={\rm rad}(e_{2}\Lambda)/{\rm soc}(e_{2}\Lambda)=\beta\Lambda/{\rm
soc}(e_{2}\Lambda)$ We will show using the Auslander-Reiten quiver that $H$ is
simple, isomorphic to $S_{1}$. We label the rows ot $\Gamma_{s}$ by
$1,2,\ldots,n$. Note that $M$ is in row $k$ if and only if $\Omega(M)$ is in
row $n-k+1$.

Recall that there is an almost split sequence

$0\to\Omega(S_{2})\to H\oplus P_{2}\to\Omega^{-1}(S_{2})\to 0$

Since $H$ is indecomposable, $\Omega(S_{2})$ has a unique successor in
$\Gamma_{s}$, and then also $S_{2}$ has a unique successor and therefore
belongs to row $1$ or to row $n$. Say $S_{2}$ is in row $1$ and then
$\Omega(S_{2})$ is in row $n$.

We know $\Omega^{2}(S_{2})$ is not simple, therefore is no other projective
attached to row $n$. Assume (for a contradiction) that $H$ is not simple. The
modules in row $n-1$ are $H$ and $\Omega^{2}(H)$ and $\Omega^{2}(H)$ is the
middle term of the other almost split sequence with end terms
$\Omega^{\pm}(S_{2})$ and hence is not simple.

With the assumption, both modules in row $n-1$ are not simple. Then the
modules in row $2$, ie $\Omega^{\pm}(H)$, are not isomorphic to
$\Omega^{\pm}(S_{1})$ and there is no projective attached to row $2$.

Now it follows by using the exactness of almost split sequences, and also the
fact that irreducible maps are injective or surjective, that the lengths of
modules in rows $2,3,\ldots$ are strictly increasing. But then there is no
space for the simple module $S_{1}$, a contradiction.

Therefore $H$ must be simple isomorphic to $S_{1}$. The socle of
$e_{2}\Lambda$ is then spanned by $\beta\alpha$, and by property (b) we have
that $\beta\rho=0$.. Property (c) shows that the socle of $e_{1}\Lambda$ is
spanned by $\alpha\beta$. It is also spanned by $\rho^{t}$ for some $t\geq 2$,
and we may assume $\alpha\beta=\rho^{t}$, after possibly rescaling an arrow.
Furthermore $\rho\alpha\beta=\rho^{t+1}=0$ and therefore $\rho\alpha\in{\rm
soc}(e_{1}\Lambda)\cap e_{1}\Lambda e_{2}=0$. These imply that $I$ is of the
stated form. $\Box$

## References

  * [1] H. Asashiba, The derived equivalence classification of representation-finite selfinjective algebras, J. Algebra 214 (1) (1999), 182-221.
  * [2] I. Assem, D. Simson, A. Skowroński, Elements of the Representation Theory of Associative Algebras 1: Techniques of Representation Theory. London Math. Soc. Student Texts. vol. 65. Cambridge University Press, Cambridge 2006.
  * [3] A. Dugas, Periodic resolutions and self-injective algebras of finite type,  J. Pure Applied Algebra 214 (2010) 990-1000
  * [4] K. Erdmann, Blocks of Tame Representation Type and Related Algebras, Springer Lecture Notes in Mathematics, 1428. 1990\.
  * [5] K. Erdmann, A. Hajduk, A. Skowyrski, Tame symmetric algebras of period four in preparation.
  * [6] K. Erdmann, A. Skowroński, Weighted surface algebras. J. Algebra 505 (2018), 490–558. arXiv 1703.02346.
  * [7] K. Erdmann, A. Skowroński, Algebras of generalized quaternion type. Adv. Math. 349 (2019), 1036–1116. arXiv 1710.09640.
  * [8] K. Erdmann, A. Skowroński, Weighted surface algebras: general version, J. Algebra 544 (2020), 170–227. Corrigendum J. Algebra 569(2021), 875-889.
  * [9] T. Holm, A. Skowroński, A. Skowyrski, Virtual mutations of weighted surface algebras. J. Algebra 619(2023), 822-859.
  * [10] A. Skowroński, A. Skowyrski, Weighted generalized surface algebras. arxiv 2106.15218
  * [11] A. Skowroński, K. Yamagata, Frobenius algebras I Basic representation theory. EMS Textbooks in Mathematics, European Mathematical Society, Zürich 2011. xii + 650 pp.


# MPCGPU: Real-Time Nonlinear Model Predictive Control
through Preconditioned Conjugate Gradient on the GPU

Emre Adabag1, Miloni Atal1∗, William Gerard1∗, Brian Plancher2 This material
is based upon work supported by the National Science Foundation (under Award
2246022). Any opinions, findings, conclusions, or recommendations expressed in
this material are those of the authors and do not necessarily reflect those of
the funding organizations.1Emre Adabag, Miloni Atal, and William Gerard are
with the School of Engineering and Applied Science, Columbia Unviersity, New
York, NY. {ea2944, ma4338<EMAIL_ADDRESS>Plancher is with Barnard
College, Columbia University, New York, NY<EMAIL_ADDRESS>authors contributed equally to this work.

###### Abstract

Nonlinear Model Predictive Control (NMPC) is a state-of-the-art approach for
locomotion and manipulation which leverages trajectory optimization at each
control step. While the performance of this approach is computationally
bounded, implementations of direct trajectory optimization that use iterative
methods to solve the underlying moderately-large and sparse linear systems,
are a natural fit for parallel hardware acceleration. In this work, we
introduce MPCGPU, a GPU-accelerated, real-time NMPC solver that leverages an
accelerated preconditioned conjugate gradient (PCG) linear system solver at
its core. We show that MPCGPU increases the scalability and real-time
performance of NMPC, solving larger problems, at faster rates. In particular,
for tracking tasks using the Kuka IIWA manipulator, MPCGPU is able to scale to
kilohertz control rates with trajectories as long as 512 knot points. This is
driven by a custom PCG solver which outperforms state-of-the-art, CPU-based,
linear system solvers by at least 10x for a majority of solves and 3.6x on
average.

##  I Introduction

Nonlinear Model Predictive Control (NMPC) is a feedback control strategy which
repeatedly solves finite horizon optimal control problems (OCP) in real time,
enabling robots to adapt to changes in their environment. This approach has
seen great recent success in applications to both locomotion and manipulation
[1, 2, 3, 4, 5].

Most implementations of NMPC leverage trajectory optimization [6] to solve the
underlying optimal control problems. Two popular classes of these algorithms
are shooting methods and direct methods. Shooting methods parameterize only
the input trajectory and use Bellman’s optimality principle [7] to iteratively
solve a sequence of smaller optimization problems [8, 9]. Direct methods
explicitly represent the states, controls, dynamics, and any additional
constraints, leading to moderately-large nonlinear programs with structured
sparsity patterns [10].

There has been historical interest in parallel strategies [11] for solving
trajectory optimization problems. This is growing increasingly important with
the impending end of Moore’s Law and the end of Dennard Scaling, which have
led to a utilization wall that limits the performance a single CPU chip can
deliver [12, 13]. Several more recent efforts have shown that significant
computational benefits are possible by exploiting the natural parallelism in
the computation of the (gradients of the) dynamics and cost functions on GPUs
and FPGAs [14, 15, 16, 17, 18, 19, 20, 21]. However, multiple-shooting and
consensus approaches to computing trajectory updates at each algorithmic
iteration [22, 23, 24, 25] have only seen modest gains when implemented on
alternative hardware platforms [26, 27].

On the other hand, direct methods naturally expose more parallelism that can
be exploited through hardware acceleration. Importantly, these approaches are
computationally dominated by the solutions of a moderately-large and sparse
linear systems. Iterative methods, like the Preconditioned Conjugate Gradient
(PCG) algorithm [28], are particularly well-suited for parallel solutions of
linear systems, as they are computationally dominated by matrix-vector
products and vector reductions [29, 30], and have shown past success in
outperforming state-of-the-art CPU implementations for solving very-large
linear systems on GPUs [31, 32].

Figure 1: At each control step, MPCGPU uses a parallel construction of the
Schur complement of the KKT system, a parallel-friendly PCG solver, and a
parallel line search to provide real-time performance through GPU
acceleration.

In this work, we introduce MPCGPU, a GPU-accelerated, real-time NMPC solver
that exploits the structured sparsity and the natural parallelism in direct
trajectory optimization (see Figure 1). At our solver’s core is a custom,
accelerated implementation of PCG tuned for the Schur complement of the KKT
systems of trajectory optimization problems.

We show that MPCGPU increases the scalability and real-time performance of
NMPC, solving larger problems, at faster rates. In particular, for tracking
tasks using the Kuka IIWA manipulator, MPCGPU is able to scale to kilohertz
rates with trajectories as long as 512 knot points. This is driven by a
custom, GPU-accelerated, PCG solver which outperforms state-of-the-art, CPU-
based, linear system solvers by at least 10x for a majority of solves and 3.6x
on average. We release our software and experiments open source at:
https://github.com/a2r-lab/MPCGPU.

##  II Related Work

There has been a significant amount of prior work developing general purpose
sparse linear system solvers on the GPU both using factorization-based
approaches [33, 34, 35, 36, 37, 38, 39, 40, 41, 42], as well as iterative
methods [29, 32, 43, 44, 45, 46, 47, 48]. There has also been work developing
and implementing Block-Cyclic-Reduction and other tree-structured methods that
are optimized for block-tridiagonal systems [49, 50, 51, 52]. These general
purpose approaches have found speedups through GPU usage, but only once the
problem size grows to more than tens if not hundreds of thousands of variables
or for instances of (linear) (power-flow) problems [32, 39, 40, 41, 42, 53].
As such, these general purpose solvers are not performant for most trajectory
optimization problems (e.g., our examples in Section V have 448 to 7,168
variables).

For the nonlinear trajectory optimization problem, evolutionary, particle-
swarm, Monte-Carlo, and other sampling based approaches have been implemented
on GPUs [54, 55, 56, 57, 58, 59, 60, 61, 62]. Most prior work on gradient-
based parallel nonlinear trajectory optimization has been fully confined to
the CPU [63, 24, 22, 23, 25, 64], relied on the CPU for many of the
computations [65, 66], focused only on the problem of optimizing BLAS
functions on the GPU [67, 68], or was limited to GPU acceleration of the
naturally parallel (gradients of the) dynamics and cost functions [14, 16, 18,
19].

There are two existing lines of work which fully implemented gradient-based
nonlinear trajectory optimization on the GPU. The first leveraged shooting
based methods and found them to not expose much natural parallelism, limiting
their performance [26, 27]. The second used a Block-Cyclic-Reduction-based
direct method to exploit the particular structure exposed by position-based
dynamics [15].

This work adds to the literature by designing a GPU-accelerated, gradient-
based, direct trajectory optimization solver for standard reduced-coordinate
dynamics [69] leveraging a custom parallel PCG solver at its core.

##  III Background

###  III-A Direct Trajectory Optimization

Trajectory optimization [6], also known as numerical optimal control, solves
an (often) nonlinear optimization problem to compute a robot’s optimal path
through an environment as a series of states (x $\in\mathbb{R}^{n}$) and
controls (u $\in\mathbb{R}^{m}$). These problems assume a discrete-time
dynamical system,

$x_{k+1}=f(x_{k},u_{k},h),\quad x_{0}=x_{s},$ (1)

with a timestep $h$, and minimize an additive cost function,

$J(X,U)=\ell_{f}(x_{N})+\sum^{N-1}_{k=0}\ell(x_{k},u_{k}).$ (2)

Direct methods for trajectory optimization form a moderately-large and sparse
nonlinear program. While there are a variety of algorithmic approaches used to
solve these problems, most methods can be reduced to a three step process
which is repeated until convergence [10, 70, 71].

_Step 1:_ Form a quadratic approximation around a nominal trajectory ($x,u$),
resulting in the quadratic program (QP):

$\begin{split}\min_{\begin{subarray}{c}\delta x_{0:N}\\\ \delta
u_{0:N-1}\end{subarray}}&\frac{1}{2}\delta x_{N}^{T}Q_{N}\delta
x_{N}+q_{N}^{T}\delta x_{N}+\\\ \sum_{k=0}^{N-1}&\frac{1}{2}\delta
x_{k}^{T}Q\delta x_{k}+q^{T}\delta x_{k}+\frac{1}{2}\delta u_{k}^{T}R\delta
u_{k}+r^{T}\delta u_{k}\\\ \text{s.t.\quad}&\delta x_{k+1}-A_{k}\delta
x_{k}-B_{k}\delta u_{k}=0\quad\forall k\in[0,N),\\\ &\delta
x_{0}=x_{s}-x_{0}\end{split}$ (3)

_Step 2:_ Compute $\delta X^{*},\delta U^{*}$ by solving the KKT system:

$\begin{split}&\begin{bmatrix}G&C^{T}\\\
C&0\end{bmatrix}\begin{bmatrix}-\delta Z\\\
\lambda\end{bmatrix}=\begin{bmatrix}g\\\ c\end{bmatrix}\\\ \end{split}$ (4)

where:

$\begin{split}\delta z_{k}&=\begin{bmatrix}\delta x_{k}&\delta
u_{k}\end{bmatrix}^{T}\quad\quad\quad\delta z_{N}=\delta x_{N}\\\
G&=\begin{bmatrix}Q_{0}&&&\\\ &R_{0}&&\\\ &&\ddots&\\\ &&&Q_{N}\\\
\end{bmatrix}\\\
g&=\begin{bmatrix}q_{0}&r_{0}&q_{1}&r_{1}&\dots&q_{N}\end{bmatrix}^{T}\\\
C&=\begin{bmatrix}I&&&&&\\\ -A_{0}&-B_{0}&I&&&\\\
&&\ddots&-A_{N-1}&-B_{N-1}&I\\\ \end{bmatrix}\\\
e_{k}&=x_{k+1}-f(x_{k},u_{k})\\\
c&=\begin{bmatrix}x_{s}-x_{0}&e_{0}&e_{1}&\dots&e_{N-1}\end{bmatrix}^{T}\end{split}$

_Step 3:_ Apply the update step, $\delta X^{*},\delta U^{*}$, while ensuring
descent on the original nonlinear problem through the use of a merit-function
and a trust-region or line-search.

###  III-B The Schur Complement Method

As the KKT system (Equation 4) is a saddle point system, one approach is to
solve it through a two step process by forming the symmetric positive definite
_Schur Complement_ , $S$, which can be used to solve for $\lambda^{*}$ and
then $\delta z^{*}$:

$\begin{split}S&=-CG^{-1}C^{T}\hskip 20.0pt\gamma=c-CG^{-1}g\\\
S\lambda^{*}&=\gamma\hskip 54.0pt\delta z^{*}=G^{-1}(g-C^{T}\lambda^{*})\\\
\end{split}$ (5)

By defining the variables $\theta$, $\phi$, and $\zeta$,

$\begin{split}\theta_{k}&=A_{k}Q_{k}^{-1}A_{k}^{T}+B_{k}R_{k}^{-1}B_{k}^{T}+Q_{k+1}^{-1}\\\
\phi_{k}&=-A_{k}Q_{k}^{-1}\\\
\zeta_{k}&=-A_{k}Q_{k}^{-1}q_{k}-B_{k}R_{k}^{-1}r_{k}+Q_{k+1}^{-1}q_{k+1},\end{split}$
(6)

$S$ and $\gamma$ take the form:

$\begin{split}S&=\begin{bmatrix}Q_{0}^{-1}&\phi_{0}^{T}&&&&\\\
\phi_{0}&\theta_{0}&\phi_{1}^{T}&&&\\\
&&\ddots&\phi_{N-2}&\theta_{N-2}&\phi_{N-1}^{T}\\\
&&&&\phi_{N-1}&\theta_{N-1}\\\ \end{bmatrix}\\\
\gamma&=c+\begin{bmatrix}Q_{0}^{-1}q_{0}&\zeta_{0}&\zeta_{1}&\dots&\zeta_{N-1}\end{bmatrix}^{T}\end{split}$
(7)

Figure 2: A high level overview of MPCGPU which: 1) in parallel on the GPU
computes $S$, $\gamma$, and $\Phi^{-1}$ and stores those values in an
optimized dense format, 2) uses our GBD-PCG solver to compute $\lambda^{*}$
and reconstructs $\delta X^{*},\delta U^{*}$ through GPU-friendly matrix-
vector multiplications and vector reductions, and 3) leverages a parallel line
search to compute the final trajectory, $X,U$. This trajectory is then passed
to the (simulated) robot and the current state of the (simulated) robot is
measured and fed back into our solver which is run again, warm-started with
our last solution.

###  III-C Iterative Methods

Iterative methods solve the problem $S\lambda^{*}=\gamma$ for a given $S$ and
$\gamma$ by iteratively refining an estimate for $\lambda$ up to some
tolerance $\epsilon$. The most popular of these methods is the conjugate
gradient (CG) algorithm which has been used for state-of-the-art results on
large-scale optimziation problems on the GPU [53, 72]. The convergence rate of
CG is directly related to the spread of the eigenvalues of $S$. Thus, a
preconditoning matrix $\Phi\approx S$ is often applied to instead solve the
equivalent problem with better numerical properties:
$\Phi^{-1}S\lambda^{*}=\Phi^{-1}\gamma$. To do so, the preconditioned
conjugate gradient (PCG) algorithm leverages matrix-vector products with $S$
and $\Phi^{-1}$, as well as vector reductions, both parallel friendly
operations (see Algorithm 1).

###  III-D Graphics Processing Units (GPUs)

Compared to a multi-core CPU, a GPU is a larger set of simpler processors,
optimized for parallel execution of identical instructions. GPUs are best at
computing regular and separable computations, over large data sets, with
limited synchronization (e.g., large matrix multiplication) [73]. Our work
uses NVIDIA’s CUDA extensions to C++ [74].

##  IV The MPCGPU Solver

In this section we describe the design of the MPCGPU solver which exploits the
sparsity and natural parallelism found in direct trajectory optimization
algorithms and iterative linear system solvers. To further promote efficient
GPU acceleration, unlike generic approaches, which require a kernel launch and
CPU-GPU synchronization for each matrix operation, MPCGPU uses only three
kernels that are asynchronously queued, resulting in only a single CPU-GPU
synchronization. We also only transfer the initial and final values between
the CPU and GPU to reduce I/O overheads.

As shown in Figure 2, our approach can be broken down into a three step
process. At each control step we first compute each block row of the Schur
complement system, $S$ and $\gamma$, as well as our preconditioner,
$\Phi^{-1}$, in parallel by taking advantage of the structured sparsity of
those matrices. Next we use our custom GPU-optimized, warm-started, PCG
solver, GBD-PCG (Algorithm 2), to compute the optimal Lagrange multipliers,
$\lambda^{*}$, and reconstruct the optimal trajectory update, $\delta
X^{*},\delta U^{*}$. Finally we leverage a parallel line search to compute the
final trajectory $X,U$ which we send to the (simulated) robot for execution
and simultaneously measure the current state of the (simulated) robot to begin
our next control step. In the remainder of this section we provide further
details on our approach. Our open-source implementation can be found at:
https://github.com/a2r-lab/MPCGPU.

1:$r=\gamma-S\lambda$

2:$\tilde{r},p=\Phi^{-1}r$

3:$\eta=r^{T}\tilde{r}$

4:for iter $i=1:\text{max\\_iter}$ do

5: $\alpha=\eta/(p^{T}Sp)$

6: $r=r-\alpha Sp$

7: $\lambda=\lambda+\alpha p$

8: $\tilde{r}=\Phi^{-1}r$

9: $\eta^{\prime}=r^{T}\tilde{r}$

10: if $\eta^{\prime}<\epsilon$ then return $\lambda$

11: $\beta=\eta^{\prime}/\eta$

12: $p=\tilde{r}+\beta p$

13: $\eta=\eta^{\prime}$

14:return $\lambda$

Algorithm 1 Preconditioned Conjugate Gradient (PCG)
($S,\Phi^{-1},\gamma,\lambda,\epsilon$) $\rightarrow$ $\lambda^{*}$

Initialization

Main Loop

###  IV-A Parallel Computation of $S,\gamma$, and $\Phi^{-1}$

To efficiently compute $S,\gamma$, and $\Phi^{-1}$ on the GPU, we need to find
a naturally parallel approach to form the values as well as an efficient data
storage format that minimizes overheads. We also need to find an effective
preconditioner that is parallel-friendly in its computation.

We first leverage the block-tridiagonal structure of the Schur complement,
$S$, as shown in Equation 7, which is for the most part independent across
timesteps, $k$, for each block-row. This pattern also extends to each block-
row of the $\gamma$ vector. To further remove the need for synchronizations,
for each $k$, we also compute the only cross-timestep quantities, $Q_{k+1}$
and $q_{k+1}$. While this results in those terms being computed twice, it
still proves to be more efficient than forcing a synchronization point between
all block-rows. To ensure efficient computation of the underlying dynamics and
kinematic quantities, we leverage the GRiD library, which was shown to
outperform state-of-the-art CPU libraries even when taking into account I/O
overheads [18].

We further parallelize across and within the many small matrix inversions and
matrix multiplications within each parallel block-row computation. Leveraging
best practices [16, 18], we also group together the various types of
mathematical operations, storing intermediate values in shared memory, and re-
ordering computations where needed.

We leverage the Symmetric Stair Preconditioner [30, 75] which is a parallel-
friendly preconditioner optimized for block-tridiagonal systems wich has an
analytical inverse and results in the block-tridiagonal matrix:

$\Phi^{-1}=\begin{bmatrix}Q_{0}&-Q_{0}\phi_{0}^{T}\theta_{0}^{-1}&\\\
-\theta_{0}^{-1}\phi_{0}Q_{0}^{-1}&\;\;\theta_{0}^{-1}\;\;&-\theta_{0}^{-1}\phi_{1}^{T}\theta_{1}^{-1}\\\
&-\theta_{1}^{-1}\phi_{1}\theta_{0}^{-1}&\theta_{1}^{-1}\\\ &&\ddots\\\
\end{bmatrix}\\\ $ (8)

The structure of $\Phi^{-1}$ also permits mostly parallel computation as only
the values of each $\theta_{k}^{-1}$ need to be shared across timesteps. Our
approach thus requires only a single global synchronization across blocks, and
allows us to store the block-tridiagonal $S$ and $\Phi^{-1}$ matrices in a
custom, compressed, dense format for increased IO bandwidth and memory
efficiency.

###  IV-B GPU Parallel PCG for Block-Tridiagonal Systems

The core of our solver is a custom GPU Parallel PCG implementation
specifically optimized for block-tridiagonal systems, GBD-PCG (Algorithm 2).
That is, we leverage the sparsity structure of $S$ and $\Phi^{-1}$ to maximize
cache usage and natural parallelism resulting in a refactored, low-latency
implementation with minimal synchronizations. These optimizations can be
leveraged in the most computationally expensive part of the algorithm, the
large matrix-vector products in lines 5, 6, and 8 of Algorithm 1, as each
element of the product depends on at most $3n_{b}$ elements from the matrix
and $3n_{b}$ elements from the vector, where $n_{b}$ is the block dimension.
We exploit this by grouping threads that access similar elements into thread
blocks and storing $S$, $\Phi^{-1}$, and all PCG iterates concurrently in
shared (cache) memory on the GPU. We also operate as many steps of the
algorithm fully in parallel as possible between the thread synchronizations
needed for the parallel reductions of scalar values on lines 6, 12, and 21 of
Algorithm 2. Similarly, we only use device memory (RAM) for those scalar
reductions and for the values of $p$ and $r$ that need to be shared between
blocks on lines 9 and 18 of Algorithm 2. This means that the choice of a
sparse preconditioner not only enables its efficient computation and memory
storage, but also reduces the number of synchronizations and amount of memory
that needs to be shared through RAM during each PCG iterate. This holistic co-
design across algorithm stages is part of the reason why MPCGPU is so
performant. Finally, we warm-start the values for $\lambda$ based on the
previous solve which we found greatly increased overall performance by
reducing the number of PCG iterations needed for convergence.

1:for block $b=0:N$ in parallel do

2: $r_{b}=\gamma_{b}-S_{b}\lambda_{b-1:b+1}$

3: Load $r_{b-1:b-1}$, $r_{b+1:b+1}$

4: $\tilde{r}_{b},p_{b}=\Phi^{-1}_{b}r_{b-1:b+1}$

5: $\eta_{b}=r_{b}^{T}\tilde{r_{b}}$

6:$\eta=\textbf{\text{ParallelReduce}}(\eta_{b})$

7:for iter $i=1:\text{max\\_iter}$ do

8: for block $b=0:N$ in parallel do

9: Load $p_{b-1}$, $p_{b+1}$

10: $\Upsilon_{b}=S_{b}p_{b-1:b+1}$

11: $\upsilon_{b}=p_{b}\Upsilon_{b}$

12: $\upsilon=\textbf{\text{ParallelReduce}}(\upsilon_{b})$

13: for block $b=0:N$ in parallel do

14: $\alpha=\eta/\upsilon$

15: $\lambda_{b}=\lambda_{b}+\alpha p_{b}$

16: $r_{b}=r_{b}-\alpha\Upsilon_{b}$

17: for block $b=0:N$ in parallel do

18: Load $r_{b-1:b-1}$, $r_{b+1:b+1}$

19: $\tilde{r}_{b}=\Phi^{-1}_{b}r_{b}$

20: $\eta^{\prime}_{b}=r^{T}_{b}\tilde{r}_{b}$

21: $\eta^{\prime}=\textbf{\text{ParallelReduce}}(\eta^{\prime}_{b})$

22: if $\eta^{\prime}<\epsilon$ then return $\lambda$

23: for block $b=0:N$ in parallel do

24: $\beta=\eta^{\prime}/\eta$

25: $p_{b}=\tilde{r}^{\prime}_{b}+\beta p_{b}$

26: $\eta=\eta^{\prime}$

27:return $\lambda$

Algorithm 2 GPU Parallel PCG for Block-Tridiagonal Systems (GBD-PCG)
($S,\Phi^{-1},\gamma,\lambda,\epsilon$) $\rightarrow$ $\lambda^{*}$

Initiali- zation

Main Loop

###  IV-C Parallel Line Search

We leverage a parallel line search, computing all possible iterates for
$\alpha\in\mathbb{A}$ in parallel and selecting the iterate with the best
value according to its L1 merit function [10]:

$M(Z;\mu)=J(Z)+\mu|c|.$ (9)

This is allows MPCGPU to evaluate all possible line search iterates in the
same amount of time as it would take the GPU to compute a single iterate under
a standard backtracking implementation. Importantly, this approach not only
reduces latency of this algorithmic step, but has also been shown to improve
the convergence of NMPC on similar whole-body trajectory tracking problems
[26].

##  V Results

In this section we present a two-part evaluation of MPCGPU through a case
study of online, dynamic, multi-goal, end-effector position tracking using
whole-body NMPC for a simulated Kuka IIWA manipulator. First, we compare the
performance of our underlying GBD-PCG iterative linear system solver with the
state-of-the-art, CPU-based, QDLDL solver [76]. Second, we show how the end-
to-end performance enabled by MPCGPU allows us to scale to long time horizons
and fast control rates. Source code accompanying this evaluation can be found
at  https://github.com/a2r-lab/MPCGPU.

###  V-A Methodology

Results were collected on high-performance workstation with a $3.2$GHz 16-core
Intel i9-12900K and a $2.2$GHz NVIDIA GeForce RTX 4090 GPU running Ubuntu
22.04 and CUDA 12.1. Code was compiled with g++11.4, and time was measured
with the Linux system call clock_gettime(), using CLOCK_MONOTONIC as the
source. Our performance analysis is drawn from 100 NMPC trials of a 10 second,
5 goal, pick-and-place circuit for a simulated Kuka IIWA-14 (see Figure 3).
Each NMPC trial consists of thousands of linear system solves which are needed
for the many iterations of the underlying trajectory optimization problem for
end-effector position tracking solved at each control step. All hyperparameter
values can be found in our open-source source code and resulted in an average
tracking error of $\sim$10cm, providing similar performance as previous
experiments with GPU accelerated NMPC [27]. In particular, we note that all
solvers used the same quadratic cost functions for each amount of knot points,
and solver-specific hyperparameter values were independently tuned for maximal
performance.

###  V-B Linear System Solver Performance

We evaluate the performance of our underlying GBD-PCG linear system solver
over its thousands of solves during each of our 100 NMPC trials running at a
500hz control rate and compare its performance to the state-of-the-art CPU-
based QDLDL solver [76] operating in the same context.

Average Solve Time: Our results show that our GPU based solver outperforms
QDLDL across most problem sizes and is only marginally slower at the smallest
problem size, obtaining as much as a 3.6x average speedup (see Figure 4). This
is driven by the GPU’s ability to leverage large scale parallelism to
gracefully scale to larger problem sizes. We note that the speedup plateaus at
256 knot points as we begin to run out of hardware resources on our specific
GPU. These results show that unlike generic approaches [32, 53], our domain-
specific co-design approach enables the GPU to outperform the CPU even on
moderately-sized linear systems (our experiments range from only 448 to 7,168
variables).

Figure 3: Screenshots from the 5 goal, pick-and-place circuit used in our
experiments. Figure 4: Average linear system solve time for the thousands of
solves within each of the 100 iterations of NMPC. As the problem scales, so
does the advantage of GBD-PCG over QDLDL achieving up to a 3.6x average
speedup. Figure 5: Cumulative density function of the linear system solve
times for a trajectory length of 128 for QDLDL and GBD-PCG under multiple
different exit tolerances, $\epsilon$. We note that the bi-modal distribution
for GBD-PCG is often much faster than the uni-modal distribution for QDLDL.
For example, for $\epsilon=1e^{-4}$, 65% of GBD-PCG solves are $\geq$10x
faster than the fastest QDLDL solve, and the slowest GBD-PCG solve is only
2.5x slower than the slowest QDLDL solve (with only 10% of solves $\geq$2x
slower).

Performance Distribution: Importantly, in most cases, the speedup is much
larger than this. This is because iterative methods have a variable runtime as
they can exit early depending upon the exit tolerance, $\epsilon$. We
demonstrate this using the 128 knot point problem as a case study in Figure 5.

We plot the distribution of solve times for QDLDL against GPU-PCG
$\epsilon=1e^{-4}$, resulting in our 1.9x average speedup in Figure 4, as well
as $\epsilon=5e^{-5}$ and $\epsilon=1e^{-5}$. QDLDL presents a uni-modal
timing distribution with almost all solves occurring between 280 and 305
$\mu$s. GBD-PCG, on the other hand, presents a bi-modal timing distribution
clustered both much faster and a little slower than QDLDL. For example, for
$\epsilon=1e^{-4}$, 65% of GBD-PCG solves are $\geq$10x faster than the
fastest QDLDL solve, and the slowest GBD-PCG solve is only 2.5x slower than
the slowest QDLDL solve (with only 10% of solves $\geq$2x slower).

Furthermore, while all values of $\epsilon$ shown in the plot were able to
successfully track the target trajectory, the less exact of a solve was
required, the more of the distribution mass shifted to being $\geq$10x faster
than QDLDL (65%, 52%, 20% for $1e^{-4},5e^{-5},1e^{-5}$ respectively).
However, when $\epsilon$ was reduced even farther, our entire NMPC controller
was unable to accurately track our target trajectory. These results present
interesting directions for future work to find ways to eliminate the second
slower mode of the solve time distribution while ensuring robust NMPC
convergence.

###  V-C End-to-End NMPC Performance

To validate efficacy for use in NMPC for robotics applications, we also
demonstrate the impact of our approach on the number of iterations of MPCGPU
we could achieve at each control step for varying control rates and trajectory
lengths. Figure 6 shows the resulting number of average trajectory
optimization solver iterations we can compute while meeting the specified
control rates and trajectory lengths using both our GBD-PCG solver as well as
QDLDL to solve the thousands of underlying linear systems.111We note that in
the QDLDL case, as the NMPC loop is running on the GPU, data needs to be
copied onto the host before executing the solve and converted into the sparse
CSR matrix format. To ensure fair comparisons and avoid overheads for
unnecessary data transfers and transformations, we implemented a variant of
our parallel Schur complement computation which directly stores data in the
CSR format expected by QDLDL.

Regardless of the linear system solver, our GPU-first approach, with both fast
parallel construction of the Schur complement and fast parallel computation of
the line search, enables trajectories as long as $128$ knot points to operate
at a $1$kHz control rate, and achieve at least $4$ iterations at a $500$Hz
control rate, for a per-iteration rate of $2$kHz. Furthermore, similar to what
we witness in the case of average linear system solve times, as the problem
gets larger and the control rate increases, our fully GPU-based approach is
increasingly performant. Highlights include our approach’s ability to scale to
$512$ knot points at a $1$kHz control rate and execute $8$ iterations for 128
knot points at a $500$Hz control rate, for a per-iteration rate of $4$kHz.
These results compare favorably to previously reported results in the
literature of about $500$hz to $1$kHz per-iteration rates for trajectories of
$30$ to $120$ knot points using state-of-the-art CPU-based [77, 78] and GPU-
based [26] solvers for similar NMPC tasks. As such, our GPU-first approach
opens up the possibility for our NMPC solver to either leverage longer horizon
trajectories, produce more optimal solutions at each control step, run at
faster control rates, or include some combination of those highly beneficial
traits.

Figure 6: Average number of trajectory optimization iterations of MPCGPU
executed at each control step for varying control rates and trajectory
lengths. Our GPU-first approach to Schur complement construction and our
parallel line search enables high control rates and long trajectories
regardless of the underlying solver. However, the improved scalability of GBD-
PCG enables our approach to scale to $512$ knot points at $1$kHz and execute
$8$ iterations for 128 knot points at $500$Hz, for a per-iteration rate of
$4$kHz.

##  VI Conclusion and Future Work

In this work, we introduce MPCGPU, a GPU-accelerated, real-time NMPC solver
built around a parallel PCG solver. MPCGPU exploits the structured sparsity
and natural parallelism in both direct trajectory optimization algorithms and
iterative linear system solvers. Our experiments show that our approach is
able to scale NMPC to larger problems, and operate it at faster rates, than is
possible with existing state-of-the-art solvers. In particular, for tracking
tasks using the Kuka IIWA manipulator, MPCGPU is able to scale to kilohertz
control rates with trajectories as long as 512 knot points. For this problem,
our GPU-based PCG solver outperforms a state-of-the-art CPU-based linear
system solver by as much as 10x for a majority of solves and 3.6x on average.

There are many promising directions for future work to improve the
functionality and usability of our approach. Most importantly, like all
iterative methods, our GPU-based PCG solver exhibits variability in its solve
times. Future work which learns when to leverage iterative methods vs.
factorization-based methods, or which learns dynamic values for
hyperparameters to reduce the worst-case runtimes, without sacrificing overall
NMPC robustness, would greatly improve average-case performance. Furthermore,
it would be interesting to explore the performance implications of adding
additional constraints either through expanding the KKT system or through
augmented Lagrangian methods [10, 79]. Finally, we would like to evaluate our
approach on physical robots at the edge using low-power GPU platforms such as
the NVIDIA Jetson [80].

## References

  * [1] F. R. Hogan, E. R. Grau, and A. Rodriguez, “Reactive planar manipulation with convex hybrid mpc,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 247–253.
  * [2] J.-P. Sleiman, F. Farshidian, M. V. Minniti, and M. Hutter, “A unified mpc framework for whole-body dynamic locomotion and manipulation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 4688–4695, 2021.
  * [3] M. Tranzatto, F. Mascarich, L. Bernreiter, C. Godinho, M. Camurri, S. Khattak, T. Dang, V. Reijgwart, J. Loeje, D. Wisth, S. Zimmermann, H. Nguyen, M. Fehr, L. Solanka, R. Buchanan, M. Bjelonic, N. Khedekar, M. Valceschini, F. Jenelten, M. Dharmadhikari, T. Homberger, P. D. Petris, L. Wellhausen, M. Kulkarni, T. Miki, S. Hirsch, M. Montenegro, C. Papachristos, F. Tresoldi, J. Carius, G. Valsecchi, J. Lee, K. Meyer, X. Wu, J. Nieto, A. Smith, M. Hutter, R. Siegwart, M. Mueller, M. Fallon, and K. Alexis, “Cerberus: Autonomous legged and aerial robotic exploration in the tunnel and urban circuits of the darpa subterranean challenge,” _arXiv preprint arXiv:2201.07067_ , 2022.
  * [4] P. M. Wensing, M. Posa, Y. Hu, A. Escande, N. Mansard, and A. Del Prete, “Optimization-based control for dynamic legged robots,” _arXiv preprint arXiv:2211.11644_ , 2022.
  * [5] S. Kuindersma, “Taskable agility: Making useful dynamic behavior easier to create,” Princeton Robotics Seminar, April 2023.
  * [6] J. T. Betts, _Practical Methods for Optimal Control Using Nonlinear Programming_ , ser. Advances in Design and Control. Society for Industrial and Applied Mathematics (SIAM), vol. 3.
  * [7] R. Bellman, _Dynamic Programming_. Dover.
  * [8] D. Q. Mayne, “A second-order gradient method of optimizing non- linear discrete time systems,” vol. 3, p. 8595.
  * [9] D. H. Jacobson and D. Q. Mayne, “Differential dynamic programming,” 1970.
  * [10] J. Nocedal and S. J. Wright, _Numerical Optimization_ , 2nd ed. Springer.
  * [11] J. T. Betts and W. P. Huffman, “Trajectory optimization on a parallel processor,” vol. 14, no. 2, pp. 431–439.
  * [12] H. Esmaeilzadeh, E. Blem, R. St. Amant, K. Sankaralingam, and D. Burger, “Dark Silicon and the End of Multicore Scaling,” in _Proceedings of the 38th Annual International Symposium on Computer Architecture_ , ser. ISCA ’11. ACM, pp. 365–376.
  * [13] G. Venkatesh, J. Sampson, N. Goulding, S. Garcia, V. Bryksin, J. Lugo-Martinez, S. Swanson, and M. B. Taylor, “Conservation Cores: Reducing the Energy of Mature Computations,” in _Proceedings of the Fifteenth Edition of ASPLOS on Architectural Support for Programming Languages and Operating Systems_ , ser. ASPLOS XV. ACM, pp. 205–218.
  * [14] T. Antony and M. J. Grant, “Rapid Indirect Trajectory Optimization on Highly Parallel Computing Architectures,” vol. 54, no. 5, pp. 1081–1091.
  * [15] Z. Pan, B. Ren, and D. Manocha, “Gpu-based contact-aware trajectory optimization using a smooth force model,” in _Proceedings of the 18th Annual ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , ser. SCA ’19. New York, NY, USA: ACM, 2019, pp. 4:1–4:12.
  * [16] B. Plancher, S. M. Neuman, T. Bourgeat, S. Kuindersma, S. Devadas, and V. J. Reddi, “Accelerating robot dynamics gradients on a cpu, gpu, and fpga,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2335–2342, 2021\.
  * [17] S. M. Neuman, B. Plancher, T. Bourgeat, T. Tambe, S. Devadas, and V. J. Reddi, “Robomorphic computing: A design methodology for domain-specific accelerators parameterized by robot morphology,” ser. ASPLOS 2021. New York, NY, USA: Association for Computing Machinery, 2021, p. 674–686. [Online]. Available: https://doi-org.ezp-prod1.hul.harvard.edu/10.1145/3445814.3446746
  * [18] B. Plancher, S. M. Neuman, R. Ghosal, S. Kuindersma, and V. J. Reddi, “GRiD: GPU-Accelerated Rigid Body Dynamics with Analytical Gradients,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, pp. 6253–6260. [Online]. Available: https://ieeexplore.ieee.org/document/9812384/
  * [19] Y. Lee, M. Cho, and K.-S. Kim, “Gpu-parallelized iterative lqr with input constraints for fast collision avoidance of autonomous vehicles,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 4797–4804.
  * [20] S. M. Neuman, R. Ghosal, T. Bourgeat, B. Plancher, and V. J. Reddi, “Roboshape: Using topology patterns to scalably and flexibly deploy accelerators across robots,” in _Proceedings of the 50th Annual International Symposium on Computer Architecture_ , ser. ISCA ’23. New York, NY, USA: Association for Computing Machinery, 2023. [Online]. Available: https://doi.org/10.1145/3579371.3589104
  * [21] Y. Yang, X. Chen, and Y. Han, “Rbdcore: Robot rigid body dynamics accelerator with multifunctional pipelines,” _arXiv preprint arXiv:2307.02274_ , 2023\.
  * [22] M. Giftthaler, M. Neunert, M. Stäuble, J. Buchli, and M. Diehl, “A Family of Iterative Gauss-Newton Shooting Methods for Nonlinear Optimal Control.” [Online]. Available: http://arxiv.org/abs/1711.11006
  * [23] F. Farshidian, E. Jelavic, A. Satapathy, M. Giftthaler, and J. Buchli, “Real-time motion planning of legged robots: A model predictive control approach,” in _2017 IEEE-RAS 17th International Conference on Humanoid Robotics_.
  * [24] D. Kouzoupis, R. Quirynen, B. Houska, and M. Diehl, “A Block Based ALADIN Scheme for Highly Parallelizable Direct Optimal Control,” in _Proceedings of the American Control Conference_.
  * [25] Y. Jiang, J. Oravec, B. Houska, and M. Kvasnica, “Parallel mpc for linear systems with input constraints,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 7, pp. 3401–3408, 2020.
  * [26] B. Plancher and S. Kuindersma, “A Performance Analysis of Parallel Differential Dynamic Programming on a GPU,” in _International Workshop on the Algorithmic Foundations of Robotics (WAFR)_.
  * [27] ——, “Realtime model predictive control using parallel ddp on a gpu,” in _Toward Online Optimal Control of Dynamic Robots Workshop at the 2019 International Conference on Robotics and Automation (ICRA)_ , Montreal, Canada, May. 2019.
  * [28] S. C. Eisenstat, “Efficient implementation of a class of preconditioned conjugate gradient methods,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 2, no. 1, pp. 1–4, 1981.
  * [29] Y. Saad, _Iterative methods for sparse linear systems_. SIAM, 2003.
  * [30] B. Plancher, “Gpu acceleration for real-time, whole-body, nonlinear model predictive control,” Ph.D. dissertation, Harvard University, Cambridge, MA, USA, April 2022.
  * [31] R. Helfenstein and J. Koko, “Parallel preconditioned conjugate gradient algorithm on GPU,” vol. 236, no. 15, pp. 3584–3590. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0377042711002196
  * [32] M. Schubiger, G. Banjac, and J. Lygeros, “Gpu acceleration of admm for large-scale quadratic programming,” _Journal of Parallel and Distributed Computing_ , vol. 144, pp. 55–67, 2020.
  * [33] J. H. Jung and D. P. O’Leary, “Cholesky decomposition and linear programming on a gpu,” _Scholarly Paper, University of Maryland_ , 2006.
  * [34] D. Yang, G. D. Peterson, and H. Li, “Compressed sensing and cholesky decomposition on fpgas and gpus,” _Parallel Computing_ , vol. 38, no. 8, pp. 421–437, 2012.
  * [35] I. E. Venetis, A. Kouris, A. Sobczyk, E. Gallopoulos, and A. H. Sameh, “A direct tridiagonal solver based on givens rotations for gpu architectures,” _Parallel Computing_ , vol. 49, pp. 101–116, 2015.
  * [36] J. D. Hogg, E. Ovtchinnikov, and J. A. Scott, “A sparse symmetric indefinite direct solver for gpu architectures,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 42, no. 1, pp. 1–25, 2016.
  * [37] X. Hu, C. C. Douglas, R. Lumley, and M. Seo, “Gpu accelerated sequential quadratic programming,” in _2017 16th International Symposium on Distributed Computing and Applications to Business, Engineering and Science (DCABES)_. IEEE, 2017, pp. 3–6.
  * [38] S. N. Yeralan, T. A. Davis, W. M. Sid-Lakhdar, and S. Ranka, “Algorithm 980: Sparse qr factorization on the gpu,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 44, no. 2, pp. 1–29, 2017.
  * [39] K. Świrydowicz, E. Darve, W. Jones, J. Maack, S. Regev, M. A. Saunders, S. J. Thomas, and S. Peleš, “Linear solvers for power grid optimization problems: a review of gpu-accelerated linear solvers,” _Parallel Computing_ , vol. 111, p. 102870, 2022.
  * [40] D. Cole, S. Shin, F. Pacaud, V. M. Zavala, and M. Anitescu, “Exploiting gpu/simd architectures for solving linear-quadratic mpc problems,” in _2023 American Control Conference (ACC)_. IEEE, 2023, pp. 3995–4000.
  * [41] S. Shin, F. Pacaud, and M. Anitescu, “Accelerating optimal power flow with gpus: Simd abstraction of nonlinear programs and condensed-space interior-point methods,” _arXiv preprint arXiv:2307.16830_ , 2023.
  * [42] F. Pacaud, S. Shin, M. Schanen, D. A. Maldonado, and M. Anitescu, “Accelerating condensed interior-point methods on simd/gpu architectures,” _Journal of Optimization Theory and Applications_ , pp. 1–20, 2023.
  * [43] J. Bolz, I. Farmer, E. Grinspun, and P. Schröoder, “Sparse Matrix Solvers on the GPU: Conjugate Gradients and Multigrid,” in _ACM SIGGRAPH 2003 Papers_ , ser. SIGGRAPH ’03. ACM, pp. 917–924. [Online]. Available: http://doi.acm.org/10.1145/1201775.882364
  * [44] H. Liu, J.-H. Seo, R. Mittal, and H. H. Huang, “Gpu-accelerated scalable solver for banded linear systems,” in _2013 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2013, pp. 1–8.
  * [45] H. Anzt, M. Gates, J. Dongarra, M. Kreutzer, G. Wellein, and M. Köhler, “Preconditioned krylov solvers on gpus,” _Parallel Computing_ , 05 2017\.
  * [46] H. Anzt, M. Kreutzer, E. Ponce, G. D. Peterson, G. Wellein, and J. Dongarra, “Optimization and performance evaluation of the idr iterative krylov solver on gpus,” _The International Journal of High Performance Computing Applications_ , vol. 32, no. 2, pp. 220–230, 2018.
  * [47] G. Flegar _et al._ , “Sparse linear system solvers on gpus: Parallel preconditioning, workload balancing, and communication reduction,” Ph.D. dissertation, Universitat Jaume I, 2019.
  * [48] M. Tiwari and S. Vadhiyar, “Strategies for efficient execution of pipelined conjugate gradient method on gpu systems,” in _International Conference on High Performance Computing_. Springer, 2022, pp. 77–89.
  * [49] L.-W. Chang, J. A. Stratton, H.-S. Kim, and W.-M. W. Hwu, “A scalable, numerically stable, high-performance tridiagonal solver using gpus,” in _SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_. IEEE, 2012, pp. 1–11.
  * [50] L.-W. Chang and W.-m. W. Hwu, _A Guide for Implementing Tridiagonal Solvers on GPUs_. Springer International Publishing, 2014, pp. 29–44. [Online]. Available: https://doi.org/10.1007/978-3-319-06548-9˙2
  * [51] A. P. Diéguez, M. Amor, and R. Doallo, “New tridiagonal systems solvers on gpu architectures,” in _2015 IEEE 22nd International Conference on High Performance Computing (HiPC)_. IEEE, 2015, pp. 85–94.
  * [52] A. Lamas Daviña and J. Roman, “Mpi-cuda parallel linear solvers for block-tridiagonal matrices in the context of slepc’s eigensolvers,” _Parallel computing_ , vol. 74, pp. 118–135, 2018.
  * [53] M. Naumov, “Incomplete-lu and cholesky preconditioned iterative methods using cusparse and cublas,” _Nvidia white paper_ , vol. 3, 2011.
  * [54] S. Heinrich, A. Zoufahl, and R. Rojas, “Real-time trajectory optimization under motion uncertainty using a GPU,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 3572–3577.
  * [55] Q. Wu, F. Xiong, F. Wang, and Y. Xiong, “Parallel particle swarm optimization on a graphics processing unit with application to trajectory optimization,” _Engineering Optimization_ , vol. 48, no. 10, pp. 1679–1692, 2016.
  * [56] G. Williams, A. Aldrich, and E. A. Theodorou, “Model predictive path integral control: From theory to parallel computation,” _Journal of Guidance, Control, and Dynamics_ , vol. 40, no. 2, pp. 344–357, 2017.
  * [57] D.-K. Phung, B. Hérissé, J. Marzat, and S. Bertrand, “Model Predictive Control for Autonomous Navigation Using Embedded Graphics Processing Unit,” vol. 50, no. 1, pp. 11 883–11 888. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S2405896317319614
  * [58] P. Hyatt and M. D. Killpack, “Real-time evolutionary model predictive control using a graphics processing unit,” in _2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids)_. IEEE, 2017, pp. 569–576.
  * [59] S. Ohyama and H. Date, “Parallelized nonlinear model predictive control on gpu,” in _2017 11th Asian Control Conference (ASCC)_. IEEE, 2017, pp. 1620–1625.
  * [60] K. M. M. Rathai, O. Sename, and M. Alamir, “Gpu-based parameterized nmpc scheme for control of half car vehicle with semi-active suspension system,” _IEEE Control Systems Letters_ , vol. 3, no. 3, pp. 631–636, 2019.
  * [61] Y. Wang, X. Luo, F. Zhang, and S. Wang, “Gpu-based model predictive control for continuous casting spray cooling control system using particle swarm optimization,” _Control Engineering Practice_ , vol. 84, pp. 349–364, 2019\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S096706611830710X
  * [62] P. Hyatt, C. S. Williams, and M. D. Killpack, “Parameterized and gpu-parallelized real-time model predictive control for high degree of freedom robots,” _arXiv preprint arXiv:2001.04931_ , 2020.
  * [63] J. V. Frasch, M. Vukov, H. J. Ferreau, and M. Diehl, “A dual newton strategy for the efficient solution of sparse quadratic programs arising in sqp-based nonlinear mpc,” _Optimization Online 3972_ , 2013.
  * [64] A. Astudillo, J. Gillis, G. Pipeleers, W. Decré, and J. Swevers, “Speed-up of nonlinear model predictive control for robot manipulators using task and data parallelism,” in _2022 IEEE 17th International Conference on Advanced Motion Control (AMC)_ , 2022, pp. 201–206.
  * [65] Y. Gang and L. Mingguang, “Acceleration of mpc using graphic processing unit,” in _Proceedings of 2012 2nd International Conference on Computer Science and Network Technology_. IEEE, 2012, pp. 1001–1004.
  * [66] N. F. Gade-Nielsen, J. B. Jørgensen, and B. Dammann, “Mpc toolbox with gpu accelerated optimization algorithms,” in _10th European workshop on advanced control and diagnosis_. Technical University of Denmark, 2012.
  * [67] Y. Huang, K. V. Ling, and S. See, “Solving Quadratic Programming Problems on Graphics Processing Unit.”
  * [68] L. Yu, A. Goldsmith, and S. Di Cairano, “Efficient Convex Optimization on GPUs for Embedded Model Predictive Control,” in _Proceedings of the General Purpose GPUs_ , ser. GPGPU-10. ACM, pp. 12–21. [Online]. Available: http://doi.acm.org/10.1145/3038228.3038234
  * [69] R. Featherstone, _Rigid Body Dynamics Algorithms_. Springer.
  * [70] A. Wächter and L. T. Biegler, “On the Implementation of an Interior-point Filter Line-search Algorithm for Large-scale Nonlinear Programming,” vol. 106, no. 1, pp. 25–57. [Online]. Available: https://doi.org/10.1007/s10107-004-0559-y
  * [71] P. E. Gill, W. Murray, and M. A. Saunders, “SNOPT: An SQP Algorithm for Large-scale Constrained Optimization,” vol. 47, no. 1, pp. 99–131.
  * [72] M. Schubiger, G. Banjac, and J. Lygeros, “GPU Acceleration of ADMM for Large-Scale Quadratic Programming,” vol. 144, pp. 55–67. [Online]. Available: http://arxiv.org/abs/1912.04263
  * [73] J. Nickolls, I. Buck, M. Garland, and K. Skadron, “Scalable Parallel Programming with CUDA,” vol. 6, pp. 40–53.
  * [74] NVIDIA, _NVIDIA CUDA C Programming Guide_ , version 9.1 ed. [Online]. Available: http://docs.nvidia.com/cuda/pdf/CUDA˙C˙Programming˙Guide.pdf
  * [75] X. Bu and B. Plancher, “Symmetric stair preconditioning of linear systems for parallel trajectory optimization,” 2023.
  * [76] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: an operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020. [Online]. Available: https://doi.org/10.1007/s12532-020-00179-2
  * [77] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, J. Carpentier, L. Righetti, S. Vijayakumar, and N. Mansard, “Crocoddyl: An efficient and versatile framework for multi-contact optimal control,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 2536–2542.
  * [78] S. Kleff, A. Meduri, R. Budhiraja, N. Mansard, and L. Righetti, “High-frequency nonlinear model predictive control of a manipulator,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 7330–7336.
  * [79] T. Howell, B. Jackson, and Z. Manchester, “Altro: A fast solver for constrained trajectory optimization,” in _Proceedings of (IROS) IEEE/RSJ International Conference on Intelligent Robots and Systems_ , November 2019, pp. 7674 – 7679.
  * [80] M. Ditty, “Nvidia orin system-on-chip,” in _2022 IEEE Hot Chips 34 Symposium (HCS)_. IEEE Computer Society, 2022, pp. 1–17.