# Image Trinarization Using a Partial Differential Equations: A Novel Approach
to Automatic Sperm Image Analysis

B. A. Jacobs<EMAIL_ADDRESS>University of Johannesburg, Department of
Mathematics and Applied Mathematics, South Africa  Vitruvian Medical
Diagnostics, YTC Building 02-01, 33 Maude Road, 208344, Singapore

###### Abstract

Partial differential equations have recently garnered substantial attention as
an image processing framework due to their extensibility, the ability to
rigorously engineer and analyse the governing dynamics as well as the ease of
implementation using numerical methods. This paper explores a novel approach
to image trinarization with a concrete real-world application of classifying
regions of sperm images used in the automatic analysis of sperm morphology.
The proposed methodology engineers a diffusion equation with non-linear source
term, exhibiting three steady-states. The model is implemented as an image
processor using a standard finite difference method to illustrate the efficacy
of the proposed approach. The performance of the proposed approach is
benchmarked against standard image clustering/segmentation methods and shown
to be highly effective.

###### keywords:

Image processing , Nonlinear PDE of parabolic type , Stability and convergence
of numerical methods , Biomedical imaging and signal processing

###### MSC:

68U10 , 35K55 , 65M12 , 92C55

††journal: Applied Mathematics and Computation

##  1 Introduction

Partial differential equations (PDEs) are a well established mathematical
modelling tool and have been proven to show enormous value in a multitude of
disciplines. More recently, PDEs have become a useful tool in image
processing. Applications in image denoising, restoration, inpainting, remote
sensing applications, document image binarization, [1, 2, 3, 4, 5, 6, 7, 8, 9,
10], are a small subset of applications that illustrate the value of using a
PDE framework to attain novel image processing algorithms. Additionally,
research on processing medical images, such as ultrasounds and MRIs [11, 12,
13].

This work builds on previous works by the author [5, 6, 8] where PDE models
were used for document image binarization, implementing global thresholding,
locally adaptive thresholding and a nonlinear system of equations which self-
governed the model parameters resulting in an unsupervised methodology. The
approach was also extended to time fractional PDEs presented in [7]. Many
other researchers have also contributed to the field using a variety of PDE
based approaches. Feng [14, 15] presents variational level set approaches to
the problem. An indirect diffusion approach to image segmentation is presented
in [16], while Guo et al. [17] show how indirect diffusion with shock filters
can be used for text binarization. Zhang et al. [18] present a nonlinear
reaction diffusion model which exhibits the Perona-Malik [1] type diffusion
(an edge preserving diffusion) coupled with a binarizing reaction term. An
interger/fractional approach is combined with an edge detection to binarize
heavily degraded document images [19, 20]. Du and He [21] propose a nonlinear
diffusion model with selective source term which selectively restores the
background or preserves the fidelity of the text.

The task of image or data trinarization has received significantly less
attention. The task aims to segment image data into three distinct classes
corresponding to some desirable segmentation. This problem can be viewed as a
base case of image segmentation, or component clustering and hence is amenable
to standard approaches. Specifically K-Means [22], K-Medoids [23], Spanning
Trees [24], Agglomerative clustering [25]. This work benchmarks the proposed
PDE approach against these standard approaches and illustrates that the fit-
for-purpose PDE approach vastly outperforms the generalist approaches.

Kido and Higo [26] present the application of water retention tests, the
results of which are amenable to data trinarization. The various phases of the
test matter are trinarized using their approach presented in [27]. This
technique is based on the frequency selection based on the image histogram
data, a natural extension of Otsu’s [28] technique for image binarization.
Another approach taken by Higo et al. [29] for this application is the
implementation of a region growing method. Müssel et al. [30] present a
multiscale approach for the binarization and trinarization of gene expression
profiles. This application is motivated by the requirement of classifying data
as ‘high’, ‘low’ or potentially ‘intermediate’, which is particularly
pertinent for biological data.

Manual assessment of sperm morphology is a subjective process, which is prone
to inter-expert variance and assessments are difficult to reproduce [31]. To
address this problem computer aided semen analysis (CASA) has been an active
research field for decades and has aimed to mitigate the subjectivity of the
manual process. A natural approach to tackling the classification task of
sperm morphology from images is by way of neural networks. Miahi et al. [31]
present a novel Genetic Neural Architecture for the automatic assessment of
sperm images. Javadi and Mirroshandel [32] provide a deep learning approach,
applying a convolutional deep neural network to the Sperm Morphology Analysis
dataset (MHSMA). Savkay et al. [33] present a convolutional neural network
approach to CASA, presenting various morphological parameter predictions. An
improved convolutional deep neural network approach is presented in [34]. A
comprehensive review of artificial intelligence in fertility is presented by
Riegler et al. [35] which covers the broader subject matter of assisted
reproductive technology, but collates key contributions to CASA using
artificial intelligence. Dai et al. [36] present a review of recent advances
in sperm analysis, illustrating that the field is an active one, garnering the
attention of many researchers.

The recent advances in CASA is largely focused on machine learning approaches.
These techniques are well established as effective tools, but do suffer from
some drawbacks typically associated with ‘black-box’ methods. Additionally,
the estimation of the acrosomal region is a key factor in the assessment of
human sperm and is inherently subjective. According to the WHO Laboratory
Manual for the Examination and Processing of Human Semen [37] guidelines for
sperm analysis this region should comprise 40% to 70% of the total sperm head.
Figures 1(a) and 1(b) illustrate an example of the sperm head and acrosomal
region as well as a schematic of the sperm anatomy.

(a) Example of Sperm Microscopy Image.

(b) Diagram illustrating sperm anatomy.

Building on the PDE based framework for image processing, the task of image
trinarization is approached with the same lens. This implementation is
motivated by the application of automatic sperm analysis from mircosopic
images. Establishing this PDE based framework as an image processing tool
gives credence to the method for tackling new data driven problems. In the
case of morphological sperm analysis, the dynamic nature of the PDE clustering
method allows for each region to dynamically transit to its final
classification. As has been shown before, this results in a method which is
robust to noise, extensible and accurate.

The following section introduces the proposed model and discusses the
methodologies used to generate the results for automated sperm segmentation
which are then discussed in Section 3. The paper concludes in Section 4 with
insights given by the results and also how the present work opens new
possibilities within the field of image processing with PDEs.

##  2 Methods

###  2.1 Model

The model, given below, is engineered to exhibit three steady states,
corresponding to background, sperm head and acrosomal region. Explicitly, the
model is defined by

$\displaystyle u_{t}$ $\displaystyle=c_{D}\nabla u+c_{S}u(u-1)(u-a)(u-b),$ (1)
$\displaystyle u(x,y,0)$ $\displaystyle=\text{Image}(x,y).$

The given model is constructed in two spatial dimensions, requiring that the
input image be grayscale. Although, the extension to RGB images is
straightforward under the PDE framework. The model is also subjected to zero-
flux boundary conditions along all boundaries of the domain. This choice of
boundary condition ensures smoothness at the boundaries and is motivated by
previous works [5, 6, 8]. Explicitly, the boundary conditions are given by,

$\displaystyle u_{x}(0,y,t)=0$ $\displaystyle u_{x}(1,y,t)=0$ $\displaystyle
u_{y}(x,0,t)=0$ $\displaystyle u_{y}(x,1,t)=0$

In engineering the dynamics of the present model, we require the model to
exhibit three stable steady states. This is achieved by the careful
construction of the source term. Figure 2 illustrate the phase space dynamics
of the source term. This structure is an extension of the Fitzhugh-Nagumo
model used in previous works [5, 6, 8]. The source term exhibits five steady
states; $u=0,u=a,u=b,u=c$ and $u=1$. Assuming that $0\leq a<b<c\leq 1$, the
stable states exist at $u=0,u=1$ and $u=b$ as illustrated in Figure 2.

Figure 2: Phase space plot of the nonlinear source term, for $a=0.25$, $b=0.5$
and $c=0.75$.

Moreover, due to the nature of the input images, specifically high
magnification with uniform illumination, there is likely minimal noise
captured in the image which is associated with the camera sensor. The
coefficients $c_{D}$ and $c_{S}$ are tunable to balance the bias of smoothing
(denoising) versus trinarization.

###  2.2 Numerical Method

The standard finite-difference implementation offers a conditionally stable,
explicit and parallelizable [6] method. Although the method is rudimentary its
ease of implementation offers an attractive option for robust experimentation.
The model is discretised as follows:

$\displaystyle u(x_{i},y_{j},t_{n})$ $\displaystyle=u(i\Delta x,j\Delta
y,n\Delta n)=u_{i,j}^{n},$ (2) $\displaystyle u_{t}(x,y,t)$
$\displaystyle=\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+\mathcal{O}(\Delta
t),$ (3) $\displaystyle u_{x}(x,y,t)$
$\displaystyle=\frac{u_{i+1,j}^{n}-u_{i-1,j}^{n}}{2\Delta
x}+\mathcal{O}(\Delta x^{2}),$ (4) $\displaystyle u_{y}(x,y,t)$
$\displaystyle=\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2\Delta
y}+\mathcal{O}(\Delta y^{2}),$ (5) $\displaystyle u_{xx}(x,y,t)$
$\displaystyle=\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta
x^{2}}+\mathcal{O}(\Delta x^{2}),$ (6) $\displaystyle u_{yy}(x,y,t)$
$\displaystyle=\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta
y^{2}}+\mathcal{O}(\Delta y^{2}).$ (7)

The computational domain is prescribed by the input image, but is artificially
normalised to $(x,y)\in[0,1]\times[0,1]$. The image dimensions, $N$ and $M$,
thus dictate $\Delta x=1/(N-1)$ and $\Delta y=1/(M-1)$. Discretising the model
in (1) yields the explicit update scheme,

$\displaystyle u_{i,j}^{n+1}=u_{i,j}^{n}$ $\displaystyle+\frac{\Delta
tc_{D}}{\Delta x^{2}}\left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right)$
$\displaystyle+\frac{\Delta tc_{D}}{\Delta
y^{2}}\left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)$
$\displaystyle-\Delta
tc_{S}u_{i,j}^{n}(u_{i,j}^{n}-1)(u_{i,j}^{n}-a)(u_{i,j}^{n}-b)(u_{i,j}^{n}-c).$
(8)

The resulting scheme is consistent with the model equation (1) as $\Delta
x,\Delta y,\Delta t\rightarrow 0$.

The boundary conditions are discretised using the central difference formulae
(4) and (5) to ensure second order accuracy at the boundaries. Explicitly we
have

$\displaystyle u_{x}(0,y,t)$ $\displaystyle=0\Rightarrow
u_{-1,j}^{n}=u_{1,j}^{n},$ (9) $\displaystyle u_{x}(1,y,t)$
$\displaystyle=0\Rightarrow u_{N+1,j}^{n}=u_{N-1,j}^{n},$ (10) $\displaystyle
u_{y}(x,0,t)$ $\displaystyle=0\Rightarrow u_{j,-1}^{n}=u_{j,1}^{n},$ (11)
$\displaystyle u_{y}(x,1,t)$ $\displaystyle=0\Rightarrow
u_{M+1,j}^{n}=u_{M-1,j}^{n}.$ (12)

For simplicity, the scheme can be written in matrix form

$\displaystyle{\bf u}^{n+1}$ $\displaystyle={\bf u}^{n+1}+\frac{cD\Delta
t}{\Delta x^{2}}{\bf A}.{\bf u}^{n}+\frac{c_{D}\Delta t}{\Delta y^{2}}{\bf
u}^{n}.{\bf A}^{T}$ $\displaystyle-c_{S}\Delta t{\bf u}^{n}\circ({\bf
u}^{n}-1)\circ({\bf u}^{n}-a)\circ({\bf u}^{n}-b)\circ({\bf u}^{n}-c),$ (13)

where $\circ$ denotes the Hadamard product and

${\bf
u}=\left(\begin{array}[]{ccccc}u_{0,0}^{n}&u_{0,1}^{n}&\ldots&u_{0,M-1}^{n}&u_{0,M}^{n}\\\
u_{1,0}^{n}&u_{1,1}^{n}&&u_{1,M-1}^{n}&u_{1,M}^{n}\\\
\vdots&&\ddots&&\vdots\\\
u_{N-1,0}^{n}&u_{N-1,1}^{n}&&u_{N-1,M-1}^{n}&u_{N-1,M}^{n}\\\
u_{N,0}^{n}&u_{N,1}^{n}&\ldots&u_{N,M-1}^{n}&u_{N,M}^{n}\\\
\end{array}\right),$ (14) ${\bf
A}=\left(\begin{array}[]{ccccccc}-2&2&0&\ldots&0&0&0\\\ 1&-2&1&&0&0&0\\\
0&1&-2&&0&0&0\\\ \vdots&&&\ddots&&&\vdots\\\ 0&0&0&&-2&1&0\\\ 0&0&0&&1&-2&1\\\
0&0&0&\ldots&0&2&-2\\\ \end{array}\right).$ (15)

###  2.3 Stability Analysis

To accommodate the nonlinearity in the stability analysis, the source term is
replaced by the worst case contribution. This leads to stability conditions
which are stricter than necessary in most cases. Nevertheless, the stability
requirement does not overly impose on the functioning the method.
Consider the source term in (1),

$f(u)=-u(u-1)(u-a)(u-b)(u-c).$ (16)

This function is maximal (in absolute value) under constraints $0\leq u\leq
1$, $0\leq a\leq 1$ and $0\leq b\leq 1$ with

$\max_{u,a,b,c}|f(u)|=\frac{256}{3125}.$ (17)

The von Neumann stability analysis is conducted by introducing the ansatz

$u_{i,j}^{n}=\zeta^{n}e^{Ii\Delta xj\Delta y\omega},$ (18)

where $\omega$ indicates an arbitrary Fourier wave number and $I=\sqrt{-1}$.
Deriving the conditions on $\Delta x$, $\Delta y$ and $\Delta t$ that ensure
the amplification factor, $\zeta$, is smaller than one in absolute value,
$|\zeta|<1$, is sufficient to guarantee the numerical method is stable.
Substituting the ansatz into equation (8) we have

$\displaystyle\zeta=1$ $\displaystyle-\frac{2\Delta tc_{D}}{\Delta
x^{2}}+\frac{2\Delta tc_{D}}{\Delta x^{2}}\frac{e^{I\Delta xj\Delta
y\omega}+e^{-I\Delta xj\Delta y\omega}}{2}$ $\displaystyle-\frac{2\Delta
tc_{D}}{\Delta y^{2}}+\frac{2\Delta tc_{D}}{\Delta y^{2}}\frac{e^{Ii\Delta
x\Delta y\omega}+e^{-Ii\Delta x\Delta y\omega}}{2}+\Delta
tc_{S}f(u_{i,j}^{n}),$ $\displaystyle=1$ $\displaystyle+\frac{2\Delta
tc_{D}}{\Delta x^{2}}\left(\cos\left(2\frac{j\Delta x\Delta
y\omega}{2}\right)-1\right)$ $\displaystyle+\frac{2\Delta tc_{D}}{\Delta
y^{2}}\left(\cos\left(2\frac{i\Delta x\Delta
y\omega}{2}\right)-1\right)+\Delta tc_{S}f(u_{i,j}^{n}),$ $\displaystyle=1$
$\displaystyle-\frac{2\Delta tc_{D}}{\Delta x^{2}}\sin^{2}\left(j\Delta
x\Delta y\omega\right)$ $\displaystyle-\frac{2\Delta tc_{D}}{\Delta
y^{2}}\sin^{2}\left(i\Delta x\Delta y\omega\right)+\Delta
tc_{S}f(u_{i,j}^{n}).$ (19)

The stability condition $|\zeta|<1$ leads to

$\displaystyle\Delta tc_{S}|f(u_{i,j}^{n})|$ $\displaystyle\leq\Delta
tc_{S}\frac{256}{3125},$ $\displaystyle<\frac{2\Delta tc_{D}}{\Delta
x^{2}}\sin^{2}\left(j\Delta x\Delta y\omega\right)+\frac{2\Delta
tc_{D}}{\Delta y^{2}}\sin^{2}\left(i\Delta x\Delta y\omega\right),$
$\displaystyle<2\Delta tc_{D}\left(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta
y^{2}}\right),$

and

$\displaystyle-2$ $\displaystyle<-\frac{2\Delta tc_{D}}{\Delta
x^{2}}\sin^{2}\left(j\Delta x\Delta y\omega\right)-\frac{2\Delta
tc_{D}}{\Delta y^{2}}\sin^{2}\left(i\Delta x\Delta y\omega\right)+\Delta
tc_{S}f(u_{i,j}^{n})$ $\displaystyle\leq-\frac{2\Delta tc_{D}}{\Delta
x^{2}}-\frac{2\Delta tc_{D}}{\Delta y^{2}}+\Delta tc_{S}f(u_{i,j}^{n})$
$\displaystyle 2\Delta tc_{D}\left(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta
y^{2}}-\frac{1}{\Delta tc_{D}}\right)$ $\displaystyle\leq\Delta
tc_{S}|f(u_{i,j}^{n})|\leq\Delta tc_{S}\frac{256}{3125}.$ (21)

The above criteria can be summarised as

$|\Delta tc_{S}\frac{128}{3125}|\leq c_{D}\left(\frac{\Delta t}{\Delta
x^{2}}+\frac{\Delta t}{\Delta y^{2}}\right).$ (22)

Since $\Delta x$ and $\Delta y$ are prescribed by the input image’s
dimensions, $\Delta t$ should be chosen so that $\frac{\Delta t}{\Delta
x^{2}}$ and $\frac{\Delta t}{\Delta y^{2}}$ are $\mathcal{O}(1)$. Then $c_{S}$
and $c_{D}$ can be chosen to balance the effect of trinarization and smoothing
respectively, while still obeying equation (22), thus ensuring the resulting
implementation is stable.

##  3 Medical Image Application

This section presents the performance results of the present method as well as
the benchmark figures for standard methodologies for image segmentation.
Performance metrics are discussed which enable the quantitative measurement of
performance of various methodologies. The application of segmentation of sperm
head and acrosomal regions is an important aspect of CASA and has been an
active field of research for many years. This work presents the novel approach
of performing this segmentation using a PDE based framework. The results
presented here illustrate that the PDE based approach allows for a PDE model
to be engineered to exhibit desirable dynamics required by the application,
and more importantly, provide a robust, provably stable and highly accurate
method for CASA.

The images of sperm obtained from microscopy shown in Figures 4(a), 5(a),
6(a), 7(a) and 8(a), illustrate the subjective nature of discerning exactly
where the acrosomal regions begin and end. As such, a robust algorithm which
provides repeatable results is of significant importance in CASA.

###  3.1 Algorithm Description

In order to improve the classification results a post processing step is
implemented, which is specific to the CASA application. Specifically this
optimisation removes extraneous information at the edges of the image, or
trimming of the sperm tail by constructing a mask from the image using
traditional image processing techniques. This post processing was done for all
clustering methods for a fair comparison of the segmentation task. The overall
algorithm is described in Algorithm 1.

Algorithm 1 Segmentation and Post-Processing Algorithm

1:$img\leftarrow\text{Input Image}$

2:$result\leftarrow clusteringMethod(img)$ $\triangleright$ from selected
clustering method

3:$result\leftarrow round(2*result)/2$ $\triangleright$ force all pixel values
to 0, 0.5 or 1

4:$mask\leftarrow closing(res,disk(20))$ $\triangleright$ image mask formed
with disk kernel

5:$result\leftarrow applyMask(result,mask)$ $\triangleright$ Mask resulting
image

###  3.2 Performance Metrics

The quantitative performance of the assessed clustering methodologies requires
a robust performance measure. The performance of each algorithm is broken down
into two binary classification tasks; measuring the number of pixels
classified as head and measuring the number of pixels classified as acrosome.
The classification results are shown in Figures 4, 5, 6, 7 and 8, where class
1 is shown in grey and class 2 is shown in black.

A natural approach to measuring the accuracy of a binary classification task
is to count the number of pixels classified as black which are also classified
as black in the ground truth image, and scale that count by the total number
of pixels in the image; in essence giving a percentage of “correctly
classified” pixels. These results are presented in Table 2 for each image and
each method. However, this is biased in the case where the number of target
pixels is much less than the total number of pixels in the entire image. For
example, if the binary ground truth image contained only one black pixel, then
a pure white image would be considered almost perfectly accurate, but in
reality the importance of capturing that one black pixel is the core task. The
overestimation of a method’s performance is evident in the high accuracies
obtained by all methods in Table 2

As such, many researchers adopt the F-measure, or $F_{1}$ score [38] as a more
meaningful measure of accuracy in binary classification tasks. This measure is
the harmonic mean of the precision and recall of the classification result,
essentially compensating for the bias in class representation. The $F_{1}$
score is calculated as

$F_{1}=2\frac{\text{precision}\times\text{recall}}{\text{precision}+\text{recall}}=\frac{\text{tp}}{\text{tp}+\frac{1}{2}\left(\text{fp}+\text{fn}\right)}.$
(23)

The $F_{1}$ scores obtained by all methods across all images are presented in
Table 3 and it can be seen that these scores correlate more meaningfully with
the qualitative results shown in Figures 4, 5, 6, 7 and 8.

###  3.3 Results

The proposed method relies on the tuning of parameters $a$, $b$ and $c$,
however in the experiments presented below, the algorithm proved to be robust
for all input images, with fixed parameters of $a=0.5$, $b=0.65$ and $c=0.7$.
As depicted in Figure 2, the value of $b$ determines the stable value to which
pixels will be driven. Due to this the value of $b$ can be fixed and $a$ and
$c$ varied to change the relative strength of the source term in driving pixel
values toward $0$, $b$ or $1$.

Typical values were selected for the remaining model parameters, namely
$c_{D}=0.01$, $c_{S}=1/\Delta t$ and $\Delta t=\Delta x\Delta y/4$ where
$\Delta x=1/N$ and $\Delta y=1/M$ with $(N,M)$ being the dimensions of the
input image. The choice of these parameters is informed by previous studies,
indicating that the relatively small value of $c_{D}$ is sufficient to
introduce classification dynamics when the image is largely free of noise. The
diffusion term encourages smoothing of noisy data and boundary regions.

The $F_{1}$ scores as a function of $a$ and $c$ are plotted in Figure 3. These
figures illustrate that near optimal results can be obtained from the present
methodology using the aforementioned fixed parameter values.

(a) Contour error plot for varying $a$ and $c$ with fixed $b=0.65$.

(b) Contour error plot for varying $a$ and $c$ with fixed $b=0.65$.

(c) Contour error plot for varying $a$ and $c$ with fixed $b=0.65$.

(d) Contour error plot for varying $a$ and $c$ with fixed $b=0.65$.

(e) Contour error plot for varying $a$ and $c$ with fixed $b=0.65$.

Figure 3: Contour plots of $F_{1}$ scores for all example images for various
values of $a$ and $c$.

Figures 4, 5, 6, 7 and 8 illustrate the classification outputs obtained using
the present method as well as four standard unsupervised image segmentation
algorithms; namely Agglomerative clustering, K-Means, K-Medoids and Spanning
Trees. The average performance results for all examples of the present method
and aforementioned methods are presented in Table 1.

(a) Input Image.

(b) Ground Truth.

(c) PDE Result.

(d) Agglomerative Clustering.

(e) K Means Clustering.

(f) K Medoids Clustering.

(g) Spanning Trees.

Figure 4: Input image, ground truth and results from a variety algorithms for
Example 1.

(a) Input Image.

(b) Ground Truth.

(c) PDE Result.

(d) Agglomerative Clustering.

(e) K Means Clustering.

(f) K Medoids Clustering.

(g) Spanning Trees.

Figure 5: Input image, ground truth and results from a variety algorithms for
Example 2.

(a) Input Image.

(b) Ground Truth.

(c) PDE Result.

(d) Agglomerative Clustering.

(e) K Means Clustering.

(f) K Medoids Clustering.

(g) Spanning Trees.

Figure 6: Input image, ground truth and results from a variety algorithms for
Example 3.

(a) Input Image.

(b) Ground Truth.

(c) PDE Result.

(d) Agglomerative Clustering.

(e) K Means Clustering.

(f) K Medoids Clustering.

(g) Spanning Trees.

Figure 7: Input image, ground truth and results from a variety algorithms for
Example 4.

(a) Input Image.

(b) Ground Truth.

(c) PDE Result.

(d) Agglomerative Clustering.

(e) K Means Clustering.

(f) K Medoids Clustering.

(g) Spanning Trees.

Figure 8: Input image, ground truth and results from a variety algorithms for
Example 5.

Table 1: Table reporting the average $F_{1}$ scores and accuracy over all images for all methods. Method |  $F_{1}$ Class 1 |  $F_{1}$ Class 2 |  Average $F_{1}$ | Accuracy Class 1 | Accuracy Class 2 | Average Accuracy
---|---|---|---|---|---|---
Agglomerate | 0.214378 | 0.708983 | 0.46168 | 0.950434 | 0.946978 | 0.948706
KMeans | 0.45513 | 0.798757 | 0.626943 | 0.962562 | 0.966829 | 0.964696
KMedoids | 0.454328 | 0.79845 | 0.626389 | 0.962536 | 0.966779 | 0.964658
PDE | 0.784823 | 0.894653 | 0.839738 | 0.975313 | 0.98833 | 0.981821
SpanningTree | 0.45042 | 0.505314 | 0.289399 | 0.942096 | 0.963594 | 0.952845

Table 2: Table reporting the accuracy for all images and all methods. Example | Method | Accuracy Class 1 | Accuracy Class 2 | Average Accuracy
---|---|---|---|---
1 | Agglomerate | 0.92285 | 0.925719 | 0.924284
2 | Agglomerate | 0.956875 | 0.961444 | 0.959159
3 | Agglomerate | 0.947506 | 0.936206 | 0.941856
4 | Agglomerate | 0.93955 | 0.922138 | 0.930844
5 | Agglomerate | 0.985706 | 0.993394 | 0.98955
1 | KMeans | 0.976012 | 0.984981 | 0.980497
2 | KMeans | 0.954438 | 0.964781 | 0.959609
3 | KMeans | 0.96795 | 0.973869 | 0.970909
4 | KMeans | 0.935137 | 0.927269 | 0.931203
5 | KMeans | 0.979019 | 0.986556 | 0.982788
1 | KMedoids | 0.976012 | 0.984981 | 0.980497
2 | KMedoids | 0.954269 | 0.964781 | 0.959525
3 | KMedoids | 0.96795 | 0.973869 | 0.970909
4 | KMedoids | 0.935137 | 0.927269 | 0.931203
5 | KMedoids | 0.979056 | 0.986306 | 0.982681
1 | PDE | 0.971231 | 0.97955 | 0.975391
2 | PDE | 0.981831 | 0.991081 | 0.986456
3 | PDE | 0.986413 | 0.996219 | 0.991316
4 | PDE | 0.95215 | 0.982944 | 0.967547
5 | PDE | 0.984938 | 0.991856 | 0.988397
1 | SpanningTree | 0.924731 | 0.980788 | 0.952759
2 | SpanningTree | 0.935456 | 0.941144 | 0.9383
3 | SpanningTree | 0.962294 | 0.968025 | 0.965159
4 | SpanningTree | 0.901563 | 0.927925 | 0.914744
5 | SpanningTree | 0.986344 | 0.995337 | 0.990841

Table 3: Table reporting the $F_{1}$ scores for all images and all methods. Example | Method |  $F_{1}$ Score Class 1 |  $F_{1}$ Score Class 2 |  Average $F_{1}$ Score
---|---|---|---|---
1 | Agglomerate | 0.126645 | 0.611055 | 0.36885
2 | Agglomerate | 0.234184 | 0.750919 | 0.492551
3 | Agglomerate | 0.000713861 | 0.61277 | 0.306742
4 | Agglomerate | 0.0724971 | 0.648119 | 0.360308
5 | Agglomerate | 0.630473 | 0.948305 | 0.789389
1 | KMeans | 0.835336 | 0.884941 | 0.860138
2 | KMeans | 0.308087 | 0.766715 | 0.537401
3 | KMeans | 0.597741 | 0.794272 | 0.696007
4 | KMeans | 0.123035 | 0.662412 | 0.392724
5 | KMeans | 0.397848 | 0.908775 | 0.653311
1 | KMedoids | 0.835336 | 0.884941 | 0.860138
2 | KMedoids | 0.307561 | 0.766715 | 0.537138
3 | KMedoids | 0.597741 | 0.794272 | 0.696007
4 | KMedoids | 0.123035 | 0.662412 | 0.392724
5 | KMedoids | 0.394361 | 0.907244 | 0.650802
1 | PDE | 0.839039 | 0.78767 | 0.813355
2 | PDE | 0.818005 | 0.92206 | 0.870033
3 | PDE |  0.870395 | 0.962067 | 0.916231
4 | PDE | 0.684471 | 0.865969 | 0.77522
5 | PDE | 0.712204 | 0.935498 | 0.823851
1 | SpanningTree | 0 | 0.803125 | 0
2 | SpanningTree | 0.574828 | 0 | 0
3 | SpanningTree | 0.501858 | 0.759383 | 0.63062
4 | SpanningTree | 0.504093 | 0 | 0
5 | SpanningTree | 0.668688 | 0.964062 | 0.816375

##  4 Conclusion

The present work extends the PDE based framework for image processing
previously established for document image binarization [5, 6, 8]. The
application of CASA raises the issue of image trinarization, whereby an image
data is classified into three distinct classes. This classification task is
performed dynamically by the carefully constructed PDE model, which is shown
to exhibit three steady states (corresponding to the three distinct classes).
This proposed model is benchmarked against several other general purpose
clustering algorithms and shown to perform exceptionally well for a variety of
input images and is robust, even with fixed parameter values.

The extensibility of the present image processing framework opens
opportunities for novel processing algorithms, for new application domains.
Specifically, the inclusion of local adaptivity [6], or automatic control of
model parameters for unsupervised methodologies [8], or edge-preserving image
denoising [39, 1]. A multitude of applications are available to this
framework, and new architectures are necessary to be developed and explored to
present novel solutions to arising problems.

## References

  * [1] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on pattern analysis and machine intelligence 12 (7) (1990) 629–639.
  * [2] P. Liu, F. Huang, G. Li, Z. Liu, Remote-sensing image denoising using partial differential equations and auxiliary images as priors, IEEE geoscience and remote sensing letters 9 (3) (2011) 358–362.
  * [3] S. Kim, Pde-based image restoration: A hybrid model and color image denoising, IEEE Transactions on Image Processing 15 (5) (2006) 1163–1170.
  * [4] F. Catté, P.-L. Lions, J.-M. Morel, T. Coll, Image selective smoothing and edge detection by nonlinear diffusion, SIAM Journal on Numerical analysis 29 (1) (1992) 182–193.
  * [5] B. Jacobs, E. Momoniat, A novel approach to text binarization via a diffusion-based model, Applied Mathematics and Computation 225 (2013) 446–460.
  * [6] B. Jacobs, E. Momoniat, A locally adaptive, diffusion based text binarization technique, Applied Mathematics and Computation 269 (2015) 464–472.
  * [7] B. Jacobs, C. Harley, Application of nonlinear time-fractional partial differential equations to image processing via hybrid laplace transform method, Journal of Mathematics 2018 (2018).
  * [8] B. Jacobs, T. Celik, Unsupervised document image binarization using a system of nonlinear partial differential equations, Applied Mathematics and Computation 418 (2022) 126806.
  * [9] S. K. Weeratunga, C. Kamath, Comparison of pde-based nonlinear anisotropic diffusion techniques for image denoising, in: Image Processing: Algorithms and Systems II, Vol. 5014, International Society for Optics and Photonics, 2003, pp. 201–212.
  * [10] T. F. Chan, J. Shen, Image processing and analysis: variational, PDE, wavelet, and stochastic methods, SIAM, 2005.
  * [11] A. Hadri, L. Afraites, A. Laghrib, M. Nachaoui, A novel image denoising approach based on a non-convex constrained pde: application to ultrasound images, Signal, Image and Video Processing 15 (5) (2021) 1057–1064.
  * [12] X.-X. Xing, Y.-L. Zhou, J. S. Adelstein, X.-N. Zuo, Pde-based spatial smoothing: a practical demonstration of impacts on mri brain extraction, tissue segmentation and registration, Magnetic Resonance Imaging 29 (5) (2011) 731–738.
  * [13] B. Lu, C. Deng, Q. Liu, J. Li, Four order adaptive pde method for mri denoising, in: 2009 3rd International Conference on Bioinformatics and Biomedical Engineering, IEEE, 2009, pp. 1–4.
  * [14] S. Feng, A novel variational model for noise robust document image binarization, Neurocomputing 325 (2019) 288–302.
  * [15] S. Feng, Effective document image binarization via a convex variational level set model, Applied Mathematics and Computation 419 (2022) 126861.
  * [16] Y. Wang, Q. Yuan, C. He, Indirect diffusion based level set evolution for image segmentation, Applied Mathematical Modelling 69 (2019) 714–722.
  * [17] J. Guo, C. He, Y. Wang, Fourth order indirect diffusion coupled with shock filter and source for text binarization, Signal Processing 171 (2020) 107478.
  * [18] X. Zhang, C. He, J. Guo, Selective diffusion involving reaction for binarization of bleed-through document images, Applied Mathematical Modelling 81 (2020) 844–854.
  * [19] U. A. Nnolim, Dynamic selective edge-based integer/fractional-order partial differential equation for degraded document image binarization, International Journal of Image and Graphics (2021) 2250030.
  * [20] U. A. Nnolim, Improved integer/fractional order partial differential equation-based thresholding, Optik 229 (2021) 166265.
  * [21] Z. Du, C. He, Nonlinear diffusion equation with selective source for binarization of degraded document images, Applied Mathematical Modelling 99 (2021) 243–259.
  * [22] J. MacQueen, et al., Some methods for classification and analysis of multivariate observations, in: Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, Vol. 1, Oakland, CA, USA, 1967, pp. 281–297.
  * [23] L. Kaufman, P. J. Rousseeuw, Finding groups in data: an introduction to cluster analysis, John Wiley & Sons, 2009.
  * [24] T. Asano, B. Bhattacharya, M. Keil, F. Yao, Clustering algorithms based on minimum and maximum spanning trees, in: Proceedings of the fourth annual symposium on Computational Geometry, 1988, pp. 252–257.
  * [25] L. Rokach, O. Maimon, Clustering methods, in: Data mining and knowledge discovery handbook, Springer, 2005, pp. 321–352.
  * [26] R. Kido, Y. Higo, Quantification of pore volume and degree of saturation in partially saturated sands with different bulk density, Japanese Geotechnical Society Special Publication 8 (6) (2020) 216–220.
  * [27] R. Kido, Y. Higo, Evaluation of distribution of void ratio and degree of saturation in partially saturated triaxial sand specimen using micro x-ray tomography, Japanese Geotechnical Society Special Publication 5 (2) (2017) 22–27.
  * [28] N. Otsu, A threshold selection method from gray-level histograms, Automatica 11 (285-296) (1975) 23–27.
  * [29] Y. Higo, F. Oka, R. Morishita, Y. Matsushima, T. Yoshida, Trinarization of $\mu$x-ray ct images of partially saturated sand at different water-retention states using a region growing method, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 324 (2014) 63–69.
  * [30] C. Müssel, F. Schmid, T. J. Blätte, M. Hopfensitz, L. Lausser, H. A. Kestler, Bitrina—multiscale binarization and trinarization with quality analysis, Bioinformatics 32 (3) (2016) 465–468.
  * [31] E. Miahi, S. A. Mirroshandel, A. Nasr, Genetic neural architecture search for automatic assessment of human sperm images, Expert Systems with Applications 188 (2022) 115937.
  * [32] S. Javadi, S. A. Mirroshandel, A novel deep learning method for automatic assessment of human sperm images, Computers in biology and medicine 109 (2019) 182–194.
  * [33] O. L. Şavkay, E. Cesur, M. E. Yalçın, V. Tavşanoğlu, Sperm morphology analysis with cnn based algorithms, in: 2014 14th International Workshop on Cellular Nanoscale Networks and their Applications (CNNA), IEEE, 2014, pp. 1–2.
  * [34] L. Prabaharan, A. Raghunathan, An improved convolutional neural network for abnormality detection and segmentation from human sperm images, Journal of Ambient Intelligence and Humanized Computing 12 (3) (2021) 3341–3352.
  * [35] M. A. Riegler, M. H. Stensen, O. Witczak, J. M. Andersen, S. Hicks, H. L. Hammer, E. Delbarre, P. Halvorsen, A. Yazidi, N. Holst, et al., Artificial intelligence in the fertility clinic: status, pitfalls and possibilities, Human Reproduction 36 (9) (2021) 2429–2442.
  * [36] C. Dai, Z. Zhang, G. Shan, L.-T. Chu, Z. Huang, S. Moskovtsev, C. Librach, K. Jarvi, Y. Sun, Advances in sperm analysis: techniques, discoveries and applications, Nature Reviews Urology 18 (8) (2021) 447–467.
  * [37] W. H. Organization, et al., Who laboratory manual for the examination and processing of human semen (2010).
  * [38] A. A. Taha, A. Hanbury, Metrics for evaluating 3d medical image segmentation: analysis, selection, and tool, BMC medical imaging 15 (1) (2015) 1–28.
  * [39] J. Guo, C. He, X. Zhang, Nonlinear edge-preserving diffusion with adaptive source for document images binarization, Applied Mathematics and Computation 351 (2019) 8–22.


# Learning outside the Black-Box:
The pursuit of interpretable models

Jonathan Crabbe
University of Cambridge
<EMAIL_ADDRESS>
&Yao Zhang
University of Cambridge
<EMAIL_ADDRESS>
William R. Zame
University of California Los Angeles
<EMAIL_ADDRESS>
&Mihaela van der Schaar
University of Cambridge
University of California Los Angeles
The Alan Turing Institute
<EMAIL_ADDRESS>
###### Abstract

Machine Learning has proved its ability to produce accurate models – but the
deployment of these models outside the machine learning community has been
hindered by the difficulties of interpreting these models. This paper proposes
an algorithm that produces a continuous global interpretation of any given
continuous black-box function. Our algorithm employs a variation of projection
pursuit in which the ridge functions are chosen to be Meijer G-functions,
rather than the usual polynomial splines. Because Meijer G-functions are
differentiable in their parameters, we can “tune” the parameters of the
representation by gradient descent; as a consequence, our algorithm is
efficient. Using five familiar data sets from the UCI repository and two
familiar machine learning algorithms, we demonstrate that our algorithm
produces global interpretations that are both highly accurate and parsimonious
(involve a small number of terms). Our interpretations permit easy
understanding of the relative importance of features and feature interactions.
Our interpretation algorithm represents a leap forward from the previous state
of the art.

##  1 Introduction

What do we need to trust the predictions of the black-box models crafted by
our Machine Learning (ML) algorithms? Although the ML community has succeeded
in generating accurate models in a very wide variety of settings, it has not
yet succeeded in convincing most practitioners to adopt these models. One
possible reason is that many practitioners will use familiar models over more
accurate models that they do not understand. For a striking example, we cite
the area of medical risk prediction, in which clinical models have remained
the standard despite the fact that ML models are demonstrably more accurate.
The explicitly stated reason for this is that an acceptable model must be both
accurate and transparent [24]; state-of-the-art clinical models are
transparent while state-of-the-art ML models are not. In responding to this
challenge, one possibility would be to focus on the design of ML models that
are themselves transparent; an alternative possibility is to make a given ML
model more transparent by providing an interpretation of that model.

As detailed in Section 5, a substantial literature in the ML community follows
the latter approach, but not entirely successfully. To understand the
challenge that we are addressing here, it is instructive to make a simple
thought experiment. Let us assume that the data follows a Cox proportional
hazards model [10]. For such a model, the hazard rate at time $t$ is given by
$H(t\mid
x_{1},\dots,x_{d})=h_{0}(t)\exp\big{(}\sum_{i=1}^{d}b_{i}x_{i}\big{)}$, where
$h_{0}(t)$ is a baseline hazard, which potentially depends on time,
$x_{1},\dots,x_{d}$ are the features and the coefficients
$b_{1},\dots,b_{d}\in\mathbb{R}$ represent the importance of each feature.
Suppose that we are interested in identifying the part of this model which
only depends on the features. This part is trivially given by
$f(x)=\exp\left(\sum_{i=1}^{d}b_{i}x_{i}\right)$. We shall here consider this
as the “black-box” that we would like to identify by using our
interpretability method. If our interpretability method produces linear
models, as it is the case for _LIME_ [37], we could obtain an estimator
$\hat{f}(x)=1+\sum_{i=1}^{d}b_{i}x_{i}$ for $f$. On the other hand, a
polynomial estimator [2] for $f$ would be of the form
$\hat{f}(x)=1+\nicefrac{{\sum_{i=1}^{d}b_{i}x_{i}}}{{1!}}+\nicefrac{{(\sum_{i=1}^{d}b_{i}x_{i})^{2}}}{{2!}}+\textit{finite
number of terms}$. It goes without saying that neither of those interpretable
models captures an exponential relationship globally, hence precluding a
global interpretation of the black-box model. Would it be possible to overcome
these limitations of the state-of-the-art algorithms for interpretability?
More precisely, is there a way to build a regressor allowing to capture a
large class of interpretable functions _globally_? Needless to say that a
positive answer to these questions would have a tremendous impact on the
interpretability landscape in Machine Learning. To tackle this problem, we
build a new approach on top of the two following corner stones.

Projection Pursuit. The projection pursuit algorithm is widely known in the
statistics literature [13, 20, 39, 22]. Projection pursuit builds an
approximation to a given function by proceeding in stages. In each stage, the
algorithm finds a direction in the feature space and a ridge function of that
direction that best approximate the residual between the given function and
the approximation constructed in the previous state. The process continues
until the desired degree of accuracy is achieved. In projection pursuit, the
ridge functions are typically polynomial splines [18]. This reliance on
polynomial splines means that prediction pursuit suffers from the very
shortcoming noted above: it permits only local interpretations and hence does
not suggest any global structure. For this reason, we should work with another
family of functions.

Meijer G-functions. Instead of polynomial splines, we use functions from the
class of Meijer G-functions. This class of functions includes all of the
familiar functions used in modeling (polynomial, exponential, logarithmic,
trigonometric and hypergeometric functions). A Meijer G-function is defined by
four non-negative integer hyperparameters and an array of real parameters. An
interesting property of these functions is that we can efficiently compute a
numerical gradient of the Meijer G-functions with respect to these parameters.
This allows to use gradient-based optimization [1].

Contribution. In this paper, we introduce a new algorithm called _Symbolic
Pursuit_ that produces an interpretable model for a given black-box. As in a
projection pursuit, our algorithm gradually adds more terms in the model until
the desired precision is achieved. In our analysis, we address a major
difficulty introduced by using G-functions as ridge functions : the
hyperparameters. As aforementioned, each G-function has four hyperparameters
to tune. Fortunately, we show that there exist a set of five hyperparameter
configurations that covers most familiar functions (polynomial, rational,
exponential, trigonometric and hypergeometric functions). By restricting to
these configurations, each G-function is optimized efficiently over a
sufficiently large class of functions. Consequently, we demonstrate that our
Symbolic Pursuit algorithm allows to produce highly accurate global models
(that we call _symbolic models_) for black-boxes with a small number of
G-functions. With our algorithm, one to five G-functions are enough to
approximate a black-box model fitting a real world dataset. This is a leap
forward from the previous state-of-the-art.

To make the paper self-contained, we start with a short presentation of the
projection pursuit algorithm and Meijer G-functions in Section 2. We address
the hyperparameter optimization problem of Meijer G-functions in Section 3.
After these theoretical considerations, we assemble all the pieces to
construct the Symbolic Pursuit algorithm in Section 4. We compare this
algorithm to the state-of-the-art interpretability methods in Section 5.
Finally, we demonstrate that our algorithm produces globally symbolic models
with outstanding approximation power and parsimonious expressions on five
real-world datasets in Section 6.

##  2 Mathematical preliminaries

We begin by recalling the projection pursuit algorithm, which we shall reshape
for our purpose in Section 4. Then we review the definition of Meijer
G-functions that we use as the building blocks of the Symbolic Pursuit
algorithm.

###  2.1 Projection Pursuit

Let $\mathcal{X}=[0,1]^{d}$ be a normalized _feature space_ ,
$\mathcal{Y}\subseteq\mathbb{R}$ be a _label space_ and
$f:\mathcal{X}\rightarrow\mathcal{Y}$ be a _black-box model_ (e.g. a neural
network). Because we view $f$ as a black-box, we do not know the true form of
$f$, but we assume we can evaluate $f(x)$ for every $x\in\mathcal{X}$; i.e. we
can query the black-box. The projection pursuit algorithm [13] constructs an
_explicit_ model $\hat{f}$ for $f$ of the form

$\hat{f}(x)=\sum_{k=1}^{K}g_{k}\left(v_{k}^{\top}x\right),$ (1)

where $v_{k}\in\mathbb{R}^{d}$ is a vector (onto which we project the feature
vector $x\in\mathcal{X}$) and each $g_{k}$ is a function belonging to a
specified class $\mathcal{F}$ of univariate functions. Because the functions
$x\mapsto g_{k}\left(v_{k}^{\top}x\right)$ are constant on the hyperplanes
normal to $v_{k}$, they are usually called _ridge functions_.

The approximation is built following an iterative process. At stage
$j\in\mathbb{N}^{*}$, we suppose that the approximation $\hat{f}_{j}$ of the
above form, but having only $j-1$ terms, has been constructed; i.e.,
$\hat{f}_{j}(x)=\sum_{k=1}^{j-1}g_{k}\left(v_{k}^{\top}x\right)$111If $j=1$
then $\hat{f}_{1}\equiv 0$.. We define the residual $r_{j}=f-\hat{f}_{j}$ and
find a vector $v_{j}$ and a function $g_{j}\in\mathcal{F}$ to minimize the
$L^{2}$ norm of the updated residual:

$(g_{j},v_{j})=\arg\min_{\mathcal{F}\times\mathbb{R}^{d}}\int_{\mathcal{X}}\
\left[r_{j}(x)-g\left(v^{\top}x\right)\right]^{2}\,dF(x)$ (2)

where $F$ is the cumulative distribution function on the feature space. This
process is continued until the $L^{2}$ norm of the residual is smaller than a
predetermined threshold. This algorithm is typically complemented with a back-
fitting strategy in practice [17].

In principle, projection pursuit would seem an excellent approach to
interpretability because each stage provides one term in the approximation and
the process can be terminated at any stage. Hence projection pursuit allows a
trade-off between the accuracy of the representation $\hat{f}$ and the
complexity of $\hat{f}$. In practice, however, the class $\mathcal{F}$ is
usually taken to be the class of polynomial splines – often of low degree
(e.g. cubic) [20]. This has the advantage of making each stage of the
algorithm computationally tractable, but the disadvantage of often leading to
representations that involve many terms. Moreover, often even the individual
terms in these representations do not have natural interpretations. These
considerations suggest looking for a different class of functions that retain
computational tractability while leading to representations that are more
parsimonious and interpretable.

###  2.2 Meijer G-functions

Is there a set of functions $\mathcal{F}$ that includes most interpretable
functions and that would allow solving the optimization problem (2) with
standard techniques? The answer to this question is yes. The set of Meijer
G-function, that we denote $\mathcal{G}$, fulfils these two requirements, as
discussed in [1]. We shall thus henceforth restrict our investigations to
$\mathcal{F}=\mathcal{G}$. Here, we briefly recall the definition of a Meijer
G-function [5].

######  Definition 2.1 (Meijer G-function [5]).

A Meijer G-function is defined by an integral along a path $\mathcal{L}$ in
the complex plane,

$G_{p,q}^{m,n}\left(\begin{tabular}[]{c|c}$a_{1},\dots,a_{p}$&\hbox{\multirowsetup$z$}\\\
$b_{1},\dots,b_{q}$\end{tabular}\right)=\frac{1}{2\pi i}\int_{\mathcal{L}}\
z^{s}\frac{\prod_{j=1}^{m}\Gamma(b_{j}-s)\prod_{j=1}^{n}\Gamma(1-a_{j}+s)}{\prod_{j=m+1}^{q}\Gamma(1-b_{j}+s)\prod_{j=n+1}^{p}\Gamma(a_{j}-s)}ds,$
(3)

where $a_{i},b_{j}\in\mathbb{R}\ ;\ \forall i=1,\dots,p\ ;\ j=1,\dots,q\ ;\
m,n,p,q\in\mathbb{N}$ with $m\leq q$ , $n\leq p$ and $\Gamma$ is Euler’s Gamma
function. The path of integration $\mathcal{L}$ is chosen so that the poles
associated to the two families of Gamma functions (one family for the $a$’s
and one family for the $b$’s) lie on different sides of $\mathcal{L}$.222We
omit some technical details and restrictions; see [5] for details. The
definition yields a complex-analytic function in the entire complex plane
$\mathbb{C}$ with the possible exception of the origin $z=0$ and the unit
circle $\\{z\in\mathbb{C}:|z|=1\\}$. We will only consider its behavior for
real $z$ in the open unit interval $(0,1)\subset\mathbb{R}$. We write
$\mathcal{G}$ for the class of all Meijer G-functions.

This set of function is interesting because it includes most familiar
functions such as exponential, trigonometric and hypergeometric functions as
particular cases [31]. For example, the negative exponential function is

$\exp(-x)=G^{1,0}_{0,1}\left(\begin{tabular}[]{c|c}---&\hbox{\multirowsetup$x$}\\\
$0$&\end{tabular}\right).$

Furthermore, we can compute a numerical gradient of these function with
respect to the real parameters $a_{i},b_{j}$, hence allowing the use of
gradient-based optimization. In the following, it will be useful to denote
$\mathcal{G}^{m,n}_{p,q}$ the set of Meijer G-functions of hyperparameters
$m,n,p,q$. The full set of Meijer G-function can thus be decomposed as

$\mathcal{G}=\bigcup_{p,q\in\mathbb{N}}\bigcup_{n=0}^{p}\bigcup_{m=0}^{q}\mathcal{G}^{m,n}_{p,q}.$
(4)

It should already be obvious at this stage that we won’t be able to search in
all of these subsets of $\mathcal{G}$ so that $m,n,p,q$ will have to be fixed
as hyperparameters. How much are we restricting the possibilities by doing so?
Is there a clever restriction that would include most familiar functions? We
will elaborate on these questions in Section 3.

##  3 Hyperparameters climbing down the trees

In this section, we address a major theoretical challenge related to the use
of Meijer G-functions: the number of hyperparameters. More precisely, we have
4 hyperparameters $(m,n,p,q)\in\mathbb{N}^{4}$ to tune for each Meijer
G-function appearing in the expansion (1). It should be obvious that we cannot
search across every subclass ${\mathcal{G}}^{m,n}_{p.q}$. Fortunately, this is
not necessary; we can choose a set $\mathbb{H}$ of 4-tuples of hyperparameters
that is large enough that the subclasses
${\mathcal{G}}^{m,n}_{p.q};(m,n,p,q)\in{\mathbb{H}}$ encompass sufficiently
many functions but small enough that searching within and across these
subclasses is computationally tractable. Our choice of $\mathbb{H}$ relies on
two following results.

First, we note that the subsets $\mathcal{G}^{m,n}_{p,q}$ of $\mathcal{G}$
associated with different values of the hyperparameters are not disjointed.
Indeed, we show that the following result holds:

######  Lemma 3.1.

For all $(m,n,p,q)\in\mathbb{N}^{4}$ such that $p,q\geq 1$:

  * •

If $m\geq 1$ : $\mathcal{G}^{m-1,n}_{p-1,q-1}\subset\mathcal{G}^{m,n}_{p,q}$

  * •

If $n\geq 1$ : $\mathcal{G}^{m,n-1}_{p-1,q-1}\subset\mathcal{G}^{m,n}_{p,q}$

###### Proof.

A detailed proof can be found in Section 1 of the supplementary material. ∎

Figure 1: The tree of inclusions starting from $\mathcal{G}^{1,1}_{2,2}$.

This important lemma tells that, during the optimization, we can explore
different values of $(m,n,p,q)$ than the one initially fixed. To illustrate,
consider the 4-tuple $(m,n,p,q)=(1,1,2,2)$. Lemma 3.1 implies that
${\mathcal{G}}^{1,1}_{2,2}$ contains both ${\mathcal{G}}^{0,1}_{1,1}$ and
${\mathcal{G}}^{1,0}_{1,1}$, and that each of these in turn contains
${\mathcal{G}}^{0,0}_{0,0}$, as in Figure 1. Therefore, starting with
$(m,n,p,q)=(1,1,2,2)$ also allows to explore
$(m,n,p,q)=(0,1,1,1),(1,0,1,1),(0,0,0,0)$ at the same time. This concept can
be generalized so that the subsets of $\mathcal{G}$ can be represented in
terms of trees of inclusion, such as the one from Figure 1. Building on this
idea of trees of inclusion, we propose to choose a clever finite set of
configurations for the hyperparameters which allows to recover most closed
form expressions.

######  Proposition 3.1.

Consider the set of Meijer G-functions of the form

$\hat{f}(z)=G_{p,q}^{m,n}\left(\begin{tabular}[]{c|c}$a_{1},\dots,a_{p}$&\hbox{\multirowsetup$s\cdot
z^{r}$}\\\ $b_{1},\dots,b_{q}$\end{tabular}\right),$ (5)

where $a_{1},\dots,a_{p},b_{1},\dots,b_{q}\in\mathbb{R}$ ; $r,s\in\mathbb{R}$
and the hyperparameters belong to the configuration set
$(m,n,p,q)\in\mathbb{H}=\left\\{(1,0,0,2),(0,1,3,1),(2,1,2,3),(2,2,3,3),(2,0,1,3)\right\\}$.
This set of function includes all the functions with the form

$f(z)=\Phi(w\cdot z^{q})\cdot z^{t},$ (6)

with $w,q,t\in\mathbb{R}$ ;
$\Phi\in\left\\{\emph{id},\sin,\cos,\sinh,\cosh,\exp,\log(1+\cdot),\arcsin,\arctan,J_{\nu},Y_{\nu},I_{\nu},\frac{1}{1+\cdot},\Gamma\right\\}$
where $J_{\nu},Y_{\nu},I_{\nu}$ are the Bessel functions and $\Gamma$ is
Euler’s Gamma function.

######  Remark.

The above set of Meijer G-function includes much more function than the one
depicted in (6). The purpose of the above proposition is only to suggest the
generality of the symbolic model that we can assemble by using these Meijer
G-functions as building blocks. Here, we note that our choice allows to cover
many exponential, trigonometric, rational and Bessel functions.

###### Proof.

A detailed proof can be found in Section 1 of the supplementary material. ∎

This theorem gives a very satisfying prescription for the hyperparameters.
Building on this, we now have a realistic restriction of $\mathcal{G}$ that we
can use for the set $\mathcal{F}$ appearing in (2). We shall henceforth use
the following notation for this restriction:

$\mathcal{G}_{\mathbb{H}}=\bigcup_{(m,n,p,q)\in\mathbb{H}}\mathcal{G}^{m,n}_{p,q}.$
(7)

Therefore, our Symbolic Pursuit algorithm will search over the set of
functions $\mathcal{F}=\mathcal{G}_{\mathbb{H}}$. All the ingredients are now
ready to build the algorithm, this is the subject of next section.

##  4 Symbolic Pursuit

In this section, we build our _Symbolic Pursuit_ algorithm. With all the
ingredients we have prepared in the previous sections, this will be
straightforward. The starting point is naturally to consider the functions
$g_{1},\dots,g_{K}$ appearing in (1) to be elements of
$\mathcal{G}_{\mathbb{H}}$. However, we have to be careful because Meijer
G-functions might not be defined when $z=0,1$. By looking at (1), we note that
the arguments of each function $g_{k}$ takes the form $v_{k}^{\top}x$ for
$k=1,\dots,K$. We shall now scale this linear combination so that it stays in
the domain $(0,1)$ of the Meijer-G function $g_{k}$.

The Cauchy-Schwartz inequality for the $l^{2}$ inner product guarantees that
$|v_{k}^{\top}x|\leq\|v_{k}\|.\|x\|$. Our normalization for
$x\in\mathcal{X}=[0,1]^{d}$ guarantees that $\|x\|\leq\sqrt{d}$. By mixing
these two ingredients, we get

$\frac{|v_{k}^{\top}x|}{\|v_{k}\|\sqrt{d}}\leq 1.$ (8)

Therefore, we can simply take the _ReLU_ of
$\tfrac{v_{k}^{\top}x}{\|v_{k}\|\sqrt{d}}$ as an admissible argument for the
G-function $g_{k}$333This is not formally true since the inequality in (8) is
not strict. However, this is not important in practice since we can perfectly
normalize the features to a closed interval included in (0,1).. In conclusion,
we make the following replacement in (1):

$v_{k}^{\top}x\xrightarrow{}\left(\frac{v_{k}^{\top}x}{\|v_{k}\|\sqrt{d}}\right)^{+}\equiv\max\left(0,\frac{v_{k}^{\top}x}{\|v_{k}\|\sqrt{d}}\right).$
(9)

Note that, in this way, the argument remains a linear combination of the
features in the region where this linear combination is positive444Which would
not be the case if we had used a sigmoid to restrict the range of the
argument..

Because our optimization problem is non-convex, it is helpful to allow our
pursuit algorithm to correct the output of previous iterations at each
iteration via a back-fitting strategy. A conventional way to implement this is
to add a weight $w_{k}$ in front of each term in the expansion (1) and
optimize this weight during the back-fitting procedure [15]. This will
guarantee that terms that are constructed at some iteration of the algorithm
and found to be largely irrelevant at some later iteration will be assigned a
small weight. We can now rewrite the symbolic model (1) as

$\hat{f}(x)=\sum_{k=1}^{K}w_{k}\cdot
g_{k}\left(\left[\frac{v_{k}^{\top}x}{\|v_{k}\|\sqrt{d}}\right]^{+}\right).$
(10)

Similarly, we can reformulate the optimization problem (2) as

$(g_{k},v_{k},w_{k})=\operatorname*{arg\,min}_{\mathcal{G}_{\mathbb{H}}\times\mathbb{R}^{d}\times\mathbb{R}}\int_{\mathcal{X}}\
\left[r_{k}(x)-w\cdot
g\left(\left[\frac{v^{\top}x}{\|v\|\sqrt{d}}\right]^{+}\right)\right]^{2}\,dF(x).$
(11)

Note that the optimization over $g_{k}$ takes into account the five
hyperparameters configuration. The objective function in (11) is solved by
trying each configuration and keeping the one associated to the smallest loss.
At each step of the projection pursuit algorithm, we shall use a back-fitting
strategy to correct all the terms that already appear in the expansion.
Keeping this in mind, the back-fitting for the term
$l\in\left\\{1,\dots,k-1\right\\}$ will consist in minimizing the residue
which excludes the contribution of term $l$ at iteration $k$

$r_{k,l}(x)\equiv f(x)-\sum_{\begin{subarray}{c}j\neq
l\end{subarray}}^{k}w_{j}\cdot
g_{j}\left(\left[\frac{v_{j}^{\top}x}{\|v_{j}\|\sqrt{d}}\right]^{+}\right).$
(12)

The pseudo code of Symbolic Pursuit is provided in Section 5 of the
supplementary material, in which we write our algorithm that solves this
optimization problem step by step. It should be stressed that, on a
mathematical ground, using Meijer G-functions in optimization problems is far
from anodyne. We elaborate on some theoretical impact of this choice on the
behavior of the loss function in Section 3 of the supplementary material.

##  5 Related work

Table 1: Symbolic Pursuit and the state-of-the-art interpretability methods.

Algorithm Feature Importance Feature Interaction Model Independent Global
Parsimonious $\text{GA}^{2}\text{M}$ [30] LIME [37] SHAP [32] DeepLIFT [40]
L2X [7] NIT [42] INVASE [45] Symb. Metamodel [1] Symbolic Pursuit

In this section, we compare our algorithm to some of the most popular state-
of-the-art interpretability methods, this discussion is summarized in Table 1.
Like most methods, our Symbolic Pursuit algorithm allows to learn about
feature importance and feature interaction for each prediction. We shall give
more details about these two points in Section 6. Unlike methods such as
_DeepLIFT_ [40] or _NIT_ [42], our method is model independent (i.e. not
specialized for a particular class of black-box), as we shall demonstrate in
Section 6. However, unlike the other methods, which provide only local
information, Symbolic Pursuit also provides global information. This comes
from the fact that Meijer-G functions allow to capture a large class of closed
form expressions globally, as we detailed in Section 1 and 3. The only other
interpretive method that provides global information is the method of Symbolic
Metamodeling [1], which also makes use of Meijer G-functions. However, this
last method fails to produce parsimonious expressions since it produces an
additive model

$\hat{f}(x_{1},\dots,x_{d})=\sum_{i=1}^{d}g_{i}(x_{i})+\sum_{i=1}^{d}\sum_{j<i}g_{ij}(x_{i}\cdot
x_{j}),$ (13)

where the $g$’s are Meijer G-function. This model is made of
$d+\nicefrac{{d(d-1)}}{{2}}=\nicefrac{{d(d+1)}}{{2}}$ terms. For ten features
$d=10$, this corresponds to $55$ terms. The size of the Symbolic metamodel
enormously increases the burden of constructing interpretations, which makes
it hard to use in most practical situations. Moreover, only one hyperparameter
configuration is explored for all of these terms since all the terms are
optimized simultaneously. In the language that we have introduced in Section
3, this means that only one hyperparameter tree can be explored for all the
G-functions $g$’s. Indeed, exploring all the hyperparameter configurations
that we have identified in Section 3 would require to solve
$|\mathbb{H}|^{\nicefrac{{d(d+1)}}{{2}}}=5^{\nicefrac{{d(d+1)}}{{2}}}$
optimization problems corresponding to each choice of hyperparameter
configuration for each $g$. This is unrealistic when $d$ is large. Therefore,
it is not possible to identify all the familiar functions simultaneously,
since these are associated to different hyperparameter trees555For instance
$x\mapsto\exp(-x)\in\mathcal{G}^{1,0}_{0,1}$ and
$x\mapsto\ln(1+x)\in\mathcal{G}^{1,2}_{2,2}$ and these two sets cannot be put
in a same tree of inclusion via Lemma 3.1..

To solve these issues, we introduced the _Symbolic Pursuit_ algorithm which
benefits from the ability to produce parsimonious expressions. More precisely,
this method encapsulates perfectly the trade-off between accuracy and
interpretability as larger (and therefore less interpretable) expressions will
typically be more accurate. Therefore, parsimony naturally translates into
stopping the optimization when a reasonable accuracy has been achieved. We
will show in Section 6 that, even for real world datasets, this does not
require a large number of terms.

##  6 Experiments

###  6.1 Synthetic data

In this section, we evaluate the performance of Symbolic Pursuit on several
synthetic datasets 666The code for Symbolic Pursuit is available at
https://bitbucket.org/mvdschaar/mlforhealthlabpub and
https://github.com/JonathanCrabbe/Symbolic-Pursuit.. In order to compare the
performance of our method against popular interpretability methods, we shall
restrict to explanations in the form of feature importance. More precisely, we
focus on a pseudo black-box $f$ for which the feature importance is known
unambiguously. Here, we consider a three dimensional linear model

$\displaystyle f\left(x_{1},x_{2},x_{3}\right)=\beta_{1}\cdot
x_{1}+\beta_{2}\cdot x_{2}+\beta_{3}\cdot x_{3}=\beta^{\top}x,$

where $\beta=(\beta_{1},\beta_{2},\beta_{3})\in\mathbb{R}^{3}$ and
$x=(x_{1},x_{2},x_{3})\in[0,1]^{3}$. Since this pseudo black-box is merely a
linear model, we expect the interpretability methods to output an importance
vector of $\beta$ for each test point. Because only relative importance
matters in practice, we assume that all the importance vectors are normalized
so that $\beta\in S^{2}$ where $S^{2}$ is the 2-sphere.

Method | $\|\beta-\hat{\beta}\|^{2}$
---|---
LIME | $0.66\pm 0.07$
SHAP | $0.66\pm 0.07$
Symbolic | $0.02\pm 0.05$
Table 2: Symbolic Pursuit on synthetic data.

In our experiment, we start by drawing a true importance vector $\beta\sim
U([1,10]^{3})$ that we normalize subsequently. Then, we build a Symbolic
model, a _LIME_ explainer [37] and a _SHAP_ explainer [32] for $f$. Finally,
we draw 30 test points $x_{test}\sim U([0,1]^{3})$ that we input to the three
interpretability methods to build an estimator $\hat{\beta}$ for $\beta$. We
reproduce this experiment 100 times and report the MSE between $\hat{\beta}$
and $\beta$ for each interpretability method in Table 2. We note similar
performances of SHAP and LIME and a significant improvement with Symbolic
Models.

###  6.2 Real world data

In this section, we evaluate the performance of Symbolic Pursuit on two
popular black-box models – a Multilayer Perceptron (MLP) and Support Vector
Machine (SVM) – applied to five UCI datasets[12] including Wine Quality Red
(Wine) and Yacht Hydrodynamics (Yacht), Boston Housing (Boston), Energy
Efficiency (Energy) and Concrete Strength (Concrete). Both models are
implemented using the scikit-learn library [6] with the default
hyperparameters. For each experiment, we split the dataset into a training set
($80\%$) and a test set ($20\%$). The training set is used to produce the
Black-Box model and the Symbolic Model (via a mixup strategy for the later, as
detailed in Section 5 of the supplementary material); the test set is used to
test the model performance. We repeat the experiment five times with different
random splits and report averages and standard deviations. Table 3 shows the
mean squared error (MSE) and $R^{2}$ of the black-box models against the true
labels, the approximation MSE and $R^{2}$ of the symbolic models against the
corresponding black-box models, the MSE and $R^{2}$ of the symbolic model
against the true labels, and the number of terms in the symbolic models
constructed for each of the five random splits.777If the symbolic model is a
good approximation of the black-box model, it will necessarily have a similar
MSE and $R^{2}$ against the true labels, as shown in the Black-box and
Symbolic column of Table 3. The precise definition of the metrics we use here
can be found in Section 6 of the supplementary material. As can be seen in
Table 3 (Symbolic vs Black-box), all the symbolic models achieve low MSE and
high $R^{2}$ with respect to the black-box model. Indeed, in nine of ten cases
the $R^{2}$ is above $0.95$, and in the tenth it is above $0.90$. Moreover, in
15/50 instances, the symbolic model has only a single term, in 36/50 instances
the symbolic model has at most two terms and in 48/50 instances the symbolic
model has at most three terms. This provides strong evidence that Symbolic
Pursuit produces interpretations that are both accurate and parsimonious, as
previously asserted.

Table 3: Symbolic Pursuit results on the five UCI datasets [12].

Models Datasets Black-box Symbolic vs Black-box Symbolic # Terms MSE $R^{2}$
MSE $R^{2}$ MSE $R^{2}$ MLP Wine 0.016 $\pm$ .002 0.179 $\pm$ .061 0.001 $\pm$
.001 0.964 $\pm$ .018 0.016 $\pm$.002 0.179 $\pm$ .054 1, 1, 1, 2, 1 Yacht
1.433 $\pm$ .681 0.426 $\pm$ .233 0.008 $\pm$ .014 0.978 $\pm$.023 1.458 $\pm$
.653 0.413 $\pm$ .225 2, 2, 2, 1, 2 Boston 0.050 $\pm$ .015 0.681 $\pm$ .063
0.004 $\pm$ .001 0.952 $\pm$ .021 0.053 $\pm$ .016 0.660 $\pm$ .066 1, 3, 3,
3, 2 Energy 0.015 $\pm$ .001 0.926 $\pm$ .012 0.002 $\pm$ .001 0.988 $\pm$
.004 0.016 $\pm$ .003 0.918 $\pm$ .013 2, 2, 2, 1, 2 Concrete 0.100 $\pm$ .006
0.538 $\pm$ .061 0.001 $\pm$ .000 0.988 $\pm$ .003 0.100 $\pm$ .005 0.533
$\pm$ .056 1, 2, 3, 2, 3 SVM Wine 0.014 $\pm$ .001 0.331 $\pm$ .026 0.001
$\pm$ .001 0.904 $\pm$ .038 0.014 $\pm$ .001 0.301 $\pm$ .039 1, 2, 1, 1, 1
Yacht 0.723 $\pm$ .179 0.555 $\pm$ .036 0.010 $\pm$ .015 0.973 $\pm$ .043
0.737 $\pm$ .197 0.547 $\pm$ .050 3, 1, 2, 2, 2 Boston 0.040 $\pm$ .005 0.740
$\pm$ .053 0.001 $\pm$ .001 0.984 $\pm$ .017 0.043 $\pm$ .006 0.724 $\pm$ .046
3, 3, 2, 1, 2 Energy 0.015 $\pm$ .002 0.928 $\pm$ .007 0.002 $\pm$ .003 0.985
$\pm$ .016 0.018 $\pm$ .002 0.913 $\pm$.006 2, 3, 2, 3, 2 Concrete 0.069 $\pm$
.005 0.676 $\pm$ .039 0.003 $\pm$ .002 0.971 $\pm$ .015 0.082 $\pm$ .011 0.623
$\pm$ .043 3, 5, 3, 1, 5

In this setup, it is possible to do algebraic manipulations on a symbolic
model to extract transparent information of the black-box model. For instance,
we could simply inspect the components of each vector $v_{k}$ to get an idea
of the feature importance. We could also differentiate the G-function $g_{k}$
with respect to the features to obtain the gradients in closed form
expression. Most importantly, it is realistic to do these operations by hand
at this stage since the expressions are short, the Meijer G-functions can be
differentiated with respect to their argument easily [31] and their arguments
are linear combinations of features. To illustrate, we use one of the
interpretations of MLP on the Wine dataset; in order not to make the task too
easy, we use the split that produces two terms in the interpretation, rather
than any of the splits that produce a single term. With our notations, $f$
denotes the MLP black-box and $\hat{f}$ denotes its symbolic model. Because
$\hat{f}$ has two terms, it can be written as

$\displaystyle\hat{f}(x)$
$\displaystyle=w_{1}g_{1}\left(\frac{v_{1}^{\top}x}{\|v_{2}\|\sqrt{d}}\right)+w_{2}g_{2}\left(\frac{v_{2}^{\top}x}{\|v_{2}\|\sqrt{d}}\right).$
(14)

By construction, $g_{1},g_{2}\in{\mathcal{G}}_{\mathbb{H}}$ but in this
instance $g_{1},g_{2}$ do not appear to have expressions in terms of familiar
functions. Despite this, we can easily extract useful information from (14) by
building local polynomial models via a Taylor expansion of $g_{1}$ and
$g_{2}$, respectively. Let $z_{j}=\tfrac{v_{j}^{\top}x}{\|v_{j}\|\sqrt{d}}$
for $j=1,2$. Using the Taylor expansion of $g_{1}$ and $g_{2}$, we produce the
first order Taylor expansion of the symbolic model $\hat{f}(x)$ around an
instance $x$ from the test set:

$\begin{split}\hat{f}_{1}(x)=\tilde{v}^{\top}x+\tilde{c}=c_{0}+c_{1,1}z_{1}+c_{1,2}z_{2}.\end{split}$
(15)

In this particular instance $x$, we find that $\tilde{c}=0.8399$,
$c_{0}=1.001$, $c_{1,1}=-0.2339$, $c_{1,2}=0.3280$; the values of
$\tilde{v},v_{1},v_{2}$ are displayed in Figure 2 (b-d).

Let us now compare this linear model offered by $\hat{f}_{1}(x)$ with a local
linear surrogate model computed with _LIME_ [37]. LIME explainer suggests that
the most important features for the MLP model at $x$ are $x_{10}$ (alcohol),
$x_{9}$ (sulphates) and $x_{8}$ (pH) (See Figure 1 of the supplementary
material). The first order interpreter $\hat{f}_{1}(x)$ in (15) agrees with
this, since $x_{8}$, $x_{9}$ and $x_{10}$ have the highest weight in
$\tilde{v}$, as shown in Figure 2 (b). However, the agreement is not perfect
since LIME also suggests that $x_{5}$ (free sulfur dioxide) and $x_{6}$ (total
sulfur dioxide) are important and our interpreter does not. Unlike LIME, our
interpreter can easily provide a second-order Talor expansion of $g_{1}$ and
$g_{2}$ to suggest important interactions between features. The second-order
Taylor expansion of $\hat{f}(x)$ has the form
$\hat{f}_{2}(x)=\hat{f}_{1}(x)+c_{2,1}z_{1}^{2}+c_{2,2}z_{2}^{2}\approx\hat{f}_{1}(x)+c_{2,1}(v_{1}^{\top}x)^{2}$.
In this case we find that $\tilde{c}=0.8399$ $c_{2,1}=0.0623$ and
$c_{2,2}=-0.0002$ so that the interactions appearing in $z_{1}^{2}$ are
important compared to the interactions in $z_{2}^{2}$. Hence we can see from
Figure 2 (c) that $x_{1}$ (fixed acidity) has important interactions with
features $x_{7}$ (density), $x_{9}$ and $x_{10}$ but not $x_{8}$, despite the
fact that $x_{8}$ itself is an important feature. This short discussion allows
to see that the major advantage of our method compared to LIME is that we can
go beyond a linear model to capture nonlinear properties of the black-box,
such as feature interactions.

Figure 2: Coefficients and feature weights in our symbolic model.

##  7 Conclusion

This paper has proposed an algorithm that produces a global interpretation of
any given continuous black-box function. Our algorithm employs a variation of
projection pursuit in which the ridge functions are chosen to be Meijer
G-functions, rather than the usual polynomial splines. A series of experiments
demonstrates that the interpretations produced by our algorithm are both
accurate and parsimonious, and that the interpretation yields more information
than is available from other methods. Because our method produces continuous
models, it may not be appropriate for the interpretation of discontinuous
black-box models such as tree-based models, at least without some adaptation
and/or qualification. Perhaps more importantly, although our interpretive
models are parsimonious, they frequently involve unfamiliar functions. We have
argued above that this is not necessarily a barrier to understanding, because
the explicit expressions can be easily used to produce local interpretations
in terms of linear functions or low-order polynomials that reveal which
features and which interactions between features are most important. In the
medical domain, for example, this information is often very valuable, but not
easily obtained.

## Broader Impact

As mentioned at the very beginning of this paper, the lack of easy
interpretation of ML models has proved a serious obstacle to their adoption –
despite their demonstrated accuracy. Because ML models are more accurate than
previous models, anything that makes ML models more widely used is likely to
have enormous positive affects in practice – simply by providing better
predictions. Even interpretations that provide only first-order information
will surely prove to be important – not least by allowing for easier
correction of measurement/recording errors. Two examples may illustrate. (1)
It is well-documented that “No-Fly” lists often flag the wrong people who
happen to have the same names as the right people, and that getting removed
from such lists can be a difficult task. (2) It is similarly well-documented
that measurement/recording errors are common in the calculation of credit
scores – but because those calculations are often quite opaque, such errors
are often left uncorrected, so some who should be approved for loans are
declined, and others who should be declined are approved .

This paper has offered a method for producing accurate and parsimonious
interpretations of black-box models. In no sense do we view this as providing
the final or best method for interpretation; indeed we view this work as
taking, along with [1], only the first few steps in a new and very promising
direction. We have already noted that, because the interpretations our
algorithm produces are continuous, our method may not be suitable for
interpretation of black-box models such as Decision Trees or Random Forests,
and may not be suitable for classification problems unless they are
transformed into regression problems by assigning probabilities instead of
decisions. In itself, this transformation is simple and unobjectionable, but
clients who expect “Buy” or “Sell” advice from their financial advisor may not
be happy with a recommendation to “Buy with probability $0.7$ and Sell with
probability $0.3$.”

Despite the fact that our symbolic models contain few terms, which is a
significant improvement compared to [1], we still have to deal with some
Meijer G-functions or hypergeometric functions that don’t reduce to familiar
expressions. If this family of functions is explicitly used in some scientific
communities, such as in physics [14], they are not likely to be deemed
interpretable by some practitioners. A possible way to deal with this issue
would be to enforce the cancellation between some zeroes and poles of the
Meijer G-functions. This can be done by adding a lasso penalty to the loss,
this approach is detailed in Section 4 of the supplementary material. We are
convinced that the Symbolic Pursuit algorithm opens up several interesting
research paths in the ML intepretability landscape. We hope that this paper
will convince the ML community to walk along them to explore this emerging
paradigm of interpretability.

## Acknowledgments

Jonathan Crabbe is supported by Aviva. Yao Zhang is supported by GSK. Mihaela
van der Schaar is supported by the Office of Naval Research (ONR), NSF
1722516.

## References

  * [1] Mihaela van der Schaar Ahmed M. Alaa.  Demystifying black-box models with symbolic metamodels.  In Neural Information Processing Systems, 2019.
  * [2] Kartik Ahuja, William Zame, and Mihaela van der Schaar.  Optimal piecewise local-linear approximations, 2018.
  * [3] George E Andrews, Richard Askey, and Ranjan Roy.  Special functions, volume 71.  Cambridge university press, 1999.
  * [4] Harry Bateman.  Higher Transcendental Functions [Volumes I-III], volume 1.  McGraw-Hill Book Company, 1953.
  * [5] Richard Beals and Jacek Szmigielski.  Meijer g–functions: A gentle introduction, 08 2013.
  * [6] Lars Buitinck, Gilles Louppe, Mathieu Blondel, Fabian Pedregosa, Andreas Mueller, Olivier Grisel, Vlad Niculae, Peter Prettenhofer, Alexandre Gramfort, Jaques Grobler, Robert Layton, Jake VanderPlas, Arnaud Joly, Brian Holt, and Gaël Varoquaux.  API design for machine learning software: experiences from the scikit-learn project.  In ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108–122, 2013.
  * [7] Jianbo Chen, Le Song, Martin J. Wainwright, and Michael I. Jordan.  Learning to explain: An information-theoretic perspective on model interpretation, 2018.
  * [8] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.
  * [9] Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis.  Modeling wine preferences by data mining from physicochemical properties.  Decision Support Systems, 47(4):547 – 553, 2009.  Smart Business Networks: Concepts and Empirical Evidence.
  * [10] D. R. Cox.  Regression models and life-tables.  Journal of the Royal Statistical Society. Series B (Methodological), 34(2):187–220, 1972.
  * [11] Persi Diaconis and Mehrdad Shahshahani.  On nonlinear functions of linear combinations.  SIAM J. Sci. Stat. Comput., 5(1):175–191, March 1984.
  * [12] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * [13] Jerome H. Friedman and Werner Stuetzle.  Projection pursuit regression.  Journal of the American Statistical Association, 76(376):817–823, 1981.
  * [14] Walter G Glöckle and Theo F Nonnenmacher.  Fox function representation of non-debye relaxation processes.  Journal of Statistical Physics, 71(3-4):741–757, 1993.
  * [15] Trevor Hastie, Jonathan Taylor, Robert Tibshirani, Guenther Walther, et al.  Forward stagewise regression and the monotone lasso.  Electronic Journal of Statistics, 1:1–29, 2007.
  * [16] Trevor Hastie and Robert Tibshirani.  Generalized additive models.  Statist. Sci., 1(3):297–310, 08 1986.
  * [17] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning: data mining, inference, and prediction.  Springer Science & Business Media, 2009.
  * [18] Trevor Hastie, Robert Tibshirani, and Martin Wainwright.  Statistical learning with sparsity: the lasso and generalizations.  CRC press, 2015.
  * [19] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2(5):359 – 366, 1989.
  * [20] Peter J Huber.  Projection pursuit.  The annals of Statistics, pages 435–475, 1985.
  * [21] Fredrik Johansson et al.  mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18), December 2013.  http://mpmath.org/.
  * [22] Lee K. Jones.  On a conjecture of huber concerning the convergence of projection pursuit regression.  Ann. Statist., 15(2):880–882, 06 1987.
  * [23] James Jordon, Jinsung Yoon, and Mihaela van der Schaar.  KnockoffGAN: Generating knockoffs for feature selection using generative adversarial networks.  In International Conference on Learning Representations, 2019.
  * [24] Michael W. Kattan, Kenneth R. Hess, Mahul B. Amin, Ying Lu, Karl G.M. Moons, Jeffrey E. Gershenwald, Phyllis A. Gimotty, Justin H. Guinney, Susan Halabi, Alexander J. Lazar, Alyson L. Mahar, Tushar Patel, Daniel J. Sargent, Martin R. Weiser, Carolyn Compton, and members of the AJCC Precision Medicine Core.  American joint committee on cancer acceptance criteria for inclusion of risk models for individualized prognosis in the practice of precision medicine.  CA: A Cancer Journal for Clinicians, 66(5):370–374, 2016.
  * [25] Pang Wei Koh and Percy Liang.  Understanding black-box predictions via influence functions, 2017.
  * [26] A. N. Kolmogorov.  On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition.  Dokl. Akad. Nauk SSSR, 114:953–956, 1957.
  * [27] Věra Kůrková.  Kolmogorov’s theorem and multilayer neural networks.  Neural Networks, 5(3):501 – 506, 1992.
  * [28] Pat Langley and Jan M Zytkow.  Data-driven approaches to empirical discovery.  Artificial Intelligence, 40(1-3):283–312, 1989.
  * [29] Zachary C. Lipton.  The mythos of model interpretability, 2016.
  * [30] Yin Lou, Rich Caruana, Johannes Gehrke, and Giles Hooker.  Accurate intelligible models with pairwise interactions.  In Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’13, page 623–631, New York, NY, USA, 2013. Association for Computing Machinery.
  * [31] Y.L. Luke.  The Special Functions and Their Approximations.  ISSN. Elsevier Science, 1969.
  * [32] Scott Lundberg and Su-In Lee.  A unified approach to interpreting model predictions, 2017.
  * [33] Telmo Menezes and Camille Roth.  Symbolic regression of generative network models.  Scientific Reports, 4, 09 2014.
  * [34] Sally C. Morton.  Interpretable exploratory projection pursuit.  In Connie Page and Raoul LePage, editors, Computing Science and Statistics, pages 470–474, New York, NY, 1992. Springer New York.
  * [35] Patryk Orzechowski, William La Cava, and Jason H. Moore.  Where are we now? a large benchmark study of recent symbolic regression methods.  In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO ’18, page 1183–1190, New York, NY, USA, 2018. Association for Computing Machinery.
  * [36] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [37] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  "why should i trust you?": Explaining the predictions of any classifier, 2016.
  * [38] Kelly Roach.  Meijer g function representations.  In Proceedings of the 1997 International Symposium on Symbolic and Algebraic Computation, ISSAC ’97, page 205–211, New York, NY, USA, 1997\. Association for Computing Machinery.
  * [39] Charles B. Roosen and Trevor J. Hastie.  Logistic response projection pursuit.  AT&T Bell Laboratories, Doc. BL011214-930806-09TM, Murray Hill, NJ, 1993.
  * [40] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences, 2017.
  * [41] David A. Sprecher.  On the structure of continuous functions of several variables.  Transactions of the American Mathematical Society, 115:340–355, 1965.
  * [42] Michael Tsang, Hanpeng Liu, Sanjay Purushotham, Pavankumar Murali, and Yan Liu.  Neural interaction transparency (nit): Disentangling learned interactions for improved interpretability.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 5804–5813. Curran Associates, Inc., 2018.
  * [43] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1. 0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [44] I-Cheng Yeh.  Modeling of strength of high-performance concrete using artificial neural networks.” cement and concrete research, 28(12), 1797-1808.  Cement and Concrete Research, 28:1797–1808, 12 1998.
  * [45] Jinsung Yoon, James Jordon, and Mihaela van der Schaar.  Invase: Instance-wise variable selection using neural networks.  In ICLR, 2019.
  * [46] Hongyi Zhang, Moustapha Cisse, Yann N. Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization, 2017.
  * [47] Xin Zhang, Armando Solar-Lezama, and Rishabh Singh.  Interpreting neural network judgments via minimal, stable, and symbolic corrections, 2018.