11institutetext: University College London
Department of Computer Science, United Kingdom 22institutetext: University of
Oxford
Oxford Robotics Institute, United Kingdom
22email<EMAIL_ADDRESS>This work was supported by the
SecondHands project, funded from the EU Horizon 2020 Research and Innovation
programme under grant agreement 643950 and by the EPSRC grants RAIN and ORCA
(EP/R026084/1, EP/R026173/1).

# 3D Pick & Mix: Object Part Blending in Joint Shape and Image Manifolds

Adrian Penate-Sanchez
Lourdes Agapito 1122 0000-0003-2876-3301 11 0000-0002-6947-1092

###### Abstract

We present 3D Pick & Mix, a new 3D shape retrieval system that provides users
with a new level of freedom to explore 3D shape and Internet image collections
by introducing the ability to reason about objects at the level of their
constituent parts. While classic retrieval systems can only formulate simple
searches such as _“find the 3D model that is most similar to the input image”_
our new approach can formulate advanced and semantically meaningful search
queries such as: _“find me the 3D model that best combines the design of the
legs of the chair in image 1 but with no armrests, like the chair in image
2”_. Many applications could benefit from such rich queries, users could
browse through catalogues of furniture and pick and mix parts, combining for
example the legs of a chair from one shop and the armrests from another shop.

###### Keywords:

Shape blending Image embedding Shape retrieval

##  1 Introduction

As databases of images and 3D shapes keep growing in size and number,
organizing and exploring them has become increasingly complex. While most
tools so far have dealt with shape and appearance modalities separately, some
recent methods [15, 10] have begun to exploit the complementary nature of
these two sources of information and to reap the benefits of creating a common
representation for images and 3D models. Once images and 3D shapes are linked
together, many possibilities open up to transfer what is learnt from one
modality to another. Creating a joint embedding allows to retrieve 3D models
based on image queries (or vice-versa) or to align images of similar 3D
shapes. However, recent retrieval methods still fall short of being flexible
enough to allow advanced queries. Crucially, they are limited to reasoning
about objects as a whole – taking a single query image (or shape) as input at
test time prevents them from combining object properties from different
inputs. 3D Pick & Mix overcomes this limitation by reasoning about objects at
the level of parts. It can formulate advanced queries such as: _“find me the
3D model that best combines the design of the backrest of the chair in image
$1$ with the shape of the legs of the chair in image $2$”_ (see Fig. 1) or
_“retrieve chairs with wheels”_. The ability to reason at the level of parts
provides users with a new level of freedom to explore 3D shape and image
datasets. Users could browse through catalogues of furniture and pick and mix
parts, combining for example the legs of a favourite chair from one catalogue
and the armrests from another (see Fig. 1).

Figure 1:  Our approach takes two (or more) inputs, either RGB images or 3D
models, and part labels such as _”legs”_ or _”backrest”_ and retrieves a shape
that combines object parts or properties from all inputs, via a cross-manifold
optimization technique.

Our system first builds independent manifold shape spaces for each object part
(see Fig. 4). A CNN-based manifold coordinate regressor is trained to map real
images of an object to the part manifolds. Our novel deep architecture jointly
performs semantic segmentation of object parts and coordinate regression for
each part on the corresponding shape manifold. This network is trained using
only synthetic data, Fig. 2 illustrates the architecture. At test time the
user provides two (or more) images (or 3D models) as input and determines
which parts they would like to pick from each (note that this only requires a
label name such as _’legs’_). Our system retrieves the model that best fits
the arrangement of parts by performing a cross-manifold optimization (see Fig.
1). The main contributions of our 3D Pick & Mix system are:

  * •

We learn embeddings (manifolds) for object parts (for instance the legs or the
armrests of chairs) which allow us to retrieve images or 3D models of objects
with similarly shaped parts.

  * •

We propose a new deep architecture that can map RGB images onto these
manifolds of parts by regressing their coordinates. Crucially, the input to
the network is simply an RGB image and the name (label) of the object part.
The CNN learns to: _(i)_ segment the pixels that correspond to the chosen
part, and _(ii)_ regress its coordinates on the shape manifold.

  * •

At query time our retrieval system can combine object parts from multiple
input images, enabled by a cross-manifold optimization technique.

##  2 Related Work

Joint 3D model/image embeddings: While most shape retrieval methods had
traditionally dealt with shape and appearance modalities separately, a recent
trend has emerged that exploits the complementary nature of appearance and
shape by creating a common representation for images and 3D models. [10]
exploits the different advantages of shape and images by using the robustness
of 3D models for alignment and pose estimation and the reliability of image
labels to identify the objects. While they do not explicitly create a joint
manifold based on shape similarity they do rely on image representations for
both modalities. Another example of 3D model/image embedding is [15] who first
builds a manifold representation of 3D shapes and then trains a CNN to
recognize the shape of the object in an image. Unlike our approach, both [10,
15], limit their representations to objects as a whole preventing the
combination of properties taken from different inputs. [23] perform shape
retrieval from sketches, words, depth maps and real images by creating a
manifold space that combines the different inputs. Since intra-class
similarity is not the main focus, most instances of the same class tend to
appear clustered. [16] learn a manifold-space metric by using triplets of
shapes where the first is similar to the third but dissimilar to the second.
Similarly to our approach, the metric space is defined based on shape and not
image similarity. [9] first generates voxel representations of the objects
present in the RGB image inputs. A shared latent shape representation is then
learnt for both images and the voxelized data. At test time RGB convolutions
and volume generation deconvolution layers are used to produce the 3D shape.

3D shape blending/mixing: Much in the line of the work presented in this
paper, there has been fruitful research in shape blending in recent years. The
_“3D model evolution“_ approach of [27] takes a small set of 3D models as
input to generate many. Parts from two models cross-over to form a new 3D
model, continuing to merge original models with new ones to generate a large
number of 3D models. In [1] new shapes are generated by interpolating and
varying the topology between two 3D models. The photo-inspired 3D modeling
method of [28] takes a single image as input, segments it into parts using an
interactive model-driven approach, then retrieves a 3D model candidate that is
finally deformed to match the silhouette of the input photo. The probabilistic
approach of [11] learns a model that describes the relationship between the
parts of 3D shapes which then allows to create an immense variety of new
blended shapes by mixing attributes from different models. The sketch driven
method of [26] edits a pre- segmented 3D shape using user-drawn sketches of
new parts. The sketch is used to retrieve a matching 3D part from a catalogue
of 3D shapes which is then snapped onto the original 3D shape to create a new
blended 3D shape. Note that the above approaches use only 3D shapes as input
for shape blending, with the exception of [28] who use a single photograph and
[26] who use sketches. However, unlike ours, neither of these approaches can
combine different input images to retrieve a shape that blends parts from each
input.

Modeling of 3D object parts: We will differentiate between 3D segmentation
approaches that seek to ensure consistency in the resulting segmentation
across different examples of the same object class (co-segmentation) and those
that seek a semantically meaningful segmentation (semantic segmentation). Some
recent examples of approaches that perform co-segmentation can be found in
[24, 8], but as we seek to describe parts that have meaning to humans we will
focus on the later. We can find examples of semantic 3D parts in approaches
like [29]. [29] provides accurate semantic region annotations for large
geometric datasets with a fraction of the effort by alternating between using
few manual annotations from an expert and a system that propagates labels to
new models. We exploit the ShapeNet annotations provided by [29] as the ground
truth part shape when constructing our joint manifold.

Recognition of 3D structure from images: The exemplar-based approach of [2]
performs joint object category detection viewpoint estimation, exploiting 3D
model datasets to render instances from different viewpoints and then learn
the combination of viewpoint-instance using exemplar SVMs. [5] uses 3D
Convolutional LSTMs to extract the 3D shape of an object from one or more
viewpoints. By using LSTM blocks that contain memory, they progressively
refine the shape of the object. [7] learn to generate a 3D point cloud from a
single RGB image, it learns purely from synthetic data. By using a point cloud
instead of a volumetric representation better definition of the details of the
shape are obtained. Their novel approach learns how to generate several
plausible 3D reconstructions from a single RGB image at test time if the
partial observation of the image is ambiguous. [22] learn to recognize the
object category and the camera viewpoint for an image using synthetically
generated images for training. This work showed that datasets of real images
annotated with 3D information were not required to learn shape properties from
images as this could be learnt from synthetically generated renderings. [21]
obtain good depth estimates for an image given a set of 3D models of the same
class.

Figure 2: Summary of the architecture of ManifoldNet, our new deep network
that takes an image as input and learns to regress the coordinates of each
object part in the different part manifolds. The architecture has 3 sections:
the first set of layers performs semantic segmentation of the image pixels
into different semantic parts (such as _”backrest”, “seat”, “armrests” or
“legs”_ in the case of chairs). The second section learns an intermediate
feature representation for manifold coordinate regression. The final section
learns to regress the shape coordinates in each of the part manifolds. We show
the nearest neighbour shapes found on the _“armrests”_ and _“legs”_ manifolds
for the depicted input image.

##  3 Overview

In this section we provide a high level overview of our 3D Pick & Mix
retrieval system. Our system requires a training stage in which: _(i)_
manifolds of 3D shapes of object parts are built (see Fig. 4) and _(ii)_ a CNN
is trained to take as input an image and regress the coordinates of each of
its constituent parts on the shape manifolds (illustrated in Fig. 2). At query
time the system receives an image or set of images as input and obtains the
corresponding coordinates on the part manifolds. If the user chooses object
parts from different images a cross-manifold optimization is carried out to
retrieve a single shape that blends together properties from different images.

Training: At training time, our method takes as input a class-specific
collection of 3D shapes (we used $ShapeNet$ [3]) for which part label
annotations are available. The first step at training time is to _learn a
separate shape manifold for each object part_ (see Fig. 4). Each shape is
represented with a Light Field descriptor [4] and characterized with a pyramid
of HoG features. The manifolds are then built using non-linear multi-
dimensional-scaling (MDS) and the $L_{2}$ norm between feature vectors as the
distance metric – in each resulting low-dimensional manifold, objects that
have similarly shaped parts are close to each other. So far these manifolds of
object parts (for instance _back-rests_ , _arm-rests_ , _legs_ , _seats_ in
the case of chairs) contain 3D shapes. The second step at training time is to
train a CNN to embed images onto each part manifold by _regressing their
coordinates_. We create a set of synthetic training images with per pixel
semantic segmentation annotations for the object parts and ground truth
manifold coordinates. The architecture of this novel CNN (which we denote
ManifoldNet and is shown in Fig. 2) has three clear parts: a set of fully
convolutional layers for semantic segmentation of the object into parts; a set
of convolutional feature extraction layers; and a set of fully connected
layers for manifold coordinate regression. This architecture can be trained
end-to-end. We give an example of the produced semantic segmentation in Fig.
3.

Retrieval: At test time, given a new query image of an unseen object,
ManifoldNet can embed it into each of the part manifolds by regressing the
coordinates. More importantly, our retrieval system can take more than one
image as input, picking different object parts from each image. Note that
ManifoldNet only needs the input images and the name of the object part that
will be used from each image. The network learns jointly to segment the image
into parts and to regress the manifold coordinates and therefore it does not
require any manual annotations as input. A cross-manifold optimization will
then take the coordinates on each of the part manifolds as input and return
the coordinates of a unique 3D shape that blends the different object parts
together. This is achieved through an energy optimization approach, described
in section 4.4.

Figure 3:  Example of the semantic segmentation performed by the first stages
of our architecture. We can see the output probabilities for each of the parts
and the background give a very strong prior of were the parts of an object can
be found. Not requiring labels for each part in the input image makes our
approach very easy to use and increases dramatically its applicability.

##  4 Methodology

###  4.1 Building shape manifolds for object parts

We choose to create an embedding space that captures the similarity between
the shape of object parts based exclusively on the 3D shapes. The reason
behind this choice is that 3D models capture a more complete, pure and
reliable representation of geometry as opposed to images that often display
occlusions, or other distracting factors such as texture or shading effects.
We then rely on our new CNN architecture to map images onto the same embedding
by regressing their coordinates on the corresponding manifolds.

####  4.1.1 Defining a smooth similarity measure between 3D shapes

Figure 4:  Shape manifold construction. The shape of each part for each 3D
model is rendered from different viewpoints and represented with a Light Field
descriptor [4]. The manifolds are then built using non-linear multi-
dimensional-scaling (MDS) and the $L_{2}$ norm between feature vectors as the
distance metric. In each low-dimensional manifold, objects that have similarly
shaped parts appear close to each other.

Shape similarity between object parts is defined as follows. Given a shape
$S_{i}$, we define its Light-field Descriptor (LfD) [4] $L_{i}$ as the
concatenation of the HoG responses [6] $L_{i}=[H_{1};H_{2};...;H_{k}]$. The
value of $k$ is fixed to $k=20$ throughout this work. The light field
descriptor $H_{k}$ for each view $k$ is defined as
$H_{k}=[H_{k}^{mid};H_{k}^{low}]\in\mathbb{R}^{2610}$. The $L_{2}$ distance
between feature vectors is then used as the similarity measure between a pair
of shapes $S_{i}$ and $S_{j}$:
$d_{ij}=\left|\left|L_{i}-L_{j}\right|\right|_{2}$ where
$L_{i}\in\mathbb{R}^{52200}$. We found that using only the mid and low
frequency parts of the HoG pyramid leads to smoother transitions in shape
similarity. For this reason we do not use the original 3 level HoG pyramids
but just the 2 higher levels of the pyramids. This allows for smooth
transitions in shape similarity between parts making the shape blending
possible. Due to most 3D models available in Internet datasets not being
watertight but only polygon soups we are required to use projective shape
features like [4]. We now build separate manifolds for each object part. Each
shape $S_{i}$ is therefore split into its constituent parts
$\forall\,S_{i}:\exists\,\\{S_{i}^{1};S_{i}^{2};...;S_{i}^{P}\\},$ where $P$
is the total number of parts and $S_{i}^{p}$ is the shape of part $p$ of
object $i$. If a part is not present in an object (for instance, chairs
without arm-rests) we set all the components of the vector $L_{i}^{p}$ to
zero, which is equivalent to computing the HoG descriptor of an empty image.

###  4.2 Building shape manifolds of parts

Using the similarity measure between the shape of object parts we use it to
construct a low dimensional representation of the shape space. In principle,
the original $L_{i}\in\mathbb{R}^{52200}$ feature vectors could have been used
to represent each shape, since distances in that space reflect well the
similarity between shapes. We reduce the dimensionality from $52,200$ to $128$
dimensions and we use non-linear Multi-Dimensional Scaling (MDS) [14] to build
the shape manifolds. We compute the distance matrix
$D^{p}\in\mathbb{R}^{n\times n}$ as $D^{p}(i,j)=d^{p}_{ij}$, were $p$ is the
index of the part and $n$ is the total number of shapes. The manifold is built
using MDS by minimizing a Sammon Mapping error [19] defined as

$E^{p}=\dfrac{1}{\sum_{i<j}D^{p}(i,j)}\sum_{i<j}\dfrac{(D^{p}(i,j)-D^{\prime
p}(i,j))^{2}}{D^{p}(i,j)},$ (1)

where $D^{p}$ is the distance matrix between shapes in the original high
dimensional feature space $L^{p}$; and $D^{\prime p}$ is the distance matrix
between shapes in the new low dimensional manifold $L^{\prime p}$. With the
different manifolds $L^{\prime p}$ for each part $p$ computed, a low
dimensional representation of shape similarity exists and all 3D shapes are
already included in it. Adding new 3D shapes to the manifold is done by
solving an optimization that minimizes the difference between the distances
between all previous shapes and the new shape in $D^{p}$ and the distances in
$D^{\prime p}$ with respect to the predicted embedding point. To understand
the shape of the produced manifolds we provide a 2D visualization in Fig. 5.
In this figure we can better understand how the manifolds relate the different
parts of an object.

Figure 5:  Two dimensional visualizations of the four low dimensional part
manifolds _“backrest”, “seat”, “armrests”_ and _“legs”_ for the chairs class.
Probabilistic PCA has been used to provide a 2D visualization of the 128
dimensional manifolds. Both images and 3D models have been represented in the
manifolds. Objects with similarly shaped parts lie close to each other on the
manifold. Several shapes are tagged in all four manifolds to show how vicinity
changes for each part. All shapes and images exist in all manifolds.

###  4.3 Learning to embed images into the shape manifolds

Building the shape manifolds $L^{\prime p}$ for each part $p$ based only on 3D
models, we have successfully abstracted away effects such as textures, colours
or materials. The next step is to train a deep neural network that can map RGB
images onto each manifold by regressing the coordinates on each part manifold
directly from RGB inputs. Crucially, the input to the network must be simply
the RGB image and the name (label) of the object part $p$ selected from that
image – for instance _embed this image of a chair into the manifold of “chair
legs”_.

We propose a novel deep learning architecture, which we call ManifoldNet, it
performs three tasks: first, it learns how to estimate the location of
different parts in the image by performing semantic segmentation, it then uses
the semantic labeling and the original input image to learn $p$ different
intermediate feature spaces for each object part and finally, $p$ different
branches of fully connected layers will learn the final image embedding into
the respective part shape manifold. The network has a general core that
performs semantic segmentation and specialized branches for each of the
manifolds in a similar fashion as [13].

Figure 6:  Detailed ManifoldNet architecture. The first set of layers take
care of assigning semantic part labels to each of the pixels in the image. The
second stage extracts intermediate features to help the manifold learning. The
final stage uses a set of fully convolutional layers to regress the manifold
coordinates.

####  4.3.1 ManifoldNet: A multi-manifold learning architecture

A summary of our new architecture is shown in Fig. 2 and a detailed
description of all the layers in Fig. 6. The common part of the architecture,
which performs the semantic segmentation, is a fully convolutional approach
closely related to [17]. The fully convolutional architecture uses a VGG-16
architecture [20] for its feature layers. A combination of FCN-32, FCN-16 and
FCN-8 is used to obtain more detailed segmentations but all sub-parts are
trained together for simplicity.

The other two parts of the architecture shown in Fig. 2 take care of: _(1)_
creating an intermediate feature space, and _(2)_ learning the manifold
embedding. The intermediate feature layers take as input the concatenation of
the original RGB image and the heat maps given as output by the semantic
segmentation layers to learn a feature representation that eases the manifold
learning task. Finally, the manifold coordinate regression module is formed by
3 fully connected layers (the first two use relu non-linearities). A dropout
scheme is used to avoid over-fitting to the data. Trained models and code are
already public but due to anonymity we cannot disclose the URL.

####  4.3.2 Details of the training strategy

The training of such a deep architecture requires careful attention. First, to
avoid vanishing gradients and to improve convergence, the semantic
segmentation layers are trained independently by using a standard cross-
entropy classification loss:

$L_{seg}=\dfrac{-1}{N}\sum_{n=1}^{N}\log{(p_{n,l_{n}})}$ (2)

where $p_{n,l_{n}}$ is the softmax output class probability for each pixel. A
batch size of only 20 is used at this stage due to memory limitations on the
GPU and the high number of weights to be trained.

When trying to train the manifold layers we found out that convergence heavily
depended on big batch sizes and many iterations. At this point we used a
learning scheme that allowed us to have bigger batch sizes during training and
faster computation of each iteration. The trick is quite simple really, we
precompute for all training images the output of the semantic layers and only
train the part branches of the network. By doing this we are training a
substantially shallower network allowing for significantly bigger batch sizes.
The network is trained by minimizing the following euclidean loss:

$L_{mani}=\dfrac{1}{2N}\sum_{i=1}^{N}{\lVert
x_{est}^{i}-x_{gt}^{i}\rVert}_{2}^{2}$ (3)

where $x_{est}^{i}$ are the manifold coordinates estimated by the network and
$x_{gt}^{i}$ are the ground truth manifold coordinates. The Euclidean loss is
chosen since the part shape manifolds are themselves Euclidean spaces. With
this good initialization of the weights we finally perform an end-to-end
training of all layers using only the final euclidean loss.

Training data and data augmentation: the training images are generated
synthetically by rendering models from $ShapeNet$ [3]. We use the 3D part
annotations on the 3D models, available from [29], to provide ground truth
values for the semantic segmentation. We generate 125 training images per
model from different poses, and a random RGB image taken from the Pascal 3D
[25] dataset is added as background. To recap, the proposed approach is
invariant to pose and manages to learn solely from rendered synthetic images.

###  4.4 Shape blending through cross-manifold optimization

Once the manifold coordinates for the different object parts have been
estimated all the information needed to blend them into a single 3D model is
available. We formulate this as a 3D shape retrieval problem: _“find the 3D
model, from the existing shapes represented in the manifolds, that best fits
the arrangement of parts”_. The user selects two (or more) images (or 3D
models) and indicates the part they wish to select from each one (note that no
annotations are needed, only the name/label of the part). The cross-manifold
optimization now finds the 3D shape in the collection that minimizes the sum
of the distances to each of the parts. In more detail – first, all manifolds
need to be normalized to allow a meaningful comparison of distances. Then,
given the set of manifold coordinates for the selected parts, a shape
prediction $b$ can be defined as the concatenation of the respective part
coordinates $b=\\{b^{1};...;b^{m}\\}$. The goal is now to retrieve a 3D model
from the shape collection whose coordinates $a=\\{a^{1};...;a^{m}\\}$ are
closest to this part arrangement by minimizing the following distance:

$B=\min_{a\in\mathbb{S}}\sum_{k=1}^{m}\lVert{a^{m}-b^{m}}\rVert$ (4)

where $\mathbb{S}$ is the set of existing shapes. Note that not all parts need
to be selected to obtain a blended shape, we define $m$ as the subset of parts
to be blended, where $m\subseteq p$. Also, notice that blending can be done by
combining any number of parts from any number of sources (shapes/images).

##  5 Results

We perform a set of qualitative and quantitative experiments to evaluate the
performance of our approach. Allthough our approach performs shape blending
from several inputs this can be understood as a retrieval task.

###  5.1 Quantitative results on Image-based Shape Retrieval

The proposed approach will be at a disadvantage when trying to retrieve whole
shapes and the same will happen when approaches that model the object as a
whole try to retrieve parts. What has been have done is a experimental
comparison of both approaches on both tasks, as a clean unbiased comparison
cannot be done on a single experiment both approaches will be used to solve
both tasks. This is possible as all compared approaches can be used as a
similarity measure between images and 3D models. By doing this, how much is
lost can be measured when either modeling the whole object or the individual
parts.

####  5.1.1 Image-based Whole Shape Retrieval

We perform whole shape image-based retrieval on the ExactMatch dataset [15].
The experiment compares against Li [15], a state-of-the-art deep volumetric
approach in Girdhar [9], HoG, AlexNet and Siamese networks. Two versions of
the proposed approach using the original three level HoG pyramid features to
build the manifold and the two level HoG manifold features that have been
shown to be better fitted for a smoother shape similarity measure. Our
approach predicts the part coordinates separately in each of the part
manifolds. The estimations of each part are then used to solve the blending
optimization and obtain a single shape prediction. The fact that the neural
network estimates the coordinates individually means that all the part co-
occurrence information that is implicitly encoded in the approaches that model
the object as a whole is lost during training, nevertheless, the proposed
approach can still yield good results that are comparable to those that model
the whole shape. It can also be seen that for the exact shape retrieval the
original three level pyramid HoG features perform better which is to be
expected. The results of these experiments are included in Fig. 7 (first
image).

Figure 7:  Left: Image-based whole shape retrieval results obtained on the
ExactMatch dataset. Center: Shape blending and retrieval from image-based part
descriptions on the ExactPartMatch dataset. Right: Matrix showing combination
results. Experiment performed by generating all possible combinations of legs
and backs from 10 test shapes.

####  5.1.2 Image-based Part Shape Blending

To test the performance of part blending and retrieval a new dataset has been
created using a subset of the shapes from the Shapenet database [3] and the
images from the ExactMatch dataset [15] to create the ExactPartMatch dataset.
The task is to find the correct 3D shape out of all the annotated 3D models
using the parts from the specified inputs. As in Shapenet many of the 3D
models are repeated many times (e.g. ikea chairs) we need to control that
there is only one correct match with the part mix. The dataset contains 845
test cases. Each test case if assigned to one of 187 hand annotated ground
truth 3D models, the same part shape combination is tested using images from
different angles, textures and lighting conditions. Each test case is a
combination of parts from two different images for which a 3$D$ model exists
in Shapenet. There are examples of all possible part combinations. The whole
dataset is already public online, the link is not disclosed due to anonymity.
The results of these experiments are included in Fig. 7 (second image). Also
in Fig. 7 (third image) we include an experiment that shows the product of
crossing 5 leg samples with 5 back samples densely, this shows the kind of
results to expect from our approach.

Figure 8:  Qualitative results extracted to show performance on different part
arrangements using parts from two different image inputs. The results show the
closest matching shapes in the Shapenet dataset. If the part is not present,
like in the second row, an $X$ in the colour of the desired part is used to
label the non-presence of that part.  Figure 9:  Qualitative results extracted
to show performance on different part arrangements using parts from three
different image inputs. The results show the closest matching shapes in the
Shapenet dataset.

For all approaches similarity is estimated in shape manifold space and then
the multi-manifold optimization is performed, which is required to obtain the
most similar shape. The approaches shown are $Ours$, with and without semantic
segmentation, $Ours$ with semantic segmentation but using the original three
level HoG features, $Li^{\prime}sSiggAsia15$ [15] and random chance. [15]
struggles to get results as good as the ones obtained by modeling the parts.
Their holistic representation enables them to better model a whole object but
loose substantial performance when trying to identify individually the parts.
In contrast our approach looses the information of part co-occurrence in
favour of being capable to model the parts individually. Also, exploiting the
semantic segmentation of the input is consistently better as it defines the
actual interest zones in the image. If considering top-5 results without
semantic segmentation $65\%$ is obtained but when using semantic segmentation
$76\%$ is obtained, which is a substantial $11\%$ improvement in performance.
It can also be seen that trying to blend shapes when using shape similarities
that do not correctly model smoothness over shape has a tremendous impact in
performance (Ours_original_HoG).

###  5.2 Qualitative results on Image-based Shape Retrieval

To further asses the quality of the results examples are shown depicting the
performance of the approach being applied to real images. In Fig. 8 and in
Fig. 9 many example images taken from the ExactPartMatch dataset detailed in
previous sections are shown. Our approach searches over the entirety of
Shapenet to find the closest matches to the input part arrangement. Many
different part arrangements are accounted for in the figure to show that the
proposed approach can capture not only the big differences but also more
subtle differences like the number of legs in the base of a specific swivel
chair, the fact that a wooden chair with a back made of bars has a round top
or a flat top, capturing the detailed shape of the interconnecting supports of
the four legs, etc.

##  6 Conclusions

An approach capable of modelling and blending object parts from images and 3D
models has been presented. This approach has demonstrated that by using a
common manifold representation very elaborate queries can be done in massive
Internet databases. It has also shown to be capable to produce accurate shape
retrieval which proves its understanding of the underlying shapes. This
provides a natural link between shape and image datasets and opens numerous
possibilities of application to similar tasks. Also, by understanding the
parts that give semantic meaning to an object volumetric approaches like [5,
12, 18, 7] could potentially address their difficulties to produce details in
their volumetric estimations by either enforcing the shape of the part
explicitly or balancing the voxel occupancy probabilities for each part
independently.

## References

  * [1] Alhashim, I., Li, H., Xu, K., Cao, J., Ma, R., Zhang, H.: Topology-varying 3d shape creation via structural blending. ACM Transactions on Graphics (2014)
  * [2] Aubry, M., Maturana, D., Efros, A., Russell, B., Sivic, J.: Seeing 3d chairs: exemplar part-based 2d-3d alignment using a large dataset of cad models. Computer Vision and Pattern Recognition (CVPR) (2014)
  * [3] Chang, A.X., Funkhouser, T., Guibas, L., Hanrahan, P., Huang, Q., Li, Z., Savarese, S., Savva, M., Song, S., Su, H., Xiao, J., Yi, L., Yu, F.: ShapeNet: An Information-Rich 3D Model Repository. Tech. Rep. arXiv:1512.03012 [cs.GR], Stanford University — Princeton University — Toyota Technological Institute at Chicago (2015)
  * [4] Chen, D.Y., Tian, X.P., Shen, Y.T., Ouhyoung, M.: On visual similarity based 3d model retrieval. Computer Graphics Forum (2003)
  * [5] Choy, C.B., Xu, D., Gwak, J., Chen, K., Savarese, S.: 3d-r2n2: A unified approach for single and multi-view 3d object reconstruction. European Conference on Computer Vision (ECCV) (2016)
  * [6] Dalal, N., Triggs, B.: Histograms of oriented gradients for human detection. In: Computer Vision and Pattern Recognition (CVPR) (2005)
  * [7] Fan, H., Su, H., Guibas, L.J.: A point set generation network for 3d object reconstruction from a single image. Computer Vision and Pattern Recognition (CVPR) (2017)
  * [8] Fish, N., van Kaick, O., Bermano, A., Cohen-Or, D.: Structure-oriented networks of shape collections. ACM Transactions on Graphics (2016)
  * [9] Girdhar, R., Fouhey, D., Rodriguez, M., Gupta, A.: Learning a predictable and generative vector representation for objects. European Conference on Computer Vision (ECCV) (2016)
  * [10] Hueting, M., Ovsjanikov, M., Mitra, N.: Crosslink: Joint understanding of image and 3d model collections through shape and camera pose variations. ACM Transactions on Graphics (Proc. SIGGRAPH Asia) (2015)
  * [11] Kalogerakis, E., Chaudhuri, S., Koller, D., Koltun, V.: A probabilistic model for component-based shape synthesis. ACM Transactions on Graphics (2012)
  * [12] Kar, A., Tulsiani, S., Carreira, J., Malik, J.: Category-specific object reconstruction from a single image. Computer Vision and Pattern Recognition (CVPR) (2015)
  * [13] Kokkinos, I.: Ubernet: Training a ‘universal’ convolutional neural network for low-, mid-, and high-level vision using diverse datasets and limited memory. In: Computer Vision and Pattern Recognition (CVPR) (2017)
  * [14] Kruskal, J.B.: Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika (1964)
  * [15] Li, Y., Su, H., Qi, C.R., Fish, N., Cohen-Or, D., Guibas, L.J.: Joint embeddings of shapes and images via cnn image purification. ACM Trans. Graph. (2015)
  * [16] Lim, I., Gehre, A., Kobbelt, L.: Identifying style of 3d shapes using deep metric learning. Computer Graphics Forum (2016)
  * [17] Long, J., Shelhamer, E., Darrell, T.: Fully convolutional networks for semantic segmentation. In: Computer Vision and Pattern Recognition (CVPR) (2015)
  * [18] Qi, C.R., Su, H., Niessner, M., Dai, A., Yan, M., Guibas, L.J.: Volumetric and multi-view cnns for object classification on 3d data. Computer Vision and Pattern Recognition (CVPR) (2016)
  * [19] Sammon, J.W.: A nonlinear mapping for data structure analysis. IEEE Transactions on Computers (1969)
  * [20] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. International Conference on Learning Representations (2015)
  * [21] Su, H., Huang, Q., Mitra, N.J., Li, Y., Guibas, L.: Estimating image depth using shape collections. ACM Transactions on Graphics (2014)
  * [22] Su, H., Qi, C.R., Li, Y., Guibas, L.J.: Render for cnn: Viewpoint estimation in images using cnns trained with rendered 3d model views. In: International Conference on Computer Vision (ICCV) (2015)
  * [23] Tasse, F.P., Dodgson, N.: Shape2vec: Semantic-based descriptors for 3d shapes, sketches and images. ACM Transactions on Graphics (2016)
  * [24] Tulsiani, S., Su, H., Guibas, L.J., Efros, A.A., Malik, J.: Learning shape abstractions by assembling volumetric primitives. Computer Vision and Pattern Recognition (CVPR) (2017)
  * [25] Xiang, Y., Mottaghi, R., Savarese, S.: Beyond pascal: A benchmark for 3d object detection in the wild. In: IEEE Winter Conference on Applications of Computer Vision (WACV) (2014)
  * [26] Xie, X., Xu, K., Mitra, N.J., Cohen-Or, D., Gong, W., Su, Q., Chen, B.: Sketch-to-Design: Context-Based Part Assembly. Computer Graphics Forum (2013)
  * [27] Xu, K., Zhang, H., Cohen-Or, D., Chen, B.: Fit and diverse: Set evolution for inspiring 3d shape galleries. ACM Transactions on Graphics (2012)
  * [28] Xu, K., Zheng, H., Zhang, H., Cohen-Or, D., Liu, L., Xiong, Y.: Photo-inspired model-driven 3d object modeling. ACM Transactions on Graphics (2011)
  * [29] Yi, L., Kim, V.G., Ceylan, D., Shen, I.C., Yan, M., Su, H., Lu, C., Huang, Q., Sheffer, A., Guibas, L.: A scalable active framework for region annotation in 3d shape collections. SIGGRAPH Asia (2016)


(ADS) and arXiv preprint repository extensively during this project and the
writing of the paper.

## Data Availability

The data products of the Ulagam simulations are publicly released and more
details can be found at https://ulagam-simulations.readthedocs.io. These
products include the density fields used in this work, alongside the requisite
scripts needed to convert them into lensing convergence fields. The halo
fields are not yet provided and will instead become public at a later time.

## References

  * Abazajian et al. (2019) Abazajian K., et al., 2019, arXiv e-prints, p. arXiv:1907.04473
  * Abbott et al. (2020) Abbott T. M. C., et al., 2020, Phys. Rev. D, 102, 023509
  * Abbott et al. (2022) Abbott T. M. C., et al., 2022, Phys. Rev. D, 105, 023520
  * Achúcarro et al. (2022) Achúcarro A., et al., 2022, arXiv e-prints, p. arXiv:2203.08128
  * Alishahiha et al. (2004) Alishahiha M., Silverstein E., Tong D., 2004, Phys. Rev. D, 70, 123505
  * Allys et al. (2020) Allys E., Marchand T., Cardoso J. F., Villaescusa-Navarro F., Ho S., Mallat S., 2020, Phys. Rev. D, 102, 103506
  * Amodeo et al. (2021) Amodeo S., et al., 2021, Phys. Rev. D, 103, 063514
  * Amon et al. (2022) Amon A., et al., 2022, Phys. Rev. D, 105, 023514
  * Anbajagane et al. (2020) Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686
  * Anbajagane et al. (2022a) Anbajagane D., Evrard A. E., Farahi A., 2022a, MNRAS, 509, 3441
  * Anbajagane et al. (2022b) Anbajagane D., et al., 2022b, MNRAS, 510, 2980
  * Anbajagane et al. (2022c) Anbajagane D., et al., 2022c, MNRAS, 514, 1645
  * Anbajagane et al. (2023a) Anbajagane D., et al., 2023a, arXiv e-prints, p. arXiv:2308.03863
  * Anbajagane et al. (2023b) Anbajagane D., et al., 2023b, arXiv e-prints, p. arXiv:2310.00059
  * Angulo & Hahn (2022) Angulo R. E., Hahn O., 2022, Living Reviews in Computational Astrophysics, 8, 1
  * Aricò et al. (2023) Aricò G., Angulo R. E., Zennaro M., Contreras S., Chen A., Hernández-Monteagudo C., 2023, arXiv e-prints, p. arXiv:2303.05537
  * Armendáriz-Picón et al. (1999) Armendáriz-Picón C., Damour T., Mukhanov V., 1999, Physics Letters B, 458, 209
  * Asgari et al. (2021) Asgari M., et al., 2021, A&A, 645, A104
  * Baldauf et al. (2015) Baldauf T., Mercolli L., Zaldarriaga M., 2015, Phys. Rev. D, 92, 123007
  * Banerjee & Abel (2023) Banerjee A., Abel T., 2023, MNRAS, 519, 4856
  * Banerjee et al. (2022) Banerjee A., Kokron N., Abel T., 2022, MNRAS, 511, 2765
  * Barreira (2022) Barreira A., 2022, J. Cosmology Astropart. Phys., 2022, 013
  * Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
  * Baumann et al. (2012) Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmology Astropart. Phys., 2012, 051
  * Bhattacharya et al. (2011) Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122
  * Blazek et al. (2019) Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Brinch Holm et al. (2023) Brinch Holm E., Herold L., Simon T., Ferreira E. G. M., Hannestad S., Poulin V., Tram T., 2023, arXiv e-prints, p. arXiv:2309.04468
  * Cabass et al. (2022a) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022a, Phys. Rev. D, 106, 043506
  * Cabass et al. (2022b) Cabass G., Ivanov M. M., Philcox O. H. E., Simonović M., Zaldarriaga M., 2022b, Phys. Rev. Lett., 129, 021301
  * Carrasco et al. (2012) Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, Journal of High Energy Physics, 2012, 82
  * Chen (2010) Chen X., 2010, Advances in Astronomy, 2010, 638979
  * Chen et al. (2023) Chen A., et al., 2023, MNRAS, 518, 5340
  * Cheng & Ménard (2021) Cheng S., Ménard B., 2021, MNRAS, 507, 1012
  * Cheung et al. (2008) Cheung C., Fitzpatrick A. L., Kaplan J., Senatore L., Creminelli P., 2008, Journal of High Energy Physics, 2008, 014
  * Chisari et al. (2018) Chisari N. E., et al., 2018, MNRAS, 480, 3962
  * Comparat et al. (2017) Comparat J., Prada F., Yepes G., Klypin A., 2017, MNRAS, 469, 4157
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * Costanzi et al. (2021) Costanzi M., et al., 2021, Phys. Rev. D, 103, 043522
  * Coulton et al. (2022) Coulton W. R., et al., 2022, arXiv e-prints, p. arXiv:2206.01619
  * Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369
  * Cui et al. (2022) Cui W., et al., 2022, MNRAS, 514, 977
  * D’Amico et al. (2022) D’Amico G., Lewandowski M., Senatore L., Zhang P., 2022, arXiv e-prints, p. arXiv:2201.11518
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * Dalal et al. (2008) Dalal N., Doré O., Huterer D., Shirokov A., 2008, Phys. Rev. D, 77, 123514
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Duffy et al. (2010) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161
  * Dvali et al. (2004a) Dvali G., Gruzinov A., Zaldarriaga M., 2004a, Phys. Rev. D, 69, 023505
  * Dvali et al. (2004b) Dvali G., Gruzinov A., Zaldarriaga M., 2004b, Phys. Rev. D, 69, 083505
  * Eifler et al. (2021) Eifler T., et al., 2021, MNRAS, 507, 1746
  * Emami et al. (2015) Emami R., Dimastrogiovanni E., Chluba J., Kamionkowski M., 2015, Phys. Rev. D, 91, 123531
  * Enqvist & Sloth (2002) Enqvist K., Sloth M. S., 2002, Nuclear Physics B, 626, 395
  * Euclid Collaboration et al. (2019) Euclid Collaboration et al., 2019, MNRAS, 484, 5509
  * Fluri et al. (2018) Fluri J., Kacprzak T., Refregier A., Amara A., Lucchi A., Hofmann T., 2018, Phys. Rev. D, 98, 123518
  * Fluri et al. (2019) Fluri J., Kacprzak T., Lucchi A., Refregier A., Amara A., Hofmann T., Schneider A., 2019, Phys. Rev. D, 100, 063514
  * Fluri et al. (2022) Fluri J., Kacprzak T., Lucchi A., Schneider A., Refregier A., Hofmann T., 2022, Phys. Rev. D, 105, 083518
  * Gatti et al. (2020) Gatti M., et al., 2020, MNRAS, 498, 4060
  * Gatti et al. (2022) Gatti M., et al., 2022, Phys. Rev. D, 106, 083509
  * Gatti et al. (2023) Gatti M., et al., 2023, arXiv e-prints, p. arXiv:2307.13860
  * Giannantonio et al. (2012) Giannantonio T., Porciani C., Carron J., Amara A., Pillepich A., 2012, MNRAS, 422, 2854
  * Giri & Schneider (2021) Giri S. K., Schneider A., 2021, J. Cosmology Astropart. Phys., 2021, 046
  * Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
  * Guth (2004) Guth A. H., 2004, in Freedman W. L., ed., Measuring and Modeling the Universe. p. 31 (arXiv:astro-ph/0404546), doi:10.48550/arXiv.astro-ph/0404546
  * Halder et al. (2021) Halder A., Friedrich O., Seitz S., Varga T. N., 2021, MNRAS, 506, 2780
  * Hartlap et al. (2007) Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399
  * Hilbert et al. (2012) Hilbert S., Marian L., Smith R. E., Desjacques V., 2012, MNRAS, 426, 2870
  * Hill et al. (2018) Hill J. C., Baxter E. J., Lidz A., Greco J. P., Jain B., 2018, Phys. Rev. D, 97, 083501
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Jarvis et al. (2004) Jarvis M., Bernstein G., Jain B., 2004, MNRAS, 352, 338
  * Jung et al. (2023) Jung G., et al., 2023, arXiv e-prints, p. arXiv:2305.10597
  * Kacprzak et al. (2023) Kacprzak T., Fluri J., Schneider A., Refregier A., Stadel J., 2023, J. Cosmology Astropart. Phys., 2023, 050
  * Kaiser & Squires (1993) Kaiser N., Squires G., 1993, ApJ, 404, 441
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Kofman (2003) Kofman L., 2003, arXiv e-prints, pp astro–ph/0303614
  * Kokron et al. (2022) Kokron N., DeRose J., Chen S.-F., White M., Wechsler R. H., 2022, MNRAS, 514, 2198
  * Krause & Hirata (2010) Krause E., Hirata C. M., 2010, A&A, 523, A28
  * Krause et al. (2021) Krause E., et al., 2021, arXiv e-prints, p. arXiv:2105.13548
  * Lamman et al. (2023) Lamman C., Tsaprazi E., Shi J., Niko Šarčević N., Pyne S., Legnani E., Ferreira T., 2023, arXiv e-prints, p. arXiv:2309.08605
  * Lee et al. (2022) Lee E., et al., 2022, MNRAS, 517, 5303
  * Lee et al. (2023) Lee M. E., Lu T., Haiman Z., Liu J., Osato K., 2023, MNRAS, 519, 573
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv e-prints, p. arXiv:1104.2932
  * Lim et al. (2021) Lim S. H., Barnes D., Vogelsberger M., Mo H. J., Nelson D., Pillepich A., Dolag K., Marinacci F., 2021, MNRAS, 504, 5131
  * LoVerde et al. (2008) LoVerde M., Miller A., Shandera S., Verde L., 2008, J. Cosmology Astropart. Phys., 2008, 014
  * Lyth & Wands (2002) Lyth D. H., Wands D., 2002, Physics Letters B, 524, 5
  * MacCrann et al. (2022) MacCrann N., et al., 2022, MNRAS, 509, 3371
  * Marian et al. (2011) Marian L., Hilbert S., Smith R. E., Schneider P., Desjacques V., 2011, ApJ, 728, L13
  * McKinney (2011) McKinney W., 2011, Python for High Performance and Scientific Computing, 14
  * Mead & Verde (2021) Mead A. J., Verde L., 2021, MNRAS, 503, 3095
  * Modi et al. (2020) Modi C., Chen S.-F., White M., 2020, MNRAS, 492, 5754
  * More et al. (2023) More S., et al., 2023, arXiv e-prints, p. arXiv:2304.00703
  * Moroi & Takahashi (2001) Moroi T., Takahashi T., 2001, Physics Letters B, 522, 215
  * Mueller et al. (2021) Mueller E.-M., et al., 2021, arXiv e-prints, p. arXiv:2106.13725
  * Munshi et al. (2022) Munshi D., Lee H., Dvorkin C., McEwen J. D., 2022, J. Cosmology Astropart. Phys., 2022, 020
  * Myles et al. (2021) Myles J., et al., 2021, MNRAS, 505, 4249
  * Pandey et al. (2022) Pandey S., et al., 2022, Phys. Rev. D, 105, 123526
  * Petri et al. (2017) Petri A., Haiman Z., May M., 2017, Phys. Rev. D, 95, 123503
  * Philcox et al. (2022a) Philcox O. H. E., Ivanov M. M., Cabass G., Simonović M., Zaldarriaga M., Nishimichi T., 2022a, Phys. Rev. D, 106, 043530
  * Philcox et al. (2022b) Philcox O. H. E., Slepian Z., Hou J., Warner C., Cahn R. N., Eisenstein D. J., 2022b, MNRAS, 509, 2457
  * Pillepich et al. (2010) Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020a) Planck Collaboration et al., 2020a, A&A, 641, A9
  * Planck Collaboration et al. (2020b) Planck Collaboration et al., 2020b, A&A, 641, A10
  * Potter et al. (2017) Potter D., Stadel J., Teyssier R., 2017, Computational Astrophysics and Cosmology, 4, 2
  * Rodríguez-Monroy et al. (2022) Rodríguez-Monroy M., et al., 2022, MNRAS, 511, 2665
  * Sasaki et al. (2006) Sasaki M., Väliviita J., Wands D., 2006, Phys. Rev. D, 74, 103003
  * Schäfer et al. (2012) Schäfer B. M., Grassi A., Gerstenlauer M., Byrnes C. T., 2012, MNRAS, 421, 797
  * Schmidt & Kamionkowski (2010) Schmidt F., Kamionkowski M., 2010, Phys. Rev. D, 82, 103002
  * Schneider et al. (2016) Schneider A., et al., 2016, J. Cosmology Astropart. Phys., 2016, 047
  * Schneider et al. (2019) Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmology Astropart. Phys., 2019, 020
  * Scoccimarro et al. (2012) Scoccimarro R., Hui L., Manera M., Chan K. C., 2012, Phys. Rev. D, 85, 083002
  * Secco et al. (2022a) Secco L. F., et al., 2022a, Phys. Rev. D, 105, 023515
  * Secco et al. (2022b) Secco L. F., et al., 2022b, Phys. Rev. D, 105, 103537
  * Sefusatti et al. (2010) Sefusatti E., Crocce M., Desjacques V., 2010, MNRAS, 406, 1014
  * Senatore et al. (2010) Senatore L., Smith K. M., Zaldarriaga M., 2010, J. Cosmology Astropart. Phys., 2010, 028
  * Sevilla-Noarbe et al. (2021) Sevilla-Noarbe I., et al., 2021, ApJS, 254, 24
  * Shao et al. (2022) Shao M., Anbajagane D., Chang C., 2022, arXiv e-prints, p. arXiv:2212.05964
  * Shirasaki et al. (2012) Shirasaki M., Yoshida N., Hamana T., Nishimichi T., 2012, ApJ, 760, 45
  * Silverstein & Tong (2004) Silverstein E., Tong D., 2004, Phys. Rev. D, 70, 103505
  * Simon et al. (2009) Simon P., Taylor A. N., Hartlap J., 2009, MNRAS, 399, 48
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Stiskalek et al. (2022) Stiskalek R., Bartlett D. J., Desmond H., Anbajagane D., 2022, MNRAS, 514, 4026
  * Sunseri et al. (2023) Sunseri J., Slepian Z., Portillo S., Hou J., Kahraman S., Finkbeiner D. P., 2023, RAS Techniques and Instruments, 2, 62
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Tegmark & Zaldarriaga (2002) Tegmark M., Zaldarriaga M., 2002, Phys. Rev. D, 66, 103508
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration 2005, arXiv e-prints, pp astro–ph/0510346
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration et al., 2018, arXiv e-prints, p. arXiv:1809.01669
  * Tinker et al. (2010) Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
  * Troxel & Ishak (2015) Troxel M. A., Ishak M., 2015, Phys. Rep., 558, 1
  * Van der Walt et al. (2011) Van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro F., et al., 2020, ApJS, 250, 2
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
  * Zhang et al. (2022) Zhang Z. J., Chang C., Larsen P., Secco L. F., Zuntz J., LSST Dark Energy Science Collaboration 2022, MNRAS, 514, 2181
  * Zürcher et al. (2021) Zürcher D., Fluri J., Sgier R., Kacprzak T., Refregier A., 2021, J. Cosmology Astropart. Phys., 2021, 028
  * Zürcher et al. (2022) Zürcher D., et al., 2022, MNRAS, 511, 2075

##  Appendix A Validation of simulations

We validate here the Ulagam simulation suite through comparisons with
theoretical predictions from other models (either based on perturbation theory
and/or from simulation-based emulators) and comparisons with the original
Quijote suite whose initial conditions are used in the Ulagam runs.

Figure 9: The density fields of a Ulagam and Quijote simulation with the same
initial condition. The slices have $500\,{\rm Mpc}/h$ thickness. The density
distribution is shown in the top right and the power spectra in the bottom
right. White circles show the location of the 10 most massive halos in each
slice. There is clear agreement between our simulations and the quijote runs.
The specific ordering of the 10 most massive halos varies slightly. The
distribution of $1+\delta$ is within $1\%$ for $\delta>1$, and grows to $10\%$
for $\delta\approx 0$. The power spectra agree to $1\%$ at z = 0 and grow to
$2-3\%$ at higher redshifts, which is within the expected deviations due to
differences in N-body solvers (Schneider et al. 2016)

Comparing Ulagam to Quijote. Figure 9 compares the properties of the $z=0$
three-dimensional density field between a single run in the two simulation
suites. While all the Ulagam simulations save lightcone maps and not three-
dimensional snapshots, we ran a single simulation that saved the snapshot
information as well, specifically to perform this comparison. The left panels
of the figure show $500\,{\rm Mpc}/h$ thick slices of the density field
projected along different axes (the cartesian axes of the projected field are
denoted in the subplot title). The top row shows the Ulagam realization, run
with Pkdgrav3, while the bottom row shows the Gadget3-based Quijote
realization. Brighter (darker) regions denote overdensities (underdensities).
A simple visual comparison shows that the structure in both runs is
consistently realized. The white circles denote the location of the top 10
most massive halos in either realization. Most of the locations are common
with some minor differences due to the mass ordering of the halos not being
completely preserved. This difference in mass-ordering is expected as the
halo-finding procedure can be sensitive to stochastic noise. While such
effects are suppressed in our comparison given both realizations start from
the same initial conditions, differences can still arise due to the use of
different gravity solvers with different numerical noise (see Schneider et al.
2016, for comparison of power spectra), and due to different FoF
implementations (Knebe et al. 2011, see their Figure 5).

The right panels of Figure 9 shows the probability distribution function of
the overdensity field, where the Ulagam and Quijote realizations are within
$1\%$ over the vast range of $\delta$ values, and the difference grows to
$10\%$ only for $\delta\approx 0$, near the tails of the distribution. The
bottom right panel shows the comparison of the power spectra for the two runs
at different redshifts. We once again find that the deviations are minimal;
they are within $1\%$ at $z=0$ and grow to $2-3\%$ at higher redshift. Such
differences, at our resolution level of $512^{3}$ particles, are consistent
with expectations from comparison studies of the Pkdgrav3 and Gadget3
(Schneider et al. 2016).

Figure 10: Residuals of matter power spectrum, $\delta P(k)=P_{\rm
fid}(k)/P_{\rm truth}(k)-1$, for five redshifts (shown in top right) where the
truth is either computed with Class \+ Halofit, Class \+ Euclid-Emu or a
lower/higher resolution run with $256^{3}$/$1024^{3}$ particles. The fiducial
resolution run uses $512^{3}$ particles. The dark (light) gray band show the
1% (5%) residuals. The power spectra match within $5\%$ up to $k=1$ Mpc-1 at z
= 0, and $k=0.3$ Mpc-1 at z = 2. Increasing the resolution of the simulation
makes the resulting power spectra agree better with both Halofit and Euclid-
Emu.

Power Spectra. We then validate the simulations against a number of
theoretical models, starting with the power spectrum of the density field.
Figure 10 shows the fractional deviation between the Ulagam suite (from the
average of five realizations) alongside four different theoretical
predictions. Two are simply the average of five lower/higher resolutions runs
from the Ulagam suite, while the other two are predictions from Halofit
(Takahashi et al. 2012) and Euclid-Emu (Euclid Collaboration et al. 2019).
Both models take the perturbation theory result from the Class boltzmann code
(Lesgourgues 2011) and modify it to model the non-linear regime. The agreement
with all models (other than the lower resolution run) is within $5\%$ up to
$k=0.3\,{\rm Mpc}^{-1}$ at $z=2$, and up to $k=1\,{\rm Mpc}^{-1}$ at $z=0$.
While the deviations with Euclid-Emu increase towards high $k$, improving the
simulation resolution to $1024^{3}$ particles results in significantly better
agreement. The Euclid-Emu model was built using Pkdgrav3 simulations, but run
in a larger volume of $L=1.9\,{\rm Gpc}$ and with $2048^{3}$ particles. The
Halofit comparison shows oscillatory residuals, particularly at $z=0$, and
this was identified in Euclid Collaboration et al. (2019, see their Figure 8)
as a systematic effect in Halofit from not accurately capturing the baryon
acoustic oscillation features.

Figure 11: The weak lensing convergence power spectra measured on the full-sky
(averaged over all 2000 simulations available at the fiducial cosmology),
compared with theoretical predictions from Class, modified by the Halofit
prescription for the non-linear regime. The dark (light) gray bands show the
1% (5%) error regime. The colored lines show the comparison for different
redshifts. The thin black lines show the cosmic variance amplitude at $z=0.4$.
The result for other redshifts is fairly similar, and so we do not show them
for brevity. The deviations are within $1\%$ for a wide range of multipoles
and redshifts, and grow towards low/high multipoles.

Lensing harmonic power spectra. Figure 11 validates the weak lensing
convergence power spectra on the full-sky. We show the comparisons for five
source planes at redshifts that span the width of the weak lensing kernel
shown in Figure 1. The theoretical model is obtained from Class, modified by
the Halofit prescription for the non-linear regime. The dark (light) bands
show the 1% (5%) error range. The harmonic spectra are the averages of all
2000 full-sky maps, measured by binning the $C_{\ell}$ in 50 log-space
multipole bins between $10<\ell<2048$. The deviations are within $1\%$ for a
wide range of $\ell$, and increase towards low and high $\ell$.

Figure 12: The halo mass function, averaged over 100 realizations of the
fiducial cosmology, compared to theoretical predictions from Watson et al.
(2013) for halos found through friends-of-friends. The agreement is 5% for
most redshifts/masses and increases to 10% for the highest mass objects at
each redshift.

Halo mass function. We also use the halo abundance in Section 4.3 to
understand the effect of different $f_{\rm NL}$. Figure 12 shows the HMF
averaged over 100 realizations, for different redshifts. The lower mass limit
of $10^{14}\,{\rm M}_{\odot}$ is determined by requiring atleast 100 particles
per halo. Overplotted is the theoretical model from Watson et al. (2013) for
FoF-based halos, computed using the Colossus package (Diemer 2018). The
agreement is 5% for most redshifts/masses and increases to 10% for the highest
mass objects at each redshift.

Figure 13: The halo bias, as a function of angular multipole scale and
redshift, estimated using Equation 4.9. We show the average result from 100
realizations at the fiducial cosmology. At large scales (low $\ell$), the bias
asymptotes to a constant value, often called the “linear bias”. The vertical
bands show the upper and lower limits from a range of bias predictions
(Pillepich et al. 2010; Tinker et al. 2010; Bhattacharya et al. 2011; Comparat
et al. 2017). The horizontal location of the bands is chosen for visualization
purposes and has no other information given the bias is scale-independent on
sufficiently large scales. At high $\ell$, there is a clear scale dependence
as found in previous work (e.g., Mead & Verde 2021).

Halo bias. Finally, we show in Figure 13 the halo bias as a function of
multipole, estimated using Equation 4.9. At large scales, where the bias
asymptotes to a constant value, we show the range of predictions from
different theoretical models as vertical lines. In specific, we compare the
models from Pillepich et al. (2010); Tinker et al. (2010); Bhattacharya et al.
(2011); Comparat et al. (2017), and all predictions are computed using
Colossus (Diemer 2018). The exact upper/lower limits used in Figure 13 are
from Comparat et al. (2017) and Tinker et al. (2010) respectively. There is a
good agreement between the predictions and the measurements across all
redshifts. Note that the bias increases with redshift and this is due to the
fixed mass cut across all redshifts. For $M>10^{14}\,{\rm M}_{\odot}$, the
halo sample at $z=1.2$ is a rarer sample than that at $z=0$, which then leads
to a higher bias. At high $\ell$, the bias grows towards small scales, as has
been found in previous work (e.g., Mead & Verde 2021).

##  Appendix B Numerical convergence of Fisher Information

A significant concern in numerical estimates of the Fisher information is the
artificial amplification of the information due to numerical noise. This has
been documented in previous studies with the Quijote suite (e.g., Coulton et
al. 2022; Jung et al. 2023). A simple test of this effect is to vary the
number of simulations used in computing either the derivatives or the
covariance matrix as used in Equation 4.1. In Figure 14 we present the
constraints on $f_{\rm NL}$, normalized by the fiducial constraints obtained
using all available simulations, as a function of the number of realizations
used to compute the derivatives and the covariance. We show only the LSST Y10
results, corresponding to the unmarginalized constraints shown in Figure 2.
The constraints for other surveys and for marginalized analyses show the same
convergence behaviors as the ones described below. The constraints from the
moments (up to 5th order) are well-converged, as an increase in 100
realizations for the derivative calculation results in less than 1% changes in
the constraints. The covariance is also well-converged. Doubling the number of
simulations used to estimate the covariance results in less than 2%
differences in the final constraints. Of all $f_{\rm NL}$ studied here,
$f_{\rm NL}^{\rm\,or,\,cmb}$ is the most poorly converged, and even then is
converged to within 10% for the combination of 2nd and 3rd moments, which is
the fiducial statistic we use in this work.

Figure 14: The change in the Fisher information for unmarginalized $f_{\rm
NL}$, for an LSST Y10-like survey, as a function of number of simulations used
to estimate the derivative (top) or the covariance (bottom). We show all
statistics presented in Figure 2. The combination of 2nd and 3rd moments — the
fiducial statistic used in this work — is converged to within 1%. The
convergence properties are similar for the other surveys and for the
marginalized $f_{\rm NL}$ constraints.

##  Appendix C Full constraints for marginalized $f_{\rm NL}$

Figure 15: The full Fisher constraints for $f_{\rm NL}^{\rm\,eq}$,
$\Omega_{\rm m}$, $\sigma_{8}$, $A_{\rm IA}$ and $\eta_{\rm IA}$. The dashed
contours show constraints over all five paramters, while solid contours are
constraints varying $f_{\rm NL}$ and cosmology parameters alone. There is only
marginal differences between the two. The 1D posterior for $f_{\rm
NL}^{\rm\,eq}$ also shows the unmarginalized constraint in the dotted lines.

For completeness, we show the full triangle plot of constraints for the
marginalized $f_{\rm NL}^{\rm\,eq}$ analysis. The $\sigma_{8}-f_{\rm
NL}^{\rm\,eq}$ plane shows an anti-correlation; increasing $f_{\rm
NL}^{\rm\,eq}$ leads to more structure on small scales, which can be partially
compensated by a reduction in $\sigma_{8}$ as that removes structures across a
wide range of scales. A similar anti-correlation is found in $\Omega_{\rm
m}-f_{\rm NL}^{\rm\,eq}$. However, this is only for the analysis of the 2nd
moments; the $\Omega_{\rm m}-f_{\rm NL}^{\rm\,eq}$ correlation is positive for
analyses combining the 2nd and 3rd moments.

We do not show the full constraints from the other parameters for brevity.
Note that in this work we do not vary all $f_{\rm NL}$ parameters at once, and
instead analyze them one at a time; the sole exception is the estimate of the
sound speed, $c_{s}$, in Section 4.2 where we vary both $f_{\rm NL}^{\rm\,eq}$
and $f_{\rm NL}^{\rm\,or,\,lss}$. The constraints for $f_{\rm NL}^{\rm\,loc}$
have qualitatively similar degeneracies in all parameter planes as the ones
described above. The degeneracies of $f_{\rm NL}^{\rm\,or,\,lss}$ and $f_{\rm
NL}^{\rm\,or,\,cmb}$ with IA and cosmological parameters are qualitatively
different, and this can be inferred from the behavior of the halo mass
function in Figure7. In particular, the sign of the correlation between
$f_{\rm NL}$ and the cosmology parameters is reversed in every plane. Figure 7
shows that the orthogonal-type $f_{\rm NL}$ reduce the number of high-density
peaks — which is also supported by the analysis of Coulton et al. (2022),
which found the change in the power spectrum is also negative — and therefore,
the sign of the correlation will be flipped.


This document contains a concise and unified reference for one of the existing mechanizations of the UD Kalman filter. The associated matrix algorithms are also included along with the corresponding references.

§ INTRODUCTION

The UD formulation of the Kalman filter does not operate on the covariance matrix $\Cov$ directly. Instead, it updates and propagates the covariance factors $\Umat$ and $\D$ ($\Cov=\Umat\D\Umat\trans$). This formulation is preferred for hardware implementation as it provides an efficient and numerically stable solution while ensuring the symmetry of the covariance matrix even for large filters <cit.>. In contrast with Potter's square root filter <cit.> that uses Cholesky factors $\Smat$, such that $\Cov = \Smat\Smat\trans $, the UD formulation does not necessarily increase the numerical precision of the filter <cit.>, but it is generally preferred for its relatively high efficiency.

The following sections summarize the UD Kalman filter, and present the details on its implementation. This independent document is also intended to be a companion paper to a more general UD filter presented in <cit.>.

§ THE UD FACTORIZATION IN KALMAN FILTERING
The UD factors can be obtained as a corollary of the Cholesky square-root factors $\Smat$ used in Potter's filter (<cit.>, Chap. 10). A closer examination of the Cholesky decomposition of an $n\times n$ matrix, reveals that $n$ square root operations are computed and that those same square-roots appear dividing each column of the matrix, motivating a pair of alternative Cholesky factors: the UD factors. To illustrate this and further see how the UD factors can be obtained, consider a $ 3 \times 3 $ covariance matrix, $ \Cov $, to be analytically factorized via Cholesky decomposition as

\begin{equation}
\Cov = \Smat\Smat\trans \
\end{equation}

with $ \Smat $ written as

\begin{equation}
\Smat = \begin{bmatrix}
\sqrt{\Cov_{11}- \frac{\Cov_{13}^2}{\Cov_{33}}-\frac{(\Cov_{12}-\frac{\Cov_{13}\Cov_{23}}{\Cov_{33}})^2}{\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}}} &
\frac{\Cov_{12}-\frac{\Cov_{13}\Cov_{23}}{\Cov_{33}}} {\sqrt{\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}}}& \frac{\Cov_{13}}{{\sqrt{\Cov_{33}}}}\\
0 		&	\sqrt{\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}}	& \frac{\Cov_{23}}{{\sqrt{\Cov_{33}}}}\\
0 & 0 & \sqrt{\Cov_{33}}
\end{bmatrix} \
\end{equation}

which can further be factorized as

\begin{equation}
\Smat	= \begin{bmatrix}
1 & \frac{\Cov_{12}-\frac{\Cov_{13}\Cov_{23}}{\Cov_{33}}} {\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}} & \frac{\Cov_{13}}{\Cov_{33}} \\
0 & 1 & \frac{\Cov_{23}}{\Cov_{33}} \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\sqrt{\Cov_{11}- \frac{\Cov_{13}^2}{\Cov_{33}}-\frac{(\Cov_{12}-\frac{\Cov_{13}\Cov_{23}}{\Cov_{33}})^2}{\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}}} & 0 & 0 \\
0 & \sqrt{\Cov_{22}-\frac{\Cov_{23}^2}{\Cov_{33}}} & 0 \\
0  & 0 & \sqrt{\Cov_{33}}
\end{bmatrix} \
\end{equation}

Defining the upper triangular matrix with ones on the diagonal as $ \Umat $, and the diagonal matrix as $\sqrt{\D} $, the square root factor $\Smat $ can now be written as

\begin{equation}
\Smat = \Umat\sqrt{\D} \
\end{equation}

which provides the alternative covariance matrix factorization:

\begin{equation}
\Cov=	 \Smat\Smat\trans = \Umat\sqrt{\D} \sqrt{\D}\Umat\trans = \Umat \D\Umat\trans \
\end{equation}

Because the square root terms are contained in the full representation of $\D$, and the filter propagates and updates the matrices $\Umat$ and $\D$, no square root operations are executed in the UD Kalman filter. Moreover, when a Kalman filter is in UD form, monitoring the covariance positive semi-definiteness only requires checking the sign of the elements of $ \D $, and enforcement of its positive semi-definiteness can be done by ensuring they are non-negative <cit.>.
A mechanization to obtain the UD factors is included here in Appendix <ref>. But many other sources are also available in the filtering literature (<cit.>, Chap. 8), (<cit.>, p. 174), (<cit.>, Chap. 7), in linear algebra books, and in matrix operations books (<cit.>, p. 563), (<cit.>, Chap. 5). Further elaboration on this topic, and other factorized variants of Kalman filter can be found in <cit.>, <cit.>.

§.§ UD Kalman Filter Propagation step

There are different alternatives to perform a time update using the UD filter. This document presents the formulation from (<cit.>, p. 176) and <cit.> as this form is common in the literature, and it does not assume a particular structure of the system state-space. Alternative UD factor propagation approaches can be found in <cit.>.

To illustrate the idea behind the propagation of the UD factors within the Extended Kalman filter (EKF) framework (<cit.>, Chap. 5), consider the following discrete nonlinear dynamic system with state vector $\x_k \in \Real^n$, known input $\uk \in \Real^r$, measurement vector  $\yktilde \in \Real^m$, process noise $ \w_k \in \Real^q$, and measurement noise $\v_k \in \Real^m$:

\begin{equation}\label{eq:nonlinear_equations}
\x_{k}=\f_{k-1}(\x_{k-1},\u_{k-1},\w_{k-1})
\end{equation}

\begin{equation}
\yktilde=\hkofxk+\vk
\end{equation}

\begin{equation}
\wk \sim \mathcal{N} (\zerovec,\Qmat_k)
\end{equation}

\begin{equation}\label{eq:measurement_cov}
\vk \sim \mathcal{N} (\zerovec,\R_k)
\end{equation}

Here, $\wk$ and $\vk $ are zero-mean Gaussian white-noise processes, with covariances $\Qmat_k=\expect[\wk\wk\trans]$ and $\R_k=\expect[\vk\vk\trans]$, respectively; the function $\hkofxk$ is the measurement model function, and all of the sub-indices denote the time instance. The UD EKF and the conventional EKF propagate states similarly via a nonlinear or linearized process model. Finally, consider the EKF discrete covariance propagation equation given as

\begin{equation}\label{eq:covariance_propagation}
\stdvec{P}[][k][-] = \Fmat_{k-1}\stdvec{P}[][k-1][+]\Fmat_{k-1}\trans+\Gmat_{k-1}\Qmat_{k-1}\Gmat_{k-1}\trans
\end{equation}

Here, $\stdvec{P}_k^- $ is the prior error covariance, $\stdvec{P}[][k-1][+]$ is the most recent updated covariance, $ \Gmat_{k-1} = \partder{\f_{k-1}}{\w}
$ is the $ n \times q $ matrix mapping the process noise from the vector $ \wk$ to the state, and $\Fmat_{k-1} = \partder{\f_{k-1}}{\x}$ is the Jacobian of the process. Both, $\Fmat_{k-1}$ and $\Gmat_{k-1}$, are evaluated at $\expect[\x_{k-1}],\ \expect[\w_{k-1}]=0,$ and $ \u_{k-1}$.

The idea of the UD propagation stage is to operate on a factorized form of $\stdvec{P}[][k][-] = \Fmat_{k-1}\stdvec{P}[][k-1][+]\Fmat_{k-1}\trans+\Gmat_{k-1}\Qmat_{k-1}\Gmat_{k-1}\trans $ composed by three factors that resemble the form $ \Cov^-_k = \Usub_k\Dsub_k\Usub_k\trans$ (where the upper bar in a quantity, i.e., $ \bar{[\,\cdot\,]}$, indicates a prior quantity). A direct attempt to factorize $\stdvec{P}[][k][-] = \Fmat_{k-1}\stdvec{P}[][k-1][+]\Fmat_{k-1}\trans+\Gmat_{k-1}\Qmat_{k-1}\Gmat_{k-1}\trans $ into the three factors, $\Wmat\Dhat\Wmat\trans$, leads to the candidate form of (with no time indices for convenience)

\begin{equation}\label{eq:P_wrong_dimensions}
\stdvec{P}^- = \Usub\Dsub\Usub\trans= \rowvec{\Fmat\Uplus,\Gmat}\begin{bmatrix}
\Dplus & \zerovec \\ \zerovec & \Qmat
\end{bmatrix}\colvec{\Uplus\trans\Fmat\trans,\Gmat\trans}= \Wmat\Dhat\Wmat\trans \
\end{equation}

where the notation $ \overset{+}{[\,\cdot\,]}$ indicates posterior quantities; in this case, the most recent posterior quantities. At this point
the propagated $\Usub$ and $\Dsub$ factors from Eq. (<ref>) cannot be obtained via a direct term-by-term comparison (i.e. $\Usub\neq\Wmat$, $\Dsub\neq\Dhat$) because the diagonal matrix $ \Dhat $ defined above is an $ (n+q)\times(n+q) $ matrix, and $ \Wmat = [\Fmat\Uplus \ \Gmat] $ is an $ n\times (n+q)$ matrix and is not upper triangular in general. Nevertheless, some work can be done on $\Wmat$ and $\Dhat$ to rewrite them in the appropriate form and dimension such that the propagated $\Usub$ and $\Dsub$ can be obtained via a direct comparison. That is, one seeks to have

\begin{equation}\label{eq:udu_equal_wdw}
\Cov^-=\Usub{\Dsub}\Usub\trans = \tilde{\Wmat}\tilde{\D}\tilde{\Wmat}\trans  \ \end{equation}

such that $ \Usub=\tilde{\Wmat} $ and $ \Dsub=\tilde{\D}$, where $\tilde{\Wmat}$ and $\tilde{\D}$ are alternative factors to be obtained. In an abuse of notation, the tilde on the $\Wmat$ and $\D$ matrices in Eq. (<ref>) are to emphasize the difference from $\Wmat$ and $\Dhat$ in Eq. (<ref>).

The upper triangular matrix $ \tilde{\Wmat} $, and the diagonal matrix $ \tilde{\D} $, can be obtained with the execution of a Weighted Modified Gram-Schmidt (WMGS) orthogonalization, given the factorized candidate form: the rows $\w_i$ of $ \Wmat$ and the elements of $\Dhat$ (<cit.>, Chap. 6). To show this, consider the WMGS procedure that consists of the following recursions:

\begin{align}
\v_n&= \w_n \  \\\label{eq:v_vectors}
\v_k &= \w_k - \sum_{j=k+1}^{n}u(k,j)\v_j \quad k = n-1,\dots,1 \ \\
u(k,j) &= \frac{\w_k\Dhat\v_j^T}{\v_j\Dhat\v_j\trans}	\quad j,k = 1,\dots, n \
\end{align}

\begin{equation}\label{eq:vDv_product}
\v_k\Dhat\v_j^T  = 0 \quad  \text{when}\  k\neq j \
\end{equation}

where the $ \v$'s are row vectors of size $ (n+q)$.
Then, by recognizing that $_k$ can be alternatively expressed as,
\begin{equation}\label{eq:alternative_expression_w}
\w_k=\v_k +  \sum_{j=k+1}^{n}u(k,j)\v_j \quad k=1, \dots,n \
\end{equation}
one can
% 	\begin{equation}
% 	\w_k^T=\v_k^T +  \sum_{j=k+1}^{n}u(k,j)\v_j^T \quad k=1, \dots,n \ ,
% 	\end{equation}
stack all $n$ vectors $$ in matrix form as
\begin{equation}
\begin{bmatrix}
\w_1\\
\vdots\\
\w_n
\end{bmatrix}=
\begin{bmatrix}
1 & u(1,2) & \dots & u(1,n) \\
0 & 1 & \ddots& \vdots\\
\vdots & \ddots & \ddots & u(n-1,n)\\
0 & \dots & \dots & 1
\end{bmatrix}
\begin{bmatrix}
\v_1\\
\vdots\\
\v_n
\end{bmatrix} \
\end{equation}
which indicates that $ $ can be constructed via the product of the  $n×n$ upper triangular matrix $$ of elements $ u(k,j) $, and the $n ×(n+q)$ orthogonal matrix $ $ as
\begin{equation}
\Wmat = \Umat\Vmat \
\end{equation}
With this expression in hand, Eq. (\ref{eq:P_wrong_dimensions}) can now be written as
\begin{equation}\label{eq:covariance_for_comparison}
\Cov^-=	\Usub{\Dsub}\Usub\trans = \Wmat\Dhat\Wmat\trans = (\Umat\Vmat)\Dhat(\Umat\Vmat)\trans = \Umat[]\ [\Vmat\Dhat\Vmat\trans]\ \Umat\trans
\end{equation}
Since the product $ $ is constructed according to Eq. (\ref{eq:vDv_product}) and $$ is $n×(n+q)$, the bracketed term $[]$ from the previous equation is diagonal of size $n×n$. At this point, a term-by-term ``equivalence'' is possible for Eq. (\ref{eq:covariance_for_comparison}) which leads to the sought alternative factors $ $ and $ $:
\begin{equation}
\Cov^-=	\Usub{\Dsub}\Usub\trans = \Umat[][] \ [\Vmat\Dhat\Vmat\trans] \ \Umat\trans= \tilde{\Wmat}\tilde{\D}\tilde{\Wmat}\trans\ \end{equation}
with $ =$ and $ =$.
Therefore, the propagated $$ and $$ factors are given by,
\begin{equation}\label{eq:finalU}
\Usub = \Umat \
\end{equation}
\begin{equation}\label{eq:finalD}
\Dsub = \Vmat\Dhat\Vmat\trans \
\end{equation}

In summary, using the candidate factors $ $ and $ $ for the execution of the WMGS routine will provide $ $ and $ $ which are used in Eqs. (\ref{eq:finalU}) and (\ref{eq:finalD}) to obtain the propagated factors $ $ and $ $. Here, the WMGS procedure will be regarded and referenced as a function, such that  WMGS($,$) indicates its execution with arguments $$ and $$.
\section{Conventional UD Kalman measurement update}
The reformulation of the standard Kalman filter into the UD Kalman update requires several steps. First, the posterior covariance $P[][][+]$ is expressed in terms of the prior covariance $P[][][-]$. To accomplish this, the Kalman gain $K[][]=P[][][-](P[][][-]+)^-1$, where $_k =_k$, is substituted into the Kalman filter update covariance equation $	P[][][+]=(I-K[][])P[][][-]$, this is
\begin{align}\label{eq:expanded_cov_update}
\stdvec{P}[][][+]&=(\stdvec{I-\stdvec{K}[][]}\Hmat)\stdvec{P}[][][-]\nonumber\\
&=\Cov^--\stdvec{P}[][][-]\Hmat\trans(\Hmat\stdvec{P}[][][-]\Hmat\trans+\R)^{-1}\Hmat\Cov^- \
\end{align}
% \begin{equation}\label{eq:Pplusplus_after_putting_UDU_terms_in}
% 	\stdvec{P}[][][+]
% 	=\Cov^--\stdvec{P}[][][-]\Hmat\trans(\Hmat\stdvec{P}[][][-]\Hmat\trans+\R)^{-1}\Hmat\Cov^-\
% \end{equation}
Then, the covariance matrix $P[][][-]$ is replaced by its UD decomposition i.e., $P[][][-] = $. In addition, because in this formulation the measurements are processed sequentially, $$ will become a row vector $_i$, and $$ a scalar, $r_i$ (the $i^th$ diagonal element of $$). Integrating these changes into Eq. (\ref{eq:expanded_cov_update}) gives
% 	\begin{flalign}
% 	&\stdvec{P}[][][++] = \UDU - \UDU\Hmat_i\trans(\Hmat_i\UDU\Hmat_i\trans + r_i)^{-1}\Hmat_i\UDU+&&\\ \nonumber
% 	&\Dmat\UDU\Hmat_i\trans(\Hmat_i\UDU\Hmat_i\trans + r_i)^{-1}\Hmat_i\UDU\Dmat \ .
% 	\end{flalign}
\begin{equation}
\stdvec{P}[][][+] = \UDU - \UDU\Hmat_i\trans(\Hmat_i\UDU\Hmat_i\trans + r_i)^{-1}\Hmat_i\UDU
\end{equation}
If we define $ = _i$ and \hlc{the scalar} $ a_i=(_i_i+ r_i)^-1 $, the previous equation can be written as
\begin{equation}\label{eq:pureUD}
\stdvec{P}[][][+] = \UDU - \UD \ \wsub a_i \ \wsub\trans \ \DU
\end{equation}
and further factorizing the matrices $ $ and $ $ leads to
\begin{equation}\label{eq:UD_partial_update_before_udu_operation}
\stdvec{P}[][][+] = \Usub[\Dsub - \Dsub\wsub a_i\wsub\trans\Dsub]\Usub\trans
\end{equation}
From Eq. (\ref{eq:UD_partial_update_before_udu_operation}), the term in brackets is symmetric and positive definite and its UD decomposition can be computed. Let the UD decomposition of this term be
\begin{equation}\label{eq:UDU_mathcal}
\mathcal{U}\mathcal{D}\mathcal{U}\trans=\Dsub - \Dsub\wsub a_i\wsub\trans\Dsub
\end{equation}
such that equation (\ref{eq:UD_partial_update_before_udu_operation}) can now be written as
\begin{equation}
\stdvec{P}[][][+] = \Usub[\mathcal{U}\mathcal{D}\mathcal{U}\trans]\Usub\trans
\end{equation}
Because the product of two upper triangular matrices with ones in the diagonal is also an upper triangular matrix with ones in the diagonal, we have
\begin{equation}
\stdvec{P}[][][+] = [\Usub\mathcal{U}]\mathcal{D}[\mathcal{U}\trans\Usub\trans]=\Uplus\Dplus\Uplus\trans
\end{equation}
In other words, the factor $$ is updated via
\begin{equation}
\Uplus = \Usub\mathcal{U}
\end{equation}
while the update of the factor $$ is equal to $𝒟$, which results from the UD decomposition of the term $[- a_i]$ as per Eq.(\ref{eq:UDU_mathcal}):
\begin{equation}
\Dplus = \mathcal{D}
\end{equation}

Alternatively to this UD update approach where we perform the UD decomposition of the term $[- a_i]$, it is common to find in the literature that the UD Kalman update is done via the \textit{modified} Turner-Agee rank-one update (\hlc{or Carlson update}) \cite{agee1972triangular}. The modified Turner-Agee update is inspired by the same idea of updating the $$ and $$ factors directly to obtain a pair of new $$ and $$, and it also uses the terms present in Eq. (\ref{eq:UD_partial_update_before_udu_operation}). However, it is a more efficient and numerically stable approach. Nonetheless, both the UD decomposition approach shown here, and the execution of the modified Turner-Agee update, will provide the updated factors $$ and $$. The execution of the modified Turner-Agee can be done implementing the algorithm from Appendix \ref{ap:modified_Agee} given the terms $$, $$,$$, and $a_i$.

\textbf{Note: }The \textit{unmodified} or standard Agee-Turner also executes a rank-one update. However, it is numerically less stable if used when the updating term is negative \cite{Grewal2001}. In particular, the negative sign of the updating term $- a_i$ in equation \ref{eq:UD_partial_update_before_udu_operation} makes the use of the standard Agee-Turner update less suitable in this case. Otherwise, the standard Agee-Turner update would have been appropriate, e.g. with $+ a_i$. The implementation of the standard Turner-Agee can be done via \ref{ap:conventional_Agee}.

Table \ref{table:UDPU_filter} summarizes the UD Kalman filter. Note that the UD decomposition needed to initialize the $$ and $$ factors is indicated as udu(), and the modified Agee-Turner as modifiedAgeeTurner(). The UD decomposition mechanization is included in Appendix \ref{ap:udu_decomposition}.

Note that if $m$ measurements are available for processing, the UD Kalman update step needs to be executed $m$ times, one execution per measurement. It is important to mention that the UD update step from Table \ref{table:UDPU_filter} assumes uncorrelated measurements, i.e., a diagonal $$, and if correlated measurements are present, i.e., if $$ is not diagonal, the practitioner can execute a decorrelation procedure prior processing the measurements. The mechanization included in Appendix \ref{ap:decorrelation} can be used to decorrelate the available measurements (diagonalize $$) so that the UD update from Table \ref{table:UDPU_filter} can still be used.

\begin{table}[h!]
\caption{U-D Kalman filter update with sequential uncorrelated measurement processing}\label{table:UDPU_filter}
\centering
\setlength{\extrarowheight}{8pt}
\begin{tabular}{ |c|c| }
\hline
\textbf{Model} & $\begin{aligned} &\x_{k}=\f_{k-1}(\x_{k-1},\u_{k-1},\w_{k-1})\\ &\yktilde=\hkofxk+\vk\\
&\wk \sim \mathcal{N} (\zerovec,\Qmat_k)\\
&\vk \sim \mathcal{N}(\zerovec,\R_k)\end{aligned}$\\
\hline
\textbf{Initialize} & $\begin{aligned}
&\xpost_0=\x_0 \\
&[\Uplus_0,\Dplus_0]= \text{udu}(\stdvec{P}[][0][+])\ \text{[Appendix \ \ref{ap:udu_decomposition}]}\\
\end{aligned}$\\
\hline
\textbf{UD Propagation} & $\begin{aligned} &\xprior_{k}=\f_{k-1}(\xpost_{k-1},\u_{k-1})\\
&\Wmat_{k-1} = \begin{bmatrix}
\Fmat_{k-1}\Uplus_{k-1} && \Gmat_{k-1}
\end{bmatrix}\ ;
\Dhat_{k-1} = \begin{bmatrix}\Dplus_{k-1} && \zerovec\\
\zerovec && \Qmat_{k-1}
\end{bmatrix}\\
&[\Usub_{k},\Dsub_{k}]=\text{WMGS}(\Wmat_{k-1},\Dhat_{k-1}) \ \text{[Section \ref{subsec:overview_UD}}] \\
%		&\Dsub= \Vmat\Dhat\Vmat\trans
\end{aligned}$\\
\hline
\textbf{UD update} &
%   \text{Given} & \ \Usub, \Dsub, \wsub \ \text{and} \ a_i\text{, do:}\\
%a_i &=(\wsub\trans\Dsub\wsub+r_i)^{-1};\ \wsub=\Usub\trans\Hmat_i\trans\\
[\Uplus_k,\Dplus_k,\Kmat_k] &=\text{ modifiedAgeeTurner}(\Usub_k, \Dsub_k, r_{i_k},  \Hmat_{i_k})\ \text{[Appendix \ \ref{ap:modified_Agee}]} \\
\xpost_k&=\xprior_{k}+\stdvec{K}[][k](\ym_k-\h(\xhat_k^-))\\
\end{aligned}$\\
\hline
\end{tabular}
\end{table}

\section*{Appendix}

\subsection{UD decomposition}\label{ap:udu_decomposition}
Given a symmetric, positive-definite, $m×m$ matrix $M$, its UD decomposition can be obtained via the following mechanization.
\begin{flalign}
\textbf{for}&\ j=m:-1:1&&\\
&\textbf{for}\ i=j:-1:1&&\\
& \quad \sigma=M(i,j)&&\\
& \quad \textbf{for}\ k=j+1:m&&\\
&\qquad \sigma=\sigma-U(i,k)D(k,k)U(j,k)&&\\
& \quad \textbf{end for}&&\\
&\quad \textbf{if}(i==j)&&\\
&\qquad D(j,j)=\sigma&&\\
&\qquad U(j,j)=1&&\\
&\qquad U(i,j)=\sigma/D(j,j)&&\\
&\quad\textbf{end if}&&\\
&\ \ \textbf{end for}&&\\
&\textbf{end for}\\
\end{flalign}

The $$ and $$ matrices constructed via this mechanization satisfy $M=$. Note that in Table \ref{table:UDPU_filter}, this algorithm is referred as the \textit{function} \text{udu()}.

\subsection{Modified Agee-Turner rank-one update}\label{ap:modified_Agee}
The UD Kalman filter uses the following algorithm to update the UD factors. Providing $$, $$, $r_i$, $_i$, the algorithm returns $$, the elements of $_j$ as [$d_1  d_2  d_n$], and the Kalman gain $$. This approach is taken from \cite{Gibbs2011a}.
\begin{flalign}
&\w = \Usub\trans\Hmat_i \ \text{where} \ \w = [w_1\ w_2 \ \hdots \ w_n]\trans&&\\
&\vchar = \Dsub\trans\w \ \text{where} \ \vchar=[\barbelow{d}_1 w_1\ \barbelow{d}_2 w_2 \ \hdots \ \barbelow{d}_n w_n]\trans \\
&\Kmat_1 = [v_1\ 0\ \hdots \ 0]\trans \\
&\alpha_1 = r_i+v_1w_1 \\
&d_1 = (r_i\barbelow{d}_1)/\alpha_1\\
&\textbf{for}\ j=2 \hdots \ n &&\\
&\ \ \ \ \alpha_j = \alpha_{j-1}+v_jw_j\\
&\ \ \ \ d_j = d_j\alpha_{j-1}/\alpha_j\\
&\ \ \ \ \lambda_j = -w_j/\alpha_{j-1}\\
&\ \ \ \ \Uplus_j = \Usub_j + \lambda_j\Kmat_{j-1}\\
&\ \ \ \ \Kmat_j = \Kmat_{j-1}+v_j\Usub_j\\
&\textbf{end for}\\
&\text{Finally compute}\nonumber\\
&\Kmat = \Kmat_{n}/\alpha&&
\end{flalign}
In this algorithm, $_j$, $_j$, and $_j$, refer to the $j^th$ column of the corresponding matrix. Note that in Table \ref{table:UDPU_filter}, this algorithm is referred as the \textit{function} \text{modifiedAgeeTurner()}.

\subsection{The Agee-Turner rank-one update}\label{ap:conventional_Agee}
This algorithm provides a mechanization to update the $n×n$ factors $$ and $$ with the information contained in $caa$ according to
\begin{equation}
\overset{+}{\stdvec{U}}\overset{+}{\stdvec{D}}\overset{+}{\stdvec{U}}\trans=\Umat\D\Umat\trans + c\stdvec{a}\stdvec{a}\trans .
\end{equation}
Where c is a positive scalar and $a$ is a column vector of size $n$. This algorithm returns the updated factors $+U $ and $+D$ given $, ,c$, and $a$. The algorithm is executed as follows:

\begin{flalign}
&C_n = c&&\\
&\textbf{for}\ j=n,n-1,\hdots,2\\
&\ \ \ \ \Uplus_{jj}=1\\
&\ \ \ \ \Dplus_{jj} = \D_{jj}+C_ja_j^2\\
&\ \ \ \ \textbf{for}\ k=1, \hdots, j-1\\
&\ \ \ \ \ \ \ a_k \coloneqq a_k-a_j\Umat_{k,j}\\
&\ \ \ \ \ \ \ \Uplus_{kj}=\Umat_{kj}+c_j a_j a_k/\Dplus_{jj}\\
&\ \ \ \ \ \textbf{end for}\\
&\ \ \ C_{j-1} = C_j\D_{jj}/\Dplus_{jj}\\
&\textbf{end for}\\
\end{flalign}
Note that all of the sub-indices are used to access the matrix elements.

\subsection{Measurement decorrelation using UD factors}\label{ap:decorrelation}
Let the UD decomposition of a non-diagonal measurement noise covariance, $ _c $, be
\begin{equation}
\R_c = \Ur\Dr\Ur\trans \
\end{equation}
and consider the measurement equation $= +$̌, being transformed by $\Uinv$ as

\begin{equation}
\stdvec{z} = \Uinv\y = \Uinv\Hmat\x+\Uinv\v \
\end{equation}

The residual (measurement minus expected measurement) can then be calculated as

\begin{equation}
\stdvec{e}=\stdvec{z}-\hat{\stdvec{z}}=\Uinv\Hmat(\stdvec{x}-\hat{\stdvec{x}})+\Uinv\v \
\end{equation}

and its corresponding covariance as

\begin{equation}
\expect[\stdvec{e}\stdvec{e}\trans]=\expect[(\Uinv\Hmat(\stdvec{x}-\hat{\stdvec{x}})+\Uinv\v)(\Uinv\Hmat(\stdvec{x}-\hat{\stdvec{x}})+\Uinv\v)\trans ] \
\end{equation}

\begin{equation}
\expect[\stdvec{e}\stdvec{e}\trans]=\Uinv\Hmat\expect[(\stdvec{x}-\hat{\stdvec{x}})(\stdvec{x}-\hat{\stdvec{x}})]\trans\Hmat\trans\Ur^{-T} +\Uinv\expect[\v\v]\trans\Ur^{-T} \
\end{equation}

\begin{equation}
\expect[\stdvec{e}\stdvec{e}\trans]=\Uinv\Hmat\Cov\trans\Hmat\trans\Ur^{-T} +\Uinv\R_c\Ur^{-T} \
\end{equation}

assuming zero mean Gaussian white noise.
But since

\begin{equation}
\Uinv\R_c  \Ur^{-T}=\Dr \
\end{equation}

and by naming $\Hmat_z=\Uinv\Hmat$, this results in

\begin{equation}
\expect[\stdvec{e}\stdvec{e}\trans] = \Hmat_z\Cov\Hmat_z\trans+\Dr \
\end{equation}

That is, the transformed equation now uses a diagonal measurement noise covariance. In order to properly use the transformed, now uncorrelated measurements, the filter simply uses $ \Uinv\Hmat $ instead of $ \Hmat $, and instead of directly computing the measurement residual with $ (\y-\Hmat\x)$ or $ (\y-\h(\xprior))$ , we now use $\Uinv(\y-\Hmat\x)$. Rather than using the correlated measurement noise $ \R_c $, the filter now will use $ \Dr $.

§ FUNDING SOURCES

This work was supported under Air Force contract FA8651-20-F-1025.