Transformers in Self-Supervised Monocular Depth Estimation with Unknown Camera Intrinsics

The advent of autonomous driving and advanced driver assistance systems necessitates continuous developments in computer vision for 3D scene understanding. Self-supervised monocular depth estimation, a method for pixel-wise distance estimation of objects from a single camera without the use of ground truth labels, is an important task in 3D scene understanding. However, existing methods for this task are limited to convolutional neural network (CNN) architectures. In contrast with CNNs that use localized linear operations and lose feature resolution across the layers, vision transformers process at constant resolution with a global receptive field at every stage. While recent works have compared transformers against their CNN counterparts for tasks such as image classification, no study exists that investigates the impact of using transformers for self-supervised monocular depth estimation. Here, we first demonstrate how to adapt vision transformers for self-supervised monocular depth estimation. Thereafter, we compare the transformer and CNN-based architectures for their performance on KITTI depth prediction benchmarks, as well as their robustness to natural corruptions and adversarial attacks, including when the camera intrinsics are unknown. Our study demonstrates how transformer-based architecture, though lower in run-time efficiency, achieves comparable performance while being more robust and generalizable.

INTRODUCTION
There have been rapid improvements in scene understanding for robotics and advanced driver assistance systems (ADAS) over the past years. This success is attributed to the use of Convolutional Neural Networks (CNNs) within a mostly encoder-decoder paradigm. Convolutions provide spatial locality and translation invariance which has proved useful for image analysis tasks. The encoder, often a convolutional Residual Network (ResNet) (He et al., 2016), learns feature representations from the input and is followed by a decoder which aggregates these features and converts them into final predictions. However, the choice of architecture has a major impact on the performance and generalizability of the task.
While CNNs have been the preferred architecture in computer vision, transformers have also recently gained traction  motivated by their success in natural language processing (Vaswani et al., 2017). Notably, they have also outperformed CNNs for object detection (Cara Equal contribution ion et al., 2020) and semantic segmentation (Zheng et al., 2021). This is also reflected in methods for monocular dense depth estimation, a pertinent task for autonomous planning and navigation, where supervised transformer-based methods (Li et al., 2020;Ranftl et al., 2021) have been proposed as an alternative to supervised CNN-based methods Aich et al., 2021). However, supervised methods require extensive RGB-D ground truth collected from costly LiDARs or multi-camera rigs. Instead, self-supervised methods have increasingly utilized concepts of Structure from Motion (SfM) with known camera intrinsics to train monocular depth and ego-motion estimation networks simultaneously (Guizilini et al., 2020;Lyu et al., 2020;Chawla et al., 2021). While transformer ingredients such as attention have been utilized for self-supervised depth estimation (Johnston and Carneiro, 2020), most methods are nevertheless limited to the use of CNNs that have localized linear operations and lose feature resolution during downsampling to increase their limited receptive field (Yang et al., 2021). We readapt modules from Dense Prediction Transformer (DPT) and Monodepth2 to be trained with appearance-based losses for self-supervised monocular depth, ego-motion, and intrinsics estimation. ductive biases allow for more globally coherent predictions with different layers attending to local and global features simultaneously . However, transformers require more training data and can be more computationally demanding (Caron et al., 2021). While multiple studies have compared transformers against CNNs for tasks such as image classification (Raghu et al., 2021;Bhojanapalli et al., 2021), no study exists that evaluates the impact of transformers in self-supervised monocular depth estimation, including when the camera intrinsics may be unknown.
In this work, we conduct a comparative study between CNN-and transformer-based architectures for self-supervised monocular depth estimation. Our contributions are as follows: • We demonstrate how to adapt vision transformers for self-supervised monocular depth estimation by implementing a method called Monocular-Transformer SfMLearner (MT-SfMLearner).
• We compare MT-SfMLearner and CNNs for their performance on the KITTI monocular depth Eigen Zhou split (Eigen et al., 2014) and the online depth prediction benchmark (Geiger et al., 2013).
• We investigate the impact of architecture choices for the individual depth and ego-motion networks on performance as well as robustness to natural corruptions and adversarial attacks.
• We also introduce a modular method that simultaneously predicts camera focal lengths and principal point from the images themselves and can easily be utilized within both CNN-and transformerbased architectures.
• We study the accuracy of intrinsics estimation as well as its impact on the performance and robustness of depth estimation.
• Finally, we also compare the run-time computational and energy efficiency of the architectures for depth and intrinsics estimation.
MT-SfMLearner provides real-time depth estimates and illustrates how transformer-based architecture, though lower in run-time efficiency, can achieve comparable performance as the CNN-based architectures while being more robust under natural corruptions and adversarial attacks, even when the camera intrinsics are unknown. Thus, our work presents a way to analyze the trade-off between the performance, robustness, and efficiency of transformer-and CNN-based architectures for depth estimation.
For supervised monocular depth estimation, Dense Prediction Transformer (DPT) (Ranftl et al., 2020) uses ViT as the encoder with a convolutional decoder and shows more coherent predictions than CNNs due to the global receptive field of transformers. TransDepth (Yang et al., 2021) additionally uses a ResNet projection layer and attention gates in the decoder to induce the spatial locality of CNNs for supervised monocular depth and surface-normal estimation. Lately, some works have inculcated elements of transformers such as self-attention (Vaswani et al., 2017) in self-supervised monocular depth estimation (Johnston and Carneiro, 2020;Xiang et al., 2021). However, there has been no investigation of transformers to replace the traditional CNN-based methods (Godard et al., 2019;Lyu et al., 2020) for self-supervised monocular depth estimation.
Moreover, self-supervised monocular depth estimation still requires prior knowledge of the camera intrinsics (focal length and principal point) during training, which may be different for each data source, may change over time, or be unknown a priori (Chawla et al., 2020). While multiple approaches to supervised camera intrinsics estimation have been proposed (Lopez et al., 2019;Zhuang et al., 2019), not many self-supervised approaches exist (Gordon et al., 2019). Therefore, we investigate the impacts of transformer architectures on self-supervised monocular depth estimation for their performance, robustness, and run-time efficiency, even when intrinsics are unknown.

METHOD
Our objective is to study the effect of utilizing vision transformers for self-supervised monocular depth estimation in contrast with the contemporary methods that utilize Convolutional Neural Networks.
Given a set of n images from a video sequence, we simultaneously train depth and ego-motion prediction networks. The inputs to the networks are a sequence of temporally consecutive RGB image triplets {I −1 , I 0 , I 1 } ∈ R H×W ×3 . The depth network learns the model f D : R H×W ×3 → R H×W to output dense depth (or disparity) for each pixel coordinate p of a single image. Simultaneously, the ego-motion network learns the model f E : R 2×H×W ×3 → R 6 to output relative translation (t x ,t y ,t z ) and rotation (r x , r y , r z ) forming the affine transformation RT 0 1 ∈ SE(3) between a pair of overlapping images. The predicted depthD and ego-motionT are linked together via the perspective projection model, that warps the source images I s ∈ {I −1 , I 1 } to the target image I t ∈ {I 0 }, with the camera intrinsics described by K. This process is called view synthesis, as shown in Figure 1. We train the networks using the appearance-based photometric loss between the real and synthesized target images, as well as a smoothness loss on the depth predictions (Godard et al., 2019).

Architecture
Here we describe Monocular Transformer Structure from Motion Learner (MT-SfMLearner), our transformer-based method for self-supervised monocular depth estimation ( Figure 1).

Depth Network
For the depth network, we readapt the Dense Prediction Transformer (DPT) (Ranftl et al., 2020) for self-supervised learning, with a DeiT-Base  in the encoder. There are five components of the depth network: • An Embed module, which is a part of the encoder, takes an image I ∈ R H×W ×3 , and converts non-overlapping image patches of size p × p into where d = 768 for DeiT-Base. This is implemented as a large p × p convolution with stride s = p where p = 16. The output from this module is concatenated with a readout token of the same size as the remaining tokens.
• The Transformer block, that is also a part of the encoder, consists of 12 transformer layers which process these tokens with multi-head selfattention (MHSA) (Vaswani et al., 2017) modules. MHSA processes inputs at constant resolution and can simultaneously attend to global and local features. The subscripts of DN refer to the transformer layer from which the respective Reassemble module takes its input (see Figure  1). Input image size is H ×W , p refers to the patch size, N p = H ·W /p 2 refers to the number of patches from the image, and d refers to the feature dimension of the transformer features.
• Four Reassemble modules in the decoder, which are responsible for extracting image-like features from the 3 rd , 6 th , 9 th , and 12 th (final) transformer layers by dropping the readout token and concatenating the remaining tokens in 2D. This is followed by pointwise convolutions to change the number of channels, and transpose convolution in the first two reassemble modules to upsample the representations (corresponding to T 3 and T 6 in Figure 1). To make the transformer network comparable to its convolutional counterparts, we increase the number of channels in the pointwise convolutions of the last three Reassemble modules by a factor of 4 with respect to DPT. The exact architecture of the Reassemble modules can be found in Table 1.
• Four Fusion modules in the decoder, based on RefineNet (Lin et al., 2017). They progressively fuse information from the Reassemble modules with information passing through the decoder, and upsample the features by 2 at each stage. Unlike DPT, we enable batch normalization in the decoder as it was found to be helpful for selfsupervised depth prediction. We also reduce the number of channels in the Fusion block to 96 from 256 in DPT.
• Four Head modules at the end of each Fusion module to predict depth at 4 scales following previous self-supervised methods (Godard et al., 2019). Unlike DPT, the Head modules use 2 convolutions instead of 3 as we found no difference in performance. For further details of the Head modules, refer to Table 2.

Ego-Motion Network
For the ego-motion network, we adapt DeiT-Base  in the encoder. Since the input to the transformer for the ego-motion network consists of two images concatenated along the channel dimension, we repeat the embedding layer accordingly.
We use a Reassemble module to pass transformer tokens to the decoder. For details on the structure of this Reassemble module, refer to Table 1. We adopt the decoder for the ego-motion network from Mon-odepth2 (Godard et al., 2019). When both depth and ego-motion networks use transformers as described above, we refer to the resulting architecture as Monocular Transformer Structure from Motion Learner (MT-SfMLearner).

Appearance-based Losses
Following contemporary self-supervised monocular depth estimation methods, we adopt the appearancebased losses and an auto-masking procedure from the CNN-based Monodepth2 (Godard et al., 2019) for the above described transformer-based architecture as well. We employ a photometric reprojection loss composed of the pixel-wise 1 distance and the Structural Similarity (SSIM) between the real and synthesized target images, along with a multi-scale edgeaware smoothness loss on the depth predictions. We also use auto-masking to disregard the temporally stationary pixels in the image triplets. Furthermore, we reduce texture-copy artifacts by calculating the total loss after upscaling the depths, predicted at 4 scales, from intermediate decoder layers to the input resolution.

Intrinsics
Accurate camera intrinsics given by are essential to self-supervised depth estimation as can be seem from Equation 1. However, the intrinsics may vary within a dataset with videos collected from different camera setups, or over a long period of time. These parameters can also be unknown for crowdsourced datasets. We address this by introducing an intrinsics estimation module. We modify the ego-motion network, which takes a pair of consecutive images as input, and learns to estimate the focal length and principal point along with the translation and rotation. Specifically, we add a convolutional path in the ego-motion decoder to learn the intrinsics. The decoder features before activation from the penultimate layer are passed through a global average pooling layer, followed by two branches of pointwise convolutions to reduce the number of channels from 256 to 2. One branch uses a softplus activation to estimate focal lengths along x and y axes as the focal length is always positive. The other branch doesn't use any activation to estimate the principal point as it has no such constraint. Note that the ego-motion decoder is the same for both the convolutional as well as transformer architectures. Consequently, the intrinsics estimation method can be modularly utilized with both architectures. Figure 1 demonstrates MT-SfMLearner with learned intrinsics.

RESULTS
In this section, we perform a comparative analysis between our transformer-based method, MT-SfMLearner, and the existing CNN-based methods for self-supervised monocular depth estimation. We also perform a contrastive study on the architectures for the depth and ego-motion networks to evaluate their impact on the prediction accuracy, robustness to natural corruptions and adversarial attacks, and the run-time computational and energy efficiency. Finally, we analyze the correctness and run-time efficiency of intrinsics predictions, and also study its impact on the accuracy and robustness of depth estimation.

Dataset
We report all results on the Eigen Split (Eigen et al., 2014) of KITTI (Geiger et al., 2013) dataset after removing the static frames as per (Zhou et al., 2017), unless stated otherwise. This split contains 39, 810 training, 4424 validation, and 697 test images, respectively. All results are reported on the per-image scaled dense depth prediction without post-processing (Godard et al., 2019), unless stated otherwise.

Training Settings
The networks are implemented in PyTorch (Paszke et al., 2019) and trained on a TeslaV100 GPU for 20 epochs at a resolution of 640 × 192 with batchsize 12. MT-SfMLearner is further trained at 2 more resolutions -416 × 128 and 1024 × 320, with batchsizes of 12 and 2, respectively for experiments in Section 4.2. The depth and ego-motion encoders are initialized with ImageNet (Deng et al., 2009) pre-trained weights. We use the Adam (Kingma and Ba, 2014) optimizer for CNN-based networks (in Sections 4.3 and 4.4) and AdamW (Loshchilov and Hutter, 2017) optimizer for transformer-based networks with initial learning rates of 1e −4 and 1e −5 , respectively. The learning rate is decayed after 15 epochs by a factor of 10. Both optimizers use β 1 = 0.9 and β 2 = 0.999.

Depth Estimation Performance
First we compare MT-SfMLearner, where both depth and ego-motion networks are transformer-based, with the existing fully convolutional neural networks for their accuracy on self-supervised monocular depth estimation. Their performance is evaluated using metrics from (Eigen et al., 2014) up to a fixed range of 80 m as shown in Table 3. We compare the methods at different input image sizes clustered into categories of Low-Resolution (LR), Medium-Resolution (MR), and High-Resolution (HR). We do not compare against methods that use ground-truth semantic labels during training. All methods assume known groundtruth camera intrinsics.
We observe that MT-SfMLearner is able to achieve comparable performance at all resolutions under error as well as accuracy metrics. This includes methods that also utilize a heavy encoder such as ResNet-101 (Johnston and Carneiro, 2020) and Pack-Net (Guizilini et al., 2020).
Online Benchmark: We also measure the performance of MT-SfMLearner on the KITTI Online Benchmark for depth prediction 1 using the metrics from (Uhrig et al., 2017). We train on an image size of 1024 × 320, and add the G2S loss (Chawla et al., 2021) for obtaining predictions at metric scale. Results ordered by their rank are shown in Table 4. The performance of MT-SfMLearner is on par with stateof-the-art self-supervised methods, and outperforms

Contrastive Study
We saw in the previous section that MT-SfMLearner performs competently on independent and identically distributed (i.i.d.) test set with respect to the stateof-the-art. However, networks that perform well on an i.i.d. test set may still learn shortcut fea-  (Geirhos et al., 2020). Since self-supervised monocular depth estimation networks concurrently train an ego-motion network (see Equation 1), we investigate the impact of each network's architecture on depth estimation. We consider both Convolutional (C) and Transformer (T) networks for depth and ego-motion estimation. The resulting four combinations of (Depth Network, Ego-Motion Network) architectures are (C, C), (C, T), (T, C), and (T, T), ordered on the basis of their increasing influence of transformers on depth estimation. To compare our transformer-based method fairly with convolutional networks, we utilize Monodepth2 (Godard et al., 2019) with a ResNet-101 (He et al., 2016) encoder. All four combinations are trained thrice using the settings described in Section 4.1 for an image size of 640 × 192. All combinations assume known ground-truth camera intrinsics.

Performance
For the four combinations, we report the best performance on i.i.d. in Table 5, and visualize the depth predictions for the same in Figure 2. The i.i.d. test set used for comparison is same as in Section 4.2.
We observe from Table 5 that the combination of transformer-based depth and ego-motion networks i.e MT-SfMLearner performs best under two of the error metrics as well as two of the accuracy metrics. The remaining combinations perform comparably on all the metrics.
From Figure 2, we observe more uniform estimates for larger objects like vehicles, vegetation, and  buildings when depth is learned using transformers. Transformers are also less affected by reflection from windows of vehicles and buildings. This is likely because of the larger receptive fields of the self-attention layers, which lead to more globally coherent predictions. On the other hand, convolutional networks produce sharper boundaries, and perform better on thinner objects such as traffic signs and poles. This is likely because of the inherent inductive bias for spatial locality present in convolutional layers.

Robustness
While the different combinations perform comparably on the i.i.d. dataset, they may differ in robustness and generalizability. Therefore, we study the impact of natural corruptions and adversarial attacks on the depth performance using the following: • Natural corruptions. Following (Hendrycks and Dietterich, 2019) and ( (Madry et al., 2018) at attack strength ε ∈ {0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0}. The gradients are calculated with respect to the training loss. Following (Kurakin et al., 2016), the perturbation is accumulated over min(ε + 4, 1.25 · ε ) iterations with a step-size of 1. When the test image is from the beginning or end of a KITTI sequence, the appearance-based loss is only calculated for the feasible pair of images.
Inspired by (Wong et al., 2020) Figure 4 further shows a general improvement in adversarial robustness when learning either depth or ego-motion with transformers instead of convolutional networks. Finally, Figure 5 shows an improvement in robustness to symmetrically flipped adversarial attacks when depth is learned using transformers instead of convolutional networks. Furthermore, depth estimation is most robust when both depth and ego-motion are learned using transformers.

RMSE
Therefore, MT-SfMLearner, where both depth and ego-motion and learned with transformers, provides the highest robustness and generalizability, in line with the studies on image classification (Paul and Chen, 2021;Bhojanapalli et al., 2021). This can be attributed to their global receptive field, which allows for better adjustment to the localized deviations by accounting for the global context of the scene.

Intrinsics
Here, we evaluate the accuracy of our proposed method for self-supervised learning of camera intrinsics and its impact on the depth estimation performance. As shown in Table 6, the percentage error for intrinsics estimation is low for both convolutional and transformer-based methods trained on an image size of 640 × 192. Moreover, the depth error as well accuracy metrics are both similar to when the ground truth intrinsics are known a priori. This is also observed in Figure 6 where the learning of intrinsics causes no artifacts in depth estimation.
We also evaluate the models trained with learned intrinsics on all 15 natural corruptions as well as on PGD and symmetrically flipped adversarial examples. We report the mean RMSE (µRMSE) across all corruptions in Table 7. The RMSE for depth estimation on adversarial examples generated by PGD method for all strengths is shown in Figure 7. The mean RMSE (µRMSE) across all attack strengths for horizontally flipped and vertically flipped adversarial ex-    amples is shown in Table 8. We observe that the robustness to natural corruptions and adversarial attacks is maintained by both architectures when the intrinsics are learned simultaneously. Furthermore, similar to the scenario with known ground truth intrinsics, MT-SfMLearner with learned intrinsics has higher robustness than its convolutional counterpart.

Efficiency
We further compare the networks on their computational and energy efficiency to examine their suitabil-  ity for real-time applications.
In Table 9, we report the mean inference speed in frames per second (fps) and the mean inference energy consumption in Joules per frame for depth and intrinsics estimation for both architectures. These metrics are computed over 10, 000 forward passes at a resolution of 640 × 192 on an NVidia GeForce RTX 2080 Ti.
Both architectures run depth and intrinsics estimation in real-time with an inference speed > 30 fps. However, the transformer-based method consumes higher energy and is computationally more de- manding than its convolutional counterpart.

CONCLUSION
This work is the first to investigate the impact of transformer architecture on the SfM inspired self-supervised monocular depth estimation. Our transformer-based method MT-SfMLearner performs comparably against contemporary convolutional methods on the KITTI depth prediction benchmark. Our contrastive study additionally demonstrates that while CNNs provide local spatial bias, especially for thinner objects and around boundaries, transformers predict uniform and coherent depths, especially for larger objects due to their global receptive field. We observe that transformers in the depth network result in higher robustness to natural corruptions, and transformers in both depth, as well as the ego-motion network, result in higher robustness to adversarial attacks. With our proposed approach to self-supervised camera intrinsics estimation, we also demonstrate how the above conclusions hold even when the focal length and principal point are learned along with depth and ego-motion. However, transformers are computationally more demanding and have lower energy efficiency than their convolutional counterparts. Thus, we contend that this work assists in evaluating the trade-off between performance, robustness, and efficiency of self-supervised monocular depth estimation for selecting the suitable architecture.


Query Processing on Tensor Computation Runtimes

The huge demand for computation in artificial intelligence (AI) is driving unparalleled investments in hardware and software systems for AI. This leads to an explosion in the number of specialized hardware devices, which are now offered by major cloud vendors. By hiding the low-level complexity through a tensor-based interface, tensor computation runtimes (TCRs) such as PyTorch allow data scientists to efficiently exploit the exciting capabilities offered by the new hardware. In this paper, we explore how database management systems can ride the wave of innovation happening in the AI space. We design, build, and evaluate Tensor Query Processor (TQP): TQP transforms SQL queries into tensor programs and executes them on TCRs. TQP is able to run the full TPC-H benchmark by implementing novel algorithms for relational operators on the tensor routines. At the same time, TQP can support various hardware while only requiring a fraction of the usual development effort. Experiments show that TQP can improve query execution time by up to 10$\times$ over specialized CPU- and GPU-only systems. Finally, TQP can accelerate queries mixing ML predictions and SQL end-to-end, and deliver up to 9$\times$ speedup over CPU baselines.

INTRODUCTION
DBMS vendors have delivered constant performance improvement for decades by evolving software to keep up with Moore's law while influencing hardware development through close relationships with manufacturers. While data volumes and demand for analytics are growing faster than ever [129], the performance improvement on CPU has slowed down [136]. However, the count of processor transistors has continued to grow over the last decade, as hardware manufacturers adopted first multi-core CPU architectures and then augmented their computing platforms with specialized components such as GPUs, FPGAs, compression and encryption chips, DSPs, and neural-network (NN) accelerators. Although DBMS builders have taken advantage of multi-core and SIMD instructions effectively [76,109,146], the explosion in the number of specialized hardware Work done while Dong, Supun, and Konstantinos were at Microsoft. This work is licensed under the Creative Commons BY-NC-ND 4.0 International License. Visit https://creativecommons.org/licenses/by-nc-nd/4.0/ to view a copy of this license. For any use beyond those covered by this license, obtain permission by emailing info@vldb.org. Copyright is held by the owner/author(s). Publication rights licensed to the VLDB Endowment. Proceedings of the VLDB Endowment, Vol. 15 components, each with different characteristics and programming abstractions, makes it challenging to support all the exciting capabilities that these new powerful devices can offer.
On the other hand, the huge demand for computation in artificial intelligence (AI) [59], combined with the market fever for AI, is driving unparalleled investments in new hardware and software for AI. Hardware makers (e.g., Intel [62], Apple [34], Xilinx [142], AMD [33]), cloud vendors (e.g., Amazon [37], Microsoft [48], Google [72]), startups (e.g., Graphcore [6], Sambanova [11], Cerebras [4]), and car companies like Tesla [135] are investing heavily in this space. Venture capitals alone are pouring nearly $2B a quarter on special hardware for AI, aiming for a market expected to exceed $200B a year by 2025 [130]. On the software side, companies and open source communities are rallying behind a small number of big efforts (e.g., PyTorch [9], TensorFlow [31], TVM [46]). The combination of investments in specialized hardware and large software communities focusing on performance allows these efforts to thrive. Our realization is that the ML community has made hardware accelerators accessible to nonspecialists (e.g., data scientists). The fact that the most popular ML frameworks are open-source, creates a virtuous cycle whereby any hardware vendor interested in the ML space must support these frameworks well to get adoption. At the same time, these large open source communities successfully tackle the labor-intensive problem of providing specialized kernels for various hardware, e.g., a month after Apple M1 was announced, TVM outperformed Apple's CoreML by 2× [134]. Hardware vendors can directly improve the kernels' performance or the hardware itself [21,22,25]. This further helps adoption since the performance improves at each new software and hardware release.
We argue that the best path forward for analytical DBMSs is to embrace this tectonic shift and take advantage of the groundswell of new hardware and software targeting AI workloads. To demonstrate the viability of this idea, we propose and prototype a new query processor which runs SQL queries atop tensor computation runtimes (TCRs) such as PyTorch, TVM, and ONNX Runtime [23]. We name our prototype Tensor Query Processor (TQP). TQP transforms a SQL query into a tensor program and executes it on TCRs. To our knowledge, TQP is the first query processor built atop TCRs. Careful architectural and algorithmic design enables TQP to: (1) deliver significant performance improvements over popular CPU-based data systems, and match or outperform custombuilt solutions for GPUs; (2) demonstrate portability across a wide range of target hardware and software platforms; and (3) achieve all the above with parsimonious and sustainable engineering effort.
The above might appear surprising as specialized hardware accelerators are notoriously hard to program, requiring much customization to extract the best performance. Furthermore, their programming abstractions differ sufficiently to make our goals of performance (G1), portability (G2), and parsimonious engineering effort (G3) seemingly hard to reconcile. However, the key is a compilation layer and a set of novel algorithms, which can map the classical database abstraction to the prevalent one in machine learning (ML), i.e., mapping relational algebra to tensor computations. This allows us to free-ride on existing labor-intensive efforts from the ML community to port and optimize TCRs across all the new specialized hardware platforms. The initial performance results are encouraging. On GPU, TQP is able to outperform open-source GPU databases in terms of query execution time. On CPU, TQP outperforms Spark [145], and it is comparable to a state-of-the-art vectorized engine, DuckDB [117], for several queries. Furthermore, when ML and SQL queries are used in concert, TQP is able to provide end-to-end acceleration for a 9× speedup over CPU baselines.
Pursuing our goals of portability and parsimonious engineering effort, we make a deliberate decision to target existing tensor APIs rather than customize lower-level operators. This decision potentially leaves some performance on the table but leads to a very sustainable long-term play, as TQP benefits from any performance enhancement and optimization added to the underlying software and hardware (e.g., [21]). To validate this proposition, we run TQP on several different hardware settings: from CPUs, to discrete GPUs, to integrated GPUs (Intel and AMD), to NN-accelerators (TPUs [72]), and web browsers. Furthermore, TQP is able to run the full TPC-H benchmark on both CPU and GPU with just around 8,000 lines of code-this is quite an achievement considering that until 2021 no GPU database was able to run all the 22 TPC-H queries [84].
Contributions. This paper makes the following core contributions: • We propose Tensor Query Processor (TQP) that comprises a collection of algorithms and a compiler stack for transforming relational operators into tensor computations. • With TQP, we demonstrate that the tensor interface of TCRs is expressive enough to support all common relational operators. • We evaluate the TQP approach extensively against state-of-theart baselines on the TPC-H benchmark.
Organization. §2 introduces some background on TCRs. §3 summarizes the challenges and the design choices we make. §4 introduces TQP, and §5 describes the algorithms used to implement several key relational operators with tensor programs. Experiments are in §6. Related works are in §7. The paper is concluded by §8.

BACKGROUND
In this section, we summarize the system support for tensor computation ( §2.1), and provide a taxonomy of the tensor operations used throughout the paper ( §2.2).

Tensor Computation Runtimes (TCRs)
The last years have witnessed an increase in the popularity of ML models based on NNs [60]. While in the heydays, these models were implemented manually in C++, data scientists now can take advantage of several open-source ML frameworks simplifying the authoring and deployment of NN models. TensorFlow [1] and PyTorch [102] are considered the most popular of such frameworks. ML frameworks follow a common architecture: at the top, they have a high-level Python API 1 where data is commonly represented as multi-dimensional arrays called tensors, while computation is expressed as a composition of tensor operations embedded into the Python language. At the lower level, they have a runtime and a dispatcher/compiler allowing to run the tensor operations over different hardware backends such as CPU, GPU, custom ASICs, and using single node execution, distributed [86], or mobile/edge [61].
Modern ML frameworks allow running computation in an interpreted mode (often referred to as eager execution), or in a compiled mode (graph execution), enabling code optimizations such as common sub-expression elimination, operator fusion, code generation [18], and removing Python dependency [137,138]. Interpreted vs. compiled execution is a popular dichotomy in query processing system implementations [75]. ML frameworks allow both modalities and we explore the trade-offs involved when using one vs. another, and the current limits of tensor compilers in §6.
We will refer to ML frameworks, runtimes [2, 23], and compilers as tensor computation runtimes (TCRs) in the rest of the paper.

Tensor Operations
TCRs provide hundreds of tensor operations. We provide a brief summary of the operators used in TQP, organized by category 2 .
Creation. This category contains all operations used to create tensors, e.g., from_numpy, fill a tensor with specific elements (zeros, ones, empty, fill, arange) or create a tensor using the same shape of another tensor (zeros_like, ones_like).
Indexing and slicing. This category involves operations for selecting one or more elements of a tensor using the square bracket notation, or using indexing (index_select), a mask (masked_select), or a range (narrow).
Reorganization. This category includes reshape, view, and squeeze that reorganize the shape of a tensor (eventually by changing only its metadata). gather, scatter reorganize the elements of a tensor using an index, while sort sorts its elements.
Comparison. eq, lt, gt, le, ge, isnan are operators in this category. Other operations are where that implements conditional statements, and bucketize that implements binary search.
Arithmetic. add, mul, div, sub, fmod, remainder are in this category. We also include logical operators such as logical_and, logical_or, negative, and shift operations.
Join. This category allows to concat or stack multiple tensors.

QUERY PROCESSING ON TCRS
In this section, we summarize the challenges ( §3.2) and the design principles we commit to ( §3.3) when building TQP. First, we show how relational operators can be implemented using tensor programs with an example ( §3.1).

Relational Operators as Tensor Programs
TCRs operate over data represented as tensors. Tensors are arrays of arbitrary dimensions containing elements of the same data type. 0d-tensors are referred to as scalars, 1d-tensors as vectors, and 2dtensors as matrices. For a tensor of dimensions, its shape is a -tuple where each element ∈ {0, 1, . . . , } specifies the size of the -dimension. For example, a matrix with 10 rows and 5 columns is a 2d-tensor of shape (10, 5). This paper only considers dense tensors where each element is explicitly stored in memory. ML practitioners implement programs (NNs) as a composition of operations over tensors. While relational operations are commonly expressed as queries in a standalone language (e.g., SQL), tensor operations are embedded in a host language (e.g., Python), which is used to implement control flows and etc. Next, we introduce examples of implementing a filter using tensors.
Let us assume that we want to implement a simple filter condition over the l_qantity column of the lineitem table: where l_qantity < 24. First, we can represent l_qantity as a 1d-tensor of floating point numbers. We can then use the lt (less than) operator to implement the filter condition (line 1 of Listing 1). lt generates a boolean mask which is then used as a parameter of the masked_select operator to generate the filtered version of the l_qantity column vector (line 2 of Listing 1).
Listing 1: Filter implementation using bitmaps. This implementation is almost identical to the Bitmap-based representation [101] of filters in vectorized query processors [110,118]. On CPU, TCRs have SIMD implementations for several condition and intersection operators. An alternative is to use indexes rather than masks. This is commonly referred to as Selection Vector representation [101,122], and can be similarly implemented using tensor operators lt, nonzero, and index_select.
Listing 2 shows another implementation. Here, we iterate over all the elements of the input tensor and use a Python conditional statement. This implementation does not take advantage of any tensor operation beyond creating the output tensor.
Listing 2: Filter implementation using Python control flow.  Table 1 shows the performance of the two implementations. The implementation using Python control flow is considerably slower. , and GPU execution of Python control flow is slower than CPU execution. This result highlights one of the design choices ( §3.3) we make in TQP: avoid the use of data-dependent code in Python.

Challenges
Implementing a query processor on TCRs requires overcoming several challenges. After all, TCRs are built for authoring and executing NN models, not relational queries.
C1: Expressivity. Relational queries can contain filters with fairly complex expressions (e.g., like, in), sub-queries, group-by aggregates, joins (e.g., natural, anti, semi, outer), etc. It is not clear whether the tensor operations currently available in TCRs are enough to support all these relational operators.
C2: Performance. Even if a relational operator is implementable using tensors, this does not automatically lead to good performance, as the example in Listing 2 suggests. In fact, it is not clear whether tensor programs can achieve good performance, beyond NNs.
C3: Data Representation. To use TCRs as execution engines, relational tables must be transformed into a tensor representation. Previous approaches have explored this challenge (e.g., [66]), but their cost of translation is not negligible. Furthermore, TCRs commonly do not support strings or date data types.
C4: Extensibility. Running relational queries over TCRs makes running a query seamlessly over different hardware (CPU, GPU, ASICs, etc.) and backends (single node, distributed, edge, web browser, etc.) possible. A single monolithic compiler architecture does not work in all situations, therefore TQP's design must be flexible enough to address all these use cases.

Design Choices
When building TQP, we embrace the following design choices.
DC1: Avoid implementing data-dependent control flow in Python. As Table 1 suggests, computation in TQP must use tensor operations as much as possible. Note that for loops and conditionals over schema elements are acceptable (e.g., loops over the columns of a table). This design choice allows us to address C2 and achieve G1.
DC2: Tensor-based columnar format for input tabular data. Relational data must be transformed into the tensor format. To do this, TQP adopts a columnar representation of tables, and considers each column in a table as a tensor. We provide more details on our data representation in §4.1. This design choice addresses C3. DC3: Adherence to TCRs' API. This design choice is required for achieving G2 and G3. In fact, if we start extending TCRs with new features and operators, eventually the system will hinter portability and increase the engineering effort because we will have to support them on any hardware. Hence, we take advantage of existing TCRs' API rather than try to extend them. With this design choice, we are also able to address C1. DC4: Extensible infrastructure allowing easy integration with relational and ML frameworks. Having a flexible infrastructure is of paramount importance since we desire to ride the wave of investments in ML. Therefore, we embrace an extensible architecture that allows different output target formats (e.g., PyTorch, ONNX), composed of a core compiler, pluggable frontends (e.g., query parser and optimizer). This design choice addresses C4.

TENSOR QUERY PROCESSOR (TQP)
In TQP, relational operators and ML models are compiled into tensor programs using a unified infrastructure, extended from Hummingbird [95,97]. Here, we focus on the relational operator part, as the ML part was described in [97].   Figure 1: TQP represents input tables in a columnar format with a 2d-tensor per column.
TQP Overview. TQP's workflow has two phases: (1) compilation: an input query is transformed into an executable tensor program; (2) execution: input data is first transformed into tensors, and then fed into the compiled program to generate the query result. Currently, TQP uses vanilla PyTorch in the compilation phase as the implementation target for the tensor programs. If necessary, PyTorch programs are lowered into different target formats for portability or performance goals. The selection of the hardware device to target is generally made in the compilation phase. Next, we first describe how TQP represents relational data using tensors ( §4.1), and then describe each phase in detail ( §4.2 and §4.3).

Data Representation
Before executing the query, TQP must convert the input (tabular) data to tensors. Databases often manage and convert data into their own proprietary format, and TQP is no different. TQP internally represents tabular data in a columnar format with virtual IDs [29], as illustrated in Figure 1. Data for each column is stored as a ( × ) tensor, where is the input number of rows, and is the length required to store the values. The translation logic is different depending on the column data type. For example, numerical columns (sid in Figure 1) can be directly represented as ( × 1) tensors. The conversion of numerical columns to tensors is often zero-copy. TQP represents date data in ( × 1) numeric tensors as the number of nanoseconds since some pre-defined epoch. In this case, (de)serialization may be required depending on the source/target date representation. Finally, TQP represents string columns using ( × ) numeric tensors, where is the maximum character length of any string for that column. Given a string, TQP stores a character per tensor column and right-pads it with 0s if its length is smaller than . We are actively working on adding support for encoded data (e.g., bit packing, run-length encoding, dictionary encoding) and more compact string representations [16].

Query Compilation
TQP's compilation phase is composed of four main layers, as shown in Figure 2  mapped into a tensor program implementation. (4) The Execution Layer ( §4.2.5), using the operator plan, generates an executor which is the program that runs on the target TCR and hardware. Next, before describing each layer in more detail, we give a quick overview of TQP's intermediate representation (IR).

Intermediate Representation (IR).
The IR is a graph-based data structure. It consists of a list of operators and variables. Each operator corresponds to a node in the graph, and it contains: (1) a list of input variables; (2) a list of output variables; (3) an alias identifying the operator type; and (4) a reference to the corresponding operator instance in the original physical plan. The latter is used to instantiate the tensor program implementing the operator. For example, to create a filter, TQP needs to access the expressions contained in the original physical operator. Edges represent data (tensors) flowing between operators. In particular, an edge connects an output variable from an operator to an input variable of another operator. A variable contains: (1) a unique identifier, and (2) the corresponding frontend column name in the original plan, which is used to translate expressions. When a variable is created, a unique identifier is generated deterministically based on information available in the graph. Variables in the IR are generated as follows. First, TQP generates a variable for each column in the input table. Then, these variables can be used as input to many operators; however, a new variable will always be created for an output of an operator. Thanks to this design: (1) properties (e.g., sorting information) can be immutably attached to columns; (2) the IR is easier to debug because variables, once defined, are never changed; and (3) TQP can detect at runtime when a column is not used anymore and safely garbage-collect it.

Parsing
Layer. The goal of the Parsing Layer is to translate input queries into TQP's internal IR. This goal is accomplished in two steps: (1) input queries are parsed, optimized, and exposed as frontend-specific physical query plans; and (2) a frontend-specific parsing logic translates the physical plan into an IR plan.
In its current version, TQP supports queries expressed as Spark SQL statements, and it uses the PySpark API to parse, optimize, and return the physical plan in a JSON format. We plan to add support for Calcite [39], DuckDB [117], and eventually Substrait [26] 3 . Then the Spark parser constructs the internal IR version of the physical plan using a DFS post-order traversal. If an unsupported operator is found in the plan, this phase will fail with an exception. The list of operators supported by the IR is extensible (DC4).

Canonicalization and Optimization
Layer. This layer implements IR graph transformations similarly to a classical rule-based optimizer. Rules are applied to the IR graph in two stages. In the first stage, canonicalization, the rules are used to eliminate any of the frontend-system idiosyncrasies in the IR graph. For example, Apache Spark returns a projection operator with no inputs for count * statements. In the second stage, optimization, rules rewrite the IR graph for obtaining better performance. While we did not explore in depth the optimization space enabled by TQP's design, we show that hand-optimized tensor programs are more efficient than the one currently generated by TQP in §6.6.

Planning Layer.
In this layer, TQP transforms the optimized IR graph into an operator plan composed of PyTorch tensor programs implementing each operator in the IR graph. In §5, we describe some operator implementations in detail. The implementation of the Planning Layer is straightforward. For each operator in the IR graph, TQP fetches the corresponding implementation containing the tensor program from a dictionary, which is then instantiated with the IR operator's reference to the frontend physical operator instance.

Execution Layer.
Here the operator plan is wrapped around a PyTorch executor object. This object is responsible for: (1) calling the tensor programs in the operator plan following a topological order; (2) wiring the output tensors generated by each program into the successive one; and (3) keeping track of tensor references to garbage collect them if not used anymore. Once the executor program is generated, TQP provides options to compile it into different target formats in addition to PyTorch interpreted execution. Currently, TQP allows lowering the executor into the TorchScript and ONNX formats, as well as to use TVM to compile it directly into machine-level code. Note that not all queries can be compiled into all formats since not all tensor operations are supported by all the target formats.

Execution
Once the executor program is generated, it can be executed over the input data. The program automatically manages (1) converting data into the tensor format; (2) data movements to/from device memory; and (3) scheduling of the operators in the selected device.
Once the data is in the proper format and on the desired device, all the operators are executed sequentially. Regarding parallelization, TQP exploits the tensor-level intra-operator parallelism provided by the TCRs. However, given the poor scalability performance ( §6.3),

OPERATOR IMPLEMENTATION IN TQP
We described how TQP uses the Planning Layer to translate relational operators in the IR graph into tensor programs. Here we provide an overview of a few program implementations. TQP provides tensor-based implementations for the following relational operators: selection, projection, sort, group-by aggregation (sortbased), natural join (hash-based and sort-based), non-equi, leftouter, left-semi, and left-anti joins. TQP supports expressions including comparison and arithmetic operations, functions on date data type, in, case, like statements, as well as aggregate expressions using sum, avg, min, max, and count aggregates (with and without distinct). Finally, TQP supports nulls, and subqueries (scalar, nested, and correlated), and predict UDF 4 [93,94]. With all the above, TQP is able to compile and execute all 22 queries of the TPC-H benchmark (C1). Interestingly, to support the full TPC-H benchmark, only the tensor operations listed in §2.2 are required, and we did not have to introduce any additional custom tensor operators (DC3). Due to space constraints, we only describe how TQP implements relational expressions with tensor operations ( §5.1), and implementations for two representative operators: join (sort-and hash-based, in §5.2 and §5.3, respectively), and group-by aggregation ( §5.4). Finally, note that the filter implementation in TQP is close to the Bitmap representation described in §3.1.

Expressions
Relational expressions such as sum(l_extendedprice * (1 -l_discount)) can be found in projection operators, filter conditions, etc. In an expression tree, each leaf node represents a column or a constant value (e.g., l_extendedprice) and each branch node represents an operator (e.g., * ). TQP keeps an internal dictionary that maps operators to their corresponding tensor operations, e.g., * to torch.mul. To implement an expression with tensor operations, TQP then performs a post-order DFS traversal on the expression tree. For each leaf node, TQP fetches (or generates) the proper column-tensor (constant value). For each internal operator, TQP retrieves the corresponding tensor operation (or a series of tensor operations) from the internal dictionary. In this way (and with the help of Python lambda functions), TQP generates a chain of tensor operations representing the evaluation of the expressions. As an example, from Q21 in TPC-H, the expressions o_orderstatus = 'f' and receiptdate > l_commitdate is implemented as torch.logical_and(torch.eq(o_orderstatus, [70]) ,torch.gt(l_receiptdate,l_commitdate)), where [70] is a 1x1 tensor storing the ASCII value for the constant 'F'.

Sort-Based Join
TQP adopts a late materialization strategy for joins, similar to the one commonly used in columnar databases [30,87]. TQP takes only the columns in the join predicate as input to the join, and the output is a set of pairs of indexes identifying the records for which the join Algorithm 1 Sort-Based Join Input: data: input columns passed as an array of tensors. Output: an array of tensors representing the join output. 1: left, right ← getJoinKeyColumns(data) ⊲ Sort join keys 2: left, leftIdx ← sort(left) 3: right, rightIdx ← sort(right) ⊲ Build histograms for the left and right key columns 4: leftHist, rightHist ← bincount(left), bincount(right) ⊲ Compute the number of rows for each pair of matching keys 5: histMul ← mul(leftHist, rightHist) ⊲ Compute the prefix sums of histograms 6: cumLeftHist ← cumsum(leftHist, dim = 0) 7: cumRightHist ← cumsum(rightHist, dim = 0) 8: cumHistMul ← cumsum(histMul, dim = 0) ⊲ Initialize the output size and output offsets 9: outSize ← cumHistMul [−1] 10: offset ← arange(outSize) ⊲ Find the bucket of matching keys to which each output belongs 11  predicate succeeds. The sort-based equi-join algorithm is shown in Algorithm 1, where, to simplify the description, we describe the case in which two integer columns are joined. With a few modifications, the algorithm is also able to support non-equi joins, left-semi joins, and outer joins. We use the typewriter font (e.g., bucketize) to denote tensor operations, and the capital font (e.g., createOutput) to denote class methods. Figure 3 further illustrates the algorithm. First, TQP sorts the join-key columns from each table (lines 1 to 3 in Algorithm 1, ➊ in Figure 3). Then, ➋, TQP builds two histograms for the join keys from left and right, respectively, i.e., TQP counts the number of occurrences for each unique join key (line 4). Then, ➌ by multiplying the values (element-wise) of the histograms (line 5), TQP computes the bucket sizes: the number of output rows for each matching join key from left and right. Afterward, TQP computes the prefix sums for the left and right histograms (➍), as well as their element-wise multiplication (➎) (lines 6 to 8). The prefix sums will be used later to retrieve, from each join output, the position in left and right. The total size of the output of the join is then computed as the last element of the prefix sum containing the bucket sizes (line 9), and ➏ TQP generates an index array (offset) of the same size (line 10). Then, ➐ TQP performs a parallel binary search on the prefix sum containing the bucket sizes to find the matching join key (bucket) to which each row in the output of the join belongs (line 11). Next, ➑ TQP computes the indexes from left and right that generate each row in the output of the join. Figure 3 shows the computation process for row 8 in the join output of the example. To compute the indexes from left and right that are part of a given offset in the output of the join, TQP first subtracts offset by the prefix sum of bucket sizes prior to the current bucket (line 12). Now offset becomes the offset in each bucket of the matching join keys. TQP then adds to the offset the previous bucket from the respective prefix sum histogram (cumLeftHist and cumRightHist, respectively), and adds the result (quotient for leftOutIdx, remainder for rightOutIdx) of offset divided by the number of join keys from right in the current bucket of matching join keys (lines 13 to 14). Finally, for each row in the join output, TQP knows which rows from left and right contributed to it. It then generates the join output (line 15, not depicted in Figure 3). It is important to note that all computations in this join implementation are achieved using tensor operations, with only minimal usage of Python code.

Hash-Based Join
The hash equi-join algorithm is shown in Algorithm 2. The definition of the input and output here is the same as in §5.2. The algorithm is similar to the classical hash join algorithm, except that the build and probe phases are interleaved and repeated as many times as the maximum number of elements that share a hash value (line 6). The algorithm is as follows: TQP first generates the indexes (line 2) and the hash values (line 3) for the left and right tables. Afterward, TQP computes a histogram over the table on which the hash table will be built (left in this case, line 4) and checks the maximum number of elements in a hash bucket (line 5). Then, TQP repeatedly builds a hash table (lines 7 and 8) and probes it (lines 11 to 14) to find matching keys (lines 15 to 17). Matching keys are accumulated across iterations (lines 18 and 19). In each iteration, TQP also keeps track of the indexes that are stored in the hash table such that they will not appear in subsequent iterations (lines 9 and 10). To achieve this, let be the hash table size; TQP appends an additional (m + 1)-th bucket to the hash table and uses it to redirect the already scattered indexes. Note that when there are no hash collisions, TQP skips the logic of lines 9 to 10 and 18 to 19. This path is therefore close to the optimal.
Compared to the sort-based join, when there are no hash collisions, this implementation is around 30% to 50% faster on CPU and 2× faster on GPU. When there are hash collisions, it is faster than the sort-based join for cases in which at most around 15 elements share a hash value; when there are more than 15 elements

Algorithm 2 Hash-Based Join
Input: data: input columns passed as an array of tensors. Output: an array of tensors representing the join output. sharing a hash value, the sort-based join is faster. We are currently working on a partitioned hash-join implementation.

Aggregation
Algorithm 3 shows the pseudocode of the aggregation implementation. First, TQP horizontally concatenates the values of the group-by columns (lines 1 and 2). TQP then sorts the values of the concatenated columns using radix sort and permutes all the input data columns according to this sorted order (lines 3 and 4). Using uniqueConsecutive, TQP eliminates all but the first key from every consecutive group of equivalent keys. Concurrently, TQP computes the inverted indexes that indicate which bucket (unique key) each row in the sorted list ends up in (line 5). Finally, with the unique key list and inverted indexes, TQP evaluates the aggregate expression for all groups. This last operation makes use of the expression generated (at compile time) as described in §5.1.
Baseline systems. Our goal is to compare TQP with state-ofthe-art query processing systems for different hardware settings. Specifically, for CPU execution, we compare TQP with Apache Spark [145] (recall that Spark and TQP share the same query plans) and DuckDB [117]: a state-of-the-art vectorized engine. Experimental setup. We use the TPC-H benchmark [49] which consists of 22 queries. We use the parameters specified in the query validation sections in [49]. We generate data at different scale factors (from 1 to 10 where 1 means 1 GB of data in total 5 ) using the dbgen tool. We load the generated data from disk into Pandas dataframes. All dataframes use the data types as specified in the benchmark, except for decimals: we use doubles for all systems since TQP does not support decimals yet. Subsequently, we register/convert each dataframe into each system's internal format, e.g., Spark dataframes for Spark 6 , PyTorch tensors for TQP, CUDA dataframes for BlazingSQL, etc., and move the data to the GPU, when applicable. We measure the total query execution time, including the time for generating the output. For each experiment, we do 10 runs where the first 5 are for warm-up. The reported numbers are median values of the last 5 runs.

Key takeaways.
(1) TQP's query execution time on CPU using a single core is better than Spark's over the same physical plans; however, (2) TQP's scalability on CPU is poor because of PyTorch lacking parallelization in some operators' implementation and its intra-operator parallelism model. (3) TQP is, in general, slower than DuckDB on CPU, but for a few queries, TQP is comparable or even better. (4) Hand-optimized plans can improve TQP's performance, which suggests that a TCR-aware query optimizer is required to achieve the best performance. (5) TQP's query execution time on GPU is usually better than both BlazingSQL's and OmnisciDB's, and TQP supports more queries than they do. (6) When ML

Single Core Execution on CPU
In this first experiment, we use a single CPU core and TPC-H at scale factor 1. The results are shown in Table 2 (under CPU). We compare Spark and DuckDB vs. TQP, using both interpreted (TQP) and compiled execution with TorchScript (TQPJ). Spark, DuckDB, and TQP can support all 22 queries. In terms of query time, TQPJ is either comparable to TQP or better. This is because TorchScript removes Python code dependency and provides optimizations not offered by vanilla PyTorch [52]. TQP outperforms Spark for most queries, sometimes by an order of magnitude (e.g., Q10, Q15, and Q22). Given that TQP uses the same physical plans as Spark, this suggests that the tensor abstraction is indeed good for executing relational queries. The practical reasons are: (1) TQP is column-oriented, while Spark is row-oriented. This makes the former better suited for analytical queries; (2) some tensor operations use SIMD instructions, while Spark does not exploit vectorization; (3) in TQP, tensor operations are implemented in C++, while Spark is Java-based; (4) Spark is designed as a scale-out system. For queries (i.e., Q1, Q13, and Q21) where TQP is slower than Spark, the reasons are: (1) TQP's left antijoin and left outer-join implementations are not optimized; (2) the performance of the uniqueConsecutive operator in PyTorch is not optimal. Finally, TQP has better performance than DuckDB only for 3 queries. For the other queries, DuckDB clearly outperforms TQP. If we exclude Q1, Q13, and Q21 (discussed above), TQP's query times are within the same order of magnitude as DuckDB's. To evaluate whether this poor performance compared with DuckDB is due to bad query plans or the tensor abstraction, we hand-code better query plans and tensor programs in §6.6 and show that TQP can match and even outperform DuckDB on CPU.

Execution on GPU
In this experiment, we evaluate the performance of TQP on GPU. The results are shown in Finally, if we compare the best CPU performance versus the best GPU ones, in general, we see that the query times on GPU are 1.5× to 48× better than the CPU ones (single core), except for Q16 where DuckDB is about 3× faster than the best-performing GPU system. This somehow counter-intuitive result is due to the fact that, at scale factor 1, GPU resources are not completely saturated. Therefore, it makes sense to explore how these systems scale with more data and more available core. This is what we explore next.

Scalability
For this and the following experiments, we select a representative set of queries: complex aggregation (Q1), joins and filters (Q2), simple filters (Q6), complex joins (Q9), simple join and aggregation (Q14), a complex mix of join, aggregation, and sub-queries (Q18).

Scaling the Number of Cores.
In this experiment, we scale the number of available CPU cores from 1 to 6 over TPC-H at scale factor 1. Figure 4a compares the scaling performance of Spark, DuckDB, and TQP. Spark has the best scalability trend lines almost for all queries. DuckDB also scales well. TQP's scaling performance is, however sub-optimal, and for some queries increasing the number of cores provides no benefits. There are two reasons: (1) PyTorch uses intra-operator parallelism , which is not as efficient as the shuffle [145] or morsel-based [85] approaches in Spark and DuckDB, respectively; (2) some PyTorch operators run on a single core (e.g., unique and unique_consecutive [116] used in aggregation). We are investigating how to overcome this limitation by adding data-parallel support to TQP leveraging PyTorch Distributed Data Parallel [24,86] or by adding parallel operator implementations.
6.3.2 Scaling the Data. In this experiment, we scale the dataset from 1 GB to 10 GB. In Figure 4b, we compare the scalability   performance of CPU implementations running over 6 cores (Spark, DuckDB), as well as GPU systems (BlazingSQL and OmnisciDB). In general, we see that TQP CPU scales the worst for almost all queries (only Spark is worst for Q6 and Q14), while GPU systems scale better than the CPU ones. For Q1, OmnisciDB provides the best performance, followed by TQP GPU. For Q2, Q14, and Q18, TQP GPU has the best performance, while for Q6, TQP GPU is comparable to OmnisciDB. Finally, for Q9, OmnisciDB has the best performance. Q9 has six joins, and OmnisciDB is able to better use the GPU resources. This query is memory-bound, and the memory bandwidth of the P100 makes it much faster on GPU than on CPU.

Cost/Performance Trade-off
We now provide a cost/performance analysis of TQP on GPU compared to a CPU-only baseline. Specifically, we select a generalpurpose (CPU-only) VM in Azure with a dollar cost similar to the cheapest VM equipped with GPU (NC4as_T4_v3), and with similar main memory size. Following these constraints, we select a D2ds_v5 with 8 CPU cores and 32GB of memory. Then we compare the performance of DuckDB on the D2ds_v5 with TQP on (1) NC4as_T4_v3 (with an NVIDIA T4 GPU, about 15% more expensive than the CPU-only machine), (2) NC6s_v2 (with an NVIDIA P100, around 4.6× more expensive than the CPU-only VM), and (3) NC6s_v3 (with an NVIDIA V100, around 6.6× more expensive than the CPU-only VM). For each GPU VM type, we show the query time speedup required to be more cost-effective than the DuckDB baseline. That is, for the T4, the speedup provided by TQP has to be more than 15% to justify the cost increase of the T4 VM compared to the DuckDB CPU baseline, 4.6× for the P100, 6.6× for the V100. The results for scale factor 10 are shown in Figure 5 for a few representative TPC-H queries. As shown, TQP on GPU is more cost-effective compared to DuckDB on the CPU-only machine: for 6 of the 6 selected queries (17 of the 21 supported queries 7 in the full TPC-H) for the T4; 5 of 6 (10 of 21 in the full TPC-H) for the P100; and 5 of 6 (9 of 21 in the full TPC-H) for the V100.

Performance Breakdown
In this experiment, we show the major contributing factors to the query execution time. TQP is integrated with TensorBoard [13], which provides performance breakdowns and makes it easy to spot bottlenecks [36]. We start by looking into which tensor operators are responsible for the majority of the execution time. Figures 6a  and 6b show the breakdown for a few selected queries on CPU and GPU, respectively. Interestingly, even if TQP uses the same algorithms on both CPU and GPU, the same query can show different operator contributions. For example, for Q1 on CPU, most of the time is spent on computing the unique elements, while on GPU, most is spent on scatter_add. This is because the quality of the operator implementations is different for CPU and GPU. Across queries, on CPU and GPU, the majority of time is also spent on different operators. On CPU, most queries are bounded by unique operators, masked_select, and indexing; on GPU, most of the time is spent on sorting, unique and nonzero. These observations suggest that: (1) the quality of kernels differs between CPU and GPU, e.g., after further investigation, we find that the GPU implementation of scatter_add is not optimal, and nonzero requires host/device synchronization [27] (however, we believe that over time the community will fix such performance issues); and (2) it might be worth investigating backend-aware tensor algorithms.
Finally, we report the GPU utilization for the same set of queries in Figure 7. As we can see, each query has different utilization characteristics. For instance, Q1 contains complex aggregation, and it spends 87% of the time on kernel execution; conversely, Q6 and Q14 are simple queries, and most of the time is spent allocating GPU memory. Finally, Q2 spends a considerable amount of time in generating the output on CPU.

Hand-Optimized Plans
Next, we study whether TQP's performance can be improved with a better optimizer able to generate better tensor programs. To understand this, we hand-optimize the tensor programs for a few selected queries similarly to what a reasonable optimizer with knowledge about cardinalities and tensor characteristics would do, e.g., avoid sorting (or computing unique) over already sorted (or unique) columns, and select better join implementations. The results are shown in Table 3, where we report the best baseline for each setting (CPU 1 and 6 cores, and GPU), and over three execution modes: interpreted PyTorch (Torch), compiled TorchScript (JIT), and compiled using TVM. TVM only supports Q6 and Q14.
If we focus on the CPU numbers first, TQP's performance is comparable to or even better than that of DuckDB's, while TQP was much slower compared to DuckDB both on single-and multi-core execution when not using the hand-optimized plans. TQP is now faster than DuckDB for all queries over 1 CPU core, and two queries over 6 CPU cores. For some queries, TQP is faster than DuckDB by a large margin, e.g., for Q6, 1-core TVM execution is 6× faster. This is because TVM uses code generation and operator fusion to minimize intermediate data materialization across operators. When scaling to 6 cores, TQP scales well only for Q14, while DuckDB scales linearly. For the other queries, TQP's query times improve by at most 2×. This again shows the limitations of PyTorch's scalability on CPU, which cannot be improved by using better tensor programs.
Finally, on GPU, we see that OmnisciDB has still better performance for Q9, although TQP's query time for Q9 on GPU improves by 4×, when using the hand-optimized plans. This is because TQP's aggregate implementation heavily uses sorting, while OmnisciDB uses hash-based implementations.

Prediction Queries
We now investigate the performance benefits of using a unified runtime for queries mixing relational and ML operators. We use prediction queries as a use case, i.e., queries embedding a trained ML model performing predictions over some input data [94]. Recall that TQP natively supports predictions of any PyTorch model (e.g., NNs), and traditional ML models through its integration with Hummingbird. Here, we join the customer and orders tables in TPC-H (scale factor 10), and train a gradient boosting tree model (with 128 trees with max depths of 8) over a mix of categorical (c_orderstatus) and numerical features (c_custkey, c_nationkey, c_acctbal, sum(o_totalprice)) after we apply one-hot encoding and feature scaling, respectively. We run a prediction query using the trained model over the query with two filter predicates added (c_mktsegment = 'building' and o_orderdate >= date '1993-10-01'). Note that this prediction query mixes ML operators (tree ensemble, one-hot encoding, scaling, and concatenation) with relational ones (join, aggregation and filtering). We compare TQP with two baselines: one where the prediction query is executed over Spark (MLlib [90] is used to build the model), and one where we use DuckDB for the relational part and scikit-learn [106] for the ML part 8 . Since TQP subsumes Hummingbird, it is able to compile both the ML and the relational operators of the query into a unified plan executable on TCRs. Figure 8 shows the result. For CPU single core, TQP is about 40% faster than Spark, while DuckDB+scikit-learn is about 7× faster than TQP. When enabling all cores, Spark and DuckDB scale much better than TQP, for the reasons described in §6.3. Finally, TQP is able to exploit GPU acceleration end-to-end, which brings a 9× improvement of query time compared to the best CPU baseline.

Overheads
Next, we evaluate the overheads of TQP for both CPU and GPU. The breakdown of the end-to-end execution with all overheads is shown in Figure 9. Note that: (1) data conversion is done once and many databases (e.g., BlazingSQL, OmnisciDB, Spark, SQL Server, etc.) requires it; (2) TQP pipelines data movement (to the GPU) with query execution (non-blocking IO), while for this experiment we explicitly make data movement blocking; (3) the machine in this experiment uses PCIe 3 which is 4× slower than the latest version, Table 3: Query execution time (in seconds) on selected TPC-H queries (scale factor 10). TQP Hand-Opt. uses hand-optimized tensor programs. We use Torch, JIT, and TVM to refer to execution using PyTorch (interpreted), TorchScript (compiled), and TVM, respectively. Bold numbers highlight the best performance for the specific setup: CPU (1 core), CPU (6 cores), or GPU.  PCIe 5; (4) query compilation can be cached, but here we report the full query compilation time as the sum of the time for the frontend database to generate the physical plan, and the time for TQP to generate the final executable tensor program.
If we focus first on the CPU side (Figure 9a), compilation and data conversion take the majority of the time only for simple queries (e.g., Q6), while for the other queries, the majority of the time is spent on the query execution. However, in the GPU case (Figure 9b), except for Q2 and Q9, the majority of the time is spent on data operations (conversion and movement) and compilation. However, in practice, as described above, these overheads are hidden (e.g., data movement using pipelining) or are one-time overheads (data conversion and query compilation). Regarding query compilation, 90% of the time is spent initializing the PyTorch models from the Spark plans, and we are currently investigating how to speed up this process. Finally, using TorchScript adds substantial compilation overheads since queries are traced using input samples.

Portability
To evaluate whether TQP can run on different hardware and software backends, we run TPC-H Query 6 with the hand-optimized plan on: (1) two integrated graphic cards, one from Intel, and one from AMD; (2) two discrete GPUs from NVIDIA (K80 and V100: the former a generation before the P100 GPU used for the experiments in the previous sections; the latter one, one generation after); (3) a custom ASIC used for NN training and inference (TPU); and (4) a web browser. We use a scale factor of 1. The results are shown in Table 4. This experiment proves the versatility of TQP. For the integrated GPUs, we use TVM to code-generate the query using Metal [35]. For the two discrete GPUs, we use vanilla PyTorch, while for the TPU, we use the XLA backend for PyTorch 9 [114]. Finally, we are able to run the query in the browser by exporting it into the ONNX format and running it in Chrome using ONNX Runtime (ORT) for WebAssembly (WASM) [96].

Engineering Effort
To demonstrate the minimal engineering effort required by TQP to run queries over different hardware, we compare the lines of code for a few relational operators (hash and sort-based joins, aggregation) across all evaluated systems. For each relational operator and each system, we use cloc [51] to count the lines of source code (excluding comment and blank lines) from the files containing the algorithmic functionality of the operator. This is admittedly a subjective process, but we believe the numbers of lines of code can roughly reflect the engineering effort required to implement relational operators in each system. Table 5 shows the results. Compared with the baselines, TQP requires significantly lower engineering effort: up to 10× less compared to CPU implementations, and 50× less compared to GPU ones. It is worth noting that TQP is able to target different hardware with the same implementation, so the engineering effort required for TQP to scale over different hardware is constant. The other baseline systems do not share this property. For instance, to run Spark on GPU (e.g., using RAPIDS [12], the same backend of BlazingSQL), we would have to add the lines of code for the GPU implementation.

RELATED WORK
Common representation for relational and ML workloads. Since the '90s [98], there have been many works trying to integrate relational queries with data science and ML workloads [15, 32, 41, 42, 45, 50, 55, 64, 67, 68, 73, 74, 79, 82, 91, 93, 107, 112, 123-125, 128, 133, 141, 143]. To our knowledge, we are the first to propose executing relational queries over TCRs. Earlier attempts tried to run a few relational operators on the TPU using TensorFlow [65]. TQP is orthogonal to previous efforts to optimize relational and tensor algebra (e.g., [67,141]), and we believe TQP can leverage them to improve its performance further. An analysis of matrix query languages can be found in [58]. Here, we focus on TCRs' tensor interface, which is more flexible than a linear algebra API.
SciDB [119,132] is a database using arrays as the base data representation. TensorDB [77] further proposes support for tensor data and decomposition operations inside databases. SciDB, TensorDB, and TQP suggest using a format closer to data science and ML to represent data. However, TQP further exploits TCRs to run both relational and ML workloads on hardware accelerators.
GPUs and hardware accelerators. Several systems have explored running relational queries over GPUs [84,88,103,104,111,127,144]. We refer readers to [105] for a recent survey. However, (a) End-to-end execution breakdown on CPU (b) End-to-end execution breakdown on GPU Figure 9: End-to-end breakdown (incl. all overheads, and w/o pipelining and caching) for selected queries at scale factor 10. Table 4: Query time (in milliseconds) of TPC-H Query 6 (scale factor 1) using the hand-optimized plan over different hardware and software backends. In parenthesis is the TCR used as well as the compilation stack (when applicable).  the majority of them focus mostly on microbenchmarks, while, to our knowledge, only RateUpDB can support the full TPC-H benchmark. TQP is able to run the TPC-H benchmark on both CPU and GPU, thanks to TCRs' flexibility to support different hardware backends. TCUDB [66] suggests using the Tensor Core Unit (TCU) of GPUs for accelerating relational operators. TCUDB requires an expensive transformation from tables to matrices and also uses low-level CUDA kernels, while TQP takes advantage of the highlevel tensor interface of TCRs. GPUs are the default hardware for running neural network models. However, there has recently been a rise in custom ASICs [4, 6,11,34,72] purposely built for ML workloads. With TQP, we propose a solution allowing us to run relational queries on any hardware supported by TCRs, since many ASICs [5,10,72] provide high-level interfaces directly through TCRs or are targetable through tensor compilers [46,83].
Query processing over heterogeneous hardware. Several recent works have started to explore query execution over heterogeneous hardware, such as CPU-GPU co-execution [44,47,57,63,108,120,121,140]. Many of them rely on OpenCL [8] to target different hardware. However, targeting a common language (or similarly a generic compiler, e.g., MLIR [83]), requires non-trivial engineering effort since each device requires proper tuning [108], algorithms, and data structures (as well as abstractions/dialects in the MLIR case). Conversely, TQP can natively run on any hardware supported by TCRs, and uses TCRs' tensor operation implementations and compilation stacks. Currently, the user has to specify which fragment of the query should run on which hardware, but we are exploring how to automate this and enable co-execution.
A trend arises recently that suggests splitting relational operators into smaller functions that can be easily composed and efficiently dispatched over heterogeneous hardware [38,80,139]. TQP fits in this trend, whereby tensor operations are sub-components.
Vectorized execution, query compilation, and columnar databases. MonetDB/X100 [43] pioneered the vectorized execution model as well as the columnar data layout [131]. TQP follows a similar design, where data is stored in a columnar format with virtual IDs [30], but each column is represented as a tensor. Recent works, such as HyPer [99] and others [92,100,126], have focused on query compilation. Nevertheless, since (1) there is no clear winner between query compilation and vectorized execution [75] ; (2) many industry-grade systems use vectorized execution because it is easier to debug and profile [40]; and (3) compiled systems start to move to vectorized execution (e.g., Spark with Photon), we evaluate TQP against a state-of-the-art vectorized engine, DuckDB [117].
On the ML systems side, TensorFlow initially embraced a compiled (graph) execution [31], while PyTorch pioneered interpreted (eager) execution [102]. Compilers [14,28,46,53,54,78,83] and optimization techniques [69][70][71] for neural networks are hot topics in the MLSys community. With TQP, we aim to ride the wave of innovation in this domain. For TQP, interpreted vs. compiled execution is just another point in the query optimization space, since TCRs allow to switch between them seamlessly.

CONCLUSION
We proposed TQP, the first system able to run relational queries on TCRs. TQP is able to take advantage of all the innovation poured into TCRs, as well as to run efficiently on any hardware devices supported by TCRs. Our experiments showed not only that TQP is capable of running the full TPC-H benchmark on TCRs, but also that TQP's performance is comparable and often superior to that of specialized CPU and GPU query processing systems.