# Flexible Modeling and Multitask Learning using Differentiable Tree Ensembles

Shibal Ibrahim MIT Department of Electrical Engineering and Computer Science.
email<EMAIL_ADDRESS>Hussein Hazimeh MIT Operations Research Center. email:
<EMAIL_ADDRESS>Rahul Mazumder MIT Sloan School of Management, Operations
Research Center and MIT Center for Statistics. email<EMAIL_ADDRESS>
(August, 2021)

###### Abstract

Decision tree ensembles are widely used and competitive learning models.
Despite their success, popular toolkits for learning tree ensembles have
limited modeling capabilities. For instance, these toolkits support a limited
number of loss functions and are restricted to single task learning. We
propose a flexible framework for learning tree ensembles, which goes beyond
existing toolkits to support arbitrary loss functions, missing responses, and
multi-task learning. Our framework builds on differentiable (a.k.a. soft) tree
ensembles, which can be trained using first-order methods. However, unlike
classical trees, differentiable trees are difficult to scale. We therefore
propose a novel tensor-based formulation of differentiable trees that allows
for efficient vectorization on GPUs. We perform experiments on a collection of
28 real open-source and proprietary datasets, which demonstrate that our
framework can lead to 100x more compact and 23% more expressive tree ensembles
than those by popular toolkits.

##  1 Introduction

Decision tree ensembles are popular models that have proven successful in
various machine learning applications and competitions (Erdman and Bates,
2016; Chen and Guestrin, 2016). Besides their competitive performance,
decision trees are appealing in practice because of their interpretability,
robustness to outliers, and ease of tuning (Hastie et al., 2009). Training a
decision tree naturally requires solving a combinatorial optimization problem,
which is difficult to scale. In practice, greedy heuristics are commonly used
to get feasible solutions to the combinatorial problem; for example CART
(Breiman et al., 1984), C5.0 (Quinlan, 1993), and OC1 (Murthy et al., 1994).
By building on these heuristics, highly scalable toolkits for learning tree
ensembles have been developed, e.g., XGBoost (Chen and Guestrin, 2016) and
LightGBM (Ke et al., 2017). These toolkits are considered a defacto standard
for training tree ensembles and have demonstrated success in various domains.

Despite their success, popular toolkits for learning tree ensembles lack
modeling flexibility. For example, these toolkits support a limited set of
loss functions, which may not be suitable for the application at hand.
Moreover, these toolkits are limited to single task learning. In many modern
applications, it is desired to solve multiple, related machine learning tasks.
In such applications, multi-task learning, i.e., learning tasks
simultaneously, may be a more appropriate choice than single task learning
(Chapelle et al., 2010; Kumar and Daumé, 2012; Han and Zhang, 2015; Crawshaw,
2020). If the tasks are sufficiently related, multi-task learning can boost
predictive performance by leveraging task relationships during training.

In this paper, we propose a flexible modeling framework for training tree
ensembles that addresses the aforementioned limitations. Specifically, our
framework allows for training tree ensembles with any differentiable loss
function, enabling the user to seamlessly experiment with different loss
functions and select what is suitable for the application. Moreover, our
framework equips tree ensembles with the ability to perform multi-task
learning. To achieve this flexibility, we build up on soft trees (Kontschieder
et al., 2015; Hazimeh et al., 2020), which are differentiable trees that can
be trained with first-order (stochastic) gradient methods. Previously, soft
tree ensembles have been predominantly explored for classification tasks with
cross-entropy loss. In such cases, they were found to lead to more expressive
and compact tree ensembles (Hazimeh et al., 2020). However, the state-of-the-
art toolkits e.g., TEL (Hazimeh et al., 2020), are slow as they only support
CPU training and are difficult to customize. Our proposed tree ensemble
learning framework supports a collection of loss functions such as
classification, regression, Poisson regression, zero-inflation models,
overdispersed distributions, multi-task learning and has seamless support for
others — the user can modify the loss function with a single line. We show
loss customization can lead to a significant reduction in ensemble sizes (up
to 20x). We also propose a careful tensor-based formulation of differentiable
tree ensembles, which leads to more efficient training (10x) on CPUs as well
as support GPU-training (20x).

We propose a novel extension for multi-task learning with tree ensembles —
this may not be readily accommodated within popular gradient boosting
toolkits. Our model can impose both ’soft’ information sharing across
tasks—soft information sharing is not possible with popular multi-task tree
ensemble toolkits e.g., RF Breiman (2001), GRF Athey et al. (2019). Our
proposed framework leads to $23\%$ performance gain (on average) and up to a
100x reduction in tree ensemble sizes on real-world datasets, offering a
promising alternative approach for learning tree ensembles.

#### Contributions

Our contributions can be summarized as follows. (i) We propose a flexible
framework for training differentiable tree ensembles with seamless support for
new loss functions. (ii) We introduce a novel, tensor-based formulation for
differentiable tree ensembles that allows for efficient training on GPUs.
Existing toolkits e.g., TEL Hazimeh et al. (2020), only support CPU training.
(iii) We extend differentiable tree ensembles to multi-task learning settings
by introducing a new regularizer that allows for soft parameter sharing across
tasks. (iv) We introduce a new toolkit (based on Tensorflow 2.0) and perform
experiments on a collection of 28 open-source and real-world datasets,
demonstrating that our framework can lead to 100x more compact ensembles and
up to $23\%$ improvement in out-of-sample performance, compared to tree
ensembles learnt by popular toolkits such as XGBoost Chen and Guestrin (2016).

##  2 Related Work

Learning binary trees has been traditionally done in three ways. The first
approach relies on greedy construction and/or optimization via methods such as
CART (Breiman et al., 1984), C5.0 (Quinlan, 1993), OC1 (Murthy et al., 1994),
TAO (Carreira-Perpinan and Tavallali, 2018). These methods optimize a
criterion at the split nodes based on the samples routed to each of the nodes.
The second approach considers probabilistic relaxations/decisions at the split
nodes and performs end-to-end learning with first order methods (Irsoy et al.,
2012; Frosst and Hinton, 2017; Lay et al., 2018). The third approach considers
optimal trees with mixed integer formulations and jointly optimize over all
discrete/continuous parameters with MIP solvers (Bennett, 1992; Bennett and
Blue, 1996; Bertsimas and Dunn, 2017; Zhu et al., 2020). Each of the three
approaches have their pros and cons. The first approach is highly scalable
because of greedy heuristics. In many cases, the tree construction uses a
splitting criterion different from the optimization objective (Breiman et al.,
1984) (e.g., gini criterion when performing classification) possibly resulting
in sub-optimal performance. The second approach is also scalable but
principled pruning in probabilistic trees remains an open research problem.
The third approach scales to small datasets with samples $N\sim 10^{4}$,
features $p\sim 10$ and tree depths $d\sim 4$.

Jointly optimizing over an ensemble of classical decision trees is a hard
combinatorial optimization problem (Hyafil and Rivest, 1976). Historically,
tree ensembles have been trained with two methods. The first method relies on
greedy heuristics with bagging/boosting: where individual trees are trained
with CART on bootstrapped samples of the data e.g., random forests (RF)
(Breiman, 2001) and its variants (Geurts et al., 2006; Athey et al., 2019); or
sequentially/adaptively trained with gradient boosting: Gradient Boosting
Decision Trees (Hastie et al., 2009) and efficient variants (Chen and
Guestrin, 2016; Ke et al., 2017; Prokhorenkova et al., 2018; Schapire and
Freund, 2012). Despite the success of ensemble methods, interesting challenges
remain: (i) RF tend to under-perform gradient boosting methods such as XGBoost
(Chen and Guestrin, 2016). (ii) The tree ensembles are typically very large,
making them a complex and hardly interpretable decision structure. Recent work
by Carreira-Perpinan and Tavallali (2018); Zharmagambetov and Carreira-
Perpiñán (2020) improve RF with local search methods via alternating
minimization. However, their implementation is not open-source. (iii) Open-
source APIs for gradient boosting are limited in terms of flexibility. They
lack support for multi-task learning, handling missing responses or catering
to customized loss functions. Modifying these APIs require significant effort
and research for custom applications.

The alternative approach for tree ensemble learning extends
probabilistic/differentiable trees and performs end-to-end learning
(Kontschieder et al., 2015; Hazimeh et al., 2020). These works build upon the
idea of hierarchical mixture of experts introduced by (Jordan and Jacobs,
1994) and further developed by (Irsoy et al., 2012; Tanno et al., 2019; Frosst
and Hinton, 2017) for greedy construction of trees. Some of these works
(Kontschieder et al., 2015; Hazimeh et al., 2020) propose using differentiable
trees as an output layer in a cascaded neural network for combining feature
representation learning along with tree ensemble learning for classification.
We focus on learning tree (ensembles) with hyperplane splits and constant leaf
nodes–this allows us to expand the scope of trees to flexible loss functions,
and develop specialized implementations that can be more efficient. One might
argue that probabilistic trees are harder to interpret and suffer from slower
inference as a sample must follow each root-leaf path, lacking _conditional
computation_ present in classical decision trees. However, Hazimeh et al.
(2020) proposed a principled way to get conditional inference in probabilistic
trees by introducing a new activation function; this allows for routing
samples through small parts of the tree similar to classical decision trees.
We refer the reader to Section 6.1 for a study on a single tree and highlight
that the a soft tree with hyperplane splits and conditional inference has
similar interpretability as that of a classical tree with hyperplane splits —
see Figure 2. Additionally, a soft tree can lead to smaller optimal depths—see
Supplemental Section S1.1.

End-to-end learning with differentiable tree ensembles appears to have several
advantages. (i) They are easy to setup with public deep learning API e.g.,
Tensorflow (Abadi et al., 2015), PyTorch (Paszke et al., 2019). We demonstrate
that with a careful implementation, the tree ensembles can perform efficient
GPU-training — this is not possible with earlier toolkits e.g., TEL (Hazimeh
et al., 2020). (ii) Tensorflow-Probability library (Dillon et al., 2017)
offers huge flexibility in modeling. For instance, they allow for modeling
mixture likelihoods, which can cater to zero-inflated data applications. This
makes it possible to combine differentiable tree ensembles with customized
losses with minimal effort. (iii) Deep learning APIs e.g., Tensorflow offer
great support for handling multi-task loss objectives without a need for
specialized algorithm. This makes it convenient to perform multi-task learning
with differentiable tree ensembles. (iv) Differentiable trees can lead to more
expressive and compact ensembles (Hazimeh et al., 2020). This can have
important implications for interpretability, latency and storage requirements
during inference.

##  3 Optimizing Tree Ensembles

We assume a supervised multi-task learning setting, with input space
$\mathcal{X}\subseteq\mathbb{R}^{p}$ and output space
$\mathcal{Y}\subseteq\mathbb{R}^{k}$. We learn a mapping
$f:\mathbb{R}^{p}\rightarrow\mathbb{R}^{k}$, from input space $\mathcal{X}$ to
output space $\mathcal{Y}$, where we parameterize function $f$ with a
differentiable tree ensemble. We consider a general optimization framework
where the learning objective is to minimize any differentiable function
$g:\mathbb{R}^{p}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$. The framework can
accommodate different loss functions arising in different applications and
perform end-to-end learning with tree ensembles. We first briefly review
differentiable tree ensembles in section 3.1. We later present a careful
tensor-based implementation of the tree ensemble in section 3.2 that makes it
possible to perform efficient training on both CPUs and GPUs. Next, we outline
the need for flexible loss modeling and present some examples that our tree
ensemble learning framework supports for learning compact tree ensembles in
Sections 4 and 5. Finally, we present a metastudy on a large collection of
real-world datasets to validate our framework.

#### Notation

We summarize our notation used in the paper. For an integer $n\geq 1$, let
$[n]:=\\{1,2,....,n\\}$. We let $1_{m}$ denote the vector in $\mathbb{R}^{m}$
with all coordinates being $1$. For matrix
$\bm{B}=((B_{ij}))\in\mathbb{R}^{M\times N}$, let the $j$-th column be denoted
by $\bm{B}_{j}:=[B_{1j},B_{2j},...,B_{Mj}]^{T}\in\mathbb{R}^{M}$ for
$j\in[N]$. A dot product between two vectors $\bm{u},\bm{v}\in R^{m}$ is
denoted as $\bm{u}\cdot\bm{v}$. A dot product between a matrix $\bm{U}\in
R^{m,n}$ and a vector $\bm{v}\in\mathbb{R}^{m}$ is denoted as
$\bm{U}\cdot\bm{v}=\bm{U}^{T}v\in\mathbb{R}^{n}$. A dot product between a
tensor $\bm{\mathcal{U}}\in\mathbb{R}^{p,m,n}$ and a vector
$\bm{v}\in\mathbb{R}^{m}$ is denoted as
$\bm{\mathcal{U}}\cdot\bm{v}=\bm{\mathcal{U}}^{T}\bm{v}\in\mathbb{R}^{p,n}$
where the transpose operation of a tensor
$\bm{\mathcal{U}}^{T}\in\mathbb{R}^{p,n,m}$ permutes the last two dimensions
of the tensor.

###  3.1 Preliminaries and Setup

We learn an ensemble of $m$ differentiable trees. Let $\bm{f}^{j}$ be the
$j$th tree in the ensemble. For easier exposition, we consider a single-task
regression or classification setting—see Section 5 for an extension to the
multi-task setting. In a regression setting $k=1$, while in multi-class
classification setting $k=C$, where $C$ is the number of classes. For an input
feature-vector $\bm{x}\in\mathbb{R}^{p}$, we learn an additive model with the
output being sum over outputs of all the trees:

$\displaystyle\bm{f}(\bm{x})=\sum_{j=1}^{m}\bm{f}^{j}(\bm{x}).$ (1)

The output, $\bm{f}(\bm{x})$, is a vector in $\mathbb{R}^{k}$ containing raw
predictions. For multiclass classification, mapping from raw predictions to
$\mathcal{Y}$ is done by applying a softmax function on the vector
$\bm{f}(\bm{x})$ and returning the class with the highest probability. Next,
we introduce the key building block of the approach: the differentiable
decision tree.

#### Differentiable decision trees for modelling $\bm{f}^{j}$

Classical decision trees perform hard sample routing, i.e., a sample is routed
to exactly one child at every splitting node. Hard sample routing introduces
discontinuities in the loss function, making trees unamenable to continuous
optimization. Therefore, trees are usually built in a greedy fashion. In this
section, we first introduce a single soft tree proposed by Jordan and Jacobs
(1994), which is utilized in Irsoy et al. (2012); Frosst and Hinton (2017);
Blanquero et al. (2021) and extended to soft tree ensembles in Kontschieder et
al. (2015); Hehn et al. (2019); Hazimeh et al. (2020). A soft tree is a
variant of a decision tree that performs soft routing, where every internal
node can route the sample to the left and right simultaneously, with different
proportions. This routing mechanism makes soft trees differentiable, so
learning can be done using gradient-based methods. Notably, Hazimeh et al.
(2020) introduced a new activation function for soft trees that allowed for
conditional computation while preserving differentiability.

Let us fix some $j\in[m]$ and consider a single tree $\bm{f}^{j}$ in the
additive model (1). Recall that $\bm{f}^{j}$ takes an input sample and returns
an output vector (logit), i.e.,
$\bm{f}^{j}:X\in\mathbb{R}^{p}\rightarrow\mathbb{R}^{k}$. Moreover, we assume
that $\bm{f}^{j}$ is a perfect binary tree with depth $d$. We use the sets
$\mathcal{I}^{j}$ and $\mathcal{L}^{j}$ to denote the internal (split) nodes
and the leaves of the tree, respectively. For any node
$i\in\mathcal{I}^{j}\cup\mathcal{L}^{j}$, we define $A^{j}(i)$ as its set of
ancestors and use the notation $\bm{x}\rightarrow i$ for the event that a
sample $\bm{x}\in\mathbb{R}^{p}$ reaches $i$.

#### Routing

Internal (split) nodes in a differentiable tree perform soft routing, where a
sample is routed left and right with different proportions. This soft routing
can be viewed as a probabilistic model. Although the sample routing is
formulated with a probabilistic model, the final prediction of the tree
$\bm{f}$ is a deterministic function as it assumes an expectation over the
leaf predictions. Classical decision trees are modeled with either axis-
aligned splits (Breiman et al., 1984; Quinlan, 1993) or hyperplane (a.k.a.
oblique) splits (Murthy et al., 1994). Soft trees are based on hyperplane
splits, where the routing decisions rely on a linear combination of the
features. Particularly, each internal node $i\in\mathcal{I}^{j}$ is associated
with a trainable weight vector $\bm{w}_{i}^{j}\in\mathbb{R}^{p}$ that defines
the node’s hyperplane split. Let $S:R\rightarrow[0,1]$ be an activation
function. Given a sample $\bm{x}\in\mathbb{R}^{p}$, the probability that
internal node $i$ routes $\bm{x}$ to the left is defined by
$S(\bm{w}_{i}^{j}\cdot\bm{x})$. Now we discuss how to model the probability
that $\bm{x}$ reaches a certain leaf $l$. Let
$[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]$ (resp.
[$i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l$]) denote the event that leaf
$l$ belongs to the left (resp. right) subtree of node $i\in\mathcal{I}^{j}$.
Assuming that the routing decision made at each internal node in the tree is
independent of the other nodes, the probability that $\bm{x}$ reaches $l$ is
given by:

$\displaystyle P^{j}(\\{x\rightarrow l\\})=\prod_{i\in
A(l)}r_{i,l}^{j}(\bm{x}),$ (2)

where $r_{i,l}^{j}(\bm{x})$ is the probability of node $i$ routing $\bm{x}$
towards the subtree containing leaf $l$, i.e.,
$r_{i,l}^{j}(x):=S(\bm{w}_{i}^{j}\cdot\bm{x})1[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]\odot(1-S(\bm{w}_{i}^{j}\cdot\bm{x}))1[i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l]$. Next, we define how the root-
to-leaf probabilities in (2) can be used to make the final prediction of the
tree.

#### Prediction

As with classical decision trees, we assume that each leaf stores a weight
vector $\bm{o}_{l}^{j}\in R^{k}$ (learned during training). Note that, during
a forward pass, $\bm{o}_{l}^{j}$ is a constant vector, meaning that it is not
a function of the input sample(s). For a sample $\bm{x}\in\mathbb{R}^{p}$, we
define the prediction of the tree as the expected value of the leaf outputs,
i.e.,

$\displaystyle\bm{f}^{j}(\bm{x})=\sum_{l\in L}P^{j}(\\{\bm{x}\rightarrow
l\\})\bm{o}_{l}^{j}.$ (3)

#### Activation Function

In soft routing, the internal nodes use an activation function $S$ in order to
compute the routing probabilities. The logistic (a.k.a. sigmoid) function is
the common choice for $S$ in the literature on soft trees (see Jordan and
Jacobs (1994); Kontschieder et al. (2015); Frosst and Hinton (2017); Tanno et
al. (2019); Hehn et al. (2019)). While the logistic function can output
arbitrarily small values, it cannot output an exact zero. This implies that
any sample $x$ will reach every node in the tree with a positive probability
(as evident from (2)). Thus, computing the output of the tree in (3) will
require computation over every node in the tree, an operation which is
exponential in tree depth. Hazimeh et al. (2020) proposed a smooth-step
activation function, which can output exact zeros and ones, thus allowing for
true conditional computation. Furthermore, the smooth-step function allows
efficient forward and backward propagation algorithms.

###  3.2 Efficient Tensor Formulation

Current differentiable tree ensemble proposals and toolkits, for example deep
neural decision
forests111https://keras.io/examples/structured_data/deep_neural_decision_forests/
(Kontschieder et al., 2015) and TEL (Hazimeh et al., 2020) model trees
individually. This leads to slow CPU-training times and makes these
implementations hard to vectorize for fast GPU training. In fact, TEL (Hazimeh
et al., 2020) doesn’t support GPU training. We propose a tensor-based
formulation of a tree ensemble that parallelizes routing decisions in nodes
across the trees in the ensemble. This can lead to 10x faster CPU training
times if the ensemble sizes are large e.g., $100$. Additionally, the tensor-
based formulation is GPU-friendly, which provides an additional $40\%$ faster
training times. See Figure 1 for a timing comparison on CPU training
without/with tensor formulation. Next, we outline the tensor-based
formulation.

Figure 1: Timing comparison on CPU training with tensor-based implementation.
GPU training with tensor-based implementation leads to an additional 40%
improvement over CPU training with tensor-based implementation.

We propose to model the internal nodes in the trees across the ensemble
jointly as a “supernodes”. In particular, an internal node
$i\in\mathcal{I}^{j}$ at depth $d$ in all trees can be condensed together into
a supernode $i\in\mathcal{I}$. We define a learnable weight matrix
$\bm{W}_{i}\in\mathbb{R}^{p,m}$, where each $j$-th column of the weight matrix
contains the learnable weight vector $\bm{w}_{i}^{j}$ of the original j-th
tree in the ensemble. Similarly, the leaf nodes are defined to store a
learnable weight matrix $\bm{O}_{l}\in\mathbb{R}^{m,k}$, where each $j$-th row
contains the learnable weight vector $\bm{o}_{l}^{j}$ in the original $j$-th
tree in the ensemble. The prediction of the tree with supernodes can be
written as

$\displaystyle\small\bm{f}(\bm{x})=\left(\sum_{l\in
L}\bm{O}_{l}\odot\prod_{i\in A(l)}\bm{R}_{i,l}\right)\cdot\bm{1}_{m}$ (4)

where $\odot$ denotes the element-wise product,
$\bm{R}_{i,l}=S(\bm{W}_{i}\cdot\bm{x})1[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]\odot(1-S(\bm{W}_{i}\cdot\bm{x}))1[i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l]\in\mathbb{R}^{m,1}$ and the
activation function $S$ is applied element-wise. This formulation of tree
ensembles via supernodes allows for sharing of information across tasks via
tensor formulation in multi-task learning — see Section 5 for more details.

###  3.3 Toolkit

Our tree ensemble learning toolkit is built in Tensorflow (TF) 2.0 and
integrates with Tensorflow-Probability. The toolkit allows the user to write a
custom loss function, and TF provides automatic differentiation. Popular
packages, such as XGBoost, require users to provide first/second order
derivatives. In addition to writing a custom loss, the user can select from a
wide range of predefined loss and likelihood functions from Tensorflow-
Probability. By relying on TF in the backend, our toolkit can easily exploit
distributed computing. It can also run on multiple CPUs or GPUs, and on
multiple platforms, including mobile platforms. Checkpointing for fault
tolerance, incremental training and warm restarts are also supported. The
toolkit will be open sourced if the paper is accepted.

##  4 Flexible loss functions

Our framework can handle any differentiable loss function. Such flexibility is
important as various applications require flexibility in loss functions beyond
what is provided by current tree ensemble learning APIs. Our framework is
built on Tensorflow, which gives us the capability to perform gradient-based
optimization on scale. This coupled with our efficient differentiable tree
ensemble formulation gives a powerful toolkit to seamlessly experiment with
different loss functions and select what is suitable for the intended
application. A few examples of flexible distributions that our toolkit
supports — due to compatibility with Tensorflow-Probability — are normal,
Poisson, gamma, exponential, mixture distributions e.g., zero-inflation models
(Lee, 2020), and compound distributions e.g., negative binomial (Yirga et al.,
2020). Other loss functions such as robustness to outliers (Barron, 2019) can
also be handled by our tree ensemble learning toolkit. To demonstrate the
flexibility of our framework, we deeply investigate two specific examples:
zero-inflated Poisson and negative binomial regression. These cannot be
handled by the popular gradient boosting APIs such as XGBoost (Chen and
Guestrin, 2016), LightGBM (Ke et al., 2017).

#### Zero-inflated Poisson Regression

Zero-inflation occurs in many applications, e.g., understanding alcohol and
drug abuse in young adults Jacobs et al. (2021), characterizing undercoverage
and overcoverage to gauge the on-going quality of the census frames (Young et
al., 2017), studying popularity of news items on different social media
platforms Moniz and Torgo (2018), financial services applications Lee (2020)
etc. Despite the prevalence of these applications, there has been limited work
on building decision tree-based approaches for zero-inflated data perhaps due
to a lack of support public APIs. Therefore, practitioners either resort to
Poisson regression with trees or simpler linear models to handle zero-inflated
responses. A Poisson model can lead to sub-optimal performance due to the
limiting equidispersion constraint (mean equals the variance). Others take a
two-stage approach (Cameron and Trivedi, 2013), where a classification model
distinguishes the zero and non-zero and a second model is used to model the
non-zero responses. This can be sup-optimal as errors in the first model can
deteriorate the performance of the second model. We employ a more well-
grounded approach by formulating the joint mixture model, where one part of
the model tries to learn the mixture proportions (zero vs non-zero) and the
other part models the actual non-zero responses. Such a mixture model permits
a differentiable loss function when both components of the model are
parameterized with differentiable tree ensembles and can be optimized with
gradient descent method in an end-to-end fashion without the need for a custom
solver. We provide an extensive study with our framework on small to large-
scale real world zero-inflated datasets and demonstrate that such flexibility
in distribution modeling can lead to significantly more compact and expressive
tree ensembles. This has large implications for faster inference, storage
requirements and interpretability.

We briefly review Poisson regression and then dive into zero-inflated Poisson
models. Poisson regression stems from the generalized linear model (GLM)
framework for modeling a response variable in the exponential family of
distributions. In general, GLM uses a link function to provide the
relationship between the linear predictors, $\bm{x}$ and the conditional mean
of the density function:

$\displaystyle g[\mathbb{E}({\textnormal{y}}|\bm{x})]=\bm{\beta}\cdot\bm{x},$
(5)

where $\bm{\beta}$ are parameters and $g(\cdot)$ is the link function. When
responses $y_{n}$ (for $n\in[N])$, are independent and identically distributed
(i.i.d.) and follow the Poisson distribution conditional on $\bm{x}_{n}$’s, we
use $log(\cdot)$ as the link function and call the model a Poisson regression
model: $log(\mu_{n}|\bm{x}_{n})=\bm{\beta}\cdot\bm{x}_{n}$. We consider more
general parameterizations with tree ensembles as given by

$\displaystyle
log(\mu_{n}|\bm{x}_{n})=f(\bm{x}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}}).$ (6)

where $f$ is parameterized with a tree ensemble as in (1) and
$\bm{\mathcal{W}},\bm{\mathcal{O}}$ are the learnable parameters in the
supernodes and the leaves of the tree ensemble. When a count data has excess
zeros, the equi-dispersion assumption of the Poisson is violated. The Poisson
model is not an appropriate model for this situation anymore. Lambert (1992)
proposed zero-inflated-Poisson (ZIP) models that address the mixture of excess
zeros and Poisson count process. The mixture is indicated by the latent binary
variable $d_{n}$ using a logit model and the density for the Poisson count
given by the log-linear model. Thus,

$\displaystyle y_{n}=\begin{cases}0,&\text{if}~{}~{}d_{n}=0\\\
y_{n}^{*},&\text{if}~{}~{}d_{n}=1\end{cases}$ (7)

where the latent indicator $d_{n}\sim Bernoulli(\pi_{n})$ with
$\pi_{n}=P(d_{n}=1)$ and $y_{n}^{*}\sim Poisson(\mu_{n})$. The mixture yields
the marginal probability mass function of the observed $y_{n}$ given as:

$\displaystyle
ZIP(y_{n}|\mu_{n},\pi_{n})=\begin{cases}(1-\pi_{n})+\pi_{n}e^{-\mu_{n}},&\text{if}~{}~{}y_{n}=0\\\
\pi_{n}e^{-\mu_{n}}\mu_{n}^{y_{n}}/y_{n}!,&\text{if}~{}~{}y_{n}=1,2,\cdots\end{cases}$
(8)

where $\mu_{n}$ and $\pi_{n}$ are modeled by

$\displaystyle log\left(\frac{\pi_{n}}{1-\pi_{n}}|\bm{x}_{n}\right)$
$\displaystyle=f(\bm{x}_{n};\bm{\mathcal{Z}},\bm{\mathcal{U}})$ (9)
$\displaystyle log(\mu_{n}|\bm{x}_{n})$
$\displaystyle=f(\bm{x}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}}).$ (10)

where $\bm{\mathcal{Z}},\bm{\mathcal{U}}$ are the learnable parameters in the
splitting internal supernodes and the leaves of the tree ensemble for the
logit model for $\pi_{n}$ and $\bm{\mathcal{W}},\bm{\mathcal{O}}$ are the
learnable parameters in the supernodes and the leaves of the tree ensemble for
the log-tree model for $\mu_{n}$ respectively. The likelihood function for
this ZIP model is given by

$\displaystyle
L(y_{n},f(\bm{x}_{n}))=\prod_{y_{n}=0}(1-\pi_{n})+\pi_{n}e^{-\mu_{n}}\prod_{y_{n}>0}\pi_{n}e^{-\mu_{n}}\mu_{n}^{y_{n}}/y_{n}!$
(11)

where $\lambda_{n}=e^{f(\bm{x}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}})}$ and
$\pi_{n}=e^{f(\bm{x}_{n},\bm{\mathcal{Z}};\bm{\mathcal{U}})}/(1+e^{f(\bm{x}_{n};\bm{\mathcal{Z}},\bm{\mathcal{U}})})$.Such
a model can be overparameterized and we observed that sharing the learnable
parameters $\bm{\mathcal{Z}}=\bm{\mathcal{W}}$ in the splitting internal
supernodes across the log-mean and logit models can lead to better out-of-
sample performance — see Section 6 for a thorough evaluation on real-world
datasets.

#### Negative Binomial Regression

An alternative distribution to zero-inflation modeling that can cater to over-
dispersion in the responses is Negative Binomial (NB) distribution. A negative
binomial distribution for a random variable y with a non-negative mean
$\mu\in\mathbb{R}_{+}$ and dispersion parameters $\phi\in\mathbb{R}_{+}$ is
given by:

$\displaystyle
NB({\textnormal{y}}|\mu,\phi)=\left(\frac{{\textnormal{y}}+\phi-1}{{\textnormal{y}}}\right)\left(\frac{\mu}{\mu+\phi}\right)^{{\textnormal{y}}}\left(\frac{\phi}{\mu+\phi}\right)^{\phi}$
(12)

The mean and variance of a random variable ${\textnormal{y}}\sim
NB({\textnormal{y}}|\mu,\phi)$ are $\mathbb{E}[{\textnormal{y}}]=\mu$ and
$\mathrm{Var}[{\textnormal{y}}]=\mu+\mu^{2}/\phi$. Recall that Poisson($\mu$)
has variance $\mu$, so $\mu^{2}/\phi>0$ is the additional variance of the
negative binomial above that of the Poisson with mean $\mu$. So the inverse of
parameter $\phi$ controls the overdispersion, scaled by the square of the
mean, $\mu^{2}$.

When the responses $y_{n}$ (for $n\in[N]$) are i.i.d, and follow NB
distribution conditioned on $\bm{x}_{n}$’s, we can use the $log(.)$ as a link
function to parameterize the log-mean and log-dispersion as linear functions
of the covariates $\bm{x}_{n}$. In our parameterization with Tree Ensembles,
we model them as given by:

$\displaystyle log(\mu_{n}|\bm{x}_{n})$
$\displaystyle=f(\bm{x}_{n};\bm{\mathcal{W}},\bm{\mathcal{O}})$ (13)
$\displaystyle log\left(\phi_{n}|\bm{x}_{n}\right)$
$\displaystyle=f(\bm{x}_{n};\bm{\mathcal{Z}},\bm{\mathcal{U}}).$ (14)

where $\bm{\mathcal{Z}},\bm{\mathcal{U}}$ are the learnable parameters in the
supernodes and the leaves of the tree ensemble for the log-mean and
$\bm{\mathcal{W}},\bm{\mathcal{O}}$ are the learnable parameters in the
supernodes and the leaves of the tree ensemble for the log-dispersion model
for $\phi_{n}$ respectively. Such a model can be overparameterized and we
observed that sharing the learnable parameters
$\bm{\mathcal{Z}}=\bm{\mathcal{W}}$ in the splitting internal supernodes
across the log-mean and log-dispersion models can lead to better out-of-sample
performance. See Section 6.4 for empirical validation on a large-scale
dataset.

##  5 Multi-task Learning with Tree Ensembles

Multi-task Learning (MTL) aims to learn multiple tasks simultaneously by using
a shared model. Unlike single task learning, MTL can achieve better
generalization performance through exploiting task relationships (Caruana,
1997; Chapelle et al., 2010). One key problem in MTL is how to share model
parameters between tasks (Ruder, 2017). For instance, sharing parameters
between unrelated tasks can potentially degrade performance. MTL approaches
for classical decision trees approaches e.g., RF (Linusson, 2013), GRF (Athey
et al., 2019) have shared weights at the splitting nodes across the tasks.
Only the leaf weights are task specific. However this can be limiting in terms
of performance, despite easier interpretability associated with the same split
nodes across tasks.

To perform flexible multi-task learning, we extend our formulation in Section
3.2 by using task-specific nodes in the tree ensemble. We consider $T$ tasks.
For easier exposition, we consider tasks of the same kind: multilabel
classification or multi-task regression. For multilabel classification, each
task is assumed to have same number of classes (with $k=C$) for easier
exposition — our framework can handle multilabel settings with different
number of classes per task. Similarly, for regression settings, $k=1$. For
multi-task zero-inflated Poisson or negative binomial regression, when two
model components need to be estimated, we set $k=2$ to predict log-mean and
logit components for zero-inflated Poisson and log-mean and log-dispersion
components for negative binomial.

We define a trainable weight tensor
$\bm{\mathcal{W}}_{i}\in\mathbb{R}^{T,p,m}$ for supernode $i\in\mathcal{I}$,
where each $t$-th slice of the tensor $\bm{\mathcal{W}}_{i}[t,:,:]$ denotes
the trainable weight matrix associated with task $t$. The prediction in this
case is given by

$\displaystyle\bm{f}(\bm{x})=\biggl{(}\sum_{l\in
L}\bm{\mathcal{O}}_{l}\odot\prod_{i\in
A(l)}\bm{\mathcal{R}}_{i,l}\biggr{)}\cdot\bm{1}_{m}$ (15)

where $\bm{\mathcal{O}}_{l}\in\mathbb{R}^{T,m,k}$ denotes the trainable leaf
tensor in leaf $l$,
$\bm{\mathcal{R}}_{i,l}=S(\bm{\mathcal{W}}_{i}\cdot\bm{x})1[l\mathrel{\text{\rotatebox[origin={c}]{225.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}i]\odot(1-S(\bm{\mathcal{W}}_{i}\cdot\bm{x}))1[i\mathrel{\text{\rotatebox[origin={c}]{315.0}{$\vrule
height=0.0pt,width=0.0pt\shortrightarrow$}}}l]\in\mathbb{R}^{T,m,1}$.

In order to share information across the tasks, our framework imposes
closeness penalty on the hyperplanes $\bm{\mathcal{W}}_{i}$ in the supernodes
across the tasks. This results in an optimization formulation:

$\displaystyle\min_{\bm{\mathcal{W}},\bm{\mathcal{O}}}\sum_{t\in
T}\sum_{\bm{x},y_{t}}g_{t}(y_{t},f_{t}(\bm{x}))+\lambda\sum_{s<t,t\in
T}\left\lVert\bm{\mathcal{W}}_{:,s,:,:}-\bm{\mathcal{W}}_{:,t,:,:}\right\rVert^{2},$
(16)

where $\bm{\mathcal{W}}\in\mathbb{R}^{\mathcal{I},T,m,p}$ denotes all the
weights in all the supernodes,
$\bm{\mathcal{O}}\in\mathbb{R}^{\mathcal{L},m,k}$ denotes all the weights in
the leaves, and $\lambda\in[0,\infty)$ is a non-negative regularization
penalty that controls how close the weights across the tasks are. The model
behaves as single-task learning when $\lambda=0$ and complete sharing of
information in the splitting nodes when $\lambda\rightarrow\infty$. Note that
when $\lambda\rightarrow\infty$, the weights across the tasks in each of the
internal supernodes become the same. This case can be separately handled more
efficiently by using the function definition in (4) for $\bm{f}(\bm{x})$
without any closeness regularization in (16). Our model can control the level
of sharing across the tasks by controlling $\lambda$. In practice, we tune
over $\lambda\in[1e-5,10]$ and select the optimal based on the validation set.
This penalty assumes that the hyperplanes across the tasks should be equally
close as we go down the depth of the trees. However this assumption maybe less
accurate as we go down the tree. Empirically, we found that decaying $\lambda$
exponentially as $\lambda/2^{d}$ with depth $d$ of the supernodes in the
ensemble can achieve better test performance.

##  6 Experiments

We study the performance of differentiable tree ensembles in various settings
and compare against the relevant state-of-the-art baselines for each setting.
The different settings can be summarized as follows: (i) Comparison of a
single soft tree with state-of-the-art classical tree method in terms of test
performance and depth. We include both axis-aligned and oblique classical tree
in our comparisons. (ii) Flexible zero-inflation models with tree ensembles.
We compare against Poisson regression with tree ensembles and gradient
boosting decision trees (GBDT). We consider test Poisson deviance and tree
ensemble compactness for model evaluation (iii) We evaluate our proposed
multi-task tree ensembles and compare them against multioutput RF, multioutput
GRF and single-task GBDT. We consider both fully observed and partially
observed responses across tasks. (iv) We also validate our tree ensemble
methods with flexible loss functions (zero-inflated Poisson and negative
binomial regression) on a large-scale multi-task proprietary dataset.

#### Model Implementation

Differentiable tree ensembles are implemented in TensorFlow 2.0 using Keras
interface.

#### Datasets

We use 27 open-source regression datasets from various domains (e.g., social
media platforms, human behavior, finance). 9 are from Mulan (Xioufis et al.,
2016), 2 are from UCI data repository (Dua and Graff, 2017), 12 are from Delve
database (Akujuobi and Zhang, 2017) and the 5 remaining are SARCOS
(Vijayakumar and Schaal, 2000), Youth Risk Behavior Survey (Jacobs et al.,
2021), Block-Group level Census data (US Census Bureau, 2021), and a financial
services loss dataset from Kaggle222https://bit.ly/3swGnTo. We also validate
our framework on a proprietary multi-task data with millions of samples from a
multi-national financial services company.

###  6.1 Studying a single tree

In this section, we compare performance and model compactness of a single tree
on 12 regression datasets from Delve database: abalone, pumadyn-family, comp-
activ (cpu, cpuSmall) and concrete.

#### Competing Methods and Implementation

We focus on two baselines from classical tree literature: CART (Breiman et
al., 1984) and Tree Alternating Optimization (TAO) method proposed by
Carreira-Perpinan and Tavallali (2018). The authors in Zharmagambetov and
Carreira-Perpiñán (2020) performed an extensive comparison of various single
tree learners and demonstrated TAO to be the best performer. Hence, we include
both axis-aligned and oblique decision tree versions of TAO in our
comparisons. Given that the authors in Carreira-Perpinan and Tavallali (2018);
Zharmagambetov and Carreira-Perpiñán (2020) do not provide an open-source
implementation for TAO, we implemented our own version of TAO. For a fair
comparison, we use binary decision trees for both axis-aligned and oblique
versions of TAO. For more details about our implementation, see Supplemental
Section S1.1.

Table 1: Test mean squared error performance of a single axis aligned and
oblique decision tree on various regression datasets.

|  Axis-Aligned |  Oblique
---|---|---
Data |  CART |  TAO |  TAO |  Soft Tree
abalone |  7.901E-03 |  8.014E-03 |  7.205E-03 |  6.092E-03
pumadyn-32nh |  9.776E-03 |  9.510E-03 |  1.146E-02 |  8.645E-03
pumadyn-32nm |  2.932E-03 |  2.741E-03 |  4.942E-03 |  1.280E-03
pumadyn-32fh |  1.324E-02 |  1.307E-02 |  1.370E-02 |  1.166E-02
pumadyn-32fm |  2.246E-03 |  2.177E-03 |  3.566E-03 |  1.602E-03
pumadyn-8nh |  1.995E-02 |  1.983E-02 |  2.027E-02 |  1.670E-02
pumadyn-8nm |  4.878E-03 |  4.584E-03 |  4.206E-03 |  2.352E-03
pumadyn-8fh |  2.182E-02 |  2.200E-02 |  2.159E-02 |  2.048E-02
pumadyn-8fm |  4.398E-03 |  4.347E-03 |  4.074E-03 |  3.543E-03
cpu |  9.655E-04 |  1.475E-03 |  1.312E-03 |  9.159E-04
cpuSmall |  1.450E-03 |  2.319E-03 |  1.171E-03 |  9.159E-04
concrete |  7.834E-03 |  6.619E-03 |  1.166E-02 |  4.139E-03

Figure 2: Classifier boundaries for CART [Left], TAO (oblique) [Middle] and
Soft tree [Right] on a synthetic dataset with
$N_{train}=N_{val}=N_{test}=2500$ generated using sklearn (Buitinck et al.,
2013). We tune for 50 trials over depths in the range $[2-4]$ for TAO
(oblique) and soft trees and $[2-10]$ for CART. Optimal depths for CART, TAO,
Soft tree are 5, 4 and 2 respectively. Test AUCs are 0.950, 0.957, and 0.994
respectively.

#### Results

We present the out-of-sample mean-squared-error performance and optimal depths
in Tables 1 and S1 (in Supplemental Section S1.1) respectively. Notably, in
all $12$ cases, soft tree outperforms all 3 baseline methods in terms of test
performance. The soft tree finds a smaller optimal depth in majority cases in
comparison with its classical counterpart i.e., oblique TAO tree — See Table
S1 in Supplemental Section S1.1. This may be due to the end-to-end learning in
a soft tree, unlike TAO that performs local search.

###  6.2 Zero-inflation

We consider a collection of real-world applications with zero-inflated data.
The datasets include (i) yrbs: nationwide drug use behaviors of high school
students as a function of demographics, e-cigarettes/mariyuana use etc.; (ii)
news: popularity of news items on social media platforms Moniz and Torgo
(2018) e.g., Facebook, Google+ as a function of topic and sentiments; (iii)
census: number of people with zero, one, or two health insurances across all
Census blocks in the ACS population as a function of housing and socio-
economic demographics; (iv) fin-services-loss: financial services losses as a
function of geodemographics, information on crime rate, weather.

#### Competing methods

We consider Poisson regression with GBDT and differentiable tree ensembles. We
also consider zero-inflation modeling with differentiable tree ensembles. We
use GBDT from sklearn Buitinck et al. (2013). For additional details about the
tuning experiments, please see Supplemental Section S1.2.

Table 2: Flexible modeling via zero-inflated Poisson for Soft Tree Ensembles
leads to better out-of-sample performance.

|  |  | GBDT | Soft Trees
---|---|---|---|---
Data | N | p | Poisson | ZIP
yrbs-cocaine | 12172 | 55 | 3.14E-02 | 3.00E-02 | 2.82E-02
yrbs-heroine | 12711 | 55 | 1.81E-02 | 1.60E-02 | 1.54E-02
yrbs-meth | 12690 | 55 | 2.38E-02 | 2.21E-02 | 2.09E-02
yrbs-lsd | 9564 | 55 | 3.51E-02 | 3.59E-02 | 3.43E-02
news-facebook | 81637 | 3 | 4.70E-03 | 4.68E-03 | 4.68E-03
news-google+ | 87495 | 3 | 5.97E-03 | 5.93E-03 | 5.93E-03
census-health0 | 220333 | 64 | 1.51E-04 | 1.28E-04 | 1.32E-04
census-health1 | 220333 | 64 | 4.63E-04 | 5.11E-04 | 4.66E-04
census-health2+ | 220333 | 64 | 2.73E-03 | 3.06E-03 | 2.72E-03
fin-services-losses | 452061 | 300 | 2.20E-03 | 2.28E-03 | 2.20E-03
#wins | - | - | 2 | 3 | 8

Table 3: Flexible modeling via zero-inflated Poisson for Soft Tree Ensembles
can lead to more compact tree ensembles, which can potentially lead to easier
interpretability

| #Trees | Depth
---|---|---
| GBDT | Soft Trees | GBDT | Soft Trees
Data | Poisson | ZIP | Poisson | ZIP
yrbs-cocaine | 575 | 77 | 13 | 4 | 2 | 3
yrbs-heroine | 1425 | 83 | 4 | 4 | 4 | 2
yrbs-meth | 1475 | 16 | 4 | 4 | 2 | 3
yrbs-lsd | 1225 | 17 | 45 | 4 | 3 | 2
news-facebook | 200 | 65 | 85 | 4 | 4 | 4
news-google+ | 750 | 81 | 74 | 4 | 4 | 4
census-health0 | 1275 | 10 | 10 | 4 | 3 | 2
census-health1 | 1275 | 17 | 17 | 4 | 2 | 2
census-health2+ | 1375 | 73 | 56 | 8 | 4 | 3
fin-services-losses | 1225 | 32 | 4 | 4 | 3 | 2
#wins | 0 | 4 | 7 | - | 5 | 8

#### Results

We present the out-of-sample Poisson deviance performance in Table 2. Notably,
tree ensembles with zero-inflated loss function leads the chart. We also
present the optimal selection of tree ensemble sizes and depths in Table 3. We
can observe that zero-inflation modeling can lead to significant benefits in
terms of model compression. Both tree ensemble sizes and depths can
potentially be made smaller, which have implications for faster inferences,
memory footprint and interpretability.

Table 4: Performance of RF, GRF and multi-task Differentiable Tree Ensembles
on 11 multi-task regression datasets with fully observed responses across
tasks.

|  |  Multi-task
---|---|---
Data |  Task |  RF |  GRF |  Soft Trees
|  1 |  2.242E-02 |  1.847E-02 |  5.383E-03
|  2 |  2.498E-02 |  1.894E-02 |  7.217E-03
|  3 |  1.127E-02 |  9.625E-03 |  1.128E-02
|  4 |  1.574E-02 |  1.504E-02 |  1.403E-02
|  5 |  2.040E-02 |  1.905E-02 |  1.182E-02
atp1d |  6 |  1.571E-02 |  1.333E-02 |  1.527E-02
|  1 |  3.244E-03 |  2.691E-03 |  8.965E-03
|  2 |  3.914E-03 |  3.931E-03 |  4.456E-03
|  3 |  1.231E-02 |  1.078E-02 |  1.059E-02
|  4 |  5.459E-03 |  6.542E-03 |  6.001E-03
|  5 |  1.842E-03 |  1.922E-03 |  3.358E-03
atp7d |  6 |  4.042E-03 |  5.303E-03 |  5.033E-03
|  1 |  4.384E-02 |  3.151E-02 |  3.345E-02
|  2 |  1.077E-02 |  4.638E-03 |  3.828E-03
sf1 |  3 |  4.252E-02 |  2.863E-02 |  2.993E-02
|  1 |  7.883E-03 |  8.807E-03 |  7.789E-03
|  2 |  3.247E-03 |  2.583E-03 |  2.206E-03
sf2 |  3 |  1.262E-03 |  3.744E-03 |  3.595E-03
|  1 |  3.027E-02 |  3.008E-02 |  2.233E-02
|  2 |  1.483E-02 |  1.405E-02 |  1.015E-02
jura |  3 |  7.896E-03 |  7.586E-03 |  6.036E-03
|  1 |  2.473E-04 |  1.865E-04 |  2.063E-04
enb |  2 |  2.830E-03 |  3.305E-03 |  1.054E-03
|  1 |  1.732E-01 |  1.396E-01 |  1.001E-01
|  2 |  1.224E-01 |  9.827E-02 |  7.368E-02
slump |  3 |  3.878E-02 |  2.944E-02 |  5.149E-03
|  1 |  3.040E-03 |  2.530E-03 |  1.794E-03
|  2 |  3.397E-03 |  3.003E-03 |  2.226E-03
|  3 |  4.178E-03 |  3.611E-03 |  2.940E-03
scm1d |  4 |  3.991E-03 |  3.376E-03 |  2.150E-03
|  1 |  4.457E-03 |  3.650E-03 |  2.198E-03
|  2 |  4.766E-03 |  3.632E-03 |  2.410E-03
|  3 |  4.892E-03 |  3.506E-03 |  2.620E-03
scm20d |  4 |  5.573E-03 |  4.072E-03 |  2.632E-03
|  1 |  3.466E-03 |  2.558E-03 |  1.730E-03
|  2 |  4.636E-03 |  4.039E-03 |  3.728E-03
bike |  3 |  5.123E-03 |  4.303E-03 |  3.822E-03
# wins |  - |  5 |  6 |  26

Table 5: Tree ensemble sizes for soft trees, RF, and GRF.

| atp1d | atp7d | sf1 | sf2 | jura | enb | slump | scm1d | scm20d | bike
---|---|---|---|---|---|---|---|---|---|---
RF | 100 | 125 | 500 | 225 | 150 | 100 | 825 | 175 | 100 | 100
GRF | 1050 | 950 | 350 | 100 | 150 | 100 | 350 | 50 | 100 | 250
Ours | 10 | 15 | 12 | 44 | 17 | 13 | 54 | 49 | 25 | 93

###  6.3 Multi-task Regression

We compare performance and model compactness of our proposed regularized
multi-task tree ensembles on 11 multi-task regression datasets from Mulan
(atp1d, atp7d, sf1, sf2, jura, enb, slump, scm1d, scm20d), and UCI data
repository (bike) and SARCOS dataset.

#### Competing Methods

We focus on 4 tree ensemble baselines from literature: single-task soft tree
ensembles, sklearn GBDT, sklearn multioutput RF Buitinck et al. (2013) and
r-grf package for GRF Athey et al. (2019). We consider two multi-task
settings: (i) All Fully observed responses for all tasks, (ii) Partially
observed responses across tasks. In the former case, we compare against RF and
GRF. In the latter case, we compare against single-task soft tree ensembles
and GBDT. Note the open-source implementations for RF and GRF do not support
partially observed responses for multi-task settings and GBDT does not have
support for multi-task setting. We refer the reader to Supplemental Section
S1.3 for tuning experiments details.

Table 6: Performance of GBDT, single-task and multi-task Soft Tree Ensembles
on 11 multi-task regression datasets with 50% missing responses per task.

|  |  Single-Task |  Multi-Task
---|---|---|---
Data |  Task |  GBDT |  Soft Trees
|  1 |  1.469E-02 |  3.091E-02 |  1.295E-02
|  2 |  7.698E-03 |  2.598E-02 |  1.137E-02
|  3 |  2.172E-02 |  2.915E-02 |  1.807E-02
|  4 |  5.905E-03 |  1.417E-02 |  9.434E-03
|  5 |  2.421E-02 |  5.631E-02 |  2.105E-02
atp1d |  6 |  5.646E-03 |  5.880E-02 |  2.724E-02
|  1 |  5.562E-02 |  6.261E-02 |  3.601E-02
|  2 |  4.033E-02 |  2.216E-02 |  1.067E-02
|  3 |  2.140E-02 |  4.989E-02 |  1.680E-02
|  4 |  1.107E-02 |  2.089E-02 |  9.926E-03
|  5 |  4.254E-02 |  6.476E-02 |  3.053E-02
atp7d |  6 |  4.195E-03 |  2.416E-02 |  , 2.199E-02
|  1 |  2.674E-02 |  2.763E-02 |  2.732E-02
|  2 |  5.825E-03 |  1.680E-02 |  5.435E-03
sf1 |  3 |  3.030E-02 |  3.704E-02 |  2.964E-02
|  1 |  1.003E-02 |  1.179E-02 |  1.009E-02
|  2 |  1.123E-02 |  6.843E-03 |  2.809E-03
sf2 |  3 |  9.271E-03 |  1.039E-02 |  8.064E-03
|  1 |  1.779E-02 |  1.934E-02 |  1.503E-02
|  2 |  1.117E-02 |  1.665E-02 |  9.304E-03
jura |  3 |  1.311E-02 |  1.514E-02 |  1.262E-02
|  1 |  1.509E-04 |  2.738E-04 |  4.167E-04
enb |  2 |  1.071E-03 |  1.227E-03 |  1.160E-03
|  1 |  1.622E-01 |  7.611E-02 |  9.485E-02
|  2 |  8.823E-02 |  1.050E-01 |  4.734E-02
slump |  3 |  8.423E-03 |  1.737E-03 |  7.744E-03
|  1 |  1.598E-03 |  1.903E-03 |  2.058E-03
|  2 |  1.952E-03 |  2.703E-03 |  2.490E-03
|  3 |  3.029E-03 |  3.194E-03 |  2.919E-03
scm1d |  4 |  2.666E-03 |  3.656E-03 |  3.272E-03
|  1 |  2.541E-03 |  2.672E-03 |  2.533E-03
|  2 |  3.640E-03 |  3.174E-03 |  3.146E-03
|  3 |  3.658E-03 |  4.015E-03 |  3.201E-03
scm20d |  4 |  3.756E-03 |  4.115E-03 |  3.670E-03
|  1 |  2.122E-03 |  2.558E-03 |  2.300E-03
|  2 |  3.680E-03 |  4.038E-03 |  3.846E-03
bike |  3 |  3.731E-03 |  4.303E-03 |  3.910E-03
|  1 |  2.582E-04 |  1.317E-04 |  1.518E-04
|  2 |  1.643E-04 |  8.310E-05 |  8.307E-05
sarcos |  3 |  3.325E-04 |  1.788E-04 |  1.933E-04
# wins |  - |  14 |  5 |  23

#### Results

We present results for fully observed response settings in Table 4 and
partially observed response settings in Table 6. In both cases, regularized
multi-task soft trees lead the charts over the corresponding baselines in
terms of out-of-sample mean squared error performance. For the fully observed
response setting, we also show tree ensemble sizes in Table 5. We see a large
reduction in the number of trees with out proposed multi-task tree ensembles.

Table 7: Out-of-sample performance of single-task and multi-task tree
ensembles with flexible loss functions for zero-inflation/overdispersion. We
evaluate performance with weighted Poisson deviance and AUC across tasks.

|  | GBDT | Soft Trees
---|---|---|---
|  | Single-task | Single-task | Multi-task
Metric | Task | Poisson | Poisson | ZIP | NB | Poisson | ZIP | NB
Poisson Deviance | 1 | 2.643E-04 | 2.624E-04 | 2.623E-04 | 2.623E-04 | 2.607E-04 | 2.605E-04 | 2.608E-04
2 | 8.029E-04 | 8.050E-04 | 8.029E-04 | 8.044E-04 | 8.022E-04 | 8.014E-04 | 8.014E-04
3 | 1.044E-03 | 1.045E-03 | 1.043E-03 | 1.042E-03 | 1.041E-03 | 1.040E-03 | 1.041E-03
AUC | 1 | 0.710 | 0.721 | 0.722 | 0.721 | 0.730 | 0.734 | 0.727
2 | 0.690 | 0.689 | 0.690 | 0.688 | 0.691 | 0.691 | 0.692
3 | 0.684 | 0.683 | 0.686 | 0.685 | 0.687 | 0.689 | 0.689

###  6.4 Large-scale multi-task data from a multinational financial services
company

We study the performance of our differentiable tree ensembles in a real-word,
large-scale multi-task setting from a multinational financial services
company. The system encompasses costs and fees for millions of users for
different products and services. The dataset has the following
characteristics: (i) It is a multi-task regression dataset with 3 tasks. (ii)
Each task has high degree of over-dispersion. (iii) All tasks are not fully
observed as each user signs up for a subset of products/services. The degree
of missing responses on average across tasks is $\sim 50\%$. (iv) Number of
features is also large ($\sim 600$).

We validate the flexibility of our end-to-end tree-ensemble learning framework
with soft trees on a dataset of 1.3 million samples. We study the following
flexible aspects of our framework: (i) Flexible loss handling with zero-
inflation Poisson regression and negative binomial regression for single-task
learning. (ii) Multi-task learning with our proposed regularized multi-task
soft tree ensembles in the presence of missing responses across tasks. (iii)
Flexible loss handling with zero-inflation Poisson/negative binomial
regression in the context of multi-task learning.

We present our results in Table 7. We can see that we achieve the lowest
Poisson deviance and highest AUC with multi-task regression via zero-inflated
Poisson/negative binomial regression.

##  7 Conclusion

We propose a flexible and scalable tree ensemble learning framework with
seamless support for new loss functions. Our framework makes it possible to
compare different loss functions or incorporate new differentiable loss
functions and train them with tree ensembles. Our proposal for tensor-based
modeling of tree ensembles allows 10x faster training on CPU than existing
toolkits and also allow for training differentiable tree ensembles on GPUs. We
also propose a novel extension for multi-task learning with trees that allow
soft sharing of information across tasks for more compact and expressive
ensembles than those by existing tree ensemble toolkits.

### Acknowledgements

We thank Denis Sai for his help with the experiments on large-scale multi-task
data. This research was supported in part, by grants from the Office of Naval
Research: ONR-N00014-21-1-2841 and award from Liberty Mutual Insurance.

## References

  * Erdman and Bates [2016] Chandra Erdman and Nancy Bates.  The low response score (LRS).  _Public Opinion Quarterly_ , 81(1):144–156, December 2016.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450342322.
  * Hastie et al. [2009] T. J. Hastie, R. J. Tibshirani, and J. Friedman.  _The Elements of Statistical Learning_.  Springer, 2 edition, 2009.
  * Breiman et al. [1984] L. Breiman, J. Friedman, C.J. Stone, and R.A. Olshen.  _Classification and Regression Trees_.  Taylor & Francis, 1984.  ISBN 9780412048418.
  * Quinlan [1993] J. Ross Quinlan.  _C4.5: Programs for Machine Learning_.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1993.  ISBN 1558602380.
  * Murthy et al. [1994] Sreerama K. Murthy, Simon Kasif, and Steven Salzberg.  A system for induction of oblique decision trees.  _J. Artif. Int. Res._ , 2(1):1–32, aug 1994.  ISSN 1076-9757.
  * Ke et al. [2017] Guolin Ke, Qi Meng, Thomas Finley, et al.  Lightgbm: A highly efficient gradient boosting decision tree.  In I. Guyon, U. V. Luxburg, S. Bengio, et al., editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Chapelle et al. [2010] Olivier Chapelle, Pannagadatta K. Shivaswamy, Srinivas Vadrevu, Kilian Q. Weinberger, Ya Zhang, and Belle L. Tseng.  Boosted multi-task learning.  _Machine Learning_ , 85:149–173, 2010.
  * Kumar and Daumé [2012] Abhishek Kumar and Hal Daumé.  Learning task grouping and overlap in multi-task learning.  In _Proceedings of the 29th International Coference on International Conference on Machine Learning_ , ICML’12, pages 1723–1730, Madison, WI, USA, 2012. Omnipress.  ISBN 9781450312851.
  * Han and Zhang [2015] Lei Han and Yu Zhang.  Learning tree structure in multi-task learning.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’15, pages 397–406, New York, NY, USA, 2015. Association for Computing Machinery.  ISBN 9781450336642.
  * Crawshaw [2020] Michael Crawshaw.  Multi-task learning with deep neural networks: A survey.  _ArXiv_ , abs/2009.09796, 2020.
  * Kontschieder et al. [2015] Peter Kontschieder, Madalina Fiterau, Antonio Criminisi, et al.  Deep neural decision forests.  In _2015 IEEE International Conference on Computer Vision (ICCV)_ , pages 1467–1475, 2015.
  * Hazimeh et al. [2020] Hussein Hazimeh, Natalia Ponomareva, Petros Mol, Zhenyu Tan, and Rahul Mazumder.  The tree ensemble layer: Differentiability meets conditional computation.  _ArXiv_ , abs/2002.07772, 2020.
  * Breiman [2001] Leo Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, 2001.
  * Athey et al. [2019] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  _The Annals of Statistics_ , 2019.
  * Carreira-Perpinan and Tavallali [2018] Miguel A. Carreira-Perpinan and Pooya Tavallali.  Alternating optimization of decision trees, with application to learning sparse oblique trees.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.
  * Irsoy et al. [2012] Ozan Irsoy, O. T. Yildiz, and Ethem Alpaydin.  Soft decision trees.  _Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012)_ , pages 1819–1822, 2012.
  * Frosst and Hinton [2017] Nicholas Frosst and Geoffrey Hinton.  Distilling a neural network into a soft decision tree, 2017.
  * Lay et al. [2018] Nathan Lay, Adam P. Harrison, Sharon Schreiber, et al.  Random hinge forest for differentiable learning.  _CoRR_ , abs/1802.03882, 2018.
  * Bennett [1992] K.P. Bennett.  _Decision Tree Construction Via Linear Programming_.  Number no. 1067 in Computer sciences technical report. University of Wisconsin-Madison, Computer Sciences Department, 1992.
  * Bennett and Blue [1996] Kristin P. Bennett and Jennifer A. Blue.  Optimal decision trees.  Technical report, R.P.I. Math Report No. 214, Rensselaer Polytechnic Institute, 1996.
  * Bertsimas and Dunn [2017] Dimitris Bertsimas and Jack Dunn.  Optimal classification trees.  _Machine Learning_ , 106:1039–1082, 2017.
  * Zhu et al. [2020] Haoran Zhu, Pavankumar Murali, Dzung Phan, Lam Nguyen, and Jayant Kalagnanam.  A scalable mip-based method for learning optimal multivariate decision trees.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 1771–1781. Curran Associates, Inc., 2020.
  * Hyafil and Rivest [1976] Laurent Hyafil and Ronald L. Rivest.  Constructing optimal binary decision trees is np-complete.  _Information Processing Letters_ , 5(1):15–17, 1976.  ISSN 0020-0190.
  * Geurts et al. [2006] Pierre Geurts, Damien Ernst, and Louis Wehenkel.  Extremely randomized trees.  _Machine Learning_ , 63(1):3–42, 2006.
  * Prokhorenkova et al. [2018] Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, et al.  Catboost: Unbiased boosting with categorical features.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18, pages 6639–6649, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * Schapire and Freund [2012] Robert E. Schapire and Yoav Freund.  _Boosting: Foundations and Algorithms_.  The MIT Press, 2012.  ISBN 0262017180.
  * Zharmagambetov and Carreira-Perpiñán [2020] Arman Serikuly Zharmagambetov and Miguel Á. Carreira-Perpiñán.  Smaller, more accurate regression forests using tree alternating optimization.  In _ICML_ , 2020.
  * Jordan and Jacobs [1994] Michael I. Jordan and Robert A. Jacobs.  Hierarchical mixtures of experts and the em algorithm.  _Neural Comput._ , 6(2):181–214, mar 1994.  ISSN 0899-7667.
  * Tanno et al. [2019] Ryutaro Tanno, Kai Arulkumaran, Daniel C. Alexander, Antonio Criminisi, and Aditya V. Nori.  Adaptive neural trees.  _ArXiv_ , abs/1807.06699, 2019.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, et al.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, et al.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, et al., editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Dillon et al. [2017] Joshua V. Dillon, Ian Langmore, Dustin Tran, et al.  Tensorflow distributions.  _CoRR_ , abs/1711.10604, 2017.
  * Blanquero et al. [2021] Rafael Blanquero, Emilio Carrizosa, Cristina Molero-Río, and Dolores Romero Morales.  Optimal randomized classification trees.  _Computers & Operations Research_, 132:105281, 2021.  ISSN 0305-0548.
  * Hehn et al. [2019] Thomas M. Hehn, Julian F. P. Kooij, and Fred A. Hamprecht.  End-to-end learning of decision trees and forests.  _International Journal of Computer Vision_ , 128:997–1011, 2019.
  * Lee [2020] Simon C. K. Lee.  Addressing imbalanced insurance data through zero-inflated poisson regression with boosting.  _ASTIN Bulletin_ , 51:27 – 55, 2020.
  * Yirga et al. [2020] Ashenafi A. Yirga, Sileshi F. Melesse, Henry G. Mwambi, and Dawit G. Ayele.  Negative binomial mixed models for analyzing longitudinal cd4 count data.  _Scientific Reports_ , 10(1):16742, 2020.
  * Barron [2019] Jonathan T. Barron.  A general and adaptive robust loss function.  In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 4326–4334, 2019.
  * Jacobs et al. [2021] Wura Jacobs, Ehikowoicho Idoko, LaTrice Montgomery, Matthew Lee Smith, and Ashley L. Merianos.  Concurrent e-cigarette and marijuana use and health-risk behaviors among u.s. high school students.  _Preventive Medicine_ , 145:106429, 2021.  ISSN 0091-7435.
  * Young et al. [2017] Derek S. Young, Andrew M. Raim, and Nancy R. Johnson.  Zero-inflated modelling for characterizing coverage errors of extracts from the us census bureau’s master address file.  _Journal of The Royal Statistical Society Series A-statistics in Society_ , 180:73–97, 2017.
  * Moniz and Torgo [2018] Nuno Moniz and Luís Torgo.  Multi-source social feedback of online news feeds.  _ArXiv_ , abs/1801.07055, 2018.
  * Cameron and Trivedi [2013] A. Colin Cameron and Pravin K. Trivedi.  _Regression Analysis of Count Data_.  Econometric Society Monographs. Cambridge University Press, 2 edition, 2013.  doi: 10.1017/CBO9781139013567.
  * Lambert [1992] Diane Lambert.  Zero-inflated poisson regression, with an application to defects in manufacturing.  _Technometrics_ , 34(1):1–14, feb 1992.  ISSN 0040-1706.
  * Caruana [1997] Rich Caruana.  Multitask learning.  _Machine Learning_ , 28(1):41–75, 1997.
  * Ruder [2017] Sebastian Ruder.  An overview of multi-task learning in deep neural networks.  _ArXiv_ , abs/1706.05098, 2017.
  * Linusson [2013] Henrik Linusson.  Multi-output random forests.  2013\.
  * Xioufis et al. [2016] Eleftherios Spyromitros Xioufis, Grigorios Tsoumakas, William Groves, et al.  Multi-target regression via input space expansion: treating targets as inputs.  _Machine Learning_ , 104:55–98, 2016.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * Akujuobi and Zhang [2017] Uchenna Akujuobi and Xiangliang Zhang.  Delve: A dataset-driven scholarly search and analysis system.  _SIGKDD Explor. Newsl._ , 19(2):36–46, nov 2017\.  ISSN 1931-0145.
  * Vijayakumar and Schaal [2000] Sethu Vijayakumar and Stefan Schaal.  Locally weighted projection regression : An o(n) algorithm for incremental real time learning in high dimensional space.  2000\.
  * US Census Bureau [2021] US Census Bureau.  Planning database, Jun 2021.
  * Buitinck et al. [2013] Lars Buitinck, Gilles Louppe, Mathieu Blondel, et al.  API design for machine learning software: experiences from the scikit-learn project.  In _ECML PKDD Workshop: Languages for Data Mining and Machine Learning_ , pages 108–122, 2013.

## Supplementary Material

##  Appendix S1 Additional details for Experimental Section 6

#### Datasets

We use a collection of 27 open-source regression datasets from various domains
(e.g., social media platforms, human behavior, financial risk data). 9 of
these are from Mulan: A Java library for multi-label learning (Mulan) [Xioufis
et al., 2016], 2 of them are from University of California Irvine data
repository (UCI) [Dua and Graff, 2017], 12 of them are from Delve database
[Akujuobi and Zhang, 2017] and the 5 remaining are SARCOS 111The original test
set has significant data leakage as noted by
https://www.datarobot.com/blog/running-code-and-failing-models/. Following
their guidance we discard the original test set and use the original train set
to generate train/validation/test splits. [Vijayakumar and Schaal, 2000],
Youth Risk Behavior
Survey222https://www.cdc.gov/healthyyouth/data/yrbs/data.htm [Jacobs et al.,
2021], Block-Group level data from US Census Planning
Database333https://www.census.gov/topics/research/guidance/planning-
databases.html[US Census Bureau, 2021], and financial services loss data from
Kaggle444https://bit.ly/3swGnTo. For scm1d and scm20d from Mulan[Xioufis et
al., 2016], we consider the first 4 tasks (out of the 16 tasks in the original
dataset). For SARCOS, we consider 3 torques for prediction (torque-3, torque-4
and torque-7; we ignore the other torques as those seem to have poor
correlations with these.)

For all datasets, we split the datasets into $64\%/16\%/20\%$
training/validation/test splits. We train the models on the training set,
perform hyperparameter tuning on the validation set and report out-of-sample
performance on the test set.

###  S1.1 Tuning parameters and optimal depths comparison for a single tree in
Section 6.1

Table S1: Optimal depth of a single axis aligned and oblique decision tree on various regression datasets. |  Axis-Aligned |  Oblique
---|---|---
Data |  CART |  TAO |  TAO |  Soft Tree
abalone |  5 |  5 |  4 |  2
pumadyn-32nh |  6 |  6 |  4 |  3
pumadyn-32nm |  8 |  8 |  5 |  4
pumadyn-32fh |  3 |  3 |  3 |  5
pumadyn-32fm |  6 |  6 |  5 |  4
pumadyn-8nh |  6 |  6 |  5 |  3
pumadyn-8nm |  9 |  8 |  7 |  5
pumadyn-8fh |  5 |  5 |  4 |  2
pumadyn-8fm |  7 |  7 |  6 |  3
cpu |  10 |  8 |  5 |  5
cpuSmall |  8 |  8 |  5 |  5
concrete |  15 |  8 |  5 |  9

#### TAO Implementation

We wrote our own implementation of the TAO algorithm proposed in Carreira-
Perpinan and Tavallali [2018]. We considered binary trees with TAO for both
axis-aligned and oblique trees. In the case of axis-aligned splits, we
initialize the tree with CART solution and run TAO iterations until there is
no improvement in training objective. In the case of oblique trees, we
initialize with a complete binary tree with random parameters for the
hyperplanes in the split nodes and use logistic regression to solve the
decision node optimization. We run the algorithm until either a maximum number
of iterations are reached or the training objective fails to improve.

#### Tuning parameters

For CART and axis-aligned TAO, we tune the depth in the range $[2-20]$. We
also optimize over the maximum number of iterations in the interval
$[20-100]$. For oblique TAO and soft tree, we tune the depth between $2-10$
and the number of iterations between $20-100$. Additionally, for soft tree, we
also tune over the learning rates $[1e-5,1e-2]$ with Adam optimizer and batch
sizes $\\{64,128,256,512\\}$. For a fair comparison, we run all $4$ methods
for $100$ trials.

#### Optimal depths

We make a comparison of optimal depths between CART, TAO (both axis-aligned
and oblique) and soft tree. The soft tree finds a smaller optimal depth in
majority cases in comparison with its classical counterpart i.e., oblique TAO
tree — See Table S1. This is hypothesized to be due to the end-to-end
optimization done by soft tree as opposed to a local search performed by the
TAO algorithm.

###  S1.2 Tuning parameters for Sections 6.2

We use HistGBDT from sklearn Buitinck et al. [2013] (GBDT in sklearn does not
support Poisson regression). We tune over depths in the range $[2-20]$, number
of trees between $50-1500$ and learning rates on the log-uniform scale in the
interval $[1e-5,1e-1]$. For differentiable tree ensembles, we tune number of
trees in the range $[2,100]$, depths in the set $[2-4]$, batch sizes
$\\{64,128,256,512\\}$, learning rates $[1e-5,1e-1]$ with Adam optimizer and
perform early stopping with a patience of $25$ based on the validation set.
For all models, we perform a random search with 1000 hyperparameter tuning
trials.

###  S1.3 Tuning parameters for Sections 6.3

We tune number of trees in the interval $[50-1500]$ for RF, GRF and GBDT. For
RF and GBDT, we also tune over depths between $2-20$. For GBDT, we tune
learning rates between $[1e-5,1e-1]$. For GRF, we also tune over min_node_size
in the set $[2-20]$ and $\alpha\in[1e-3,1e-1]$. For single-task and multi-task
trees, we tune over depths $[2-4]$, number of trees $[5-100]$, batch sizes
$\\{64,128,256,512\\}$, epochs $[20-500]$, Adam learning rates $[1e-5,1e-2]$.
We also optimize over the regularization penalty for multi-task soft decision
trees $[1e-5,1e1]$. All single-task models (soft tree ensembles, GBDT) are
tuned for 1000 trials per task. All multi-task models (RF, GRF, multi-task
soft trees) are tuned for 1000 trials in total.

###  S1.4 Tuning parameters for Sections 6.4

We use GBDT from XGBoost [Chen and Guestrin, 2016], where we tune number of
trees in the interval $[50-1500]$, depths between $2-20$ and learning rates
between $[1e-4,1e-0]$. For single-task and multi-task trees, we tune over
depths $[2-4]$, number of trees $[5-100]$, batch sizes $\\{64,128,256,512\\}$,
epochs $[20-200]$, Adam learning rates $[1e-5,1e-2]$. We also optimize over
the regularization penalty for multi-task soft decision trees $[1e-5,1e1]$.
All single-task models (soft tree ensembles, GBDT) with Poisson, Zero-
Inflated-Poisson, Negative Binomial are tuned for 1000 trials per task. All
multi-task soft-trees are tuned for 1000 trials in total.


# Moiré Gravity and Cosmology

Alireza Parhizkar  Joint Quantum Institute, University of Maryland, College
Park, MD 20742, USA  Victor Galitski  Joint Quantum Institute, University of
Maryland, College Park, MD 20742, USA

###### Abstract

The vacuum catastrophe is a fundamental puzzle, where the observed scales of
the cosmological constant are many orders of magnitude smaller than the
natural scales expected in the theory. This work proposes a new “bi-world”
construction that may offer an insight into the cosmological constant problem.
The model includes a $(3+1)$-dimensional manifold with two different
geometries and matter fields residing on them. The diffeomorphism invariance
and causality highly constrain the two metrics to be conformally related,
$\eta_{\mu\nu}=\phi^{2}g_{\mu\nu}$. This reduces the theory to a standard
single-world description, but introduces a new inherently geometrical “moiré
field,” $\phi$. Interestingly, the moiré field has the character of both a
dilaton and Higgs field familiar in the conventional theory. Integrating out
the moiré field naturally gives rise to the Starobinsky action and
inflationary dynamics. In the framework of the Friedmann-
Lemaitre–Robertson–Walker solution, we reduce an effective action for the
moiré field to that of a particle moving in a Mexican hat potential. The
equations of motion are then solved numerically and the moiré field is shown
to approach a Mexican-hat minimum in an oscillatory fashion, which is
accompanied by the decay of the Hubble parameter. Under additional reasonable
assumptions, the vacuum energy asymptotically approaches zero in the end of
inflationary evolution. The physics presented here shares similarities with
the moiré phenomena in condensed matter and elsewhere, where two similar
structures superimposed upon give rise to a superstructure with low emergent
energy scales compared to the native theories.

Introduction—As ubiquitous phenomena, moiré patterns are not hard to come by.
Generally when two layers of grids combine, like overlapping fabrics or when a
digital photo of a pixel screen is watched through another such screen, an
additional larger pattern emerges. When the original layers are close enough,
the moiré pattern becomes more than just an optical illusion. For example, in
bilayer graphene the moiré pattern dictates the interlayer hopping process
Parhizkar and Galitski (2021); Balents (2019). This concept is generalizable
to continuous models as in Parhizkar and Galitski (2021) where we can define
moiré physics to be: The emergence of large length scales, or synonymously
small energy scales, due to coalition of two underlying theories with small
length or large energy characteristics.

Moiré physics as a conceptual tool is potentially utilizable in many different
contexts. Here we are going to investigate its possible presence in
gravitational systems and its implication in cosmology.

Bi-world—By definition, for a moiré pattern to appear, we first need two more-
or-less similar systems as the underlying structures. In case of gravity those
are metric systems and therefore we are combining two curved space-times, and
potentially the matter fields residing on them, with each other. Such a
universe that embraces two worlds is described most generally by an action
$S=S_{g}+S_{\eta}+S_{c}+S_{a}$ where $S_{g}$ and $S_{\eta}$ are built only out
of fields that live solely on either $g$-world or $\eta$-world, while $S_{c}$
contains the cross terms consisted of these fields. But there is also the
possibility of amphibian fields that live on both worlds, the physical
description of which is encapsulated in $S_{a}$. If not through $S_{c}$ and
$S_{a}$ the g-world and $\eta$-world could not talk and would have been
completely independent of each other.

Naturally, we now have to integrate over all fields contributing to $S$. We
usually take $S_{\eta}$ to be a copy of $S_{g}$ and then $S_{c}$ will contain
various coupling terms between fields in g-world and their counterparts in
$\eta$-world. In this fashion the path-integral looks like the following,

$\displaystyle
I=\int\mathcal{D}g_{\mu\nu}\mathcal{D}\eta_{\mu\nu}\mathcal{D}\left[\bar{\Psi}^{g},\Psi^{g}\right]\mathcal{D}\left[\bar{\Psi}^{\eta},\Psi^{\eta}\right]e^{iS}\,,$
(1)

where we have excluded gauge fixing and Faddeev-Popov fields, and have used
fermions as a representative for all matter fields to avoid clutter where it
is not necessary. Superscripts $\eta$ and $g$ label the native world of each
field. There remains, hence, no ambiguity as to which metric or vielbein
should be used for raising and lowering indices of the fields.

An example of a typical action containing intra-world and inter-world bilinear
terms is the following,

$\tilde{S}=\int
d^{4}x\left[\bar{\Psi}^{g}\Psi^{g}+\bar{\Psi}^{g}\Psi^{\eta}+\bar{\Psi}^{\eta}\Psi^{g}+\bar{\Psi}^{\eta}\Psi^{\eta}\right]\,.$
(2)

Variables $\Psi^{g,\eta}$ are defined to be carrying a forth root of the
corresponding metric determinant as $\Psi^{g}\equiv\sqrt[4]{|g|}\psi^{g}$ with
$\psi^{g}$ being independent of metric. Writing the action in terms of
$\psi^{g,\eta}$ makes the diffeomorphism invariance manifest. The weight half
density variables $\Psi^{g,\eta}$ and their consequent measures
$\mathcal{D}\Psi^{g,\eta}$ are suitable for curved background considerations
(cf. Hawking (1977); Fujikawa  _et al._ (2004); Fujikawa (1983)). Looking at
the above we see that along with the cross-terms in $S_{c}$ comes the
crossbreed volume element $\sqrt[4]{g\eta}d^{4}x$. Consequently, out of the
bi-world construction a new inter-world vacuum energy density will emerge
which couples to gravity as $2\bar{\Lambda}\sqrt[4]{g\eta}$ in the effective
form of $S_{c}$.

Interestingly, $2\bar{\Lambda}\sqrt[4]{g\eta}$ is the on-shell value of
$\bar{\Lambda}\sqrt{|g|}e^{\chi}+\bar{\Lambda}\sqrt{|\eta|}e^{-\chi}$ when
solving for the auxiliary scalar $\chi$, which describes the simultaneous
concentration and dilution of vacuum energy density in g-world and
$\eta$-world respectively. The saddle-point occurs at
$e^{\chi}=\sqrt[4]{\eta/g}$. The ratio of the metric determinants, which are
tensor densities of the same weight, is a scalar under general coordinate
transformations: $\eta$ transforms in exactly the same way as $g$ does and for
their ratio the resulting Jacobians cancel.

Terms that describe a vacuum energy density, in Einstein’s equations are
proportional to metric, i.e. $\rho g_{\mu\nu}$. Since metric is locally
Lorentz invariant, boosting observers does not result in them measuring a
different vacuum energy density, therefore they remain unable to find a
preferred inertial frame only by looking at a stationary vacuum. Terms in
Einstein’s equations that are proportional to metric, come from those terms in
the action that depend on metric only through its determinant. With just one
metric, the only such term we can construct, while preserving diffeomorphism
invariance, is the four volume $\int d^{4}x\sqrt{|g|}$. With two metrics,
however, we will have a new scalar, $\sqrt[4]{\eta/g}$, and now we can use any
power of this scalar to construct various terms that can describe vacuum as
explained above. Let us define $\phi^{2}\equiv\sqrt[4]{\eta/g}$ for future
convenience.

If either metric, or both, are dynamical we can expect a kinematic term for
$\phi$ to appear in the action. Therefore, $\phi$ becomes a suitable candidate
for a field that describes the vacuum and possibly its dynamical relaxation.
But $\phi$ is more than that. If we factor out a square root of metric
determinant from terms in $S_{c}$ we can write $S_{c}$ seemingly with a
regular volume element. Then, considering the fact that fields hop from one
world to another, we can anticipate the following few instances in $S_{c}$,

$\displaystyle S_{c}=\int$ $\displaystyle
d^{4}x\sqrt{|g|}\Big{[}m\phi^{2}\bar{\psi}^{g}\psi^{\eta}+w\phi^{2}G^{\eta}_{\mu}G^{g\mu}+$
$\displaystyle\alpha\left(\phi^{2}+\phi^{-2}\right)\epsilon^{\mu\nu\rho\sigma}F^{g}_{\mu\nu}F^{\eta}_{\rho\sigma}+g\leftrightarrow\eta\Big{]}+\cdots\,,$
(3)

where $m$, $w$ and $\alpha$ are coupling parameters which symbolically
represent inter-world hopping between different fermionic fields, gauge fields
and field strengths respectively. Also
$\epsilon^{\mu\nu\rho\sigma}\equiv\varepsilon^{\mu\nu\rho\sigma}/\sqrt{|g|}$
is the Levi-Civita tensor with $\varepsilon^{\mu\nu\rho\sigma}=\pm 1$ being
the totally anti-symmetric symbol. Therefore, $\phi$ has a Higgs
characteristic as well. Other possible terms to include in $S_{c}$ are
variations of $\int
d^{4}x\sqrt{|g|}\phi^{2}\bar{\psi}^{\eta}\not{D}^{g}\psi^{g}$ with
$\not{D}^{g}$ being the Dirac operator on the curved geometry given by
$g_{\mu\nu}$. This of course will modify the chiral anomaly and the non-
conservation of the chiral current
$\nabla^{g}_{\mu}j^{\mu}_{5}=\mathcal{A}_{5}$ will have additional
contributions that are coupled to $\phi$ such as $\phi^{2}F\wedge
F\subset\mathcal{A}_{5}$ which describes a scalar boson decay into vector
bosons carried through a triangle of fermions Shifman  _et al._ (1979). It is
also worth noting that the first term in the second line of Eq. 3 is a
naturally arising axionic term Wilczek (1978); Weinberg (1978); Sikivie
(1983), however, a detailed and rigorous investigation of all implications of
the moiré field $\phi$ will be considered elsewhere. Below, we instead focus
on its potential cosmological consequences.

Bi-causality—On a manifold with only one metric, say $g_{\mu\nu}$, it is
always possible to transform the coordinates to one in which the metric is
locally Minkowskian where the line element is expressed as
$ds^{2}_{g}=-dt^{2}+\mathbf{dx}^{2}$. This of course is not generally possible
for two metrics: We can not always find a local basis in which the line
elements match as $ds^{2}_{g,\eta}=-dt^{2}+\mathbf{dx}^{2}$. (See Fig. 1)
Among other things, this means that in general the light cones of g-world and
$\eta$-world do not coincide. On the whole, the gauge classes of a bi-metric
$(g_{\mu\nu},\eta_{\mu\nu})$ are smaller than the product of their separate
gauge classes.

Figure 1: Events that affect $P$, namely its past $J^{-}(P)$, reside in the
hatched regions. Horizontal hatching designates the past of $P$ according to
$g_{\mu\nu}$ while vertical hatching is according to $\eta_{\mu\nu}$. (a) The
two pasts and futures of $P$ do not match in general. The region
$\xi\equiv\Sigma\cap J^{+}_{g}(P)\setminus J^{+}_{\eta}(P)$ is then
effectively forbidden for excitations of $\hat{a}(x)$. Therefore, in the space
of all configurations of $g_{\mu\nu}$ and $\eta_{\mu\nu}$ those with
coinciding light cones, as in (b), are statistically favored by an exponential
factor of the volume of $\xi$.

Metric is responsible both for setting up the causal structure of spacetime
and for determining the four-volume element at each point on it. The latter is
given by the metric determinant and the former by the remaining degree of
freedom in four dimensions. An observer can, in principle, reconstruct the
causal structure by gathering which measurements performed at different points
of spacetime affect each other. Measurements with a space-like separation must
be independent, and hence we expect, the commutators $[a(x),a(y)]$ with
$(y-x)^{2}>0$ to vanish Peskin (2018). So microcausality ensures that
superluminal propagations are exponentially suppressed; but superluminal with
respect to which metric?

Consider, for example, an amphibian field $a(x)$. Using the local commutator
data $[a(x),a(y)]$ at a point $P$, segregating disconnected events to $P$ from
the rest, we should be able to sketch the light cones of both g-world and
$\eta$-world at $P$ simultaneously. This means the two metrics must eventually
end up sharing their causal structure. Therefore, one degree of freedom is
effectively thrown away and what remains is a causal structure and two metric
determinants. Let us look at the probability amplitude of creating $a$ at $P$
and finding it at a later time anywhere,

$\displaystyle\int_{V}\langle
a(x)a(P)\rangle\equiv\frac{\int_{V}\int\mathcal{D}g_{\mu\nu}\mathcal{D}\eta_{\mu\nu}\mathcal{D}a\,e^{iS}a(x)a(P)}{\int\mathcal{D}g_{\mu\nu}\mathcal{D}\eta_{\mu\nu}\mathcal{D}a\,e^{iS}}\,.$
(4)

Since the causal structures described by the two metrics do not necessarily
match, $P$ in general has two futures and two pasts, each given by one of the
metrics. The path-integral is invariant under coordinate transformations, so
we can always locally choose a coordinate where the fastest light cone is on
$t=\pm|\mathbf{x}|$. The integral $\int_{V}$ is over all points on an
arbitrary space-like hypersurface $\Sigma$ that avoids crossing the past of
$P$. But since the propagator is exponentially suppressed for all $x$ that sit
outside any of the two light cones, the integral effectively only covers the
events inside the narrowest future light cone. Therefore, when we are
integrating over all $g_{\mu\nu}$ and $\eta_{\mu\nu}$, those configurations of
the two metrics which do not share causal structures have, in effect, less
accessible volume and hence exponentially less microstates in comparison to
those with coinciding light cones. The light cones are sketched by
$ds^{2}_{g,\eta}=0$. Thus, when the causal structure is set, there remains a
conformal degree of freedom, as in
$\phi^{2}ds^{2}_{g}=\phi^{2}(-dt^{2}+\mathbf{dx}^{2})=ds^{2}_{\eta}$ where the
light cones of both worlds are coinciding on $t=\pm|\mathbf{x}|$.
Consequently, the contributing configurations are expected to approximately be
those with conformally related metrics

$\eta_{\mu\nu}=\sqrt[4]{\frac{\eta}{g}}\,g_{\mu\nu}=\phi^{2}g_{\mu\nu},$

which directs us to rewrite the geometrical part of path-integration (1) in
terms of the moiré field, $\phi$, as follows

$\int\mathcal{D}g_{\mu\nu}\mathcal{D}\eta_{\mu\nu}e^{iS[g_{\mu\nu},\eta_{\mu\nu}]}\approx\int\mathcal{D}g_{\mu\nu}\mathcal{D}\phi\,e^{iS[g_{\mu\nu},\phi^{2}g_{\mu\nu}]}.$
(5)

Inflation and dynamic relaxation— Let us now, for the sake of simplicity,
resort to considering only two dynamical metrics with the usual Einstein-
Hilbert action each, and minimize the appearance of matter fields to three
accumulated vacuum energy densities,

$\displaystyle S=\int d^{4}x\Big{[}$
$\displaystyle\sqrt{|g|}\left(R^{g}-2\Lambda^{g}\right)$ (6)
$\displaystyle+\sqrt{|\eta|}\left(R^{\eta}-2\Lambda^{\eta}\right)-4\bar{\Lambda}\sqrt[4]{g\eta}\Big{]}\,,$

where $R^{g,\eta}$ are the Ricci scalars built out of each metric,
$\Lambda^{g,\eta}$ the intra-world cosmological constants and $\bar{\Lambda}$
the inter-world cosmological constant coming from the tunneling terms.
Classically, it is enough to solve the Einstein’s equations derived from
taking the variation of the above action. $\delta S/\delta g^{\mu\nu}$ yields,

$R^{g}_{\mu\nu}-\frac{1}{2}R^{g}g_{\mu\nu}+\Lambda^{g}+\bar{\Lambda}\sqrt[4]{\frac{\eta}{g}}=0\,,$
(7)

and the equations of motion for $\eta_{\mu\nu}$ are given by
$g\leftrightarrow\eta$. With $\Lambda^{g}$ and $\bar{\Lambda}$ being
constants, the contracted Bianchi identity makes sure that
$\phi^{2}=\sqrt[4]{\eta/g}$ is constant as well. This constant will be
dictated by the equations of motion and the symmetries imposed on the desired
solution. For instance, for the case of $\Lambda^{g}=\Lambda^{\eta}$ any value
of $\phi^{2}$ other than one will result in different vacuum energy densities
for $g_{\mu\nu}$ and $\eta_{\mu\nu}$ which for an expanding universe means a
time varying ratio of metric determinants, in contradiction with a changeless
$\phi$. Thus in that case $\phi^{2}=1$ classically. However, we will see that
the quantum effects discussed in the previous section will allow $\phi$ to
vary. But we expect such quantum gravitational deviations from the classical
theory only at huge energy or temperature scales. Therefore, the effective
theory should descend to the classical regime in the low-energy limit (which
is equivalent to the long-time limit in the context of cosmology).

Following (5) we can change the action (6) to its effective form by
substituting $\eta_{\mu\nu}\rightarrow\phi^{2}g_{\mu\nu}$ everywhere,

$\displaystyle S=\int d^{4}x\sqrt{|g|}\Big{[}$
$\displaystyle\left(1+\phi^{2}\right)R+6k\partial_{\mu}\phi\partial^{\mu}\phi$
(8)
$\displaystyle\quad-2\Lambda^{\eta}\phi^{4}-2\Lambda^{g}-4\bar{\Lambda}\phi^{2}\Big{]}\,.$

Here we have omitted the surface terms and $k$ represents the effect of
unstated parts of the action (6) on the kinetic term for $\phi$. We have also
dropped the superscript on $R\equiv R^{g}$. The above action has one metric
$g_{\mu\nu}$ and one scalar $\phi^{2}$ as its dynamical degrees of freedom
which are going to be integrated over in the path-integral approach.
Therefore, apart from its similarities to Higgs and a vacuum order parameter,
$\phi^{2}$ also resembles a dilaton field.

The action (6), and consequently the above action (8), are devised for Neumann
boundary conditions Krishnan  _et al._ (2016). This fact particularly means
that in the effective action the initial condition for $\phi$ is set by fixing
its derivative on a generic Cauchy surface, so that we can begin in an “slow-
roll” regime where we can disregard the kinetic term. Then the saddle-point
approximation, using the on-shell value of
$\phi^{2}\doteq(R-4\bar{\Lambda})/4\Lambda^{\eta}$, will give us a
renormalized Starobinsky action:

$S\doteq\int
d^{4}x\sqrt{|g|}\left[R-2\Lambda^{g}+\frac{\left(R-4\bar{\Lambda}\right)^{2}}{4\Lambda^{\eta}}\right]\,.$
(9)

The Starobinsky action Starobinsky (1980); Vilenkin (1985) and in general
higher curvature terms can arise from quantum gravity considerations Davies
(1977). The recurrence of the above, therefore, might be a hint towards the
bi-world construction having a similar origin.

To achieve a minimally coupled theory of gravity out of (8) we conformally
transform our (now only) metric by $g_{\mu\nu}\rightarrow
g_{\mu\nu}/(1+\phi^{2})$ to get,

$\displaystyle S=\int d^{4}x\sqrt{|g|}\Bigg{[}R$
$\displaystyle-6\frac{(1-k)\phi^{2}-k}{(1+\phi^{2})^{2}}\partial_{\mu}\phi\partial^{\mu}\phi$
$\displaystyle\
\quad-2\frac{\Lambda^{\eta}\phi^{4}+\Lambda^{g}+2\bar{\Lambda}\phi^{2}}{(1+\phi^{2})^{2}}\Bigg{]}\,,$
(10)

where the dynamics governing $\phi$ is more clear. Since we are interested in
the cosmological consequences here we will simplify the model by choosing the
Friedmann-Lemaître–Robertson–Walker metric $ds^{2}\equiv
g_{\mu\nu}dx^{\mu}dx^{\nu}=-dt^{2}+a^{2}(t)d\mathbf{x}^{2}$ as an ansatz to
the equations of motion. These are,

$\displaystyle\ddot{\phi}=-3H\dot{\phi}+\frac{\phi}{(1-k)\phi^{2}-k}\Bigg{[}\frac{(1-k)\phi^{2}-k-1}{1+\phi^{2}}\dot{\phi}^{2}$
$\displaystyle-\frac{2}{3}\frac{(\Lambda^{\eta}-\bar{\Lambda})\phi^{2}-(\Lambda^{g}-\bar{\Lambda})}{1+\phi^{2}}\Bigg{]}\,,$
(11)

with $H\equiv\dot{a}/a$ being the Hubble parameter. Depending only on time,
$\phi(t)$ now has the dynamics of a single particle rolling on a potential.
The first and second terms on the right hand side are friction exerted on this
particle while the last one is the conservative force which is always zero at
$\phi^{2}=(\Lambda^{\eta}-\bar{\Lambda})/(\Lambda^{g}-\bar{\Lambda})$. The
potential corresponding to this force is a Mexican hat except for $k=0$ where
the crown of the hat is uplifted to infinity. The constraint part of
Einstein’s equations is given by,

$3H^{2}=3\frac{(1-k)\phi^{2}-k}{(1+\phi^{2})^{2}}\dot{\phi}^{2}+\frac{\Lambda^{\eta}\phi^{4}+\Lambda^{g}+2\bar{\Lambda}\phi^{2}}{(1+\phi^{2})^{2}}\,,$
(12)

while the dynamical part is,

$2\dot{H}+3H^{2}=-3\frac{(1-k)\phi^{2}-k}{(1+\phi^{2})^{2}}\dot{\phi}^{2}+\frac{\Lambda^{\eta}\phi^{4}+\Lambda^{g}+2\bar{\Lambda}\phi^{2}}{(1+\phi^{2})^{2}}\,,$
(13)

which combine to simply give,

$\dot{H}=-3\frac{(1-k)\phi^{2}-k}{(1+\phi^{2})^{2}}\dot{\phi}^{2}\,.$ (14)

Although not necessary for our purposes, hereon we adopt the symmetric
situation where the two worlds are similar whence the desired destiny is for
them to eventually reach equilibrium at $\phi^{2}=1$. According to (14), for
all values of $k$ below $1/2$ the Hubble parameter will be ever decreasing up
to the point of equilibrium. Therefore, we can imagine the story of the
universe beginning with the two worlds either already out of equilibrium or
kicked out of it, while according to equation (12) the state of non-
equilibrium, $\phi^{2}\neq 1$, is accompanied with a huge vacuum energy
density. All this however will be relaxed when the worlds eventually reach
equilibrium with the help of the friction forces. An example has been plotted
in Fig. 2. At equilibrium $\phi$ remains constant and hence the classical
theory (7) is reinstated.

Figure 2: Evolutions of the moiré field $\phi$ and the Hubble parameter $H$ in
units of Planck time $t_{P}$. Initially the worlds are away from equilibrium
by $\phi=1.5$. Here $k$ is set to zero while the cosmological constants are
set to be $\Lambda^{g,\eta}=-\bar{\Lambda}=M_{P}^{2}$. $H$ (the plot on top in
blue) starts from $0.22M_{P}$ and approaches zero as time goes by while $\phi$
starts from $1.5$ and oscillates around $\phi=1$ henceforward. Whenever the
rate of change in $\phi$ is small the universe experiences a period of
stationary $H$.

Having two similar worlds implies that their characteristic parameters
(coupling constants, cosmological constants, etc.) are the same, e.g.
$\Lambda^{\eta}=\Lambda^{g}\approx M_{P}^{2}$. This resembles a system of two
overlapping grids with characteristic lengths of $L_{P}$. The zero of the
conservative force in (11) happens at $\phi^{2}=1$. As discussed before, the
classical limit further requires this point of equilibrium to be stable which
is met only if $\bar{\Lambda}$ has a different sign than $\Lambda^{g,\eta}$.
On the other hand, the cross terms in $S_{c}$ are built out of the same fields
and parameters appearing in $S^{g,\eta}$. There is an inter-world hopping term
from g-world to $\eta$-world in addition to its conjugate which describes the
reverse hopping. Therefore the magnitude of inter-world vacuum energy density
must be twice that of the intra-world ones in the action (6). We can reinforce
this picture by looking at a few possible diagrams of the gravitating bi-world
vacuum Fig. 3

Figure 3: Coupling of external gravitons (wiggly lines) to vacuum loops (solid
lines) of intra-world vacua (on the left) and inter-world vacua (on the
right).

where for each vacuum loop in a single world there is a counterpart that
passes through both worlds. All in all, we have
$\Lambda^{g,\eta}=-\bar{\Lambda}\approx M_{P}^{2}$. This combination yields
the desired history for the universe (Fig. 2) starting from an inflationary
period and reaching a small vacuum energy density that is continuously
approaching zero at late times.

Conclusion—Although the classical limit in a way ensures it, a robust and
detailed consideration of why the cosmological constants have such a relation
requires a solid knowledge of the origin of the bi-world in the proposed
scenario. Discounting the classical limit we are still left with some degree
of maneuver. For example, we know that the inter-world quantum vacuum is
constructed out of the same fields and processes as the intra-world quantum
vacuum. Since there are no new scales introduced in the theory we expect the
vacuum energy densities to be of equal magnitude, which leaves us merely with
deciding what sign must be attributed to inter-world vacuum energy density. It
is worth noting that $\mathbb{Z}^{g,\eta}_{2}$ invariant terms in $S_{g,\eta}$
have counterparts in $S_{c}$ which do not remain invariant under either
$\mathbb{Z}^{g}_{2}$ or $\mathbb{Z}^{\eta}_{2}$ separately. Thus one might
naively decide that the sign choice is arbitrary and/or a property attributed
to a symmetry broken universe. In any case, the negative sign fits our picture
more naturally: If the worlds were lying on membranes separated by some
distance as in Fig. 3 the tunneling terms would have exerted an attractive or
repulsive force between them and their distance would have settled on the
point of equilibrium between these forces. Albeit the separation in Fig. 3 is
only for illustration purposes. There is no separation between the membranes
in the bi-world structure which in membrane picture translates to an always
attractive force between the membranes. This in turn means that the
corresponding potential for the inter-layer force appears in the Lagrangian as
$\mathcal{L}\subset-V(l)$ with $l$ being the separation between the membranes
and $\frac{\text{d}^{2}V}{\text{d}l^{2}}>0$. For terms in $S_{c}$ to result in
such a potential they must appear with a minus sign themselves which results
in a negative inter-world vacuum energy density given the intra-world ones are
positive.

The theory presented here is experimentally falsifiable also through
considering the detailed behaviour of the moiré field $\phi$. For example, the
axion term in (3) when integrated out will modulate the effective moiré field
potential in an oscillatory manner Peccei and Quinn (1977); Abbott (1985) near
the extrema, which possibly leaves a trace to be found by early universe
physics observations.

Lastly, it should be mentioned that the construction briefly sketched in this
Letter differs from theories of bimetric or massive gravity Hassan and Rosen
(2012, 2011); Baccetti  _et al._ (2012) in several respects of which we will
emphasise on the following: The metrics couple to each other and to matter in
a different manner as we have here; the current construction generally
describes a universe containing two worlds rather than only two metrics, it
includes matter fields the effects of which are crucial and measurable at
least through cosmological observations; there is an emergent inherently
geometric moiré field $\phi$ that plays multiple roles within the theory;
gravity remains massless here; and the foundational philosophy is different
promising a feasible origin story.

###### Acknowledgements.

This work originated from a project supported by the Templeton Foundation.
This research was also supported by the Simons Foundation. A.P. wishes to
thank Batoul Banihashemi and Ted Jacobson for useful discussions.

## References

  * Parhizkar and Galitski (2021) A. Parhizkar and V. Galitski, (2021), arXiv:2108.04252 [cond-mat.mes-hall] .
  * Balents (2019) L. Balents, SciPost Phys. 7, 48 (2019).
  * Hawking (1977) S. W. Hawking, in  _Euclidean quantum gravity_ (World Scientific, 1977) pp. 114–129.
  * Fujikawa  _et al._ (2004) K. Fujikawa, K. Fujikawa, H. Suzuki, _et al._ ,  _Path integrals and quantum anomalies_ , 122 (Oxford University Press on Demand, 2004).
  * Fujikawa (1983) K. Fujikawa, Nuclear Physics B 226, 437 (1983).
  * Shifman  _et al._ (1979) M. A. Shifman, A. I. Vainshtein, M. B. Voloshin, and V. I. Zakharov, Sov. J. Nucl. Phys. 30, 711 (1979).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Sikivie (1983) P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * Peskin (2018) M. E. Peskin,  _An introduction to quantum field theory_ (CRC press, 2018).
  * Krishnan  _et al._ (2016) C. Krishnan, K. P. Kumar, and A. Raju, Journal of High Energy Physics 2016, 43 (2016).
  * Starobinsky (1980) A. Starobinsky, Physics Letters B 91, 99 (1980).
  * Vilenkin (1985) A. Vilenkin, Phys. Rev. D 32, 2511 (1985).
  * Davies (1977) P. Davies, Physics Letters B 68, 402 (1977).
  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Abbott (1985) L. Abbott, Physics Letters B 150, 427 (1985).
  * Hassan and Rosen (2012) S. F. Hassan and R. A. Rosen, Journal of High Energy Physics 2012, 1 (2012).
  * Hassan and Rosen (2011) S. F. Hassan and R. A. Rosen, Journal of High Energy Physics 2011, 1 (2011).
  * Baccetti  _et al._ (2012) V. Baccetti, P. Martin-Moruno, and M. Visser, Classical and Quantum Gravity 30, 015004 (2012).