# Non-convex SVM for cancer diagnosis based on morphologic features of tumor
microenvironment

Sean Kent0000-0001-8697-9069  Menggang
<EMAIL_ADDRESS>[ Department of
Statistics, University of Wisconsin – Madison  Department of Biostatistics and
Medical Informatics, University of Wisconsin – Madisonpresep=, ]e2

###### Abstract

The surroundings of a cancerous tumor impact how it grows and develops in
humans. New data from early breast cancer patients contains information on the
collagen fibers surrounding the tumorous tissue—offering hope of finding
additional biomarkers for diagnosis and prognosis—but poses two challenges for
typical analysis. Each image section contains information on hundreds of
fibers, and each tissue has multiple image sections contributing to a single
prediction of tumor vs. non-tumor. This nested relationship of fibers within
image spots within tissue samples requires a specialized analysis approach.

We devise a novel support vector machine (SVM)-based predictive algorithm for
this data structure. By treating the collection of fibers as a probability
distribution, we can measure similarities between the collections through a
flexible kernel approach. By assuming the relationship of tumor status between
image sections and tissue samples, the constructed SVM problem is non-convex
and traditional algorithms can not be applied. We propose two algorithms that
exchange computational accuracy and efficiency to manage data of all sizes.
The predictive performance of both algorithms is evaluated on the collagen
fiber data set and additional simulation scenarios. We offer reproducible
implementations of both algorithms of this approach in the R package
`mildsvm`.

Breast cancer,

Support Vector Machines,

Weakly supervised,

Functional data,

###### keywords:

and

##  1 Introduction

The tumor microenvironment is the ecosystem that surrounds and interacts with
a tumor. The composition of the tumor microenvironment includes the immune
cells, stromal cells, blood vessels, and extracellular matrix (ECM) (Anderson
and Simon, 2020). It is now widely accepted that the tumor microenvironment
plays an important role in supporting tumor growth and progression (Anderson
and Simon, 2020). Early in tumor growth, a dynamic and reciprocal relationship
develops between cancer cells and components of the tumor microenvironment
that supports cancer cell survival, local invasion, and metastatic
dissemination (Baghban et al., 2020). Therefore, investigation of the tumor
microenvironment can help discover noteworthy characteristics (e.g.
biomarkers), leading to better cancer diagnosis and prognosis. Expanding
literature on the tumor microenvironment has even identified new targets
within it for therapeutic intervention (Bejarano, Jordāo and Joyce, 2021).

In this article, we analyze morphological features of the ECM for breast
cancer diagnosis. Breast tumor formation is associated with a stromal
response, termed the desmoplastic reaction, characterized by amplified
collagen matrix deposition and stromal cell recruitment and activation,
thereby promoting tumor progression (Arendt et al., 2010; Zeltz et al., 2020).
Because increased cell numbers and increased collagen are sources of contrast
within the mammogram, they are difficult to distinguish; traditional
clinically proven methods, such as radiography and ultrasound imaging, do not
have the resolution to differentiate the tumor from collagen at the cellular
level. Recent imaging development, known as second harmonic generation (SHG)
imaging (Chen et al., 2012), can accurately capture over 20 collagen
features—including alignment, curvature, width, length, closeness to nearby
fibers, and other features of ECM composition. These morphological features
can be used as early biomarkers for cancer diagnosis and prognosis (Conklin et
al., 2018).

A recent data set with these morphological features contains a nested
structure that poses a few challenges. The data includes slides of breast
tissue from ductal carcinoma in situ (DCIS) patients that are either tumor
tissue or adjacent, normal tissue. For each tissue slide, between one and
eight spots (median 5) are selected for SHG imaging (Figure 1). Within each
spot, SHG imaging captures features on hundreds of collagen fibers. Figure 1
shows detailed counts of the number of fibers per spot across the data set. An
example of two subject slides, provided in Figure 2 for two informative
features, illustrates two challenges in creating a tumor prediction algorithm.
The first is how to summarize collagen relationships within each spot from the
hundreds of fibers. Figure 2 shows that the mean and standard deviation of the
distance to the eighth closest fiber may be predictors; for other features,
correlations and quantiles may also be relevant. The second challenge is how
to relate the summarized spot features to the prediction on the whole tissue
slide where our training labels occur since each has many spots with
information.

Figure 1: Number of spots and fibers in the data. A. Counts of the number of
spots selected per tissue, with 5 spots being the most common. B. Counts of
the number of collagen fibers per spot, with 450 being the most common. Figure
2: Example data from two subjects. Each small point represents the
measurements on individual fibers within the tumor microenvironment, while
large points (outlined in black) summarize the mean for each chosen spot. Both
axes are on the log scale. A. Color distinguishes subjects while shape
distinguishes whether the tissue sample was tumorous or not. We observe
subject 1 with some spots close to the spots from subject 2, and other spots
far from the cluster. B. Only subject 1 is shown. Color distinguishes the
individual spots for this subject, showing distributional differences in both
location and shape among the five spots chosen.

We propose to tackle these challenges by connecting approaches from previous
literature that address one of the two difficulties and constructing a non-
convex support vector machine (SVM) learning algorithm. First, we capture the
pooled information from the many fibers per spot by treating the collection as
a sample from a multivariate probability distribution. This distribution
summarizes the relationships between collagen fibers rather than looking
individually within them. We leverage growing literature on learning from
probability distributions in works such as Muandet et al. (2012, 2017); Lin et
al. (2013), which is part of the broader category of functional data analysis
(Ramsay, 2006). Second, we assume that for each tumor slide, at least one of
the spots has signs of cancer. Conversely, among non-tumor slides, no spots
have tumor features. This assumption creates weak supervision between the
labels at the spot level and the predictions at the tissue sample level
(Campanella et al., 2019; Li et al., 2013). Our approach combines both ideas
and minimizes a hinge loss function on the data while maximizing the margin of
error for classification, resulting in a non-convex type of SVM.

Our proposed methodology introduces two algorithmic solutions to the non-
convex SVM. One heuristically moves between selecting important spots and then
solving a reduced problem on the selected spots. This approach has roots in
several optimization approaches, including the Expectation-Maximization
algorithm (Moon, 1996), the Convex Concave Procedure (Yuille and Rangarajan,
2003), and a similar non-convex SVM for weakly supervised data (Andrews,
Tsochantaridis and Hofmann, 2003). The other approach solves the problem
directly by modifying its form and using efficient algorithms for problems
involving integer variables (Gurobi Optimization, 2021). In doing so, we
contribute new ideas for approximating non-linear kernels when learning from
data with probability distributions. Both methods are applied to the
motivating data on tumor microenvironment features, and they perform favorably
without requiring the researcher to select the meaningful collagen features a
priori. We also provide an R package to solve the non-convex SVM problem for
data with a similar structure.

The remainder of this paper contains the following: a detailed account of our
methodology—including a formalized data set structure, non-convex SVM
derivations, and two computational algorithms—in Section 2; predictive
performance on simulated data sets in Section 3; an empirical evaluation on
the motivating collagen feature data set in Section 4; and a discussion of
methodological considerations, implications, and future research directions in
Section 5.

##  2 Methodology

In this section, we formulate a predictive model from data structured with
weakly-supervised outcome labels and distributional instances—as seen in the
motivating collagen fiber data set. We will transform the joint goal of
minimizing training loss and maximizing the error margin into a non-convex SVM
optimization problem. We then develop two separate algorithms to solve the SVM
problem: one with a heuristic approach iteratively solving a convex sub-
problem and another directly solving the problem after a feature map
approximation for probability distributions.

###  2.1 Data set structure

Consider a data set with $n$ distributional instances. Let $\mathcal{P}$
denote the set of all probability measures $\mathbb{P}$ over the data domain
$\mathcal{X}\subseteq\mathbb{R}^{d}$ for $d\in\mathbb{N}$. Formally, we say
that each instance $i$ has information contained by a probability measure
$\mathbb{P}_{i}$ on the data domain $\mathcal{X}$. In practice, this
distribution could manifest in the data as a sample of size $r_{i}$ drawn from
$\mathbb{P}_{i}$, denoted as
$\hat{\mathbb{P}}_{i}=\\{x_{i,j}\\}_{j=1}^{r_{i}}$, where
$x_{i,j}\in\mathcal{X}$. Collectively, we can refer to this data set as
$\\{\mathbb{P}_{i}\\}_{i=1}^{n}$ or $\\{\hat{\mathbb{P}}_{i}\\}_{i=1}^{n}$. In
our motivating data set, each distributional instances is a spot with samples
from the individual collagen fibers.

Next, we consider weakly-supervised outcome labels within the data set. Let
$\mathcal{I}$ represent the set of all top-level groups in the data, with each
group $I\in\mathcal{I}$ having a corresponding label of interest given by
$Y_{I}\in\\{-1,1\\}$. Each top-level group is a collection of spots $i\in I$
and can thus be thought of as an index set.111The index sets $I\in\mathcal{I}$
together create a partition of the spots $\\{1,2,\dots,n\\}$ into mutually
exclusive groups that include all spots:
$\bigcup_{I\in\mathcal{I}}I=\\{1,\dots,n\\}$ and $I\cap J=\emptyset$ for all
pairs $I,J\in\mathcal{I}$. The data set information for each top-level group
is denoted by $B_{I}=\\{\mathbb{P}_{i}:i\in I\\}$. Each probability
distributions $\mathbb{P}_{i}$ within $I$ has a potentially unobserved binary
label $y_{i}\in\\{-1,1\\}$, and, collectively, they contribute to the top-
level label $Y_{I}$. This nested structure results in a weak supervision
between the spots and the outcome, and the term has been coined in previous
works (Campanella et al., 2019; Li et al., 2013). In our data, the top-level
groups correspond to tissue slides, and labels refer to tumor (+1) and non-
tumor (-1).

###  2.2 Non-convex SVM formulation

Our goal is to find a function $F$ such that $F(B_{I})\approx Y_{I}$ using the
collective information from $\\{\mathbb{P}_{i}:i\in I\\}$.222We use $\approx$
to indicate closeness with respect to some loss function First, we assume that
the spot level outcomes $y_{i}$ relate to the slide-level group outcomes
$Y_{I}$ as follows: a positive slide label $(Y_{I}=1)$ is observed if and only
if at least one spot label in it is positive $(\exists i\in I:y_{i}=1)$.
Mathematically, this is equivalent to the following assumption:

######  Assumption 2.1.

$Y_{I}=\max_{i\in I}y_{i}\quad\forall I\in\mathcal{I}$.

Assumption 2.1 works well for the motivating data set on collagen fibers using
tumor for the positive label. We assume that if a single spot on the tissue
sample is tumorous, the specimen is tumorous. Conversely, non-tumor samples
must have all spots with non-tumorous features. This assumption could be
violated if no tumorous areas were selected as spots in a tumor tissue sample.
However, the procedure for determining spots had experts focus on tumorous
areas, making this violation unlikely in practice.

Working at the spot level, we suppose that a linear form is sufficient to
classify labels after a possible non-linear transformation $\phi$, so that
$f(\mathbb{P}_{i})=\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b$ and
$f(\mathbb{P}_{i})\approx y_{i}$. For the moment, consider $\boldsymbol{w}$
and $\phi$ as finite dimensional. Then following Assumption 2.1, a natural
form for $F$ is $F(B_{I})=\max_{i\in I}f(\mathbb{P}_{i})$. The hyperplane
$F(B_{I})=0$ divides the slide-level groups into two regions which are
classified as positive and negative, and in a similar way
$f(\mathbb{P}_{i})=0$ divides the instances into positive and negative
regions.

Following an SVM approach (Burges, 1998) we seek to 1) minimize the hinge loss
function $\ell(y,x)=\max(0,1-yx)$ between labels $Y_{I}$ and predictions
$F(B_{I})$ for all $I$ and 2) maximize the margin around the hyperplane
created by $f(\mathbb{P}_{i})=0$, which is equivalent to minimizing
$\left\lVert\boldsymbol{w}\right\rVert^{2}$. Hinge loss is used because if
$\ell\left(Y_{I},\max_{i\in
I}\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)=0$, then
$Y_{I}\max_{i\in
I}\left(\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)\geq 1$,
i.e. group labels fall outside the margin of the classifier and are correctly
classified; otherwise, the hinge loss increases proportionally to the distance
from that margin. Larger margins are beneficial because they offer the hope of
strong performance on unseen data. We balance these two objectives with a
tuneable hyper-parameter $C$:

$\min_{\boldsymbol{w},b}\left\\{\frac{1}{2}\left\lVert\boldsymbol{w}\right\rVert^{2}+C\sum_{I\in\mathcal{I}}\ell\left(Y_{I},\max_{i\in
I}\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)\right\\}$

We can re-write this objective in a more common form by using auxiliary
variables $\boldsymbol{\xi}=\left[\xi_{I}\right]_{I\in\mathcal{I}}$ defined by
$\xi_{I}=\ell\left(Y_{I},\max_{i\in
I}\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)$. From the
definition of $\ell$, $\xi_{I}$ is the smallest non-negative number such that
$Y_{I}\max_{i\in
I}\left(\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)\geq
1-\xi_{I}$, and so we have the following optimization problem:

$\begin{split}&\min_{\boldsymbol{w},b,\boldsymbol{\xi}\geq
0}\frac{1}{2}\left\lVert\boldsymbol{w}\right\rVert^{2}+C\sum_{I\in\mathcal{I}}\xi_{I}\\\
\textrm{s.t. \quad}&Y_{I}\max_{i\in
I}\left(\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)\geq
1-\xi_{I},\enskip\forall I\in\mathcal{I}\end{split}$ (1)

where $\boldsymbol{\xi}\geq 0$ indicates that $\xi_{I}\geq 0\enskip\forall I$.
We refer to (1) as the primal problem. It is a non-convex mixed-integer
quadratic programming (MIQP) problem. The mixed-integer nature can be revealed
by explicitly considering which $i\in I$ yields the maximum value of
$\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b$ for each of the slide-
level groups in the constraint of (1).

Related research has studied non-convex SVMs under several different contexts.
One substantial branch added a non-convex penalty to the hinge loss to induce
sparsity for variable selection (Zhang et al., 2006; Laporte et al., 2013;
Zhang et al., 2016; Liu, Yao and Li, 2016; Guan et al., 2020). Another
introduced non-convexity through the ramp loss to reduce the influence of
outliers (Ertekin, Bottou and Giles, 2010; Zhao et al., 2021). Our approach in
Problem (1) differs because non-convexity results from the $\max$ operator in
the constraint, instead of from a change in the loss function. Andrews,
Tsochantaridis and Hofmann (2003) arrive at a similar formulation when
considering standard vector data instead of distributional instances. In both
the semi- and weakly-supervised settings, Li et al. (2013) explore a more
general non-convex, MIQP SVM problem through a convex relaxation approach.

###  2.3 A heuristic algorithm from dual SVM problem

To calculate predictive models from data structured with weakly-supervised
outcome labels and distributional instances, we need algorithms that can solve
the optimization problem in (1). Our first algorithm, which we will call mi-
smm (heuristic), relies on finding the dual to a convex sub-problem of (1) and
iteratively solving that problem using an appropriate kernel function on
distributional instances.

####  2.3.1 The convex sub-problem and dual SVM

While Problem (1) has an intuitive setup, it is often inconvenient to work
with. If $\boldsymbol{w}$ and $\phi$ are infinite-dimensional or the form of
$\phi$ is unknown, solving (1) is intractable. One way of resolving the
intractability is by using Wolfe duality theory to find an equivalent problem
on which the “kernel trick” can be applied (Boyd and Vandenberghe, 2004).
However, many of the guarantees within duality theory do not apply to non-
convex problems.

One solution is to find a convex sub-problem within (1). Consider a function
$s$ that selects one spot $i$ from each slide-level group $I$ having positive
label $Y_{I}=1$, so that $s(I)=i$. An optimal $s$ will choose
$s(I)=\arg\max_{i\in
I}\left(\langle\boldsymbol{w},\phi(\mathbb{P}_{i})\rangle+b\right)$. For any
fixed (potentially non-optimal) $s$, we can solve the convex SVM problem:

$\begin{split}&\min_{\boldsymbol{w},b,\boldsymbol{\xi}\geq
0}\frac{1}{2}\left\lVert\boldsymbol{w}\right\rVert^{2}+C\sum_{I\in\mathcal{I}}\xi_{I}\\\
\textrm{s.t.
\quad}&Y_{I}\left(\langle\boldsymbol{w},\phi(\mathbb{P}_{s(I)})\rangle+b\right)\geq
1-\xi_{I},\enskip\forall I\in\mathcal{I}\end{split}$ (2)

On the convex SVM problem (2), we can follow standard duality theory by
calculating a Lagrangian function and evaluating the Karush–Kuhn–Tucker
conditions (Boyd and Vandenberghe, 2004) to find the dual problem, which has a
similar form as that of Andrews, Tsochantaridis and Hofmann (2003):

$\begin{split}\max_{\boldsymbol{\alpha}}&\sum_{i\in
E}\alpha_{i}-\frac{1}{2}\sum_{i\in E}\sum_{j\in
E}\alpha_{i}\alpha_{j}Y_{B(i)}Y_{B(j)}\langle\phi(\mathbb{P}_{i}),\phi(\mathbb{P}_{j})\rangle\\\
\textrm{s.t. \quad}\forall I\in\mathcal{I}\text{ : }&Y_{I}=+1\text{ and
}0\leq\alpha_{i}\leq C,\enskip\forall i\in I,\\\ \text{or }&Y_{I}=-1\text{ and
}0\leq\sum_{i\in I}\alpha_{i}\leq C,\\\ &\sum_{i\in
E}\alpha_{i}Y_{B(i)}=0,\quad\text{and}\quad\alpha_{i}\geq 0,\enskip\forall
i\in E\end{split}$ (3)

where $E=\left(\bigcup_{I:Y_{I}=-1}\\{i\in
I\\}\right)\cup\left(\bigcup_{I:Y_{I}=+1}s(I)\right)$ is the effective set of
indexes consisting of all negative spots and the selected spot from each
positive group, and where $B(i)$ is the slide associated with spot $i$ (i.e.
$B(i)=I$ for $i\in I$). Problem (3) is a convex quadratic program that is
solvable by standard QP solvers or by chunking algorithms like the sequential
minimal optimization (SMO) algorithm (Platt, 1998; Boser, Guyon and Vapnik,
1992).

Similar to the standard SVM case, we construct the classifier based on
$\alpha_{i}$ as $f(\mathbb{P})=\sum_{i\in
E}\alpha_{i}Y_{B(i)}\langle\phi(\mathbb{P}),\phi(\mathbb{P}_{i})\rangle+b$.
Given the final $h$, we compute predictions on new slides $B_{I^{\prime}}$ as
$\textrm{sign}\left(\max_{i\in I^{\prime}}h(\mathbb{P}_{i})>0\right)$,
although the threshold can also be modified away from 0 to tailor the
sensitivity or specificity of the classifier.

The constant $b$ is determined by solving the following equation for a given
spot $i$, or by taking the average over all such equations, as has been
suggested for SVM (Burges, 1998):

$\begin{gathered}Y_{B(i)}\left(\sum_{j\in
E}\alpha_{j}Y_{B(j)}K(\mathbb{P}_{i},\mathbb{P}_{j})+b\right)=1\\\
\textrm{s.t. }0<\alpha_{i}<C\enskip\text{if }Y_{B(i)}=+1\quad\text{or}\quad
0<\sum_{j\in B(i)}\alpha_{j}<C\enskip\text{if }Y_{B(i)}=-1\end{gathered}$ (4)

Problem (1) is equivalent to the minimum of all choices of selection functions
$s$ of Problem (3). The last step to solving the intractability in (1) when
$w$ and $\phi$ are infinite-dimensional is to utilize the “kernel trick” by
replacing $\langle\phi(\mathbb{P}),\phi(\mathbb{Q})\rangle$ with a kernel
function $K(\mathbb{P},\mathbb{Q})$. While kernels on vector data are commonly
seen, a few works consider kernels appropriate for distributional instances.
The following sub-section describes the support measure machine (SMM) kernel,
an approach by Muandet et al. (2012) and Muandet et al. (2017), that offers
flexibility in choosing the underlying non-linearity similar to kernels on
vector data.

####  2.3.2 SMM kernel

The approach from Muandet et al. (2012) and Muandet et al. (2017) involves a
family of kernels defined on probability distributions, which they use to
learn a classifier from the SVM dual problem.

Given a reproducing kernel $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$
associated with a reproducing kernel Hilbert space (RKHS) $\mathcal{H}$, we
can define $K:\mathcal{P}\times\mathcal{P}\to\mathbb{R}$ by

$K(\mathbb{P}_{i},\mathbb{P}_{j})=\mathbb{E}_{X\sim\mathbb{P}_{i},Y\sim\mathbb{P}_{j}}\left[k(X,Y)\right]=\int\int
k(x,y)d\mathbb{P}_{i}(X)d\mathbb{P}_{j}(Y)$

They show under standard assumptions that $K$ is a positive definite kernel on
$\mathcal{P}$ and that $K$ works as an inner-product of so-called “mean maps”.
Let $\mu:\mathcal{P}\to\mathcal{H}$ denote the mean map:
$\mu(\mathbb{P})=\int_{\mathcal{X}}k(x,\cdot)d\mathbb{P}(x)$. It follows that
$K(\mathbb{P},\mathbb{Q})=\langle\mu(\mathbb{P}),\mu(\mathbb{Q})\rangle_{\mathcal{H}}$.

In some cases, this kernel can be computed in closed form (Muandet et al.,
2012). When this is not possible, we can approximate
$K(\mathbb{P}_{i},\mathbb{P}_{j})$ when the distributions manifest as the
empirical distributions $\hat{\mathbb{P}}_{i}$ and $\hat{\mathbb{P}}_{j}$ of
random samples $\\{x_{i,l}\\}_{l=1}^{r_{i}}$ and
$\\{z_{j,m}\\}_{m=1}^{r_{j}}$, respectively:

$K_{\textrm{emp}}(\hat{\mathbb{P}}_{i},\hat{\mathbb{P}}_{j})=\frac{1}{r_{i}\cdot
r_{j}}\sum_{l=1}^{r_{i}}\sum_{m=1}^{r_{j}}k(x_{i,l},z_{j,m})$

Once the kernel $K(\mathbb{P}_{i},\mathbb{P}_{j})$ is computed for all
combinations $i,j$, it can be plugged into an SVM dual problem such as (3).

####  2.3.3 Algorithm

Since the dual problem (3) relies on a fixed selection function $s$, we can
use it in a heuristic algorithm that updates the model at each step based on a
solution to the convex SVM problem (3) and subsequently updates the selection
function $s$. This idea is similar to that of the Expectation-Maximization
(E-M) algorithm (Moon, 1996), the Convex Concave Procedure (Yuille and
Rangarajan, 2003), and a non-convex SVM with weak labels studied previously
(Andrews, Tsochantaridis and Hofmann, 2003). We describe the algorithm in
detail in Algorithm 1.

1

Input: training data $\left\\{Y_{B(i)},\mathbb{P}_{i}\right\\}_{i=1}^{n}$,
embedding kernel $k$.

2 initialize _ $s(I)$ as a random selector function, $\forall I$ such that
$Y_{I}=1$ _

3 Pre-compute individual kernel matrix entries
$K(\mathbb{P}_{i},\mathbb{P}_{j})$ for all combinations $i,j$ in the data

4 while  _(selector variables $s(I)$ have changed)_ do

5 $E\leftarrow\left(\bigcup_{I:Y_{I}=-1}\\{i\in
I\\}\right)\cup\left(\bigcup_{I:Y_{I}=+1}s(I)\right)$

6 Solve for $\\{\alpha_{i}\\}_{i\in E}$ in Problem (3) based on the instances
within $E$

7 Compute $b$ from (4), or an average over all eligible equations

8 $h(\mathbb{P}_{i})\leftarrow\sum_{j\in
E}\alpha_{j}Y_{B(j)}K(\mathbb{P}_{i},\mathbb{P}_{j})+b$ for all instances
$\mathbb{P}_{i}$ in positive groups

9 $s(I)\leftarrow\arg\max_{i\in I}h(\mathbb{P}_{i})\;\forall I$ such that
$Y_{I}=1$

10

11 end while

return SMM classifier $h(\mathbb{P})=\sum_{j\in
E}\alpha_{j}Y_{B(j)}K(\mathbb{P},\mathbb{P}_{j})+b$

Algorithm 1 mi-smm (heuristic)

This algorithm can be improved by only computing entries of
$K(\mathbb{P}_{i},\mathbb{P}_{j})$ as needed in the set $E$, which will
increase the speed slightly. Additionally, one could perform Algorithm 1 over
numerous random initial selections $s$ and combine the outputs by ensembling
or choosing the highest performing version from a hold-out data set. It is
also common to limit the number of times selector variables can change to
avoid getting stuck in the while block forever.

Even though Algorithm 1 and similar approaches empirically show good
performance, there is no guarantee that it will find the globally optimal
solution to (1) (Andrews, Tsochantaridis and Hofmann, 2003; Yuille and
Rangarajan, 2003). The algorithm can get stuck in local minimum, and it is
difficult to know when this occurs.

###  2.4 A direct algorithm from MIQP formulation

As an alternative to the heuristic approach in Section 2.3, we can solve for
the MI-SMM classifier in a more direct fashion. We call this approach mi-smm
(miqp). For simplicity of construction, assume first that the form of $\phi$
in Problem (1) maps explicitly to a real-valued feature vector:
$\phi(\mathbb{P}_{i})=z_{i}\in\mathbb{R}^{p}$. We define index sets of the
spots based on outcome labels as $\mathcal{I}_{+}=\\{I:Y_{I}=1\\}$ and
$\mathcal{I}_{-}=\\{I:Y_{I}=-1\\}$. Using a well-known optimization trick
(Lemma provided in the Supplementary Materials), Problem (1) is equivalent to

$\begin{split}&\min_{\boldsymbol{w},b,\boldsymbol{\xi}\geq
0,\boldsymbol{\zeta}\in\\{0,1\\}}\frac{1}{2}\left\lVert\boldsymbol{w}\right\rVert^{2}+C\sum_{I}\xi_{I}\\\
\textrm{s.t. \quad}&-\langle\boldsymbol{w},z_{i}\rangle-b\geq
1-\xi_{I},\enskip\forall i\in I,\enskip\forall I\in\mathcal{I}_{-}\\\
&+\langle\boldsymbol{w},z_{i}\rangle+b\geq
1-\xi_{I}-L\cdot\zeta_{I,i},\enskip\forall i\in I,\enskip\forall
I\in\mathcal{I}_{+}\\\ &\sum_{i\in I}\zeta_{I,i}\leq|I|-1,\enskip\forall
I\in\mathcal{I}_{+}\end{split}$ (5)

provided that $L$ is a sufficiently large scalar. In the above formulation,
$|I|$ is the cardinality of $I$,
$\boldsymbol{\zeta}=\left[\zeta_{I,i}\right]_{I\in\mathcal{I}_{+},i\in I}$,
and $\boldsymbol{\zeta}\in\\{0,1\\}$ means
$\zeta_{I,i}\in\\{0,1\\}\enskip\forall I\in\mathcal{I}_{+}$, $i\in I$. This
unpacks the $\max$ constraint for positively labeled slides through the
introduction of binary variables $\zeta_{I,i}$. For a given slide-level group
$I$, not all of $\zeta_{I,i}$ can be 1, which ensures that
$\langle\boldsymbol{w},z_{i}\rangle+b\geq 1-\xi_{I}$ for at least one spot $i$
in each slide.

####  2.4.1 Feature map approximation

An unfortunate aspect of the formulation in (5) is that it is still non-
convex, leaving out the option of using duality theory to find an equivalent
formulation that directly uses the kernel $K(\mathbb{P},\mathbb{Q})$, as was
done in Section 2.3.1. When the form of $\phi$ does not map to a real-valued
feature vector, Problem (5) is intractable. To get around this challenge, we
leverage the theory developed on feature map approximations (Minh, Niyogi and
Yao, 2006; Rahimi and Recht, 2008; Vedaldi and Zisserman, 2012) to find a
real-valued function $\tilde{\phi}$ that approximates a kernel $k$ through
inner products:

$k(x,z)\approx\langle\tilde{\phi}(x),\tilde{\phi}(z)\rangle\quad\forall
x,z\in\mathcal{X}$

There is little discussion of feature map approximations on probability
distributions. Lemma 2.2 below provides some intuition for the SMM kernel in a
special case. The empirical kernel $K_{emp}$ summarizes the underlying feature
map by taking the mean of that map over the distribution samples, and non-
linearity is introduced by the embedding kernel $k$ and its feature map
$\phi$. We provide a short proof in the Supplementary Material since Muandet
et al. (2012) do not describe this result.

######  Lemma 2.2.

Let $k:X\times X\to\mathbb{R}$ from a RKHS admit a feature map $\phi:X\to
F\subset\mathbb{R}^{p}$ so that $k(x,z)=\langle\phi(x),\phi(z)\rangle$. Then

$K_{\textrm{emp}}(\hat{\mathbb{P}}_{i},\hat{\mathbb{P}}_{j})=\langle\bar{\phi}(\hat{\mathbb{P}}_{i}),\bar{\phi}(\hat{\mathbb{P}}_{j})\rangle$

where
$\bar{\phi}(\hat{\mathbb{P}}_{i})=\frac{1}{r_{i}}\sum_{l=1}^{r_{i}}\phi(x_{i,l})$
and
$\bar{\phi}(\hat{\mathbb{P}}_{j})=\frac{1}{r_{j}}\sum_{m=1}^{r_{j}}\phi(x_{j,m})$
are the empirical means of the random samples evaluated on the feature map.

Give a feature map approximation $\tilde{\phi}$ for the underlying embedding
kernel $k$, Lemma 2.2 provides a natural extension to a feature map
approximation for $K$:

$K_{emp}(\hat{\mathbb{P}}_{i},\hat{\mathbb{P}}_{j})\approx\left\langle\frac{1}{r_{i}}\sum_{l=1}^{r_{i}}\tilde{\phi}(x_{i,l}),\frac{1}{r_{j}}\sum_{m=1}^{r_{j}}\tilde{\phi}(x_{j,m})\right\rangle,$

for the random samples $\\{x_{i,l}\\}_{l=1}^{r_{i}}$ and
$\\{x_{j,m}\\}_{m=1}^{r_{j}}$ from empirical distributions
$\hat{\mathbb{P}}_{i}$ and $\hat{\mathbb{P}}_{j}$, respectively.
Computationally, this is much faster than approximating $K$ as
$\frac{1}{r_{i}\cdot
r_{j}}\sum_{l=1}^{r_{i}}\sum_{m=1}^{r_{j}}k(x_{i,l},x_{j,m})$, since it
requires calculating $k(x_{i,l},x_{j,m})$ from
$\langle\tilde{\phi}(x_{i,l}),\tilde{\phi}(x_{j,m})\rangle$ for all
combinations of samples.

In this paper, we use the Nyström feature map (Williams and Seeger, 2001) to
approximate the embedding kernel $k$ because of its simplicity and excellent
performance. Yang et al. (2012) provide a detailed description and intuition
for this feature map, and Chatalic et al. (2022) have recently applied Nyström
methods to kernels on probability distributions. We modify the Nyström
algorithm by taking a stratified sub-sample from the concatenation of all
samples $\left\\{\\{x_{i,l}\\}_{l=1}^{r_{i}}\right\\}_{i=1}^{n}$ instead of a
simple sub-sample, ensuring each slide-level group is represented almost
uniformly (details follow in Algorithm 2).

####  2.4.2 Algorithm

Algorithm 2 connects the Nyström feature map procedure, SMM kernel
approximation, and MIQP from Problem (5) into the direct approach, mi-smm
(miqp). It requires a MIQP solver and two integer constants, $m_{1}$ and
$m_{2}$, for the kernel rank approximation and sub-sample size, respectively.
The choice of $m_{1}$ and $m_{2}$ can be made based on computational
feasibility or chosen via cross-validation.

Input: kernel $k$, training data
$\left\\{Y_{B(i)},\\{x_{i,l}\\}_{l=1}^{r_{i}}\right\\}_{i=1}^{n}$, and
$m_{1},m_{2}\in\mathbb{N}$ s.t. $m_{1}\leq m_{2}\leq n$.

1 begin create the Nyström feature map

2 Take $m_{2}$ stratified sub-samples from
$\left\\{\\{x_{i,l}\\}_{l=1}^{r_{i}}\right\\}_{i=1}^{n}$ as
$\hat{x}_{1},\dots,\hat{x}_{m_{2}}$

3
$\hat{K}\leftarrow\left[k\left(\hat{x}_{i},\hat{x}_{j}\right)\right]_{m_{2}\times
m_{2}}$

4 Find eigenvalue, eigenvector pairs of $\hat{K}$ as
$(\hat{\lambda}_{j},\hat{\mathbf{v}}_{j}),j=1,\dots,m_{2}$ ranked by
decreasing order of eigenvalues.

5
$\hat{V}\leftarrow\left(\hat{\mathbf{v}}_{1},\dots,\hat{\mathbf{v}}_{m_{1}}\right)$
and
$\hat{D}\leftarrow\textrm{diag}(\hat{\lambda}_{1},\dots,\hat{\lambda}_{m_{1}})$,
using $m_{1}$ out of $m_{2}$ eigen-pairs.

6 Return the feature map
$\tilde{\phi}(x)=\hat{D}^{-1/2}\hat{V}^{\top}\left(k\left(x,\hat{x}_{1}\right),\ldots,k\left(x,\hat{x}_{m_{2}}\right)\right)^{\top}$

7

8 end

9

10$z_{i}\leftarrow\frac{1}{r_{i}}\sum_{l=1}^{r_{i}}\tilde{\phi}(x_{i,l})$, for
$i=1,\dots,n$

11 return $(\boldsymbol{w},b)$ from Problem (5) via an MIQP solver

Algorithm 2 mi-smm (miqp)

The advantage of Algorithm 2 is that it is guaranteed to find a global
solution to (1) under the feature map approximation. In contrast to Algorithm
1, it is potentially much slower, as MIQP problems have exponential
computational complexity and MIQP problems are NP-complete (Pia, Dey and
Molinaro, 2017). Despite this, solvers are improving their speed and can solve
small problems efficiently through heuristic modifications of the branch and
bound algorithm (Lazimy, 1982; Gurobi Optimization, 2021).

##  3 Simulations

###  3.1 Simulation data

To evaluate our proposed approach, we consider four distinct simulation
scenarios that contain data with weakly supervised labels and a distributional
instance structure. Each scenario is parametrized by three quantities that
control the size of the data: number of top-level groups (20, 50), number of
instances per group (3, 6), and the sample size from each distributional
instance (20, 50, 100), for 12 total combinations. To construct the data, we
first generate instance labels $y_{i}$ as i.i.d. Bernoulli random variables
with $P(y_{i}=1)=0.15$ for all instances. Then the group label $Y_{I}$ is
determined from the instance label through Assumption 2.1. The samples
$\\{x_{i,j}\\}_{j=1}^{r_{i}}$ are generated depending on the four distinct
scenarios, where $\mathcal{X}=\mathbb{R}^{10}$, as described in Table 1.

Table 1: Summary of Simulations Setup

| Differing Components |
---|---|---
Setting | When $y_{i}=+1$ | When $y_{i}=-1$ | Remaining Components
Student $t$ vs. Normal | $(x_{i,j}^{(1)},\dots,x_{i,j}^{(5)})\sim\textrm{MVT}{\textsuperscript{1}}\left(\nu=3,\boldsymbol{0}_{5},\frac{1}{3}\boldsymbol{I}_{5}\right)$ | $(x_{i,j}^{(1)},\dots,x_{i,j}^{(5)})\sim\textrm{MVN}{\textsuperscript{2}}\left(\boldsymbol{0}_{5},\boldsymbol{I}_{5}\right)$ | $\textrm{MVN}\left(\boldsymbol{0}_{5},\boldsymbol{I}_{5}\right)$
Covariance differences | $(x_{i,j}^{(1)},x_{i,j}^{(2)})\sim\mathrm{MVN}\left(\mathbf{0}_{2},\left(\begin{smallmatrix}{}1&-0.5\\\ -0.5&1\end{smallmatrix}\right)\right)$ | $(x_{i,j}^{(2)},x_{i,j}^{(3)})\sim\mathrm{MVN}\left(\mathbf{0}_{2},\left(\begin{smallmatrix}{}1&0.5\\\ 0.5&1\end{smallmatrix}\right)\right)$ | $\textrm{MVN}\left(\boldsymbol{0}_{8},\boldsymbol{I}_{8}\right)$
Mean differences | $(x_{i,j}^{(1)},\dots,x_{i,j}^{(5)})\sim\mathrm{MVN}\left(\mathbf{0.2}_{5},\boldsymbol{I}_{5}\right)$ | $(x_{i,j}^{(1)},\dots,x_{i,j}^{(5)})\sim\mathrm{MVN}\left(\mathbf{0}_{5},\boldsymbol{I}_{5}\right)$ | $\textrm{MVN}\left(\boldsymbol{0}_{5},\boldsymbol{I}_{5}\right)$
Large covariance differences | $(x_{i,j}^{(1)},\dots,x_{i,j}^{(5)})\sim\mathrm{MVN}\left(\mathbf{0}_{5},\mathbf{\Sigma}_{S4}\right)$ | $(x_{i,j}^{(6)},\dots,x_{i,j}^{(10)})\sim\mathrm{MVN}\left(\mathbf{0}_{5},\mathbf{\Sigma}_{S4}\right)$ | $\textrm{MVN}\left(\boldsymbol{0}_{5},\boldsymbol{I}_{5}\right)$

  * We use $\boldsymbol{0}_{p}$ to denote a length $p$ vector with all elements being 0 (with similar notation for other scalars) and $\boldsymbol{I}_{p}$ to denote the $p\times p$ identity matrix.

  * 1

$\textrm{MVT}\left(\nu,\boldsymbol{\delta},\mathbf{\Sigma}\right)$ is the
multivariate $t$ distribution with degrees of freedom $\nu\geq 3$, mean
$\boldsymbol{\delta}$, covariance matrix $\frac{\nu}{\nu-2}\mathbf{\Sigma}$,
and density
$f(x;\nu,\boldsymbol{\delta},\mathbf{\Sigma})=\frac{\Gamma\left(\frac{v+k}{2}\right)}{\Gamma\left(\frac{v}{2}\right)\sqrt{(v\pi)^{k}|\mathbf{\Sigma}|}}\cdot\left(1+\frac{(x-\boldsymbol{\delta})^{t}\mathbf{\Sigma}^{-1}(x-\boldsymbol{\delta})}{v}\right)^{-\frac{v+k}{2}}$.

  * 2

$\textrm{MVN}\left(\boldsymbol{\mu},\mathbf{\sigma}\right)$ is the
multivariate normal distribution with mean $\boldsymbol{\mu}$, covariance
matrix $\mathbf{\Sigma}$, and density
$f\left(x;\boldsymbol{\mu},\mathbf{\Sigma}\right)=\frac{1}{\sqrt{(2\pi)^{k}|\mathbf{\Sigma}|}}\cdot\exp\left(-\frac{1}{2}(x-\boldsymbol{\mu})^{\mathrm{T}}\boldsymbol{\Sigma}^{-1}(x-\boldsymbol{\mu})\right)$.

  * 3

$\mathbf{\Sigma}_{S4}$ is a $5\times 5$ covariance matrix with 1 on the
diagonal and 0.5 in every other entry.

###  3.2 Methods compared

Both algorithms for our proposed non-convex SVM solution, mi-smm (heuristic)
and mi-smm (miqp), were tested on the simulation scenarios and compared to two
ad hoc approaches. The comparison methods simplify the data to a structure
that can be solved by previous approaches.

si-smm: The first approach, which we call Single-Instance Support Measure
Machines (si-smm), relies on imputing the spot label as the slide label
$\tilde{y_{i}}=Y_{I}$. For the weakly-supervised setting, this results in the
spot-level data set of $\\{\tilde{y}_{i},\mathbb{P}_{i}\\}_{i=1}^{n}$, on
which we can train an SMM classifier $h$ (Muandet et al., 2012). A new slide
$B_{I^{\prime}}=\\{\mathbb{P}_{i}:i\in I^{\prime}\\}$ can be classified as
$\hat{Y}_{I^{\prime}}=\max_{i\in I^{\prime}}h(\mathbb{P}_{i})$, following
Assumption 2.1. This technique has previously been considered for data with
the weakly supervised label structure, but not for distributional instance
data (Ray and Craven, 2005; Alpaydın et al., 2015).

mi-svm: The second approach focuses on summarizing each distributional
instance into a feature vector. We define a function
$\varphi:\mathcal{P}\to\mathbb{R}^{p}$ which summarizes a distribution into a
feature vector: $\varphi(\mathbb{P}_{i})=z_{i}$. For example, this could
include the mean, standard deviation, or correlation of $\mathbb{P}_{i}$.
Applying $\varphi$ results in a group of feature vectors
$B_{I}^{\varphi}=\\{z_{i}:i\in I\\}$ corresponding to each original slide-
level group $B_{I}$. We apply a similar non-convex SVM algorithm to the data
set $\\{Y_{I},B_{I}^{\varphi}\\}$, called mi-svm from Andrews, Tsochantaridis
and Hofmann (2003).

For this approach, we consider three sets of summary functions on the instance
distributions $\hat{\mathbb{P}}_{i}=\\{x_{i,j}\\}_{j=1}^{r_{i}}$ where
$x_{i,j}\in\mathcal{X}\subseteq\mathbb{R}^{d}$:

  * •

$\varphi_{\textrm{univ}1}=\left[\textrm{mean}(\hat{\mathbb{P}}_{i}),\textrm{SD}(\hat{\mathbb{P}}_{i})\right]\in\mathbb{R}^{2d}$,
which includes the sample mean and standard deviation of $\mathbb{P}_{i}$.

  * •

$\varphi_{\textrm{univ}2}=\left[\textrm{skew}(\hat{\mathbb{P}}_{i}),\textrm{kurt}(\hat{\mathbb{P}}_{i}),\textrm{Q1}(\hat{\mathbb{P}}_{i}),\textrm{Q3}(\hat{\mathbb{P}}_{i})\right]\in\mathbb{R}^{4d}$,
which includes the sample skew, sample kurtosis, and the 25th and 75th
percentile of $\mathbb{P}_{i}$.

  * •

$\varphi_{\textrm{cor}}$: 2-way correlations. All 2-way correlations between
each pair of features are computed using the sample Pearson correlation
coefficient, so that
$\varphi_{\textrm{cor}}(\hat{\mathbb{P}}_{i})\in\mathbb{R}^{p}$ where
$p=\binom{d}{2}$.

From these three sets of summary statistics, we evaluate four distinct
combinations. We name each combination based on the concatenated set of
included functions: mi-svm (univ1), mi-svm (univ1, univ2), mi-svm (univ1,
cor), mi-svm (univ1, univ2, cor). All four utilize the first set of univariate
features $\varphi_{\textrm{univ}1}$, while the second set of univariate
features $\varphi_{\textrm{univ}2}$ and the multivariate features
$\varphi_{\textrm{cor}}$ are toggled on or off.

We use the Gaussian kernel $k(x,z)=\exp\left(-\frac{\left\lVert
x-z\right\rVert^{2}}{2\sigma^{2}}\right)$ in si-smm and mi-smm as the
embedding kernel, and in all mi-svm approaches to introduce non-linearity.

###  3.3 Implementation details

Each combination of data size and data scenario is evaluated similarly. We
generate 100 data sets to serve as training data. For training, we employ a
5-fold cross-validated grid-search to determine the optimal choice of $C$ and
$\sigma$ in the SVM optimization problem and Gaussian kernel, respectively.
The parameter $C$ is weighted inversely proportional to the number of positive
groups and the number of negative instances to help alleviate potential class
imbalance. A corresponding testing data set of 500 groups is generated for
each training data set. We evaluate model performance on the predicted top-
level group labels from the area under the receiver operating characteristic
curve (AUROC).

On each data set, we apply the seven distinct methods described previously.
Since there are 10 covariates for each data point, the mi-svm methods have 20,
60, 65, and 105 total features for uinv1, uinv1+uinv2, uinv1+cor,
uinv1+uinv2+cor, respectively. For the mi-smm (miqp), we take $m_{1}=240$
samples and use rank $m_{2}=240$ for the Nyström kernel approximation.

###  3.4 Performance evaluation

Figure 3 compares the AUROC performance of each method for the data sets
generated with 50 groups, 3 instances per group, and 50 samples from each
instance distribution. Overall, mi-smm (miqp) and si-smm show strong
performance in most scenarios. In Scenario 3, where the distributions differ
in their means, si-smm outperforms mi-smm (miqp) by a fair margin, but it is
less consistent in Scenarios 2 and 4.

Figure 3: Scatter plot of AUROC vs method for four simulation scenarios when
there are 100 top-level groups, 3 instance per group, and 50 distribution
samples per instance. The points in each subplot represent performance on 100
individually simulated combinations of training and testing data sets. The
width of each cluster of points closely corresponds to the number of points
with the given performance. A black line shows the mean AUROC within each
scenario and method.

Figure 4 summarizes the results across the 1,200 data sets (12 sample size
combinations $\times$ 100 data sets) for each of the simulation scenarios by
providing the average classifier AUROC rank (1-7). The trends illustrated are
similar to those in Figure 3. The performance of the mi-svm methods depends
heavily on the scenario. They perform better when the summary features
correspond to the scenario differences. However, it’s not sufficient to throw
in all features since mi-svm (univ1, univ2, cor) is not the strongest
performing of the mi-svm methods in any data scenario. The full results,
including plots of AUROC for each scenario and sample size combination, are
available in the Supplementary Materials.

Figure 4: Average ranks of test AUROC across all data sets. Each panel
represents one simulation scenario. Ranks range from 1-7, with a rank close to
1 indicating the best performance.

##  4 Analysis of collagen features data

Using our non-convex SVM methods, mi-smm, proposed in Section 2, we analyze
the collagen features data. This section includes an overview of the key
variables and an evaluation of how well our methodology predicts the tumor
label compared to other approaches.

###  4.1 Data set overview

244 DCIS patients have a total of 348 tissue sample slides, which represent
the top-level groups. We seek to classify the 235 tumor tissue slides
($Y_{I}=1$) vs. the 113 normal tissue slides ($Y_{I}=-1$). Instances are
represented through the multiple spots chosen on each image (min 1, max 8), as
shown in Figure 1. These spots yield a weakly-supervised label data structure.
We assume that tumorous tissue samples have at least one spot with features
indicating tumor, while non-tumorous samples have no spots indicating tumor,
which aligns with Assumption 2.1.

Each spot was processed using the second harmonic generating (SHG) imaging
(Chen et al., 2012) to find morphological features on numerous collagen fibers
(min 73, max 1092). Since there are many collagen fibers within each spot, and
the relationship between the fibers may be predictive of tumor tissue, we
consider each spot as an empirical probability distribution, where fibers are
sampled from the spot. Figure 1 shows that the number of fibers is commonly
between 350 – 750.

This SHG process creates 20 features such as length, width, curvature, and
measurements of alignment and density. Table 2 provides an overview of all
features and compares the mean of subject-level averages between patients
having tumorous and non-tumorous slides. Many features show significant
differences in the two groups, most notably the distance of one fiber to other
nearby fibers, the angle to the nearest relative boundary, and the box
alignment. Some features are heavily skewed, which we log-transform. All are
centered and scaled in training.

Table 2: Comparison of subject-level averages of morphologic features. | Mean |
---|---|---
Variable | Non-tumor | Tumor | P value
Total length | 57.97 | 58.63 | 0.0570
End-to-end length | 52.80 | 53.51 | 0.0345
Width | 5.55 | 5.46 | 0.0023
Nearest relative boundary angle | 8.33 | 10.44 | ¡0.0001
Curvature | 0.91 | 0.92 | 0.0783
Distance to nearest 2 | 24.30 | 25.14 | ¡0.0001
Distance to nearest 4 | 31.33 | 32.55 | ¡0.0001
Distance to nearest 8 | 42.00 | 43.91 | ¡0.0001
Distance to nearest 16 | 58.13 | 61.30 | ¡0.0001
Mean nearest distance | 38.94 | 40.73 | ¡0.0001
Std nearest distance | 14.80 | 15.83 | ¡0.0001
Alignment of nearest 2 | 0.76 | 0.77 | 0.0077
Alignment of nearest 4 | 0.64 | 0.65 | 0.0014
Alignment of nearest 8 | 0.56 | 0.58 | 0.0015
Alignment of nearest 16 | 0.51 | 0.53 | 0.0069
Mean nearest alignment | 0.62 | 0.63 | 0.0031
Std nearest alignment | 0.17 | 0.17 | 0.0035
Box alignment 32 | 0.87 | 0.88 | 0.0086
Box alignment 64 | 0.68 | 0.70 | ¡0.0001
Box alignment 128 | 0.54 | 0.57 | ¡0.0001

###  4.2 Predictive performance

We compare the predictive performance of our proposed non-convex SVM framework
to the comparison methods from Section 3.2. We use 10 replications of a
10-fold cross-validation procedure to evaluate each method (Kim, 2009). In
each training fold, a 5-fold cross-validated grid-search determines the
optimal choice of $C$ and $\sigma$ in the SVM optimization problem and
Gaussian kernel, respectively. Test performance is combined across the testing
folds through averaging. For the Nyström approximation in mi-smm (miqp), we
take $m_{1}=$ 360 samples and use rank $m_{2}=$ 360\. Due to the data set
size, the mi-smm (miqp) algorithm takes considerable time with state-of-the-
art optimization solvers (Gurobi Optimization, 2021); we accept the best
solution after 10 minutes.

Table 3: Summary of predictive performance. | AUROC
---|---
Method | Mean | SD
MI-SMM (MIQP) | 0.768 | 0.023
MI-SMM (HEURISTIC) | 0.695 | 0.017
SI-SMM | 0.695 | 0.025
MI-SVM (UNIV1) | 0.761 | 0.013
MI-SVM (UNIV1, UNIV2) | 0.772 | 0.016
MI-SVM (UNIV1, COR) | 0.765 | 0.015
MI-SVM (UNIV1, UNIV2, COR) | 0.759 | 0.018

The AUROC performance for each method on the fiber-feature data set is
summarized in Table 3. The top methods on this data are mi-svm (univ1, univ2),
mi-smm (miqp), and mi-svm (univ1, cor), although all of the mi-svm methods do
well with an AUROC near 0.77. The mean and standard deviation of features from
$\varphi_{\textrm{univ}1}$ are adequate to predict tumor or normal tissue;
however, this performance can be slightly improved with additional univariate
and multivariate features. Importantly, mi-smm (miqp) is comparable with these
methods without requiring a potentially arbitrary choice in these features.

##  5 Discussion

mi-smm, a non-convex, maximum-margin-based learning approach, was constructed
to classify tumorous tissue vs. non-tumorous tissue from a data set of
collagen fiber features within the tumor microenvironment. It directly handles
both aspects that make the data set unique: the weakly-supervised label
structure stemming from the multiple spots per subject and the distributional
instance structure stemming from many fibers per spot. In addition, mi-smm
allows for non-linear transformations of the data through an embedding kernel
and two separate algorithms that trade off computational efficiency for
optimality guarantees. In its entirety, it offers flexibility to adapt to
similarly-structured, unseen data sets.

Our approach shows strong performance across all of the experimental settings
we tested in this work. The other ad-hoc approaches worked well in either the
experimental simulations or the motivating data set, but not both, showing
that they were more data-dependent. In contrast, mi-smm was more adaptable to
the data sets we tested on. Further testing on additional data is needed to
make broad conclusions, but this shows a promising start.

This work outlines the first combination of weakly-supervised labels and
distributional instances that we are aware of, but data sets with similar
structures may likely surface in the future. Thinking of instances as
distributions is justified when we observe a large sample of random events
from a unique identity. As the size and number of data sets grow in the
future, examples of this structure may become more common. Such examples could
include: the behavior of a single user on an internet platform across many
days; weights from trained inner layers of a neural network; or images under a
sample of scaling, translation, and rotation transformations (see Muandet et
al., 2012). The weakly-supervised setting has already occurred in previous
works (Dietterich, Lathrop and Lozano-Pérez, 1997; Andrews, Tsochantaridis and
Hofmann, 2003; Ray and Craven, 2005; Campanella et al., 2019; Chen et al.,
2022). It spans a host of applications, such as drug discovery, tumor
classification from whole slide images, and text document classification. We
hope future researchers can use our methodology on data with both
distributional instances and weakly-supervised label structure.

One limitation of the motivating data set is the assumption of how tissue
slide labels relate to the spot labels. If, for a given tumorous tissue
sample, no spots are observed that contain the features of cancer, this could
violate our assumption and provide incorrect label information to the training
algorithm. These violations are unlikely in practice because of the expert
selection of significant spots. Another limiting factor is the computation in
large data sets. While mi-smm (miqp) works better than mi-smm (heuristic), the
computational cost of solving the MIQP problem directly is high, as the number
of integer combinations multiplies rapidly with an increasing number of
instances. Setting a time limit on this algorithmic approach worked in
practice on the collagen fiber data set, despite the large optimality gap
observed. Alternatively, the heuristic algorithm for mi-smm is typically quite
fast, even for large data sets.

All approaches described in this work are available in an open-source R
package, `mildsvm`.333Available at github.com/skent259/mildsvm The package is
designed for building models with the familiar formula interface or through
data frames. It also includes useful tools for SVM-based methods, including
methods for the kernel feature approximation of Section 2.4.1, calculation of
kernels on probability distributions (Muandet et al., 2012), and summarization
of distributional instances to features. We believe this work can provide a
basis to expand the literature on non-convex SVMs, potentially through
additional kernels that work on probability distribution data or alternative
non-convex SVM problems targeting different assumptions of data structure.

Research reported in this publication was partially supported by the National
Cancer Institute of the National Institutes of Health under Award Number
T32LM012413 and the University of Wisconsin Carbone Cancer Center Support
Grant P30 CA014520. The content is solely the responsibility of the authors
and does not necessarily represent the official views of the National
Institutes of Health.

Supplement to “Non-convex SVM for cancer diagnosis based on morphologic
features of tumor microenvironment” The supplementary material contains proofs
of all theoretical results and additional simulation results (pdf file). An
associated R package mildsvm is currently available on GitHub
(github.com/skent259/mildsvm) and will be made available on the Comprehensive
R Archive Network (CRAN). Full scripts to recreate the simulation results and
data analysis are provided in the following GitHub repository for
reproducibility: github.com/skent259/mildsvm-sims; the motivating data set is
excluded from the public repository but can be made available upon request.

## References

  * Alpaydın et al. (2015) [author] Alpaydın, EthemE., Cheplygina, VeronikaV., Loog, MarcoM. and Tax, David M. J.D. M. J. (2015). Single- vs. multiple-instance classification. Pattern Recognition 48 2831–2838. 10.1016/j.patcog.2015.04.006
  * Anderson and Simon (2020) [author] Anderson, Nicole M.N. M. and Simon, M. CelesteM. C. (2020). The tumor microenvironment. Current Biology 30 R921–R925. 10.1016/j.cub.2020.06.081
  * Andrews, Tsochantaridis and Hofmann (2003) [author] Andrews, StuartS., Tsochantaridis, IoannisI. and Hofmann, ThomasT. (2003). Support Vector Machines for Multiple-Instance Learning. In Advances in Neural Information Processing Systems 15 (S.S. Becker, S.S. Thrun and K.K. Obermayer, eds.) 577–584. MIT Press.
  * Arendt et al. (2010) [author] Arendt, Lisa M.L. M., Rudnick, Jenny A.J. A., Keller, Patricia J.P. J. and Kuperwasser, CharlotteC. (2010). Stroma in breast development and disease. Seminars in Cell & Developmental Biology 21 11–18. 10.1016/j.semcdb.2009.10.003
  * Baghban et al. (2020) [author] Baghban, RoghayyehR., Roshangar, LeilaL., Jahanban-Esfahlan, RanaR., Seidi, KhaledK., Ebrahimi-Kalan, AbbasA., Jaymand, MehdiM., Kolahian, SaeedS., Javaheri, TaherehT. and Zare, PeymanP. (2020). Tumor microenvironment complexity and therapeutic implications at a glance. Cell Communication and Signaling 18 59\. 10.1186/s12964-020-0530-4
  * Bejarano, Jordāo and Joyce (2021) [author] Bejarano, LeireL., Jordāo, Marta J. C.M. J. C. and Joyce, Johanna A.J. A. (2021). Therapeutic Targeting of the Tumor Microenvironment. Cancer Discovery 11 933–959. 10.1158/2159-8290.CD-20-1808
  * Boser, Guyon and Vapnik (1992) [author] Boser, Bernhard E.B. E., Guyon, Isabelle M.I. M. and Vapnik, Vladimir N.V. N. (1992). A training algorithm for optimal margin classifiers. In Proceedings of the fifth annual workshop on Computational learning theory 144–152. Association for Computing Machinery. 10.1145/130385.130401
  * Boyd and Vandenberghe (2004) [author] Boyd, StephenS. and Vandenberghe, LievenL. (2004). 5\. Duality. In Convex Optimization 1 ed. 215–287. Cambridge University Press, Cambridge, UK ; New York.
  * Burges (1998) [author] Burges, Christopher J. C.C. J. C. (1998). A Tutorial on Support Vector Machines for Pattern Recognition. Data Mining and Knowledge Discovery 2 121–167. 10.1023/A:1009715923555
  * Campanella et al. (2019) [author] Campanella, GabrieleG., Hanna, Matthew G.M. G., Geneslaw, LukeL., Miraflor, AllenA., Werneck Krauss Silva, VitorV., Busam, Klaus J.K. J., Brogi, EdiE., Reuter, Victor E.V. E., Klimstra, David S.D. S. and Fuchs, Thomas J.T. J. (2019). Clinical-grade computational pathology using weakly supervised deep learning on whole slide images. Nature Medicine 25 1301–1309. 10.1038/s41591-019-0508-1
  * Chatalic et al. (2022) [author] Chatalic, AntoineA., Schreuder, NicolasN., Rudi, AlessandroA. and Rosasco, LorenzoL. (2022). Nyström Kernel Mean Embeddings. arXiv:2201.13055 [cs, stat].
  * Chen et al. (2012) [author] Chen, XiyiX., Nadiarynkh, OlegO., Plotnikov, SergeyS. and Campagnola, Paul J.P. J. (2012). Second harmonic generation microscopy for quantitative analysis of collagen fibrillar structure. Nature Protocols 7 654–669. 10.1038/nprot.2012.009
  * Chen et al. (2022) [author] Chen, JunhuaJ., Zeng, HaiyanH., Zhang, ChongC., Shi, ZhenweiZ., Dekker, AndreA., Wee, LeonardL. and Bermejo, InigoI. (2022). Lung cancer diagnosis using deep attention-based multiple instance learning and radiomics. Medical Physics 49 3134–3143. 10.1002/mp.15539
  * Conklin et al. (2018) [author] Conklin, Matthew W.M. W., Gangnon, Ronald E.R. E., Sprague, Brian L.B. L., Gemert, Lisa VanL. V., Hampton, John M.J. M., Eliceiri, Kevin W.K. W., Bredfeldt, Jeremy S.J. S., Liu, YumingY., Surachaicharn, NuntidaN., Newcomb, Polly A.P. A., Friedl, AndreasA., Keely, Patricia J.P. J. and Trentham-Dietz, AmyA. (2018). Collagen Alignment as a Predictor of Recurrence after Ductal Carcinoma In Situ. Cancer Epidemiology and Prevention Biomarkers 27 138–145. 10.1158/1055-9965.EPI-17-0720
  * Dietterich, Lathrop and Lozano-Pérez (1997) [author] Dietterich, Thomas G.T. G., Lathrop, Richard H.R. H. and Lozano-Pérez, TomásT. (1997). Solving the multiple instance problem with axis-parallel rectangles. Artificial Intelligence 89 31–71. 10.1016/S0004-3702(96)00034-3
  * Ertekin, Bottou and Giles (2010) [author] Ertekin, SeydaS., Bottou, LeonL. and Giles, C LeeC. L. (2010). Nonconvex online support vector machines. IEEE Transactions on Pattern Analysis and Machine Intelligence 33 368–381. Publisher: IEEE.
  * Guan et al. (2020) [author] Guan, LeiL., Sun, TaoT., Qiao, Lin-boL.-b., Yang, Zhi-huiZ.-h., Li, Dong-shengD.-s., Ge, Ke-shiK.-s. and Lu, Xi-chengX.-c. (2020). An efficient parallel and distributed solution to nonconvex penalized linear SVMs. Frontiers of Information Technology & Electronic Engineering 21 587–603. Publisher: Springer.
  * Gurobi Optimization (2021) [author] Gurobi Optimization, LLCL. (2021). Mixed-Integer Programming (MIP) - A Primer on the Basics.
  * Kim (2009) [author] Kim, Ji-HyunJ.-H. (2009). Estimating classification error rate: Repeated cross-validation, repeated hold-out and bootstrap. Computational Statistics & Data Analysis 53 3735–3745. 10.1016/j.csda.2009.04.009
  * Laporte et al. (2013) [author] Laporte, LéaL., Flamary, RémiR., Canu, StéphaneS., Déjean, SébastienS. and Mothe, JosianeJ. (2013). Nonconvex regularizations for feature selection in ranking with sparse SVM. IEEE Transactions on Neural Networks and Learning Systems 25 1118–1130. Publisher: IEEE.
  * Lazimy (1982) [author] Lazimy, RafaelR. (1982). Mixed-integer quadratic programming. Mathematical Programming 22 332–349. 10.1007/BF01581047
  * Li et al. (2013) [author] Li, Yu-FengY.-F., Tsang, Ivor WI. W., Kwok, James TJ. T. and Zhou, Zhi-HuaZ.-H. (2013). Convex and scalable weakly labeled SVMs. Journal of Machine Learning Research 14\.
  * Lin et al. (2013) [author] Lin, H. T.H. T., Lee, S.S., Bui, N.N. and Honavar, V.V. (2013). Learning Classifiers from Distributional Data. In 2013 IEEE International Congress on Big Data 302–309. 10.1109/BigData.Congress.2013.47
  * Liu, Yao and Li (2016) [author] Liu, HongchengH., Yao, TaoT. and Li, RunzeR. (2016). Global solutions to folded concave penalized nonconvex learning. Annals of statistics 44 629\. Publisher: NIH Public Access.
  * Minh, Niyogi and Yao (2006) [author] Minh, Ha QuangH. Q., Niyogi, ParthaP. and Yao, YuanY. (2006). Mercer’s theorem, feature maps, and smoothing. In International Conference on Computational Learning Theory 154–168. Springer.
  * Moon (1996) [author] Moon, T. K.T. K. (1996). The expectation-maximization algorithm. IEEE Signal Processing Magazine 13 47–60. 10.1109/79.543975
  * Muandet et al. (2012) [author] Muandet, KrikamolK., Fukumizu, KenjiK., Dinuzzo, FrancescoF. and Schölkopf, BernhardB. (2012). Learning from distributions via support measure machines. Advances in neural information processing systems 25\.
  * Muandet et al. (2017) [author] Muandet, KrikamolK., Fukumizu, KenjiK., Sriperumbudur, BharathB. and Schölkopf, BernhardB. (2017). Kernel Mean Embedding of Distributions: A Review and Beyond. Foundations and Trends® in Machine Learning 10 1–141. 10.1561/2200000060
  * Pia, Dey and Molinaro (2017) [author] Pia, Alberto DelA. D., Dey, Santanu S.S. S. and Molinaro, MarcoM. (2017). Mixed-integer quadratic programming is in NP. Mathematical Programming 162 225–240. 10.1007/s10107-016-1036-0
  * Platt (1998) [author] Platt, JohnJ. (1998). Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines.
  * Rahimi and Recht (2008) [author] Rahimi, AliA. and Recht, BenjaminB. (2008). Random Features for Large-Scale Kernel Machines. In Advances in Neural Information Processing Systems 20 1177–1184. Curran Associates, Inc.
  * Ramsay (2006) [author] Ramsay, J. O.J. O. (2006). Functional Data Analysis. In Encyclopedia of Statistical Sciences American Cancer Society. 10.1002/0471667196.ess3138
  * Ray and Craven (2005) [author] Ray, SoumyaS. and Craven, MarkM. (2005). Supervised versus multiple instance learning: an empirical comparison. In Proceedings of the 22nd international conference on Machine learning. ICML ’05 697–704. Association for Computing Machinery, New York, NY, USA. 10.1145/1102351.1102439
  * Vedaldi and Zisserman (2012) [author] Vedaldi, AndreaA. and Zisserman, AndrewA. (2012). Efficient Additive Kernels via Explicit Feature Maps. IEEE Transactions on Pattern Analysis and Machine Intelligence 34 480–492. 10.1109/TPAMI.2011.153
  * Williams and Seeger (2001) [author] Williams, ChristopherC. and Seeger, MatthiasM. (2001). Using the Nyström Method to Speed Up Kernel Machines. In Advances in Neural Information Processing Systems 13 682–688. MIT Press.
  * Yang et al. (2012) [author] Yang, TianbaoT., Li, Yu-fengY.-f., Mahdavi, MehrdadM., Jin, RongR. and Zhou, Zhi-HuaZ.-H. (2012). Nyström Method vs Random Fourier Features: A Theoretical and Empirical Comparison. In Advances in Neural Information Processing Systems 25 476–484.
  * Yuille and Rangarajan (2003) [author] Yuille, A. L.A. L. and Rangarajan, AnandA. (2003). The Concave-Convex Procedure. Neural Computation 15 915–936. 10.1162/08997660360581958
  * Zeltz et al. (2020) [author] Zeltz, CédricC., Primac, IrinaI., Erusappan, PugazendhiP., Alam, JahedulJ., Noel, AgnesA. and Gullberg, DonaldD. (2020). Cancer-associated fibroblasts in desmoplastic tumors: emerging role of integrins. Seminars in Cancer Biology 62 166–181. 10.1016/j.semcancer.2019.08.004
  * Zhang et al. (2006) [author] Zhang, Hao HelenH. H., Ahn, JeongyounJ., Lin, XiaodongX. and Park, CheolwooC. (2006). Gene selection using support vector machines with non-convex penalty. bioinformatics 22 88–95. Publisher: Oxford University Press.
  * Zhang et al. (2016) [author] Zhang, XiangX., Wu, YichaoY., Wang, LanL. and Li, RunzeR. (2016). Variable selection for support vector machines in moderately high dimensions. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 78 53–76. Publisher: Wiley Online Library.
  * Zhao et al. (2021) [author] Zhao, JieJ., Xu, YitianY., Xu, ChangC. and Wang, TingT. (2021). A two-stage safe screening method for non-convex support vector machine with ramp loss. Knowledge-Based Systems 228 107250\. Publisher: Elsevier.