# Modality-invariant Visual Odometry for Embodied Vision

Marius Memmel1 Roman Bachmann2 Amir Zamir2
1University of Washington 2Swiss Federal Institute of Technology (EPFL)
https://vo-transformer.github.io Work done on exchange at EPFL

###### Abstract

Effectively localizing an agent in a realistic, noisy setting is crucial for
many embodied vision tasks. Visual Odometry (VO) is a practical substitute for
unreliable GPS and compass sensors, especially in indoor environments. While
SLAM-based methods show a solid performance without large data requirements,
they are less flexible and robust w.r.t. to noise and changes in the sensor
suite compared to learning-based approaches. Recent deep VO models, however,
limit themselves to a fixed set of input modalities, _e.g_., RGB and depth,
while training on millions of samples. When sensors fail, sensor suites
change, or modalities are intentionally looped out due to available resources,
_e.g_., power consumption, the models fail catastrophically. Furthermore,
training these models from scratch is even more expensive without simulator
access or suitable existing models that can be fine-tuned. While such
scenarios get mostly ignored in simulation, they commonly hinder a model’s
reusability in real-world applications. We propose a Transformer-based
modality-invariant VO approach that can deal with diverse or changing sensor
suites of navigation agents. Our model outperforms previous methods while
training on only a fraction of the data. We hope this method opens the door to
a broader range of real-world applications that can benefit from flexible and
learned VO models.

##  1 Introduction

Artificial intelligence has found its way into many commercial products that
provide helpful digital services. To increase its impact beyond the digital
world, personal robotics and embodied AI aims to put intelligent programs into
bodies that can move in the real world or interact with it [15]. One of the
most fundamental skills embodied agents must learn is to effectively traverse
the environment around them, allowing them to move past stationary
manipulation tasks and provide services in multiple locations instead [40].
The ability of an agent to locate itself in an environment is vital to
navigating it successfully [12, 64]. A common setup is to equip an agent with
an RGB-D (RGB and 0pt) camera and a GPS+Compass sensor and teach it to
navigate to goals in unseen environments [2]. With extended data access
through simulators [40, 47, 57, 39, 28], photo-realistic scans of 3D
environments [28, 7, 46, 58, 56], and large-scale parallel training, recent
approaches reach almost perfect navigation results in indoor environments
[55]. However, these agents fail catastrophically in more realistic settings
with noisy, partially unavailable, or failing RGB-D sensor readings, noisy
actuation, or no access to GPS+Compass [6, 64].

Figure 1: An agent is tasked to navigate to a goal location using RGB-D
sensors. Because GPS+Compass are not available, the location is inferred from
visual observations only. Nevertheless, sensors can malfunction, or
availability can change during test-time (indicated by $\sim$), resulting in
catastrophic failure of the localization. We train our model to react to such
scenarios by randomly dropping input modalities. Furthermore, our method can
be extended to learn from multiple arbitrary input modalities, _e.g_., surface
normals, point clouds, or internal measurements.

Visual Odometry (VO) is one way to close this performance gap and localize the
agent from only RGB-D observations [2], and deploying such a model has been
shown to be especially beneficial when observations are noisy [12, 64].
However, those methods are not robust to any sensory changes at the test-time,
such as a sensor failing, underperforming, or being intentionally looped out.
In practical applications [43], low-cost hardware can also experience serious
bandwidth limitations, causing RGB (3 channels) and 0pt(1 channel) to be
transferred at different rates. Furthermore, mobile edge devices must balance
battery usage by switching between passive (_e.g_., RGB) and active (_e.g_.,
LIDAR) sensors depending on the specific episode. Attempting to solve this
asymmetry by keeping separate models in memory, relying on active sensors, or
using only the highest rate modality is simply infeasible for high-speed and
real-world systems. Finally, a changing sensor suite represents an extreme
case of sensor failure where access to a modality is lost during test-time.
These points demonstrate the usefulness of a certain level of modality
invariance in a VO framework. Those scenarios decrease the robustness of SLAM-
based approaches [32] and limit the transferability of models trained on RGB-D
to systems with only a subset or different sensors.

We introduce “optional” modalities as an umbrella term to describe settings
where input modalities may be of limited availability at test-time. Figure 1
visualizes a typical indoor navigation pipeline, but introduces uncertainty
about modality availability (_i.e_. at test-time, only a subset of all
modalities might be available). While previous approaches completely neglect
such scenarios, we argue that explicitly accounting for “optional” modalities
already _during training_ of VO models allows for better reusability on
platforms with different sensor suites and trading-off costly or unreliable
sensors during test-time. Recent methods [12, 64] use Convolution Neural
Network (ConvNet) architectures that assume a constant channel size of the
input, which makes it hard to deal with multiple ”optional” modalities. In
contrast, Transformers [51] are much more amenable to variable-sized inputs,
facilitating the training of models that can optionally accept one or multiple
modalities [4].

Transformers are known to require large amounts of data for training from
scratch. Our model’s data requirements are significantly reduced by
incorporating various biases: We utilize multi-modal pre-training [30, 17, 4],
which not only provides better initializations but also improves performance
when only a subset of modalities are accessible during test-time [4].
Additionally, we propose a token-based action prior. The action taken by the
agent has shown to be beneficial for learning VO [64, 35] and primes the model
towards the task-relevant image regions.

We introduce the Visual Odometry Transformer (VOT), a novel modality-agnostic
framework for VO based on the Transformer architecture. Multi-modal pre-
training and an action prior drastically reduce the data required to train the
architecture. Furthermore, we propose explicit modality-invariance training.
By dropping modalities during training, a single VOT matches the performance
of separate uni-modal approaches. This allows for traversing different sensors
during test-time and maintaining performance in the absence of some training
modalities.

We evaluate our method on point-goal navigation in the Habitat Challenge 2021
[1] and show that VOT outperforms previous methods [35] with training on only
5% of the data. Beyond this simple demonstration, we stress that our framework
is modality-agnostic and not limited to RGB-D input or discrete action spaces
and can be adapted to various modalities, _e.g_., point clouds, surface
normals, gyroscopes, accelerators, compass, etc. To the best of our knowledge,
VOT is the first widely applicable modality-invariant Transformer-based VO
approach and opens up exciting new applications of deep VO in both simulated
and real-world applications. We make our code available at
github.com/memmelma/VO-Transformer.

Figure 2: The Visual Odometry Transformer architecture for RGB-D input. Image
patches are turned into tokens through modality-specific linear projections
before a fixed positional embedding is added to them. We pass an action token
that embeds the action  taken by the agent as we find it acts as a strong
prior on the VO problem. An MLP-head  then estimates the VO parameters
$\widehat{\beta},\widehat{\bm{\xi}}$, i.e., translation and rotation of the
agent, from the output token. By randomly dropping either RGB or 0ptduring
training, the Transformer backbone  becomes modality-agnostic, allowing it to
deal with a subset of these input modalities during test-time without losing
performance. When more modalities are available during training, other
modality-specific linear projections can be added to process the additional
information.

##  2 Related Work

SLAM- vs Learning-based Navigation: Simultaneous Localization and Mapping
(SLAM) approaches decompose the navigation task into the components of
mapping, localization, planning, and control [49]. These methods rely on
explicit visual feature extraction and, therefore, fail in realistic settings
with noisy observations [64], while learning-based methods are more robust to
noise, ambiguous observations, and limited sensor suites [27, 32]. However,
learning-based methods require an order of magnitude more data, _e.g_.,
available through simulation [40]. To deal with the large data requirements,
SLAM\- and learning-based methods can be combined [8, 9, 5, 63, 61, 48, 11].

Visual Odometry for Realistic Indoor Navigation: While most VO methods
estimate an agent’s pose change from more than two frames [52, 53] or optical
flow [66], subsequent frames in indoor environments share almost no overlap
and contain many occlusions due to the large displacement caused by the
discrete action space [64]. Datta _et al_. [12] propose to estimate the pose
change from consecutive frames via a ConvNet architecture and decouple
learning the VO from the task-specific navigation policy to allow for
retraining modules when dynamics change or the actuation experiences noise.
Zhao _et al_. [64] improve the model’s robustness to observation and actuation
noise through geometric invariance losses [54], separate models for moving and
turning, pre-process observations, and introduce dropout [44]. Finally,
Partsey _et al_. [35] explore the need for explicit map building in autonomous
indoor navigation. They apply train- and test-time augmentations and
concatenate an action embedding similar to Zhao _et al_. [64] to the extracted
visual features. A trend is to exploit simulators to gather large datasets (1M
[64], 5M [35]). While this is a reasonable progression, it is infeasible to
re-train the VO model whenever dynamics or sensor configurations change.

Multi-modal Representation Learning: The availability of multi-modal or
pseudo-labeled [4] data [13, 34, 59, 65, 16, 38], _e.g_., depth, video, and
audio, makes it possible to learn feature-rich representations over multiple
modalities. Together with Transformer’s [51] ability to process a token
sequence of arbitrary length, this leads to general-purpose architectures that
can handle various modalities [23] like video, images, and audio [30] or
single-view 3D geometry [17]. In particular, Multi-modal Multi-task Masked
Autoencoder (MultiMAE) [4] is a multi-modal pre-training strategy that
performs masked autoencoding [19] with RGB, 0pt, and Semantic Segmenation
(SemSeg). We show that fine-tuning a pre-trained MultiMAE model can
significantly increase VO performance using only 5% of the training data
amount of previous methods [35].

##  3 Proposed Method

###  3.1 Preliminaries

In the realistic PointGoal Navigation task [2], an agent spawns at a random
position in an unseen environment and is given a random goal location
$\bm{g}_{t}$ relative to its starting position. At each time step $t$ of an
episode, the agent perceives its environment through observations $\bm{o}_{t}$
and executes an action $a_{t}$ from a set of discrete actions (move fwd
$0.25m$, turn left and right by $30\degree$). The stop action indicates the
agent’s confidence in having reached the goal. Because the relative goal
position $\bm{g}_{t}$ is defined at the beginning of each episode, it has to
be updated throughout the episode as the actions change the agent’s position
and orientation. Following [12, 64], we update $\bm{g}_{t}$ through an
estimate of the agent’s coordinate transformation. With access to GPS+Compass,
computing this transformation is trivial. However, since those sensors are
unavailable, we estimate the transformation from the agent’s subsequent
observations $\bm{o}_{t},\bm{o}_{t+1}$ and update the estimated relative goal
position $\widehat{\bm{g}}_{t}$. When taking an action $a_{t}$, the agent’s
coordinate system $C_{t}$ transforms into $C_{t+1}$. Because the agent can
only navigate planarly in the indoor scene, we discard the 3rd dimension for
simplicity. We define the estimated transformation as $\bm{\widehat{H}}\in
SE(2)$, with $SE(2)$ being the group of rigid transformations in a 2D plane
and parameterize it by the estimated rotation angle
$\widehat{\beta}\in\mathbb{R}$ and estimated translation vector
$\bm{\widehat{\xi}}\in\mathbb{R}^{2}$:

$\displaystyle\bm{\widehat{H}}=\begin{bmatrix}\widehat{R}&\bm{\widehat{\xi}}\\\
0&1\end{bmatrix},\quad\widehat{R}=\begin{bmatrix}[r]\cos(\widehat{\beta})&-\sin(\widehat{\beta})\\\
\sin(\widehat{\beta})&\cos(\widehat{\beta})\end{bmatrix}\in SO(2).$

(1)

We then learn a VO model $f_{\phi}$ with parameters $\phi$ predicting
$\widehat{\beta},\bm{\widehat{\xi}}$ from observations
$\bm{o}_{t},\bm{o}_{t+1}$:
$\widehat{\beta},\bm{\widehat{\xi}}=f_{\phi}(\bm{o}_{t},\bm{o}_{t+1})$
Finally, we transform $\widehat{\bm{g}}_{t}$ in coordinate system $C_{t}$ to
the new agent coordinate system $C_{t+1}$ by
$\widehat{\bm{g}}_{t+1}=\bm{\widehat{H}}\cdot\widehat{\bm{g}}_{t}$.

###  3.2 Visual Odometry Transformer

Model Architecture:  When facing “optional” modalities, it is not yet clear
how systems should react. Options range from constructing an alternative
input, _e.g_., noise [29], to falling back on a model trained without the
missing modalities, to training the network with placeholder inputs [31].
Besides these, recent approaches depend on a fixed set of modalities during
train- and test-time due to their ConvNet-based backbone. Transformer-based
architectures can process a variable number of input tokens and can be
explicitly trained to accept fewer modalities during test-time while observing
multiple modalities throughout training [51, 4]. Furthermore, the
Transformer’s global receptive field could be beneficial for VO, which often
gets solved with correspondence or feature matching techniques [41]. We,
therefore, propose the Visual Odometry Transformer (VOT), a multi-modal
Transformer-based architecture for VO.

Visual Odometry Estimation: To estimate the VO parameters, we pass the encoded
Action Token ($[ACT]$) token to a prediction head. We use a two-layer Multi-
layer Perceptron (MLP) with learnable parameters $\psi$ composed into
$\bm{W_{0}}\in\mathbb{R}^{d\times d_{h}},\bm{b}_{0}\in\mathbb{R}^{d_{h}}$, and
$\bm{W_{1}}\in\mathbb{R}^{d_{h}\times 3},\bm{b}_{1}\in\mathbb{R}^{3}$ with
token dimensions $d=768$, and hidden dimensions $d_{h}=d/2$. A Gaussian Error
Linear Unit (GELU) [21] acts as the non-linearity between the two layers. The
VO model can then be defined as a function
$f_{\phi,\psi}(\bm{o}_{t},\bm{o}_{t+1},a_{t})$ taking as input the action
$a_{t}$ and the observations $\bm{o}_{t},\bm{o}_{t+1}$ corresponding to either
RGB, 0pt, or RGB-D and predicting the VO parameters
$\widehat{\beta},\bm{\widehat{\xi}}$. Simplifying the backbone as
$b_{\phi}(\bm{o}_{t},\bm{o}_{t+1},a_{t})$ that returns extracted visual
features $\bm{v}_{t\rightarrow t+1}\in\mathbb{R}^{1\times d}$, and governed by
parameters $\phi$, the resulting model is:

$\begin{split}b_{\phi}(\bm{o}_{t},\bm{o}_{t+1},a_{t})&=\bm{v}_{t\rightarrow
t+1}\\\
\mathrm{MLP}_{\psi}(\bm{v})&=\mathrm{GELU}(\bm{v}\bm{W_{0}}+\bm{b}_{0})\bm{W_{1}}+\bm{b}_{1}\\\
f_{\phi,\psi}(\bm{o}_{t},\bm{o}_{t+1},a_{t})&=\mathrm{MLP}_{\psi}(b_{\phi}(\bm{o}_{t},\bm{o}_{t+1},a_{t}))=\widehat{\beta},\bm{\widehat{\xi}}\\\
\end{split}$ (2)

Action Prior: The action $a_{t}$ taken by the agent to get from $\bm{o}_{t}$
to $\bm{o}_{t+1}$ is a powerful prior on the VO parameters. To provide this
information to the model, we embed the action using an embedding layer [36].
This layer acts as a learnable lookup for each action, mapping it to a fixed-
size embedding. With the embedding size equal to the token dimensions, we can
create an $[ACT]$ and pass the information directly to the model (_cf_. Figure
2). In contrast to [64, 35], we pass the token directly to the encoder instead
of concatenating it to the extracted features. This practice conditions the
visual feature extraction on the action and helps ignore irrelevant parts of
the image. Note that this approach is not limited to discrete actions but
tokens could represent continuous sensor readings like accelerometers,
gyroscopes, and compasses, allowing for flexible deployment, _e.g_., in
smartphones or autonomous vehicles [43].

Explicit Modality-invariance Training: Explicitly training the model to be
invariant to its input modalities is one way of dealing with missing sensory
information during test-time. To enforce this property, we drop modalities
during training to simulate missing modalities during test-time. Furthermore,
this procedure can improve training on less informative modalities by
bootstrapping model performance with more informative ones. For example, RGB
is more prone to overfitting than 0ptbecause the model can latch onto spurious
image statistics, _e.g_. textures. Training on RGB-only would likely cause the
model to latch onto those and converge to local minima, not generalizing well
to unseen scenes. By increasing the amount of 0ptobservations seen during
training, the model learns to relate both modalities, acting as
regularization. We model this notion as a multinomial distribution over
modality combinations (here: RGB, 0pt, RGB-D) with equal probability. For each
batch, we draw a sample from the distribution to determine on which
combination to train.

##  4 Experimental Evaluation

###  4.1 Setup

Simulation: We use the AI Habitat (Habitat) simulator for data collection and
model evaluation, following the Habitat PointNav Challenge 2020 [1]
specifications. The guidelines define an action space of fwd (move forward
$0.25m$), left (turn left by $30\degree$), right (turn right by $30\degree$),
and stop (indicate the agent reached its goal), and include a sensor suite of
RGB-D camera, and GPS+Compass (not used in the realistic PointGoal navigation
task). The RGB observations get returned into a $[0,255]$ range while the
0ptmap is scaled to $[0,10]$. Both sensors are subject to noise, _i.e_., noisy
actuations [33] and observations [10]. Furthermore, collision dynamics prevent
sliding, a behavior that allows the agent to slide along walls on collision.
Cosmetic changes bring the simulation closer to the LoCoBot [18], a low-cost
robotic platform with an agent radius of $0.18m$ and height of $0.88m$. An
optical sensor resolution of $341\times 192$ (width $\times$ height) emulates
an Azure Kinect camera. An episode is successful if the agent calls stop in a
radius two times its own, _i.e_., $0.36m$, around the point goal and does so
in $T=500$ total number of time steps. By specification, the 3D scenes loaded
into Habitat are from the Gibson [57] dataset, more precisely Gibson-4+ [40],
a subset of 72 scenes with the highest quality. The validation set contains 14
scenes, which are not part of the training set.

Dataset: For training VOT, we collect a training- and a validation dataset.
Each set consists of samples containing the ground truth translation
$\bm{\xi}$ and rotation parameters $\beta$ retrieved from a perfect
GPS+Compass sensor, observations $\bm{o_{t}},\bm{o}_{t+1}$, and taken action
$a_{t}$. We keep samples where the agent collides with its environment as the
transformations strongly differ from standard behavior [64]. The collection
procedure follows Zhao _et al_. [64] and is performed as: 1) initialize the
Habitat simulator and load a scene from the dataset, 2) place the agent at a
random location within the environment with a random orientation, 3) sample a
navigable PointGoal the agent should navigate to, 4) compute the shortest path
and let the agent follow it, and 5) randomly sample data points along the
trajectory. We collect $250\text{\,}\mathrm{k}{}$ observation-transformation
pairs from the training and $25\text{\,}\mathrm{k}{}$ from the validation
scenes of Gibson-4+, which is significantly less than comparable methods
($1\text{\,}\mathrm{M}{}$ [64], $5\text{\,}\mathrm{M}{}$ [35]). Furthermore,
we apply data augmentation during training to the left and right actions by
horizontally flipping the observations and computing the inverse
transformation.

Drop: –

Drop: RGB

Drop: 0pt

VOT-B

VOT-B w/ inv.

Figure 3:  Top-down map of the agent navigating the Cantwell scene [58] from
start ( ) to goal ( ). The plot shows the shortest path ( ), the path taken by
the agent ( ), and the ”imaginary” path the agent took, _i.e_., its VO
estimate ( ). We evaluate the model without RGB or 0pt(Drop) to determine
performance when modalities are missing. As expected, the VOT relies heavily
on both modalities, causing the estimation to drift when either RGB or 0ptis
unavailable (top row). The localization error accumulates over the course of
the trajectory and causes the true and imaginary path to diverge, resulting in
failure to complete the episodes. Training a VOT to be modality-invariant (VOT
w/ inv.) removes those reliances and leads to success even when modalities are
missing (bottom row).

Loss Function: Our loss function is the $L2$-norm between the ground truth VO
parameters and their estimated counterparts. We further add the geometric
invariance losses $\mathcal{L}_{inv}$ proposed by Zhao _et al_. [64] and use
the Adam [26] optimizer ($\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=1e^{-8}$) to
minimize the resulting loss function
$\mathcal{L}=\|\bm{\xi}-\bm{\widehat{\xi}}\|^{2}_{2}+\|\beta-\widehat{\beta}\|^{2}_{2}+\mathcal{L}_{inv}$.

Method | Drop | $S\uparrow$ | SPL$\uparrow$ | SSPL$\uparrow$ | $d_{g}\downarrow$
---|---|---|---|---|---
VOT RGB | – | 59.3 | 45.4 | 66.7 | 66.2
VOT 0pt | – | 93.3 | 71.7 | 72.0 | 38.0
[12] | – | 64.5 | 48.9 | 65.4 | 85.3
VOT | – | 88.2 | 67.9 | 71.3 | 42.1
VOT w/ inv. | – | 92.6 | 70.6 | 71.3 | 40.7
[12] | RGB | 0.0 | 0.0 | 5.4 | 398.7
VOT | RGB | 75.9 | 58.5 | 69.9 | 59.5
VOT w/ inv. | RGB | 91.0 | 69.4 | 71.2 | 37.0
[12] | 0pt | 0.0 | 0.0 | 5.4 | 398.7
VOT | 0pt | 26.1 | 20.0 | 58.7 | 148.1
VOT w/ inv. | 0pt | 60.9 | 47.2 | 67.7 | 72.1
Table 1: Results for dropping modalities during test-time. Training a VOT to
be modality-invariant (w/ inv.) leads to no performance drop in comparison to
a VOT trained on a single modality (VOT RGB, VOT 0pt). This shows that a
single VOT w/ inv. can replace multiple modality-dependent counterparts.
Previous approaches [12, 64, 35] become inapplicable, converging to a Blind
behavior. Metrics reported as $e^{-2}$. Bold indicates best results.

Pre-training: Pre-training is a well-known practice to deal with the large
data requirements of Vision Transformers (ViTs) [14, 60], especially in a VO
setting where data is scarce [14, 45, 25]. We use the pre-trained MultiMAE
(RGB \+ 0pt+ SemSeg) made publicly available by Bachmann _et al_. [3]. Since
SemSeg is unavailable in our setting, we discard the corresponding projection
layers.

Training Details: We follow prior work [12, 64, 35] and train our navigation
policy and VO model separately before jointly evaluating them on the
validation set. In contrast to [12, 64], we do not fine-tune the navigation
policy on the trained VO models as it has shown minimal navigation performance
gains in [64] and was abandoned in [35].

We train all models, including baselines, for 100 epochs with 10 warm-up
epochs that increase the learning rate linearly from $0.0$ to $2e^{-4}$, and
evaluate the checkpoints with the lowest validation error. We further find
gradient norm clipping [62] (max gradient norm of $1.0$) to stabilize the
training of VOT but to hurt the performance of the ConvNet baselines. The
training was done with a batch size of 128 on an NVIDIA V100-SXM4-40GB GPU
with automatic mixed-precision enabled in PyTorch [36] to reduce memory
footprint and speed up training. Our backbone is a ViT-B [14] with a patch
size of $16\times 16$ and 12 encoder blocks with 12 Multi-head Attention (MHA)
heads each, and token dimensions 768. To encode the input into tokens, we use
a 2D sine-cosine positional embedding and separate linear projection layers
for each modality. Note that if additional modalities are available, our model
can be extended by adding additional linear input projections or fine-tuning
existing ones [4]. Finally, we pass all available tokens to the model and
resize each observation to $160\times 80\times c$ (width $\times$ height
$\times$ channels c) and concatenate modalities along their height to
$160\times 160\times c$ to reduce computation. We keep a running mean and
variance to normalize RGB and 0ptto zero mean and unit variance.

Figure 4: Absolute difference between ground truth translations
$\bm{\xi_{x}},\bm{\xi_{z}}$ and rotation angle $\bm{\beta}$ to their estimated
counterparts $\bm{\widehat{\cdot}}$. We compare Zhao et al. [64] (Table 2, 2)
to the VOT (Table 2, 13). Our model estimates fwd translation along the
$z$-axis (middle), left, right along $z$-, $x$-axis (left, middle), and the
turning angle $\beta$ (right) more accurately than the baseline. We
successfully capture the displacements caused by the noisy actuation with an
average error (over both axis $x$, $z$) of
$0.25\text{\,}\mathrm{c}\mathrm{m}{}$ (fwd),
$0.7\text{\,}\mathrm{c}\mathrm{m}{}$ (right), and
$0.65\text{\,}\mathrm{c}\mathrm{m}{}$ (left).

Evaluation Metrics: Anderson _et al_. [2] propose the Success weighted by
(normalized inverse) Path Length (SPL) to evaluate agents in a PointGoal or
ObjectGoal navigation setting. A crucial component of this metric is the
success of an episode (success $S=1$, failure $S=0$). With $l$ the shortest
path distance from the starting position and $p$ the length of the path taken
by the agent, the SPL over $N$ episodes is defined as
$\text{SPL}=\frac{1}{N}\sum_{i=0}^{N-1}S^{(i)}\frac{l^{(i)}}{\max(p^{(i)},l^{(i)})}$.
While SPL depends on the success of an episode, [12] propose the Soft Success
Path Length (SSPL) that provides a more holistic view of the agent’s
navigation performance. The authors replace the binary success $S$ of an
episode with a soft value consisting of the ratio between the (geodesic)
distances to target upon start $d_{init}$ and termination of an episode
$d_{g}$. The resulting metric is then
$\text{SSPL}=\frac{1}{N}\sum_{i=0}^{N-1}\left(1-d_{g}^{(i)}/d_{init}^{(i)})\right)\frac{l^{(i)}}{\max\left(p^{(i)},l^{(i)}\right)}$
The closer the agent gets to the goal, the higher the SSPL, even if the
episode is unsuccessful. This softening allows distinguishing agents that fail
to complete a single or multiple episodes but move significantly close to the
goal from ones that move away from it. Without access to GPS+Compass, SSPL
becomes significantly more important as an agent might call stop prematurely
due to inaccurate localization. We report the SPL, SSPL, success $S$, and
(geodesic) distance to goal on termination $d_{g}$ on the validation scenes of
Gibson-4+ with decimals truncated.

Navigation Policy: Similar to prior work [12, 64, 35], we replace the
GPS+Compass with our VO model to estimate the relative goal position, which
serves as the input to a pre-trained navigation policy. We use the same pre-
trained policy as Zhao _et al_. [64] for our experiments, which was trained
using a goal position updated by ground truth localization. The policy
architecture consists of a Long Short-Term Memory (LSTM) [22] with two
recurrent layers that process 1) a 512-dimensional encoding of the agent’s
observations $\bm{o_{t}}$ (here: 0pt), 2) a 32-dimensional embedding of the
previous action, and 3) a 32-dimensional embedding of the updated relative
goal position. The observation encoding gets obtained by passing the
observations $\bm{o}_{t}$ through a ResNet-18 [20] backbone, flattening the
resulting feature map to dimensionality 2052, and projecting it to
dimensionality 512 with a fully-connected layer. Finally, the output of the
LSTM is fed through another fully-connected layer to produce a distribution
over the action space and a value function estimate. The policy was trained
using DDPO [55], a distributed version of Proximal Policy Optimization (PPO)
[42].

| Method | Observations | Pre-train | $[ACT]$ | $S\uparrow$ | SPL$\uparrow$ | SSPL$\uparrow$ | $d_{g}\downarrow$ | $\mathcal{L}_{train}\downarrow$ | $\mathcal{L}_{val}\downarrow$
---|---|---|---|---|---|---|---|---|---|---
1 |  [64] (separate) | RGB-D |  |  | 22.4 | 13.8 | 31.5 | 305.3 | 0.125 | 0.186
2 |  [64] (unified) | RGB-D |  | ✔ | 64.5 | 48.9 | 65.4 | 85.3 | 0.264 | 0.420
3 | Blind | – |  |  | 0.0 | 0.0 | 5.4 | 398.7 | 48.770 | 47.258
4 | VOT-B | RGB |  |  | 27.1 | 21.2 | 57.7 | 177.0 | 0.735 | 1.075
5 | VOT-B | 0pt |  |  | 43.2 | 32.0 | 59.3 | 122.5 | 0.441 | 0.644
6 | VOT-B | RGB-D |  |  | 47.3 | 36.3 | 61.2 | 119.7 | 1.256 | 1.698
7 | Blind | – |  | ✔ | 13.3 | 10.0 | 46.3 | 251.8 | 1.637 | 1.641
8 | VOT-B | RGB |  | ✔ | 42.0 | 32.3 | 62.7 | 107.0 | 0.043 | 0.571
9 | VOT-B | 0pt |  | ✔ | 76.1 | 58.8 | 69.2 | 60.7 | 0.017 | 0.113
10 | VOT-B | RGB-D |  | ✔ | 72.1 | 55.6 | 68.5 | 64.4 | 0.019 | 0.129
11 | VOT-B | RGB | ✔ |  | 54.5 | 41.3 | 65.2 | 69.9 | 0.056 | 0.347
12 | VOT-B | 0pt | ✔ |  | 83.2 | 63.4 | 69.1 | 49.9 | 0.079 | 0.205
13 | VOT-B | RGB-D | ✔ |  | 85.7 | 65.7 | 69.7 | 56.1 | 0.021 | 0.060
14 | VOT-B | RGB | ✔ | ✔ | 59.3 | 45.4 | 66.7 | 66.2 | 0.003 | 0.280
15 | VOT-B | 0pt | ✔ | ✔ | 93.3 | 71.7 | 72.0 | 38.0 | 0.004 | 0.044
16 | VOT-B | RGB-D | ✔ | ✔ | 88.2 | 67.9 | 71.3 | 42.1 | 0.004 | 0.051
17 | VOT-B w/ inv. | RGB-D | ✔ | ✔ | 92.6 | 70.6 | 71.3 | 40.7 | 0.008 | 0.094
| oracle | GPS+Compass | – | – | 97.8 | 74.8 | 73.1 | 29.9 | – | –

Table 2: Ablation study of architecture design and input modalities. We
further investigate pre-training with MultiMAE [4] in models 11-14. Losses
$\mathcal{L}$, Success $S$, SPL, SSPL, and $d_{g}$ reported as $e^{-2}$. Bold
indicates best results.

###  4.2 Dealing With Optional Modalities

We evaluate the models’ robustness to missing modalities by randomly dropping
access to one of the training modalities. This setup probes VOT for
dependencies on the input modalities, which directly influence the downstream
performance under limited access. In case of sensor malfunctioning, _e.g_.,
only a single modality might be available, a ConvNet’s failure is
predetermined as it requires a fixed-size input. If not given, the system
converges to a Blind behavior, exemplified in Table 1. Limiting access to
modalities reveals VOT’s dependency on 0pt. Dropping RGB barely decreases
performance, while dropping 0ptcauses the localization to fail more
drastically. Comparing the true agent localization and its ”imaginary”,
_i.e_., VO estimate, it becomes clear why. Figure 3 shows how the errors
accumulate, causing the true location to drift away from the estimate. While
the effect is less drastic when dropping RGB, the agent still fails to reach
the goal.

Training VOT with the proposed invariance training (w/ inv.), _i.e_., sampling
RGB for 20%, 0ptfor 30%, and RGB-D for 50% of the training batches, eliminates
this shortcoming. Removing RGB now only decreases the success rate by $1.6\%$,
while removing 0ptalso leads to a stronger performance. This observation
suggests that RGB is less informative for the VO task than 0pt. Especially
when navigating narrow passages, RGB might consist of uniform observations,
_e.g_., textureless surfaces like walls, making it hard to infer the
displacement, unlike 0ptwhich would still provide sufficient geometric
information (_cf_. Figure 3). However, this information asymmetry only leads
to a decline in the metrics that are sensitive to subtle inconsistencies in
the localization, _i.e_., $S$, and SPL. Inspecting the SSPL, the drop of
$-3.5$ is less drastic. Explicit modality-invariance training keeps VOT-B
(RGB-D) from exploiting this asymmetry and matches the performance of VOT-B
(RGB) when 0ptis dropped during test-time Tab. 1.

Rank | Participant team | S | SPL | SSPL | $d_{g}$
---|---|---|---|---|---
1 |  MultiModalVO (VOT) (ours) | 93 | 74 | 77 | 21
2 | VO for Realistic PointGoal [35] | 94 | 74 | 76 | 21
3 | inspir.ai robotics | 91 | 70 | 71 | 70
4 | VO2021 [64] | 78 | 59 | 69 | 53
5 | Differentiable SLAM-net [24] | 65 | 47 | 60 | 174

Table 3: Habitat Challenge 2021. Results for the Point Nav Test-Standard Phase
(test-std split) retrieved on 05-Nov-2022.

###  4.3 Quantitative Results

We compare our approach to Zhao _et al_. [64] in terms of downstream
navigation performance, _i.e_., the VO model as GPS+Compass replacement for a
learned navigation agent. We use the same publicly available navigation policy
for both approaches and the published VO models of the baseline [64]. Using
only 25% of the training data, VOT improves performance by $S+12.3$, SPL
$+9.7$, SSPL $+2.0$ (_cf_. Table 2 15) and $S+7.2$, SPL $+5.7$, SSPL $+1.3$
(_cf_. Table 2 16). When training the baseline on our smaller data set (_cf_.
Table 2 2, unified, ResNet-50), this improvement increases to $S+29.8$, SPL
$+22.8$, SSPL $+6.6$ (_cf_. Table 2 15) and $S+23.7$, SPL $+19.0$, SSPL $+5.9$
(_cf_. Table 2 16).

To capture the raw VO performance detached from the indoor navigation task, we
inspect the absolute prediction error in Figure 4. We differentiate between
translation $\bm{\xi}$ in $x$\- and $y$\- direction ($\xi_{x}$, $\xi_{y}$),
and taken action. VOT is accurate up to $0.36\text{\,}\mathrm{c}\mathrm{m}{}$
(fwd), $1.04\text{\,}\mathrm{c}\mathrm{m}{}$ (right),
$1.05\text{\,}\mathrm{c}\mathrm{m}{}$ (left) in $x$\- direction and
$0.20\text{\,}\mathrm{c}\mathrm{m}{}$ (fwd),
$0.41\text{\,}\mathrm{c}\mathrm{m}{}$ (right),
$0.38\text{\,}\mathrm{c}\mathrm{m}{}$ (left) in $z$-direction. Note how the
baseline struggles to capture $\xi_{z}$, corresponding to the forward-moving
direction $z$ when taking the fwd action.

Given the results in Table 2, we advise using VOT trained on 0pt-only when
access is assumed, as the difference to using GPS+Compass is a mere $S-4.5$,
SPL $-3.1$, SSPL $-1.1$. When ”optional” modalities are needed, _e.g_., they
are expected to change during test-time, invariance training should be used.
Trained on RGB-D, this setup also reaches GPS+Compass like performance with
differences of only $S-5.2$, SPL $-4.2$, SSPL $-1.8$.

$o_{t}$ (RGB)

$o_{t+1}$ (RGB)

$o_{t}$ (0pt)

$o_{t+1}$ (0pt)

(a) Ground truth action: left

$o_{t}$ (RGB)

$o_{t+1}$ (RGB)

$o_{t}$ (0pt)

$o_{t+1}$ (0pt)

(b) Ground truth action: fwd

Figure 6:  Attention maps of the last attention layer of VOT (_cf_. Table 2
13). Brighter color indicates higher ( ) and darker color lower ( ) weighting
of the image patch. The VOT learns to focus on regions present in both time
steps $t,t+1$, _i.e_., outer image regions for turning left, and center
regions for moving fwd. Artifacts of the Gibson dataset get ignored (_cf_.
Figure 6(b)).

###  4.4 Ablation Study

We identify the impact of different input modalities and model design choices
in our ablation study (_cf_. Table 2). Without observations, the Blind VO
model cannot update the goal position. This means the agent can only act
without goal-related feedback, resulting in a $0\%$ success rate.

Extending the model with our proposed $[ACT]$ token allows it to surpass the
Blind performance. Able to update the relative goal position, the agent
reaches an SSPL of $46.3$, but due to the actuation noise, it calls stop
correctly only $13.3\%$ of the time. Access to RGB or 0ptallows the VO model
to adjust to those unpredictable displacements. While the RGB and
0ptobservations correlate with the $[ACT]$ token, they also contain
information about the noisy actuation. Vice versa, $[ACT]$ disambiguates
corner cases where the visual observations do not provide explicit information
about the underlying action. For instance, a fwd action colliding with a wall
might be hard to distinguish from a noisy left turning less than $30\degree$
[64].

Our results show that MultiMAE pre-training provides useful multi-modal
features for VO that fine-tuned outperform the ConvNet baselines. In addition,
these features are complementary to the $[ACT]$ prior, together achieving
state-of-the-art results. We conclude that the $[ACT]$ prior biases the model
towards the mean of the corresponding transformation, while the pre-training
supports the learning of the additive actuation noise.

Training separate models for each modality reveals that 0ptis a more
informative modality than RGB for VO. We assume this to be a direct result of
its geometric properties, _i.e_., the 3D structure of the scene. We find that
training VOT on noisy RGB even hurts the localization. The model overfits the
visual appearance of the scenes and is unable to generalize to unseen ones. In
turn, 0ptdoes not suffer from this issue as it only contains geometric
information.

###  4.5 Action-conditioned Feature Extraction

We show what image regions the model attends to by visualizing the attention
maps of the last MHA-layer (_cf_. Table 2 16) corresponding to the $[ACT]$
token in Figure 6. To reduce the dimensionality of the visualization, we fuse
the heads’ weights via the $max$ operator and align the attention maps with
the input images. We normalize the maps to show the full range of the color
scheme.

We find that passing different actions to VOT primes it to attend to
meaningful regions in the image. When passed turning actions left or right,
VOT focuses on regions present at both time steps. This makes intuitive sense,
as a turning action of $30\degree$ strongly displaces visual features or even
pushes them out of the agent’s field of view. A similar behavior emerges for a
fwd action which leads to more attention on the center regions, _e.g_., the
walls and the end of a hallway (_cf_. Figure 6(b)). These results are
particularly interesting as the model has no prior knowledge about the VO task
but learns something about its underlying structure.

###  4.6 Habitat Challenge 2021 PointNav

We compare our approach (_cf_. Table 2 16) to several baselines submitted to
the Habitat Challenge 2021 benchmark in Table 3. Using the same navigation
policy as Partsey _et al_. [35], VOT achieves the highest SSPL and on par SPL
and $d_{g}$ training on only 5% of the data. These results clearly show that
reusability doesn’t come with a price of lower performance and that scaling
data requirements doesn’t seem to be the answer to solving deep VO.

###  4.7 Limitations

In our work, we separate the VO model from the navigation policy and only
focus on the modality-invariance of the former, neglecting that the navigation
policy expects 0ptas input [12, 64, 35]. Designing policies to be modality-
invariant is subject to future research. Assuming an accurate sensor failure
detection when dropping modalities, additionally, is an idealized setup.
Furthermore, our experiments in the Habitat’s simulator limit the available
modalities to RGB-D. Even though SemSeg has shown to be beneficial for some VO
applications [50, 37], there is no specific sensor for it. However, SemSeg
could be estimated from RGB. While our experiments focus on discrete actions
and RGB-D, our architecture could be adapted to continuous actions and other
sensor types. However, training might become more difficult due to a lack of
pre-trained weights.

##  5 Conclusions

We present Visual Odometry Transformers for learned Visual Odometry. Through
multi-modal pre-training and action-conditioned feature extraction, our method
is sample efficient and outperforms current methods trained on an order of
magnitude more data. With its modality-agnostic design and modality-invariance
training, a single model can deal with different sensor suites during training
and can trade-off subsets of those during test-time.

## References

  * [1] Abhishek Kadian*, Joanne Truong*, Aaron Gokaslan, Alexander Clegg, Erik Wijmans, Stefan Lee, Manolis Savva, Sonia Chernova, and Dhruv Batra.  Are We Making Real Progress in Simulated Environments? Measuring the Sim2Real Gap in Embodied Visual Navigation.  In arXiv preprint arXiv:1912.06321, 2019.
  * [2] Peter Anderson, Angel Chang, Devendra Singh Chaplot, Alexey Dosovitskiy, Saurabh Gupta, Vladlen Koltun, Jana Kosecka, Jitendra Malik, Roozbeh Mottaghi, Manolis Savva, et al.  On evaluation of embodied navigation agents.  In arXiv preprint arXiv:1807.06757, 2018.
  * [3] Roman Bachmann and David Mizrahi.  Multimae: Multi-modal multi-task masked autoencoders.  https://github.com/EPFL-VILAB/MultiMAE, 2022.
  * [4] Roman Bachmann, David Mizrahi, Andrei Atanov, and Amir Zamir.  Multimae: Multi-modal multi-task masked autoencoders.  In European Conference on Computer Vision, 2022.
  * [5] Shehroze Bhatti, Alban Desmaison, Ondrej Miksik, Nantas Nardelli, N Siddharth, and Philip HS Torr.  Playing doom with slam-augmented deep reinforcement learning.  In arXiv preprint arXiv:1612.00380, 2016.
  * [6] Lukas Biewald.  Habitat challenge 2020.  https://eval.ai/web/challenges/challenge-page/580/leaderboard/1631, 2020\.
  * [7] Angel Chang, Angela Dai, Thomas Funkhouser, Maciej Halber, Matthias Niessner, Manolis Savva, Shuran Song, Andy Zeng, and Yinda Zhang.  Matterport3d: Learning from rgb-d data in indoor environments.  In International Conference on 3D Vision, 2017.
  * [8] Devendra Singh Chaplot, Dhiraj Gandhi, Saurabh Gupta, Abhinav Gupta, and Ruslan Salakhutdinov.  Learning to explore using active neural slam.  In International Conference on Learning Representations, 2020\.
  * [9] Devendra Singh Chaplot, Ruslan Salakhutdinov, Abhinav Gupta, and Saurabh Gupta.  Neural topological slam for visual navigation.  In Conference on Computer Vision and Pattern Recognition, 2020\.
  * [10] Sungjoon Choi, Qian-Yi Zhou, and Vladlen Koltun.  Robust reconstruction of indoor scenes.  In Conference on Computer Vision and Pattern Recognition, 2015\.
  * [11] Jan Czarnowski, Tristan Laidlow, Ronald Clark, and Andrew J Davison.  Deepfactors: Real-time probabilistic dense monocular slam.  In IEEE Robotics and Automation Letters, 2020.
  * [12] Samyak Datta, Oleksandr Maksymets, Judy Hoffman, Stefan Lee, Dhruv Batra, and Devi Parikh.  Integrating egocentric localization for more realistic point-goal navigation agents.  In Conference on Robot Learning, 2021.
  * [13] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In Conference on Computer Vision and Pattern Recognition, 2009\.
  * [14] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale.  In International Conference on Learning Representations, 2021\.
  * [15] Jiafei Duan, Samson Yu, Hui Li Tan, Hongyuan Zhu, and Cheston Tan.  A survey of embodied ai: From simulators to research tasks.  In IEEE Transactions on Emerging Topics in Computational Intelligence, 2022.
  * [16] Ainaz Eftekhar, Alexander Sax, Jitendra Malik, and Amir Zamir.  Omnidata: A scalable pipeline for making multi-task mid-level vision datasets from 3d scans.  In International Conference on Computer Vision, 2021.
  * [17] Rohit Girdhar, Mannat Singh, Nikhila Ravi, Laurens van der Maaten, Armand Joulin, and Ishan Misra.  Omnivore: A single model for many visual modalities.  In Computer Vision and Pattern Recognition, 2022.
  * [18] Abhinav Gupta, Adithyavairavan Murali, Dhiraj Prakashchand Gandhi, and Lerrel Pinto.  Robot learning in homes: Improving generalization and reducing dataset bias.  In Advances in Neural Information Processing Systems, 2018.
  * [19] Kaiming He, Xinlei Chen, Saining Xie, Yanghao Li, Piotr Dollár, and Ross Girshick.  Masked autoencoders are scalable vision learners.  In Computer Vision and Pattern Recognition, 2022.
  * [20] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Conference on Computer Vision and Pattern Recognition, 2016\.
  * [21] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  In arXiv preprint arXiv:1606.08415, 2016.
  * [22] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  In Neural Computation, 1997.
  * [23] Andrew Jaegle, Sebastian Borgeaud, Jean-Baptiste Alayrac, Carl Doersch, Catalin Ionescu, David Ding, Skanda Koppula, Daniel Zoran, Andrew Brock, Evan Shelhamer, Olivier J Henaff, Matthew Botvinick, Andrew Zisserman, Oriol Vinyals, and Joao Carreira.  Perceiver IO: A general architecture for structured inputs & outputs.  In International Conference on Learning Representations, 2022\.
  * [24] Peter Karkus, Shaojun Cai, and David Hsu.  Differentiable slam-net: Learning particle slam for visual navigation.  In Computer Vision and Pattern Recognition, 2021.
  * [25] Alex Kendall, Matthew Grimes, and Roberto Cipolla.  Posenet: A convolutional network for real-time 6-dof camera relocalization.  In International Conference on Computer Vision, 2015.
  * [26] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In International Conference on Learning Representations, 2015\.
  * [27] Noriyuki Kojima and Jia Deng.  To learn or not to learn: Analyzing the role of learning for navigation in virtual environments.  In arXiv preprint arXiv:1907.11770, 2019.
  * [28] Eric Kolve, Roozbeh Mottaghi, Winson Han, Eli VanderBilt, Luca Weihs, Alvaro Herrasti, Daniel Gordon, Yuke Zhu, Abhinav Gupta, and Ali Farhadi.  Ai2-thor: An interactive 3d environment for visual ai.  In arXiv preprint arXiv:1712.05474, 2017.
  * [29] Yu-Jhe Li, Jinhyung Park, Matthew O’Toole, and Kris Kitani.  Modality-agnostic learning for radar-lidar fusion in vehicle detection.  In Conference on Computer Vision and Pattern Recognition, 2022.
  * [30] Valerii Likhosherstov, Anurag Arnab, Krzysztof Choromanski, Mario Lucic, Yi Tay, Adrian Weller, and Mostafa Dehghani.  Polyvit: Co-training vision transformers on images, videos and audio.  In arXiv preprint arXiv:2111.12993, 2021.
  * [31] Marius Memmel, Camila Gonzalez, and Anirban Mukhopadhyay.  Adversarial Continual Learning for Multi-domain Hippocampal Segmentation.  In Domain Adaptation and Representation Transfer, and Affordable Healthcare and AI for Resource Diverse Global Health, 2021.
  * [32] Dmytro Mishkin, Alexey Dosovitskiy, and Vladlen Koltun.  Benchmarking classic and learned navigation in complex 3d environments.  In arXiv preprint arXiv:1901.10915, 2019.
  * [33] Adithyavairavan Murali, Tao Chen, Kalyan Vasudev Alwala, Dhiraj Gandhi, Lerrel Pinto, Saurabh Gupta, and Abhinav Gupta.  Pyrobot: An open-source robotics framework for research and benchmarking.  In arXiv preprint arXiv:1906.08236, 2019.
  * [34] Pushmeet Kohli Nathan Silberman, Derek Hoiem and Rob Fergus.  Indoor segmentation and support inference from rgbd images.  In European Conference on Computer Vision, 2012.
  * [35] Ruslan Partsey, Erik Wijmans, Naoki Yokoyama, Oles Dobosevych, Dhruv Batra, and Oleksandr Maksymets.  Is mapping necessary for realistic pointgoal navigation?  In Conference on Computer Vision and Pattern Recognition, 2022\.
  * [36] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems, 2019.
  * [37] Noha Radwan, Abhinav Valada, and Wolfram Burgard.  Vlocnet++: Deep multitask learning for semantic visual localization and odometry.  In IEEE Robotics and Automation Letters, 2018.
  * [38] Mike Roberts, Jason Ramapuram, Anurag Ranjan, Atulit Kumar, Miguel Angel Bautista, Nathan Paczan, Russ Webb, and Joshua M. Susskind.  Hypersim: A photorealistic synthetic dataset for holistic indoor scene understanding.  In International Conference on Computer Vision, 2021.
  * [39] Manolis Savva, Angel X Chang, Alexey Dosovitskiy, Thomas Funkhouser, and Vladlen Koltun.  Minos: Multimodal indoor simulator for navigation in complex environments.  In arXiv preprint arXiv:1712.03931, 2017.
  * [40] Manolis Savva, Abhishek Kadian, Oleksandr Maksymets, Yili Zhao, Erik Wijmans, Bhavana Jain, Julian Straub, Jia Liu, Vladlen Koltun, Jitendra Malik, et al.  Habitat: A platform for embodied ai research.  In Conference on Computer Vision and Pattern Recognition, 2019\.
  * [41] Davide Scaramuzza and Friedrich Fraundorfer.  Visual odometry [tutorial].  In IEEE Robotics & Automation Magazine, 2011.
  * [42] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms.  In arXiv preprint arXiv:1707.06347, 2017.
  * [43] Siddhartha S Srinivasa, Patrick Lancaster, Johan Michalove, Matt Schmittle, Colin Summers, Matthew Rockett, Joshua R Smith, Sanjiban Choudhury, Christoforos Mavrogiannis, and Fereshteh Sadeghi.  Mushr: A low-cost, open-source robotic racecar for education and research.  In arXiv preprint arXiv:1908.08031, 2019.
  * [44] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  In The Journal of Machine Learning Research, 2014.
  * [45] Andreas Steiner, Alexander Kolesnikov, Xiaohua Zhai, Ross Wightman, Jakob Uszkoreit, and Lucas Beyer.  How to train your vit? data, augmentation, and regularization in vision transformers.  In arXiv preprint arXiv:2106.10270, 2021.
  * [46] Julian Straub, Thomas Whelan, Lingni Ma, Yufan Chen, Erik Wijmans, Simon Green, Jakob J Engel, Raul Mur-Artal, Carl Ren, Shobhit Verma, et al.  The replica dataset: A digital replica of indoor spaces.  In arXiv preprint arXiv:1906.05797, 2019.
  * [47] Andrew Szot, Alexander Clegg, Eric Undersander, Erik Wijmans, Yili Zhao, John Turner, Noah Maestre, Mustafa Mukadam, Devendra Singh Chaplot, Oleksandr Maksymets, et al.  Habitat 2.0: Training home assistants to rearrange their habitat.  In Advances in Neural Information Processing Systems, 2021.
  * [48] Zachary Teed and Jia Deng.  Droid-slam: Deep visual slam for monocular, stereo, and rgb-d cameras.  In Advances in Neural Information Processing Systems, 2021.
  * [49] Sebastian Thrun.  Probabilistic robotics.  In Communications of the ACM, 2002.
  * [50] Abhinav Valada, Noha Radwan, and Wolfram Burgard.  Deep auxiliary learning for visual localization and odometry.  In IEEE International Conference on Robotics and Automation, 2018\.
  * [51] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is All You Need.  In Advances in Neural Information Processing Systems, 2017.
  * [52] Sen Wang, Ronald Clark, Hongkai Wen, and Niki Trigoni.  Deepvo: Towards end-to-end visual odometry with deep recurrent convolutional neural networks.  In IEEE International Conference on Robotics and Automation, 2017\.
  * [53] Sen Wang, Ronald Clark, Hongkai Wen, and Niki Trigoni.  End-to-end, sequence-to-sequence probabilistic visual odometry through deep neural networks.  In The International Journal of Robotics Research, 2018.
  * [54] Xiangwei Wang, Daniel Maturana, Shichao Yang, Wenshan Wang, Qijun Chen, and Sebastian Scherer.  Improving learning-based ego-motion estimation with homomorphism-based losses and drift correction.  In International Conference on Intelligent Robots and Systems, 2019.
  * [55] Erik Wijmans, Abhishek Kadian, Ari Morcos, Stefan Lee, Irfan Essa, Devi Parikh, Manolis Savva, and Dhruv Batra.  Dd-ppo: Learning near-perfect pointgoal navigators from 2.5 billion frames.  In International Conference on Learning Representations, 2020\.
  * [56] Yi Wu, Yuxin Wu, Georgia Gkioxari, and Yuandong Tian.  Building generalizable agents with a realistic and rich 3d environment.  In International Conference on Learning Workshops, 2018.
  * [57] Fei Xia, William B Shen, Chengshu Li, Priya Kasimbeg, Micael Edmond Tchapmi, Alexander Toshev, Roberto Martín-Martín, and Silvio Savarese.  Interactive gibson benchmark: A benchmark for interactive navigation in cluttered environments.  In IEEE Robotics and Automation Letters, 2020.
  * [58] Fei Xia, Amir R Zamir, Zhiyang He, Alexander Sax, Jitendra Malik, and Silvio Savarese.  Gibson env: Real-world perception for embodied agents.  In Conference on Computer Vision and Pattern Recognition, 2018\.
  * [59] Amir R. Zamir, Alexander Sax, William Shen, Leonidas J. Guibas, Jitendra Malik, and Silvio Savarese.  Taskonomy: Disentangling task transfer learning.  In Conference on Computer Vision and Pattern Recognition, 2018\.
  * [60] Xiaohua Zhai, Alexander Kolesnikov, Neil Houlsby, and Lucas Beyer.  Scaling vision transformers.  In Computer Vision and Pattern Recognition, 2022.
  * [61] Huangying Zhan, Chamara Saroj Weerasekera, Jia-Wang Bian, and Ian Reid.  Visual odometry revisited: What should be learnt?  In IEEE International Conference on Robotics and Automation, 2020\.
  * [62] Jingzhao Zhang, Tianxing He, Suvrit Sra, and Ali Jadbabaie.  Why gradient clipping accelerates training: A theoretical justification for adaptivity.  In International Conference on Learning Representations, 2020\.
  * [63] Jingwei Zhang, Lei Tai, Ming Liu, Joschka Boedecker, and Wolfram Burgard.  Neural slam: Learning to explore with external memory.  In arXiv preprint arXiv:1706.09520, 2017.
  * [64] Xiaoming Zhao, Harsh Agrawal, Dhruv Batra, and Alexander G Schwing.  The surprising effectiveness of visual odometry techniques for embodied pointgoal navigation.  In International Conference on Computer Vision, 2021.
  * [65] Bolei Zhou, Hang Zhao, Xavier Puig, Tete Xiao, Sanja Fidler, Adela Barriuso, and Antonio Torralba.  Semantic understanding of scenes through the ade20k dataset.  In International Journal of Computer Vision, 2019.
  * [66] Ran Zhu, Mingkun Yang, Wang Liu, Rujun Song, Bo Yan, and Zhuoling Xiao.  Deepavo: Efficient pose refining with feature distilling for deep visual odometry.  In Neurocomputing, 2022.


# Parallel two-stage reduction to Hessenberg-triangular form

Thijs Steel  Raf Vandebril

###### Abstract

We present a two-stage algorithm for the parallel reduction of a pencil to
Hessenberg-triangular form. Traditionally, two-stage Hessenberg-triangular
reduction algorithms achieve high performance in the first stage, but struggle
to achieve high performance in the second stage. Our algorithm extends
techniques described by Karlsson et al. [9] to also achieve high performance
in the second stage. Experiments in a shared memory environment demonstrate
that the algorithm can outperform state-of-the-art implementations.

##  1 Introduction

We consider the parallel reduction of a matrix pencil to Hessenberg-triangular
form. For a pencil $(A,B)\,,\,A,B\in\mathbb{C}^{n\times n}$, this consists of
finding unitary matrices $Q$ and $Z$, a Hessenberg matrix $H$ and an upper
triangular matrix $T$ so that

$Q(H,T)Z^{*}=(A,B).$

The most common use for such a decomposition is as a preprocessing step for
the QZ algorithm, which solves generalized eigenvalue problems and was
introduced by C.B. Moler and G.W. Stewart [10]. The goal of this paper is to
introduce and evaluate a new Hessenberg-triangular reduction algorithm that
scales well on multicore machines.

A common way to parallelize dense linear algebra algorithms is to formulate
them so that as many of the flops as possible are within large matrix-matrix
multiplications. Implementations of these algorithms can then rely on highly
optimized implementations of matrix-matrix multiplications that perform well
in parallel. In 2008, B. Kågström , D. Kressner, E.S. Quintana-Ortí, and G.
Quintana-Orti [7] introduced an algorithm for the reduction of a pencil to
Hessenberg-triangular form. Their algorithm performs at least 60% of its
operations via matrix-matrix multiplications. This is a significant
improvement over the original reduction by Moler and Stewart, but it is still
problematic when we consider parallelization. If we rely only on the
parallelization of the matrix-matrix multiplications, then 40% of the work
will not be parallelized.

A possible solution is to use two-stage algorithms. In the first stage, the
pencil is reduced to r-Hessenberg-triangular form. A pencil $(H,T)$ is in
r-Hessenberg-triangular form if $T$ is upper triangular and
$H_{i,j}=0\,\forall i>j+r$, i.e. $H$ has at most $r$ nonzero subdiagonals.
Efficiently parallelizing a reduction to r-Hessenberg-triangular form is much
easier. A parallel algorithm for this reduction was first introduced by K.
Dackland and B. Kågström [6] and improved upon by Kågström et al. [7]. In the
second stage, the pencil in r-Hessenberg-triangular form is then reduced to
Hessenberg-triangular form. To the extent of our knowledge, no effective
parallelizations of the second stage are present in the literature.
Additionally, the combination of the two stages requires approximately 50%
more flops than the one-stage algorithms. Despite these disadvantages,
Kågström et al. report that the two-stage approach can achieve good
performance.

Similar two-stage algorithms exist for the reduction of a matrix to Hessenberg
form. In particular, L. Karlsson and B. Kågström [8, 9] were able to
efficiently parallelize the second stage for Hessenberg matrices. The main
contribution of this paper is to extend their work to Hessenberg-triangular
reductions.

Another solution can be found in the work of B. Adlerborn, L. Karlsson, and B.
Kågström [1]. They also realized that relying on the parallelization of
matrix-matrix multiplications is not an effective strategy. Instead of
switching to a two-stage algorithm, they parallelize each step of the one-
stage reduction of Kågström et al. [7] with great care.

The rest of this paper is organized as follows. First, we recall the blocked
reduction to r-Hessenberg-triangular form by Kågström et al. [7] and discuss
ways to parallelize it in Section 2. Second, we introduce a new parallel
algorithm to reduce r-Hessenberg-triangular pencils to Hessenberg-triangular
form in Section 3. Finally, we evaluate the performance of the parallel two-
stage algorithm in Section 4.

##  2 Reduction to r-Hessenberg-triangular form

###  2.1 WY form of Householder reflectors

In the algorithms we will present later, we will often need to apply a
sequence of Householder reflectors to a matrix. Applying these reflectors one
at a time is usually not optimal. For better cache performance, we can use the
WY representation of this sequence. If $Q=Q_{r}Q_{r-1}\dots Q_{1}$ is the
product of these reflectors, then we can factorize $Q$ as

$Q=I-WY^{*},$

where $Q$ is an $n\times n$ matrix and $W$ and $Y$ are $n\times r$ matrices.
An algorithm to compute the matrices $W$ and $Y$ is due to C. Bischof and C.
Van Loan [4]. Multiplying $Q$ with an $n\times m$ matrix $A$ takes roughly the
same amount of flops as applying the reflectors individually, but it can be
performed using just two matrix-matrix multiplications instead of $r$
applications of a reflector. If a highly optimized matrix-matrix
multiplication is available, using the WY representation is typically much
more efficient.

In the rest of this text, when we refer to something as a block reflector, we
mean that it is a sequence of reflectors that is represented using its WY
representation.

###  2.2 Blocked stage one

Essentially, the blocked stage one algorithm (reduction to r-Hessenberg-
triangular form) is a panel reduction. In each iteration, the algorithm
reduces a panel of $n_{b}$ columns of $A$ using block reflectors from the left
while using block reflectors from the right to preserve the structure of $B$.
The key factor that makes this algorithm a panel reduction is that only
submatrices of $A$ and $B$ are updated while calculating the reflectors. After
these have been calculated, the rest of the updates can be applied efficiently
using matrix-matrix multiplications. Let us consider one iteration of the
algorithm, where we want to reduce columns $j$ through $j+n_{b}-1$. Figure 1
illustrates this iteration. The rest of this subsection explains the details
of the algorithm.

The initial state of the pencil is shown in Figure 1(a). The simplest way to
calculate the reflectors from the left would be to take a QR factorization of
the block $A(j+n_{b}:n,j:j+n_{b}-1)$. The orthogonal factor of this QR
factorization reduces an entire panel in $A$. However, applying this block
reflector would create a large amount of fill-in in $B$. To avoid this fill-
in, we split $A(j+n_{b}:n,j:j+n_{b}-1)$ into several blocks of size
$pn_{b}\times n_{b}$ and take the QR factorization of those blocks, with $p$
an implementation parameter. The reflectors can be calculated while only
updating the panel. Afterward, the reflectors can be accumulated and applied
to the rest of the pencil using matrix-matrix multiplications. The reduced
blocks are illustrated in Figure 1(b). If $p$ is larger, then the block
reflectors are larger and fewer in number, which is typically more efficient.
However, if $p$ is large, it will also lead to a large amount of fill-in.
Kågström et al. [7] report that modest values between 5 and 12 are usually
optimal.

After the multiplications from the left, the pencil is in a state shown in
Figure 1(c). Now that the block column in $A$ has been reduced, we need to
remove the fill-in in $B$ without perturbing the reduced columns of $A$. The
simplest way to reduce $B$ back to upper triangular form would be to calculate
the RQ factorizations of all the subblocks that were filled in and apply their
orthogonal factors from the right. It was noted by Kågström et al. [7] that
fully reducing $B$ back to upper triangular form is not necessary, only the
first $n_{b}$ columns of the subblock need to be reduced.

Reducing a selection of columns using reflectors applied from the right is
nontrivial. Usually, a reflector applied from the right reduces a selection of
rows, not columns. To perform this task, opposite Householder reflectors can
be employed. As first noted by D. Watkins [12], an opposite Householder
reflector can reduce a column when applied from the right. To use this
technique, we take the LQ factorization of the first $n_{b}$ rows of the
orthogonal factor of the RQ factorization of the subblock. If we apply the
orthogonal factor of this LQ factorization to $B$ from the right it will
reduce the first $n_{b}$ columns of the subblock. For a detailed analysis, we
refer to Kågström et al. [7].

The way to generate the opposite reflectors may seem strange. We employ
opposite reflectors to avoid the RQ decomposition of the subblock, but to
generate the opposite reflectors we need to calculate the RQ decomposition
anyway. The advantage becomes clear when we consider the cost of applying
reflectors to the rest of the matrix. The RQ factorization uses $pn_{b}$
reflectors, whereas the LQ factorization only uses $n_{b}$ reflectors. The
full reduction from the right is illustrated in Figure 1(d). Note that these
multiplications do not affect the structure of $A$.

Figure LABEL:sub@subfig:blockedstage1afterrightmul shows the final state of
the pencil in this iteration. The subdiagonal blocks in $B$ are such that they
will not interfere with subsequent iterations. During each iteration, the
blocks in $B$ move down by $n_{b}$ positions eventually leading them off the
edge of the pencil.

Finally, we consider the number of flops required for stage one. If we include
the flops required to update $Q$ and $Z$, Algorithm 1 requires
$\frac{28p+14}{3(p-1)}n^{3}+O(n^{2})$ flops. In our implementation, we have
chosen $p=8$, which results in a cost of $11.33\,n^{3}+O(n^{2})$ flops.

(a) Initial state of the pencil at the start of a panel reduction.

(b) Reducing a panel using multiplications from the left.

(c) The pencil after the multiplications from the left.

(d) Reducing the fill in

(e) After fill in has been reduced

Figure 1: Illustration of the reduction of one panel in stage 1 (Algorithm 1),
with parameters: $n=30$, $n_{b}=4$, $p=3$. The blocks indicated in gray will
be reduced. The updates will affect the parts of the pencil indicated in light
gray (the slightly darker gray shows where these updates overlap).

###  2.3 Parallel stage one

Because it is a blocked algorithm, Algorithm 1 already has some potential for
parallelization. The application of a block reflector to a matrix uses several
matrix-matrix multiplications, which can usually be parallelized well. This is
essentially the parallelization strategy used by Dackland et al. [6]. However,
for typical choices of $r$ and $p$, the dimensions of these matrices are too
small to effectively parallelize them. Here, we present a different
parallelization of stage one.

We will be utilizing two levels of parallelization. First, we will identify
several large-grained tasks. Second, we will split some of the large-grained
tasks into fine-grained tasks.

Algorithm 1 operates in the following way: Calculate a block reflector and
then apply this block reflector to the relevant matrices. For the updates from
the left, an important realization is that the sequence of block reflectors
depends only on the panel that is being reduced. We can introduce more
parallelism by generating the sequence of reflectors as one task (the $G_{L}$
task) and then applying this sequence to $A$, $B$, and $Q$ in parallel (the
$L_{A}$, $L_{B}$ and $L_{Q}$ tasks). For the updates from the right, the
sequence of reflectors cannot be fully calculated without updating $B$. We can
still split the work into the generation of the sequence of reflectors, which
includes applying the sequence to $B$ (the $G_{R}$ task) and the application
of this sequence to $A$ and $Z$ (the $R_{A}$ and $R_{Z}$ tasks). Figure 2
illustrates the dependencies between the tasks.

$G_{L}$$L_{A}$$L_{B}$$L_{Q}$$G_{R}$$R_{A}$$R_{Z}$ Figure 2: Dependency graph
of parallel stage one.

These large-grained tasks allow for some parallelism, but at most three tasks
can be executed in parallel this way. To obtain an efficient algorithm, the
tasks must be split into smaller tasks. The application tasks: $L_{Q}$,
$L_{A}$ and $L_{B}$ apply a sequence of block reflectors to their respective
matrices. To update a certain column of these matrices, the tasks only require
that specific column, so $L_{Q}$, $L_{A}$ and $L_{B}$ can easily be
parallelized by splitting the matrices into several column slices and updating
the slices in parallel. The same trick can be applied to $R_{Z}$ and $R_{A}$
by splitting the matrices into row slices instead of column slices. The
distribution is illustrated in Figure 3. We emphasize that such a distribution
of the matrices is probably not well suited to a distributed memory
implementation, because the cores constantly switch between operating on rows
and columns of the matrices. We also note that $L_{B}$ can have some load-
balancing issues. Because of the upper-triangular structure of $B$, the
different subtasks of $L_{B}$ consist of a different amount of flops. This
could be solved by letting the tasks that affect fewer rows process more
columns to equalize the number of flops. In our implementation, we chose to
let the dynamic scheduler handle these load imbalances.

In a way, this parallelization is similar to the simple parallelization of the
matrix-matrix multiplication mentioned at the start of this section. The
difference is that instead of distributing the application of a single block
reflector, we distribute the application of a sequence of block reflectors.
This results in the same amount of parallelism, but there are fewer
synchronization points. Finally, we consider the parallelization of the
generation tasks: $G_{L}$ and $G_{R}$. $G_{L}$ is a small task and is not
worth parallelizing. $G_{R}$ is the hardest task to parallelize. Calculating
the next block reflector requires applying the previous block reflector to
$B$, so the parallelization we used for $R_{A}$ is not applicable here. Only
the simple parallelization of the matrix-matrix multiplications is possible.

$T_{1}$$T_{2}$$T_{3}$$T_{4}$ (a) $L_{A}$

$T_{1}$$T_{2}$$T_{3}$$T_{4}$ (b) $R_{A}$

$T_{1}$$T_{2}$$T_{3}$$T_{4}$ (c) $L_{B}$

Figure 3: Illustration of the distribution of the $L_{A}$, $R_{A}$ and $L_{B}$
tasks of the parallel phase 1 algorithm to different subtasks
$T_{1},T_{2},\dots$. The dark gray areas indicate the blocks that are reduced.
The light gray areas indicate the parts of the matrix that are affected by the
multiplications.

##  3 Reduction to Hessenberg-triangular form

To reduce a pencil in r-Hessenberg-triangular form to Hessenberg-triangular
form, we adapt the reduction algorithm for r-Hessenberg matrices by Karlsson
et al. [8, 9] to r-Hessenberg-triangular pencils. We start by explaining the
unblocked algorithm and will explain how to delay and efficiently apply
updates later.

###  3.1 Unblocked stage two

Just like in stage one, we will reduce entries in $A$ with reflectors applied
from the left, while eliminating the fill-in in $B$ with reflectors applied
from the right. The algorithm consists of multiple sweeps. During one sweep,
we will reduce one column in $A$ and eliminate the fill-in. At the start of
sweep $j$, the first $j-1$ columns are already in Hessenberg-triangular form.
The first step in a sweep is to apply a reflector $\hat{Q}_{0}^{j}$ from the
left so that $\hat{Q}_{0}^{j}A(j+1:j+r,j)$ is reduced. This will create fill-
in in $B$ so that there is a block in $B(j+1:j+r,j+1:j+r)$. This first step is
illustrated in Figure 4(b). Next, we apply an opposite reflector
$\hat{Z}_{0}^{j}$ so that the first column of
$B(j+1:j+r,j+1:j+r)\hat{Z}_{0}^{j}$ is reduced. This second step is
illustrated in Figure 4(d). Notice that $\hat{Z}_{0}^{j}$ creates some fill-in
in $A$. We can eliminate this fill-in using $\hat{Q}_{1}^{j}$ and eliminate
the resulting fill-in in $B$ using $\hat{Z}_{1}^{j}$. As Figures 4(e) and 4(f)
illustrate, this leads to a bulge-chasing algorithm where eliminating the
fill-in (sometimes called a bulge) leads to more fill-in further down the
matrix. This process is repeated until the bulge is chased off the end of the
matrix. The full process is shown in Algorithm 2.

(a) Initial state of the pencil

(b) Reducing a column in $A$

(c) New state of the pencil

(d) Reducing fill-in in $B$

(e) Chase down fill-in

(f) Chase down fill-in

Figure 4: Illustration of Algorithm 2.

Note that just like in stage one $B$ is block-upper-triangular (and not upper
triangular) after the first sweep. Similarly, $A$ is not in r-Hessenberg-
triangular form. These blocks do not interfere with subsequent sweeps and move
down as more columns are reduced. Traditionally, stage two for Hessenberg-
triangular matrices has been implemented using Givens rotations, a description
of this algorithm can be found in the work of Kågström et al. [7]. Givens
rotations avoid the RQ factorizations required for calculating the opposite
reflectors, but they (asymptotically) require slightly more flops.

If we include the flops required to update $Q$ and $Z$, Algorithm 2 requires
$10n^{3}+O(n^{2})$ flops. Note that the $r^{2}n^{2}$ cost of the RQ
factorizations can be significant if $r$ is large even if it is asymptotically
negligible. In our implementation, we usually choose $r=16$. The original
reduction algorithm of Moler and Stewart requires only $14n^{3}+O(n^{2})$
flops, while the two-stage reduction requires $21.33n^{3}+O(n^{2})$ flops.
That is an increase of more than $40\%$. In the right circumstances,
parallelization will compensate for the increased computational cost.

###  3.2 Blocked stage two

Our blocked stage two algorithm has a simple premise: update as few entries as
possible while generating the reflectors and then use the additional freedom
to reorder the operations to achieve higher efficiency when fully applying the
updates. We will start by explaining how to reorder and apply the updates and
end by showing how to generate the reflectors. These are respectively the
application and generate phases of the blocked algorithm.

#### Applying the reflectors

Assume for the moment that we can somehow generate the reflectors for $q$
consecutive columns without updating the entire matrix. Even in this case,
efficiently applying the reflectors is nontrivial. The unblocked variant
applies the reflectors in the following order

$\hat{Q}_{0}^{j},\hat{Z}_{0}^{j},\hat{Q}_{1}^{j},\hat{Z}_{1}^{j},\dots,\hat{Q}_{0}^{j+1},\hat{Z}_{0}^{j+1},\hat{Q}_{1}^{j+1},\hat{Z}_{1}^{j+1},\dots$

An important realization is that if we already have all the reflectors we can
apply all the $\hat{Z}$ first and apply the $\hat{Q}$ later because of
associativity. This leads to the following order

$\hat{Z}_{0}^{j},\hat{Z}_{1}^{j},\dots,\hat{Z}_{0}^{j+1},\hat{Z}_{1}^{j+1},\dots,\hat{Q}_{0}^{j},\hat{Q}_{1}^{j},\dots,\hat{Q}_{0}^{j+1},\hat{Q}_{1}^{j+1},\dots$

One might hope that because the updates from the left and right are now
separated, applying them will be more efficient. Unfortunately, there is not
much cache reuse. The entries affected by $\hat{Q}_{k}^{j}$ do not overlap
with the entries affected by $\hat{Q}_{k+1}^{j}$, so there is almost no cache
reuse if the reflectors are applied in this order. We have some additional
freedom in choosing the order. Reflectors that do not overlap commute. This
ultimately means that as long as a reflector $\hat{Q}_{k}^{j}$ is preceded by
$\hat{Q}_{k}^{j-1}$ and $\hat{Q}_{k+1}^{j-1}$ (and similarly for $\hat{Z}$) we
can choose any order we want. A reflector $\hat{Q}_{k}^{j}$ affects $r$
columns (or rows in the case of $\hat{Z}_{k}^{j}$), $r-1$ of which are shared
with $\hat{Q}_{k}^{j+1}$, so there is much more potential for cache reuse.
Grouping the reflectors by $k$ instead of $j$ is likely more efficient.
Consider the following sequence

$\hat{Z}_{l-1}^{j},\hat{Z}_{l-1}^{j+1},\dots,\hat{Z}_{l-2}^{j},\hat{Z}_{l-2}^{j+1},\dots,\hat{Q}_{l-1}^{j},\hat{Q}_{l-1}^{j+1},\dots,\hat{Q}_{l-2}^{j},\hat{Q}_{l-2}^{j+1},\dots,$

where $l$ is the maximal value of $k$. It is easy to see that the reflectors
are now grouped by $k$ and that the previously mentioned constraint is not
violated. This technique was initially discovered by C. Bishof, X. Sun, and B.
Lang [3] for tridiagonal reductions. Processing the reflectors in this order
also allows us to use WY representations. Each group of $q$ reflectors
belonging to the same $k$ can be accumulated into a block reflector so that we
can use matrix-matrix multiplications.

#### Generating the reflectors

Now that we know how to efficiently apply the reflectors, we need a way to
calculate the reflectors for several sweeps while updating as few entries as
possible. This amounts to the following task. Generate the reflectors for the
columns $j=j_{1}:j_{1}+q-1$, while only applying the reflectors
$\hat{Q}_{k}^{j}$ to a range of columns $c_{1A}(k,j):c_{2A}(k,j)$ of $A$ and
columns $c_{1B}(k,j):c_{2B}(k,j)$ of $B$ and applying the reflectors
$\hat{Z}_{k}^{j}$ to a range of rows $r_{1A}(k,j):r_{2A}(k,j)$ of $A$ and rows
$r_{1B}(k,j):r_{2B}(k,j)$ of $B$.

Let us start by deriving expressions for $A$. In order to generate the
reflector $\hat{Q}_{k}^{j+1}$,

$r_{1A}(k,j)\leq j_{1}+kr+1-r$ (1)

and in order to generate the reflector $\hat{Q}_{k+1}^{j}$,

$r_{2A}(k,j)\geq\min(j+(k+2)r,n).$ (2)

We must also respect the constraints that were previously mentioned,
$Z_{k}^{j}$ must be preceded by $Z_{k}^{j-1}$ and $Z_{k+1}^{j-1}$, this means
that

$r_{1A}(k,j)\geq r_{1A}(k,j-1)$ (3)

and

$r_{1A}(k,j)\geq r_{1A}(k-1,j-1)$

Setting $r_{1A}(k,j_{1}+q-1)=j_{1}+kr+1-r$ and combining that with equation
(3.2) we get that $r_{1A}(k,j_{1}+q-2)=j_{1}+kr+1-2r$. We can continue this
logic to find that $r_{1A}(k,j_{1}+q-3)=j_{1}+kr+1-3r$ and we ultimately get

$\begin{split}r_{1A}(k,j)&=j_{1}+1+\max(0,kr-r-(j_{1}+q-1-j)r)\\\
r_{2A}(k,j)&=\min(j+(k+2)r,n).\end{split}$ (4)

For $Q_{k}^{j}$, the range is simpler, we can simply set

$c_{1A}(k,j)=j+\max(0,(k-1)r+1)$

and

$c_{2A}(k,j)=j_{1}+q-1+\max(0,(k-1)r+1)$

These are the minimal ranges required for generating the reflectors and it
does not need to be expanded to account for overlapping updates. These two
ranges were first presented by L. Karlsson et al. [8] for Hessenberg matrices.

A similar derivation for $B$ results in:

$\begin{split}r_{1B}(k,j)&=j_{1}+1+\max(0,kr-r-(j_{1}+q-1-j)r)\\\
r_{2B}(k,j)&=\min(j+(k+1)r,n)\end{split}$ (5)

and

$\begin{split}c_{1B}(k,j)&=j+kr+1\\\
c_{2B}(k,j)&=\min(j_{1}+q-1+(k+1)r,n).\end{split}$ (6)

Algorithm 3 generates the reflectors for a sequence of columns while only
updating a small part of $A$ and $B$. We invite the reader to study the
differences between Algorithm 2 and Algorithm 3. Algorithm 4 applies the rest
of the updates using reordered reflectors and block reflectors. Figures 5 and
6 show the parts of the pencil that are updated from the left and right while
generating the reflectors.

(a) $A$ after Sweep 1

(b) $A$ after Sweep 2

(c) $A$ after Sweep 3

(d) $B$ after Sweep 1

(e) $B$ after Sweep 2

(f) $B$ after Sweep 3

Figure 5: Parts of the matrices $A$ and $B$ that are updated from the left
during the generate phase. In light gray, gray and dark gray respectively:
parts updated by the reflectors from sweeps 1, 2 and 3. Note that during each
sweep, we apply the reflectors of the previous sweeps to one extra column.
This cannot be done earlier because it would interfere with the updates from
the right.

(a) $A$ after sweep 1

(b) $A$ after sweep 2

(c) $A$ after sweep 3

(d) $B$ after sweep 1

(e) $B$ after sweep 2

(f) $B$ after sweep 3

Figure 6: Parts of the matrices $A$ and $B$ that are updated from the right
during the generate phase. In light gray, gray and dark gray respectively: the
parts updated by the reflectors from sweeps 1, 2 and 3. Note that the parts of
the matrices that get updated are much larger than for the left updates.

###  3.3 Parallel stage two

The parallelization strategy of stage two follows the same general principle
as stage one. We split the generation of the reflectors and the updates of the
different matrices into different tasks. What is special here is that we also
split the application tasks into lookahead tasks and smaller application
tasks. During these lookahead tasks, we update $A$ and $B$ just enough so that
the $O(rq)$ band required for the generation of the reflectors is fully
updated. This allows us to overlap the updating of the matrices with the
generation of the reflectors. This is important because the generation of the
reflectors is neither negligible nor easily parallelizable. Figure 7
illustrates the dependencies between the tasks of different iterations.

Next, we parallelize the large tasks in the same way as in stage one. For both
the lookahead and update tasks, we split the matrices into either column or
row slices and apply the updates in parallel. Figure 8 shows this
distribution.

$\text{Generate}^{(1)}$Lookahead $A^{(1)}$Lookahead $B^{(1)}$Update
$Q^{(1)}$Update $Z^{(1)}$Update $A^{(1)}$Update
$B^{(1)}$$\text{Generate}^{(2)}$Lookahead $A^{(2)}$Lookahead $B^{(2)}$Update
$Q^{(2)}$Update $Z^{(2)}$Update $A^{(2)}$Update $B^{(2)}$ Figure 7: Dependency
graph of parallel stage two. The superscript indicates which iteration the
task belongs to.

(a) full updates from the right and generation of reflectors

(b) full updates from the left and generation of reflectors

(c) lookahead updates from the right

(d) lookahead updates from the left

Figure 8: Distribution of the update, lookahead and generate tasks of phase 2
into multiple subtasks. In dark gray, the generate task, this task cannot be
parallelized effectively and is not split into subtasks. The full update tasks
are shown in light gray. They encompass most of the flops of the algorithm and
are split into the most subtasks. The lookahead updates are shown in gray.
These tasks encompass fewer flops so they are split into larger blocks.

##  4 Experiments

We compare the following implementations

  * •

ParaHT: The algorithm presented in this paper

  * •

IterHT: Iterative reduction algorithm by T. Steel and R. Vandebril [11]

  * •

DGGHD3: The blocked Hessenberg-triangular reduction by Kågström et al. [7] as
implemented in LAPACK 3.9 [2].

  * •

HouseHT: Householder reflector based reduction by Bujanovic, L. Karlsson and
D. Kressner [5].

The tests are performed on a machine with two Intel Xeon E-2697 v3 CPUs (14
cores each) that share 128GB of memory. A single node of this CPU has a
theoretical peak double precision flop rate of 20.8 Gflops.

ParaHT, IterHT, and DGGHD3 are written in Fortran and compiled with GNU
Fortran version 7.5.0 with optimization flags `-O3 -march=native`. HouseHT is
written in C++ and compiled with GNU C Compiler version 7.5.0 with
optimization flags `-O3 -march=native`. The code is linked with MKL version
2019.0.1.

For HouseHT, the parameter $n_{b}$ is set to 64. For ParaHT, the parameters
$r$, $p$, and $q$ are set to 16, 8, and 8 respectively. These parameters were
tuned by running the algorithm on randomly generated pencils of different
sizes and choosing the set of parameters that performed the best on average.

The experiments will focus on the runtime of the algorithms. We will not go
into detail about the accuracy of the methods. Each of the algorithms that we
tested always produces results that have relative backward errors on the order
of the machine precision.

#### Tests on random pencils

To evaluate the general performance of our algorithm, we run it for randomly
generated matrices of varying sizes and compare it with other implementations.
After generating the random pencil, we take a QR factorization of $B$ so that
it is upper triangular. Note that a random matrix is usually well conditioned,
which is important for two of the algorithms we compare to. IterHT and to a
lesser degree HouseHT are iterative algorithms that require more iterations if
the pencil has many infinite eigenvalues. Pencils of which the matrix $B$ has
a small condition number avoid these extra iterations.

Figure 9(a) shows the runtime of several algorithms relative to the sequential
runtime of LAPACK for a varying number of threads. Because of the extra flops,
our algorithm is significantly slower than the other algorithms on one core.
As the number of threads increases, the effective parallelization makes up for
the increased cost. We also note that HouseHT and IterHT achieve their highest
parallel speedup when using 14 threads. In other tests that do not explicitly
vary the number of threads, we limit HouseHT and IterHT to 14 threads to get a
fair comparison.

Figure 9(b) shows the speedup our algorithm achieves over LAPACK, HouseHT and
IterHT. We achieve a speedup of 2 over HouseHT. Against LAPACK, our algorithm
is slightly slower for small matrices but achieves a speedup of 4 for large
matrices. IterHT is faster than our algorithm in most cases. As mentioned
before, the runtime of IterHT strongly depends on the number of infinite
eigenvalues present in the pencil and the condition of $B$. In most of the
tests, because $B$ is well-conditioned, it only needs a single iteration. For
$n=8000$ and $n=16000$, it needs two iterations and in those cases, our
algorithm is slightly faster than IterHT.

Figure 10 shows the parallel performance of ParaHT in more detail. First, note
that most of the runtime of the algorithm is spent inside phase 2 despite
phase 1 requiring slightly more flops. This indicates that phase 2 is less
efficient. Second, the parallel speedup of phase 1 and phase 2 (and
consequently, the full algorithm) is very similar. The issues that make phase
2 less efficient on a single core persist when executing the phases in
parallel. Finally, we can also see that the matrices need to be quite large to
achieve good speedups. For $n=1000$, our algorithm only achieves a parallel
speedup of around 2, but for $n=8000$, our algorithm achieves a parallel
speedup of 10.

124814280.512468Number of threadsParallel speedupLAPACKHouseHTIterHTParaHT (a)
Parallel speedup

5001000200040008000160000.51248$n$SpeedupLAPACKHouseHTIterHT (b) Speedup for
varying n

Figure 9: On the left, the parallel speedup (fraction of single-threaded
runtime of LAPACK and multi-threaded runtime of the relevant algorithm) for a
randomly generated pencil of size 8000. On the right, the speedup ParaHT
achieves over other algorithms for varying pencil sizes.

12481428$0.3$$0.4$$0.5$$0.6$$0.7$Number of threadsRelative runtimephase 1phase
2 (a) Runtime of phase 1 and 2 relative to the total runtime

12481428124810Number of threadsParallel speedup1000200040008000 (b) Parallel
speedup of ParaHT

12481428124810Number of threadsParallel speedup1000200040008000 (c) Parallel
speedup of phase 1

12481428124810Number of threadsParallel speedup1000200040008000 (d) Parallel
speedup of phase 2

Figure 10: Parallel speedup and relative runtime of our algorithm and its
phases.

#### Tests on saddle point problems

Following Bujanovic et al. [5], we also test our algorithm on saddle point
problems. These are pencils of the form:

$(A,B)=\left(\begin{bmatrix}X&Y\\\ Y^{*}&0\end{bmatrix},\begin{bmatrix}I&0\\\
0&0\end{bmatrix}\right),$

where $I$ is the identity matrix, $Y$ is a random matrix and $X$ is a random
symmetric positive definite matrix. By choosing the dimension of $X$ and the
identity matrix, we can control the number of infinite eigenvalues in the
pencil. In our experiments, the dimensions have been chosen so that $25\%$ of
the eigenvalues of the pencil are infinite. Such a pencil is particularly
difficult for HouseHT and IterHT.

We note that while pencils with a large number of infinite eigenvalues exist
in practice, it is often possible to deflate these eigenvalues prior to the
Hessenberg-triangular reduction. LAPACK can detect infinite eigenvalues based
on the sparsity structure of the pencil and a different deflation technique we
presented in a previous paper [11] can deflate eigenvalues that are not
exactly infinite; but can be considered numerically infinite in finite
precision.

Figure 11 shows the speedup of ParaHT for saddle point problems of varying
sizes. As expected, the speedup over LAPACK has not changed significantly,
because the runtimes of LAPACK and ParaHT do not depend on the number of
infinite eigenvalues. The speedup over HouseHT is much larger than before;
because HouseHT must do more iterative refinement to deal with the ill-
conditioned matrix. IterHT is even more sensitive to the presence of a large
number of infinite eigenvalues and fails to converge.

5001000200040008000160001248$n$SpeedupLAPACKHouseHT Figure 11: The speedup
ParaHT achieves over other algorithms for saddle point pencils of varying
size. IterHT is not listed because it failed to converge within 10 iterations
of iterative refinement.

##  5 Conclusion

We have presented a parallel algorithm for
Hessenberg-triangular reduction. The algorithm consists of two stages: a
reduction to r-Hessenberg-triangular form and a reduction to Hessenberg-
triangular form. Two-stage reductions typically require more operations than
one-stage reductions but are easier to parallelize and our new formulation is
no different. Experiments in a shared memory environment have shown that our
algorithm is accurate and can outperform state-of-the-art sequential
algorithms using parallel BLAS. Future work will need to determine whether our
algorithm can also outperform parallelized versions of these algorithms. As an
extra advantage, our algorithm does not suffer a performance penalty when
infinite eigenvalues are present in the pencil.

## Acknowledgements

We thank Bujanovic, Karlsson, and Kressner for providing us with their
implementation of HouseHT to compare against.

## References

  * [1] Björn Adlerborn, Lars Karlsson, and Bo Kgström.  Distributed one-stage Hessenberg-triangular reduction with wavefront scheduling.  SIAM Journal on Scientific Computing, 40(2):C157–C180, 2018.
  * [2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen.  LAPACK Users’ Guide.  Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.
  * [3] Christian Bischof, Xiaobai Sun, and Bruno Lang.  Parallel tridiagonalization through two-step band reduction.  In Proceedings of IEEE Scalable High Performance Computing Conference, pages 23–27. IEEE, 1994.
  * [4] Christian Bischof and Charles Van Loan.  The WY representation for products of Householder matrices.  SIAM Journal on Scientific and Statistical Computing, 8(1):s2–s13, 1987.
  * [5] Zvonimir Bujanovic, Lars Karlsson, and Daniel Kressner.  A Householder-based algorithm for Hessenberg-triangular reduction.  SIAM Journal on Matrix Analysis and Applications, 39(3):1270–1294, 2018.
  * [6] Krister Dackland and Bo Kgström.  A Scalapack-style algorithm for reducing a regular matrix pair to block Hessenberg-triangular form.  In International Workshop on Applied Parallel Computing, pages 95–103. Springer, 1998.
  * [7] Bo Kgström, Daniel Kressner, ES Quintana-Ortí, and Gregorio Quintana-Orti.  Blocked algorithms for the reduction to Hessenberg-triangular form revisited.  BIT Numerical Mathematics, 48(3):563–584, 2008.
  * [8] Lars Karlsson and Bo Kgström.  Efficient reduction from block Hessenberg form to Hessenberg form using shared memory.  In International Workshop on Applied Parallel Computing, pages 258–268. Springer, 2010.
  * [9] Lars Karlsson and Bo Kgström.  Parallel two-stage reduction to Hessenberg form using dynamic scheduling on shared-memory architectures.  Parallel Computing, 37(12):771–782, 2011.
  * [10] C. B. Moler and G. W. Stewart.  An algorithm for generalized matrix eigenvalue problems.  SIAM Journal on Numerical Analysis, 10(2):241–256, 1973.
  * [11] Thijs Steel and Raf Vandebril.  A novel, blocked algorithm for the reduction to Hessenberg-triangular form.  The Electronic Journal of Linear Algebra, Accepted for publication, 2023.
  * [12] David S Watkins.  Performance of the QZ algorithm in the presence of infinite eigenvalues.  SIAM Journal on Matrix Analysis and Applications, 22(2):364–375, 2000.

## Appendix

0: $n\times n$ pencil $(A,B)$, with $B$ upper triangular

10: decomposition $(A_{orig},B_{orig})=Q(A,B)Z^{*}$, with $(A,B)$ in
r-Hessenberg-triangular form

1: $Q=I$

2: $Z=I$

3: for $j=1:n_{b}:n-2$ do

4: $j_{1}=j$

5: $j_{2}=\min(n,j+n_{b}-1)$

6: $n_{blocks}=\lceil(n-nb-j+1)/((p-1)n_{b})\rceil$

7: for $k=n_{blocks}-1:-1:0$ do

8: $i_{1}=j+n_{b}+k(p-1)n_{b}$

9: $i_{2}=\min(n,i_{1}+p*n_{b}-1)$

10: $A(i_{1}:i_{2},j_{1}:j_{2})=\hat{Q}_{k}^{j}\hat{R}$ {QR factorization}

11: $A(i_{1}:i_{2},j_{1}:j_{2})=\hat{R}$

12: $A(i_{1}:i_{2},j_{2}+1:n)=\hat{Q}_{k}^{j*}A(i_{1}:i_{2},j_{2}+1:n)$

13: $B(i_{1}:i_{2},i_{1}:n)=\hat{Q}_{k}^{j*}B(i_{1}:i_{2},i_{1}:n)$

14: $Q(1:n,i_{1}:i_{2})=Q(1:n,i_{1}:i_{2})\hat{Q}_{k}^{j}$

215: end for

16: for $i=j+n_{b}+(n_{blocks}-1)(p-1)n_{b}:-(p-1)n_{b}:j+n_{b}$ do

17: $i_{1}=i$

318: $i_{2}=\min(n,i+p*n_{b}-1)$

19: $B(i_{1}:i_{2},i_{1}:i_{2})=\tilde{R}\tilde{Q}${RQ factorization}

20: $\tilde{Q}(1:\min(n_{b},i_{2}-i_{1}+1),1:i_{2}-i_{1}+1)=L\hat{Z}${LQ
factorization}

21: $A(1:n,i_{1}:i_{2})=A(1:n,i_{1}:i_{2})\hat{Z}$

22: $B(1:i_{2},i_{1}:i_{2})=B(1:i_{2},i_{1}:i_{2})\hat{Z}$

23: $Z(1:n,i_{1}:i_{2})=Z(1:n,i_{1}:i_{2})\hat{Z}$

24: end for

25: end for

Algorithm 1 Blocked reduction to r-Hessenberg-triangular form [6]

0: $n\times n$ pencil $(A,B)$ in r-Hessenberg-triangular form

10: decomposition $(A_{orig},B_{orig})=Q(A,B)Z^{*}$, with $(A,B)$ in
Hessenberg-triangular form

1: $Q=I$

2: $Z=I$

3: for $j=1:n-2$ do

4: $n_{blocks}=1+\lfloor\frac{n-j-2}{r}\rfloor$

5: for $k=0:n_{blocks}-1$ do

6: $j_{b}=j+\max(0,(k-1)r+1)$

7: $i_{1}=j+kr+1$

8: $i_{2}=\min(j+(k+1)r,n)$

9: $i_{3}=\min(j+(k+2)r,n)$

10: Generate a reflector $\hat{Q}_{k}^{j}=I-\tau_{l}v_{l}v_{l}^{*}$ that
reduces $A(i_{1}:i_{2},j_{b})$

11: $A(i_{1}:i_{2},j_{b}:n)=\hat{Q}_{k}^{j}A(i_{1}:i_{2},j_{b}:n)$

12: $B(i_{1}:i_{2},i_{1}:n)=\hat{Q}_{k}^{j}B(i_{1}:i_{2},i_{1}:n)$

13: $Q(:,i_{1}:i_{2})=Q(:,i_{1}:i_{2})\hat{Q}_{k}^{j}$

14: $B(i_{1}:i_{2},i_{1}:i_{2})=\tilde{R}\tilde{Q}${RQ factorization}

15: Generate a reflector $\hat{Z}_{k}^{j}=I-\tau_{r}v_{r}v_{r}^{*}$ that
reduces $\tilde{Q}(1,1:i_{2}-i_{1}+1)$

16: $A(1:i_{3},i_{1}:i_{2})=A(1:i_{3},i_{1}:i_{2})\hat{Z}_{k}^{j}$

17: $B(1:i_{2},i_{1}:i_{2})=B(1:i_{2},i_{1}:i_{2})\hat{Z}_{k}^{j}$

18: $Z(:,i_{1}:i_{2})=Z(:,i_{1}:i_{2})\hat{Z}_{k}^{j}$

19: end for

20: end for

Algorithm 2 Reduction to Hessenberg-triangular form without blocking

0: $n\times n$ pencil $(A,B)$ in r-Hessenberg-triangular form, where $A$ is in
Hessenberg form in columns $1:j_{1}-1$.

1: for $j=j_{1}:j_{1}+q-1$ do

2: $n_{blocks}=2+\lfloor\frac{n-j-1}{r}\rfloor$

3: for $k=0:n_{blocks}-1$ do

4: $j_{b}=j+\max(0,(k-1)r+1)$

5: $i_{1}=j+kr+1$

6: $i_{2}=\min(j+(k+1)r,n)$

7: $i_{3}=\min(j+(k+2)r,n)$

18: $i_{4}=j_{1}+1+\max(0,(k+j-j_{1}-q+2)r)$

9: for $\hat{j}=j_{1}:j-1$ do

10: $\hat{i_{1}}=\hat{j}+kr+1$

11: $\hat{i_{2}}=\min(\hat{j}+(k+1)r,n)$

12: if $\hat{i_{2}}-\hat{i_{1}}\geq 1$ then

13:
$A(\hat{i_{1}}:\hat{i_{2}},j_{b})=\hat{Q}_{k}^{\hat{j}}A(\hat{i_{1}}:\hat{i_{2}},j_{b})$

14: if $i_{1}+r-1\leq n$ then

15:
$B(\hat{i_{1}}:\hat{i_{2}},i_{1}+r-1)=\hat{Q}_{k}^{\hat{j}}B(\hat{i_{1}}:\hat{i_{2}},i_{1}+r-1)$

16: end if

17: end if

18: end for

19: Generate a reflector $\hat{Q}_{k}^{j}=I-\tau_{l}v_{l}v_{l}^{*}$ that
reduces $A(i_{1}:i_{2},j_{b})$

20: $A(i_{1}:i_{2},j_{b})=\hat{Q}_{k}^{j}A(i_{1}:i_{2},j_{b})$

21: $B(i_{1}:i_{2},i_{1}:i_{2})=\hat{Q}_{k}^{j}B(i_{1}:i_{2},i_{1}:i_{2})$

22: $B(i_{1}:i_{2},i_{1}:i_{2})=\tilde{R}\tilde{Q}${RQ factorization}

23: Generate a reflector $\hat{Z}_{k}^{j}=I-\tau_{r}v_{r}v_{r}^{*}$ that
reduces $\tilde{Q}(1,1:i_{2}-i_{1}+1)$

24: $A(i_{4}:i_{3},i_{1}:i_{2})=A(i_{4}:i_{3},i_{1}:i_{2})\hat{Z}_{k}^{j}$

25: $B(i_{4}:i_{2},i_{1}:i_{2})=B(i_{4}:i_{2},i_{1}:i_{2})\hat{Z}_{k}^{j}$

26: end for

27: end for

Algorithm 3 Generate phase of stage two

0: $n\times n$ pencil $(A,B)$ in r-Hessenberg-triangular form, where $A$ is in
Hessenberg form in columns $1:j_{1}-1$.

1: $n_{blocks}=1+\lfloor\frac{n-j_{1}-2}{r}\rfloor$

2: for $k=n_{blocks}-1:-1:0$ do

3: $i_{5}=j_{1}+1+\max(0,(k-q+2)r)$

4: for $j=j_{1}+1:j_{1}+q-1$ do

5: $i_{1}=j+kr+1$

6: $i_{2}=\min(j+(k+1)r,n)$

7: $i_{4}=j_{1}+\max(0,(k+j-j_{1}-q+2)r)$

8: $A(i_{5}:i_{4},i_{1}:i_{2})=A(i_{5}:i_{4},i_{1}:i_{2})\hat{Z}_{k}^{j}$

9: $B(i_{5}:i_{4},i_{1}:i_{2})=B(i_{5}:i_{4},i_{1}:i_{2})\hat{Z}_{k}^{j}$

10: end for

11:
$\hat{Z}_{k}=\hat{Z}_{k}^{j_{1}}\hat{Z}_{k}^{j_{2}}\dots\hat{Z}_{k}^{j_{1}+q-1}${Form
compact WY}

12: $i_{1}=j_{1}+kr+1$

13: $i_{2}=\min(j_{1}+q-1+(k+1)r,n)$

14: $i_{5}=j_{1}+\max(0,(k-q+2)r)$

15: $A(1:i_{5},i_{1}:i_{2})=A(1:i_{5},i_{1}:i_{2})\hat{Z}_{k}$

16: $B(1:i_{5},i_{1}:i_{2})=B(1:i_{5},i_{1}:i_{2})\hat{Z}_{k}$

17: $Z(:,i_{1}:i_{2})=Z(:,i_{1}:i_{2})\hat{Z}_{k}$

18: end for

19: for $k=n_{blocks}-1:-1:0$ do

20:
$\hat{Q}_{k}=\hat{Q}_{k}^{j_{1}}\hat{Q}_{k}^{j_{2}}\dots\hat{Q}_{k}^{j_{1}+q-1}${Form
compact WY}

21: $i_{1}=j_{1}+kr+1$

22: $i_{2}=\min(j_{1}+q-1+(k+1)r,n)$

23: $i_{5}=j_{1}+q-1+\max(0,(k-1)r+1)$

24: $i_{6}=j_{1}+q+(k+1)r-1$

25: $A(i_{1}:i_{2},i_{5}:n)=\hat{Q}_{k}^{*}A(i_{1}:i_{2},i_{5}:n)$

26: $B(i_{1}:i_{2},i_{6}:n)=\hat{Q}_{k}^{*}B(i_{1}:i_{2},i_{6}:n)$

27: $Q(:,i_{1}:i_{2})=Q(:,i_{1}:i_{2})\hat{Q}_{k}$

28: end for

Algorithm 4 Application phase of stage two