# Distilling Interpretable Models into Human-Readable Code

Walker Ravina∗, Ethan Sterling, Olexiy Oryeshko, Nathan Bell, Honglei Zhuang,
Xuanhui Wang, Yonghui Wu, Alexander Grushetsky  Google, Mountain View, CA, USA
walkerravina, esterling, olexiy, nathanbell, hlz, xuanhui, yonghui,
<EMAIL_ADDRESS>
###### Abstract.

The goal of model distillation is to faithfully transfer teacher model
knowledge to a model which is faster, more generalizable, more interpretable,
or possesses other desirable characteristics. Human-readability is an
important and desirable standard for machine-learned model interpretability.
Readable models are transparent and can be reviewed, manipulated, and deployed
like traditional source code. As a result, such models can be improved outside
the context of machine learning and manually edited if desired. Given that
directly training such models is difficult, we propose to train interpretable
models using conventional methods, and then distill them into concise, human-
readable code.

The proposed distillation methodology approximates a model’s univariate
numerical functions with piecewise-linear curves in a localized manner. The
resulting curve model representations are accurate, concise, human-readable,
and well-regularized by construction. We describe a piecewise-linear curve-
fitting algorithm that produces high-quality results efficiently and reliably
across a broad range of use cases. We demonstrate the effectiveness of the
overall distillation technique and our curve-fitting algorithm using four
datasets across the tasks of classification, regression, and ranking.

Model distillation; human readable; piecewise-linear curves

∗Corresponding author

††ccs: Computing methodologies Machine learning approaches

##  1\. Introduction

Interpretable models are critical for high-stakes decision-making scenarios
(Rudin, 2018) such as guiding bail or parole decisions, assessing loan
eligibility, and guiding medical treatment decisions. In these cases, the
explanation of a model’s output (_e.g_. individual feature contributions)
should be examinable and understandable, to ensure transparency,
accountability, and fairness of the outcomes.

To achieve intrinsic interpretability, univariate functions are widely used in
interpretable models. In the classic Generalized Additive Models (GAMs)
(Hastie and Tibshirani, 1986), the model is a sum of univariate shape
functions,

$M=f_{0}+f_{1}(x_{1})+f_{2}(x_{2})+f_{3}(x_{3})+\dots+f_{n}(x_{n}).$

where $x_{i}$’s are $n$ features and $f_{i}$’s are the shape functions. Such a
model is simple but often less accurate than a model with feature
interactions. Recently, Lou _et al_. (Lou et al., 2013) showed that adding a
limited number of pairwise feature interactions allows GAM-style additive
models to capture a significant fraction of the accuracy of a fully-
interacting model. In many cases of interest, such feature interactions are
intuitively captured with products of univariate functions,

$g_{1}(c_{1})\cdot f_{1}(x_{1})+g_{2}(c_{2})\cdot f_{2}(x_{2})+\dots,$

or products of groups of features,

$(g_{1,1}(c_{1})+g_{1,2}(c_{2}))\cdot
f_{1}(x_{1})+(g_{2,1}(c_{1})+g_{2,2}(c_{2}))\cdot f_{2}(x_{2})+\dots,$

where the magnitude of one function (_i.e_. $f_{i}$) is modulated by a
function (_i.e_. $g_{i}$ or $g_{i,j}$) of another ”context” feature (_i.e_.
$c_{i}$) (Zhuang et al., 2021). In other cases, the interaction amongst
features is adequately approximated by additive models of univariate functions
nested within univariate functions,

$\displaystyle f(x_{1},x_{2},x_{3})\approx$ $\displaystyle\
g_{1}(f_{1,1}(x_{1})+f_{1,2}(x_{2})+f_{1,3}(x_{3}))\ +$ $\displaystyle\
g_{2}(f_{2,1}(x_{1})+f_{2,2}(x_{2})+f_{2,3}(x_{3}))+\dots,$

where the outer function $g_{i}$ captures nonlinear behavior (Chen et al.,
2018). Indeed, the Kolmogorov–Arnold representation theorem (Kolmogorov, 1957;
Wikipedia, 2020c) guarantees that every continuous multivariate function of
$n$ inputs can be represented as a sum of $2n$ such terms,

$f(x_{1},\dots,x_{n})=\sum_{i=0}^{2n}g_{i}\left(\sum_{j=1}^{n}f_{i,j}(x_{j})\right).$

In practice a single outer function is often sufficient, yielding an
interpretable model.

Figure 1. Shape plots for numerical features from models learned on the COMPAS
dataset. The GAM forest model is shown in blue dots, while its distillation
into the two-segment curve model is shown in orange lines (both map to the
left Y axis). Cumulative distribution functions of the corresponding signals
are shown in grey (right Y axis).

In the classic GAM models, splines are used as shape functions (Hastie and
Tibshirani, 1986). Another commonly used shape function is piecewise-linear
functions (Wikipedia, 2020f). These representations contain a small number of
variables (_e.g_. knots) and thus are concise and human-readable. However,
directly optimizing such representations often yields less accurate models
than alternative model representations. For example, Lou _et al_. (Lou et al.,
2012) showed that learning spline GAMs is less accurate than learning bagged
boosted decision forest GAMs. Our experiments show similar results for
directly optimizing GAMs composed of piecewise-linear curves using Stochastic
Gradient Descent (SGD) methods. Broadly speaking, the model representations
using decision forest GAMs have the advantage during model optimization, but
the resultant models are not human-readable. This is the case even when there
exists a simpler model with a concise, human-readable form that provides
comparable accuracy.

Inspired by the model distillation work in which relatively small decision
forests or neural networks can be distilled from much larger ensembles, but
not trained directly from data, to match the accuracy of complex models
(Buciluundefined et al., 2006; Hinton et al., 2015), we propose to distill
interpretable models into readable representations in a separate process after
model optimization. This decouples the initial, learned model representation
from the final, published model representation. For example, the proposed
distillation methodology can be applied to additive models trained using
bagged boosted decision trees (Lou et al., 2012), as well as additive neural
nets (Agarwal et al., 2020; Zhuang et al., 2021).

In this paper, we describe a technique for distilling models composed of
univariate components into human readable representations, in particular, the
piecewise-linear curves described in Section 2.2. The output of our
distillation technique is illustrated in Listing 1 and Figure 1, which show
textual and graphical representations of piecewise-linear curves obtained by
applying our approach to a decision forest GAM trained on the COMPAS dataset
(described in Section 2.1). The distilled model is a concise representation of
the decision forest GAM model and is converted to human-readable source code.

Listing 1: Code for a distilled COMPAS model

⬇

score = sum([

PWLCurve("age", [(18, 3.13), (21, 0.5914),

(46, -0.7206)], fx="log"),

PWLCurve("priors_count", [(0, -0.8415), (1, -0.4452), (38, 2.146)],
fx="log1p"),

PWLCurve("length_of_stay", [(0, -0.1855),

(3, -0.04099), (4, 0.2443)], fx="log1p"),

EnumCurve("c_charge_degree", {1: 0.0198, 2: -0.0384}),

## ... other features ...

])

From here on, we will use ”curves” to refer to piecewise-linear curves, ”curve
models” to refer to models where each component is a curve, and ”code” to
refer to the textual representations of curve models or curves.

The rest of this paper is structured as follows. After presenting the
preliminaries in Section 2, we elaborate on the benefits of using curve models
in Section 3. We then describe the localized distillation process in Section 4
and piecewise-linear approximation algorithm, sometimes referred to as
segmented regression (Wikipedia, 2020f), for creating curve models in Section
5. Lastly, we present experimental results on for datasets: COMPAS, FICO,
MSLR-WEB30K, and CWS in Section 6 and conclude the paper in Section 7.

##  2\. Preliminaries

Throughout the paper, we will use the data sets used in this paper as concrete
examples to explain our methods. Thus, we first describe them in this section.
We also give the formal definition of piecewise-linear-curves in this section.

###  2.1. Data Sets

We used the following four datasets to represent different settings:
classification, regression, and ranking. The first three are publicly
available.

  * •

The COMPAS dataset111https://github.com/propublica/compas-analysis is the
result of a ProPublica investigation (Angwin et al., 2016) into possible
racial bias of the proprietary COMPAS model score for defendants in Broward
county, Florida. The dataset has been studied extensively in the context of
bias, fairness, and interpretability (Tan et al., 2018; Dressel and Farid,
2018; Kleinberg, 2018; Chouldechova, 2017). Labels are binary and indicate
whether recidivism occurred for an individual within a time period. We use
area under the receiver operating characteristic curve (AUC-ROC) to measure
classifier accuracy. COMPAS has 6 features and four of them are used as
examples in this paper: age, priors_count, length_of_stay, and
c_charge_degree.

  * •

The FICO dataset (FIC, 2018) is composed of real-world anonymized credit
applications along with risk scores. Labels are a risk score for an
individual. We use root mean square error (RMSE) to measure regressor
accuracy. FICO has 24 features and we use two features as examples in our
paper: MSinceMostRecentDelq, Months Since Most Recent Delinquency;
PercentTradesWBalance, Percent Trades with Balance.

  * •

The MSLR-WEB30K dataset (Qin and Liu, 2013) is a widely used learning-to-rank
benchmark dataset. Labels are per document relevance judgements. We use
normalized discounted cumulative gain at $k=5$ (NDCG@5) to measure ranker
accuracy. MSLR-WEB30K is significantly larger both in number of features (136)
and number of training examples (~2,000,000 per cross validation fold). We use
it to compare our curve approximation algorithm to pwlf (Jekel and Venter,
2019), a publicly available alternative, on the basis of accuracy, robustness
and efficiency. We use two features as examples in our paper: feature_0011,
Body stream length; feature_0128, Inlink number.

  * •

The Chrome Web Store (CWS) dataset is a private and anonymized dataset
originating from Chrome Web Store logs. Each query corresponds to a visit to
the Chrome Web Store. The items within each query were the ones shown to the
user. Labels correspond to user actions such as clicking, installing, or no
action whatsoever. We again use NDCG@5 to measure ranker accuracy. A similar,
but distinct datset from the Chrome Web Store was previously studied by Zhaung
_et al_. (Zhuang et al., 2021). Unlike in that previous work, in this instance
we do not utilize query level ”context” features, instead using only 14 item
level features. The queries are also distinct.

In each case, we distill a decision forest GAM and evaluate the accuracy of
the distilled curve models. The COMPAS and FICO datasets represent high-stakes
domains (Rudin, 2018) in which the benefits of curve models, discussed below,
are particularly compelling. FICO, MSLR-WEB30K, and CWS have been previously
studied in the context of interpretability (Agarwal et al., 2020; Zhuang et
al., 2021; Lou et al., 2013; Chen et al., 2018). Furthermore, the results from
MSLR-WEB30K demonstrate that the accuracy of this approach is not limited to
small datasets.

###  2.2. Piecewise-Linear Curves

A piecewise linear curve (PWLCurve) is defined by a list of control points
$S=[(x_{k},y_{k})]_{k=1}^{K}$ through which the curve must pass. Between
control points, output $y$ values are determined by performing linear
interpolation between neighboring control points. Beyond the leftmost or
rightmost control points, output values are capped to the $y_{k}$-value of the
neighboring control point. More formally, assuming $x_{k}$’s are ordered,
_i.e_. $x_{k}<x_{k+1}$, the definition of a piecewise linear curve can be
described as:

$PWL(x;S)=\begin{cases}y_{1}&\text{if }x<x_{1},\\\
\frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}}(x-x_{k})+y_{k}&\text{if }x_{k}\leq x\leq
x_{k+1},\\\ y_{K}&\text{if }x>x_{K}.\end{cases}$

In most cases of interest 5 or 6 control points, defining 4 or 5 interior
segments, is sufficient to capture the desired behavior.

We allow for an optional $x$-transformation, specified with the fx argument,
to fit curves to data with different scales. When an $x$-transformation is
present it is applied to the input value and $x$-values of all the control
points, and then linear interpolation is performed in the transformed space.
We support identity (default), log, log1p and symlog1p transformations. Here
symlog1p is defined as sgn(x) * log1p(abs(x)) and is suitable for highly-
variable features that take on both positive and negative values.

Univariate categorical functions are represented by EnumCurve, which directly
maps input values to outputs using a discrete mapping.

##  3\. Background & Motivation

Interpretable models are critical for high-stakes decisions (Rudin, 2018) and
provide many advantages over more complex model structures (Caruana et al.,
2015; Du et al., 2019). In this section we explain how distilling
interpretable models into curve models reinforces these benefits and addresses
a variety of real-world engineering challenges. Here, one underlying theme is
that distilling models into human-readable source code _reduces a novel
machine learning problem to an established software engineering problem with
an abundance of existing solutions_.

###  3.1. Greater Transparency

A model is transparent if it provides a textual or graphical representation
that enables its behavior to be understood comprehensively (Ustun and Rudin,
2014). One way in which the proposed method provides greater transparency is
by simplifying graphical depictions of a model while retaining its essential
characteristics. It is often argued, implicitly or explicitly, that the shape
plots of an interpretable model are an _exact description_ of the model and
therefore provide a reliable way to understand the model. While this claim is
narrowly true, it is misleading in general. Unless given specific guidance,
humans will naturally discount certain fine-grained details of the plots when
developing an understanding of the model. By distilling interpretable models
to a concise representation, we discard extraneous characteristics and reduce
the mental effort necessary to understand the model. For example, it is not
immediately obvious what understanding an individual should derive from the
shape plots of the feature_0011 (body stream length), and feature_0128 (inlink
number) features in the initially-learned MSLR-WEB30K model, shown in Figure
2. Indeed, different individuals may derive qualitatively different
understandings from these graphical depictions. However, given the additional
knowledge that the distilled curve model represented by the overlaid curves in
Figure 2 has nearly identical accuracy, an observer can make much stronger
inferences about the model’s essential characteristics. Interpretability can
be increased even further by imposing monotonicity constraints. We discuss the
effect of such constraints in Section 6.4.

Figure 2. Shape plots for a GAM forest model (in blue dots) and 5 segment
curve distillation (in orange lines) for the MSLR-WEB30K dataset.

Clearly when distillation yields a simpler model with comparable accuracy we
would say the distillation process has succeeded. However, instances where
distillation yields a model with inferior accuracy warrant further
investigation because the apparent ”failure” can often be attributed to
essential characteristics of the teacher model that were not successfully
transferred to the student _because they violate a prescribed notion of human-
interpretability_. We examine one representative example of this phenomenon in
Section 6.2. While a complete discussion of this principle is beyond the scope
of this paper, we note that the idea can be viewed as an extension of the use
of structural constraints to define “interpretable” models, just now applied
to the structure of individual functions in the model. Under this policy, if
the accuracy of a candidate model cannot be reproduced using a predefined
class of expressive, “human-scale” functions (_e.g_. curves with a small
number of truncated control points) its transparency would be called into
question.

###  3.2. Constructive Regularization

The proposed method can also be viewed as a post-hoc regularization process
that is completely compatible with, and complementary to, optimization-based
regularization techniques (_e.g_. L1/L2 penalties or monotonicity
constraints). In the context of regularization, our emphasis on conciseness is
aligned with the minimum description length principle (Wikipedia, 2020d) for
model selection. Ustun and Rudin (Ustun and Rudin, 2014) applied similar
reasoning to motivate linear models with small, integer-valued weights. The
constrained description length of curves provides limited capacity for
capturing idiosyncratic behavior. As a result, curve distillation successfully
removes aberrations from teacher model functions. This regularization effect
can be seen in Figure 1 and Figure 2. The fewer segments the greater the
effect. To find the most concise curve model we can repeatedly apply the
proposed method with decreasing number of control points. Naturally, the
optimality of this approach is subject to the limitations of our localized
distillation methodology (see Section 4) and curve approximation algorithm
(see Section 5). While it is difficult to directly compare models with
different functional representations, comparing the length and readability of
their corresponding code is instructive.

One practical advantage of curve-based regularization is that regularity is
enforced by construction and the complexity of individual curves is readily
apparent and quantifiable. Therefore, organizations that adopt curve models
can set objective guidelines about model complexity that developers can
anticipate when submitting model candidates for approval. Such guidelines can
specify the maximum number of curve segments, maximum number of significant
digits per curve control point, or monotonicity of the curve. Similar to the
use of nothing-up-my-sleeve numbers in cryptography (Wikipedia, 2020e), curve
models enable developers to preemptively address suspicions about potential
weaknesses and constructively prove the robustness of a given model candidate.
In general, standardizing development around curve models is a straightforward
way for organizations to systematically enforce best practices, defend against
common mistakes and pitfalls, and expedite model verification and approval.
The accessible, readable nature of curve models enables organization members
beyond engineers (_e.g_. executives, product managers, etc.) to participate in
this approval process.

###  3.3. Readable, Editable Code

Curve model code can be read, reviewed, merged, and versioned like
conventional source code. An example model for the COMPAS dataset is shown in
Listing 1. One can understand how a curve model would behave under novel or
extremal conditions by mentally “evaluating” the model under hypothetical
“what if?” scenarios without the need for additional tools. Subjecting models
to a traditional source code review process facilitates a more rigorous
examination of the model’s characteristics and greater accountability than is
possible with non-readable models. Indeed, conducting “model review” through
source code review ensures that the candidate model itself - not some
separate, potentially inconsistent description or artifact of the model or how
it was trained - is the subject of review. In the event that undesirable model
behavior is discovered, the model’s code may be directly edited to correct
such issues. For example, in the case of the COMPAS model a user may wish to
deliberately cap the contribution of features such as priors_count and
length_of_stay features for legitimate policy reasons not captured by
classification metrics such as AUC-ROC. The contribution of other features can
be entirely removed. Agarwal _et al_. (Agarwal et al., 2020) discussed how
such an approach of training with biased features and then removing them can
potentially be better than simply training without biased features. This
approach can prevent the model from extracting bias through other features
which are correlated with biased ones.

Model transparency is essential in the context of high-stakes decisions
(Rudin, 2018) arising in criminal justice, finance, health care, and other
areas. Providing the complete source of the model in simple, portable, human-
readable code makes the models transparent. Compared to human-readable models
produced by CORELS (Angelino et al., 2017), which are expressed in
universally-understandable if-then language, curve models sacrifice
accessibility for greater expressiveness and general-purpose application.

###  3.4. Collaborative Model Development

Curve distillation is compatible with any algorithm or modeling technique that
results in univariate functions. In the experiments section we apply the
proposed technique to decision forest GAMs on several datasets. Previous work
(Zhuang et al., 2021) applied the proposed technique to GAMs learned via
neural networks, as well as similar neural networks with limited interactions
via multiplicative pairs. Organizing collaborative development around curve
models enables engineers to apply a plurality of different tools, techniques,
or platforms to optimize components of a (potentially large-scale) model.
Engineers are free to choose a modeling approach that maximizes their
productivity, similarly to how engineers use multiple IDEs, code formatters,
or linters to collaboratively develop software. Curve distillation can be
viewed as a “format conversion” tool that translates an arbitrary and
potentially exotic model representation into a fixed, agreed-upon vocabulary
of human-readable building blocks.

###  3.5. Straightforward Deployment

Curve models are fast-to-evaluate and straightforward to deploy. Since
evaluation requires minimal computation - just a handful of floating point
operations per curve - curve models are well-suited for performance-critical
applications. Curves are a portable, platform-agnostic representation that can
be natively supported in a variety of languages or systems with little effort.
For example, Listing 2 shows a C++ implementation of a COMPAS model with 2
segments. In general, curve models are straightforward to deploy because they
offer a multitude of integration options. Curves can be embedded in
configuration files, passed via CGI parameters, manually embedded into complex
applications in a piecemeal fashion, systematically translated to a target
representation, or evaluated by existing runtime systems with a few
incremental extensions.

Listing 2: A COMPAS model as a C++ function

⬇

double COMPAS(double age, double priors_count,

double length_of_stay, int charge_degree,

// ... other features ...

) {

static auto age_curve = PWLCurve({{18, 3.13},

{21, 0.5914}, {46, -0.7206}}, "log");

static auto priors_count_curve = PWLCurve(

{{0, -0.8415}, {1, -0.4452}, {38, 2.146}},"log1p");

static auto length_of_stay_curve = PWLCurve(

{{0, -0.1855}, {3, -0.04099}, {4, 0.2443}}, "log1p");

static auto charge_degree_curve = EnumCurve({{1, 0.0198}, {2, -0.0384}});

// ... other features ...

return (age_curve.Eval(age) +

priors_count_curve.Eval(priors_count) +

length_of_stay_curve.Eval(length_of_stay) +

charge_degree_curve.Eval(charge_degree) +

// ... other features ...

);

}

##  4\. Localized Distillation

Our distillation process takes two inputs: a teacher model containing one or
more univariate functions, and a representative dataset (generally the
training data). Our method differs from conventional distillation techniques
in that we (1) distill each univariate function in isolation and (2) optimize
for mean squared error (MSE) when approximating each univariate function.
Specifically, each univariate function in the teacher model is evaluated on
the dataset to produce representative $(x,y)$ example pairs. For discrete
categorical features we create a mapping where each unique $x$ is mapped to
the mean $y$. For numerical features, we produce a PWLCurve using the
approximation algorithm described in Section 5. If the teacher model contains
univariate functions nested within other univariate functions we replace the
source functions in a bottom-up fashion. Otherwise, all non-nested functions
can be approximated in parallel. The final model is constructed by replacing
each original univariate function with its PWLCurve approximation.

Conventionally, model distillation involves a global optimization using the
same (or at least similar) objective to the original teacher model training.
This objective may differ from a point-wise MSE objective. For example,
ranking objectives often have pair-wise definitions. Why then do we advocate a
localized optimization using a MSE objective in all circumstances? The primary
answer is that, in the context of interpretable models, there is substantial
value in maintaining a strong one-for-one correspondence between each source
function and target function. Notably, this allows us to visualize each shape
function in the teacher model against its corresponding curve replacement.
Additionally, we can attribute distillation failures - data instances where
the curve model is less accurate than the teacher model - to specific
univariate functions, and to take remedial actions. For example, in the Figure
5 we can immediately tell that the shape function of $x_{1}$ was not well-
approximated by a curve. In the experiments section we show that the
meaningful behavior of nearly all shape functions can be accurately captured
by curves with three to five segments. Furthermore, when the meaningful
behavior is not captured, it is generally due to inherently non-interpretable
behavior being lost.

While a global optimization approach (_i.e_. optimizing the parameters of all
curves in the target model simultaneously) using a problem-specific metric
might produce a more accurate result, it is computationally more expensive and
would lack the same one-to-one correspondence with the teacher model, making
distillation failures more difficult to diagnose. Furthermore, if higher
accuracy is desired, the output of the proposed distillation process can be
used to initialize a global optimization of the curve model’s parameters.

##  5\. Piecewise-Linear Curve Approximation

Given a univariate numerical function $f(x)\rightarrow y$, our goal is to
produce a PWLCurve $c(x)\rightarrow y$ that faithfully approximates $f(x)$ by
minimizing the $MSE(c(x),f(x))$ over sample data. Clearly, the accuracy of the
overall distillation method depends critically on the accuracy of the
individual curve approximations - _i.e_. how much metric loss is incurred when
each $c(x)$ is substituted for the corresponding $f(x)$ in the trained model.

Additionally, the practical success of the methodology also depends on the
robustness and efficiency of the approximation algorithm. To enable systematic
use of curve distillation in model training pipelines, the approximation
algorithm must run with minimal configuration. Complex hyperparameters pose a
significant barrier to entry. We have designed pwlfit, our piecewise linear
approximation algorithm, so that in practice users only need to consider the
num_segments and mono (monotonicity) parameters. While num_segments=5 segments
and mono=False is sufficient to get high accuracy (as demonstrated by our
experiments), it is desirable to investigate whether the model can be further
simplified with fewer segments or with monotonicity restrictions. To
facilitate such investigations it is important that distillation runs quickly
(_e.g_. less than 1 second per function) which enables interactive analysis
via Jupyter notebooks (Kluyver et al., 2016) or other tools. These practical
considerations have informed various decisions in the design of pwlfit. In
particular, we prefer an algorithm which quickly and reliably yields high
accuracy results with minimal configuration to one which sacrifices either of
these practical considerations for marginal gains in accuracy.

In this section we will describe the salient characteristics and noteworthy
features of pwlfit. We invite interested readers to consult the publicly-
available source code of pwlfit (Sterling and Ravina, 2019), for additional
details.

###  5.1. Algorithm

Given a list of $(x,y,weight)$ points and a desired number of segments $k$, we
search for a PWLCurve to minimize mean squared error, MSE. A PWLCurve with $k$
segments is characterized by its $k+1$ control points – a set of $x$-knots and
their corresponding $y$-knots. Given only the $x$-knots, we can solve a linear
least squares expression for the optimal $y$-knots and the resulting error.
Since we don’t know the correct $x$-knots, we search through the space of
possible $x$-knots and solve a least squares expression at each step to
calculate the error.222pwlf (Jekel and Venter, 2019) implements a similar
approach. We will compare with it in our experiments.

####  5.1.1. Initial Downsampling

For performance, we randomly downsample large datasets to approximately one
million points before fitting. We downsample to reduce the cost of sorting,
which dominates the runtime for large data. This downsampling imposes a
negligible quality loss. To further reduce runtime, we discretize the search
space for $x$-knots. We choose num_samples $x$-values from the data, spaced
equally by cumulative weight, and search over the combinations of $x$-knots
from that sampled set of candidates. Using the default 100 samples, our
candidates are the $x$-values at $(0\%,1.01\%,\dots,98.9\%,100\%)$ of the
cumulative weight.

####  5.1.2. Knot Discretization

For data with many repeated $x$-values, some of our candidates will be
duplicates. For example, $55\%$ of the values in the length_of_stay feature in
the COMPAS data set are 0 or 1. In such cases, we iteratively resample at
higher rates (such as $0\%,0.505\%,1.01\%$, etc.) until we collect a suitable
number of distinct candidates, never exceeding the specified num_samples
parameter.

####  5.1.3. Condensation

To minimize the cost of each linear least squares step, we condense the data
using a novel technique described in Appendix B. Given num_samples candidate
knots, we condense the full data into two synthetic points per adjacent pair
of candidates, for a total of 2 * (num_samples - 1) synthetic points. For any
function that’s linear between each adjacent pair of candidate $x$-knots,
which is guaranteed by our choice of discrete candidate $x$-knots, these
condensed points perfectly recreate the loss of that function over the full
data set. We run our linear least squares solver on the condensed points
instead of the full data set, which reduces our cost per solve from
${\mathcal{O}(\text{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\footnotesize{\@listingGroup{ltx_lst_identifier}{num\textunderscore
points}}}}}})$ to
${\mathcal{O}(\text{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\footnotesize{\@listingGroup{ltx_lst_keyword}{\color[rgb]{0.0,0.27,0.13}{n}um\textunderscore
samples}}}}}})$. This is purely a performance optimization, with no quality
implications.

Figure 3. Candidate $x$-knots (red vertical lines) and derived condensed
points (pink large dots) on the age piece of a COMPAS GAM forest model (blue
dots). For visual clarity, this illustration considers only five $x$-knot
candidates.

####  5.1.4. Global Optimization via Greedy Search

After discretization, the solution space consists of
${{\text{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\footnotesize{\@listingGroup{ltx_lst_keyword}{\color[rgb]{0.0,0.27,0.13}{n}um\textunderscore
samples}}}}}}\choose{\text{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\footnotesize{\@listingGroup{ltx_lst_keyword}{\color[rgb]{0.0,0.27,0.13}{n}um\textunderscore
segments}}}}}}+1}$ $x$-knot combinations, which is still too large for an
exhaustive search. To make the search tractable we use a greedy search
heuristic that optimizes one $x$-knot at a time. Specifically, at each step of
the process we evaluate the error associated with each candidate $x$-knot, and
keep the candidate that yields the least error.

With this approach, we optimize in two stages. We begin with a single $x$-knot
as our solution, and greedily add the best remaining candidate $x$-knot until
our solution consists of (num_segments + 1) $x$-knots. Then we cycle through
our solution, removing one $x$-knot at a time and replacing that $x$-knot with
the best remaining candidate $x$-knot, which could be the same $x$-knot that
we just removed. We continue this cycle of iterative improvements until our
solution converges, or until we’ve exceeded the maximum number of iterations
(defaulting to 10 iterations).

####  5.1.5. Slope Constraints & Monotonicity

pwlfit can impose a minimum and/or maximum slope on the solution via bounded
least squares. Instead of solving the least squares expression directly for
the $y$-knots, we solve it for the deltas between adjacent $y$-knots. Then we
impose a min/max slope by bounding the deltas. Slope restrictions can be used
to limit the spikiness of curves, but we primarily use them to impose
monotonicity. For example, specifying min_slope=0 restricts to monotonically
non-decreasing functions while max_slope=0 restricts to monotonically non-
increasing functions. Specifying a min_slope greater than 0 or a max_slope
less than 0 restricts to strictly increasing or decreasing functions,
respectively.

pwlfit can deduce the direction of monotonicity by applying isotonic
regression (Wikipedia, 2020b) to the condensed points. We fit an increasing
and a decreasing isotonic regression, and use the direction that minimizes
mean squared error. The user can override this behavior by specifying the
direction explicitly or by disabling monotonicity entirely.

####  5.1.6. Input Transformations

pwlfit can also interpolate in a transformed $x$-coordinate space instead of
the original space, as a simple form of feature engineering. pwlfit transforms
the $x$-values before learning the curve. Specifically, pwlfit will choose a
candidate $x$-transformation, fx, among log, log1p, or symlog1p based on the
range of the $x$-values and then proceed with that transformation if it
increases the Pearson correlation between fx and $y$ by a noticeable amount
over the identity transformation. Alternatively, the user can specify any
strictly increasing 1D transform or specify the identity transform to disable
transformation.

##  6\. Experiments

###  6.1. Distillation Accuracy

Table 1 and Figure 4(a) show the results obtained from experiments on the
different datasets. A complete set of results can be found in Table 2 in
Appendix A. The results of applying our distillation technique with our
piecewise-linear approximation algorithm are presented as pwlfit. We present
results from using various numbers of segments with and without a monotonicity
restriction and otherwise default parameters. In all cases we truncated the
control points to four significant digits. We also present several additional
reference points to provide context.

  * •

SGD: We directly learn the curves with the Adadelta(Zeiler, 2012) optimizer.
We initialize the $y$ values of the control points as zeros. For the $x$
values of the control points we use the quantiles for numerical features
(_e.g_. 0%, 50%, 100% for a three point, two segment curve) or all unique
values for categorical features. We then apply Adadelta to optimize the $y$
values. Simultaneously optimizing $x$ and $y$ values was also attempted, but
the results were always worse than optimizing $y$ values alone.

  * •

NAM: Neural Additive Models (NAMs) (Agarwal et al., 2020) is another method
for learning interpretable models proposed by Agarwal _et al_. We present
their result for reference where applicable.

  * •

Interacting forest: We train a bagged, boosted decision forest allowing
feature interactions to demonstrate the accuracy of a non-interpretable, high-
complexity ”black box” model.

  * •

GAM forest: We train a bagged boosted decision forest GAM by restricting each
tree to use only one feature. This model is also the source model for our
distillation technique.

  * •

pwlf: We apply our distillation technique using an alternative piecewise-
linear approximation algorithm pwlf(Jekel and Venter, 2019).

On each dataset we used five fold cross validation and present the metric mean
and sample standard deviation across folds. We used three different metrics to
evaluate accuracy: AUC-ROC, RMSE, and NDCG@5 for the three different tasks of
classification, regression, and ranking. Further details on our
experimentation setup can be found in Appendix A and further details on the
datasets, labels, and metrics can be found in Preliminaries 2.1.

Table 1. Metrics across datasets. For AUC-ROC and NDCG@5 higher is better, and for RMSE lower is better. MSLR-WEB30K and CWS were not used by the NAM paper and are omitted from that row. Metric values are the mean from five fold cross validation $\pm$ the sample standard deviation. Model | COMPAS (AUC-ROC) | FICO (RMSE) | MSLR-WEB30K (NDCG@5) | CWS (NDCG@5)
---|---|---|---|---
Interacting forest | $0.742\pm 0.011$ | $3.128\pm 0.092$ | $0.485\pm 0.002$ | $0.461\pm 0.003$
GAM forest | $0.741\pm 0.013$ | $3.495\pm 0.109$ | $0.442\pm 0.002$ | $0.460\pm 0.004$
NAM | $0.741\pm 0.009$ | $3.490\pm 0.081$ |  |
pwlfit num_segments=5, mono=False | $0.743\pm 0.012$ | $3.494\pm 0.096$ | $0.441\pm 0.002$ | $0.454\pm 0.002$
pwlfit num_segments=5, mono=True | $0.743\pm 0.013$ | $3.693\pm 0.101$ | $0.437\pm 0.003$ | $0.452\pm 0.003$
pwlf num_segments=5 | $0.743\pm 0.012$ | $3.503\pm 0.096$ | $0.433\pm 0.003$ | $0.454\pm 0.003$
SGD num_segments=5 | $0.741\pm 0.010$ | $3.643\pm 0.097$ | $0.405\pm 0.002$ | $0.448\pm 0.004$
SGD num_segments=20 | $0.738\pm 0.011$ | $3.499\pm 0.117$ | $0.419\pm 0.003$ | $0.455\pm 0.003$

(a) Top level metrics across datasets. For AUC-ROC and NDCG@5 higher is
better, for RMSE lower is better.

(b) Per fit metrics for MSLR-WEB30K Fold 1. Note that for 4 segments pwlf has
extreme outliers for RMSE vs the source submodel.

Figure 4. Comparisons of different methods across the 4 datasets.

Our results show that applying our distillation technique with 4-5 segments
with pwlfit produces models which are as accurate as both the source GAM
forest and NAM models for all datasets except CWS where a small gap remains.
We investigate this accuracy gap in detail in Section 6.2 below. In the case
of the COMPAS dataset these models are as accurate as full complexity models.
Applying our technique with pwlf produces competitive results, albeit less
accurate on the MSLR-WEB30K dataset. By contrast, the results show that
learning curves directly via SGD is less general. On the FICO and CWS datasets
more segments are required to achieve accuracy comparable to the GAM forest
models. On the MSLR-WEB30K dataset the accuracy is inferior even with many
more segments.

The consistent accuracy of applying our distillation approach with pwlfit on
these four datasets and three separate tasks (classification, regression,
learning to rank) demonstrates that the process is not sensitive to either the
specific data or the top level objective being used.

###  6.2. Distillation Failures

In Section 3.1 we explained how distillation yielding a model with inferior
accuracy warrants further investigation because the purported ”failure” can
often be attributed to essential yet non-interpretable characteristics of the
teacher model not transferring to the student model. The accuracy gap observed
on the CWS dataset is an example of this phenomenon. Figure 5 shows the worst
two fits from the CWS dataset. The plots have been redacted to maintain the
privacy of the dataset. For each plot it is clear that the original teacher
submodel had some non-interpretable behavior which was lost during
distillation. This is most evident for feature $x_{1}$, the worst offender,
where the output is highly erratic. If the original teacher submodel is not
distilled for these two features then the accuracy gap between the original
teacher model and 5 segment non-monotonic distillation drops from 0.0059 to
0.003 (_i.e_. ~50% of the gap is recovered).

To identify the above two failures we applied the following method.

  * •

Begin with the original teacher model. For each submodel compute the metric
delta against the teacher model from distilling only that submodel and no
others.

  * •

Perform the above on each cross validation fold using the validation set and
average the metric deltas across folds.

  * •

Sort the features by their associated metric delta to determine the worst
distillations.

Figure 5. The worst curve distillations from the CWS dataset using 5 segments.

###  6.3. Efficiency & Robustness

The experiments of the previous section showed that pwlfit more accurately
distills the source model across datsets than pwlf. We also found on the MSLR-
WEB30K dataset that pwlfit is more efficient and robust than pwlf. Figure 4(b)
shows per fit metrics from the first fold of the MSLR-WEB30K dataset as the
number of segments varies without monotonicity. The top plot shows the time in
seconds, as measured on a ThinkStation P520, to fit each of the 136 submodels
of the source GAM forest. We find that pwlfit is faster in the average case as
the number of segments increases, and has a narrower distribution. The bottom
plot shows the RMSE of each fit against the 136 submodels of the source GAM
forest. We again find that pwlfit performs favorably in the average case with
a narrower runtime distribution.

It’s worth noting that pwlf by default does not perform any downsampling. For
the MSLR-WEB30K dataset running pwlf without any downsampling was
prohibitively expensive. For all of our experiments we ran pwlf with a pre-
processing downsample to a random 1000 examples. We found this to be a fair
point for balancing speed and quality when comparing to pwlfit. It is of
course possible with both algorithms to modify the number of samples used to
strike a different trade-off between run time and accuracy.

###  6.4. Monotonicity

As discussed in Section 5, pwlfit can fit monotonic curves with automatic
direction detection. Figure 4(a) compares curve models fit with and without
monotonicity constraints (automatically inferring the direction) across
datasets. For the COMPAS dataset monotonic and non monotonic models are
comparably accurate, while for FICO, MSLR-WEB30K, and CWS, non-monotonic
models are more accurate.

Monotonicity with respect to appropriate features is desirable for
interpretable models. In these cases a monotonic model may be preferable to a
non-monotonic one, even if it is less accurate. For example, Figure 6 compares
monotonic and non-monotonic 5 segment curve models on the FICO dataset for the
MSinceMostRecentDelq, and PercentTradesWBalance features. Given the semantic
meaning of these features, it is desirable from a transparency and incentives
standpoint for the model output to be monotonic with respect to each of them.

Figure 6. Curve distillations on the FICO dataset using 5 segments without
monotonicity (orange) and with monotonicity (green).

##  7\. Conclusion

We have introduced a novel method for distilling interpretable models into
human-readable code using piecewise-linear curves and demonstrated its
efficacy on four datasets. We have shown that curve models match or outperform
the accuracy achieved by other additive models. On smaller datasets, curve
models match the accuracy of more complex models, like interacting decision
forests. Our localized distillation methodology is applicable to any model
containing univariate numerical functions and is straightforward to implement
using the publicly-available pwlfit(Sterling and Ravina, 2019) library.

We have explained how curve model distillation reinforces interpretability and
addresses a variety of real-world engineering challenges. Curve models are 1)
transparent, 2) well-regularized, 3) easy to analyze for presence of biases or
other fairness issues, and 4) can be directly edited or improved outside the
context of machine learning to fix the aforementioned fairness issues.
Distilling models into human-readable code allows one to address novel machine
learning problems using well-established software engineering methods. Curve
models can be improved by multiple contributors in parallel, reviewed, and
made to systematically follow best practices. Curve models are well-suited for
production applications, since they can be natively supported in many
languages, are easy to deploy, and fast to evaluate.

###### Acknowledgements.

We thank Vytenis Sakenas, Jaime Fernandez del Rio, Benoit Zhong, and Petr
Mitrichev for their support and providing the algorithms and optimization
infrastructure used in our experiments. We also thank Paul Heymann, Diego
Federici, Mike Bendersky, Paul Haahr, and Petr Mitrichev for their helpful
feedback and detailed reviews. Lastly, we thank Xinyu Qian, Janelle Lee, Po
Hu, and Chary Chen for preparing the CWS data set for our experiments.

## References

  * (1)
  * FIC (2018) 2018\.  FICO Explainable Machine Learning Challenge.  https://community.fico.com/s/explainable-machine-learning-challenge.
  * Agarwal et al. (2020) Rishabh Agarwal, Nicholas Frosst, Xuezhou Zhang, Rich Caruana, and Geoffrey E. Hinton. 2020.  Neural Additive Models: Interpretable Machine Learning with Neural Nets.  arXiv:cs.LG/2004.13912
  * Angelino et al. (2017) Elaine Angelino, Nicholas Larus-Stone, Daniel Alabi, Margo Seltzer, and Cynthia Rudin. 2017\.  Learning Certifiably Optimal Rule Lists. In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ _(KDD ’17)_.
  * Angwin et al. (2016) Julia Angwin, Jeff Larson, Surya Mattu, and Lauren Kirchner. 2016.  Machine bias: There’s software used across the country to predict future criminals.  _And it’s biased against blacks. ProPublica_ 23 (2016).
  * Buciluundefined et al. (2006) Cristian Buciluundefined, Rich Caruana, and Alexandru Niculescu-Mizil. 2006.  Model Compression. In _Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ _(KDD ’06)_. 535–541.
  * Caruana et al. (2015) Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Marc Sturm, and Noemie Elhadad. 2015\.  Intelligible Models for HealthCare: Predicting Pneumonia Risk and Hospital 30-Day Readmission. In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ _(KDD ’15)_. 1721–1730.
  * Chen et al. (2018) Chaofan Chen, Kangcheng Lin, Cynthia Rudin, Yaron Shaposhnik, Sijia Wang, and Tong Wang. 2018\.  An Interpretable Model with Globally Consistent Explanations for Credit Risk.  arXiv:cs.LG/1811.12615
  * Chouldechova (2017) Alexandra Chouldechova. 2017\.  Fair Prediction with Disparate Impact: A Study of Bias in Recidivism Prediction Instruments.  _Big Data_ 5, 2 (Jun 2017), 153–163.  https://doi.org/10.1089/big.2016.0047
  * Dressel and Farid (2018) Julia Dressel and Hany Farid. 2018.  The accuracy, fairness, and limits of predicting recidivism.  _Science Advances_ 4, 1 (2018).  https://doi.org/10.1126/sciadv.aao5580
  * Du et al. (2019) Mengnan Du, Ninghao Liu, and Xia Hu. 2019.  Techniques for interpretable machine learning.  _Commun. ACM_ 63, 1 (Dec 2019), 68–77.  https://doi.org/10.1145/3359786
  * Hastie and Tibshirani (1986) Trevor Hastie and Robert Tibshirani. 1986.  Generalized Additive Models.  _Statist. Sci._ 1, 3 (08 1986), 297–310.
  * Hinton et al. (2015) Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. 2015\.  Distilling the Knowledge in a Neural Network.  arXiv:stat.ML/1503.02531
  * Jekel and Venter (2019) Charles F. Jekel and Gerhard Venter. 2019.  _pwlf: A Python Library for Fitting 1D Continuous Piecewise Linear Functions_.  https://github.com/cjekel/piecewise_linear_fit_py
  * Kleinberg (2018) Jon Kleinberg. 2018\.  Inherent Trade-Offs in Algorithmic Fairness. In _Abstracts of the 2018 ACM International Conference on Measurement and Modeling of Computer Systems_ _(SIGMETRICS ’18)_.
  * Kluyver et al. (2016) Thomas Kluyver, Benjamin Ragan-Kelley, Fernando Pérez, Brian Granger, Matthias Bussonnier, Jonathan Frederic, Kyle Kelley, Jessica Hamrick, Jason Grout, Sylvain Corlay, Paul Ivanov, Damián Avila, Safia Abdalla, and Carol Willing. 2016\.  Jupyter Notebooks – a publishing format for reproducible computational workflows. In _Positioning and Power in Academic Publishing: Players, Agents and Agendas_ , F. Loizides and B. Schmidt (Eds.). IOS Press, 87 – 90.
  * Kolmogorov (1957) A. K. Kolmogorov. 1957\.  On the Representation of Continuous Functions of Several Variables by Superposition of Continuous Functions of One Variable and Addition.  _Doklady Akademii Nauk SSSR_ 114 (1957), 369–373.
  * Lou et al. (2012) Yin Lou, Rich Caruana, and Johannes Gehrke. 2012.  Intelligible Models for Classification and Regression. In _Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ _(KDD ’12)_. 150–158.
  * Lou et al. (2013) Yin Lou, Rich Caruana, Johannes Gehrke, and Giles Hooker. 2013\.  Accurate Intelligible Models with Pairwise Interactions. In _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ _(KDD ’13)_. 623–631.
  * Qin and Liu (2013) Tao Qin and Tie-Yan Liu. 2013.  Introducing LETOR 4.0 Datasets.  _CoRR_ abs/1306.2597 (2013).  http://arxiv.org/abs/1306.2597
  * Qin et al. (2010) Tao Qin, Tie-Yan Liu, and Hang Li. 2010.  A general approximation framework for direct optimization of information retrieval measures.  _Information retrieval_ 13, 4 (2010), 375–397.
  * Rudin (2018) Cynthia Rudin. 2018\.  Stop Explaining Black Box Machine Learning Models for High Stakes Decisions and Use Interpretable Models Instead.  arXiv:stat.ML/1811.10154
  * Sterling and Ravina (2019) Ethan Sterling and Walker Ravina. 2019.  _pwlfit: A Piecewise-Linear Curve Fitting Library_.  https://github.com/google/pwlfit
  * Tan et al. (2018) Sarah Tan, Rich Caruana, Giles Hooker, and Yin Lou. 2018\.  Distill-and-Compare: Auditing Black-Box Models Using Transparent Model Distillation. In _Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society_ _(AIES ’18)_. 303–310.
  * Ustun and Rudin (2014) Berk Ustun and Cynthia Rudin. 2014.  Methods and Models for Interpretable Linear Classification.  arXiv:stat.ME/1405.4047
  * Wikipedia (2020a) Wikipedia. 2020a.  Bhatia–Davis inequality.  http://en.wikipedia.org/w/index.php?title=Bhatia%E2%80%93Davis%20inequality&oldid=875899600.  [Online; accessed 03-September-2020].
  * Wikipedia (2020b) Wikipedia. 2020b.  Isotonic regression.  http://en.wikipedia.org/w/index.php?title=Isotonic%20regression&oldid=989717822.  [Online; accessed 30-November-2020].
  * Wikipedia (2020c) Wikipedia. 2020c.  Kolmogorov–Arnold representation theorem.  http://en.wikipedia.org/w/index.php?title=Kolmogorov%E2%80%93Arnold%20representation%20theorem&oldid=964097101.  [Online; accessed 10-August-2020].
  * Wikipedia (2020d) Wikipedia. 2020d.  Minimum description length.  http://en.wikipedia.org/w/index.php?title=Minimum%20description%20length&oldid=965620302.  [Online; accessed 12-August-2020].
  * Wikipedia (2020e) Wikipedia. 2020e.  Nothing-up-my-sleeve number.  http://en.wikipedia.org/w/index.php?title=Nothing-up-my-sleeve%20number&oldid=972510276.  [Online; accessed 12-August-2020].
  * Wikipedia (2020f) Wikipedia. 2020f.  Segmented regression.  http://en.wikipedia.org/w/index.php?title=Segmented%20regression&oldid=910888930.  [Online; accessed 10-August-2020].
  * Zeiler (2012) Matthew D. Zeiler. 2012\.  ADADELTA: An Adaptive Learning Rate Method.  _CoRR_ abs/1212.5701 (2012).  arXiv:1212.5701 http://arxiv.org/abs/1212.5701
  * Zhuang et al. (2021) Honglei Zhuang, Xuanhui Wang, Michael Bendersky, Alexander Grushetsky, Yonghui Wu, Petr Mitrichev, Ethan Sterling, Nathan Bell, Walker Ravina, and Hai Qian. 2021\.  Interpretable Ranking with Generalized Additive Models. In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_ _(WSDM ’21)_. to appear.

##  Appendix A Experimental Details

Table 2. Metrics across datasets. For AUC-ROC and NDCG@5 higher is better, and for RMSE lower is better. MSLR-WEB30K and CWS were not used by the NAM paper and are omitted from that row. Metric values are the mean from five fold cross validation $\pm$ the sample standard deviation. Model | COMPAS (AUC-ROC) | FICO (RMSE) | MSLR-WEB30K (NDCG@5) | CWS (NDCG@5)
---|---|---|---|---
Interacting forest | $0.742\pm 0.011$ | $3.128\pm 0.092$ | $0.485\pm 0.002$ | $0.461\pm 0.003$
GAM forest | $0.741\pm 0.013$ | $3.495\pm 0.109$ | $0.442\pm 0.002$ | $0.460\pm 0.004$
NAM | $0.741\pm 0.009$ | $3.490\pm 0.081$ |  |
pwlfit num_segments=1, mono=False | $0.741\pm 0.008$ | $3.781\pm 0.105$ | $0.432\pm 0.003$ | $0.449\pm 0.003$
pwlfit num_segments=1, mono=True | $0.741\pm 0.008$ | $3.780\pm 0.105$ | $0.432\pm 0.003$ | $0.450\pm 0.003$
pwlfit num_segments=2, mono=False | $0.742\pm 0.011$ | $3.617\pm 0.101$ | $0.438\pm 0.002$ | $0.452\pm 0.003$
pwlfit num_segments=2, mono=True | $0.742\pm 0.011$ | $3.713\pm 0.103$ | $0.435\pm 0.003$ | $0.451\pm 0.003$
pwlfit num_segments=3, mono=False | $0.743\pm 0.010$ | $3.536\pm 0.099$ | $0.440\pm 0.002$ | $0.453\pm 0.002$
pwlfit num_segments=3, mono=True | $0.743\pm 0.011$ | $3.691\pm 0.100$ | $0.437\pm 0.002$ | $0.451\pm 0.003$
pwlfit num_segments=4, mono=False | $0.742\pm 0.012$ | $3.505\pm 0.094$ | $0.441\pm 0.002$ | $0.454\pm 0.003$
pwlfit num_segments=4, mono=True | $0.742\pm 0.012$ | $3.691\pm 0.101$ | $0.437\pm 0.003$ | $0.452\pm 0.004$
pwlfit num_segments=5, mono=False | $0.743\pm 0.012$ | $3.494\pm 0.096$ | $0.441\pm 0.002$ | $0.454\pm 0.002$
pwlfit num_segments=5, mono=True | $0.743\pm 0.013$ | $3.693\pm 0.101$ | $0.437\pm 0.003$ | $0.452\pm 0.003$
pwlf num_segments=2 | $0.742\pm 0.014$ | $3.728\pm 0.101$ | $0.428\pm 0.004$ | $0.451\pm 0.003$
pwlf num_segments=3 | $0.743\pm 0.012$ | $3.565\pm 0.108$ | $0.434\pm 0.004$ | $0.453\pm 0.003$
pwlf num_segments=4 | $0.744\pm 0.012$ | $3.498\pm 0.099$ | $0.436\pm 0.004$ | $0.453\pm 0.003$
pwlf num_segments=5 | $0.743\pm 0.012$ | $3.503\pm 0.096$ | $0.433\pm 0.003$ | $0.454\pm 0.003$
SGD num_segments=1 | $0.728\pm 0.008$ | $4.349\pm 0.059$ | $0.352\pm 0.002$ | $0.435\pm 0.003$
SGD num_segments=2 | $0.734\pm 0.007$ | $4.028\pm 0.095$ | $0.382\pm 0.001$ | $0.443\pm 0.004$
SGD num_segments=3 | $0.742\pm 0.008$ | $3.887\pm 0.080$ | $0.394\pm 0.002$ | $0.449\pm 0.003$
SGD num_segments=4 | $0.742\pm 0.010$ | $3.742\pm 0.110$ | $0.403\pm 0.003$ | $0.449\pm 0.004$
SGD num_segments=5 | $0.741\pm 0.010$ | $3.643\pm 0.097$ | $0.405\pm 0.002$ | $0.448\pm 0.004$
SGD num_segments=6 | $0.742\pm 0.010$ | $3.583\pm 0.105$ | $0.408\pm 0.002$ | $0.452\pm 0.002$
SGD num_segments=7 | $0.741\pm 0.010$ | $3.604\pm 0.098$ | $0.408\pm 0.002$ | $0.452\pm 0.003$
SGD num_segments=8 | $0.741\pm 0.010$ | $3.561\pm 0.111$ | $0.411\pm 0.003$ | $0.449\pm 0.004$
SGD num_segments=9 | $0.741\pm 0.011$ | $3.522\pm 0.103$ | $0.414\pm 0.002$ | $0.449\pm 0.003$
SGD num_segments=10 | $0.741\pm 0.010$ | $3.544\pm 0.117$ | $0.415\pm 0.002$ | $0.452\pm 0.005$
SGD num_segments=15 | $0.740\pm 0.011$ | $3.499\pm 0.101$ | $0.419\pm 0.003$ | $0.454\pm 0.003$
SGD num_segments=20 | $0.738\pm 0.011$ | $3.499\pm 0.117$ | $0.419\pm 0.003$ | $0.455\pm 0.003$

###  A.1. Cross Validation

We performed 5-fold cross validation on all datasets.

  * •

COMPAS & FICO: The datasets were split into 5 equal parts. Each part was used
once as a test set (20%) with the remaining parts as the training set (80%).
We used the same random folds as in the NAM paper (Agarwal et al., 2020). No
validation set was used given the small size of the data. Instead we used out
of bag evaluation wherever a validation set would be used (see below).

  * •

MSLR-WEB30K: We used the predefined folds and partitions from the original
dataset. For each fold it allocates 60% for training 20% for validation and
20% for testing.

  * •

CWS: We used a dataset of 60,000 queries and 2,690,439 items with an average
of ~44 items per query. The dataset was split into 5 equal parts. Each part
was used once as a test set. Of the remaining parts 80% was used as training
and 20% as validation. Overall this resulted in 64% for training, 16% for
validation and 20% for test for each fold.

###  A.2. Ensemble Learning

For both SGD and tree models, we trained ensembles with 56 bags using a bag
fraction of $\frac{7}{8}$ to produce the random subsets. For MSLR-WEB30K and
CWS, queries were randomly divided into bags. For the other datasets,
individual examples were randomly divided into bags. When applying our
distillation technique, we distilled the ensembles into a single PWLCurve per
feature. When learning the curves directly via SGD, we averaged the learned
$y$-coordinate values across bags to obtain the final model.

###  A.3. Loss Functions

We trained SGD and tree models using log-loss for the COMPAS dataset, mean
squared error (MSE) for the FICO dataset, and ranking loss (ApproxNDCG (Qin et
al., 2010)) for the MSLR-WEB30K and CWS datasets.

###  A.4. Hyper-parameters

For the COMPAS, and FICO datasets hyper-parameters were tuned using out of bag
evaluation on the training set of the first fold. For MSLR-WEB30K and CWS, we
used the validation sets of the first fold.

  * •

SGD: We tuned the batch size in {128, 256, 512, 1024, 4096}. We used the
Adadelta (Zeiler, 2012) optimizer and tuned a sufficient maximum number of
steps for convergence. No other parameters were tuned.

  * •

Interacting forest: We trained depth 5 trees using an internal boosted forest
algorithm. We tuned a sufficient maximum number of steps for convergence. No
other parameters were tuned.

  * •

GAM forest: We trained depth 3 trees restricted to using a single feature with
an internal boosted forest algorithm. We tuned a sufficient maximum number of
steps for convergence. No other parameters were tuned.

In all cases, we trained models for the tuned maximum number of steps and then
truncated models after training. Truncation used a confidence-based truncation
algorithm which attempts to select the earliest step for which no later step
provides a confident win. This algorithm was run on the validation set if
present or otherwise utilized out of bag evaluation.

###  A.5. Code

The GitHub repository for pwlfit (Sterling and Ravina, 2019) contains several
Jupyter notebooks applying our distillation technique and performing the
analyses shown in this paper. Please reference the v0.2.0 release to get the
accompanying data files and appropriate version of the Jupyter notebooks.

##  Appendix B Linear Condense

Linear condensing is a data optimization designed to reduce the runtime
complexity of our piecewise-linear curve fitting.

###  B.1. Motivation/Overview

pwlfit picks a set of candidate $x$-knots and searches through combinations of
those $x$-knots. For each combination considered, it solves a linear least
squares expression for the ideal $y$-knots, calculates the resulting squared
error, and prefers the combination that yields the lowest error.

Each solve is linear in the size of input, which is slow for large data. We
could downsample to save compute at the cost of accuracy. Instead, we
introduce a technique to save compute at no cost in accuracy. We condense the
data into $\mathcal{O}(\\#candidates)$ synthetic points. These synthetic
points perfectly recreate the true squared error over the full data for every
PWLCurve that will be considered. We then optimize over the synthetic points
instead of the real points.

This is possible because we know the candidate $x$-knots ahead of time. A
PWLCurve defined on those $x$-knots will always be linear between any adjacent
$x$-knots in the set of candidates. As we show in the theorem, we can condense
arbitrarily many points down to two points such that linear fits are the same
on those two points as on the full set. In the corollary, we apply this
process separately between each pair of candidate $x$-knots, producing two
points between each pair. Together, the squared error of such a PWLCurve is
the same on those synthetic points as it is on the full data set. (Up to a
constant that we safely ignore because it’s the same for each PWLCurve.)

###  B.2. Definitions

For convenience, we take standard definitions and specialize them for weighted
2D points.

###### Definition B.1.

Let a ‘point’ refer to a real-valued triple of the form $(x,y,weight)$ with
$weight>0$.

###### Definition B.2.

Let a ‘line’ refer to a function of the form $f(x)=mx+b$ for
$m,b,x\in\mathbb{R}$.

###### Definition B.3.

For any function $f:\mathbb{R}\to\mathbb{R}$ and finite point set $P$, define
the squared error $SE(f,P)$ as the sum of $(f(x)-y)^{2}\cdot weight$ for each
point in $P$. If $P$ is empty, we consider the squared error to be $0$.

###### Definition B.4.

For any finite point set $P$, define the ‘best fit line’ $bestfitline(P)$ as
the line $L$ that minimizes $SE(L,P)$. In the degenerate case where multiple
lines minimize $SE$, let the best fit line be the solution with zero slope,
and if multiple solutions have zero slope, let the best fit line be the
solution with a zero $y$-intercept.

There are two degenerate cases that require tie-breaking. If the point set is
empty, every line has the same squared error, so our definition chooses
$f(x)=0$ as the best fit line. If the point set is nonempty but all its points
have the same $x$, then any line with the correct value at $x$ will minimize
the squared error, so our definition choose the horizontal line.

###  B.3. Theorem

###### Theorem B.5.

Given a set of points $P$, we can construct a set $P^{\prime}$ of two or fewer
points such that

1\. $min_{x}(P)<=min_{x}(P^{\prime})<=max_{x}(P^{\prime})<=max_{x}(P)$, and

2\. For any line $L$, $SE(L,P)=SE(L,P^{\prime})+SE(bestfitline(P),P)$.

######  Remark.

These properties are desirable because (2) allows us to compute the squared
error of $M$ lines over a data set of $N$ points in $\mathcal{O}(N+M)$ instead
of the naive $\mathcal{O}(NM)$, and (1) allows us to extend this property from
lines to a useful class of piecewise-linear curves in the corollary.

Note that the points in $P^{\prime}$ are constructed, rather than chosen from
$P$. The construction of $P^{\prime}$ is implemented in pwlfit (Sterling and
Ravina, 2019) as linear_condense.linear_condense.

###### Proof.

Let $X$, $Y$, and $W$ represent the $x,y$, and $weight$ values of $P$,
respectively. We dismiss the trivial case where $P$ is empty; in that case, an
empty $P^{\prime}$ satisfies the requirements. Likewise, we dismiss the case
where $min(X)=max(X)$ since $P^{\prime}=\\{Centroid(P)\\}$ fulfills our
desired properties. With those cases resolved, we assume for the rest of this
proof that $min(X)<max(X)$.

####  B.3.1. Reframe the Coordinate System

To begin, we reframe the coordinate system such that the origin is the
centroid of $P$ and $y=0$ is the best fit line. (This simplifies the math.) We
ensure that the shift of coordinates is reversible and preserves the squared
error.

$Centroid(P)=(X\cdot W/sum(W),Y\cdot W/sum(W))$. We translate the coordinate
frame by this centroid so that, under the new coordinates,
$Centroid(P)=(0,0)$. After translation, $X\cdot W=0$ and $Y\cdot W=0$.

Additionally, we skew the coordinate system by the slope of the best fit line:
we replace $Y$ with $Y-X\cdot slope(bestfitline(P))$. With the centroid at the
origin, the slope of the best fit line is $Covariance(X,Y,W)/Variance(X,W)$ =
$(XY\cdot W)/(XX\cdot W)$. After skewing this slope to 0, $XY\cdot W$ = 0.

Under the new coordinate frame, $SE(bestfitline(P),P)=SE(y=0,P)=Y^{2}\cdot W$.

We will determine $P^{\prime}$ under this new coordinate system. Afterwards,
we can easily convert $P^{\prime}$ back to the original coordinate system by
reversing the skew and the translation.

####  B.3.2. Squared Error of an arbitrary line

We will express $SE(line,P)$ as $SE(bestfitline(P),P)$ plus leftover terms.
From that, we will derive a $P^{\prime}$ such that $SE(line,P^{\prime})$
equals those leftover terms.

For an arbitrary line $y=mx+b$,

$SE(y=mx+b,P)=(mX+b-Y)^{2}\cdot W=(m^{2}X^{2}+2bmX-2mXY+b^{2}-2bY+Y^{2})\cdot
W.$

In our coordinate frame, $X\cdot W=0$, $Y\cdot W=0$, and $XY\cdot W=0$. So
$SE(y=mx+b,P)=(m^{2}X^{2}+b^{2}+Y^{2})\cdot W.$

$Y^{2}\cdot W=SE(bestfitline(P),P)$. Therefore,

$SE(y=mx+b,P)=m^{2}X^{2}\cdot W+b^{2}\cdot W+SE(bestfitline(P),P).$
$m^{2}X^{2}\cdot W+b^{2}\cdot W=SE(y=mx+b,P)-SE(bestfitline(P),P).$

####  B.3.3. Squared error over $P^{\prime}$

$SE(y=mx+b,P^{\prime})=SE(y=mx+b,P)-SE(bestfitline(P),P)$

for all lines $y=mx+b$ $\iff$ $(mX^{\prime}+b-Y^{\prime})^{2}\cdot
W^{\prime}=m^{2}X^{2}\cdot W+b^{2}\cdot W$ for all lines $y=mx+b$.

The above equation can be viewed as a quadratic polynomial in the two
variables $m$ and $b$. To hold for all values of $m$ and $b$, the coefficients
of each $m^{c}b^{d}$ must be equal on both sides of the equation. Then the
equation holds iff:

1\. $X^{\prime 2}\cdot W^{\prime}=X^{2}\cdot W$, and

2\. $X^{\prime}\cdot W^{\prime}=0$, and

3\. $sum(W)=sum(W^{\prime})$, and

4\. $Y^{\prime}\cdot W^{\prime}=0$, and

5\. $Y^{\prime 2}\cdot W^{\prime}=0$, and

6\. $X^{\prime}Y^{\prime}\cdot W^{\prime}=0$.

(5) $\iff$ $Y^{\prime}=0$, which also guarantees (4) and (6). We will use 1-3
to derive a satisfactory $X^{\prime}$ and $W^{\prime}$.

####  B.3.4. Deriving $X^{\prime}$ and $W^{\prime}$

We’ve determined that $Y^{\prime}=0$.

Let $X^{\prime}:=(x_{1},x_{2})$ and $W^{\prime}:=(w_{1},w_{2})$. Without loss
of generality, let $x_{1}$ ¡= $x_{2}$. Then, to satisfy our squared error
expression, it’s necessary and sufficient that:

1\. $x_{1}^{2}w_{1}+x_{2}^{2}w_{2}=X^{2}\cdot W$, and

2\. $x_{1}w_{1}+x_{2}w_{2}=0$, and

3\. $w_{1}+w_{2}=sum(W)$.

Because we have three equations in four unknowns, we cannot directly solve for
$x_{1},x_{2},w_{1},w_{2}.$ To produce a fourth equation, we choose the
constraint that $x_{1}/x_{2}$ = $min(X)/max(X)$. This choice will simplify the
math, and will ensure that $min(X)<=x_{1}<=x_{2}<=max(X)$.

With this fourth equation, we solve the simultaneous equations to produce:

$x_{1}=-stddev(X,W)\sqrt{-min(X)/max(X)}$

$x_{2}=stddev(X,W)\sqrt{max(X)/-min(X)}$.

$w_{1}=sum(W)\cdot max(X)/(max(X)-min(X))$

$w_{2}=sum(W)\cdot-min(X)/(max(X)-min(X))$.

Note that, because the centroid is zero, $min(X)<0<max(X)$, so these
expressions are all defined. (The denominators are never 0 and values beneath
the square roots are never negative.)

$P^{\prime}={(x_{1},0,w_{1}),(x_{2},0,w_{2})}$ satisfies our requirements.

####  B.3.5. Verify that $min(X)<=x_{1}<=x_{2}<=max(X)$

We wanted $P^{\prime}$ to satisfy the the squared error expression, which it
does, and also have its x-values bounded by the x-values of $P$, which we
prove now. Let $\mu:=E(X,W)$, the expected value of $X$ weighted by $W$, which
is equivalent to the x-value of $Centroid(P)$. By the Bhatia–Davis inequality
(Wikipedia, 2020a),

$stddev(X,W)^{2}<=(\mu-min(X))(max(X)-\mu)$. (This inequality is equivalent to
the observation that the standard deviation of a distribution is maximized
when all the xs are at the extremes – i.e. equal to min(X) or max(X).)

Since $\mu$ is zero for $P$, $stddev(X,W)^{2}<=-min(X)max(X)$.

$x_{1}^{2}=stddev(X,W)^{2}\cdot(-min(X)/max(X))<=-min(X)max(X)\cdot(-min(x)/max(X))=min(X)^{2}.$

$x_{1}<0$ and $min(X)<0$, so $x_{1}^{2}<=min(X)^{2}\implies min(X)<=x_{1}$.
The proof that $x_{2}<=max(X)$ is similar.

Therefore $min(X)<=x_{1}<=x_{2}<=max(X)$, as desired. ∎

###  B.4. Corollary

###### Corollary B.6.

Given a set of points $P$ and a set of x-knots $K$, we can construct a set of
points $P^{\prime}$ with $|P^{\prime}|<=2(|K|-1)$ such that, for any PWLCurve
$C$ whose x-knots are elements of $K$, $SE(C,P)=SE(C,P^{\prime})+c$, where $c$
is a constant determined exclusively by $P$ and $K$ that’s the same for every
$C$.

Note that the points in $P^{\prime}$ are constructed, rather than chosen from
$P$. The construction of $P^{\prime}$ is implemented in pwlfit (Sterling and
Ravina, 2019) as linear_condense.condense_around_knots.

####  B.4.1. Preprocess $P$ by clamping

Let $k:=|K|$, and consider $K$ in sorted order. Piecewise-linear curves are
constant for input values that exceed the range of their x-knots, so $C$ is
constant for $x<=min(K)=K[0]$ and for $x>=max(K)=K[k-1]$.

Therefore we can clamp the x-values of $P$ to $[K[0],K[k-1]]$ without altering
$SE(C,P)$. We do so as a preprocess.

####  B.4.2. Partition $P$ by $K$

To construct $P^{\prime}$ from $P$, we first partition $P$ by $K$ into $k+1$
disjoint pieces labeled $P_{0}$, $P_{1}$, …, $P_{k}$.

\- $P_{0}$ := $\\{x\in P|x<K[0]\\}$.

\- $P_{i}$ := $\\{x\in P|K[i-1]<=x<K[i]\\}$ for $1<=i<=k-2$.

\- $P_{k-1}$ := $\\{x\in P|K[k-2]<=x<=K[k-1]\\}$.

\- $P_{k}$ := $\\{x\in P|K[k-1]<x\\}$.

Because we clamped $P$, $P_{0}$ and $P_{k}$ are empty. Therefore
$\bigcup_{i=1}^{k-1}P_{i}=P$.

A PWLCurve is linear between consecutive control points, so $C$ is linear over
each $P_{i}$.

####  B.4.3. Convert each partition into two points

From the theorem, for each $P_{i}$, we can produce a two-point set
$P_{i}^{\prime}$ with
$min_{x}(P_{i})<=min_{x}(P_{i}^{\prime})<=max_{x}(P_{i}^{\prime})<=max_{x}(P_{i})$,
such that for any line $L$,

$SE(L,P_{i})=SE(L,P_{i}^{\prime})+SE(bestfitline(P_{i}),P_{i})$. $C$ is linear
over each $P_{i}$, so

$SE(C,P_{i})=SE(C,P_{i}^{\prime})+SE(bestfitline(P_{i}),P_{i})$.

####  B.4.4. Recombine partitions

Let $P^{\prime}:=\bigcup_{i=1}^{k-1}P_{i}^{\prime}$. Each $P_{i}^{\prime}$
consists of two points, so $P^{\prime}$ consists of $2(|K|-1)$ points.

$\displaystyle SE(C,P)$ $\displaystyle=\sum_{i=1}^{k-1}SE(C,P_{i})$
$\displaystyle={\sum_{i=1}^{k-1}(SE(C,P_{i}^{\prime})+SE(bestfitline(P_{i}),P_{i}))}$
$\displaystyle={SE(C,P^{\prime})+\sum_{i=1}^{k-1}SE(bestfitline(P_{i}),P_{i})}.$

$\sum_{i=1}^{k-1}SE(bestfitline(P_{i}),P_{i})$ is determined by $P$ and $K$,
and is therefore the same for every $C$. Therefore we’ve proven the corollary.


$\displaystyle=\frac{({2\pi}t)^{2}}{2}\text{Tr}\Big{[}\lambda^{z}(f_{T}-f_{S})\left(\hat{\rho}_{T}\hat{\rho}_{S}+t^{2}\hat{\rho}_{T}(G_{SS}g_{T}\hat{\rho}_{S}+\hat{\rho}_{S}g_{T}^{\dagger}G_{SS}^{\dagger})+t^{4}\hat{\rho}_{T}G_{SS}g_{T}\hat{\rho}_{S}g_{T}^{\dagger}G_{SS}^{\dagger}\right)+t^{2}\lambda^{z}\hat{\rho}_{T}[f_{T},G_{SS}]\hat{\rho}_{T}G^{\dagger}_{SS}\Big{]},$
(67)

where again the trace includes an integral over frequency. It now remains only
to calculate $G_{SS}$, which can be done easily with the help of the last
Dyson equation in (60):

$G_{SS}=(g^{-1}_{SS}-t^{2}g_{T})^{-1}.$ (68)

The first group of terms in (67) vanishes when $\hat{\rho}_{S}$ does; these
terms are therefore responsible for normal tunneling processes that are only
active outside of the superconducting gap. The last term on the other hand
does not depend on the superconducting DOS, and can lead to a current even for
biases inside the gap; this is therefore the term responsible for the in-gap
Andreev current.

The conductance is obtained by differentiating (67) with respect to $V$,
yielding

$\displaystyle\frac{dI}{dV}$
$\displaystyle=\frac{({2\pi}t)^{2}}{2}\int_{\mathbb{R}}d\omega\,f^{\prime}(\omega-V)\text{Tr}\left[\hat{\rho}_{T}\hat{\rho}_{S}+t^{2}\hat{\rho}_{T}\left(G_{SS}g_{T}\hat{\rho}_{S}+\hat{\rho}_{S}g_{T}^{\dagger}G_{SS}^{\dagger}\right)+t^{4}\hat{\rho}_{T}G_{SS}g_{T}\hat{\rho}_{S}g_{T}^{\dagger}G_{SS}^{\dagger}+2t^{2}\hat{\rho}_{T}G_{SS,a}\hat{\rho}_{T}G^{\dagger}_{SS,a}\right](\omega),$
(69)

where $G_{SS,a}$ is the anomalous part of $G_{SS}$ (viz. the part of $G_{SS}$
off-diagonal in Nambu indices), and where we have now written out the
frequency integral explicitly.

####  D.2.1 Momentum non-conserving tunneling

We first consider the limit of ‘incoherent’ tunneling, where the tunneling
process does not conserve in-plane momentum. Therefore the tip Green’s
function is, treating the tip as being in the wide band limit with a constant
density of states $\rho_{T}$ and remembering the projection onto the SVL
subspace,

$g_{T}(\omega)=\int_{\mathbb{R}}d\xi\,\frac{\rho_{T}}{\omega+i\eta-\lambda^{z}\xi}\mathscr{P}=-\pi
i\hat{\rho}_{T}\lambda^{0}\mathscr{P}$ (70)

where

$\mathscr{P}\equiv\mathcal{P}\oplus\mathcal{P}^{T},$ (71)

with the $\oplus$ in Nambu space.

On the other hand, the free Green’s function of the SC in the incoherent limit
(a matrix in spin and Nambu space) is, reading off from (25) (and continuing
to let $\ket{\eta},\ket{\eta^{\prime}}$ be real),

$\displaystyle g_{S}(\omega)$
$\displaystyle=\sum_{\tau,\tau^{\prime}}\int_{\mathbf{k}}\left[\frac{\mathscr{P}(\omega+\lambda^{z}(\xi_{\alpha^{z}\mathbf{k}}-\tau^{z}\delta_{\alpha^{z}\mathbf{k}}))-\widehat{\Delta}_{\theta}\lambda^{+}-\widehat{\Delta}_{\theta}^{\dagger}\lambda^{-}}{(\omega-\alpha^{z}\delta_{\alpha^{z}\mathbf{k}})^{2}-E^{2}_{\alpha^{z}\mathbf{k}}}\right]_{\tau\tau^{\prime}}$
(72) $\displaystyle=-\frac{1}{{4\pi}}\int_{\theta}\Big{(}\sum_{s=\pm
1}\frac{1+s\lambda^{z}}{2}\left(\Upsilon_{\theta,s}(\omega-s\delta_{\theta_{s}})\ket{\eta}\langle\eta|+\Upsilon_{\theta,-s}(\omega+s\delta_{\theta_{-s}})\ket{\eta^{\prime}}\langle\eta^{\prime}|\right)$
$\displaystyle\qquad-\Delta_{\theta}\lambda^{+}(\Upsilon_{\theta,-}\ket{\eta^{\prime}}\langle\eta|-p\Upsilon_{\theta,+}\ket{\eta}\langle\eta^{\prime}|)-\Delta_{\theta}^{*}\lambda^{-}(-p\Upsilon_{\theta,+}\ket{\eta^{\prime}}\langle\eta|+\Upsilon_{\theta,-}\ket{\eta}\langle\eta^{\prime}|)\Big{)},$

where $\alpha^{z}\equiv\lambda^{z}\tau^{z}$ and we have defined

$\Upsilon_{\theta,\pm}\equiv\frac{m_{\theta_{\pm}}}{\sqrt{|\Delta_{\theta}|^{2}-(\omega\mp\delta_{\theta})^{2}}},$
(73)

with the notation $\theta_{+}=\theta,\theta_{-}=\theta+\pi$.

We will now use (72) to argue that either a) TBG is actually an $s$-wave SC
(which for the reasons explained in the main text we are inclined to disfavor)
or b) the tunneling experiments of [22] are in fact rather well in the
momentum-conserving tunneling regime (which given the experimental setup is
perhaps slightly surprising).

To argue this, all we need to do is note the rather strong zero-bias peak in
$dI/dV$ observed at $\mathbf{B}=0$ in [22]; we claim that (72) cannot
generically reproduce such a peak. Indeed, setting $\delta_{\theta}=0$ (and
specializing to the case of AF SVL), we have

$\displaystyle g_{S}(\omega)$
$\displaystyle=-\int\frac{d\theta}{{2\pi}}\,\Upsilon_{\theta,+}(\omega+(\lambda^{+}\Delta_{\theta}-\lambda^{-}\Delta^{*}_{\theta})is^{y}).$
(74)

Since $m_{\theta}$ is the mass in the $K$ valley, $C_{3}$ symmetry imposes
$m_{\theta}=m_{\theta+2\pi/3}$, and so as in (49) we parametrize

$m_{\theta}=m(1+\zeta\cos(3\theta)),$ (75)

with $|\zeta|<1$ controlling the extent of the angular anisotropy. Now at
$T=0$ the zero bias conductance only depends on $g_{S}(0)$, which is (taking
$\Delta_{\theta}\in\mathbb{R}$ wolog)

$\displaystyle g_{S}(0)$
$\displaystyle=\int\frac{d\theta}{{2\pi}}\,m_{\theta}\,{\rm
sgn}(\Delta_{\theta})\lambda^{y}s^{y}$ (76)
$\displaystyle=m\int\frac{d\theta}{{2\pi}}\,\left(\frac{1+p}{2}+\zeta\frac{1-p}{2}\cos(3\theta)\right){\rm
sgn}(\Delta_{\theta}).$

This means that the strength of the Andreev reflection signal is proportional
to $\zeta$ if $p=-1$, and is controlled by $m\int d\theta\,{\rm
sgn}(\Delta_{\theta})$ if $p=+1$. On the other hand, the conductance in the
$V\rightarrow\infty$ limit is determined only by $\int
d\theta\,m_{\theta}/{2\pi}=m$. This leads to the normalized ZBC being
depressed from the maximal value of 2 it takes in the s-wave case to something
smaller, with the suppression being either by an integral like $\zeta\int
d\theta\cos(3\theta){\rm sgn}(\Delta_{\theta})$ in the $p=-1$ case, or $\int
d\theta\,{\rm sgn}(\Delta_{\theta})$ in the $p=1$ case.

(a)

(b)

Figure 6:  $dI/dV$ curves for $\Delta_{\theta}=\Delta\cos(\theta)$, normalized
to their respective $V\rightarrow\infty$ normal-state values. Here
$\mathsf{t}$ is the dimensionless tunneling strength defined by
$\mathsf{t}\equiv t\sqrt{m\pi\rho_{T}}$. In terms of BTK theory,
$\mathsf{t}=1$ corresponds to the case of a transparent barrier, $Z=0$ (and
$\mathsf{t}>1$ corresponds to an attractive barrier, with $Z<0$).

The net result of these effects is to produce tunneling curves that are in
general significantly suppressed near zero bias. For a p-wave order parameter
$\Delta_{\theta}=\Delta\cos(\theta)$, the normalized conductance is shown in
figure 6, for an isotropic effective mass ($\zeta=0$, left) and a large
$C_{3}$ anisotropic mass ($\zeta=0.7$, right). The conductance in the two
cases looks similar, but the anisotropy $\zeta\neq 0$ is responsible for a
(very small) nonzero ZBC in the latter case.

(a)

(b)

Figure 7:  Same as in figure 6, but with
t$\Delta_{\theta}=\Delta(\cos(2\theta)+\varsigma\cos(4\theta))$. The multiple
singularities are due to the fact that the gap possesses more than one
harmonic.

A nonzero ZBC is obtained in the $p=1$ case only by working with an order
parameter for which $\int d\theta\,{\rm sgn}(\Delta_{\theta})\neq 0$. A simple
example of such an order parameter is
$\Delta_{\theta}=\Delta(\cos(2\theta)+\varsigma\cos(4\theta))$, which has 4
nodes if $|\varsigma|<1$. Examples of the normalized conductance for
$\varsigma=0.25,0.75$ are shown in figure 7. The nonzero value of $\varsigma$
allows the ZBC to be finite (though still suppressed relative to its normal
state value), but the fact that $\Delta_{\theta}$ includes multiple harmonics
leads to multiple coherence peaks at larger biases, which are not observed in
experiment. Even with $\varsigma=0.75$, we are clearly a long way away from
producing a peak in the ZBC.

####  D.2.2 Momentum conserving tunneling

Given the issues with incoherent tunneling we now consider the case of
coherent tunneling, where the tunneling between the tip and the sample
conserves in-plane momentum in the TBG BZ.101010One could also consider an
intermediate case, where the tunneling is incoherent in the TBG BZ but
diagonal in valley; this leads to the same problems as in the valley-
incoherent case as they originate from the integral over $\theta$, rather than
the sum over valleys. To treat this case carefully, we would work with
momentum-resolved Green’s functions, and include a single integral over
momentum in the trace appearing in (69). We will instead do something slightly
simpler: for a given point on the Fermi surface at angle $\theta$, we will
take the tunneling to conserve momenta normal to $\widehat{\bm{\theta}}$
(along the Fermi surface), but not parallel to $\widehat{\bm{\theta}}$ (normal
to the Fermi surface). This treatment nevertheless captures the essential way
in which coherence affects the conductance, since in this treatment there is
no interference between the signs of $\Delta_{\theta}$ at different angles
$\theta$.

(a)

(b)

Figure 8:  Tunneling conductance in the coherent tunneling model for a
$p$-wave gap ($d$-wave is essentially the same), normalized to its normal-
state value. The mass anisotropy parameter is fixed as $\zeta=1/4$, and
$\mathsf{t}=t\sqrt{m\pi\rho_{T}}$ is the dimensionless tunneling strength as
in figure 6. Left: zero magnetic field. Right: Field of magnitude
$B_{\,\parallel\,}=0.2T$ and $\theta_{B}=\pi/4$, $\gamma=0$, and with the
parameters in (50) fixed at $\delta_{0}=4\delta_{3}=0.15$meV/T.

In this framework, the appropriate Green’s function (a matrix in valley, spin,
and Nambu space) to use when computing $dI/dV$ is

$\displaystyle g_{S}(\omega,\theta)$
$\displaystyle=-\frac{1}{{4\pi}}\Big{(}\Upsilon_{\theta,\alpha^{z}}(\omega-\alpha^{z}\delta_{\theta_{\alpha^{z}}})\mathscr{P}-\lambda^{+}\Delta_{\theta}\left(\tau^{-}\Upsilon_{\theta,-}\ket{\eta^{\prime}}\langle\eta|-p\tau^{+}\Upsilon_{\theta,{+}}\ket{\eta}\langle{\eta^{\prime}}|\right)-\lambda^{-}\Delta_{\theta}^{*}\left(-p\tau^{-}\Upsilon_{\theta,+}|\eta^{\prime}\rangle\langle\eta|+\tau^{+}\Upsilon_{\theta,-}|\eta\rangle\langle\eta^{\prime}|\right)\Big{)},$
(77)

where $\theta$ is integrated over as part of the trace in (69). The
conductance computed in this framework gives results much more in accordance
with experiment, an example of which is shown in figure 8. As an in-plane
magnetic field modifies the conductance essentially by way of it acting as a
chemical potential for quasiparticles, it has the effect of flattening out the
conductance curves near zero bias, as well as splitting the coherence peak at
low tunneling strength. In the right panel of figure 8 we have taken
$B_{\,\parallel\,}=0.2$T, which in the experiment of [14] is rather
comfortably below the minimum $B_{c,{\,\parallel\,}}$ across the majority of
the superconducting dome. The $\theta_{B}$-dependence of $dI/dV$ at low bias
is essentially the same as that of the $\theta_{B}$-dependence of
$\rho_{S}(0)$ studied above.