# An interpretable multiple kernel learning approach for the discovery of
integrative cancer subtypes

Nora K. Speicher 1 and Nico Pfeifer 1,2, 111to whom correspondence should be
addressed
1Department of Computational Biology and Applied Algorithmics, Max Planck
Institute for Informatics, Saarland Informatics Campus Saarbrücken, Germany
2Methods in Medical Informatics,
Department of Computer Science, University of Tübingen, Germany

Abstract: Due to the complexity of cancer, clustering algorithms have been
used to disentangle the observed heterogeneity and identify cancer subtypes
that can be treated specifically. While kernel based clustering approaches
allow the use of more than one input matrix, which is an important factor when
considering a multidimensional/manifold disease like cancer, the clustering
results remain hard to evaluate and, in many cases, it is unclear which piece
of information had which impact on the final result. In this paper, we propose
an extension of multiple kernel learning clustering that enables the
characterization of each identified patient cluster based on the features that
had the highest impact on the result. To this end, we combine feature
clustering with multiple kernel dimensionality reduction and introduce FIPPA,
a score which measures the feature cluster impact on a patient cluster.
Results: We applied the approach to different cancer types described by four
different data types with the aim of identifying integrative patient subtypes
and understanding which features were most important for their identification.
Our results show that our method does not only have state-of-the-art
performance according to standard measures (e.g., survival analysis), but,
based on the high impact features, it also produces meaningful explanations
for the molecular bases of the subtypes. This could provide an important step
in the validation of potential cancer subtypes and enable the formulation of
new hypotheses concerning individual patient groups. Similar analysis are
possible for other disease phenotypes.
Availability: Source code for rMKL-LPP with feature clustering and the
calculation of fFIPPA scores is available at
https://github.molgen.mpg.de/nora/fFIPPA.
Contact<EMAIL_ADDRESS>
##  1 Introduction

Despite the growing amount of multidimensional data available for cancer
patients, the optimal treatment for a new patient remains in many cases
unknown. In this setting, clustering algorithms are increasingly helpful to
identify cancer subtypes, i.e., homogeneous groups of patients that could
benefit from the same treatment. Most clustering algorithms, e.g., k-means
clustering (Hartigan, 1975) or spectral clustering (von Luxburg, 2007)
generate a hard clustering, i.e., each sample is assigned to one cluster,
which works well when the clusters are separated. When analyzing cancer data,
samples might lie between two or more clusters because they show
characteristics of each of them. In these cases, fuzzy or soft clustering
methods give additional information, since they assign to each patient a
membership probability for each cluster instead of a binary cluster
assignment. For the purpose of soft-clustering, fuzzy c-means (Dunn, 1974;
Bezdek, 1981) extends the prominent k-means clustering method. To handle the
increasing diversity in the data describing the same sample, fuzzy c-means can
be combined with multiple kernel learning approaches, either directly as shown
by Hsin-Chien Huang et al. (2012) or via application to the result of multiple
kernel based dimensionality reduction (Speicher and Pfeifer, 2015).

However, after the process of calculating kernel matrices for the input data,
the subsequent data integration and the clustering of the samples,
interpreting the results is difficult and the crucial features usually remain
obscure. This lack of interpretability also hinders the validation of the
identified groups in cases where the ground truth is not known.

We base our approach on the assumption that sample clusters are not defined by
all features but rather by a subset of features which is also the motivation
for the technique of biclustering or co-clustering, i.e., the simultaneous
clustering of samples and features (Hartigan, 1972). This assumption becomes
even more realistic with the increasing number of features within the data
types, e.g., for DNA methylation or gene expression data, that are being
measured nowadays for a cancer sample. Therefore, our approach combines
feature clustering with subsequent data integration and sample clustering,
thereby providing several advantages. First, the feature clustering identifies
groups of homogeneous features, as typically, even for patients in the same
cluster, not all measured features will behave similarly. Consequently, the
sample similarities obtained from each of these feature groups are less noisy
than the similarities one would obtain based on all the features. Second, due
to the learned kernel weights we are able to identify which feature clusters
were especially influential for each sample clustering allowing the further
characterization of the identified sample clusters.

#### Related work

Data integration has been combined with clustering algorithms in different
ways. One approach that is regularly used for cancer data is the so-called
late integration, where clusterings are combined after they have been
retrieved for each data source independently (The Cancer Genome Atlas, 2012).
However, signals that might be spread over distinct data types could be too
weak in each data source to have an impact on the individual clustering and
therefore, might not be present in the final result. To circumvent this
problem, different approaches integrate the data at an earlier time point,
this way performing intermediate integration.

iCluster creates the clustering using a joint latent variable model in which
sparsity of the features is imposed by using an $l_{1}$ norm or a variance-
related penalty (Shen et al., 2012). The authors applied their approach to a
dataset of 55 glioblastoma patients described by three data types. In a
different approach, Wang et al. (2014) extract the sample similarities from
each data type and generate sample similarity networks. These can then be
fused iteratively by using a message-passing algorithm until convergence to a
common integrated network.

Another possibility for data integration is multiple kernel learning whose
strength lies in its wide applicability in different scenarios, e.g., to
integrate data from arbitrary sources or to integrate distinct representations
of the same input data. Therefore, it has been combined with a number of
unsupervised algorithms in the last years. Yu et al. (2012) developed a
multiple kernel version of k-means clustering. The optimization of the weights
for the conic sum of kernels is a non-convex problem tackled by a bi-level
alternating minimization procedure. A different approach for multiview k-means
combines non-linear dimensionality reduction based on multiple kernels with
subsequent clustering of the samples in the projection space (Speicher and
Pfeifer, 2015). Kumar et al. (2011) extend spectral clustering such that it
allows the clustering based on multiple inputs. Here, data integration is
achieved by co-regularization of the kernel matrices which enforces the
cluster structures in the different matrices to agree. Moving one step
further, Gönen and Margolin (2014) introduced localized kernel k-means, which,
instead of optimizing one weight per kernel, uses per-sample weights to have
more flexibility to account for sample-specific characteristics. Some of the
algorithms described have been applied to cancer data and could show
improvements, such as increased robustness to noise in individual data sources
and increased flexibility of the model (Gönen and Margolin, 2014; Speicher and
Pfeifer, 2015). However, they do not allow the extraction of importances
directly related to specific features or the interpretation of the result on
the basis of the applied algorithm, which is especially in unsupervised
learning an important step to be able to associate clusters with their
specific characteristics.

A number of approaches aim at increasing performance and interpretability by
incorporating prior knowledge. Multiview biclustering (Le Van et al., 2016)
integrates copy number variations with Boolean mutation data, which are
transformed into ranked matrices using a known network of molecular
interactions between the genes. The subtype identification problem is then
formulated as a rank matrix factorization of the ranked data. On breast cancer
data, the method could refine the established subtypes which are based on
PAM50 classification. Approaches similar to ours concerning the subdivision
into feature groups before applying multiple kernel learning were presented by
Sinnott and Cai (2018) for supervised survival prediction of cancer patients
and by Rahimi and Gönen (2018) for discriminating early- from late-stage
cancers, showing a better interpretability compared to standard approaches.
However, in both cases the feature groups are identified using prior knowledge
(i.e., gene groups reflecting the memberships to pre-defined biological
pathways) and the known outcome of interest is used to train the models, which
is not the case for unsupervised clustering.

##  2 Approach

Due to the flexibility of the kernel functions, methods that can handle
multiple input kernel matrices allow the integration of very distinct data
types, such as numerical, sequence, or graph data, as well as the integration
of different preprocessings of the input data. Additionally, the individual
up- or downweighting of kernel matrices can account for differences in the
quality or relevance of the kernel matrices. However, the features within one
data type can be very heterogeneous and might not all be equally important to
a meaningful distance measure between the patients. Moreover, for non-linear
kernels in general, it is not straightforward how to identify feature
importances related to the result. Here, we propose a procedure that combines
feature clustering with sample clustering based on multiple kernels, thereby
increasing the potential to interpret the result without losing the power of
the multiple kernel learning approach.

Figure 1: Overview of the presented approach for two input data matrices.
First, feature clustering is performed, here with $C=3$. Each feature cluster
gives rise to one kernel matrix, which are integrated using rMKL-LPP. This
method optimizes one weight for each kernel matrix and a projection matrix,
leading to a low-dimensional representation of the samples. Using these
learned coordinates, the samples are clustered using fuzzy c-means. Finally,
the feature clusters, kernel weights per feature cluster, and the patient
clusters are used to calculate FIPPA scores, which describe feature cluster
impact on a patient cluster.

As illustrated in Figure 1, we cluster the features of each data type using
k-means such that we can generate one kernel matrix based on each feature
cluster. The kernel matrices are then integrated using a multiple kernel
learning approach, here regularized multiple kernel learning for locality
preserving projection. Based on the low-dimensional projection, we cluster the
samples using fuzzy c-means. Our approach provides two advantages compared to
standard procedures. First, increasing the homogeneity of the features by
identifying feature clusters can reduce the noise in each kernel matrix, this
way, a signal that is generated by only a few features can still significantly
influence the final result if these features behave very similarly over a
subset of samples. Second, the availability of the feature clusters and
respective kernel weights allows for increased interpretability of the
identified patient clusters as each cluster can be traced back to those groups
of features that had the highest influence on the sample similarities. To our
knowledge, this is the first extension of a multiple kernel clustering
algorithm towards integrative biclustering, where biclustering describes the
simultaneous clustering of features and samples, or of rows and columns of a
data matrix (Hartigan, 1972).

##  3 Materials and methods

###  3.1 Regularized multiple kernel learning for locality preserving
projections (rMKL-LPP)

Given a set of $M$ kernel matrices $\\{K_{1},...,K_{M}\\}$ that are based on
$N$ samples. Multiple kernel learning describes the optimization of a specific
weight $\beta$ for each kernel matrix, such that an ensemble kernel matrix
$\mathbb{K}$ can be calculated via:

$\mathbb{K}=\sum_{m=1}^{M}\beta_{m}K_{M},~{}~{}\text{with}~{}~{}\beta_{m}\geq
0.$ (1)

Regularized multiple kernel learning for locality preserving projections
(Speicher and Pfeifer, 2015) optimizes these kernel weights together with a
projection matrix such that local neighborhoods are preserved optimally when
reducing the dimensionality of the data. The projection matrix
$A\in\mathbb{R}^{p\times N}$ and the kernel weights
$\boldsymbol{\beta}\in\mathbb{R}^{M}$ are learned via the minimization problem

$\displaystyle\underset{\boldsymbol{A},\boldsymbol{\beta}}{\mbox{minimize
}}~{}$
$\displaystyle\sum_{i,j=1}^{N}\|\boldsymbol{A}^{T}\mathbb{K}^{(i)}\boldsymbol{\beta}-\boldsymbol{A}^{T}\mathbb{K}^{(j)}\boldsymbol{\beta}\|^{2}w_{ij}$
(2) subject to
$\displaystyle\sum_{i=1}^{N}\|\boldsymbol{A}^{T}\mathbb{K}^{(i)}\boldsymbol{\beta}\|^{2}d_{ii}=const.$
(3) $\displaystyle\|\boldsymbol{\beta}\|_{1}=1$ (4)
$\displaystyle\beta_{m}\geq 0,m=1,2,...,M.$ (5)

with

$\mathbb{K}^{(i)}=\begin{pmatrix}K_{1}(1,i)&\cdots&K_{M}(1,i)\\\
\vdots&\ddots&\vdots\\\ K_{1}(N,i)&\cdots&K_{M}(N,i)\\\
\end{pmatrix}\in\mathbb{R}^{N\times M}.$ (6)

The projection matrix $A$ consists of $p$ vectors
$[\boldsymbol{\alpha_{1}}\cdots\boldsymbol{\alpha_{p}}]$ that project the data
into $p$ dimensions. The weight matrix $W$ and the diagonal matrix $D$ are
determined by locality preserving projections using the $k$-neighborhood
$\mathcal{N}_{k}(j)$ of a sample $j$ as follows

$\displaystyle w_{ij}$
$\displaystyle=\begin{cases}1,~{}~{}\text{if}~{}i\in\mathcal{N}_{k}(j)\vee
j\in\mathcal{N}_{k}(i)\\\ 0,~{}~{}\text{else}\end{cases}$ (7)
$\displaystyle\text{and }~{}~{}~{}~{}~{}d_{ii}$ $\displaystyle=\sum_{i\neq
j}w_{ij}.$ (8)

The optimization of $A$ and $\beta$ is performed using coordinate descent. $A$
is optimized by solving a generalized eigenvalue problem, whereas the optimal
$\beta$ is identified using semidefinite programming.

###  3.2 Fuzzy c-means clustering (FCM)

Given a data matrix $X\in\mathbb{R}^{N\times l}$ describing $N$ samples
$x_{i}$ with $l$ features, FCM identifies $C$ cluster centers $v_{c}$ and
assigns the degrees of memberships $u_{i,c}$ for each sample $i$ and cluster
$c$ by minimization of the following objective function (Dunn, 1974; Bezdek,
1981):

$\displaystyle J(U,v)~{}=~{}$
$\displaystyle\sum_{i=1}^{N}\sum_{c=1}^{C}u_{i,c}^{f}\lVert
x_{i}-v_{c}\rVert^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(9) subject to $\displaystyle\sum_{c=1}^{C}u_{i,c}=1~{}~{}~{}~{}~{}\forall i$
$\displaystyle u_{i,c}\geq 0~{}~{}~{}~{}~{}\forall i,c$
$\displaystyle\sum_{i=1}^{N}u_{i,c}>0~{}~{}~{}~{}~{}\forall c.$

The resulting $U$ is an $N\times C$ matrix of the degrees of cluster
memberships $u_{i,c}$, which depends on $f\geq 1$, a parameter of the method
that controls the degree of fuzzification. Choosing $f=1$ results in the hard
clustering of the samples (i.e., $u_{i,c}\in\\{0,1\\}$); choosing
$f\rightarrow\infty$ results in uniform cluster probabilities
$u_{i,c}=\nicefrac{{1}}{{C}}$ for all $i,j$.

###  3.3 Increased interpretability due to simultaneous clustering of features
and samples

The impact of each feature cluster $m\in\\{1,..,M\\}$ on each identified
sample cluster $c\in\\{1,...,C\\}$ (FIPPAc,m) can be calculated based on the
kernel weights $\beta$ and the kernel matrices $K_{m}$ using the following
equation

$\text{FIPPA}_{c,m}=\dfrac{1}{|c|^{2}}\sum_{x_{i},x_{j}\in
c}\dfrac{\beta_{m}K_{m}[i,j]}{\mathbb{K}[i,j]},$ (10)

with $\mathbb{K}$ being the ensemble kernel matrix. Here, we suppose a hard
clustering of the patients, which can either be generated using a hard
clustering algorithm or using the modal class of a fuzzy clustering. However,
fuzzy clustering provides additional information concerning the reliability of
the cluster assignment for each sample. Using these probabilities can make the
results more robust given that some samples might have an ambiguous signature
and therefore lie between two or more clusters. Given a probability
$p_{c}(x_{i})=p(x_{i}\in c)$ that patient $x_{i}$ belongs to cluster $c$, the
fuzzy FIPPA can be calculated as follows:

$\text{fFIPPA}_{c,m}=\dfrac{1}{|N|^{2}}\sum_{i,j=1}^{N}p_{c}(x_{i})p_{c}(x_{j})\dfrac{\beta_{m}K_{m}[i,j]}{\mathbb{K}[i,j]}.$
(11)

The incorporation of the joint probability $p(x_{i}\in c\land x_{j}\in
c)=p_{c}(x_{i})p_{c}(x_{j})$ replaces the selection of sample pairs performed
in Equation 10 but still allows to downweight the effects of samples that are
uncertain or unlikely to belong to $c$.

Moreover, we can separate this overall score into a positive part
($\text{fFIPPA}^{+}_{k,m}$) that leads to high intra-cluster similarity and a
negative part that leads to high inter-cluster ($\text{fFIPPA}^{-}_{k,m}$)
dissimilarity. For this purpose, we define

$\mathbb{K}^{+}=\sum_{m=1}^{M}\beta_{m}K_{m}^{+}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}\mathbb{K}^{-}=\sum_{m=1}^{M}\beta_{m}K_{m}^{-}$
(12)

with $K_{m}^{+}$ being the positive part of the matrix $K_{m}$ (all negative
values set to zero) and vice versa for $K_{m}^{-}$. For this step, all kernel
matrices need to be centered in the feature space, such that the mean of each
matrix is equal to zero. When combining Formula 12 with Formula 11, the
calculation of the positive and negative fFIPPA is given by

$\displaystyle\text{fFIPPA}^{+}_{c,m}$
$\displaystyle=\dfrac{1}{|N|^{2}}\sum_{i,j=1}^{N}p_{c}(x_{i}\land
x_{j})\dfrac{\beta_{m}K^{+}_{m}[i,j]}{\mathbb{K^{+}}[i,j]},~{}~{}~{}~{}\text{and}$
$\displaystyle\text{fFIPPA}^{-}_{c,m}$
$\displaystyle=\dfrac{1}{|N|^{2}}\sum_{i,j=1}^{N}p_{c}(x_{i}\oplus
x_{j})\dfrac{\beta_{m}K^{-}_{m}[i,j]}{\mathbb{K^{-}}[i,j]}.$ (13)

Besides the fact that the two scores are based on the positive and the
negative part of the kernel matrices, respectively, the main difference
between them is the probability factor for each summand. For the positive
fFIPPA, the joint probability $p_{c}(x_{i}\land x_{j})$ is used to generate
high influence for pairs where both partners are likely to belong to cluster
$c$, while for the negative fFIPPA, the exclusive or, defined by

$p_{c}(x_{i}\oplus
x_{j})=(p_{c}(x_{i})+p_{c}(x_{j})-2p_{c}(x_{i})p_{c}(x_{j})),$ (14)

results in an increased factor for pairs of samples of which exactly one has a
high probability for $c$.

The fFIPPA scores calculated allow the identification of feature clusters that
contribute more than average to the similarity of the samples within a sample
cluster and the dissimilarity of the samples in two different clusters,
thereby, revealing the underlying basis of the generated clustering.

###  3.4 Materials

We applied our approach to six different cancer data sets generated by The
Cancer Genome Atlas which were downloaded from the UCSC Xena browser (Goldman
et al., 2017). The cancer types covered are breast invasive carcinoma (BRCA),
lung adenocarcinoma (LUAD), head and neck squamous cell carcinoma (HNSC),
lower grade glioma (LGG), thyroid carcinoma (THCA) and prostate adenocarcinoma
(PRAD). For each cancer patient, we used DNA methylation, gene expression
data, copy number variations, and miRNA expression data for clustering. DNA
methylation was mapped from methylation sites to gene promoter and gene body
regions using RnBeads (Assenov et al., 2014). Using GeneTrail (Stöckel et al.,
2016), we mapped the miRNAs to their target genes (i.e., the gene that is
regulated by the miRNA) such that the features of each data type were genes.
Due to the high number of features available, we performed the whole analysis
on the 10% of the features with the highest variance from each data type. For
the cluster evaluation, we further leveraged the survival times of the
patients for all cancer types except PRAD and THCA. For these two cancer
types, survival analysis would not be informative due to the low number of
events.

##  4 Results and discussion

###  4.1 Parameter selection

When applying our approach to a data set, we need to choose the number of
feature clusters per data type as well as the number of patient clusters. For
our experimental validation, we set both parameters to the same value
($c\in\\{2,...,6\\}$). Feature clustering was performed using k-means before
generating the kernel matrices using the Gaussian radial basis kernel
function. The parameter of the kernel $\gamma$ was chosen dependent on the
number of features $d$ in the respective feature set based on the rule of
thumb $\gamma=\frac{1}{2d^{2}}$ (Gärtner et al., 2002). We generated three
kernels per data type by multiplying this $\gamma$ with a factor
$f_{\gamma}\in\\{0.5,1,2\\}$ and only used the one kernel matrix providing the
highest variance in the first $d$ principal components. The number of
neighbors for rMKL-LPP was set to 9 and the the number of dimensions $d$ to 5,
as explained in our previous work (Speicher and Pfeifer, 2015). The
fuzzification degree $f$ of the soft-clustering algorithm was set to the
default value of $2$ in concordance with Dunn (1974). If necessary for the
subsequent analysis, we assigned each patient to its modal cluster (i.e., the
cluster with the highest probability), otherwise, we used the cluster
membership probabilities returned by FCM.

###  4.2 Robustness

Since each kernel matrix in this approach is constructed on the basis of a
feature cluster, slight changes in these initial clusters will propagate
through the algorithm and might have an effect on the final result. We
repeated the complete approach 50 times with different seeds to analyze the
robustness of the final patient clustering. Figure 2 shows the pairwise
similarities measured by the Rand index between the final patient clusterings.
When including all patients according to their modal class, we observe high
reproducibility for all cancer types except LUAD, for which the Rand index is
approx. 0.85. When using the class probabilities for each sample to exclude
patients where the prediction has a low confidence, here defined as being more
than one standard deviation lower than the mean, we observe that the cluster
assignments for the remaining patients are more stable for all cancer types
including LUAD.

Figure 2: Robustness of FC+rMKL-LPP in 50 repetitions. Dark grey boxes indicate Rand indices on the basis of all patients, light grey boxes indicate Rand indices calculated without low confidence samples, i.e., patients with a modal class probability $p_{c}$ more than one standard deviation smaller than the mean of the modal class probabilities. Table 1: For each cancer type and method, we report the most significant results when varying the number of clusters from two to six. Average kLPP stands for kernel locality preserving projection on the (unweighted) average kernel, rMKL-LPP for the standard multiple kernel learning approach with one kernel per data type, and FC + rMKL-LPP represents the proposed approach for which the reported p-values are the median of 50 runs. All three approaches are combined with subsequent fuzzy c-means clustering. $C$ indicates the chosen number of clusters. Cancer | average kLPP | rMKL-LPP | FC + rMKL-LPP
---|---|---|---
|  p-value |  $C$ |  p-value |  $C$ |  p-value |  $C$
BRCA |  3.7E-2 |  6 |  7.3E-2 |  6 |  5.0E-2 |  4
HNSC |  1.4E-3 |  6 |  1.4E-3 |  6 |  9.96E-3 |  5
LGG |  $<$1.0E-16 |  3:6 |  $<$1.0E-16 |  3:6 |  $<$1.0E-16 |  3:6
LUAD |  0.15 |  2 |  2.9E-2 |  2 |  3.1E-2 |  6

###  4.3 Survival analysis

To evaluate if a patient clustering could be clinically relevant, we performed
survival analysis on the results with $C\in\\{2,...,6\\}$ clusters for each
cancer type. To assess the influence of the feature clustering step, we also
generated patient clusterings with two established methods: kLPP with the
average kernel, and rMKL-LPP without feature clustering, which has been shown
to perform well in comparison to other data integration methods (Rappoport and
Shamir, 2018). To obtain a fair comparison, the same kernel matrices have been
used for all approaches.

The p-value for the log-rank test (Hosmer Jr. et al., 2011) indicates if there
is at least one group of patients among the identified ones with a significant
difference in the survival rate compared to the other groups. While the number
of clusters is reflected in the degrees of freedom that are used in the test,
we still corrected for the number of tests performed due to the variation of
$C$ using Benjamini-Hochberg correction. The results of the survival analysis
shown in Table 1 indicate that on the four cancer types where the number of
events allows for a survival analysis, our approach has a comparable
performance to the established methods and is able to find biologically
meaningful groups of patients. Moreover, despite the additional flexibility
which is achieved by the feature clustering, the optimal number of clusters
does not increase in comparison to less complex methods. However, as not all
patients receive an optimal treatment, the survival data can be biased.
Therefore, we need additional means to evaluate and interpret the underlying
structure of the identified clusterings.

###  4.4 Interpretation

As described in Section 3.3 the weights for each feature cluster in
combination with the kernel matrices allow us to calculate the fFIPPA scores.
Figure 3(a) visualizes the positive fFIPPA scores for the BRCA clustering with
$c=4$. Here, we can see clear differences between the patient clusters. While
for all of them, copy number variations have the strongest impact, the exact
contributions of the data types vary (e.g. the fFIPPA of gene expression,
especially feature cluster 3, is strongest in patient cluster 3). When
comparing to the negative fFIPPAs (cf. Figure 3(b)), we observe in general
similar patterns. Differences still occur, e.g., in patient cluster 2 and 4,
one can see that the fFIPPA of copy number variations on the intra-cluster
similarity is larger than their impact on inter-cluster dissimilarity.

(a) fFIPPA scores for intra-cluster similarity of the patient clusters.

(b) fFIPPA scores for inter-cluster dissimilarity of the patient clusters.

Figure 3: Positive and negative fFIPPA of each feature cluster and patients
cluster for BRCA with $C=4$. Each row represents one patient cluster, i.e., it
is derived using the probabilities for this specific cluster. The width of a
bar segment is defined by $\text{fFIPPA}_{c,m}^{+}$ or
($\text{fFIPPA}_{c,m}^{-}$) for the respective patient cluster $c$ and feature
cluster $m$.

Observing the impact of each data type is also possible with traditional
multiple kernel learning, however, in our approach, we can additionally
analyze the impact of the individual feature clusters. To gain further
insights into the underlying biological mechanisms, we used the positive and
negative fFIPPA scores to identify which feature clusters contributed more
than average to a patient cluster, i.e., we initialized two lists of feature
clusters, the first one explaining high dissimilarity to other patient
clusters, the second one explaining high similarity within the patient
cluster. After identification of these high-fFIPPA feature clusters, we
filtered the associated features to keep only those that behave homogeneously
over the patients of the specific cluster. For the homogeneity filtering, we
introduced a homogeneity score, counting in how many patients of the
respective patient cluster the gene was constantly expressed (i.e., either
over- or underexpressed, whichever happened more often). Based on these
scores, we obtained a distribution representing the homogeneity of the
features within a cluster. For further analysis, we kept every gene with a
score larger than the mean of this distribution and distinguish between over-
and underexpressed genes. This way, we generated four lists of features for
each patient cluster describing a) overexpressed genes leading to high intra-
cluster similarity, b) underexpressed genes leading to high intra-cluster
similarity, c) overexpressed genes leading to high inter-cluster dissimilarity
and d) underexpressed genes leading to high inter-cluster dissimilarity.

In these lists, we found genes with known implications to cancer formation or
progression. Some of the genes were related to several cancer types (e.g.
PTEN, which regulates the AKT/PKB signaling pathway and thereby represents a
tumor suppressor in various tumors (Cantley and Neel, 1999) was present in
lists for BRCA, THCA, and PRAD) or to several clusters of the same cancer type
(e.g.SOX2 for HNSC, which has been shown to control tumor initiation in
squamous-cell carcinomas (Boumahdi et al., 2014)). Others were only found for
one specific cluster of a cancer type (e.g., among the four clusters of BRCA,
only one gene list included GRPR, a G-protein coupled receptor that was shown
to influence the viability of breast cancer cells (Cornelio et al., 2013)).
From a biological perspective, this finding is consistent with our expectation
that the lists are not disjoint, since cancer does not form due to an
individual mutation but rather as a consequence of a number of different
aberrations, which can be shared between different subtypes or even cancer
types.

Besides looking at the lists at the level of individual genes, we also
inferred biological functions that are related to these gene lists using
GeneTrail (Stöckel et al., 2016). We applied over-representation analysis
which tests, if significantly more genes in the list are associated to a
specific category than would be expected based on a background gene list. As
background all the genes were used that entered the clustering and we tested
for enrichment of KEGG pathways and gene ontology categories (GO).

As there was some variation in the cluster assignments, one would also expect
some variation in the enriched terms. To quantify the stability of the result,
we compared the identified terms of the best to the terms of the median
result, where best and median was defined according to the survival p-value.
The similarity of the GO terms was determined using the relevance score
Schlicker et al. (2006), for the similarity between two sets of KEGG pathways
we used the Jaccard index. The average similarities correlate to the observed
stabilities of the cluster assignments: The value for LUAD is slightly lower
than for BRCA and HNSC (0.48 vs. 0.58 and 0.68, respectively), however, for
all three cancer types the scores indicate that similar functions or pathways
are found. For LGG, the results were stable such that the optimum and median
were equivalent, for PRAD and THCA, the best and median result could not be
identified due to the lack of adequate survival data.

Based on the intersection of the optimal and median result, Table 2 shows a
subset of the findings for HNSC for Cluster 3, which has the poorest prognosis
(median survival time 1079 days), and Cluster 4, which has the best prognosis
(median survival time $>4241$ days). The terms identified are very different
for the two clusters. In Cluster 3, we obtained enriched terms that are mainly
associated with muscle cells, more specifically, their development or their
activation. Skeletal muscle invasion seems to be correlated to the recurrence
of the disease (Chandler et al., 2011), which agrees with our observation that
patients having a distinct signature in the genes involved in these functions
have a poor survival.

Table 2: Significant GO terms and KEGG pathways associated with high intra-cluster similarity identified by over-representation analysis for patient clusters 3 and 4 for HNSC. For patient cluster 3, the terms are based on inactive genes, for patient cluster 4 on active genes. $c=3$ |  GO-Cellular component (p-value $<1.0E-3$)
---|---
|  contractile fiber
|  myofibril
|  sarcomere
|  I band
|  myosin II complex
|  muscle myosin complex
|  A band
|  myosin complex
|  Z disc
|  myofilament
|  GO - Biological Process (p-value $<1.0E-3$)
|  myofibril assembly
|  striated muscle cell development
|  muscle filament sliding&actin myosin filament sliding
|  striated muscle contraction
|  actomyosin structure organization
|  muscle cell development
$c=4$ |  GO -Biological process (p-value $<1.0E-3$)
|  regulation of peptidyl serine phosphorylation of STAT protein
|  natural killer cell activation involved in immune response
|  serine phosphorylation of STAT protein
|  response to exogenous dsRNA
|  natural killer cell activation
|  positive regulation of peptidyl serine phosphorylation
|  regulation of type I interferon mediated signaling pathway

For patient cluster 4, we identified terms related to the phosphorylation of a
STAT protein which regulates as a part of the JAK/STAT signaling pathway cell
growth and has been shown to play a role in tumor formation and progression
(Song and Grandis, 2000). On the other hand, the identified terms mainly refer
to functions of the immune system, including the regulation of type I
interferon mediated signaling pathway which is a regulation of the innate
immune response. Despite of the controversial role of inflammation and the
immune system in cancer (Coussens and Werb, 2002), its activation could help
eliminate cancer cells and therefore, be beneficial for the respective patient
group. For HNSC, the terms associated to high inter-cluster dissimilarity and
high intra-cluster similarity differ only for Cluster 1, where we identify
terms, such as B cell differentiation and proliferation, that are associated
to high dissimilarity to other clusters but not to high intra-cluster
similarity.

We also applied our approach to PRAD and THCA where survival analysis is not
meaningful because only a small number of patients die from the disease. For
both cancer types, we set $C$ to 4 and could identify different significant
terms for some patient clusters. For PRAD, we identified amongst others a
number of terms related to the activity of olfactory receptors, which have
been shown to participate in the process of tumor cell proliferation and
apoptosis (Chen et al., 2018). Several hits for THCA, e.g., DNA deamination,
indicate an influence of epigenetic regulation mechanisms, which is supported
by current literature (Rodríguez-Rodero et al., 2017).

##  5 Conclusion

The approach presented in this paper provides a step towards a better
interpretabilty of multiple kernel clustering. For a few kernel functions it
is possible to extract feature importances, however, our new procedure is
independent of the kernel function used and can consequently be applied in
combination with any kernel function. The application to biological data sets
shows that the approach gives state-of-the-art performance with respect to
survival analysis. Additionally, the weights for the feature clusters can be
used to extract meaningful information describing the individual sample
clusters identified, thereby, leading to comprehensive insights on individual
subtypes which allows to generate hypotheses that can be specifically tested.
The initial screening for overlap between the found information and biological
literature shows very promising results.

In our current approach, the feature clusters influence the patient clusters
but not vice versa. As shown on biological examples, these feature groups
already provide valuable information. However, to obtain feature clusters,
that are even more specific to the patient clusters, simultaneous optimization
of feature and patient clusters could be beneficial. Also are the feature
clusters currently restricted to be non-overlapping and to include every
feature. In many real-world scenarios we might want to a more flexible
approach that is relaxing these constraints.

## References

  * Assenov et al. (2014) Assenov, Y., Müller, F., Lutsik, P., Walter, J., Lengauer, T., and Bock, C. (2014).  Comprehensive analysis of DNA methylation data with RnBeads.  Nature Methods, 11(11), 1138–1140.
  * Bezdek (1981) Bezdek, J. C. (1981).  Pattern Recognition with Fuzzy Objective Function Algorithms.  Plenum Press.
  * Boumahdi et al. (2014) Boumahdi, S., Driessens, G., Lapouge, G., Rorive, S., Nassar, D., Le Mercier, M., Delatte, B., Caauwe, A., Lenglez, S., Nkusi, E., Brohée, S., Salmon, I., Dubois, C., del Marmol, V., Fuks, F., Beck, B., and Blanpain, C. (2014).  SOX2 controls tumour initiation and cancer stem-cell functions in squamous-cell carcinoma.  Nature, 511(7508), 246–250.
  * Cantley and Neel (1999) Cantley, L. C. and Neel, B. G. (1999).  New insights into tumor suppression: PTEN suppresses tumor formation by restraining the phosphoinositide 3-kinase/AKT pathway.  Proceedings of the National Academy of Sciences of the United States of America, 96(8), 4240–5.
  * Chandler et al. (2011) Chandler, K., Vance, C., Budnick, S., and Muller, S. (2011).  Muscle invasion in oral tongue squamous cell carcinoma as a predictor of nodal status and local recurrence: just as effective as depth of invasion?  Head and neck pathology, 5(4), 359–63.
  * Chen et al. (2018) Chen, Z., Zhao, H., Fu, N., and Chen, L. (2018).  The diversified function and potential therapy of ectopic olfactory receptors in non-olfactory tissues.  Journal of Cellular Physiology, 233(3), 2104–2115.
  * Cornelio et al. (2013) Cornelio, D. B., DE Farias, C. B., Prusch, D. S., Heinen, T. E., Dos Santos, R. P., Abujamra, A. L., Schwartsmann, G., and Roesler, R. (2013).  Influence of GRPR and BDNF/TrkB signaling on the viability of breast and gynecologic cancer cells.  Molecular and clinical oncology, 1(1), 148–152.
  * Coussens and Werb (2002) Coussens, L. M. and Werb, Z. (2002).  Inflammation and cancer.  Nature, 420(6917), 860–7.
  * Dunn (1974) Dunn, J. C. (1974).  A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters.  Journal of Cybernetics, 3(3), 32–57.
  * Gärtner et al. (2002) Gärtner, T., Flach, P. A., Kowalczyk, A., and Smola, A. J. (2002).  Multi-instance kernels.  In Proc. 19th International Conf. on Machine Learning, pages 179–186. Morgan Kaufmann.
  * Goldman et al. (2017) Goldman, M., Craft, B., Zhu, J., and Haussler, D. (2017).  The UCSC Xena system for cancer genomics data visualization and interpretation [abstract].  Proceedings of the American Association for Cancer Research Annual Meeting 2017; Cancer Research, 77(13 Supplement), 2584–2584.  Available from http://xena.ucsc.edu [Accessed: 2017/10/24].
  * Gönen and Margolin (2014) Gönen, M. and Margolin, A. A. (2014).  Localized Data Fusion for Kernel k-Means Clustering with Application to Cancer Biology.  In Advances in Neural Information Processing Systems 27 (NIPS 2014), pages 1305–1313.
  * Hartigan (1972) Hartigan, J. A. (1972).  Direct Clustering of a Data Matrix.  Journal of the American Statistical Association, 67(337), 123–129.
  * Hartigan (1975) Hartigan, J. A. (1975).  Clustering Algorithms.  John Wiley & Sons.
  * Hosmer Jr. et al. (2011) Hosmer Jr., D. W., Lemeshow, S., and May, S. (2011).  Applied Survival Analysis: Regression Modeling of Time to Event Data.  Wiley.
  * Hsin-Chien Huang et al. (2012) Hsin-Chien Huang, Yung-Yu Chuang, and Chu-Song Chen (2012).  Multiple Kernel Fuzzy Clustering.  IEEE Transactions on Fuzzy Systems, 20(1), 120–134.
  * Kumar et al. (2011) Kumar, A., Rai, P., and Daume, H. (2011).  Co-regularized multi-view spectral clustering.  In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 1413–1421. Curran Associates, Inc.
  * Le Van et al. (2016) Le Van, T., van Leeuwen, M., Carolina Fierro, A., De Maeyer, D., Van den Eynden, J., Verbeke, L., De Raedt, L., Marchal, K., and Nijssen, S. (2016).  Simultaneous discovery of cancer subtypes and subtype features by molecular data integration.  Bioinformatics, 32(17), i445–i454.
  * Rahimi and Gönen (2018) Rahimi, A. and Gönen, M. (2018).  Discriminating early- and late-stage cancers using multiple kernel learning on gene sets.  Bioinformatics, 34(13), i412–i421.
  * Rappoport and Shamir (2018) Rappoport, N. and Shamir, R. (2018).  Multi-omic and multi-view clustering algorithms: review and cancer benchmark.  Nucleic Acids Research.
  * Rodríguez-Rodero et al. (2017) Rodríguez-Rodero, S., Delgado-Álvarez, E., Díaz-Naya, L., Martín Nieto, A., and Menéndez Torre, E. (2017).  Epigenetic modulators of thyroid cancer.  Endocrinología, Diabetes y Nutrición, 64(1), 44–56.
  * Schlicker et al. (2006) Schlicker, A., Domingues, F. S., Rahnenführer, J., and Lengauer, T. (2006).  A new measure for functional similarity of gene products based on Gene Ontology.  BMC Bioinformatics, 7(1), 302.
  * Shen et al. (2012) Shen, R., Mo, Q., Schultz, N., Seshan, V. E., Olshen, A. B., Huse, J., Ladanyi, M., and Sander, C. (2012).  Integrative subtype discovery in glioblastoma using iCluster.  PloS ONE, 7.
  * Sinnott and Cai (2018) Sinnott, J. A. and Cai, T. (2018).  Pathway aggregation for survival prediction via multiple kernel learning.  Statistics in Medicine.
  * Song and Grandis (2000) Song, J. I. and Grandis, J. R. (2000).  STAT signaling in head and neck cancer.  Oncogene, 19(21), 2489–2495.
  * Speicher and Pfeifer (2015) Speicher, N. K. and Pfeifer, N. (2015).  Integrating different data types by regularized unsupervised multiple kernel learning with application to cancer subtype discovery.  Bioinformatics, 31(12), i268–i275.
  * Stöckel et al. (2016) Stöckel, D., Kehl, T., Trampert, P., Schneider, L., Backes, C., Ludwig, N., Gerasch, A., Kaufmann, M., Gessler, M., Graf, N., Meese, E., Keller, A., and Lenhof, H.-P. (2016).  Multi-omics enrichment analysis using the GeneTrail2 web service.  Bioinformatics, 32(10), 1502–1508.
  * The Cancer Genome Atlas (2012) The Cancer Genome Atlas (2012).  Comprehensive molecular portraits of human breast tumors.  Nature, 490, 61–70.
  * von Luxburg (2007) von Luxburg, U. (2007).  A tutorial on spectral clustering.  Statistics and Computing, 17(4), 395–416.
  * Wang et al. (2014) Wang, B., Mezlini, A. M., Demir, F., Fiume, M., Tu, Z., Brudno, M., Haibe-Kains, B., and Goldenberg, A. (2014).  Similarity network fusion for aggregating data types on a genomic scale.  Nature Methods, 11(3), 333–337.
  * Yu et al. (2012) Yu, S., Tranchevent, L.-C., Liu, X., Glanzel, W., Suykens, J. A. K., De Moor, B., and Moreau, Y. (2012).  Optimized Data Fusion for Kernel k-Means Clustering.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(5), 1031–1039.


# Data-Driven Learning of $3$-Point Correlation Functions as Microstructure
Representations

Sheng Cheng  Yang Jiao  Yi Ren  Computer Science, Arizona State University,
Tempe AZ 85287, United States  Materials Science and Engineering, Arizona
State University, Tempe AZ 85287, United States  Aerospace and Mechanical
Engineering, Arizona State University, Tempe AZ 85287, United States
<EMAIL_ADDRESS>
###### Abstract

This paper considers the open challenge of identifying complete, concise, and
explainable quantitative microstructure representations for disordered
heterogeneous material systems. Completeness and conciseness have been
achieved through existing data-driven methods, e.g., deep generative models,
which, however, do not provide mathematically explainable latent
representations. This study investigates representations composed of three-
point correlation functions, which are a special type of spatial convolutions.
We show that a variety of microstructures can be characterized by a concise
subset of three-point correlations, and the identification of such subsets can
be achieved by Bayesian optimization. Lastly, we show that the proposed
representation can directly be used to compute material properties based on
the effective medium theory.

###### keywords:

Quantitative microstructure representation, higher-order spatial correlations,
heterogeneous material reconstruction , Bayesian optimization

††journal: Acta Materialia

##  1 Introduction

This paper considers the open challenge of identifying quantitative
microstructure representations for disordered heterogeneous material systems.
Microstructures of such a material system can be captured as images drawn from
a 2D or 3D random field. We consider representations as the encoding of random
fields in a finite-dimensional Euclidean space, and ideally, this encoding
should be complete and concise. To explain, a representation is complete if
there exists a decoder, i.e., a mapping from the representation space to the
space of random fields, such that the decoded random field matches with the
sample distribution of authentic microstructures used to compute the
representation. And the representation is concise if it has minimal dimensions
within all complete representations. Both attributes are desired in a
materials design context: Completeness ensures that the representation can be
used for property prediction as it captures the morphological features of the
material system. Conciseness enables efficient search in a low-dimensional
representation space for new microstructures that potentially achieve
properties beyond existing limits. Nonetheless, these two attributes do not
guarantee explainability of predictive models built upon the representation,
e.g., for structure-property prediction. And the lack of explainability
hampers knowledge extraction, i.e., the natural-language understanding of how
properties of a material system are affected by its morphological features.
This is what differentiates our study from existing data-driven attempts at
learning microstructure representations: We seek representations that are
mathematically explainable, in addition to be complete and concise.

Figure 1: An example of successful modeling and reconstruction of a ceramic
matrix composite material [(a), original and (b), reconstruction] using the
two-point correlation function $S_{2}$ and an example where $S_{2}$ is not
sufficient to capture to key features of particle-reinforced composite [(c),
original and (d), reconstruction].

To provide a background, a large body of work on representation learning
exists at the intersection of materials science and machine learning. Most of
these propose either complete but non-explainable, or explainable but
incomplete representations. The former include purely data-driven generative
models, e.g., restricted Boltzmann machines [1], variational autoencoders
(VAE) [2], and generative adversarial networks (GAN) [3, 4], where a concise
and near-complete representation is learned through microstructure samples,
yet the encoders of which are composed of general-purpose neural networks and
are non-explainable. The alternative approach is to use statistical models
with more explainable, as opposed to learnable, architectures. These include
Gaussian random fields [5], geometric descriptors (e.g., grain/particle size
and shape distribution) [6, 7, 8, 9], spectral density functions [10, 11, 12],
and $n$-point correlation functions [13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28], among others. The encoding process of these methods
(e.g., statistical inference and geometric characterization) are manually
defined and can be explained. Yet due to the manual definitions, these
representations often have limited degrees of freedom to be learned to
approach completeness. To elaborate, Gaussian random fields have been used to
characterize interpenetrating bi-phase morphologies (e.g., rising from the
Spinodal decomposition process). This model assumes that the fluctuations in
the system can be completely specified by second order statistics, which does
not always hold for heterogeneous materials such as particle-reinforced
composites. Geometric descriptors are typically material-specific, which are
derived for microstructures composed of well-defined geometrical objects
(e.g., particles/grains). The spectral density function is the Fourier
transform of the two-point correlation function in the wave-number space,
which allows one to access the large-scale correlations (via small wave-number
analysis), yet have difficulty at accurately encoding small-scale
morphological features that are crucial for structure-property prediction.
Lastly, the $n$-point correlation functions encode the occurrence
probabilities of specific $n$-point configurations in the microstructure [29].
The set of correlation functions up to infinite orders fully characterizes a
random field [29, 30], and is therefore asymptotically complete. While it is
empirically shown that some material systems can be represented by concise
sets of correlation functions, e.g., metallic alloys, ceramic matrix
composites, and certain porous systems [31, 32, 33, 34, 35, 36] (Fig. 1a and
b), there is currently a lack of systematic tools for choosing a concise and
nearly complete set of correlations for any particular material system.

Within this context, our study investigates the feasibility and efficacy of
computational tools for identifying material-specific correlation functions as
complete, concise, and explainable representations of disordered heterogeneous
material systems.

To elaborate on the challenge, the number of correlation functions grows
exponentially as $d^{n-1}$ along the correlation order $n$ and the
dimensionality $d$ of the microstructure sample (usually 2 or 3). This causes
limited applications of higher-order correlation functions ($n\geq 3$) to
material systems [37, 38]. On the other hand, high-order correlation functions
are often required to represent disordered heterogeneous material systems [39,
40, 41] (Fig. 1c,d). So far the choice of representative correlations relies
on human intuition [42, 43].

As a first step towards automated selection of correlations, this paper
empirically verifies the following hypothesis: A concise and nearly-complete
representation composed of 3-point correlations can be learned from a dataset
of microstructure samples, for which 2-point correlations do not compose a
complete representation. To the authors’ best knowledge, this study is the
first to demonstrate the feasibility of representation learning using 3-point
correlations. The following technical contributions lay the foundation for
computationally tractable learning of higher-order correlations:

  * 1.

(i) We show that the encoder, i.e., the computation of $n$-point correlations
from microstructure samples, can be modeled as a special convolutional neural
network, allowing leverage of GPU computing and an existing deep learning
framework.

  * 2.

(ii) We show that the decoder, i.e., the reconstruction of microstructures,
can be formulated as a topology optimization problem, allowing leverage of the
state-of-the-art gradient-based algorithm which drastically improves the
reconstruction efficiency from standard stochastic methods such as simulated
annealing.

  * 3.

(iii) Building on top of (i) and (ii), we show that Bayesian Optimization can
effectively identify a concise subset of correlations as a nearly-complete
representation.

  * 4.

(iv) We demonstrate that the representations resulted from (iii) can be
directly utilized in combination of the effective medium theory to make
accurate and explainable predictions of physical properties of heterogeneous
material systems, which further supports the claim that they are nearly
complete.

The rest of the paper is organized as follows: Sec. 2 describes the models and
algorithms for learning correlation functions and demonstrates using a toy
case. In Sec. 3, the proposed method is applied to five material systems:
porous materials (i.e., sandstone), metal-ceramic composites, metallic alloys,
concrete microstructure, and particle-reinforced composites, which possess a
variety of distinct structural features, including both well separated compact
domains with different degrees of size and shape dispersion and
interpenetrating phase morphology with varying degrees of clustering. In Sec.
4, we show that the learned representations can provide reasonable estimates
of the physical properties of the material systems via the effective medium
theory. In Sec. 5, we provide concluding remarks and future directions.

##  2 Methods

To overview, our method includes three components: (i) An encoder that
computes $n$-point correlations via image convolution; (ii) a decoder that
performs gradient-based reconstruction of microstructures; and (iii) a search
algorithm for selecting a small and nearly-complete set of $n$-point
correlations. In the subsequent discussion, we describe each component in
details, which are further illustrated through a toy case of learning 3-point
correlations to represent an isosceles triangle.

###  2.1 Notations and preliminaries

Notations: The following notations will be used throughout the paper: The
indicator function $\mathbbm{1}(\textbf{c})$ returns 1 for all elements of the
Boolean vector c that are true and 0 otherwise. $\delta(x,a)=1$ if $x=a$ or
otherwise 0. 1 is a vector with all components equal to 1. $||\cdot||_{p}$ is
the $l_{p}$ norm, and $|\cdot|$ measures the size of a countable set.
$\textbf{C}(\cdot,\cdot;l,l^{\prime}):\mathcal{A}\times\mathcal{B}\rightarrow\mathbb{R}^{l\times
l}$ denotes the image convolution that discretizes the input $a\in\mathcal{A}$
and the filter $b\in\mathcal{B}$ as $l$-by-$l$ and
$l^{\prime}$-by-$l^{\prime}$ matrices, respectively.

For ease of exposition, we will focus on 2D bi-phase microstructures, although
the presented method can be extended to 3D multi-phase microstructures. A
microstructure sample is an indicator function
$y:\mathcal{X}\rightarrow\\{0,1\\}$ where $\mathcal{X}=[0,L]^{d}$ is the
support and $L$ is its linear size. We denote by $\mathcal{Y}$ the space of
$y$, and $\tilde{\mathcal{Y}}=\\{y_{i}\\}_{i=1}^{N}\subset\mathcal{Y}$ a
finite dataset that contains $N$ microstructure samples of the same material
system. The empirical data distribution is denoted by
$p_{\tilde{\mathcal{Y}}}$.

Figure 2: Visualization of the 3-point correlation set
$\mathcal{C}_{3}=\\{c_{1},c_{2},...\\}$, the polytopes
$\\{p_{1},p_{2},...p_{9}\\}$ of $c_{1}$, and the configuration-specific
vertices $\\{{\bf x}_{1}^{c_{1}},...,{\bf x}_{3}^{c_{1}}\\}$

$n$-point correlation functions: Given $n$, an $n$-point configuration space
$\mathcal{C}_{n}\in\mathbb{R}^{m}$ can be constructed, where $m=d^{n-1}$ is
the degrees of freedom of the $n$-point configuration. To explain using 2D
examples ($d=2$): When $n=2$, $\mathcal{C}_{2}=[0,\infty)\times[0,\pi)$
contains all distance-angle pairs that define line segments in
$\mathbb{R}^{2}$; when $n=3$, $\mathcal{C}_{3}=[0,\infty)^{3}\times[0,\pi)$
contains all distances-angle tuples that define triangles. Note that
$\mathcal{C}_{n}$ is a subspace of $\mathcal{C}_{n+1}$. Therefore
$\mathcal{C}_{\infty}$ contains all correlation configurations. Each
correlation configuration $c\in\mathcal{C}_{\infty}$ defines a set of
polytopes $\mathcal{P}_{c}$ and each polytope $p\in\mathcal{P}_{c}$ is defined
by a set of vertices $\\{{\bf x}_{1}^{p},...,{\bf x}_{n}^{p}\\}$ translated
from a configuration-specific set $\\{{\bf x}_{1}^{c},...,{\bf x}_{n}^{c}\\}$
(Fig. 2). Given $y$, an $n$-point correlation
$s:\mathcal{C}_{\infty}\times\mathcal{Y}\rightarrow[0,1]$ for the phase of
interest (assumed to be $y(\cdot)=1$) is defined as

$s(c,y)=\frac{1}{|\mathcal{P}_{c}|}\sum_{\\{{\bf x}_{1}^{p},...,{\bf
x}_{n}^{p}\\}\in\mathcal{P}_{c}}\mathbbm{1}\left(y({\bf x}_{1}^{p})=y({\bf
x}_{2}^{p})=...=y({\bf x}_{n}^{p})=1\right),$ (1)

which is the probability that a uniformly selected polytope
$p\in\mathcal{P}_{c}$ has vertices all belonging to the phase of interest
[30]. The order of $s$ is determined by $c$.

###  2.2 Computing $n$-point correlations through image convolution

We now show that $n$-point correlations can be computed through image
convolution, and therefore a standard convolutional neural network
architecture can be adopted. Let convolution kernels parameterized by an
$n$-point configuration $c$ be a function
$w_{c}:\mathcal{X}\rightarrow\\{0,1\\}$:

$w_{c}(x)=\sum_{x_{i}\in\\{{\bf x}_{1}^{c},...,{\bf
x}_{n}^{c}\\}}\delta(x,x_{i}).$ (2)

We can now compute the correlation as

$s(c,y)=\frac{\textbf{1}^{T}\mathbbm{1}\left(\textbf{C}(y,w_{c};l,l^{\prime})\geq
n\right)\textbf{1}}{l^{2}}.$ (3)

Since the indicator function is not differentiable, we propose the following
approximation:

$\hat{s}(c,y)=\frac{\textbf{1}^{T}\left(\sigma(\textbf{C}(y,w_{c};l,l^{\prime})-n+1)\right)\textbf{1}}{l^{2}},$
(4)

where $\sigma(x)=x$ if $x>0$ and 0 otherwise, and has differentiable
implementations [44]. Note that for bi-phase microstructures, the
approximation is exact. Lastly, we denote by
$\hat{s}(c,\tilde{\mathcal{Y}})=\mathbb{E}_{Y\sim
p_{\tilde{\mathcal{Y}}}}\left[\hat{s}(c,Y)\right]$ the mean correlation of the
dataset with respect to $c$.

###  2.3 Gradient-based microstructure reconstruction

Given a set of correlations $\hat{\mathcal{C}}\in\mathcal{C}_{\infty}$, a
dataset $\tilde{\mathcal{Y}}$, and target correlations
$\\{\hat{s}(c,\tilde{\mathcal{Y}})\\}_{c\in\hat{\mathcal{C}}}$, the
reconstruction problem is to find a microstructure $y$ with correlations close
to the target. With discretization of $\mathcal{X}$ and a focus on bi-phase
microstructures, $y$ can be considered as a $l$-by-$l$ binary matrix. This
leads to the following topology optimization problem:

$\min_{y\in\\{0,1\\}^{l^{2}}}\quad
L_{\hat{\mathcal{C}}}(y):=\sum_{c\in\hat{\mathcal{C}}}\left(\hat{s}(c,y)-\hat{s}(c,\tilde{\mathcal{Y}})\right)^{2}.$
(5)

We solve Eq. (5) using gradient descent with the following techniques.
Relaxation to a continuous problem: Note that Eq. (5) is combinatorial when
$y$ is binary. Therefore it is standard to solve a relaxed problem by
considering $y$ as a soft step function of another variable $u$:

$y(u;\beta)=\frac{\tanh(\beta/2)+\tanh(\beta(u-1/2))}{2\tanh(\beta/2)},$ (6)

where $u\in\mathbb{R}^{l\times l}$. We can then perform gradient-based search
using the gradient
$\nabla_{u}L_{\hat{\mathcal{C}}}(y(u;\beta))=\nabla_{y}L_{\hat{\mathcal{C}}}(y)\nabla_{u}y(u;\beta)$.
Avoiding gradient vanishing: The latter part of the gradient,
$\nabla_{u}y(u;\beta)$, vanishes quickly as $u$ moves away from 0 and this
issue exacerbates with larger $\beta$. On the other hand, large $\beta$ is
necessary for $y$ to approach a binary solution. A plausible solution to this
dilemma is to gradually increase $\beta$ during the search [45]. Our
experiments suggest that this scheme for changing $\beta$ is data dependent
(see Sec. 3). Removing artifacts in reconstructions: When solving topology
optimization problems, gradient descent can lead to local and artificial
checkerboard patterns in the solution. To address this issue, we apply a
Gaussian filter $w_{g}$ to the gradient at every search step, where $w_{g}$ is
a $3$-by-$3$ Gaussian kernel with a mean of 0 and a standard deviation of 3.
The standard deviation is set based on experiments.

###  2.4 A goodness metric for choosing correlations

Let $p_{\hat{\mathcal{C}}}$ be the distribution of reconstructed
microstructures resulted from Eq. (5) given a correlation set
$\hat{\mathcal{C}}$. The randomness is due to the random initialization of
$u$. Here we discuss the procedure for choosing a concise set
$\hat{\mathcal{C}}$ so that $p_{\hat{\mathcal{C}}}$ is representative of
$p_{\mathcal{Y}}$. To start, we propose the following problem:

$\min_{\hat{\mathcal{C}}\subset\mathcal{C}}\quad
L(\hat{\mathcal{C}})=\phi(p_{\hat{\mathcal{C}}},p_{\tilde{\mathcal{Y}}}),$ (7)

where $\phi(\cdot,\cdot)$ measures a distance between two distributions. We
note that morphological consistency between $p_{\hat{\mathcal{C}}}$ and
$p_{\tilde{\mathcal{Y}}}$ is of higher priority than matching between the
distributions, i.e., we consider $\hat{\mathcal{C}}$ a good representation as
long as samples from $p_{\hat{\mathcal{C}}}$ are morphologically similar to
$\tilde{\mathcal{Y}}$, even when the support of $p_{\hat{\mathcal{C}}}$ does
not fully overlap with that of $p_{\tilde{\mathcal{Y}}}$. With this insight,
we choose to define $\phi$ on the full correlation set $\mathcal{C}_{\infty}$:

$\phi(p_{1},p_{2})=\mathbb{E}_{y_{1}\sim p_{1},y_{2}\sim
p_{2}}\left[\frac{1}{|\mathcal{C}_{\infty}|}\sum_{c\in\mathcal{C}_{\infty}}(\hat{s}(c,y_{1})-\hat{s}(c,y_{2}))^{2}\right].$
(8)

However, Eq. (8) is ill-defined since $\mathcal{C}_{\infty}$ is an infinite-
dimensional space, and approximating $\mathcal{C}_{\infty}$ with large $n$ is
intractable. In fact, if we discretize each dimension of $\mathcal{C}_{n}$ by
$m$ levels, there will be $m^{d^{(n-1)}}$ correlation values to compute for
every $y$. As practical solutions, in this paper we will introduce task-
dependent and computationally tractable surrogates (denoted by $\hat{\phi}$))
as our best effort to approximating $\mathcal{C}_{\infty}$ (see Sec. 3 for
details).

###  2.5 Search of representations

Since computing $p_{\hat{\mathcal{C}}}$ relies on the reconstruction, which is
challenging to be differentiated, we propose to solve Eq. (7) using grid
search and Bayesian optimization (BO). For both algorithms, we fix the subset
size $|\hat{\mathcal{C}}|$. Grid search enumerates over all combinations of
correlations using a discretized $\mathcal{C}_{n}$; BO actively samples such
combinations without full enumeration. Both algorithms suffer from the curse
of dimensionality. To alleviate this challenge in the material cases in Sec.
3, we define $\hat{\mathcal{C}}$ as drawn from a distribution on
$\mathcal{C}_{\infty}$, and search only for the low-dimensional parameters of
this distribution.

###  2.6 A toy case

Here we use a toy case to walk through the algorithmic settings for the
encoding, decoding, and Bayesian optimization. It should be noted that same
details, including the parameterization of $\hat{\mathcal{C}}$ and the choice
of $\hat{\phi}$ are case-dependent. Case setup: Let the dataset be 2D binary
images, where each image $y\in\\{0,1\\}^{12\times 12}$ contains a filled
isosceles right triangle. The goal is to find a subset of correlations that
reconstruct the triangles. Since correlations are invariant to linear
translation, it is sufficient to keep a single image as the dataset (Fig. 3a).
While primitive, this case demonstrates the need for higher-order correlations
for reconstruction: Fig. 3b shows the inferior reconstruction result using the
full set of 2-point correlations. Reconstruction: We initialize a
deterministic reconstruction process by setting the initial guess to
$y(\cdot)=0$, and increase $\beta$ from 7 to 15 with an interval of 2. $\beta$
is increased for every 2000 iterations of gradient descent on the relaxed
version of Eq. (5). The step size for gradient descent is set to 2000 for
$\beta\leq 11$ and 1000 for larger $\beta$. Parameterization of the
correlation space: For the toy case, we choose to use a subset of three
3-point correlations for this toy study:
$\hat{\mathcal{C}}=\\{c_{1},c_{2},c_{3}\\}$, with the following
parameterization. First, we set $c_{1}$ as the volume fraction (i.e., with all
three vertices overlap). For $c_{2}$ and $c_{3}$, we fix one vertex and allow
the other two to move in the horizontal and vertical directions, respectively,
within a 18-by-18 box. This yields a 4-dimensional continuous and bounded
space where $\hat{\mathcal{C}}$ will be searched from. Criterion for choosing
correlations: For the toy case, we replace $\mathcal{C}_{\infty}$ with
$\mathcal{C}_{3}$ in Eq. (8), based on the empirical evidence to be presented
that 3-point correlations are sufficient for reconstructing the triangle.
Bayesian optimization: We start the search with 100 samples from the four
dimensional search space where the samples are drawn using a Latin Hypercube
sampler. Each sample defines a correlation set from which $\hat{\phi}$ can be
computed. With this data, a Gaussian process model can be built on the search
space, which allows prediction of $\hat{\phi}$ values and their uncertainties
for any correlation set. With this, we search for a new correlation set that
has the maximum expected improvement in $\hat{\phi}$ from the current samples.
Once sampled, this new data point is added to the existing dataset to update
the Gaussian process model. This process continues for 150 iterations.

Results and remarks: Fig. 3b-e summarizes the resultant $\hat{\mathcal{C}}$,
the corresponding reconstruction, and the BO loss history. Due to the
probabilistic nature of the algorithm, loss variances are estimated over five
independent BO trials. The results suggest that a subset of three 3-point
correlations is nearly complete for reconstructing the given sample. To show
that the resultant representation is concise, we compare the reconstruction
results from using two, three, and four 3-point correlations in Fig. 3d,f. The
correlations are chosen by BO in all cases. This toy case demonstrates that a
concise subset of 3-point correlations can be used to reconstruct images. In
the next section, we apply this idea to material microstructures.

Figure 3: Bayesian optimization for finding a concise representation of a
filled isosceles right triangle in (a) using 3-point correlations. (b)
Reconstruction fails with a full set of 2-point correlations. (c) Bayesian
optimized representations visualized as triangles (along with the volume
fraction at the origin) in the 3-point correlation landscape spanned by the
polytope parameters. (d) Resultant reconstruction based on (c). (e) BO
convergence. (f) Bayesian optimized reconstructions with two and four
correlations.

##  3 Experimental Results

We test five material systems, including sandstone, particle-reinforced
composites, concrete micro-structure, metal-ceramic composites and metallic
alloys. These material systems possess a variety of distinct structural
features, including both well separated compact domains with different degrees
of size and shape dispersion and interpenetrating phase morphology with
varying degrees of clustering. Our experiments suggest that a rich variety of
morphologically consistent microstructure samples can be reconstructed using
the proposed method, even when a single microstructure sample is used for
computing the representations. Therefore, all experiments will be based on a
single microstructure sample from each of the material systems. We will focus
on sandstone to elaborate on the algorithmic settings.

Figure 4: Sandstone sample

Case setup: A sandstone microstructure sample $y\in\\{0,1\\}^{128\times 128}$
is composed of the rock (shown in black) and the pore phase (shown in white).
Fig. 4 shows 2D slices of a 3D sandstone sample obtained via X-ray tomography
[46]. It can be seen that the sedimented rock grains possess a wide shape and
size dispersion, leading to complex morphology of the pore space. The material
possesses a statistically homogeneous and isotropic microstructure so that
3-point correlation information extracted from the 2D images are also
representative of the full 3D system [38].

Reconstruction: We initialize with $\beta=5$ and double its value up to 20
after a fixed amount of gradient descent steps. This amount is set to 50k for
$\beta\leq 13$ and 10k otherwise. The step size for each step is set to 2000.
The initial guess is randomized by sampling each pixel value from a Bernoulli
distribution parameterized by the volume fraction of the given microstructure.

Parameterization of the correlation space: For each correlation, we fix one
vertex at the center and allow the other two to move in a 41-by-41 box with a
stride of 2. To define a search space of correlations, we hypothesize that
both the number and the length scales of correlations affect the
reconstruction quality. Specifically, given a maximum length scale, more
correlations leads to better reconstruction, yet we expect that the benefit of
adding more correlations diminishes. On the other hand, for a fixed number of
correlations, there exists an optimal maximum length scale within which the
correlations are representative. With these hypotheses, we propose to search
on a 2D grid with the maximum correlation length scale as
$l\in\\{13,21,29,37\\}$ and the number of correlations as
$k\in\\{\frac{1}{64}K,\frac{1}{32}K,\frac{1}{16}K,\frac{1}{10}K,\frac{1}{8}K\\}$
where $K=7381$ is the total number of correlations when the maximum length
scale is 21. These parameter ranges are empirically chosen. Experimental
results suggest that larger $l$s and $k$s contribute minimally to
reconstruction improvements. Each grid point $(k,l)$ defines a distribution of
representations, each of which is composed of $k$ 3-point correlations which
are uniformly sampled from the set of correlations with a maximum length scale
of $l$.

Criterion for choosing correlations: From experimental results, we propose to
evaluate the reconstruction results using the pore-size distribution function
$P(r;y)$ as a surrogate. Specifically, $P(r;y)dr$ provides the probability
that a spherical “cavity” of radius $r$ can be entirely inserted into the
phase of interest centered at a randomly selected point in that phase of $y$
[30]. $P(r;y)$ was originally introduced to quantify the void phase in porous
materials and was subsequently generalized as a generic statistical descriptor
to quantify disordered heterogeneous materials. Since the porous phase is
hosting all transport processes such as fluid flow and chemical diffusion in
the porous media, the function $P(r;y)$ is shown to quantitatively related to
a variety of transport properties including effective diffusivity, fluid
permeability and mean-survival time of chemicals [30]. Moreover, $P(r;y)$
essentially provides a spherical measure of the clustering information in the
phase of interest and thus encodes information on higher-order correlations
[47]. In this study, we set $r\in\\{0,...,63\\}$ pixels and define

$\hat{\phi}(p_{1},p_{2})=\mathbb{E}_{y_{1}\sim p_{1},y_{2}\sim
p_{2}}\left[\frac{1}{64}\sum_{r\in\\{0,...,63\\}}\left|P(r;y_{1})-P(r;y_{2})\right|\right].$
(9)

Figure 5: (a) Reconstructions of sandstone with different numbers of
correlations $k$ and maximum correlation length scales. Sample means and
standard deviations of pore-size distribution function values are reported
based on 25 randomly drawn correlation sets from the corresponding $(k,l)$.
Unit: $\times 10^{-2}$. (b-d) The pore-size distribution, $l_{1}$ loss, and
weighted $l_{1}$ loss on the $(k,l)$ grid. (e) Convergence of Bayesian
optimization. (f, g) Reconstructions by random initial guesses and random
correlation sets drawn from the optimal $(k,l)$, respectively. (h) Correlation
landscape and the optimal maximum length scale $l$.

Bayesian Optimization: Since we have a relatively small candidate set, BO is
started with five configurations of $(k,l)$. For each configuration, we
compute $\hat{\phi}$ following Eq. 9 with five draws of correlations from the
uniform distribution specified by $(k,l)$. BO terminates with nine samples.
See Fig. 5e for BO convergence, where loss variances are estimated over five
independent BO trials. A full grid search is performed to verify the efficacy
of BO. See Fig. 5a for reconstruction results and pore-size values from
sampled correlations for all $(k,l)$ configurations and the same initial
guess. Standard deviations are estimated from five correlation sets drawn
independently based on the corresponding $(k,l)$. The Bayesian optimized
configuration is consistent with the grid search.

Remarks: (1) How much is the variance in reconstruction quality due to the
random initialization and the random sampling of correlations? We draw 25
initial guesses of $u$ from independent Bernoulli distributions parameterized
by the volume fraction of the given sandstone sample. The sample variance in
reconstruction quality is 3.31e-3. Fig. 5f shows five sample reconstructions.
Since the variance is small, we do not use multiple random initialization for
reconstruction during BO.

Fig. 5a reports variances in reconstruction quality for all $(k,l)$
configurations due to random sampling of correlations. Variances are estimated
based on 25 random draws of correlations for each grid point using the same
initialization of $u$. While the variances are not negligible, they are small
enough so that an optimal configuration ($k=738$ and $l=21$) can be identified
with statistical significance. We demonstrate reconstruction from five
randomly drawn correlation sets for the optimal configuration in Fig. 5g.

(2) Is pore-size distribution a good surrogate of $\phi$? Visual comparison in
Fig. 5a supports the use of pore-size distribution: The optimal configuration
achieves reconstruction quality no worse than others. We provide a second
evidence by comparing the pore-size results with those from alternatives that
are based on the full set of 3-point correlations with a maximum length scale
of 41. Minimal variations in correlations are observed beyond this length
scale. With minor abuse of notation, we will denote this set by
$\mathcal{C}_{3}$, and 3-point correlation vector of an input $y$ by
$c(y)=[c_{1}(y),...,c_{|\mathcal{C}_{3}|}(y)]$, where $c_{i}(y)$ is the
correlation of $y$ according to the $i$th element of $\mathcal{C}_{3}$. The
alternatives are based on the $l_{1}$ norm (Fig. 5c) and the weighted $l_{1}$
norm of $c(y)-c(y^{\prime})$, where $y$ and $y^{\prime}$ are the
microstructure sample and its reconstruction, respectively. We notice that the
sizes of subsets of $c_{i}(y)$ corresponding to different length scales
increases along the length scale, i.e., there are more correlations evaluated
on large triangles than small ones in $c(y)$. Since correlations
asymptotically approach zero as the length scale increases, $c(y)$ will
increasingly be dominated by zero elements when the maximum length scale
increases, making reconstruction qualities less distinguishable between two
choices of $k$ when using $l_{1}$ norm as a metric (See Fig. 5c). To address
this, we introduce a weighted norm as
$\sum_{i=1}^{|\mathcal{C}_{3}|}w_{i}|c_{i}(y)-c_{i}(y^{\prime})|$, where
$w_{i}$ is the proportion of the subset of correlations of the length scale
where $c_{i}$ belongs. This leads to an optimal representation consistent with
that from using the pore-size distribution function (Fig. 5d). However, the
choice of reweighting is not grounded in theory and can potentially be case
dependent.

(3) Is the chosen representation explainable? From the grid search results, we
conjecture that there is a critical length scale within which correlations are
representative, because when we choose larger $l$s, some of the randomly drawn
correlations will have larger length scales that are not representative, thus
lowering the reconstruction quality (see rows for $l=29$ and $37$ in Fig.
5a,b). To better understand this, we examine the 3-point correlation landscape
over a range of length scales, as shown in Fig. 5h. For visualization purpose,
we parameterize the three points of a polytope $p$ using
$r_{1}=||\textbf{x}_{1}^{p}-\textbf{x}_{2}^{p}||_{2}$,
$r_{2}=||\textbf{x}_{1}^{p}-\textbf{x}_{3}^{p}||_{2}$, and set the angle
$\theta=0\deg$ at vertex $x_{1}$. Each grid point in Fig. 5h represents the
average correlation value over all orientations of the polytope defined by
$r_{1}$ and $r_{2}$. We see that all variations in correlations occur within a
length scale, which corresponds well to the optimized one. The same
characteristic length scale can also been seen on correlation landscapes with
$\theta=45\deg$ and $\theta=90\deg$.

Figure 6: (a) Reconstructions of new materials using the optimal $(k,l)$ from
sandstone (b) Increasing the maximum length scale for concrete improves its
reconstruction quality (measured by the pore-size distribution function)

(4) Can the chosen representation be generalized to other material systems? We
now test whether the learned representation for sandstone can be applied to
other binary-phase material systems. From (3), we hypothesize that the
representation may fail on a new material system if its microstructure has a
different characteristic length scale from sandstone. Results on metal-ceramic
composite, Pb-Sn alloy, concrete, and particle-reinforced composite are
summarized in Fig. 6a. The reconstruction quality is acceptable for all
material systems except for concrete. We conjecture that this is because
metal-ceramic composite, Pb-Sn alloy, and particle-reinforced composite
possess similar morphology features as those for sandstone, e.g., distinct
local (short-range) correlations mainly determined by the shape and size of
local structural elements (e.g., particles) and their short-range ordering,
which then quickly decay to the average asymptotic values. Concrete, on the
other hand, has larger characteristic length scale than that of sandstone.
This hypothesis is empirically tested in Fig. 6b where we show that increasing
$l$ improves the reconstruction quality of concrete.

##  4 Property Predictions via Optimized Correlations

In this section, we show that the optimized representations can directly
provide reasonable estimates of the physical properties of the material
systems. This is achieved by employing the effective medium theory, in
particular, the so-called strong contrast expansion (SCE) formalism which
allows one to analytically express material properties such as elastic modulus
and thermal/electrical conductivity as an infinite series of integrals of
$n$-point correlation functions convoluted with the proper field kernels [48,
49], given the properties of individual material phases. It has been shown
that the SCE formalism possess superior convergence behavior such that
truncation of the infinite series at relatively lower order $n$ can already
yield very accurate estimates of the physical properties of interest for many
heterogeneous material systems [50].

Here, we focus on mechanical properties, including the overall (or effective)
bulk modulus $K_{e}$ and the shear modulus $G_{e}$, which are well defined for
all of the material systems investigated here. In particular, the 3-point SCE
approximation obtained by truncating the SCE series at $n=3$ for the effective
bulk modulus $K_{e}$ is given by

$\phi_{2}\frac{\kappa_{21}}{\kappa_{e1}}=1-\frac{(d+2)(d-1)G_{1}\kappa_{21}\mu_{21}}{d(K_{1}+2G_{1})}\phi_{1}\chi$
(10)

Similarly, the 3-point SCE approximation obtained by truncating the SCE series
at the $n=3$ for the effective shear modulus $G_{e}$ is given by

$\phi_{2}\frac{\mu_{21}}{\mu_{e1}}=1-\frac{2G_{1}\kappa_{21}\mu_{21}}{d(K_{1}+G_{1})}\phi_{1}\chi-\frac{2G_{1}\kappa_{21}\mu_{21}}{d(K_{1}+2G_{1})}\phi_{1}\chi-\frac{1}{2d}[\frac{dK_{1}+(d-2)G_{1}}{K_{1}+2G_{1}}]^{2}\mu_{21}^{2}\phi_{1}\eta_{2}$
(11)

where $\phi_{1}$ and $\phi_{2}$ are respectively the volume fraction of the
black (“matrix”) and the white (“inclusion”) phases; $d$ is the spatial
dimension of the material system; $K_{p}$ and $G_{p}$ are respectively the
bulk and shear modulus of phase $p$; the scalar parameters $\kappa_{pq}$ and
$\mu_{pq}$ ($p,q=1,2,e$) are respectively the bulk and shear modulus
polarizability, defined as

$\kappa_{pq}=\frac{K_{p}-K_{q}}{K_{q}+\frac{2(d-1)}{d}G_{q}}$ (12)

and

$\mu_{pq}=\frac{G_{p}-G_{q}}{G_{q}}\frac{1}{1+\frac{\frac{d}{2}K_{q}+\frac{(d+1)(d-2)}{d}G_{1}}{K_{q}+2G_{q}}};$
(13)

and $\chi$ and $\eta_{2}$ are the microstructural parameters associated with
the phases of interest incorporating the 3-pt correlations $s_{3}$ as well as
2-pt correlation correlations $s_{2}$ (which are certain subset of $s_{3}$),
i.e.,

$\chi=\frac{9}{2\phi_{1}\phi_{2}}\int_{0}^{\infty}\frac{dr}{r}\int_{0}^{\infty}\frac{ds}{s}\int_{-1}^{1}d(\cos\theta)P_{2}(\cos\theta)[s_{3}(r,s,t)-\frac{s_{2}(r)s_{2}(t)}{\phi_{2}}]$
(14)

and

$\eta_{2}=\frac{5}{21}\chi+\frac{150}{7\phi_{1}\phi_{2}}\int_{0}^{\infty}\frac{dr}{r}\int_{0}^{\infty}\frac{ds}{s}\int_{-1}^{1}d(\cos\theta)P_{4}(\cos\theta)[s_{3}(r,s,t)-\frac{s_{2}(r)s_{2}(t)}{\phi_{2}}],$
(15)

where $t=(r^{2}+s^{2}-2rs\cos\theta)^{1/2}$, and $P_{2}$ and $P_{4}$ are
respectively the Legendre polynomials of order two and four, i.e.,

$P_{2}(x)=\frac{1}{2}(3x^{2}-1),\quad
P_{4}(x)=\frac{1}{8}(35x^{4}-30x^{2}+3).$ (16)

Materials | $K_{1}$ | $G_{1}$ | $K_{2}$ | $G_{2}$
---|---|---|---|---
Sandstone (rock & void) | 35.8 | 10.2 | 0 | 0
Metal-ceramic Composite (SiC & Al) | 180 | 110 | 68 | 24
Pb-Sn Alloy | 42 | 14 | 45.8 | 4.9
Concrete (cement & rock) | 32 | 9 | 38 | 11.5
Particle-reinforced Composite (Al & SiC) | 68 | 24 | 180 | 110
Table 1: Individual phase properties ($K_{i}$ and $G_{i}$) for the heterogeneous material systems investigated in this work. The unit for the bulk and shear modulus are GPa. Materials | $K_{L}$ | $\bar{K}_{e}$ | $K_{e}$ | $K_{U}$ | $G_{L}$ | $\bar{G}_{e}$ | $G_{e}$ | $G_{U}$
---|---|---|---|---|---|---|---|---
Sandstone | 0 | 16.7$\pm$0.9 | 17.3 | 19.3 | 0 | 5.1$\pm$0.8 | 5.5 | 7.1
Metal-ceramic Composite | 132.2 | 134$\pm$3.8 | 136.9 | 141.0 | 68.4 | 70$\pm$1.7 | 74.3 | 78.8
Pb-Sn Alloy | 43 | 41.5$\pm$0.8 | 43.02 | 43.1 | 10.3 | 9.6$\pm$0.7 | 10.4 | 10.8
Concrete | 33.8 | 31.4$\pm$0.7 | 34.0 | 34.4 | 9.8 | 9.1$\pm$0.4 | 9.9 | 10.2
Particle-reinforced Composite | 86.9 | 87$\pm$2.9 | 89.2 | 96.3 | 37.9 | 38$\pm$1.5 | 38.8 | 44.7
Table 2: Comparison of the estimated elastic properties ($\bar{K}_{e}$ and
$\bar{G}_{e}$) for the heterogeneous material systems using the reduced set of
3-point representations to the corresponding Hashin-Shtrikman bounds ($K_{U}$,
$K_{L}$ and $G_{U}$, $G_{L}$), as well as the ground truth values ($K_{e}$ and
$G_{e}$) computed from SCE using the complete set of 3-point correlations. The
unit for the bulk and shear modulus are GPa.

The elastic modulus $K_{e}$ and $G_{e}$ for different material systems are
respectively estimated using Eq. (10)) and Eq. (11) by computing the 3-point
microstructural parameters $\chi$ and $\eta_{2}$ using the optimized
distribution of the 3-point representations obtained from the previous
section. The properties of individual phases for each material system are
given in Tab. 1. For each material system, 5 sets of 3-point statistics are
independently sampled based on the optimized distribution, and used to compute
the estimated properties using the corresponding SCE. The 2-point correlations
are computed by using the same set of distances defined by edge length of the
sampled triangular kernels. The averaged $K_{e}$ and $G_{e}$ are give in Tab.
2, which are compared to Hashin-Shtrikman bounds [51], as well as the “ground
truth” values, which are computed from SCE using the complete set of 3-point
functions. It can be seen from the table that the estimated elastic properties
using the reduced set of 3-point representations are within the theoretical
bounds and can provide reasonable estimates of the actual material properties.

We note again that purely data-driven models (e.g., representation via GAN or
VAE and representation-property mapping via supervised learning) can also be
used for property prediction [2]. However, such methods do not provide
explainable representations and require extra compute for supervised learning.

##  5 Conclusion and Discussions

We introduced a method to learn concise, complete, and explainable
representations of complex heterogeneous material systems based on 3-point
correlations. The key components of our method include a convolutional network
architecture for efficient computation of 3-point correlations, and an
algorithm for gradient-based microstructure reconstruction. Bayesian
optimization or grid search is subsequently applied to obtain the optimal
subset or distribution of 3-point correlations based on physics-inspired
reconstruction quality metrics. The utility of our method is demonstrated in
detail via the quantification, modeling and reconstruction of heterogeneous
material systems, each with distinct morphological features and degrees of
clustering. Moreover, we showed that the learned representations, in
combination with the effective medium theory, can be used to estimate elastic
properties of these material systems.

Although the example material systems used here are statistically homogeneous
and isotropic two-phase systems, we note that our general procedure can be
readily generalized to anisotropic, inhomogeneous and multi-phase material
systems. This is because neither the convolution-based encoder nor the
gradient-based decoder depends on the isotropy and homogeneity assumptions.
However, Bayesian optimization would be challenging as the kernel space for
the $n$-point correlations might not be easily characterized by a single
parameter distribution. To address this challenge, a possible solution is to
formulate the encoder-decoder process as a differential network, which allows
direct optimization on the kernel space. In our future work, we will explore
this approach, as well as constructing spatial-temporal correlations to
represent dynamics of evolving material systems.

## Data availability

The source codes and data in this work are available here.

## Acknowledgement

This work is supported by the National Science Foundation, Division of
Material Research under grant NO. 2020277 (AI Institute: Planning: Novel
Neural Architectures for 4D Materials Science)

## References

  * [1] R. Cang, Y. Xu, S. Chen, Y. Liu, Y. Jiao, M. Yi Ren, Microstructure representation and reconstruction of heterogeneous materials via deep belief network for computational material design, Journal of Mechanical Design 139 (7) (2017).
  * [2] R. Cang, H. Li, H. Yao, Y. Jiao, Y. Ren, Improving direct physical properties prediction of heterogeneous materials from imaging data via convolutional neural network and a morphology-aware generative model, Computational Materials Science 150 (2018) 212–221.
  * [3] Z. Yang, X. Li, L. Catherine Brinson, A. N. Choudhary, W. Chen, A. Agrawal, Microstructural materials design via deep adversarial learning methodology, Journal of Mechanical Design 140 (11) (2018).
  * [4] X. Li, Y. Zhang, H. Zhao, C. Burkhart, L. C. Brinson, W. Chen, A transfer learning approach for microstructure reconstruction and structure-property predictions, Scientific reports 8 (1) (2018) 1–13.
  * [5] A. P. Roberts, Statistical reconstruction of three-dimensional porous media from two-dimensional images, Physical Review E 56 (3) (1997) 3203.
  * [6] S. E. Wilding, D. T. Fullwood, Clustering metrics for two-phase composites, Computational materials science 50 (7) (2011) 2262–2272.
  * [7] P. G. Callahan, Quantitative characterization and comparison of precipitate and grain shape in nickel-base superalloys using moment invariants, Ph.D. thesis, Carnegie Mellon University (2012).
  * [8] M. Wang, J. Williams, L. Jiang, F. De Carlo, T. Jing, N. Chawla, Three dimensional (3d) microstructural characterization and quantitative analysis of solidified microstructures in magnesium-based alloys, Metallography, Microstructure, and Analysis 1 (1) (2012) 7–13.
  * [9] S. Ratanaphan, Y. Yoon, G. S. Rohrer, The five parameter grain boundary character distribution of polycrystalline silicon, Journal of materials science 49 (14) (2014) 4938–4945.
  * [10] A. Iyer, R. Dulal, Y. Zhang, U. F. Ghumman, T. Chien, G. Balasubramanian, W. Chen, Designing anisotropic microstructures with spectral density function, Computational Materials Science 179 (2020) 109559.
  * [11] U. Farooq Ghumman, A. Iyer, R. Dulal, J. Munshi, A. Wang, T. Chien, G. Balasubramanian, W. Chen, A spectral density function approach for active layer design of organic photovoltaic cells, Journal of Mechanical Design 140 (11) (2018).
  * [12] D. Chen, S. Torquato, Designing disordered hyperuniform two-phase materials with novel physical properties, Acta Materialia 142 (2018) 152–161.
  * [13] S. Niezgoda, D. Fullwood, S. Kalidindi, Delineation of the space of 2-point correlations in a composite material system, Acta Materialia 56 (18) (2008) 5285–5292.
  * [14] A. Cecen, T. Fast, S. R. Kalidindi, Versatile algorithms for the computation of 2-point spatial correlations in quantifying material structure, Integrating Materials and Manufacturing Innovation 5 (1) (2016) 1–15.
  * [15] A. Choudhury, Y. C. Yabansu, S. R. Kalidindi, A. Dennstedt, Quantification and classification of microstructures in ternary eutectic alloys using 2-point spatial correlations and principal component analyses, Acta Materialia 110 (2016) 131–141.
  * [16] H. Okabe, M. J. Blunt, Pore space reconstruction using multiple-point statistics, Journal of Petroleum Science and Engineering 46 (1-2) (2005) 121–137.
  * [17] D. T. Fullwood, S. R. Niezgoda, S. R. Kalidindi, Microstructure reconstructions from 2-point statistics using phase-recovery algorithms, Acta Materialia 56 (5) (2008) 942–948.
  * [18] Y. Jiao, F. Stillinger, S. Torquato, Modeling heterogeneous materials via two-point correlation functions: Basic principles, Physical Review E 76 (3) (2007) 031110.
  * [19] Y. Jiao, F. Stillinger, S. Torquato, Modeling heterogeneous materials via two-point correlation functions. ii. algorithmic details and applications, Physical Review E 77 (3) (2008) 031135.
  * [20] A. Hajizadeh, A. Safekordi, F. A. Farhadpour, A multiple-point statistics algorithm for 3d pore space reconstruction from 2d images, Advances in water Resources 34 (10) (2011) 1256–1267.
  * [21] P. Tahmasebi, M. Sahimi, Cross-correlation function for accurate reconstruction of heterogeneous media, Physical review letters 110 (7) (2013) 078002.
  * [22] P. Tahmasebi, A. Hezarkhani, M. Sahimi, Multiple-point geostatistical modeling based on the cross-correlation functions, Computational Geosciences 16 (3) (2012) 779–797.
  * [23] H. Xu, Y. Li, C. Brinson, W. Chen, A descriptor-based design methodology for developing heterogeneous microstructural materials system, Journal of Mechanical Design 136 (5) (2014).
  * [24] K. M. Gerke, M. V. Karsanina, Improving stochastic reconstructions by weighting correlation functions in an objective function, EPL (Europhysics Letters) 111 (5) (2015) 56002.
  * [25] M. V. Karsanina, K. M. Gerke, Hierarchical optimization: Fast and robust multiscale stochastic reconstructions with rescaled correlation functions, Physical review letters 121 (26) (2018) 265501.
  * [26] J. Feng, Q. Teng, X. He, X. Wu, Accelerating multi-point statistics reconstruction method for porous media via deep learning, Acta Materialia 159 (2018) 296–308.
  * [27] Y. Gao, Y. Jiao, Y. Liu, Efficient high-dimensional material reliability analysis with explicit voxel-level stochastic microstructure representation, Applied Mathematical Modelling 91 (2021) 1117–1140.
  * [28] Y. Gao, Y. Jiao, Y. Liu, Ultra-efficient reconstruction of 3d microstructure and distribution of properties of random heterogeneous materials containing multiple phases, Acta Materialia 204 (2021) 116526.
  * [29] S. Torquato, G. Stell, Microstructure of two-phase random media. i. the n-point probability functions, The Journal of Chemical Physics 77 (4) (1982) 2071–2077.
  * [30] S. Torquato, Random heterogeneous materials: microstructure and macroscopic properties, Appl. Mech. Rev. 55 (4) (2002) B62–B63.
  * [31] Y. Jiao, E. Padilla, N. Chawla, Modeling and predicting microstructure evolution in lead/tin alloy via correlation functions and stochastic material reconstruction, Acta Materialia 61 (9) (2013) 3370–3377.
  * [32] E.-Y. Guo, N. Chawla, T. Jing, S. Torquato, Y. Jiao, Accurate modeling and reconstruction of three-dimensional percolating filamentary microstructures from two-dimensional micrographs via dilation-erosion method, Materials Characterization 89 (2014) 33–42.
  * [33] Y. Jiao, N. Chawla, Modeling and characterizing anisotropic inclusion orientation in heterogeneous material via directional cluster functions and stochastic microstructure reconstruction, Journal of Applied Physics 115 (9) (2014) 093511.
  * [34] S. Chen, H. Li, Y. Jiao, Dynamic reconstruction of heterogeneous materials and microstructure evolution, Physical Review E 92 (2) (2015) 023301.
  * [35] S. Chen, A. Kirubanandham, N. Chawla, Y. Jiao, Stochastic multi-scale reconstruction of 3d microstructure consisting of polycrystalline grains and second-phase particles from 2d micrographs, Metallurgical and Materials Transactions A 47 (3) (2016) 1440–1450.
  * [36] Y. Xu, S. Chen, P.-E. Chen, W. Xu, Y. Jiao, Microstructure and mechanical properties of hyperuniform heterogeneous materials, Physical Review E 96 (4) (2017) 043301.
  * [37] S. R. Kalidindi, S. R. Niezgoda, A. A. Salem, Microstructure informatics using higher-order statistics and efficient data-mining protocols, Jom 63 (4) (2011) 34–41.
  * [38] H. Malmir, M. Sahimi, Y. Jiao, Higher-order correlation functions in disordered media: Computational algorithms and application to two-phase heterogeneous materials, Physical Review E 98 (6) (2018) 063317.
  * [39] C. J. Gommes, Y. Jiao, S. Torquato, Density of states for a specified correlation function and the energy landscape, Physical review letters 108 (8) (2012) 080601.
  * [40] C. J. Gommes, Y. Jiao, S. Torquato, Microstructural degeneracy associated with a two-point correlation function and its information content, Physical Review E 85 (5) (2012) 051140.
  * [41] Y. Jiao, F. H. Stillinger, S. Torquato, Geometrical ambiguity of pair statistics. ii. heterogeneous media, Physical Review E 82 (1) (2010) 011106.
  * [42] P.-E. Chen, W. Xu, N. Chawla, Y. Ren, Y. Jiao, Hierarchical n-point polytope functions for quantitative representation of complex heterogeneous materials and microstructural evolution, Acta Materialia 179 (2019) 317–327.
  * [43] P.-E. Chen, W. Xu, Y. Ren, Y. Jiao, Probing information content of hierarchical n-point polytope functions for quantifying and reconstructing disordered systems, Physical Review E 102 (1) (2020) 013305.
  * [44] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer, Automatic differentiation in pytorch (2017).
  * [45] J. Wu, N. Aage, R. Westermann, O. Sigmund, Infill optimization for additive manufacturing—approaching bone-like porous structures, IEEE transactions on visualization and computer graphics 24 (2) (2017) 1127–1140.
  * [46] H. Li, P.-E. Chen, Y. Jiao, Accurate reconstruction of porous materials via stochastic fusion of limited bimodal microstructural data, Transport in Porous Media 125 (1) (2018) 5–22.
  * [47] Y. Jiao, F. Stillinger, S. Torquato, A superior descriptor of random textures and its predictive capacity, Proceedings of the National Academy of Sciences 106 (42) (2009) 17634–17639.
  * [48] S. Torquato, Effective stiffness tensor of composite media—i. exact series expansions, Journal of the Mechanics and Physics of Solids 45 (9) (1997) 1421–1448.
  * [49] S. Torquato, Exact expression for the effective elastic tensor of disordered composites, Physical review letters 79 (4) (1997) 681.
  * [50] H. Li, S. Singh, C. S. Kaira, J. C. Mertens, J. J. Williams, N. Chawla, Y. Jiao, Microstructural quantification and property prediction using limited x-ray tomography data, JOM 68 (8) (2016) 2288–2295.
  * [51] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, Journal of the Mechanics and Physics of Solids 11 (2) (1963) 127–140.