Roman Klokov
email: roman.klokov@inria.fr

Jakob Verbeek
email: jakob.verbeek@inria.fr

Edmond Boyer
email: edmond.boyer@inria.fr

Probabilistic Reconstruction Networks for 3D Shape Inference from a Single Image

We study end-to-end learning strategies for 3D shape inference from images, in particular from a single image. Several approaches in this direction have been investigated that explore different shape representations and suitable learning architectures. We focus instead on the underlying probabilistic mechanisms involved and contribute a more principled probabilistic inference-based reconstruction framework, which we coin Probabilistic Reconstruction Networks. This framework expresses image conditioned 3D shape inference through a family of latent variable models, and naturally decouples the choice of shape representations from the inference itself. Moreover, it suggests different options for the image conditioning and allows training in two regimes, using either Monte Carlo or variational approximation of the marginal likelihood. Using our Probabilistic Reconstruction Networks we obtain single image 3D reconstruction results that set a new state of the art on the ShapeNet dataset in terms of the intersection over union and earth mover's distance evaluation metrics. Interestingly, we obtain these results using a basic voxel grid representation, improving over recent work based on finer point cloud or mesh based representations.

Introduction

The overwhelming success of convolutional neural networks on image data [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF][START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] instigated the exploration of CNNs for other problems, in particular in 3D visual computing. 3D CNNs for shapes represented with uniform voxel grids have been investigated for recognition [START_REF] Maturana | VoxNet: A 3D convolutional neural network for real-time object recognition[END_REF][START_REF] Wu | 3D ShapeNets: A deep representation for volumetric shapes[END_REF] and generative modelling tasks [START_REF] Brock | Generative and discriminative voxel modeling with convolutional neural networks[END_REF][START_REF] Wu | Learning a probabilistic latent space of object shapes via 3D generative-adversarial modeling[END_REF]. For 3D shape inference, initial works [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF][START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF] successfully demonstrated the ability of 3D CNNs to produce coherent voxelized shapes given single images. This task has since gained a significant attention, as a result of its vast application field and despite its challenging ill-posed nature.

Further exploring CNNs in this context, recent works have investigated beyond straightforward adaptions of 2D CNNs to 3D voxel grids, notably to overcome the cubic complexity in time and memory associated with it. For instance, sparse representations of large voxel grid have been proposed to reduce complexity while allowing for finer shape details [START_REF] Graham | 3D semantic segmentation with submanifold sparse convolutional networks[END_REF][START_REF] Tatarchenko | Octree generating networks: Efficient convolutional architectures for high-resolution 3D outputs[END_REF]. Other more scalable shape representations suitable for recognition and generation tasks have been investigated, including rendered images [START_REF] Su | Multi-view convolutional neural networks for 3D shape recognition[END_REF], geometry images [START_REF] Sinha | Deep learning 3D shape surfaces using geometry images[END_REF], 2D depth maps [START_REF] Soltani | Synthesizing 3D shapes via modeling multi-view depth maps and silhouettes with deep generative networks[END_REF], point clouds [START_REF] Klokov | Escape from cells: Deep Kd-networks for the recognition of 3D point cloud models[END_REF][START_REF] Qi | Pointnet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Qi | Pointnet++: Deep hierarchical feature learning on point sets in a metric space[END_REF], and graphs [START_REF] Monti | Geometric deep learning on graphs and manifolds using mixture model CNNs[END_REF][START_REF] Verma | Feastnet: Feature-steered graph convolutions for 3D shape analysis[END_REF]. Importantly, these representations Conditioning between 2D and 3D tensors is achieved by means of FiLM [START_REF] Perez | FiLM: Visual reasoning with a general conditioning layer[END_REF] layers. The inference network q ψ is only used during training for variational inference.

come with specific network architectures and loss functions suited to the corresponding data structures.

The variety of approaches complicates comparisons and identification of the sources of improved performance. In particular, most works do not decouple the problems inherent to the task, and mix new shape representations, along with the associated network architectures and loss functions, with different image conditioning schemes, different probabilistic formulations of the shape prediction task, and in some cases the use of additional training data. This leads to possibly difficult or even unfair comparisons, due to the inability to confidently determine the source of improvements in new models, and emphasizes the need for a more systematic approach.

In this paper we look at the single image 3D shape inference through the prism of a family of generic probabilistic latent variable models, which we term Probabilistic Reconstruction Networks (PRN), see Figure 1. The formalism encompassing these models naturally decouples different aspects of the problem including the shape representation, the image conditioning, and the usage of latent shape space for the shape prediction. It also allows to categorize previous models for this task by their structural properties. Without loss of generality, we use voxel grids as shape representations and focus our attention on other aspects. We systematically analyze the impact of several design choices: (i) the dependency structure between the input image, the latent shape variable and the output variables; (ii) the effectiveness of training the model using Monte Carlo sampling or variational inference to approximate the log-likelihood; (iii) the effectiveness of a deterministic version of the model that suppresses any uncertainty associated with the latent variable; (iv) the effect of jointly learning the shape reconstruction model along with a generative shape model, which share their latent shape space.

For our experiments we use the ShapeNet dataset for the single image 3D shape reconstruction. We obtain excellent single image 3D reconstruction results with our Probabilistic Reconstruction Networks, setting new state-of-the-art results in terms of the IoU and EMD performance metrics. Interestingly, our results improve over recent works based on pointcloud and mesh-based shape representations.

In summary, our contributions are:

• a generic latent variable model for the single image 3D shape reconstruction;

• exploration of modeling options in a systematic and comparable manner;

• new state-of-the-art single image reconstruction results on the ShapeNet dataset. In the following we first introduce our generic latent variable model in Section 2, which is then used to review and categorize previous work on shape inference from a single image in Section 3. We then present our experimental results and comparisons in Section 4.

Probabilistic framework for 3D shape reconstruction

Below we present our generic latent variable model in Section 2.1, and detail the network architectures used for our experiments in Section 2.2. Although proposed probabilistic model is agnostic to the underlying shape representation, we present it using voxel grid representation in accordance to our experimental setup.

Latent variable model for single image 3D shape reconstruction

We consider a shape v as a uniform voxel occupancy grid of a predefined resolution. Our task is to predict the shape v given an input image i, i.e. to model p (v|i). While images and occupancy grids live in different spaces, both are representations of an underlying object that has a 3D shape and an appearance. Using this observation, we assume that an observed shape v has a latent parametrization z within a latent shape space of lower dimension, that captures shape variations. We then define our latent variable model for single image 3D shape reconstruction as:

p (v|i) = z p θ (v|z, i) p φ (z|i) dz, (1)
where φ and θ are parameters of the model. This model consists of two modules: an image informed latent variable prior p φ (z|i), and a decoder p θ (v|z, i) that predicts the shape v given the image and the latent variables. Being a generic latent variable model, it allows us to decouple and study different aspects of the single image 3D shape reconstruction task. Image conditioning options. In Eq. (1) both latent variable prior and decoder are conditioned on the input image. If we drop the dependence on i from one module, we maintain dependence of the shape v on the image i, and obtain two alternative models:

p (v|i) = z p θ (v|z) p φ (z|i) dz, (2)
p (v|i) = z p θ (v|z, i) p φ (z) dz. (3)
In the first case, we omit the image conditioning in the decoder p θ (v|z) and assume that the image conditioned prior p φ (z|i) is sufficient to obtain a valid reconstruction by the decoder. This dependency structure corresponds to an assumption of conditional independence of i and v given z. In the second case, we leave the image dependence in the decoder p θ (v|z, i) but use an unconditional prior for the latent variable p φ (z). This corresponds to an assumption that i and z are a-priori independent. If we drop the image conditioning in both components, the model becomes a generative latent variable shape model:

p (v) = p θ (v|z) p κ (z) dz.
Latent variable sampling during training. Due to the non-linear dependencies in the integral in the models defined in equations ( 1)-( 3), exact computation of the log-likelihood and its gradient is intractable. We consider two alternative approaches to overcome this difficulty. The first is to use a Monte Carlo approximation, as e.g. in [START_REF] Gordon | Meta-learning probabilistic inference for prediction[END_REF],

ln p(v|i) ≈ ln 1 M M ∑ m=1 p θ (v|z m ), z m ∼ p φ (z|i), (4)
where we make use of the re-parametrization trick [START_REF] Kingma | Auto-encoding variational Bayes[END_REF][START_REF] Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF] to back-propagate the gradient of the log-likelihood w.r.t. φ through the sampling from p φ (z|i).

Alternatively, we can use the variational inference framework [START_REF] Kingma | Auto-encoding variational Bayes[END_REF][START_REF] Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF] to obtain a training signal based on more informed samples from the latent variable. We introduce a variational approximate posterior q ψ (z|v, i), which is learned jointly with the prior and decoder by the maximization of the variational lower bound on the log-likelihood:

L (φ , ψ, θ , v, i) = IE q ψ (z|v,i) [log p θ (v|z, i)] -D KL q ψ (z|v, i) ||p φ (z|i) ≤ ln p(v|i). (5)
To evaluate this bound and its gradient, we sample from q ψ (z|v, i). Since the samples are conditioned on the shape, unlike the Monte Carlo approximation case, we expect this approach to be more sample efficient. On the one hand, image conditioning in the posterior may be omitted, since the posterior is already conditioned on the shape information. On the other hand, conditioning on the image may in principle further improve the accuracy of the approximate posterior, as it is based on more information. Deterministic shape model. We can obtain deterministic versions of the presented models by considering z as a function of φ and, if required, i. Although this simplifies the models, it also discards an important property. Typically, p(v|z) is factorized, with each voxel occupancy being modelled as an independent Bernoulli, e.g. [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF][START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF]. With this factorization, any shape ambiguity given a single image cannot be modelled properly, since self-occluded parts of the shape lack structure in the prediction. A latent variable that conditions the factorized distribution can be used to induce dependencies among the voxel occupancies to reflect the structured ambiguity resulting from partially observed shapes.

In case of the variational training we also use a deterministic posterior, and substitute the KL-divergence in Eq. ( 5) with a suitable similarity measure. For example, the L 2 -norm of the difference between the output of image encoder p φ and the posterior q ψ . Merging unconditional generation with reconstruction. The models presented above can be trained along with a generative shape model, that may be of interest on its own, or used to regularize the conditional model. To achieve this, we consider a variational bound on the log-likelihood of an unconditional generative model:

L (κ, ψ, θ , v) = IE q ψ (z|v) [log p θ (v|z)] -D KL q ψ (z|v) ||p κ (z) ≤ ln p(v), (6)
where p(v) = p θ (v|z)p κ (z)dz. Looking at Eq. ( 5) and Eq. ( 6), we observe that if we omit image conditioning from the decoder and the posterior in Eq. ( 5), and share their parameters in both conditional and unconditional models, we can obtain a unified training objective:

L (κ, φ , ψ, θ , v, i) = IE q ψ (z|v) [log p θ (v|z)] - 1 2 D KL q ψ (z|v) ||p κ (z) - 1 2 D KL q ψ (z|v) ||p φ (z|i) , (7)
which is the average of the lower bounds on the marginal likelihood on the shape v with and without conditioning on the image. The term corresponding to the unconditional likelihood can be viewed as a regularization that encourages z from the entire latent space to correspond to realistic shapes, instead of just the z in the support of the conditional distributions p(z|i) on latent coordinates given an input image. Note that in the generative model neither the decoder p θ (v|z) nor the encoder q ψ (z|v) are conditioned on the image. These can thus be shared with the single image reconstruction model if the latter is not conditioned on the image in these components.

Network architectures

Although our probabilistic model is not representation specific, we focus on the voxel grid shape representation, and leave the comparison to alternative representations for the future work. Thus, we implement the different conditional distributions in our model as 2D and 3D CNNs that output the parameters of distributions on the latent variable z, or on the voxel occupancies. In particular,

• The unconditional Gaussian latent prior p κ (z) is characterized by κ that consists of means and variances for all latent dimensions implemented as trainable parameters. • The image conditioned prior p φ (z|i) is a 2D CNN, consisting of six blocks of pairs of convolutions: a standard and a strided one, interleaved with batch normalization and point-wise non-linearities, followed by two fully-connected layers. It processes input images into the means and variances of a factored Gaussian on z. • The shape conditioned variational posterior q ψ (z|v, i) is a 3D CNN consisting of an initial 3D convolution and a series of four modified residual blocks [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] He | Identity mappings in deep residual networks[END_REF], each using an additional 1×1 convolution instead of identity and feature map concatenation instead of summation and each followed by 2×2×2 spatial average pooling. Final convolutional features are fed to two additional fully-connected layers. This encoder processes input shapes into the means and variances of a factored Gaussian on z. • The 3D deconvolutional decoder p θ (v|z, i) is mirrored from the approximate posterior q ψ , with the pooling being substituted by the 2×2×2 upscaling by trilinear interpolations, producing the parameters of Bernoulli distributions on the voxel occupancies from a latent variable input. Image conditioning in the two latter modules is inspired by the FiLM conditioning mechanism [START_REF] Perez | FiLM: Visual reasoning with a general conditioning layer[END_REF]. Intermediate 2D feature maps from the first five convolutional blocks of the image encoder p φ (z|i) are averaged spatially, transformed by two additional fully-connected layers into weights and biases, that are used to scale the according five intermediate batchnormalized 3D feature maps in the shape encoder, the latent variable decoder, or both. Instead of affine transformation used in FiLM, we use non-negative scaling weights by predicting them in logarithmic scale, in our experiments this resulted in more stable training and slightly better results. See Figure 1 for a schematic overview of the model architecture.

To ensure fair comparison between the variations of the model we use identical architectures for every component concurring in different models, except for additional fullyconnected layers associated with the different image conditioning options. When we include an unconditional generative shape model and optimize Eq. ( 7), we share the decoder p θ (v|z) and the variational posterior q ψ (z|v) between the conditional and unconditional models. Exact architectures, training procedures and their hyperparameters are available on the implementation code page. q (z|s)

[5] p (s|z) p (z|i) q (z|s) [START_REF] Mandikal | 3D-LMNet: Latent embedding matching for accurate and diverse 3D point cloud reconstruction from a single image[END_REF] Table 1: Overview of how related work fits into our probabilistic reconstruction framework.

Related work

In this section we review related work on single image 3D shape reconstruction, and relate it to our generic latent variable model presented in the previous section. Representations for 3D shape inference. The majority of works studying the inferencebased single image 3D shape reconstruction introduce new shape representations and suitable neural network architectures. The seminal works by Choy et al. [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF] and Girdhar et al. [START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF] used 3D CNNs to predict voxel occupancy grids. To reduce the computational complexity of the voxel grid representation Tatarchenko et al. [START_REF] Tatarchenko | Octree generating networks: Efficient convolutional architectures for high-resolution 3D outputs[END_REF] proposed an architecture to process octrees computed on top of the voxel grids. Richter and Roth [START_REF] Richter | Matryoshka Networks: Predicting 3d geometry via nested shape layers[END_REF] proposed to use a set of six depth maps to represent voxel grids and to combine a series of such sets in a nested manner to model non-trivial shapes. Su et al. [START_REF] Su | A point set generation network for 3D object reconstruction from a single image[END_REF] combined 2D CNNs with fully-connected networks to output point clouds, which are learned by optimizing Chamfer distance or differentiable approximation of earth mover's distance. Wang et al. [START_REF] Wang | Pixel2Mesh: Generating 3D mesh models from single RGB images[END_REF] applied graph-convolutional networks [START_REF] Bronstein | Geometric deep learning: going beyond Euclidean data[END_REF] to the mesh-based shape representation. Shin et al. [START_REF] Shin | Pixels, voxels, and views: A study of shape representations for single view 3d object shape prediction[END_REF] proposed to predict multiple depth maps and according silhouettes and fuse them into meshes by post-processing with Poisson reconstruction algorithm.

Although related in their overall goals, these approaches are difficult to compare since they use target shapes approximated to different degrees. Ideally, a fair comparison across shape representations should be performed while maintaining the same level of granularity across representations, and for different levels, since it is possible that some representations work better for rough shape reconstruction, while other are best for detailed reconstruction. Image-shape consistency and additional data. Another significant stream of works originates from the idea of ensuring consistency between input data and target shapes. Initial work by Yan et al. [START_REF] Yan | Perspective Transformer Nets: Learning single-view 3D object reconstruction without 3D supervision[END_REF] introduced the consistency between 3D shapes and their silhouettes produced by different viewpoints in a form of a loss function. Similar ideas were investigated by Wiles and Zisserman [START_REF] Wiles | SilNet: Single-and multi-view reconstruction by learning from silhouettes[END_REF]. Tulsiani et al. [START_REF] Tulsiani | Multi-view supervision for single-view reconstruction via differentiable ray consistency[END_REF] expanded this approach by the use of differentiable ray tracing in the loss function, ensuring correspondence of inferred voxelized shapes to foreground segmentation masks and depth images. Wu et al. [START_REF] Wu | MarrNet: 3D shape reconstruction via 2.5D sketches[END_REF] developed the idea even further and introduced a two-step reconstruction framework. The first part of the model is trained to infer 2.5D shape sketches (unions of segmentation, depth and normal maps) from images, while the second is separately trained to predict shapes from 2.5D sketches. Both components are then fine-tuned, using reprojection consistency.

Henderson and Ferrari [START_REF] Henderson | Learning to generate and reconstruct 3d meshes with only 2d supervision[END_REF] proposed a probabilistic framework for image generation conditioned on a latent shape variable and an additional latent variable for the shape pose. This framework was used to train an underlying 3D mesh generator with the help of differentiable rendering of 3D meshes into images. Differentiable point clouds [START_REF] Insafutdinov | Unsupervised learning of shape and pose with differentiable point clouds[END_REF] Relations to our framework. In addition to the work discussed above, given similarities between VAEs and GANs [START_REF] Hu | On unifying deep generative models[END_REF], our work is also related to adversarial approaches involving shape discriminators [START_REF] Smith | Improved adversarial systems for 3d object generation and reconstruction[END_REF][START_REF] Wu | Learning a probabilistic latent space of object shapes via 3D generative-adversarial modeling[END_REF][START_REF] Wu | Learning shape priors for single-view 3D completion and reconstruction[END_REF]. In Table 1 we organize related work in terms of how it fits into our generic probabilistic reconstruction framework, abstracting away from the implementation of various components.

We see that most previous works use a dependency structure where the latent variable is inferred from the image, and the shape decoder only depends on the latent variable and not on the image. Moreover, most works rely on deterministic models, except for GANbased approaches of [START_REF] Smith | Improved adversarial systems for 3d object generation and reconstruction[END_REF][START_REF] Wu | Learning a probabilistic latent space of object shapes via 3D generative-adversarial modeling[END_REF], the point cloud based approaches of [START_REF] Mandikal | 3D-LMNet: Latent embedding matching for accurate and diverse 3D point cloud reconstruction from a single image[END_REF][START_REF] Su | A point set generation network for 3D object reconstruction from a single image[END_REF], and the mesh based method of [START_REF] Henderson | Learning to generate and reconstruct 3d meshes with only 2d supervision[END_REF]. Finally, only TL-Networks [START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF] and 3D-LMNet [START_REF] Mandikal | 3D-LMNet: Latent embedding matching for accurate and diverse 3D point cloud reconstruction from a single image[END_REF] make use of the variational inference for shape modelling.

Experiments

In this section we present the experimental setup, our quantitative and qualitative evaluation results, as well as our analysis of these results.

Dataset, evaluation metrics, and experimental details

Dataset. We evaluate PRNs on a subset of the ShapeNet dataset [START_REF] Chang | ShapeNet: An information-rich 3D model repository[END_REF] introduced by Choy et al. [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF]. It contains about 44k 3D shapes from 13 major categories of ShapeNet dataset represented as voxel grids of resolution 32 3 , as well as renderings from 24 different randomized viewpoints as 137 2 images. Following Choy et al., we use 80% of the shapes from each category for training and remaining shapes for testing. Evaluation metrics. We evaluate using the standard intersection-over-union (IoU) metric [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF], which averages the per-category IoU metric between the inferred shape and the ground-truth voxel representation. In addition, to allow comparison to recent work based on point-cloud and mesh based representations, we also report the Chamfer distance (CD) and earth mover's distance (EMD), computed using the code of [START_REF] Sun | Pix3D: Dataset and methods for single-image 3D shape modeling[END_REF], where we removed the square root from the distance computations in CD to make it comparable to the related work. In particular, each ground truth and predicted voxel grid is mapped to a point cloud by sampling the surface induced using the marching cubes algorithm [START_REF] Lewiner | Efficient implementation of marching cubes' cases with topological guarantees[END_REF]. We then compute the CD and EMD on the resulting point clouds.

Training and evaluating. When using either Monte Carlo approximation or a variational objective function, we always use a single sample to compute the gradients during training. A unified training protocol is used for all the models: all components are trained simultaneously (contrary to [START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF][START_REF] Mandikal | 3D-LMNet: Latent embedding matching for accurate and diverse 3D point cloud reconstruction from a single image[END_REF][START_REF] Wu | Learning shape priors for single-view 3D completion and reconstruction[END_REF]), with the AMSGrad [START_REF] Reddi | On the convergence of Adam and beyond[END_REF] optimizer with decoupled weight decay regularization [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF] with step-like scheduling for learning rate and weight decay parameter, and restarts of gradient moments accumulation at the beginning of each step.

During testing, we use a deterministic approach. In particular, we take the means of the conditional distributions rather than samples from them. We found this to significantly improve the reconstruction quality, compared to sampling.

Experimental results

Evaluation of PRN variants. To explore the various possibilities of our general latent variable model, we consider three options to condition on the image: (i) using the latent space to carry all image information: p(v|z)p(z|i), (ii) using additional conditioning of the decoder on the image: p(v|z, i)p(z|i), and (iii) using an uninformed prior on the latent variable: p(v|z, i)p(z). To train the models we either approximate the integral in Eq. ( 1) directly with Monte Carlo samples from the prior on z, or with a variational lower bound and samples from the variational posterior. We also test a deterministic model, where the distribution on the latent variable is replaced by a deterministic function. Finally, we consider the option to train the model jointly with an an unconditional generative model. In Table 2 we present the results, using two thresholds τ on the voxel occupancy probability: the neutral 0.5, and following [START_REF] Choy | 3D-R2N2: A unified approach for single and multi-view 3d object reconstruction[END_REF] the looser 0.4 which overall leads to improved IoU scores.

In case of the Monte Carlo approximation (top three rows), additional image conditioning in the decoder improves the results. Conditioning both the latent variable prior and the decoder on the image achieves best results, suggesting that these different pathways to use the image are complementary.

The use of variational training consistently improves the results over the Monte Carlo approximation. In this case, the additional image conditioning of the decoder or the variational posterior, see line five and six, is not effective and even somewhat reduces the performance. This is contrary to the results obtained with Monte Carlo sampling; in the latter case the sampling inefficiency is probably partially compensated by the additional conditioning pathway. Variational inference leads to more accurate samples, which obviates the need for the additional image conditioning (at least for the chosen mechanism of the additional image conditioning).

We also consider a deterministic variant of our best performing model, which resembles the TL-network of [START_REF] Girdhar | Learning a predictable and generative vector representation for objects[END_REF]. The results show that probabilistic handling of the latent variable reduces overfitting in the model and leads to IoU of 2.5 points higher. Finally, we also tested the training with a joint generative shape model, which TL-Networks used as pretraining. Although we did not observe a significant effect due to the joint training with a generative shape model, it does offer additional functionality by being able to sample shapes, Table 3: Comparison of PRN to the state-of-the-art. All results are taken from the original papers, except for †, which were provided in [START_REF] Su | A point set generation network for 3D object reconstruction from a single image[END_REF]. Pixel2Mesh additionally uses camera information and surface normals during training.

or compute their likelihoods under the model.

Comparison to the state-of-the-art. In Table 3 we compare to earlier state-of-the-art approaches. All methods use the same input images, but use slightly different image preprocessing: 3D-R2N2 uses random cropping, 3D-LMNet central cropping, AtlasNet crops and resizes, while PSGN and Pixel2Mesh resize the image. As OGN, we use original images, but also add a grey-scale version of each image as a fourth input channel.

In Table 3 we report Chamfer distance computed for 1024 predicted and 1024 ground truth points as it was done in most of the related work. Although they all used different numbers of predicted and ground truth points for training, The authors of [START_REF] Groueix | A papier-mâché approach to learning 3D surface generation[END_REF][START_REF] Mandikal | 3D-LMNet: Latent embedding matching for accurate and diverse 3D point cloud reconstruction from a single image[END_REF][START_REF] Wang | Pixel2Mesh: Generating 3D mesh models from single RGB images[END_REF] explicitly state this protocol for evaluation, while there is no information about it in the text of [START_REF] Su | A point set generation network for 3D object reconstruction from a single image[END_REF]. Judging from the code of [START_REF] Su | A point set generation network for 3D object reconstruction from a single image[END_REF], we assume that, in their case, the metric was obtained for 1024 predicted and 16384 ground truth points. For the reference, we recomputed CD under the same protocol and obtained a better value: 3.90. This shows that the metric is affected by a negative bias, which decreases with the increasing number of evaluated points, and underlines the need for unified evaluation protocol.

Our PRN obtains excellent results, and significantly improves over previous state-of-theart results in terms of IoU and EMD, including methods based on point cloud and mesh representations. Point-based approaches use loss functions based on the Chamfer distance, and so naturally perform well in terms of this metric, but this does not per se transfer to better performance in the other metrics. In our case, we do not explicitly optimize for either of these metrics, relying on the binary cross-entropy for the voxel occupancies instead, and yet obtain best results in two of the three metrics (with only one competitor being better in terms of the third metric).

Qualitative reconstruction results. In Figure 2 we provide a selection of qualitative reconstruction results. We show results for the models in rows one, four and seven in Table 2, i.e. with the p(v|z)p(z|i) dependency structure, using Monte Carlo (PRN MC) and variational (PRN var.) training, and the deterministic version of the latter (PRN var. det.). We show four examples where the variationally trained model is the best, and one case where it is the worst. Overall, variationally trained model output fewer failed reconstructions, as well as more detailed reconstructions, compared to more failures and over simplified reconstructions from the model trained with Monte Carlo. For reference, the average IoU score of the variationally trained model is 66.2 (median 69.9), which corresponds to a fairly accurate reconstruction level, in particular given the challenging nature of the task.

Conclusion

In this paper we presented Probabilistic Reconstruction Networks, a generic probabilistic framework for 3D shape inference from single image. This framework naturally decouples different aspects of the problem, including the shape representation, the image conditioning structure, and the usage of the latent shape space. In our experiments with voxel-grid shape representations, we systematically explored the impact of image conditioning, Monte Carlo vs. variational likelihood approximation for training, the stochastic nature of the latent variable, and joint training with a generative shape model. We obtained single image shape reconstruction results that surpass the previous state of the art in terms of the IoU and EMD performance metrics, and outperform recent work based on point-cloud and mesh-based shape representation.

Given the interpretation of the inference-based reconstruction as an instance of conditional generation, future work includes further adaptation of the generative modelling approaches to the task, as well as the investigation of different shape representations within the proposed framework.

A Histograms of IoU for selected methods

In Figure 3 we provide a histogram of the IoU scores obtained across the shapes in the ShapeNet test set using three of the model evaluated in Table 2: • Using p(v|z)p(z|i), with Monte Carlo training (Table 2, row 1).

• Using p(v|z)p(z|i), with variational training (Table 2, row 4).

• Using p(v|z)p(z|i), with deterministic modeling (Table 2, row 7). For each shape in the test set there are 24 views, giving a total of about 210k shape inferences. The histograms show that the variational learning approach leads to more accurate reconstructions, leading to the largest number of reconstructed shapes in the last three bins for shape with > 70% IoU. For all other bins of less accurate results, the variational method has the smallest number of shapes.

Compared to the deterministic model, Monte Carlo training leads to more accurate reconstructions, but also to more very poor reconstructions.

B Visualization of shape reconstruction results

In this section we provide visualization of additional shape reconstruction results, similar to the ones presented in Section 4. Contrary to them, we put randomly sampled examples here.

For each example in Figure 4, we show from left to right:

• the input image;

• ground-truth shape;

• inferred shape with Monte Carlo training;

• inferred shape with variational training;

• inferred shape with deterministic model. These shape inference approaches correspond to rows one, four, and seven of Table 2.

Figure 1 :

 1 Figure 1: Probabilistic Reconstruction Networks for 3D shape inference from a single image. Arrows show the computational flow through the model, dotted arrows show optional image conditioning. Conditioning between 2D and 3D tensors is achieved by means of FiLM [23] layers. The inference network q ψ is only used during training for variational inference.

Figure 2 :

 2 Figure 2: Qualitative reconstruction results for three variants of PRNs.

Figure 3 :

 3 Figure 3: Histogram of IoU values on the ShapeNet test set for the Monte Carlo, variational, and deterministic model. See text for details. Bins of size 10 from 0 to 10, then 10 to 20, etc.

Figure 4 :

 4 Figure 4: Reconstruction results for random input images from the test set.

Table 2 :

 2 closed the Evaluation results for variations of PRN. Monte Carlo training uses samples from the unconditional or image-informed prior, while variational training relies on samples from the shape-conditioned approximate posterior. We report IoU under two occupancy probability thresholds τ. consistency loop between inferred point clouds and input images, by rendering point clouds as images and minimizing loss between such renderings and input images.Similarly to the previous class of models, these methods also enrich available training data by considering different forms of additional data: camera information, 2.5D sketches, etc. This, again, makes comparison problematic, since it is not always clear what the impact of the additional training data is.

	Dependencies	Sampling Deterministic Shape model IoU↑ (0.5) IoU↑ (0.4)
	p (v|z) p (z|i)	p (z|i)	63.7	65.0
	p (v|z, i) p (z|i)	p (z|i)	64.6	65.6
	p (v|z, i) p (z)	p (z)	64.0	65.0
	p (v|z) p (z|i)	q (z|v)	65.9	66.2
	p (v|z, i) p (z|i)	q (z|v)	64.8	65.3
	p (v|z) p (z|i)	q (z|v, i)	65.4	65.8
	p (v|z) p (z|i)	q (z|v)	63.4	63.7
	p (v|z) p (z|i)	q (z|v)	65.6	66.1

https://github.com/Regenerator/prns

Acknowledgements

This work has been partially supported by the grant "Deep in France" (ANR-16-CE23-0006) and LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).


Julia Chris Elrod
email: elrodc@gmail.com

Franc ¸ois Févotte

Accurate and Efficiently Vectorized Sums and Dot Products in

Keywords: Vectorization, Compensated Algorithms, Julia programming language

This paper presents an efficient, vectorized implementation of various summation and dot product algorithms in the Julia programming language. These implementations are available under an open source license in the AccurateArithmetic.jl Julia package.

Besides naive algorithms, compensated algorithms are implemented: the Kahan-Babuška-Neumaier summation algorithm, and the Ogita-Rump-Oishi simply compensated summation and dot product algorithms. These algorithms effectively double the working precision, producing much more accurate results while incurring little to no overhead, especially for large input vectors.

This paper also tries and builds upon this example to make a case for a more widespread use of Julia in the HPC community. Although the vectorization of compensated algorithms is no particularly simple task, Julia makes it relatively easy and straightforward. It also allows structuring the code in small, composable building blocks, closely matching textbook algorithms yet efficiently compiled.

I. INTRODUCTION

Computing the dot product of two vectors and, perhaps to a lesser extent, summing the elements of a vector, are two very common basic building blocks for more complex linear algebra algorithms. As such, any change in their performance is likely to affect the overall performance of scientific computing codes; any change in their accuracy is likely to induce a loss of reproducibility in overall computed results.

Such a lack of reproducibility is becoming more and more common as the performance of supercomputers increase, due to a combination of multiple factors. First, an increase in computational power allows for ever more expensive models (such as, for example, finer grids for the discretization of partial differential equations), which in turn generates larger vectors, whose summation is often more ill-conditioned. Second, the order in which operations are performed in modern systems is far from being deterministic and reproducible: vectorization, parallelism, and compiler optimizations are only a few examples of possible causes for unexpected changes in the execution order of a given program. Since floatingpoint (FP) arithmetic is not associative, such non-determinism in execution order entails non-reproducibilities in computed results. For these reasons, accurate summation and dot product algorithms are sometimes useful.

However, increasing the algorithms accuracy at the price of an increase in their computing cost is not acceptable in all cases. For programs that already use double precision (64bit FP numbers), increasing the precision is not an option, since quadruple precision (128-bit FP numbers) remains too expensive for most purposes 1 . On the other hand, hitting the memory wall means that many algorithms (and in particular BLAS1-like algorithms such as summation and dot product) are limited by the memory bandwidth in a wide range of vector lengths on modern systems. In turn, this means that CPUs are idling for a fraction of the computing time, which could instead be used to perform additional operations for free, as long as these operations do not involve additional memory accesses. Much work has therefore been devoted to computing sums and dot products in a way that is both fast and accurate, leading to efficient algorithms implemented in a variety of programming languages.

The work presented here comes from the observation2 that, as of November 2018, no such algorithm was currently implemented in an efficient way in the Julia programming language. This paper describes an attempt made to implement accurate summation and dot product algorithms in Julia, in a way that is as efficient as possible, and in particular makes use of SIMD vector units. All the code described here is available under an open source license in the AccurateArithmetic.jl package. This paper will also try to make the argument that writing such efficient implementations of compensated algorithms is a relatively easy and straightforward task thanks to the nice properties of the Julia language.

Part II of this paper takes the viewpoint of computational sciences and discusses some aspects related to the use of vector instructions in fast implementations of summation algorithms. Part III then takes the point of view of floatingpoint arithmetic, and reviews a subclass of accurate summation and dot product algorithms: compensated algorithms. Nothing particularly new is presented in these parts, except maybe some specific details about interactions between vectorization and compensation (paragraph III-C), that have to be dealt with in order to gain the full benefit of both worlds: speed and Algorithm 1 Scalar, naive summation 1: {Initialization:} 2: a ← 0 3: {Loop on vector elements:} 4: for i ∈ 1 : N do 5: a ← a ⊕ x i 6: end for 7: return a accuracy. Part IV gives some insight about the Julia implementation of such compensated algorithms, and discusses the performance obtained on some test systems. These results are summarized in part V, along with some closing notes and outlooks.

II. THE COMPUTATIONAL SCIENCE POINT OF VIEW

For the sake of simplicity, let us only discuss here the case of the summation algorithm:

sum(x) = N i=1 x i ,
where x ∈ R N . Since a dot product can also be expressed as a summation, no loss of generality is entailed by this choice:

dot(x, y) = N i=1 x i y i , where x ∈ R N and y ∈ R N .
Algorithm 1 presents the most naive possible algorithm for a summation, in its scalar version. An accumulator a stores partial sums. During each loop iteration, a new vector element x i is added to this accumulator. Note that this accumulation is performed using the computer arithmetic addition operation ⊕, which operates on FP numbers, rounds its result, and is therefore only an approximation of the mathematical + operator which operates on real numbers. In the remainder of this paper, similar notations will be used for other operators, such as denoting the FP subtraction.

In Algorithm 2, this algorithm is refined to make use of vector instructions. In the following, the system will be assumed to efficiently implement SIMD instructions on vectors (hereafter referred to as "packs" so that the "vector" term can unambiguously refer to mathematical vectors such as the operand x) of width W . For x86_64 architectures (which are the main target of the algorithms described here), the AVX extension for example supports efficiently operating on 256bit registers holding W = 4 double-precision FP numbers. All variables denoted in bold font represent such packs, and are meant to be stored and manipulated as SIMD registers. The main loop in Alg. 2 (lines 4-8) now operates on a SIMD accumulator, in which W vector elements are added at once at each iteration. The number of iterations is therefore divided by W . In case the vector length is not a multiple of W , a Algorithm 2 Vectorized, naive summation 1: {Initialization:} 2: a ← 0 3: {Loop on full packs:} 4: for j ∈ 1 : N W do 5:

i ← W (j -1) + 1

6: p ← (x i , x i+1 , x i+2 , x i+3 ) 7:
a ← a ⊕ p 8: end for a ← a ⊕ x i 14: end for 15: return a remainder has to be accounted for. There are multiple ways of doing this; in this example, the remaining elements are added in a scalar loop. Before doing this, the SIMD accumulator has to be reduced to a scalar: this is the purpose of the vsum function, which sums all W components of a SIMD register:

vsum(a) = a[1] ⊕ a[2] ⊕ . . . ⊕ a[W ],
where a[i] represents the i-th element in a SIMD pack a.

Although such an implementation is already much more complex than its scalar version, it is not enough to obtain optimal performance on modern hardware. Because of the loopcarried dependency of the accumulator and instruction latency in superscalar hardware architectures, full performance can only be obtained if multiple instructions can be running at the same time. This means breaking that sequential dependency between one iteration and the next; a traditional way to do that consists in partially unrolling the main loop. An illustration of this technique is given in Algorithm 3, where the loop is unrolled twice (U = 2). A consequence of this technique is that there are now as many SIMD accumulators as the level of partial unrolling U (2 in this case). These accumulators need to be summed together and reduced before proceeding with the accumulation of remaining scalar elements.

If the number of packs N W is not a multiple of U , an additional type of remainder needs to be handled: the full packs that were not handled in the partially unrolled loop. In the algorithm presented here, SIMD accumulator a 1 is used to accumulate these.

Figure 1 illustrates the gains to be expected from such a technique: for various vector sizes, the total summation time (normalized by the vector size) is plotted against the unrolling level. The parameter given in the x-axis is the base-2 logarithm of the unrolling level U . The figure shows that a value log 2 (U ) = 3 seems to be an optimal trade-off (on Algorithm 3 Vectorized, naive summation with partial unrolling (U = 2 in this example) 1: {Initialization:} 2: a 1 ← 0 3: a 2 ← 0 4: {Loop on full packs, unrolled twice:} 5: for j ∈ 1 : N 2 W do 6:

i 1 ← 2 W (j -1) + 1 7: p 1 ← (x i1 , x i1+1 , x i1+2 , x i1+3 ) 8: a 1 ← a 1 ⊕ p 1 9: i 2 ← 2 W (j -1) + W + 1 10: p 2 ← (x i2 , x i2+1 , x i2+2 , x i2+3 ) 11:
a 2 ← a 2 ⊕ p 2 12: end for 13: {Loop on remaining full packs:}

14: for j ∈ 2 N 2 W + 1 : N W do 15: i ← W (j -1) + 1 16: p ← (x i , x i+1 , x i+2 , x i+3 ) 17: a 1 ← a 1 ⊕ p 18: end for 19: {Reduction of SIMD accumulators:} 20: a 1 ← a 1 ⊕ a 2 21: a ← vsum(a 1 ) 22: {Loop on remaining elements:} 23: for j ∈ W N W + 1 : N do 24:
a ← a ⊕ x i 25: end for 26: return a this particular system): this corresponds to unrolling the loop U = 2 3 = 8 times. Depending on the vector size, the gains achieved can be as large as a 7× speed-up w.r.t. the version without any unrolling (log 2 (U ) = 0).

III. THE ARITHMETIC POINT OF VIEW

A. Accuracy and condition number

The naive summation algorithm presented above produces results of relatively poor quality when used on ill-conditioned vectors. Throughout this paper, the computer arithmetic will be assumed to follow the IEEE-754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF]. Under this hypothesis, the accuracy of the naive summation algorithm is given by

snaive (x) - N i=1 x i = K N N i=1 |x i |,
where s(x) denotes the result of naive summation on vector x ∈ R N , K N is a constant depending only on the vector size N , and is the machine epsilon for the FP format used in the computation. All examples in this paper use double precision (known as binary64 in IEEE-754), for which = 2 -53 1.11 × 10 -16 . In other words, the relative error on the naive summation results verifies

ε naive (x) = snaive (x) - N i=1 x i N i=1 x i = K N C(x) , (1)
where the condition number of the summation was introduced:

C(x) = N i=1 |x i | N i=1 x i .
Similar results can be obtained for the dot product. Following the line of reasoning introduced in [START_REF] Ogita | Accurate sum and dot product[END_REF], Figure 2 plots the relative error as a function of the condition number, in a log-log scale. Randomly generated input vectors of sizes ranging from 100 to 110 are produced by the algorithm described in [START_REF] Ogita | Accurate sum and dot product[END_REF] to feature condition numbers ranging from 1 to 10 45 . In the figure, errors lower than are arbitrarily set to ; conversely, when the relative error is more than 100% (i.e no digit is correctly computed anymore), the error is capped there in order to avoid affecting the scale of the graph too much. In Figure 2a, the "pairwise" line describes the accuracy of Julia's standard summation function (Base.sum), which uses a pairwise summation algorithm; the "naive" line corresponds to a naive summation implementation similar to the one described in Alg. 3 and available as AccurateArithmetic.sum_naive. Both algorithms behave as described by Equation (1): they start losing accuracy as soon as the condition number increases, computing only noise when the condition number exceeds 1/ 10 16 . A very similar behavior is observed in Figure 2b for the two naive dot product implementations: "blas" denotes Julia's default LinearAlgebra.dot function, which internally uses OpenBLAS's ddot implementation [START_REF] Wang | Augem: automatically generate high performance dense linear algebra kernels on x86 cpus[END_REF]; "naive" denotes a naive dot product implementation

• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • • • • • •• •• • •• ••• • •• • • • • ••• • • • • •• • • • • • • •• • • ••• • •••• ••• •• • • •• •• • •••• • • • • • •• • • • ••• • • •• • • • • • + + + + + +
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • •• • • • • • • • ••• ••• •• • ••• •• • • •• • •• • • • • • • • •• •• •• • • • •• •• • • • • • ••• • • •• • • •• • • •• • •• •• • • • • • • •• • • • •• • • • • Condition number 10

B. Accurate sum and dot product

As said above, much work has been devoted to developing more accurate and/or reproducible summation and dot product algorithms, of which an extensive survey would not fit the format of this document. For the sake of simplicity, let us consider here that a large fraction of accurate summation algorithms leverage a common and very simple idea: use a larger accumulator to store partial sums. In its crudest variant, this might consist in using mixed precision, where a 64-bit FP number accumulates partial results in the summation of 32-bit floats. More elaborate algorithms consider longer fixed-point accumulators (superaccumulators), capable of representing the result of every possible sum of FP numbers of a given format; this is for example the case in ExBLAS [START_REF] Iakymchuk | ExBLAS: Reproducible and accurate BLAS library[END_REF]. The xsum algorithm [START_REF] Radford | Fast exact summation using small and large superaccumulators[END_REF] is a variation on this concept. Algorithms such as those of ReproBLAS [START_REF] Demmel | Parallel reproducible summation[END_REF] use smaller accumulators, which focus on making the result of the summation independent of the order of the summands, but also feature increased accuracy w.r.t. naive summation.

In this work, we focus on a family of so-called "compensated algorithms", where the error introduced by each accumulation is computed using "error-free transforms" (EFTs), and re-introduced later in the computation. As an example, Alg. 4 describes the fast_two_sum algorithm [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF], which takes two FP numbers (a, b) as inputs and returns two FP numbers (x, e) such that

x = a ⊕ b,

x + e = a + b.

As such, e can be considered as the round-off error introduced in the FP addition a + b. Note that in the particular case of fast_two_sum, operands a and b have to be sorted by decreasing magnitude in order for the EFT to be correct.

The compensated summation algorithm due to Kahan, Babuška and Neumaier (KBN) [START_REF] Babuška | Numerical stability in mathematical analysis[END_REF], [START_REF] Kahan | Further remarks on reducing truncation errors[END_REF], [START_REF] Neumaier | Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen[END_REF] is probably the most famous in this category. A scalar implementation of the KBN algorithm in Julia is available in the KahanSummation.jl package. Ogita, Rump and Oishi showed in [START_REF] Ogita | Accurate sum and dot product[END_REF] that the KBN algorithm could be reinterpreted as a "cascaded summation" (Alg. 6) based on the fast_two_sum EFT. Since fast_two_sum inputs need to be sorted by magnitude in order for the transform to be error-free, Ogita, Rump and Oishi propose to build a similar cascaded summation based on the two_sum EFT [START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF] (Alg. 5), whose practical implementation does not require any test or branch to sort the inputs and should therefore be more suitable for efficient vectorization. This two_sum-based cascaded summation will be referred to as the ORO algorithm in the following. Compensated summation algorithms built with this technique (as well as compensated dot product algorithms built in a similar way) can be shown to produce results that are as accurate as if they had been computed using a naive algorithm at twice the working precision, rounding the final result to the nearest working precision value:

ε comp (x) = K N C(x) 2 .
(

) 2
For the sake of completeness, let us mention that EFTs can also be used to build Floating-Point Expansions (FPEs), in which a higher-precision number is represented as the (unevaluated) sum of multiple FP numbers. FPEs of size 2, also known as double-double [START_REF] Bailey | A fortran 90-based multiprecision system[END_REF], are at the core of the implementation of the XBLAS library [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] whose algorithms (σ, e) ← EFT_sum(σ, x i )

6:
τ ← τ ⊕ e 7: end for 8: return σ ⊕ τ are very similar in spirit to Ogita, Rump and Oishi's cascaded summation, except in the way lower-order terms are handled.

C. Vectorized implementation of compensated algorithms

The remainder of this paper concentrates on the efficient, vectorized implementation of compensated algorithms such as those introduced by Ogita, Rump and Oishi. Looking at the cascaded summation algorithm 6, one can observe the similarity to a naive summation (Alg. 1):

• the single accumulator a (Alg. 1, line 2) is replaced by an EFT-based accumulator (σ, τ ) (Alg. 6, line 2) • the accumulation operation itself (Alg. 1, line 5) is replaced by the use of an EFT to accumulate high-order terms in σ and compute an error (Alg. 6, line 5) followed by a simple accumulation of low-order terms in τ (line 6). Nevertheless, a few questions arise when trying to vectorize such an algorithm. Looking at the efficiently vectorized version of the naive summation (Alg. 3), it appears that there are quite a few additional places where intermediate results are summed: (acc+pack) lines 8, 11 and 17: addition of a pack of vector elements into a SIMD accumulator; (acc+acc) line 20: addition of a SIMD accumulators into another SIMD accumulator; (vsum(acc)) line 21: sum of all elements of a SIMD accumulator, yielding the initial value of a scalar accumulator; (acc+elem) line 24: addition of a (scalar) vector element to a scalar accumulator. While summing an ill-conditioned vector, a catastrophic cancellation might occur in any of these sums, so that special care must be taken to compensate for errors at all these locations so that the overall accuracy of the compensated algorithm is not lost. The (acc+elem) accumulation type is

(σ, τ ) ← (0, 0) for i ∈ 1 : W do (σ, e) ← EFT_sum(σ, σ[i]) τ ← (τ ⊕ τ [i]
) ⊕ e end for return (σ, τ ) exactly the same as what appears in a standard (scalar) cascaded summation, and should therefore be handled in the way described by Alg. 6. The (acc+pack) accumulation is the straightforward SIMD generalization of this operation, and can be handled in the same way, provided that vectorized EFT implementations are available.

In order to keep the "spirit" of the cascaded summation algorithm, we propose to handle the two remaining accumulation types (acc+acc) and (vsum(acc)) using algorithms 7 and 8 respectively, in such a way that high-order terms are always summed using an EFT, whose error term is simply accumulated with the low-order terms.

The soundness of this approach can be assessed by looking at the relative error of the results produced by vectorized implementations of the compensated algorithms. Looking at Figure 2a, it is clear that both compensated algorithms behave as predicted, with relative errors following (2): KBN and ORO give the exact same results, and only start losing accuracy for condition numbers in the order of 1/ 10 16 , computing meaningless results when the condition number reaches 1/ 2 10 32 . The compensated dot product algorithm labeled "dot oro" on figure 2b re-uses the same vectorized cascaded summation algorithm, using an additional two_prod EFT to compute error terms for the elementwise product of the two vectors, as proposed in [START_REF] Ogita | Accurate sum and dot product[END_REF]. Again, this compensated dot product implementation behaves accordingly to the theory, with relative errors following (2).

IV. JULIA IMPLEMENTATION & BENCHMARKING

A. Implementation in Julia

The full Julia implementation of all algorithms described here (naive and compensated versions of the summation and dot product) is available under an open source license in the AccurateArithmetic.jl package 3 . Any reader interesting in the details can therefore look at the code, and only some interesting features will be highlighted here. While doing so, we will also try and make the case that the Julia programming language is well suited to easily develop efficient implementation of such algorithms. In particular, it should be noted that the actual implementation closely follows the algorithmic description made here, and that all individual algorithms can be implemented in their own separate functions, that the Julia compiler will then assemble in an efficient way. This makes it easy to maximize code re-use, as well as to structure the code as building blocks, which will later allow for an easy exploration of new or unusual combinations of elementary components.

As a first example, Listing 1 gives the Julia implementation of the two_sum EFT (Alg. 5). This implementation, which is a nearly direct translation from the textbook algorithm, works for any type of inputs, including SIMD packs. The SIMDops.@explicit macro at the beginning arranges for elementary operators such as + orto be vectorized when they operate on SIMD packs, using the SIMDPirates.jl package 4 . Special care is taken to ensure that no further optimization is performed (such as fusing multiplications and additions) which would potentially break the EFT: this is the meaning of @explicit.

As a second, more complicated example, Listing 2 shows the partially unrolled loop corresponding to lines 4-12 in Alg. 3. Again, although the problem at hand is much more complicated, the Julia implementation closely follows the algorithmic description. The addition of a SIMD pack to an accumulator has been abstracted away in an add! function, which can have specific implementations (methods, in the Julia terminology) for each accumulator type. This way, the same function implementing Alg. 3 can be used as a building block for all kinds of summation algorithms (naive or compensated). Following the same kind of logic, the mechanism loading SIMD packs is designed in such a way that it can handle a 1-tuple containing the input vector in the case of the summation, but also load packs from two input vectors contained in a 2-tuple in the case of the dot product.

B. Performance evaluation

Although the implementation makes use of small, composable functional units, the Julia compiler specializes and optimizes the program for the specific combination of building blocks being used, so as to emit very efficient code.

Figure 3 plots, for various implementations of the sum and dot product, the run time (renormalized per element) against the vector size. These results have been obtained on a system featuring an Intel Core i5 6200U, which is a typical mobile AVX2-enabled Skylake-generation CPU. It can be observed that our naive implementation performs a bit better than Julia's standard pairwise summation for small vector sizes, but almost identically afterwards. Once vectors start having significant sizes (more than a few hundreds of elements), both implementations are memory bound, as expected of a typical BLAS1 operation. This explains why significant decreases in the performance can be observed when the vector can't fit into the L2 cache (around 30k elements, or 256kB on the test system) or the L3 cache (around 400k elements, or 3MB on the test system). Our explicitly vectorized naive summation remains memory bound in the L1 cache, where it gains a speed-up of about 2× with respect to Julia's automatically vectorized pairwise summation.

The ORO algorithm, when implemented with a suitable unrolling level (2 2 = 4 unrolled iterations in this case), is CPU-bound when vectors fit inside the cache. However, when vectors are too large to fit into the L3 cache, the implementation becomes memory-bound again, which means that it achieves the same performance as the standard summation. In other words, the improved accuracy is free for sufficiently large vectors. For smaller vectors, the accuracy comes with a slow-down no larger than 3× for vectors which fit in the L2 cache or 2.5× in the L3 cache.

The KBN algorithm performs less well on this system: although the fast_two_sum EFT uses fewer operations than two_sum, the need for a test sorting the arguments by magnitude makes it less suitable for a vectorized implementation on this system.

The same kind of conclusions can be reached for the dot product. In order to make an apples-to-apples comparison, the use of threads has been disabled in OpenBLAS. Naive implementations are always memory bound (the L1 cache boundary can be seen clearly in this case). The ORO compensated implementation incurs a slow-down no larger than 5× in the L1 cache, which decreases to around 2.5× in the L2 cache, 1.7× in the L3 cache, and no overhead outside the cache.

The same tests can be performed on other type of architectures, and lead to different conclusions. For example, Figure 4 presents the results obtained on a system powered by an Intel Xeon Silver 4110 CPU, of the same Skylake generation, but supporting the AVX512 extension. On this system, it is interesting to note that the KBN algorithm performs better than ORO. This can probably be attributed to this generation of CPUs better handling of the branches in fast_two_sum, thanks to the operands masking capabilities introduced in AVX512. Compensated summation algorithms remain memory bound much longer, with KBN achieving the same performance as naive algorithms even in the L3 cache. With respect to Julia's pairwise summation, a 4× slow-down is observed in the L1 cache, which decreases to 2.5× in the L2 cache. A 5 -10% overhead is observed outside the cache.

As for the dot product, the OpenBLAS implementation seems to be a bit more optimized than our naive implemen-tation, with a 10 -20% difference in speed overall. The compensated dot product implementation achieves the same performance as the naive one in the L3 cache, while incurring a 3.5× slow-down in the L2 cache and a 7% overhead outside the cache.

V. CLOSING NOTES

In summary, the work presented here tries to bring into the Julia ecosystem efficient, vectorized implementations of stateof-the-art compensated summation algorithms. All the code described here is available as an open source Julia package on github: https://github.com/JuliaMath/AccurateArithmetic.jl

The family of algorithms implemented here is described in [START_REF] Ogita | Accurate sum and dot product[END_REF]: the Kahan-Babuška-Neumaier summation algorithm, and the simply compensated Ogita-Rump-Oishi summation and dot product algorithms. In accordance to the theory, these algorithms are observed to effectively double the floating-point precision (yielding, for double-precision inputs, a result as accurate as if it had been naively computed using FP numbers with a 106-bit mantissa).

From a performance viewpoint, these implementations are shown to incur no significant penalty with respect to naive algorithms, as soon as input vectors are larger than the L3 cache. For smaller vectors fitting in the L2 cache, the slowdown is no larger than approximately 3×, the precise value depending on the architecture. As observed by Ogita, Rump and Oishi, the branches of the KBN algorithm make it a bit less efficient than the ORO algorithm on AVX2 architectures. However, benchmarks on AVX512 systems suggest that this might not always be the case, so that compensated algorithms based on the fast_two_sum EFT might have to be reconsidered.

The Julia implementation provided here is very compact and modular, which makes it a very suitable platform for further experiments. Future work might include implementations of other types of accurate summation and dot product algorithms, such as those based on fixed-point-like accumulators and used in other software packages such as ExBLAS, XBLAS, ReproBLAS or xsum. Alternatively, the compensated summation and dot product algorithms implemented here could be used as a basis to implement accurate versions of other BLAS functions like matrix decompositions (QR, Cholesky) or solving triangular linear systems.

An other potential direction for future work would be to focus on improving the accuracy of single-or half-precision algorithms (for example via the use of mixed precision) or to develop similar implementations specifically targeting GPUs.

9 :

 9 {Reduction of SIMD accumulator:} 10: a ← vsum(a) 11: {Loop on remaining elements:} 12: for j ∈ W N W + 1 : N do 13:

Fig. 1 :

 1 Fig. 1: Performance of a vectorized, naive summation as a function of the unrolling level. Computing times are normalized by the vector size.

Fig. 2 :

 2 Fig. 2: Comparison of the accuracy of various summation and dot product algorithms, as a function of the condition number

Algorithm 5 Algorithm 6 3 :

 563 two_sum error-free transform Require: (a, b) ∈ F 2 Ensure: x = a ⊕ b and x + e = a + b x ← a ⊕ b y ← x a e ← (a (x y)) ⊕ (b y) Cascaded summation. EFT_sum can be either fast_two_sum (which gives the KBN algorithm) or two_sum (which gives the ORO summation algorithm) 1: {Initialization:} 2: (σ, τ ) ← (0, 0) {Loop on vector elements:} 4: for i ∈ 1 : N do 5:

Algorithm 7 Algorithm 8

 78 Addition of two EFT-based SIMD accumulators In: two EFT-based SIMD accumulators a 1 = (σ 1 , τ 1 ) and a 2 = (σ 2 , τ 2 ) Out: EFT-based SIMD accumulator representing a 1 ⊕ a 2 (σ, e) ← EFT_sum(σ 1 , σ 2 ) τ ← (τ 1 ⊕ τ 2 ) ⊕ e return (σ, τ ) Summation of an EFT-based accumulator In: EFT-based SIMD accumulator a = (σ, τ ) of width W Out: EFT-based scalar accumulator a = vsum(a)

Listing 1 two_sum implementation in Julia 1 function 2 SIMDops

 12 two_sum(a::T, b::T) where {T}

Listing 2 1 # 6 for j in 1 : 8 # Unroll U times 9 11 # 12 # 15 #

 216189111215 Partial unrolling of the vectorized loop N÷($W * $U) 7 Base.Cartesian.@nexprs $U i -> begin 10 Load a pack of width W and type T starting at position (px+offset) 13 p_i = vload.(Vec{$W,$T}, px.+offset) 14 Add pack p_i into accumulator acc_i 16 add!(acc_i, p_i...) 17 18 # Update offset for the next pack 19 offset += $W * $sizeT

Fig. 3 :

 3 Fig. 3: Comparison of the performance of various summation and dot product implementations (Intel Core i5 6200U)

Fig. 4 :

 4 Fig. 4: Comparison of the performance of various summation and dot product implementations (Intel Xeon Silver 4110)

at least on x86-64 architectures, where quadruple precision is implemented in software

https://discourse.julialang.org/t/floating-point-summation/17785

https://github.com/JuliaMath/AccurateArithmetic.jl

https://github.com/chriselrod/SIMDPirates.jl

ACKNOWLEDGEMENTS

The authors feel indebted to the members of the Julia discourse community 5 and in particular Jeffrey Sarnoff and Laurent Plagne. This work would not have started without their encouragement. Their insightful comments and help greatly contributed to increase the quality of this paper.

APPENDIX REPRODUCIBILITY RESOURCES

The AccurateArithmetic.jl Julia package implements all algorithms described in this paper. It also contains everything needed to reproduce the results shown in this paper. Below is a full list of the actions to be taken by anyone wanting to try and reproduce the results on their own Linux system. Commands prefixed with a "sh>" prompt are to be entered in a shell; commands prefixed with a "julia>" prompt are to be entered in a Julia interactive session. The Julia programming language should have been downloaded and installed beforehand from: https://julialang.org/downloads/ Any version above v1.0.0 should work; the results presented in this paper have been obtained using Julia v1.1.1.

Be aware that step (4) in this procedure might take a few hours to complete. Afterwards, all measurements should be available as JSON files in the AccurateArithmetic.jl/test directory. After step ( 5), all figures showed in this paper should be available as PDF files in the same directory.

Should anyone actually run this procedure, the authors would really like to get a copy of their results, so as to gain more insight about the behavior of compensated algorithms on a variety of systems.

# (1) Get the repository sh> cd /tmp sh> git clone --branch paper-correctness-2019 \ https://github.com/JuliaMath/AccurateArithmetic.jl.git sh> cd AccurateArithmetic.jl/test # (2) Install required dependencies # and run the test suite sh> julia --project julia> using Pkg julia> pkg"instantiate" julia> pkg"test" julia> exit() # (3) Additional dependencies for performance tests sh> julia --project= julia> using Pkg julia> pkg"add BenchmarkTools Plots Printf Statistics Test" julia> exit() # (4) Run the performance tests sh> julia --project -O3 -L perftests.jl -e 'run_tests()' # (5) Plot the graphs # (this step can be run on another machine, provided # that all JSON files resulting from the previous # step are copied) sh> julia --project -L perfplots.jl -e 'plot_results()'


Semi-supervised Learning using Deep Generative Models and Metric Embedding Auxiliary Task

Jhosimar George Arias Figueroa jariasf03@gmail.com

Motivation

In many real world applications, it is relatively easy to acquire large amounts of unlabeled data. However, their corresponding labels often requires the effort of a human expert and special devices. This bottleneck in the labeling results in small amounts of labeled data and large amounts of unlabeled data. Therefore, being able to utilize the surplus unlabeled data is desirable.

Background

Variational Autoencoder Evidence Lower Bound (ELBO)

Approximate p θ (z|x) with variational posterior q ϕ (z|x), define lower bound on log p(x)

Intuition

▶ Our encoder transforms the data into feature representations, from these representations we infer a Gaussian latent variable to provide more robust set of features and a discrete latent variable to provide category assignments.

▶ We take advantage of the small amount of labeled data to perform the category assignments for the unlabeled data.

▶ The decoder uses the gaussian representations and category information for reconstruction.

Probabilistic Model

Inference Model Generative Model

Evidence Lower Bound (ELBO)

The first term is the reconstruction loss (L R ). The second and third terms work as regularizers of the categorical (L C ) and Gaussian (L G ) distributions respectively.

Auxiliary Task Loss Assignment Loss (L A )

▶ In order to assign labels to the unlabeled data, we use the Distance-Weighted k-nearest neighbor (k-NN) on the feature space.

▶ Once we obtained the assignments for the unlabeled data, we use the Cross-entropy loss.

Metric Embedding Loss (L M )

▶ The goal of this loss is to regularize the feature space, in such a way that the distance between features that belong to the same category should be smaller than features of different categories.

▶ We perform metric learning via the lifted structured loss. LatinX in AI Workshop @ ICML 2019, Long Beach, CA, USA.

L M = 1 2| P | ∑ (i,j)∈ P max(0, J (f i , f j )) 2 J (f i , f j ) = D(f i , f j ) + max { max (i,k)∈ N α -D(f i , f k ), max
https://github.com/jariasf/semisupervised-vae-metric-embedding

  proposed model. The encoder and decoder are given by Convolutional Neural Networks. The red diamonds represent the reconstruction loss (L R ), categorical loss (L C ), gaussian loss (L G ), assignment loss (L A ) and metric embedding loss (L M ).Quantitative EvaluationMNIST SVHNPerformance of our probabilistic model after adding the assignment (L A ) and metric embedding (L M ) loss at different epochs on the MNIST and SVHN datasets.Semi-supervised test error (%) on MNIST, SVHN and CIFAR-10 datasets for randomly and evenly distributed labeled data. Colored rows denote methods that use Bayesian approaches.