# Conformal prediction for frequency-severity modeling

Helton Graziadei 111Corresponding author ✉<EMAIL_ADDRESS>
School of Applied Mathematics, Getulio Vargas Foundation, Rio de Janeiro,
Brazil

Paulo C. Marques F.
Insper Institute of Education and Research, São Paulo, Brazil

Eduardo F. L. de Melo
School of Applied Mathematics, Getulio Vargas Foundation, Rio de Janeiro,
Brazil
SUSEP - Superintendence of Private Insurance, Rio de Janeiro, Brazil
UERJ - State University of Rio de Janeiro, Brazil

Rodrigo S. Targino
School of Applied Mathematics, Getulio Vargas Foundation, Rio de Janeiro,
Brazil

(July 2023)

###### Abstract

We present a nonparametric model-agnostic framework for building prediction
intervals of insurance claims, with finite sample statistical guarantees,
extending the technique of split conformal prediction to the domain of two-
stage frequency-severity modeling. The effectiveness of the framework is
showcased with simulated and real datasets. When the underlying severity model
is a random forest, we extend the two-stage split conformal prediction
procedure, showing how the out-of-bag mechanism can be leveraged to eliminate
the need for a calibration set and to enable the production of prediction
intervals with adaptive width.

Keywords: Frequency-severity modeling; Prediction sets; Two-stage split
conformal prediction; Random forests; Two-stage out-of-bag conformal
prediction.

##  1 Introduction

The statistical modeling of insurance claims is a crucial task of the property
and casualty insurance industry. An essential ingredient in this modeling
process is the two-stage approach, encompassing a frequency model and a
severity model. At the first stage, a frequency model predicts the number of
claims, while, at the second stage, a severity model predicts the average
financial impact or size of a claim, given that it has occurred. Together,
these two models map relevant predictors such as the policyholder’s age,
geographical location, and claim history, to the response variables describing
the frequency and severity of the claims. This classic approach, known as the
frequency-severity model, has been instrumental in the process of risk
categorization, premium calculation, and, in a broader context, risk
quantification of business portfolios for specific industry segments [1, 2].

More specifically, for insurance pricing, the form of the two-stage model
varies according to the characteristics of the available data[1]. A common
approach is to incorporate the modeling of two assumed independent variables:
the frequency of claims, representing the number of claims within a given
portfolio, and the severity (cost) associated with a claim. The pure premium
is then calculated multiplying the predicted frequency and severity of a
claim. Another approach is to model the so called aggregated losses, defined
as the sum of claims in subsets of the underlying portfolio. Our results are
applicable to both contexts.

Traditionally, generalized linear models [3] have been the standard choice for
both frequency and severity stages. A count regression model, typically based
on the Poisson or the negative binomial distributions, is employed for the
frequency variable, while a continuous regression model is designed for the
severity component, commonly using the gamma or the log-normal distributions.
In recent years, we have witnessed the inception of a paradigm shift, with the
insurance industry increasingly gravitating towards the use of modern machine
learning techniques [4, 5] for the two modeling stages.

However, the plethora of currently available machine learning algorithms share
a common limitation: their focus is on pointwise predictions, without an
appropriate accounting of the confidence in their forecasts. In light of this,
the aim of this paper is to explore a model agnostic, universally applicable
framework, suitable for frequency-severity modeling tasks, which is capable of
producing prediction intervals with statistical guarantees.

In Section 2 we discuss the general process of frequency-severity modeling.
Section 3 presents the key ideas behind split conformal prediction [6, 7, 8],
showing how this technique can be extended to the two-stage frequency-severity
scenario. The framework is exemplified in Section 4 using simulated and real
datasets, including a novel dataset of crop insurance claims in Brazil. In
Section 5 we show how to extend the two-stage split conformal prediction
procedure when the underlying predictive models are random forests, exploiting
the out-of-bag mechanism as a way to eliminate the need for a calibration set,
as well as a means to produce prediction intervals of adaptive width. We state
our conclusions in Section 6.

### Related work

To the best of the authors knowledge, this is the first application of
conformal prediction techniques to the process of frequency-severity modeling.

##  2 Frequency-severity modeling

Let $(X_{1},D_{1},Y_{1}),\dots,(X_{n},D_{n},Y_{n})$ be a sample of independent
and identically distributed triplets, in which, for the $i$-th sample unit,
$X_{i}\in\mathbb{R}^{p}$ is a vector of predictors, $D_{i}\in\mathbb{N}$
denotes the corresponding number of incurred claims, and $Y_{i}\in\mathbb{R}$
is the severity of the corresponding claims, following the definition given in
[4, 5]. More specifically, in Section 4.2, we analyse a motor third party
liability dataset from the Belgium insurance market, in which each row
corresponds to a policyholder. In this case, $D_{i}$ is the number of incurred
claims, while $Y_{i}$ stands for the average claim cost of a claim. The
structure of the data is illustrated in Table 1.Furthermore, our results are
also applicable to cases in which there are some level of data aggregation.
For instance, in Section 4.3, a soybean crop insurance dataset from Brazil is
analyzed, in which each row is an aggregation of all contracts at the
municipality level, with $D_{i}$ representing the number of claims and $Y_{i}$
being the average claim cost of within a municipality.

Policy | Claim Count ($D_{i}$) | Total Claim Amount | Severity ($Y_{i}$) | Age ($X_{i,1}$) | Sex ($X_{i,2}$)
---|---|---|---|---|---
1 | $2$ | $2,\\!852.76$ | $1,\\!426.38$ | 24 | male
2 | 0 | $0$ | $0$ | 64 | female
3 | $1$ | $1,\\!618.00$ | $1,\\!618.00$ |  $5$0 | male
Table 1: Structure of the motor third party liability dataset discussed in
Section 4.2, showing the first three policies and the values of two
predictors: $X_{i,1}$ is the age of the policyholder (in years), and $X_{i,2}$
is the sex of the policyholder.

We assume a two-stage data generating process of the form

$D_{i}\mid X_{i}\sim F,\qquad\qquad Y_{i}\mid X_{i},D_{i}\sim G,$

in which the severity distribution $G$ is such that $Y_{i}=0$, whenever
$D_{i}=0$. The total claim cost is predicted using two regression models, a
first one for the frequency component and a second one for the severity
component. Among the many available options, Poisson regression and gamma
regression, respectively, are commonly employed to model the two referred
components. Modern machine learning models, such as random forests [9], have
been recently applied in this two-stage scenario, with remarkable improvements
in terms of predictive performance [4, 5].

##  3 Two-stage split conformal prediction

Conformal prediction [6, 10, 7, 11, 8, 12] was developed in the late 1990s as
a nonparametric way to quantify the confidence in the forecasts made by
general predictive models through the construction of prediction intervals
with finite sample coverage guarantees 222The coverage guarantees refers to
the probability of the prediction interval will contain the observed value of
the response variable within the test set. It does not relate to the type of
coverage provided by an insurance contract..

Figure 1: An illustration of the split conformal procedure for simple linear
regression. The six gray and five black points represent the training and
calibration samples, respectively. The gray line is the least squares line
fitted to the training sample. The size of the black segments are the values
of the calibration conformity scores. The nominal miscoverage level
$\alpha=40\%$, so that $\lceil(1-\alpha)((\text{calibration sample
size})+1)\rceil=4$. A future predictor $x_{n+1}$ and the corresponding
prediction interval is depicted in the figure.

In a regression context, we have a sample of independent and identically
distributed333This assumption can be weakened by assuming only exchangeability
of the pairs. pairs $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$, in which
$X_{i}\in\mathbb{R}^{d}$ is a vector of predictors and $Y_{i}\in\mathbb{R}$ is
a response variable. In the split conformal prediction scheme [10, 8], we
randomly split the sample units indexes into two disjoint sets $I_{1}$ and
$I_{2}$, with $I_{1}\cup I_{2}=\\{1,\dots,n\\}$, corresponding to training and
calibration samples, with sizes $n_{1}$ and $n_{2}$, respectively. Some
regression model $\hat{\mu}$ is trained using the pairs in the training
sample, and this model is used to compute the conformity scores
$R_{i}=|Y_{i}-\hat{\mu}(X_{i})|$ in the calibration set, i.e., for $i\in
I_{2}$. Let $R_{(1)}<R_{(2)}<\dots<R_{(n_{2})}$ denote the ordered calibration
conformity scores.

The key idea of split conformal prediction is that, due to the assumed
distributional symmetry of the data, the conformity score
$R_{n+1}=|Y_{n+1}-\hat{\mu}(X_{n+1})|$ for a future random pair
$(X_{n+1},Y_{n+1})$ is ranked uniformly among the ordered calibration
conformity scores, implying that (assuming that there are no ties among the
conformity scores)

$1-\alpha\leq P(Y_{n+1}\in C^{(1-\alpha)}(X_{n+1}))\leq
1-\alpha+\frac{1}{n_{2}+1},$

for every training and calibration sample sizes, and any prescribed nominal
miscoverage level $0<\alpha<1$, in which the conformal prediction set

$C^{(1-\alpha)}(X_{n+1})=[\,\hat{\mu}(X_{n+1})-\hat{r},\hat{\mu}(X_{n+1})+\hat{r}\,],$

with $\hat{r}=R_{(\lceil(1-\alpha)(n_{2}+1)\rceil)}$, where $\lceil t\rceil$
denotes the smallest integer greater than or equal to the real number $t$ (see
the proof of Theorem 2 in the supplementary materials of [8] for the formal
details).

It is worth emphasizing that property $(\star)$ holds for high dimensional
predictors ($d\gg 1$), no matter what model or algorithm is used to construct
$\hat{\mu}$. All these combined features account for the so-called
universality of the split conformal prediction framework. Figure 1 provides
some intuition on the split conformal prediction procedure in a simple case.

Algorithm 1 Two-stage split conformal prediction

1:Dataset $\\{(x_{i},d_{i},y_{i})\\}_{i=1}^{n}$, training and calibration
indexes $I_{1}$ and $I_{2}$, future vector of predictors
$x_{n+1}\in\mathbb{R}^{p}$, and nominal miscoverage level $0<\alpha<1$.

2:

3:function PredictionInterval($\\{(x_{i},d_{i},y_{i})\\}_{i=1}^{n}$, $I_{1}$,
$I_{2}$, $x_{n+1}$, $\alpha$)

4: Train frequency model $\hat{\mu}$ from $\\{(x_{i},d_{i})\\}_{i\in I_{1}}$

5: if $\hat{\mu}(x_{n+1})=0$ then

6: return $\emptyset$

7: end if

8: $I^{+}_{2}\leftarrow\\{i\in I_{2}:\hat{\mu}(x_{i})>0\\}$

9: $n^{+}_{2}\leftarrow|I^{+}_{2}|$

10: Train severity model $\hat{\psi}$ from $\\{(x_{i},d_{i},y_{i}):i\in
I_{1},d_{i}>0\\}$

11: for $i\in I^{+}_{2}$ do

12: $r_{i}\leftarrow|y_{i}-\hat{\psi}(x_{i},\hat{\mu}(x_{i}))|$

13: end for

14: $\hat{r}\leftarrow r_{(\lceil(1-\alpha)(n^{+}_{2}+1)\rceil)}$

15: return
$[\,\max\,\\{\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))-\hat{r},0\\},\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))+\hat{r}\,]$

16:end function

Extension of split conformal prediction to the two-stage frequency-severity
scenario consists in maintaining the distributional symmetries of the original
procedure while accounting for the presence of the different models on each
one of the two stages.

We start with the random sample
$(X_{1},D_{1},Y_{1}),\dots,(X_{n},D_{n},Y_{n})$ described in Section 2. The
data is split into training and calibration samples as described above. A
frequency model $\hat{\mu}$ is built using the information in the training
sample units $(X_{i},D_{i})$ with $i\in I_{1}$. After this, a severity model
$\hat{\psi}$ is built from $(X_{i},D_{i},Y_{i})$, using only those training
sample units $i\in I_{1}$ for which the observed frequency $D_{i}>0$. Using
the frequency model $\hat{\mu}$, we predict the frequency $\hat{\mu}(X_{i})$
for the calibration sample units $i\in I_{2}$. Now, the severity model
$\hat{\psi}$ is used to predict the total claim cost
$\hat{\psi}(X_{i},\hat{\mu}(X_{i}))$ only for the calibration sample units
$i\in I_{2}$ such that $\hat{\mu}(X_{i})>0$. We denote the indexes of the
sample units in this calibration sub-sample by $I^{+}_{2}$. After these steps,
we compute the conformity scores
$R_{i}=|Y_{i}-\hat{\psi}(X_{i},\hat{\mu}(X_{i}))|$ for the calibration sample
units $i\in I^{+}_{2}$.

As in the original split conformal procedure, the assumed distributional
symmetry of the data sequence implies, for a future random triplet
$(X_{n+1},D_{n+1},Y_{n+1})$ such that $\hat{\mu}(X_{n+1})>0$, the conformity
score $R_{n+1}=|Y_{n+1}-\hat{\psi}(X_{n+1},\hat{\mu}(X_{n+1}))|$ is uniformly
ranked among the ordered calibration conformity scores, and the coverage
property $(\star)$ holds for the conformal prediction set

$C^{(1-\alpha)}(X_{n+1})=[\,\max\,\\{\hat{\psi}(X_{n+1},\hat{\mu}(X_{n+1}))-\hat{r},0\\},\hat{\psi}(X_{n+1},\hat{\mu}(X_{n+1}))+\hat{r}\,],$

if in $(\star)$ we replace $n_{2}$ by $n^{+}_{2}$ (the number of calibration
sample units in $I^{+}_{2}$), and define
$\hat{r}=R_{(\lceil(1-\alpha)(n^{+}_{2}+1)\rceil)}$. A programmatic
description of this two-stage split conformal prediction procedure is given in
Algorithm 1.

##  4 Synthetic and real datasets

###  4.1 Synthetic data

Figure 2: Frequency and severity distributions for the synthetic dataset.
Figure 3: Prediction intervals for the synthetic dataset with nominal
miscoverage level $\alpha=5\%$. The black dots are the observed severity
values.

Following the structure of the two-stage data generating process described in
Section 2, we simulate a size $10,000$ dataset with $10$ independent
predictors $X_{i,1},X_{i,2},\dots X_{i,10}$, each one with $10\times U[0,1]$
distribution. The frequency variable $D_{i}$ is drawn from a mixture, with
symmetric weights, comprising a point mass at zero and a
Poisson($\lambda_{i}$) distribution, with rate parameter
$\lambda_{i}=e^{0.01\times X_{1,i}}$. The purpose of this mixture is to mimic
the behavior of the zero-inflated frequencies commonly found in insurance
claims. The severity variable $Y_{i}$ is sampled conditionally: if $D_{i}=0$,
then $Y_{i}=0$, otherwise, we draw $Y_{i}$ from an exponential distribution
with expectation $4\times e^{X_{i,2}}+\sin(X_{i,3}\times X_{i,4})+5\times
X_{i,5}^{3}$. Figure 2 depicts the distributions of claim frequency ($D_{i}$)
and severity ($Y_{i}$). Claims with zero frequency ($D_{i}=0$) make up a
significant portion ($67.8\%$) of the data.

The synthetic dataset is partitioned into training, calibration, and test
samples, with sizes $5,\\!000$, $2,\\!500$, and $2,\\!500$, respectively. We
use a random forest with $1,\\!000$ trees as the predictive model for the
claims frequency. For the severity model, we consider two alternatives: a
gamma regression and a random forest with $1,\\!000$ trees. In terms of
predictive performance, considering only the positive severity predictions on
the test set, the mean absolute error using the gamma regression and the
random forest are $14,\\!196.73$ and $6,\\!189.51$, respectively.

Using a nominal miscoverage level $\alpha=5\%$, Figure 3 shows the prediction
intervals for $50$ random test sample units. For the option using the gamma
regression severity model, the prediction intervals width is $81,\\!551.56$,
while the random forest severity model gives prediction intervals with width
equal to $29,\\!166,26$. The observed coverages of the prediction intervals
considering the full test sample are $94.93\%$ and $94.87\%$, for the gamma
regression and the random forest severity models, respectively.

###  4.2 Motor third party liability in Belgium

Figure 4: Frequency and severity distributions for the MTPL dataset. Figure 5:
Prediction intervals for the MTPL dataset with nominal miscoverage level
$\alpha=5\%$. The black dots are the observed severity values.

In this example, we use a motor third party liability (MTPL) portfolio dataset
from a Belgium insurance company, related to the year 1997, also analysed by
[5, 13, 14]. This dataset contains information about $163,\\!212$
policyholders. The proportion of zero claims in this dataset is approximately
$88.88\%$. Figure 4 depicts the distributions of claim frequency and severity.
A description of the variables is given in Table 2.

Variable | Description
---|---
Coverage | Type of coverage
Fuel | Type of the motor fuel of the vehicle
Sex | Gender of the policyholder: male or female
Use | Main use of the vehicle: private or work
Fleet | Indicates whether the vehicle is part of a fleet
Ageph | Age of the policyholder (in years)
Power | Horsepower of the motor vehicle
Bm | Bonus malus scale (lower value indicates good claim history)
Lat | Latitude of the center of the municipality where the policyholder resides
Long | Longitude of the center of the municipality where the policyholder resides
NClaims | The number of claims
Amount | The total amount claimed by the policyholder (in euros)
Expo | The duration of the contract (proportion of the year)
Table 2: Variables in the MTPL dataset (reproduced from [5]).

The MTPL dataset is randomly partitioned into training, calibration, and test
samples, according to the proportions $70\%$, $15\%$, and $15\%$,
respectively. Again, we use a random forest with $1,\\!000$ trees as the
predictive model for the claims frequency. As before, we consider two
alternatives for the severity model: a gamma regression and a random forest
with $1,\\!000$ trees. For the severity predictions on the test set, the mean
absolute error using the gamma regression and the random forest are
$1,\\!389.59$ and $628.57$ euros, respectively.

Using a nominal miscoverage level $\alpha=5\%$, Figure 5 shows the prediction
intervals for $50$ random test sample units. For the scenario using the gamma
regression severity model, the prediction intervals average width is
$3,\\!153.94$ euros, while the random forest severity model gives prediction
intervals with average width equal to $1,\\!568.74$ euros. The observed
coverages of the prediction intervals are $96.18\%$ and $97.13\%$, for the
gamma regression and the random forest severity models, respectively.

###  4.3 Crop insurance in Brazil

Now, we consider data related to crop insurance in Brazil. The data was
compiled from the Brazilian Agricultural Ministry public records444The raw
data is publicly available at:
https://dados.agricultura.gov.br/dataset/sisser3, encompassing crop insurance
policies within a time frame ranging from March 2016 to February 2022. All
policies belong to the Rural Grant Program (Programa de Subvenção Rural, PSR),
which is maintained by the Brazilian federal government in order to subsidize
a percentage of the insurance premiums for farmers.

In our analysis we focus on policies related to soybean crops, due to the
peculiarities of its agricultural cycles, and its relevance for the Brazilian
economy. The dataset is aggregated at the municipality level. We consider the
sample unit as the municipality in a given year crop. Therefore, the second
stage component actually refers to an aggregated loss within each
municipality. Furthermore, municipalities without at least one policy during a
specific agricultural year were excluded from the dataset.

The dataset incorporates relevant climate variables with monthly cumulative
precipitation (in millimeters) and monthly average temperatures (in Celsius)
which were obtained from version 4 of the Climatic Research Unit gridded Time
Series [15]. The detailed list of variables of this dataset is presented in
Table 3.

Variable | Description
---|---
Municipality | Name of the municipality
Year | Agricultural year, numbered from 1 to 6
Latitude | Latitude of the municipality
Longitude | Longitude of the municipality
AWC | Available water capacity of the municipality
Soil | Predominant type of soil found in the municipality
TempPC1 | First principal component of temperature
TempPC2 | Second principal component of temperature
PrecPC1 | First principal component of precipitation
PrecPC2 | Second principal component precipitation
PrecPC3 | Third principal component precipitation
PrecPC4 | Fourth principal component of precipitation
Claims | Number of claims
RelativeLoss | Measure of relative loss, ratio of total claims amounts to number of claims
Table 3: Variables in the crop insurance dataset. Figure 6: Frequency and
severity distributions for the crop insurance dataset.

The distributions of frequency and relative loss by municipality is given in
Figure 6. The proportion of claims in this dataset is approximately $32.57\%$.
This substantially high proportion is due to the fact that soybean crops are
highly sensitive to climate factors, especially droughts, which happened in
the considered time frame.

We randomly split the crop insurance dataset into training, calibration, and
test samples, according to the proportions $70\%$, $15\%$, and $15\%$,
respectively. Again, we use a random forest with $1,\\!000$ trees as the
predictive model for the claims frequency. As before, we consider two
alternatives for the severity model: a gamma regression and a random forest
with $1,\\!000$ trees. For the positive relative loss predictions on the test
set, the mean absolute error using the gamma regression and the random forest
are $142,\\!988.00$ and $105,\\!796.80$, respectively.

Using a nominal miscoverage level $\alpha=5\%$, Figure 7 shows the prediction
intervals for $50$ random test sample units. For the scenario using the gamma
regression severity model, the prediction intervals average width is
$719,\\!859.50$, while the random forest relative loss model gives prediction
intervals with width equal to $473,\\!635.00$. The observed coverages of the
prediction intervals are $96.40\%$ and $94.42\%$, for the gamma regression and
the random forest relative loss models, respectively.

Figure 7: Prediction intervals for the crop insurance dataset with nominal
miscoverage level $\alpha=5\%$. The black dots are the observed severity
values.

##  5 Out-of-bag extension

The two-stage split conformal prediction has two undesired features. Firstly,
the need to have a separate calibration sample reduces the amount of data
available to train the frequency and the severity models. Secondly, the
prediction intervals for future observations have exactly the same width when
the lower limits are non negative numbers. When a random forest, or any other
predictive model relying on bagging [16], is used for the severity stage, the
out-of-bag mechanism [9, 17] can be exploited to address both issues: a
calibration sample is no longer needed, and the prediction intervals begin to
have adaptive width. To understand this extension of Algorithm 1 we need to
review a few details about the inner workings of a random forest, in the
context of regression.

We are in the usual supervised learning regression context [18, 19, 20], and a
training sample $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ is available, for which
$x_{i}\in\mathbb{R}^{p}$ is a vector of predictors, and $y\in\mathbb{R}$ is
the response variable. In the 1980s, Breiman and coauthors [21] developed the
classification and regression trees (CART) algorithm, which recursively
partitions the space of predictor variables, greedily looking for splits in
the training data which minimize a quadratic cost function. If we resample the
training set, drawing $n$ observations uniformly with replacement, this
bootstrap sample [22] can be used to train a tall regression tree using the
CART algorithm. By repeating this process $B$ times, we get the regression
trees $\hat{\psi}^{(1)},\dots,\hat{\psi}^{(B)}$, which can be averaged to
produce the aggregated regression function

$\hat{\psi}(\;\cdot\;)=\frac{1}{B}\sum_{j=1}^{B}\hat{\psi}^{(j)}(\;\cdot\;).$

This general aggregation process of regression functions trained from
bootstrap samples, known as bagging [16], can be further optimized by
uniformly drawing without replacement a random subset of predictors when
deciding each split of each regression tree in the ensemble, resulting in the
definition of a random forest [9].

From the classical bias-variance trade-off perspective [18], the effectiveness
of random forests in making accurate predictions can be attributed to two
factors. Firstly, by averaging numerous tall regression trees, which have low
bias and high variance, random forests effectively reduce variance. Secondly,
Breiman’s randomized split decision mechanism further enhances the performance
by breaking up correlations between predictions made by individual trees in
the ensemble.

It is easy to prove that this use of bootstrapping in the training process of
a random forest implies that each training observation is not used when
cultivating approximately $e^{-1}\approx 36.8\%$ of the trees in the ensemble.
Letting $\mathcal{O}_{i}\subset\\{1,2,\dots,B\\}$ denote the indexes of the
trees for which the $i$-th training observation stayed out of the
corresponding bootstrap samples – stayed “out-of-bag”, in Breiman’s
terminology – the regression trees
$\\{\hat{\psi}^{(j)}:j\in\mathcal{O}_{i}\\}$ form a random subforest which can
be used to extend Algorithm 1 in the following manner.

Any frequency model $\hat{\mu}$ is trained from
$\\{(x_{i},d_{i})\\}_{i=1}^{n}$. Let $I^{+}$ denote the indexes of the $n^{+}$
training sample units whose observed frequency $d_{i}$ are not zero. For the
severity model, grow a random forest $\hat{\psi}=\\{\hat{\psi}^{(j)}\\}_{j\in
B}$, with $B$ trees, using the training sub-sample
$\\{(x_{i},d_{i},y_{i})\\}_{i\in I^{+}}$. Compute

$\hat{y}_{i}=\frac{1}{|\mathcal{O}_{i}|}\sum_{j\in\mathcal{O}_{i}}\hat{\psi}^{(j)}(x_{i},\hat{\mu}(x_{i}))\quad\text{and}\quad\hat{\delta}_{i}=\frac{1}{|\mathcal{O}_{i}|}\sum_{j\in\mathcal{O}_{i}}\left|\hat{\psi}^{(j)}(x_{i},\hat{\mu}(x_{i}))-\hat{y}_{i}\right|,$

for $i\in I^{+}$, in which $|\mathcal{O}_{i}|$ denotes the number of trees
whose indexes are in $\mathcal{O}_{i}$. For $i\in I^{+}$, compute the weighted
absolute residuals $r_{i}=\left|y_{i}-\hat{y}_{i}\right|/\hat{\delta_{i}}$.
For a nominal miscoverage level $0<\alpha<1$, letting
$\hat{r}=r_{(\lceil(1-\alpha)(n^{+}\,+\,1)\rceil)}$, the corresponding
prediction interval for a future vector of predictors $x_{n+1}$ is given by

$[\,\max\,\\{\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))-\hat{r}\times\hat{\delta}_{n+1},0\\},\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))+\hat{r}\times\hat{\delta}_{n+1}\,],$

in which
$\hat{\delta}_{n+1}=\frac{1}{B}\sum_{j=1}^{B}\left|\hat{\psi}^{(j)}(x_{n+1},\hat{\mu}(x_{n+1}))-\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))\right|$.
A more programmatic description of this two-stage out-of-bag conformal
prediction procedure is given in Algorithm 2.

Algorithm 2 Two-stage out-of-bag conformal prediction

1:Dataset $\\{(x_{i},d_{i},y_{i})\\}_{i=1}^{n}$, number $B$ of trees used to
train the severity random forest, future vector of predictors
$x_{n+1}\in\mathbb{R}^{p}$, and nominal miscoverage level $0<\alpha<1$.

2:

3:function PredictionInterval($\\{(x_{i},d_{i},y_{i})\\}_{i=1}^{n}$, $B$,
$x_{n+1}$, $\alpha$)

4: Train frequency model $\hat{\mu}$ from $\\{(x_{i},d_{i})\\}_{i=1}^{n}$

5: if $\hat{\mu}(x_{n+1})=0$ then

6: return $\emptyset$

7: end if

8: $I^{+}\leftarrow\\{i\in\\{1,\dots,n\\}:d_{i}>0\\}$

9: $n^{+}\leftarrow|I^{+}|$

10: Train severity random forest $\hat{\psi}=\\{\hat{\psi}^{(j)}\\}_{j\in B}$
from $\\{(x_{i},d_{i},y_{i})\\}_{i\in I^{+}}$

11: for $i\in I^{+}$ do

12: $\mathcal{O}_{i}\leftarrow\emptyset$

13: for $j\leftarrow 1\textrm{ to }B$ do

14: if $i$-th sample unit not in the $j$-th bootstrap sample then

15: $\mathcal{O}_{i}\leftarrow\mathcal{O}_{i}\cup\\{j\\}$

16: end if

17: end for

18:
$\hat{y}_{i}\leftarrow\frac{1}{|\mathcal{O}_{i}|}\sum_{j\in\mathcal{O}_{i}}\hat{\psi}^{(j)}(x_{i},\hat{\mu}(x_{i}))$

19:
$\hat{\delta}_{i}\leftarrow\frac{1}{|\mathcal{O}_{i}|}\sum_{j\in\mathcal{O}_{i}}\left|\hat{\psi}^{(j)}(x_{i},\hat{\mu}(x_{i}))-\hat{y}_{i}\right|$

20: $r_{i}\leftarrow\left|y_{i}-\hat{y}_{i}\right|/\hat{\delta_{i}}$

21: end for

22: $\hat{r}\leftarrow r_{(\lceil(1-\alpha)(n^{+}\,+\,1)\rceil)}$

23:
$\hat{\delta}_{n+1}\leftarrow\frac{1}{B}\sum_{j=1}^{B}\left|\hat{\psi}^{(j)}(x_{n+1},\hat{\mu}(x_{n+1}))-\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))\right|$

24: return
$[\,\max\,\\{\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))-\hat{r}\times\hat{\delta}_{n+1},0\\},\hat{\psi}(x_{n+1},\hat{\mu}(x_{n+1}))+\hat{r}\times\hat{\delta}_{n+1}\,]$

25:end function

In order to elucidate our methodology, we randomly partition the crop
insurance dataset into training and test samples, according to the proportions
$70\%$ and $30\%$, respectively. For both the frequency and relative loss
models, we employ a random forest with $1,\\!000$ trees. Setting a nominal
miscoverage level $\alpha=5\%$, Figure 8 shows the out-of-bag conformal
prediction intervals for $100$ random sample units. The median width of the
prediction intervals is $581,\\!607.70$, whereas the observed coverage
obtained is $96.84\%$.

Figure 8: The out-of-bag prediction intervals for the crop insurance dataset
with nominal miscoverage level $\alpha=10\%$. The black dots are the observed
relative loss values.

##  6 Concluding remarks

The increasing use of machine learning algorithms for the predictive modeling
of frequency-severity tasks demands the development of a comprehensive
framework to assess the confidence in the corresponding forecasts. It is
desirable that this framework avoids relying on strong assumptions about the
underlying data generation process and the specific features of the predictive
algorithms selected for the modeling task. Additionally, it should generate
prediction intervals with reliable finite sample coverage guarantees. This
paper contributes in this direction, proposing conformal prediction procedures
suitable for dealing with the two-stage nature of frequency-severity modeling,
codified in the more general Algorithm 1 and the more specialized Algorithm 2,
which exhibit different features in terms of data usage and properties of the
generated prediction intervals. The general performance exhibited by the
procedures developed in the paper motivates the future investigation of
conformal prediction ideas applied to other actuarial science problems.

## Data availability

Open source R [23] codes and data for all the examples in the paper are openly
available at https://github.com/heltongraziadei/conformal-fs.

## References

  * [1] E. W. Frees, Regression modeling with actuarial and financial applications.  Cambridge University Press, 2009.
  * [2] P. Shi, X. Feng, and A. Ivantsova, “Dependent frequency–severity modeling of insurance claims,” Insurance: Mathematics and Economics, vol. 64, pp. 417–428, 2015.
  * [3] C. E. McCulloch, “Generalized linear models,” Journal of the American Statistical Association, vol. 95, no. 452, pp. 1320–1324, 2000.
  * [4] X. Su and M. Bai, “Stochastic gradient boosting frequency-severity model of insurance claims,” PloS one, vol. 15, no. 8, p. e0238000, 2020.
  * [5] R. Henckaerts, M.-P. Côté, K. Antonio, and R. Verbelen, “Boosting insights in insurance tariff plans with tree-based machine learning methods,” North American Actuarial Journal, vol. 25, no. 2, pp. 255–285, 2021.
  * [6] V. Vovk, A. Gammerman, and C. Saunders, “Machine-learning applications of algorithmic randomness,” in Proceedings of the Sixteenth International Conference on Machine Learning, ICML ’99, (San Francisco, CA, USA), pp. 444–453, Morgan Kaufmann Publishers Inc., 1999.
  * [7] V. Vovk, A. Gammerman, and G. Shafer, Algorithmic learning in a random world.  Springer Science & Business Media, 2005.
  * [8] J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman, “Distribution-free predictive inference for regression,” Journal of the American Statistical Association, vol. 113, no. 523, pp. 1094–1111, 2018\.
  * [9] L. Breiman, “Random Forests,” Machine Learning, vol. 45, no. 1, pp. 5–32, 2001.
  * [10] H. Papadopoulos, K. Proedrou, V. Vovk, and A. Gammerman, “Inductive confidence machines for regression,” in Machine Learning: ECML 2002 (T. Elomaa, H. Mannila, and H. Toivonen, eds.), (Berlin, Heidelberg), pp. 345–356, Springer Berlin Heidelberg, 2002.
  * [11] G. Shafer and V. Vovk, “A tutorial on conformal prediction,” Journal of Machine Learning Research, vol. 9, no. 12, pp. 371–421, 2008.
  * [12] M. Fontana, G. Zeni, and S. Vantini, “Conformal prediction: A unified review of theory and new challenges,” Bernoulli, vol. 29, no. 1, pp. 1–23, 2023\.
  * [13] M. Denuit and S. Lang, “Non-life rate-making with bayesian gams,” Insurance: Mathematics and Economics, vol. 35, no. 3, pp. 627–647, 2004.
  * [14] N. Klein, M. Denuit, S. Lang, and T. Kneib, “Nonlife ratemaking and risk management with bayesian generalized additive models for location, scale, and shape,” Insurance: Mathematics and Economics, vol. 55, pp. 225–249, 2014\.
  * [15] I. Harris, T. J. Osborn, P. Jones, and D. Lister, “Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset,” Scientific data, vol. 7, no. 1, p. 109, 2020.
  * [16] L. Breiman, “Bagging predictors,” Machine Learning, vol. 24, no. 2, pp. 123–140, 1996.
  * [17] P. C. Marques F., “Confidence intervals for the random forest generalization error,” Pattern Recognition Letters, vol. 158, pp. 171–175, 2022.
  * [18] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, 2nd ed., 2009.
  * [19] K. P. Murphy, Machine Learning: A Probabilistic Perspective.  MIT Press, 2012.
  * [20] C. M. Bishop, Pattern Recognition and Machine Learning.  Springer, 2006.
  * [21] L. Breiman, J. H. Friedman, R. A. Olshen, and C. J. Stone, Classification and Regression Trees.  Routledge, 1984.
  * [22] B. Efron, “Bootstrap methods: another look at the jackknife,” The Annals of Statistics, vol. 7, no. 1, pp. 1–26, 1979.
  * [23] R Core Team, R: a language and environment for statistical computing.  R Foundation for Statistical Computing, Vienna, Austria, 2017.


# AutoWeird: Weird Translational Scoring Function Identified by Random Search

Hansi Yang
Department of Computer Science and Engineering
Hong Kong University of Science and Technology
Clear Water Bay, Hong Kong
<EMAIL_ADDRESS>
&Yongqi Zhang
4Paradigm Inc.
Beijing, China
<EMAIL_ADDRESS>
&Quanming Yao
Department of Electronic Engineering
Tsinghua University
Beijing, China
<EMAIL_ADDRESS>

###### Abstract

Scoring function (SF) measures the plausibility of triplets in knowledge
graphs. Different scoring functions can lead to huge differences in link
prediction performances on different knowledge graphs. In this report, we
describe a weird scoring function found by random search on the open graph
benchmark (OGB). This scoring function, called AutoWeird, only uses tail
entity and relation in a triplet to compute its plausibility score.
Experimental results show that AutoWeird achieves top-1 performance on ogbl-
wikikg2 data set, but has much worse performance than other methods on ogbl-
biokg data set. By analyzing the tail entity distribution and evaluation
protocol of these two data sets, we attribute the unexpected success of
AutoWeird on ogbl-wikikg2 to inappropriate evaluation and concentrated tail
entity distribution. Such results may motivate further research on how to
accurately evaluate the performance of different link prediction methods for
knowledge graphs.

##  1 Introduction

Scoring function (SF), which measures the plausibility of triplets in
knowledge graphs (KGs), plays an important role to learn from KGs [13, 7]. To
fairly evaluate different SFs, open graph benchmark (OGB) [5] sets up
leaderboards to evaluate the different models on large-scale graphs. Among the
two leaderboards for link prediction task in KG, we observe that the leading
models on ogbl-wikikg2 have similar forms by modeling embeddings as
translations in vector spaces. To analyze the success of these SFs and try to
design a better SF, we are motivated to use automated machine learning
(AutoML) techniques [14, 17, 18] to search SFs in this type.

Based on the translational models [1, 12, 15, 2, 16], we first set up a search
space that can cover all these translational models as special cases. Then, we
conduct random search in this space to explore new designs of scoring
functions. After evaluating $10$ SFs, we identified a weird SF which does not
contain head entity embedding but outperforms the existing models. This is
abnormal as triplets in KGs are formed by head entity, relation and tail
entity. By analyzing the distribution of tail entities, we guess the success
of such a weird SF may attribute to the concentrated distribution and the
evaluation protocol with random negative sampling.

This tech-report is organized as follows. First, we introduce the search
space, search algorithm and the training objective in Section 2. Then in
Section 3, we show the the top-1 performance achieved by the weird
translation-based SF on ogbl-wikikg2 and our analysis on the searched results.

##  2 Method

The scoring function used here is found by automated machine learning
techniques (AutoML) [6, 14]. We first describe the search space where we find
the scoring function (SF) in Section 2.1. And our search algorithm is
introduced in Section 2.2.

###  2.1 Search space

Our search space is motivated by translation-based models [1, 12, 15, 2, 16]
on the ogbl-wikikg2 leaderboard
111https://ogb.stanford.edu/docs/leader_linkprop/#ogbl-wikikg2. Our search
space is designed as follows:

  * •

We set up the entity embeddings with two parts ${\bm{e}}_{0},{\bm{e}}_{1}$,
relation embeddings with three parts ${\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2}$.
All these vectors have same dimension $d$, i.e.
${\bm{e}}_{0},{\bm{e}}_{1},{\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2}\in\mathbb{R}^{d}$.

  * •

The scoring function is given by combining the the different components, i.e.,
$s(h,r,t)=-\|f({\bm{e}}_{0}^{h},{\bm{e}}_{1}^{h},{\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{e}}_{0}^{t},{\bm{e}}_{1}^{t})\|_{1}$,
where ${\bm{e}}_{0}^{h},{\bm{e}}_{1}^{h}$ denote the embeddings for head
entity $h$, ${\bm{e}}_{0}^{t},{\bm{e}}_{1}^{t}$ denote the embeddings for tail
entity $t$, ${\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2}$ denote the embeddings for
relation $r$, and $\|\cdot\|_{1}$ is the $\ell_{1}$ norm.

  * •

The function $f(\cdot)$ outputs a vector in $\mathbb{R}^{d}$, and is computed
by adding/subtracting first-order terms
(${\bm{e}}_{0}^{h},{\bm{e}}_{1}^{h},{\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{e}}_{0}^{t},{\bm{e}}_{1}^{t}$)
and/or second-order terms (element-wise product between any of these two
vectors), e.g., $\bm{e}_{0}^{h}\circ\bm{e}_{1}^{t}$ and
$\bm{e}_{0}^{h}\circ\bm{r}_{1}$, for entity and relation embeddings.

We can see that many translation-based SFs that are already on the OGB
leaderboard [5] can all be covered in this search space. Some examples are
listed in Table 1.

Table 1: Example knowledge graph embedding models that can be covered by our search space. Model | Scoring Function
---|---
TransE [1] | $-\|{\bm{e}}_{0}^{h}-{\bm{e}}_{0}^{t}+{\bm{r}}_{0}\|$
InterHT [12] | $-\|{\bm{e}}_{0}^{h}\circ{\bm{e}}_{1}^{t}-{\bm{e}}_{1}^{h}\circ{\bm{e}}_{0}^{t}+{\bm{r}}_{0}\|$
TripleRE [15] | $-\|{\bm{e}}_{0}^{h}\circ{\bm{r}}_{1}-{\bm{e}}_{0}^{t}\circ{\bm{r}}_{2}+{\bm{r}}_{0}\|$
PairRE [2] | $-\|{\bm{e}}_{0}^{h}\circ{\bm{r}}_{1}-{\bm{e}}_{0}^{t}\circ{\bm{r}}_{2}\|$
TranS [16] | $-\|{\bm{e}}_{0}^{h}\circ{\bm{e}}_{1}^{t}-{\bm{e}}_{1}^{h}\circ{\bm{e}}_{0}^{t}+{\bm{r}}_{0}+{\bm{e}}_{0}^{h}\circ{\bm{r}}_{1}+{\bm{e}}_{0}^{t}\circ{\bm{r}}_{2}\|$
AutoWeird | $-\|-{\bm{e}}_{1}^{t}\circ{\bm{r}}_{2}+{\bm{e}}_{0}^{t}\circ{\bm{r}}_{0}+{\bm{e}}_{0}^{t}\circ{\bm{r}}_{2}-{\bm{r}}_{0}\|$

This search space covers a large number of different scoring functions. Note
that with 7 different vectors
${\bm{e}}_{0}^{h},{\bm{e}}_{1}^{h},{\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2},{\bm{e}}_{0}^{t},{\bm{e}}_{1}^{t}$,
the total number of first-order and second-order terms is $7+7\times 7=56$.
And for each of these 56 terms, we have three possible choices for its
coefficient: 1, -1 and 0, where 0 coefficient means that it does not appear in
$f(\cdot)$. Therefore, the total number of possible scoring functions can be
as large as $3^{56}\approx 5.23\times 10^{26}$.

###  2.2 Search algorithm

To find new scoring functions from this search space, we consider the
following bi-level optimization problem:

$\displaystyle\min_{\bm{\alpha}}\mathcal{M}_{\text{val}}\left(\\{{\bm{e}}(\bm{\alpha}),{\bm{r}}(\bm{\alpha})\\}\right)\quad\text{s.t.}\quad\\{{\bm{e}}(\bm{\alpha}),{\bm{r}}(\bm{\alpha})\\}=\min_{\\{{\bm{e}},{\bm{r}}\\}}\mathcal{L}_{\text{tr}}(\\{{\bm{e}},{\bm{r}}\\};s_{\bm{\alpha}}),$

where $\bm{\alpha}$ denotes the hyper-parameters that specify the SF
$s_{\bm{\alpha}}$, $\mathcal{M}_{\text{val}}$ refers to the validation
performance (e.g., MRR on the validation set), and $\mathcal{L}_{\text{tr}}$
refers to the training loss. We use the self-adversarial negative sampling
loss [9], which is defined as follows:

$\displaystyle\mathcal{L}(\\{{\bm{e}},{\bm{r}}\\};s_{\bm{\alpha}})=-\log\
\sigma(\gamma-s_{\bm{\alpha}}(h,r,t))-\sum_{i=1}^{n}\frac{1}{k}\log\
\sigma(s_{\bm{\alpha}}(h^{\prime}_{i},r,t^{\prime}_{i})-\gamma),$ (1)

where $\gamma$ is a fixed margin, $\sigma$ is the sigmoid function, and
$(h^{\prime}_{i},r,t^{\prime}_{i})$ is the $i$-th negative triplet.

For simplicity, we use random search to solve this problem. We first restrict
the $f(\cdot)$ function to be computed by 5 first-order or second-order terms,
and randomly choose 4 terms among all 56 terms. Then for each term, we
randomly assign a coefficient 1 or -1 for it. The final function $f(\cdot)$ is
computed by summing over these terms multiplied by their coefficients. We
repeat the random search for 10 times. The searched SF, referred as AutoWeird,
is also shown in Table 1. Different with all previous scoring functions, it
only uses information from the relation ($\bm{r}_{0},\bm{r}_{2}$) and tail
entity ($\bm{e}_{0}^{t}$ and $\bm{e}_{1}^{t}$) to compute the score.
Intuitively, it should have bad performances as it does not consider the head
entity. However, as will be demonstrated in experimental parts, it can achieve
state-of-the-art performances on some KG data sets.

##  3 Experiments

In this section, we will first introduce the dataset. Then we will introduce
implementation details, experimental results, and the post-analysis on the
results.

###  3.1 Dataset and metric

The ogbl-wikikg2 [5] is a large KG dataset extracted from Wikidata [11]. It
contains a set of triplet edges, capturing the different types of relations
between entities in the world. There are $2,500,604$ entities, $535$ relation
types and $17,137,181$ triplets in the whole data set.

We use original dataset splitting configuration. It splits the triplets
according to time, simulating a realistic KG completion scenario that aims to
fill in missing triplets that are not present at a certain timestamp. The
training set contains 16,109,182 triplets, the validation set contains 429,456
triplets, and the test set contains 598,543 triplets.

Following the official guideline, we evaluate the KG embedding performance by
predicting new triplet edges according to the training edges. The evaluation
metric follows the standard filtered metric [1] which is widely used in KG.
Specifically, each test triplet edges are corrupted by replacing its head or
tail with randomly-sampled 500 negative entities, while ensuring the resulting
triplets do not appear in KG. The goal is to rank the true head (or tail)
entities higher than the negative entities, which is measured by Mean
Reciprocal Rank (MRR).

###  3.2 Implementation details

In our experiments, the Adam [8] optimizer is used with a learning rate of
$0.0005$. The batch size of the model is set to 512. To prevent overfitting,
we use dropout technique, and set it to 0.1. The negative sampling size is set
to 128. And the dimension $d$ of each embedding vector
${\bm{e}}_{0},{\bm{e}}_{1},{\bm{r}}_{0},{\bm{r}}_{1},{\bm{r}}_{2}\in\mathbb{R}^{d}$
is set to 200. The maximum number of training steps is 300,000. We validate
the model every 20 thousand steps. The number of anchors for NodePiece are
20,000. And $\gamma$ in the loss function is set to 6. These hyper-parameters
are selected according to the performance on the validation set.

###  3.3 Results

The results on ogbl-wikikg2 are shown in Table 2. By only using information
from relation and tail entity, the AutoWeird model still achieves an MRR of
0.7368 on validation set and 0.7368 on test set, which outperforms all
previous models on the ogbl-wikikg2 dataset. The revised version of AutoWeird
model also has a performance better than all the previous methods. These
counter-intuitive results demonstrate some hidden problems on the evaluation
of KG embedding models. Hence, we conduct a post-analysis in the next
subsection.

Table 2: Results on the ogbl-wikikg2 dataset. Model | #Params | #Dims | Test MRR | Valid MRR
---|---|---|---|---
TransE [1] | 1251M | 500 | 0.4256 $\pm$ 0.0030 | 0.4272 $\pm$ 0.0030
RotatE [9] | 1250M | 250 | 0.4332 $\pm$ 0.0025 | 0.4353 $\pm$ 0.0028
PairRE [2] | 500M | 200 | 0.5208 $\pm$ 0.0027 | 0.5423 $\pm$ 0.0020
AutoSF [17] | 500M | - | 0.5458 $\pm$ 0.0052 | 0.5510 $\pm$ 0.0063
ComplEx [10] | 1251M | 250 | 0.5027 $\pm$ 0.0027 | 0.3759 $\pm$ 0.0016
TripleRE [15] | 501M | 200 | 0.5794 $\pm$ 0.0020 | 0.6045 $\pm$ 0.0024
ComplEx-RP [3] | 250M | 50 | 0.6392 $\pm$ 0.0045 | 0.6561 $\pm$ 0.0070
AutoSF + NodePiece [4] | 6.9M | - | 0.5703 $\pm$ 0.0035 | 0.5806 $\pm$ 0.0047
TripleREv2 + NodePiece [15] | 7.3M | 200 | 0.6582 $\pm$ 0.0020 | 0.6616 $\pm$ 0.0018
InterHT + NodePiece [12] | 19.2M | 200 | 0.6779 $\pm$ 0.0018 | 0.6893 $\pm$ 0.0015
TripleREv3 + NodePiece [15] | 36.4M | 200 | 0.6866 $\pm$ 0.0014 | 0.6955 $\pm$ 0.0008
TranS + NodePiece [16] | 19.2M | 200 | 0.6882 $\pm$ 0.0019 | 0.6988 $\pm$ 0.0006
AutoWeird + NodePiece | 19.2M | 200 | 0.7353 $\pm$ 0.0006 | 0.7362 $\pm$ 0.0006
EntOccur | - | - | 0.4419 | 0.4467
Table 3: Results on the ogbl-biokg dataset. Model | #Params | #Dims | Test MRR | Valid MRR
---|---|---|---|---
TransE [1] | 188M | 500 | 0.7452 $\pm$ 0.0004 | 0.7456 $\pm$ 0.0003
RotatE [9] | 188M | 250 | 0.7989 $\pm$ 0.0004 | 0.7997 $\pm$ 0.0002
ComplEx [10] | 188M | 250 | 0.8095 $\pm$ 0.0007 | 0.8105 $\pm$ 0.0001
PairRE [2] | 188M | 200 | 0.8164 $\pm$ 0.0005 | 0.8172 $\pm$ 0.0005
AutoSF [17] | 93.8M | - | 0.8309 $\pm$ 0.0008 | 0.8317 $\pm$ 0.0007
TripleRE [15] | 470M | 200 | 0.8348 $\pm$ 0.0007 | 0.8360 $\pm$ 0.0006
ComplEx-RP [3] | 188M | 1000 | 0.8492 $\pm$ 0.0002 | 0.8497 $\pm$ 0.0002
AutoBLM+KGBench [18] | 192M | - | 0.8536 $\pm$ 0.0003 | 0.8548 $\pm$ 0.0002
AutoWeird | 376M | 500 | 0.6755 $\pm$ 0.0004 | 0.6752 $\pm$ 0.0003
EntOccur | - | - | 0.3387 | 0.3387

###  3.4 Post-analysis: why AutoWeird is effective?

To analyze why AutoWeird is effective, we first count the occurrence of tail
entities after adding inverse relations to the data set. Specifically, for
each entity, we count how many times it appears as the tail entity 222 For
simplicity, we regard the head prediction as tail prediction by adding inverse
relations, i.e., $(t,r_{\text{inv}},h)$ for $(h,r,t)$ in the datasets. in the
training set. As is shown in Figure 1(a), a few number of entities has more
than one million occurrences in the training set. Such a concentrated
distribution also exists in the test triplets (Figure 1(b)), leading to a
strong bias to the high occurrence entities.

Note that the evaluation protocol in ogbl-wikikg2 is to rank the positive tail
entities over a set of 500 randomly sampled negative entities. Hence, the bias
will become more severe during evaluation. As in Figure 1(a) and 1(b), if we
only randomly sample $500$ entities from $10^{6}$ entities, the majority of
them will has low occurrence.

(a) Training set of ogbl-wikikg2.  (b) Testing set of ogbl-wikikg2.  (c)
Training set of ogbl-biokg.  (d) Testing set of ogbl-biokg.

Figure 1: Tail entity occurrence distribution for different KG data sets.

Based on these observations, we design a simple method, named as EntOccur, by
using the occurrence as scores. For each triplet $(h,r,t)$ in the training
set, we count the occurrence of entities $t$ for each relation $r$. Then the
score of $(h,r,t)$ in inference is directly set as the occurrence of $t$ with
respect to the specific relation $r$. As in Table 2, such a simple model
EntOccur outperforms some embedding models like TransE and RotatE.

We further evaluate the searched SF AutoWeird and the simple solution EntOccur
on the other KG link prediction dataset ogbl-biokg [5]. By comparing the
performance of AutoWeird and EntOccur with the results on the leaderboard
333https://ogb.stanford.edu/docs/leader_linkprop/#ogbl-biokg in Table 3,
AutoWeird and EntOccur are much worse. By comparaing ogbl-biokg with ogbl-
wikikg2, we find that

  * •

the distributions of tail entity occurrence in ogbl-biokg in Figure 1(c) and
1(d) are not as concentrated as those of ogbl-wikikg2 in Figure 1(a) and 1(b);

  * •

the number of entities in ogbl-biokg is $93,773$, much smaller than ogbl-
wikikg2 with $2,500,604$;

  * •

the negative samples for evaluation in ogbl-biokg have the same type with the
positive tail entity.

Hence, the negative entities in ogbl-biokg have higher correlation with the
positive tail entity than ogbl-wikikg2, alleviating the bias towards entities
with high occurrence. As such, we attribute the success of AutoWeird in ogbl-
wikikg2 to the inappropriate evaluation protocol.

In conclusion, these results show that there is a relatively strong
correlation between relations and a few high occurrence entities in ogbl-
wikikg2. And the correlation may lead to the weird scoring function, AutoWeird
in Table 1. A better way to compare the different models may be (1) using the
full set of negative entities, or (2) using some highly correlated negative
entities, or (3) using the larger number of negative entities rather than
$500$ randomly sampled ones. However, due to the time complexity in evaluating
millions of negative entities for each triplet in ogbl-wikikg2, we leave it as
a future work.

##  4 Conclusion

In this paper, we propose a novel search space for the scoring function in
knowledge graph embedding model. We find a novel scoring function from this
search space by random search, which only depends on tail entity and relation.
Empirical results on the ogbl-wikikg2 dataset demonstrate that the searched
scoring function is able to achieve state-of-the-art results among all
previous models. These strange results should motivate us to consider
potential problems on the evaluation of KG embedding models.

## References

  * [1] A. Bordes, N. Usunier, A. Garcia-Duran, J. Weston, and O. Yakhnenko.  Translating embeddings for modeling multi-relational data.  Advances in neural information processing systems, 26, 2013.
  * [2] L. Chao, J. He, T. Wang, and W. Chu.  Pairre: Knowledge graph embeddings via paired relation vectors.  arXiv preprint arXiv:2011.03798, 2020.
  * [3] Y. Chen, P. Minervini, S. Riedel, and P. Stenetorp.  Relation prediction as an auxiliary training objective for improving multi-relational graph representations.  arXiv preprint arXiv:2110.02834, 2021.
  * [4] M. Galkin, J. Wu, E. Denis, and W. L. Hamilton.  Nodepiece: Compositional and parameter-efficient representations of large knowledge graphs.  arXiv preprint arXiv:2106.12144, 2021.
  * [5] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  arXiv preprint arXiv:2005.00687, 2020.
  * [6] F. Hutter, L. Kotthoff, and J. Vanschoren.  Automated machine learning: methods, systems, challenges.  Springer Nature, 2019.
  * [7] S. Ji, S. Pan, E. Cambria, P. Marttinen, and S. Y. Philip.  A survey on knowledge graphs: Representation, acquisition, and applications.  IEEE Transactions on Neural Networks and Learning Systems, 2021\.
  * [8] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [9] Z. Sun, Z.-H. Deng, J.-Y. Nie, and J. Tang.  Rotate: Knowledge graph embedding by relational rotation in complex space.  arXiv preprint arXiv:1902.10197, 2019.
  * [10] T. Trouillon, C. R. Dance, Éric Gaussier, J. Welbl, S. Riedel, and G. Bouchard.  Knowledge graph completion via complex tensor factorization.  Journal of Machine Learning Research, 18(130):1–38, 2017.
  * [11] D. Vrandečić and M. Krötzsch.  Wikidata: a free collaborative knowledgebase.  Communications of the ACM, 57(10):78–85, 2014.
  * [12] B. Wang, Q. Meng, Z. Wang, D. Wu, W. Che, S. Wang, Z. Chen, and C. Liu.  Interht: Knowledge graph embeddings by interaction between head and tail entities.  arXiv preprint arXiv:2202.04897, 2022.
  * [13] Q. Wang, Z. Mao, B. Wang, and L. Guo.  Knowledge graph embedding: A survey of approaches and applications.  TKDE, 29(12):2724–2743, 2017.
  * [14] Q. Yao and M. Wang.  Taking human out of learning applications: A survey on automated machine learning.  arXiv preprint arXiv:1810.13306, 2018.
  * [15] L. Yu, Z. Luo, D. Lin, H. Liu, and Y. Deng.  Triplere: Knowledge graph embeddings via triple relation vectors.  arXiv preprint arXiv:2112.0095, 2021.
  * [16] X. Zhang, Q. Yang, and D. Xu.  Trans: Transition-based knowledge graph embedding with synthetic relation representation.  arXiv preprint arXiv:2204.08401, 2022.
  * [17] Y. Zhang, Q. Yao, W. Dai, and L. Chen.  Autosf: Searching scoring functions for knowledge graph embedding.  In 2020 IEEE 36th International Conference on Data Engineering (ICDE), pages 433–444. IEEE, 2020.
  * [18] Y. Zhang, Q. Yao, and J. T. Kwok.  Bilinear scoring function search for knowledge graph learning.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022\.


# A Comparison of Virtual Analog Modelling Techniques for Desktop and Embedded
Implementations

###### Abstract

We develop a virtual analog model of the Klon Centaur guitar pedal circuit,
comparing various circuit modelling techniques. The techniques analyzed
include traditional modelling techniques such as nodal analysis and Wave
Digital Filters, as well as a machine-learning technique using recurrent
neural networks. We examine these techniques in the contexts of two use cases:
an audio plug-in designed to be run on a consumer-grade desktop computer, and
a guitar pedal-style effect running on an embedded device. Finally, we discuss
the advantages and disdvantages of each technique for modelling different
circuits, and targeting different platforms.

##  1 Introduction

The Klon Centaur is an overdrive guitar pedal designed by Bill Finnegan in the
early 1990’s, that has developed cult acclaim amongst guitarists [1]. The
circuit is notable for producing “transparent distortion” [2], a term used to
describe the way the pedal seems to add distortion to a guitar’s sound without
otherwise affecting the tone. While the original manufacturing run of the
pedal ended in 2004, many “clones” of the pedal have been produced by other
manufacturers, adding to its cult following.

Circuit modelling is typically broken down into “white-box” and “black-box”
approaches [3]. A “white-box” approach uses knowledge of the internal
mechanisms of the circuit, often modelling the physical interactions of the
electrical components. Popular white-box methods include nodal analysis [4],
Port-Hamiltonian analysis [5], Wave Digital Filters [6, 7], and nonlinear
state space analysis [8].

“Black-box” circuit modelling methods generally use measurements taken from
the circuit being modelled and attempt to model the response of the circuit
without knowledge of the internal workings of the system. Traditional black-
box techniques include impulse response measurements [9] and extensions
thereof, including the Weiner-Hammerstein method [3]. Recently, researchers
have begun using machine learning methods for black-box modelling. Damskägg
et. al. model several guitar distortion circuits, using a WaveNet style
architecture to generate an output signal sample-by-sample [10]. Parker et.
al. use deep fully-connected networks to approximate nonlinear state-space
solutions for the Korg MS-20 filter circuit, effectively a “grey-box” approach
[11]. Finally, Wright et. al. use a recurrent neural network to model the
behavior of guitar distortion circuits with control parameters [12].

The structure of the paper will be as follows: in §2 we give background
information on circuit modelling using nodal analysis and Wave Digital
Filters. §3 describes the use of recurrent neural networks for circuit
modelling, and outlines the model and training process used for emulating the
“Gain Stage” circuit from the Klon Centaur. In §4 we discuss the real-time
implementation of a complete emulation of the Klon Centaur pedal using the
methods outlined in the previous sections. §5 shares the results of Klon
Centaur emulation as well as recommendations for circuit modelling using the
methods discussed here.

##  2 Traditional Circuit Modelling Techniques

First, we examine the use of traditional circuit modelling techniques,
specifically nodal analysis and Wave Digital filters, using sub-circuits from
the Klon Centaur as examples.

###  2.1 Nodal Analysis

The process for creating a digital model of a circuit using nodal analysis is
as follows:

  1. 1.

Convert the circuit into the Laplace domain.

  2. 2.

Form a Laplace domain transfer function of the circuit.

  3. 3.

Use a conformal map to transform the circuit into the digital domain.

As an example circuit, we examine the Tone Control circuit from the Klon
Centaur (see fig. 1).

$-$$+$+4.5VR22VinR21RV2C14R24R23Vout Figure 1: Klon Centaur Tone Control
Circuit

The first step is to convert the circuit into the Laplace domain, using the
Laplace variable $s=j\omega$. The impedances for each principle circuit
component: resistors ($Z_{R}$), capacitors ($Z_{C}$), and inductors ($Z_{L}$),
are as follows:

$Z_{R}=R,\quad Z_{C}=\frac{1}{Cs},\quad Z_{L}=Ls$ (1)

From there, using linear circuit theory, one can construct a Laplace Domain
transfer function for the circuit. Note that this assumes an ideal operational
amplifier operating in its linear region. For more information on this
process, see [13]. For the tone control circuit, the Laplace domain transfer
function can be written as:

$\frac{V_{out}(s)}{V_{in}(s)}={\scriptscriptstyle\frac{C_{14}\left(\frac{1}{R_{22}}+\frac{1}{R_{21}+R_{v2b}}\right)s+\frac{1}{R_{22}}\left(\frac{1}{R_{21}+R_{v2b}}+\frac{1}{R_{23}+R_{v2a}}\right)}{C_{14}\left(\frac{1}{R_{23}+R_{v2a}}+\frac{1}{R_{24}}\right)s+\frac{-1}{R_{24}}\left(\frac{1}{R_{21}+R_{v2b}}+\frac{1}{R_{23}+R_{v2a}}\right)}}$
(2)

Note that we refer to the the section of potentiometer $R_{v2}$ that is above
the wiper as $R_{v2a}$, and the section below as $R_{v2b}$, and that we ignore
the DC offset created by the $4.5$V voltage source at the positive terminal of
the op-amp.

Next we use a conformal map to transform the transfer function from the
Laplace domain to the z-plane where it can be implemented as a digital filter.
The most commonly used conformal map is the bilinear transform, defined as

$s\leftarrow\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$ (3)

Where $T$ is the sample period of the digital system. For more information on
the use of the bilinear transform to digitize an analog system, see [14]. The
resulting filter is known as a “high-shelf” filter, that accentuates high
frequency content in the signal. The resulting frequency response of the
digital model, validated against the response of the analog circuit is shown
in fig. 2.

Figure 2:  Tone control frequency response at various values of the Treble
parameter, comparing the responses of the analog filter with the digital
model.

####  2.1.1 Advantages and Limitations

The advantages of nodal analysis are that the circuit model is simple and
computationally efficient. The model can be constructed with minimal knowledge
of circuit theory, and a basic understanding of digitial signal processing.
The main disdvantage is that nodal analysis cannot be used for nonlinear
circuits, though it can be extended to model this class of circuits through
modified nodal analysis (MNA) [15]. Another disadvantage of nodal analysis-
based methods is that (typically) large portions of the system need to be
recomputed when a circuit element is changed, such as a potentiometer. While
this computation is fairly simple in the example shown here, it can become
vastly more difficult for more complex systems.

###  2.2 Wave Digital Filters

The Wave Digital Filter (WDF) formalism allows circuits to be modelled in
modular and flexible manner. Originally developed by Alfred Fettweis in the
1970’s [6], WDFs have recently gained popularity in modelling audio circuits,
and have been extended to model a wider class of circuits [7]. The WDF
formalism defines each circuit element as a port with some characteristic
resistance $R_{0}$, and uses wave variables passing through each port, rather
than the typical voltage and current variables. The incident wave at a certain
port is defined as:

$a=v+R_{0}i$ (4)

where $v$ is the voltage across the port, and $i$ is the current passing
through the port. The reflected wave is similarly defined as:

$b=v-R_{0}i$ (5)

A Wave Digital Filter defines circuit elements (resistors, capacitors,
inductors, etc.) in the wave domain, and allows the elements to be connected
by series and parallel adaptors also defined in the wave domain. The full
derivation of these WDF elements is given in [6] and [7].

C3R7R19$-$$+$4.5VC16 Figure 3: Klon Centaur Feed-Forward Network 1 Circuit
$V_{in}$$\mathcal{S}_{1}$$C_{3}$$\mathcal{S}_{2}$$R_{7}$$\mathcal{P}_{1}$$C_{16}$$\mathcal{S}_{3}$$R_{19}$$V_{4.5}$
Figure 4: WDF tree for the Klon Centaur Feed-Forward Network 1 Circuit.
$\mathcal{S}$ and $\mathcal{P}$ nodes refer to series and parallel adaptors
respectively.

Once each circuit element and adaptor has been defined, they are connected
together in a structure often referred to as a WDF tree. As an example, we
examine the “feed-forward network 1” from the Klon Centaur circuit (see fig.
3). The corresponding WDF tree is shown in fig. 4. Simulation results compared
to the analog reference are shown in fig. 5.

Figure 5: WDF simulation results compared to the analog reference.

####  2.2.1 Advantages and Limitations

The primary advantage of the Wave Digital approach is its modularity. The
ability to construct a circuit model with each circuit component treated
completely independently in the digital domain opens up many interesting
possibilities for circuit prototyping, modelling circuit-bent instruments, and
more. Additionally, this modularity allows each circuit component to be
discretized separately, even using different conformal maps, which can improve
model behavior for certain classes of circuits (see [3]). Finally, the
separability of components means that when a component is changed (e.g. a
potentiometer), the component change is propagated so that only components
with behavior that depends on the impedance of the changed component need to
be recomputed.

The main disdvantage of WDFs is their difficulty in handling circuits with
complex topologies or multiple nonlinearities. While the recent addition of
$\mathcal{R}$-type adaptors to the Wave Digital formalism [7] has begun to
make these circuits tractable, the WDF models of these types of circuits are
significantly more computationally complex. Further, the use of
$\mathcal{R}$-type adaptors can somewhat compromise the modularity that makes
WDFs advantageous in the first place.

##  3 Recurrent Neural Network Model

While several styles of machine-learning based models have been used for
modelling analog audio circuitry [10, 11, 16], we choose the recurrent neural
network approach developed in [12] as our starting point. Using a recurrent
neural network (RNN) allows the for a significantly smaller neural network
than would be possible with a traditional deep neural network or convolutional
neural network, meaning that the network can be evaluated much faster for
real-time use, while maintaining a smaller memory footprint (an important
advantage on embedded platforms). Additionally, recurrent neural networks are
a sensible candidate for modelling distortion circuits, particularly circuits
with stateful behavior, given the fact that recurrent network building blocks,
such as gated recurrent units, themselves resemble audio distortion effects
and can directly be used as such [17].

In the following paragraphs, we outline the use of an RNN for modelling the
gain stage circuit from the Klon Centaur pedal. While our model is similar to
the model used in [12], it differs in some notable ways. For instance, the
model described in [12] accepts the values of control parameters to the
circuit as inputs to the RNN, however we were unable to successfully train a
network in this fashion. Instead, we construct separate networks for five
different values of the “Gain” parameter, and fade between the outputs of the
networks in real-time in the final implementation of the model. Other
differences are outlined further below.

###  3.1 Model Architecture

The model architecture described in [12] consists of a single recurrent layer
followed by a fully connected layer consisting of a single “neuron” (see fig.
6). In our models, we use a recurrent layer made up of 8 Gated Recurrent
Units. For training, all models are implemented in Python using the Keras
framework [18].

Input $x[n]$Recurrent LayerCurrent State $h[n]$$z^{-1}$Previous State
$h[n-1]$Fully Connected LayerOutput $y[n]$ Figure 6: RNN Architecture.

####  3.1.1 Recurrent Layer

Recurrent layers are typically comprised of one of two types of recurrent
units: Long Short-Term Memory units (LSTMs) or Gated Recurrent Units (GRUs).
For this application, we choose to use GRUs [19] since they require fewer
operations, allowing for faster computation, and since they requre fewer
weights, thereby allowing the model to have a smaller memory footprint. The
GRU consists of three “gates”: the update gate $z[n]$, reset gate $r[n]$, and
the new gate $c[n]$. These gates are used to compute the cell’s current output
$h[n]$ from its current input $x[n]$ and previous output $h[n-1]$ as follows:

$z[n]=\sigma(W_{z}x[n]+U_{z}h[n-1]+b_{z})$ (6)
$r[n]=\sigma(W_{r}x[n]+U_{r}h[n-1]+b_{r})$ (7) $c[n]=\tanh(W_{c}x[n]+r[n]\circ
U_{c}h[n-1]+b_{c})$ (8) $h[n]=z[n]\circ h[n-1]+(1-z[n])\circ c[n]$ (9)

Where $W_{z},W_{r},W_{c}$ are the kernel weights for each gate,
$U_{z},U_{r},U_{c}$ are the recurrent weights for each gate, and
$b_{z},b_{r},b_{c}$ are the biases for each gate. Note that as the inputs and
outputs to the GRU layer may be vectors, all products in the above equations
are assumed to be standard matrix-vector products, except those Hadamard
products denoted $\circ$. $\sigma(x)$ refers to the sigmoid function
$\sigma(x)=\frac{1}{1+e^{-x}}$.

####  3.1.2 Fully Connected Layer

A fully connected layer computes an output vector $y[n]$ from input vector
$x[n]$ as follows:

$y[n]=\alpha(Wx[n]+b)$ (10)

Where $W$ is the kernel weights, $b$ is the layer bias, and $\alpha(x)$ is the
layer activation. In our model, we use no activation, i.e., $\alpha(x)=x$.

###  3.2 Training Data

Our dataset consists of $\sim 4$ minutes of electric guitar recordings, from a
variety of electric guitars including a Fender Stratocaster and a Gibson Les
Paul. The guitars are recorded “direct” meaning that the recorded signal is
equivalent to the signal received by the pedal coming directly from the
guitar. Recordings were made using a Focusrite Scarlett audio interface at
44.1 kHz. Note that this sample rate is very low compared to that used for
other neural network models of nonlinear audio effects (e.g. [10, 12]). This
sample rate was chosen because the embedded hardware on which the final model
was implemented processes audio at this sample rate. The recordings were then
separated into segments of 0.5 seconds each, resulting in a total of 425
segments.

Since the original Klon Centaur pedal is quite expensive ($>1500$ USD), we
used a SPICE simulation of the Centaur circuit in order to obtain a “ground
truth” reference dataset. The reference dataset measures the output voltage of
the summing amplifier from the circuit at five different values for the “Gain”
potentiometer.

Figure 7:  Comparison of predicted output of the model against reference
output shown in the time domain (above) and frequency domain (below). The
frequency domain plot uses frequency band smoothing using $1/24$ octave bands
for improved clarity. Figure 8:  Training accuracy for the RNN, shown over
epochs.

###  3.3 Training

We trained our models on 400 of the 425 audio samples, saving 25 samples for
validation. Training was performed using the Adam optimizer [20], with an
initial learning rate of $2\text{\times}{10}^{-3}$. Each model was trained for
500 epochs; each training session ran for $\sim 8$ hours. Similar to [12], we
use an error-to-signal ratio (ESR) as the loss function for our models. For a
signal of lenght $N$, ESR is defined as:

$\mathcal{E}_{ESR}=\frac{\sum_{n=0}^{N-1}|y[n]-\hat{y}[n]|^{2}}{\sum_{n=0}^{N-1}|y[n]|^{2}}$
(11)

where $y[n]$ is the reference output, and $\hat{y}[n]$ is the predicted output
of the network.

####  3.3.1 Training Results

For each model, the trained network achieved a validation ESR of less than
$2\%$. Training and validation accuracies are shown in table 1. The training
accuracy over epochs is shown in fig. 8. Results comparing the output of the
network to the reference output are shown in fig. 7. Note that the high
frequency response of the RNN output is slightly damped compared to the
reference.

Gain Parameter | Training ESR | Validation ESR
---|---|---
0.0 | 0.50 | 0.70
0.25 | 0.51 | 0.57
0.5 | 0.57 | 0.50
0.75 | 0.70 | 0.67
1.0 | 1.63 | 1.72
Table 1:  Training and validation accuracies given in error-to-signal ratio
percentages for each trained RNN model.

###  3.4 Advantages and Limitations

The recurrent neural network is a flexible and powerful black-box modelling
tool for stateful nonlinear systems. The main limitation of the RNN model is
its computational complexity for large models, mostly due to the fact that the
$\tanh$ and sigmoid functions required by the recurrent layer can be costly to
compute. Further, it can be difficult to include control parameters in the
model, a persistent challenge with black-box approaches. Finally, the
recurrent neural network cannot be used at arbirtary sample rates, and must be
trained at the same sample rate that is used for processing.

##  4 Implementation

In order to compare the virtual analog methods described above, we construct
two emulations of the Klon Centaur circuit: one emulation using traditional
circuit modelling methods (non-ML implementation), and a second using a
recurrent neural network (ML implementation). The Centaur circuit can be
broken down into four separable parts (see fig. 10):

  1. 1.

Input Buffer

  2. 2.

Gain Stage

  3. 3.

Tone Control

  4. 4.

Output Buffer

Due to their relative simplicity and linearity, in both emulations the input
buffer, output buffer, and tone control circuits were modelled using nodal
analysis. The “Gain Stage” circuit can be further broken down into six
(mostly) separable parts (see fig. 11):

  1. 1.

Feed-Forward Network 1 (FF-1)

  2. 2.

Feed-Forward Network 2 (FF-2)

  3. 3.

Pre-Amp Stage

  4. 4.

Amplifier Stage

  5. 5.

Clipping Stage

  6. 6.

Summing Amplifier

In the ML implementation, we treat the Gain Stage as a black box with a single
user-facing control (the “Gain” control). The RNN model is designed to
completely replace the Gain Stage in the circuit model. In the non-ML
implementation, we use nodal analysis to model the amplifier stage, and
summing amplifier circuits. For FF-2 and the clipping stage, we use a wave
digital filter. Since FF-1 and the pre-amp circuit share a capacitor, we
construct a joint WDF model of these two circuits, using the voltage output
from the pre-amp circuit as the input to the amplifier stage, and the current
output from FF-1 (summed with the current outputs of FF-2 and the clipping
stage) as the input to the summing amplifier.

Figure 9:  Audio plugin implementation of the Klon Centaur circuit model. Note
controls for “Gain”, “Tone”, and “Level” analogous to the original circuit, as
well as the “Traditional/Neural” parameter to control whether the emulation
uses the traditional circuit model, or the RNN model. Figure 10:  Full circuit
schematic for the Klon Centaur guitar pedal with different circuit sections
outlined. Adapted from [2]. Figure 11:  Circuit schematic for the gain stage
from the Klon Centaur pedal, with the two Feed-Forward networks highlighted.
Adapted from [2].

###  4.1 Audio Plugin

Digital audio effects are often implemented as audio plugins that can be used
by mixing engineers, producers, and musicians in a consumer digital audio
workstation (DAW) software. Common plugin formats include the Avid Audio
Extension (AAX), Steinberg’s Virtual Studio Technology (VST), and Apple’s
Audio Unit (AU) for desktop use, as well as Apple’s Audio Unit v3 (AUv3) for
mobile use. The JUCE C++ framework111https://github.com/juce-framework/JUCE is
commonly used to create cross-platform, cross-format plugins.

As a demonstration of the two circuit emulations, we construct an audio plugin
containing both models, allowing the user to switch between the two models for
comparison. The plugin is implemented using JUCE/C++, along with a real-time
Wave Digital Filter
library222https://github.com/jatinchowdhury18/WaveDigitalFilters for the WDF
models. For computing the output of the RNN models, we have implemented a
custom inferencing engine in C++, with two modes, one using the Eigen linear
algebra library [21], the second using only the C++ standard library. In the
future, we plan to add a third mode that uses the Tensorflow Lite
library.333https://www.tensorflow.org/lite/

###  4.2 Embedded Implementation

Digital audio effects are sometimes implemented on embedded devices for use in
stage performances, often in the form of a guitar pedal, or synthesizer
module. Deploying an audio effect on an embedded device can be difficult, due
to the constraints in processing power and memory availability. Further, in
order to achieve a more expressive performance, musicians often prefer effects
that add minimal latency to the signal, meaning that the embedded
implementation must be able to run with a very small buffer size.

We chose the Teensy 4.0 microcontroller as our embedded platform, since it
contains a reasonably powerful floating point processor at a relatively low
price point. The Teensy can be purchased along with an Audio Shield, which
provides 16-bit stereo audio input/output at 44.1 kHz sampling rate. The
Teensy has gained popularity in the audio community due to the Teensy Audio
Library444https://www.pjrc.com/teensy/td_libs_Audio.html that contains useful
audio DSP functionality, as well as the Faust programming language which
allows audio effects and synthesizers made in Faust to be exported for use on
the Teensy [22]. The Teensy 4.0 with the audio shield can be purchased for 35
USD.

The Teensy implementation is writteen in C++ using the Teensy Audio Library,
along with the same WDF library as the audio plugin, and the standard library
mode of the same RNN inferencing engine. The emulation can be compiled to use
either the ML or non-ML implementation. Variables in the code can be connected
to potentiometers or push-buttons to control model parameters in real-time.

Figure 12:  Teensy microcontroller implementation.

##  5 Results

The results of the real-time implementations described above can best be seen
through audio performance examples. To that end, we provide video examples of
both implementations being used in real-time on a guitar input being performed
live. These examples can be seen on
YouTube.555https://www.youtube.com/playlist?list=PLrcXtWXbPsj11cNBamVyMmDcWY1SXZHvz
From subjective listening, the ML and non-ML implementations sound very
similar, although the ML implementation has slightly damped high frequencies,
as predicted by the results of model training (see fig. 7). The high frequency
damping is slightly more noticeable when the audio input is something other
than a guitar, e.g. drums. These issues could likely be alleviated by training
on a more diverse set of audio, and possibly by adjusting the loss function to
weight high frequencies more.

###  5.1 Performance

We also evaluate the computational performance of the emulations. For real-
time performance it is important to have fast computational performance in
order to reduce audio latency. In table 2, we show the compute time per second
of audio processed of the various models at different input block sizes. Note
that at all block sizes, the ML implementation outperforms the non-ML
implementation. Performance evaluation was completed using a 2017 Dell
Precision laptop with a 2.9 GHz Intel Core i7 processor.

Block Size | NonML Time | ML Time
---|---|---
8 | 0.0723437 | 0.0528792
16 | 0.0703079 | 0.0510437
32 | 0.0652856 | 0.0511147
64 | 0.0662835 | 0.0502434
128 | 0.0666593 | 0.0495194
256 | 0.0696844 | 0.0480298
512 | 0.0669037 | 0.0477946
1024 | 0.060816 | 0.0488841
2048 | 0.0695175 | 0.0488309
4096 | 0.0623839 | 0.0472191
Table 2:  Benchmark results comparing processing speed of the audio plugin
implementation using ML processing vs. non-ML processing. Speed is measured in
compute time per second of audio processed.

###  5.2 Recommendations

From the process of implementing the circuit emulations described above, we
provide the following recommendations for circuit modellers:

  * •

For simple, linear circuits, nodal analysis is the easiest and most performant
circuit modelling method.

  * •

When modularity is important, prefer Wave Digital Filters. This modularity can
refer to the circuit topology, the components in the circuit, or the way in
which the components are discretized.

  * •

For complex nonlinear systems, particularly systems with multiple nonlinear
elements, or stateful nonlinear topologies, consider using recurrent neural
networks.

  * •

Small RNNs can outperform more complex circuit modelling methods while still
maintaining model accuracy.

  * •

While handling control parameters with RNNs can be difficult, this can be
acceptably solved by training multiple models for different values of the
control parameter and fading between them in real time.

##  6 Conclusion

We have constructed two emulations of the Klon Centaur guitar pedal circuit,
using circuit modelling techniques including nodal analysis, wave digital
filters, and recurrent neural networks. We described and compared the
advantages and limitations of each method, and showed how they can be used
together to achieve good results. We implemented the circuit emulations in the
form of an audio plugin and guitar-pedal style effect embedded on a Teensy
microcontroller. Finally, we provided recommendations for utilising different
circuit modelling methods for different types of circuits, and for different
platforms. The code for both implementations, as well as the model training,
is open source and can be found on
GitHub.666https://github.com/jatinchowdhury18/KlonCentaur

In future works, we would like to extend the RNN framework to be able to
implement larger networks in real-time. Specifically, the Differentiable
Digital Signal Processingg (DDSP) library from Google’s Magenta project
implements complex audio effects including timbral transfer, dereverberation,
and more, using an auto-encoder that contains two 512-unit GRUs, along with
several other complex operations [23]. Being able to implement the DDSP auto-
encoder for use on real-time signals would be a powerful tool for musicians
and audio engineers.

##  7 Acknowledgments

The author would like to thank Pete Warden and the EE292 class at Stanford
University for inspiring this project, as well as Julius Smith, Kurt Werner,
and Jingjie Zhang for assistance with Wave Digital Filter modelling. Thanks as
well to the Center for Computer Research in Music and Acoustics (CCRMA) for
providing computing resources.

## References

  * [1] West Warren,  “Builder Profile: Klon’s Bill Finnegan,”  Premier Guitar.
  * [2] ElectroSmash,  “Klon centaur analysis,” https://www.electrosmash.com/klon-centaur-analysis.
  * [3] François Germain,  Non-oversampled physical modeling for virtual analog simulations,  Ph.D. thesis, Stanford University, June 2019.
  * [4] D.T. Yeh,  Digital Implementation of Musical Distortion Circuits by Analysis and Simulation,  Ph.D. thesis, Stanford Univeristy, June 2009.
  * [5] Antoine Falaize and Thomas Hélie,  “Passive guaranteed simulation of analog audio circuits: A Port-Hamiltonian approach,”  Applied Sciences, vol. 6(10), pp. 273, Sept. 2016.
  * [6] A. Fettweis,  “Wave digital filters: Theory and practice,”  Proceedings of the IEEE, vol. 74, no. 2, pp. 270–327, Feb. 1986\.
  * [7] Kurt James Werner,  Virtual Analog Modeling of Audio Circuitry Using Wave Digital Filters,  Ph.D. thesis, Stanford Univeristy, June 2016.
  * [8] Martin Holters and Udo Zölzer,  “A generalized method for the derivation of non-linear state-space models from circuit schematics,”  in Proc. of the 23rd European Signal Processing Conference (EUSIPCO), Sept. 2015, pp. 1078–1082.
  * [9] Julius O. Smith,  Spectral Audio Signal Processing,  http:http://ccrma.stanford.edu/ jos/sasp///ccrma.stanford.edu/~jos/sasp/, accessed 2020-5-22,  online book, 2011 edition.
  * [10] Eero-Pekka Damskägg, Lauri Juvela, and Vesa Välimäki,  “Real-time modeling of audio distortion circuits with deep learning,”  in Proc. of the 16th Sound and Music Computing Conference, May 2019\.
  * [11] Julian D. Parker, Fabián Esqueda, and André Bergner,  “Modelling of nonlinear state-space systems using a deep neural network,”  in Proc. of the 22nd Int. Conference on Digital Audio Effects (DAFx-19), Sept. 2019.
  * [12] Alec Wright, Eero-Pekka Damskägg, and Vesa Välimäki,  “Real-time black-box modelling with recurrent neural networks,”  in Proc. of the 22nd Int. Conference on Digital Audio Effects (DAFx-19), Sept. 2019.
  * [13] Edward W. Maby,  Solid State Electronic Circuits, vol. 4,  1 edition, 2014.
  * [14] Julius O. Smith,  Physical Audio Signal Processing,  http:http://ccrma.stanford.edu/ jos/pasp///ccrma.stanford.edu/~jos/pasp/, accessed 2020-5-22,  online book, 2010 edition.
  * [15] Chung-Wen Ho, A. Ruehli, and P. Brennan,  “The modified nodal approach to network analysis,”  IEEE Transactions on Circuits and Systems, vol. 22, no. 6, pp. 504–509, 1975.
  * [16] Marco A. Martínez Ramirez and Joshua D. Reiss,  “Modeling of nonlinear audio effects with end-to-end deep neural networks,”  arXiv e-prints, p. arXiv:1810.06603, Oct. 2018.
  * [17] Jatin Chowdhury,  “Complex nonlinearities for audio signal processing,” https://ccrma.stanford.edu/~jatin/papers/Complex_NLs.pdf, Feb. 2020.
  * [18] François Chollet et al.,  “Keras,” https://github.com/fchollet/keras, 2015.
  * [19] Kyunghyun Cho, Bart van Merrienboer, Çaglar Gülçehre, Fethi Bougares, Holger Schwenk, and Yoshua Bengio,  “Learning phrase representations using RNN encoder-decoder for statistical machine translation,”  CoRR, vol. abs/1406.1078, 2014.
  * [20] Diederik P. Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  CoRR, vol. abs/1412.6980, 2015.
  * [21] Gaël Guennebaud, Benoît Jacob, et al.,  “Eigen v3,” http://eigen.tuxfamily.org, 2010.
  * [22] Romain Michon, Yann Orlarey, Stéphane Letz, and Dominique Fober,  “Real time audio digital signal processing with Faust and the Teensy,”  in Proc. of the 16th Sound and Music Computing Conference, May 2019\.
  * [23] Jesse Engel, Lamtharn (Hanoi) Hantrakul, Chenjie Gu, and Adam Roberts,  “DDSP: Differentiable digital signal processing,”  in International Conference on Learning Representations, 2020.