# Quadratic Unconstrained Binary Optimisation via Quantum-Inspired Annealing

Joseph Bowles  ICFO - Institut de Ciencies Fotoniques, The Barcelona Institute
of Science and Technology, Castelldefels (Barcelona) 08860, Spain  Alexandre
Dauphin  ICFO - Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, Castelldefels (Barcelona) 08860, Spain  Patrick
Huembeli  Institute of Physics, École Polytechnique Fédérale de Lausanne
(EPFL), CH-1015 Lausanne, Switzerland  José Martinez  Quside Technologies SL,
Carrer d’Esteve Terradas, 1, 08860 Castelldefels, Barcelona, Spain  Antonio
Acín  ICFO - Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, Castelldefels (Barcelona) 08860, Spain  ICREA -
Institucio Catalana de Recerca i Estudis Avançats, Lluis Companys 23, 08010
Barcelona, Spain

###### Abstract

We present a classical algorithm to find approximate solutions to instances of
quadratic unconstrained binary optimisation. The algorithm can be seen as an
analogue of quantum annealing under the restriction of a product state space,
where the dynamical evolution in quantum annealing is replaced with a
gradient-descent based method. This formulation is able to quickly find high-
quality solutions to large-scale problem instances, and can naturally be
accelerated by dedicated hardware such as graphics processing units. We
benchmark our approach for large scale problem instances with tuneable
hardness and planted solutions. We find that our algorithm offers a similar
performance to current state of the art approaches within a comparably simple
gradient-based and non-stochastic setting.

##  I Introduction

Combinatorial optimisation is a class of optimisation problems that has
applications in nearly all areas of society. Such problems involve searching
for an optimal object amongst an often enormous but finite range of potential
candidates, and are notoriously difficult to solve. One type of combinatorial
optimisation, called _quadratic unconstrained binary optimisation_ (QUBO) [1,
2] involves searching for a bitstring that minimises a quadratic function of
its elements. QUBO problems have recently attracted considerable attention,
largely because they can be naturally tackled by quantum computers by first
mapping the problem to the energy minimisation of a classical Ising model, and
then promoting this system to a quantum Ising model [3, 4, 5, 6, 7, 8, 9, 10].
By exploiting phenomena such as superposition and entanglement, the hope is
that these quantum algorithms provide faster or higher-quality solutions than
their classical counterparts. Much of the focus of showing a quantum advantage
versus classical optimisation has consequently been focused around the class
of QUBO problems; for example, the D-Wave quantum computer [11, 12]
exclusively solves this type of optimisation problem.

At the same time, the interest in QUBO problems has inspired new classical
algorithms [13, 14, 15, 16, 17, 18] and corresponding optimization devices
[19, 20, 21, 22, 23, 24]. Since these algorithms can be run on digital logic,
they can often handle orders of magnitude more variables than current quantum
computers, and as such will serve as valuable classical benchmarks as quantum
computers increase in size. Furthermore, since it is known that many hard
optimization problems can be mapped to the QUBO setting [25, 26, 27, 28, 29],
these classical solvers may also lead to improvements in classical
optimization heuristics in general. An important question is thus: what types
of classical algorithm are best suited to solving large-scale QUBO problems? A
standard approach is to use algorithms such as simulated annealing or
population annealing [30], since these algorithms are naturally suited to the
discrete nature of QUBO problems. Although they perform well in general, the
fast implementation of these algorithms for large-scale problems is limited by
hardware (for example, Fujitsu’s digital annealer [31] is currently limited to
8192 variables due to limitation of CPU cache), and it is not clear how best
to parallelize these algorthims since in their standard versions they rely on
sequential parameter updates. A more recent method involves classically
simulating the dynamical evolution of a physical system. This is generally
done using continuous degrees of freedom (such as position and momentum) whose
energy is related to the corresponding Ising Hamiltonian. Two recently
introduced methods, Toshiba’s simulated bifurcation (SB) [13, 14], and the
simulated coherent Ising machine (SIM-CIM) [18] are of this form. The original
SB algorithm [13] is based on a simulation of classical nonlinear Hamiltonian
system of oscillators, and was later developed into two modified versions of
the algorithm [14]. SIM-CIM [18] is a classical simulation of the optical
neural network called the coherent ising machine [9]. Importantly, both these
methods are suited to solving large-scale QUBO problems of tens or even
hundreds of thousands of variables, owing to the possibility of GPU
acceleration of the computational demanding part of these algorithms. Another
interesting approach is to map the QUBO problem to the optimisation of a
classical neural network. This has been recently investigated via the use of
graph neural networks [32], autoregressive neural networks [33], and neural-
network quantum states [34, 35].

In this work we explore an alternative, quantum-inspired classical algorithm
for QUBO problems. The algorithm is inspired by quantum annealing [8, 36], and
is called _local quantum annealing_ (LQA). In a standard quantum annealing
algorithm, one evolves a multi-qubit quantum state through a Schrödinger
evolution that is generated by a time-dependent Hamiltonian, where the ground
state solution of the final Hamiltonian encodes the solution to the problem.
In LQA, we use the same Hamiltonian to define a time-dependent cost function
that we iteratively optimise via a momentum-accelerated gradient descent based
approach. In order to make the optimisation tractable, the cost-function is
optimised over a subset of product quantum states that is guaranteed to
contain the problem solution. In this way, the system stays in a low energy
(product) state throughout the optimisation, which can be seen as an
approximation of the annealing process. Since we use a gradient descent
approach on the energy landscape, the method is not however equivalent to
simulating a type of adiabatic quantum evolution in which the system stays in
the global minimum at all times, as is the case for quantum annealing.
Surprisingly however, for small systems the method is often able to find the
global minimum via the route defined by the gradient descent procedure, and
for larger systems gives a method for finding good approximate solutions in
very short time.

We note that this approach is reminiscent but not equivalent to those studied
in [37], which propose to update the parameters via a dynamical physical
evolution defined by the energy of the system. We comment more on this in
section II. Similarly to SB and SIM-CIM, the computationally expensive step in
the algorithm corresponds to a matrix-vector multiplication which can be
accelerated by dedicated hardware such as GPUs and FPGAs, and so the algorithm
is naturally parallelizable and requires approximately the same resources per
optimisation step as SB and SIM-CIM. Unlike simulated bifurcation however, our
approach is purely gradient based, and unlike SIM-CIM, it does not consume
randomness during the optimisation. We note that our approach shares some
similarities with the molecular dynamics part of the hybrid quantum annealing
algorithm presented in [38], where it is possible that the recursive form of
the leap-frog algorithm used there plays a similar role to gradient descent in
our approach.

We benchmark our algorithm against the three alternative versions of the SB
algorithm [14] and SIM-CIM [18]. Since all these algorithms require the same
computational resources per optimisation step this makes for a fair and
uncontroversial comparison in terms of the solution quality per optimisation
step. We focus most of our benchmarking around recently developed methods of
planted solutions [39, 40, 41], which provide constructions of QUBO problems
with tuneable hardness whose global optimal solution is known. In our opinion
this provides a more complete picture than benchmarks that focus on a single
problem, or an arbitrary sets of problems whose hardness is unknown. We find
that LQA provides comparable solution quality to SB and SIM-CIM over all
problem instances and thus opens an alternative route to large-scale QUBO
optimisation via a purely gradient based deterministic algorithm. We also
believe that the relative simplicity of our approach paves the way for a
number of potential improvements or extensions that we discuss at the end of
the article.

An implementation of our algorithm in PyTorch is available as a Python
package; see https://github.com/josephbowles/localquantumannealing.

Figure 1: Plot of the cost function value as a function of $t$ for a fully
connected 20-spin problem with randomly chosen weights in the interval
$[-1,1]$. The dashed straight lines show the cost of the global solutions for
the times $t=0$ and $t=1$, evaluated at different values of $t$ (i.e. the
energy crossing of the two ground states if the two Hamiltonian $H_{x}$ and
$H_{z}$ would commute). The lower dashed line shows the minimum of the cost
function of the product state ansatz (8) obtained via an intensive basin-
hopping algorithm. Here, we consider a vanilla gradient descent approach
(mom=0) and a momentum-assisted approach (mom=0.99). One sees that the use of
momentum is particularly effective at helping the system leave the initial
minimum, since this point becomes a saddle point for some value of $t$. This
helps the momentum assisted approach stay closer to the global minimum
throughout the optimisation, and reach the global minimum.

##  II Methods

Formally, a QUBO optimisation problem is one of the form

$\displaystyle\min_{\mathbf{x}\in\\{0,1\\}^{n}}\quad\mathbf{x}^{T}\;Q\;\mathbf{x}+\mathbf{x}^{T}\;\mathbf{a}$
(1)

where $Q$ is a $n\times n$ real symmetric matrix such that $Q_{ii}=0\;\forall
i$, and $\mathbf{a}$ is an $n\times 1$ real vector. By defining the $\pm 1$
valued variables $s_{i}=2x_{i}-1$, and corresponding vector $\mathbf{s}$, the
problem is equivalent (up to a problem dependent constant) to

$\displaystyle\min_{\mathbf{s}\in\\{+1,-1\\}^{n}}\quad\mathbf{s}^{T}\;J\;\mathbf{s}+\mathbf{s}^{T}\;\mathbf{b}$
(2)

where $J=Q/4$, $\mathbf{b}=(\mathbf{a}+Q\mathbf{1})/2$, and $\mathbf{1}$ is
the vector of ones. Note that the second term in (2) can be incorporated into
the first at the expense of adding one extra variable with fixed value $+1$.
From hereon we therefore consider problems in the above form with
$\mathbf{b}=0$ so that our problem becomes

$\displaystyle\min_{\mathbf{s}\in\\{+1,-1\\}^{n}}\quad\mathbf{s}^{T}\;J\;\mathbf{s}.$
(3)

In this form, the optimisation is equivalent to minimising the energy of a
classical Ising Hamiltonian $J$ with no bias field, where the variables
$s_{i}$ are interpreted as classical spin values.

The minimum (3) is equivalently obtained as the ground state of the quantum
Ising Hamiltonian

$\displaystyle H_{z}=\sum_{ij}J_{ij}\sigma_{z}^{(i)}\sigma_{z}^{(j)}$ (4)

where $\sigma_{z}^{(i)}$ denotes the Pauli z matrix applied to the
$i^{\text{th}}$ qubit of an $n$ qubit system. That is, one way to find the
minimum (2) is via

$\displaystyle\min_{\ket{\psi}\in\mathbb{C}^{2^{n}}}\bra{\psi}H_{z}\ket{\psi}$
(5)

where $\ket{\psi}$ is a normalised quantum state. Since $H_{z}$ is diagonal in
the z basis, this minimum energy will be attained by one (or more) of the z
basis product states. In a quantum annealing algorithm, one attempts to solve
(5) by considering a time-dependent Hamiltonian such as

$\displaystyle H(t)=tH_{z}\gamma-(1-t)H_{x}$ (6)

where

$\displaystyle H_{x}=\sum_{i}\sigma_{x}^{(i)}.$ (7)

Here $\gamma>0$ controls the relative strength of $H_{z}$ contribution to the
energy. The quantum state of the system is initially prepared in the state
$\ket{+}^{\otimes n}$, which is the ground state of $H(0)$, and the system is
evolved by changing the Hamiltonian from $t=0$ to $t=1$. If this evolution is
done slowly enough, the adiabatic theorem guarantees that the state will stay
in the ground state throughout the evolution and the global solution to the
problem will therefore be obtained.

Our algorithm is inspired from the quantum annealing approach. Since a
classical simulation of quantum annealing is likely impossible due to the
exponential memory requirements needed to store the quantum state vector, we
do not consider the set of all quantum states, but restrict ourselves to
product states of the form

$\displaystyle\ket{\boldsymbol{\theta}}=\ket{\theta_{1}}\otimes\ket{\theta_{2}}\otimes\cdots\otimes\ket{\theta_{n}}$
(8)

where

$\displaystyle\ket{\theta_{i}}=\cos\frac{\theta_{i}}{2}\ket{+}+\sin\frac{\theta_{i}}{2}\ket{-}.$
(9)

We therefore have

$\displaystyle\bra{\theta_{i}}\sigma_{z}\ket{\theta_{i}}=\sin\theta_{i}\;,\quad\quad\bra{\theta_{i}}\sigma_{x}\ket{\theta_{i}}=\cos\theta_{i}.$
(10)

Considering these states, the minimisation (5) is equivalent to

$\displaystyle\min_{\boldsymbol{\theta}}\bra{\boldsymbol{\theta}}H_{z}\ket{\boldsymbol{\theta}}$
$\displaystyle=\min_{\boldsymbol{\theta}}\sum_{ij}J_{ij}\sin\theta_{i}\sin\theta_{j}$
(11)

In analogy with the quantum annealing algorithm we define a time-dependent
cost function corresponding to the energy of the system at time $t$:

$\displaystyle\mathcal{C}(t,\boldsymbol{\theta})$
$\displaystyle=\bra{\boldsymbol{\theta}}tH_{z}\gamma-(1-t)H_{x}\ket{\boldsymbol{\theta}}$
$\displaystyle=t\gamma\sum_{ij}J_{ij}\sin\theta_{i}\sin\theta_{j}-(1-t)\sum_{i}\cos\theta_{i}$
(12)

In our algorithm, we further parameterise the values $\theta_{i}$ as
$\theta_{i}=\frac{\pi}{2}\tanh{(w_{i})}$ where $w_{i}\in\mathbb{R}$ so that
$w_{i}\rightarrow\pm\infty\implies\ket{\theta_{i}}\rightarrow\ket{0},\ket{1}$.
The cost (II) therefore becomes

$\displaystyle\mathcal{C}(t,\mathbf{w})=t\gamma\;\mathbf{z}^{T}J\;\mathbf{z}\;-(1-t)\mathbf{x}^{T}\cdot\mathbf{1},$
(13)

with $\mathbf{z}=(\sin(\frac{\pi}{2}\tanh
w_{1}),\cdots,\sin(\frac{\pi}{2}\tanh w_{n}))$,
$\mathbf{x}=(\cos(\frac{\pi}{2}\tanh w_{1}),\cdots,\cos(\frac{\pi}{2}\tanh
w_{n}))$. The gradient of (13) with respect to the parameters $\mathbf{w}$ is

$\displaystyle\nabla_{\mathbf{w}}\;\mathcal{C}(\mathbf{w},t)=\frac{\pi}{2}[t\gamma(2J\mathbf{z})\circ\mathbf{x}+(1-t)\mathbf{z}]\circ
a(\mathbf{w}),$ (14)

where $\circ$ denotes element-wise multiplication of vectors and
$a(\cdot)=1-\tanh^{2}(\cdot)$ is the derivative of the $\tanh$ function and
$a(\cdot)$ acts element-wise.

Figure 2: Benchmarking using the $K_{2000}$ Max-Cut problem introduced in
[18]. This problem is equivalent to solving a 2000 spin fully connected Ising
problem with $J_{ij}=\pm 1$ chosen uniformly random. Here, we compare LQA to
the original simulated bifurcation algorithm [13] and the simulated coherent
Ising mahcine [18] (top panel), and to the more recent ballistic and discrete
adaptations of the simulated bifurcation algorithm [14] (right). We perform
100 trials for each algorithm, each with 5000 optimisation steps. The thick
lines show the average values over the 100 trials as a function of the
optimisation step. The shaded regions correspond to the max and min cut values
over all trials. The histograms plot the final obtained cut values at the end
of the optimisation over the 100 trials. For LQA, we chose an initial step
size of 1, and set $\gamma=0.1$. For SB, SBB and SBD the step size is set to
0.5, 1.25, 1.25 respectively, with other parameters being set as recommended
in [14]. For SIM-CIM, a step size of 1 was chosen and the noise parameter
$A_{n}$ set to 0.25. In all algorithms, the initial parameters were set as
0.1*rand, where rand is a random list with elements in the range [-1,1].

Figure 3: Top panels: Benchmark results for the tile planting problems defined
in [39], for 11 problem instances (x-axis) with increasing hardness. The
number of spins in each problem is 1024. For each instance, 500 trials with
500 steps are performed. The solid lines plot the average relative error to
the known global minimum (defined as $|\frac{C-C_{0}}{C_{0}}|$ where $C_{0}$
is the global minimum and $C$ is the value obtained in a given trial) over the
500 trials. The upper and lower coloured markers denote the maximum and
minimum values obtained over the 500 trials. The shaded regions correspond to
1 standard deviation around the mean. Lower panel: Analogous benchmarking
results for the Wishart problems defined in [40] which feature and easy-hard-
easy problem hardness transition (note that the global minima is obtained for
the last instance for some algorithms). Here, 10 problem instance were chosen
and 500 trials were performed each with 500 optimisation steps. The initial
step size for LQA was set to $2$ for both benchmarks. For SIM-CIM the step
size was set to $1$ and $A_{n}$ to 0.25 for both benchmarks. For tile
planting, the step sizes for SB, SBB and SBD were set to $0.5$, $1.25$,
$1.25$. For the Wishart planting, the step sizes for SB, SBB and SBD depend on
the problem instance due to instabilities in the optimisaiton for the first
two problem instances; see appendix A for a table with the precise stepsizes.
The parameter initilisation were set as in figure 2.

In order to obtain an approximate solution to the QUBO problem, one may use a
gradient descent routine on $\mathcal{C}(t,\mathbf{w})$. More precisely, in
each step $i$ of the algorithm, one updates the parameters $\mathbf{w}$
according to the gradient of $\mathcal{C}(t,\mathbf{w})$, where $t=i/N$ and
$N$ is the total number of steps. A spin configuration can be obtained at any
time via $\mathbf{s}=\text{sign}(\mathbf{w})$. Here, one may chose from a
plethora of strategies that exist for gradient-based optimisation, such as
momentum-assistance and adaptive step sizes. For example, the standard
momentum-assisted gradient descent technique uses an additional velocity
vector $\boldsymbol{\nu}$ (typically initialised as the zero vector) and a
momentum parameter $\mu\in[0,1]$, and corresponds to the parameter update (at
time t)

$\displaystyle\boldsymbol{\nu}\leftarrow\mu\boldsymbol{\nu}-\eta\nabla\mathcal{C}_{\mathbf{w}}(\mathbf{w},t)$
(15) $\displaystyle\mathbf{w}\leftarrow\mathbf{w}+\boldsymbol{\nu},$ (16)

where $\eta$ is the step-size of the gradient descent. Momentum can help to
accelerate the gradient descent in flat regions of optimisation, and is widely
used in the training of machine learning models. Another common technique is
the Adam update [42], which uses the momentum technique and further adapts the
step size after each update. In the following pseudocode we denote a generic
update as $\mathbf{w}\leftarrow
g(\mathbf{w},\nabla\mathcal{C}_{\mathbf{w}}(\mathbf{w},t))$, where $g$ may
incorporate additional information such as momentum and step size.

Data: $J\in\mathbb{R}^{n\times n}$: symmetric Ising matrix;
$\mathbf{w}_{0}\in\mathbb{R}^{n}$: initial weights; $N$: total steps

initialization $\mathbf{w}=\mathbf{w}_{0}$;

for $i=1,\cdots,N$ do

$\mathbf{w}\leftarrow
g(\mathbf{w},\nabla\mathcal{C}_{\mathbf{w}}(\mathbf{w},i/N))$

end for

return sign$(\mathbf{w})$

Algorithm 1 Local quantum annealing

The effect of this algorithm is to continuously push the parameters toward a
local minimum of the time-evolving cost function. One may therefore hope to
mimic the evolution of a quantum annealing algorithm, in which the system
stays in the instantaneous ground state of the time-dependent Hamiltonian
throughout the optimisation. For this reason, it is best to initialise the
parameters close to zero, since the ground state of
$\mathcal{C}(\mathbf{w},0)$ is given by $\ket{+}^{\otimes n}$.

A couple of important points are in order here. First, note that the
computationally expensive part of the algorithm is the matrix multiplication
$J\mathbf{z}$ that appears in the gradient calculation (14). This can be
accelerated by dedicated hardware such as GPUs or FPGAs, which for large
problems results in a tremendous performance boost with respect to using a
CPU. Secondly, we perform only a single parameter update for each value of
$t$, rather than waiting for convergence to a local minima before stepping
$t$, which requires a much larger optimisation time and typically leads to
similar results. Even for large optimisation times, the method is not
guaranteed to converge to the globally optimal solution, since the
optimisation is performed over a much smaller space of product quantum states,
and there is no corresponding guarantee that following a locally optimal
minimum throughout the optimisation will lead to the ground state of the final
Hamiltonian (due to e.g. a first order phase transition). Nevertheless, one
may expect that the behaviour approximates the quantum annealing behaviour to
some extent.

In figure 1 we plot the evolution of the cost, given by
$\mathcal{C}(\mathbf{w},t)$ for a simple problem involving 20 spins. Here we
also investigate the use of momentum-assisted gradient descent. Although the
system does not stay in the global minimum throughout the optimisation (lower
dashed curve), the momentum-assisted approach eventually finds the ground
state of the system via a different path. The pure gradient-descent algorithm
has difficultly leaving its initial parameters $\ket{+}^{\otimes{n}}$, and
more frequently does not converge to the global solution. This can be
attributed to the fact that the initial local minimum becomes a saddle point
from which it is very slow to escape without momentum assistance. We therefore
suspect that momentum-assistance is vital in achieving good performance, which
we have found to generally be the case.

##  III Results

Here we present a number of benchmarking results. We compare our algorithm
against the three variants of simulated bifurcation (the original algorithm
(SB) [13], ballistic simulated bifurcation (SBB) [14], and discrete simulated
bifurcation (SBD) [14]), and the simulated coherent Ising machine (SIM-CIM)
[18]. All algorithms were implemented in pytorch on a standard laptop CPU. For
LQA, the Adam gradient descent method [42] was adopted for the parameter
updates.

We focus on three types of problem. The first is the $K_{2000}$ Max-Cut
problem (see figure 2). This is a benchmark introduced in [9] that has been
tested on both simulated bifurcation and SIM-CIM which exploits the Max-Cut
problem mapping to QUBO (see e.g. [43]). Here, larger values correspond to
higher quality solutions. We find that LQA achieves the highest mean value
over 100 optimisation trials of 5000 steps each, with SBD achieving the best
value over all trials. The distribution of final values varies significantly
over the algorithms, and LQA and SBB both feature a small variance with a
strong peak in a single solution. We note that this behaviour is sensitive to
hyperparameter choice (one can achieve a larger variance at the cost of a
lower mean value by adjusting the initial step size) and thus seems not to be
a generic feature of LQA.

Next, we consider performance with respect to problems with planted solutions
(see figure 3). To generate these problems we make use of the Chook package
[41]. We consider two such classes: the first are the 2D ‘tile planting’
problems (figure 3 top panels) introduced in [39]. These problems are defined
on a 2D lattice and although in principle solvable in polynomial time due to
their planar nature, are often challenging for heuristic solvers. Here we
consider 1024 spin problems that are constructed from the C2 and C3 tiles (see
[39] for more details). We generate 11 problem instances of this type, where
the probability to use a C3 tile is increased linearly from 0 to 1 over the
problem instances. It has been observed that this results in progressively
harder problems, which is reproduced in our results. For these problems, the
best results were obtained by SBB, with LQA giving similar results to SB.

We also study the class of Wishart planted solutions (figure 3 lower panels)
defined in [40], which generate fully connected problems with an easy-hard-
easy transition. Here we consider 500 spins problems. We find a very similar
performance among all algorithms, with SBB and LQA performing the best in the
hard central problem region. SBD performed almost identically to SIM-CIM so we
do not show these results here for clarity. SB has some problems with
stability for the initial problem instance which we suspect could be remedied
with appropriate hyper-parameter tuning. For problem instances at the latter
easy tail of the sequence, it is known that the probability for local
optimisation methods to find the ground state solution by chance increases
significantly [40]; the relatively few trials in which the algorithms find the
global minimum in these cases may thus be more an indication of luck for this
particular batch of trials rather than a feature of the algorithm itself.

##  IV Discussion

Our method is reminiscent of but different to the approaches suggested in [37,
44], which use a similar product state ansatz and Hamiltonian to mimik the
effects of quantum annealing. These approaches are based on a dynamical
(physical) evolution of the parameters $\boldsymbol{\theta}$ under the
Hamiltonian (6). In its most basic form, our method uses a simple gradient
descent update, which results in a different (un-physical) trajectory of the
parameters $\boldsymbol{\theta}$. Under the momentum-assisted update (15) and
in the continuous-time limit (i.e. for infinitesimal step size), it is known
that a physical interpretation of our approach is possible [45]: namely, as a
dynamical evolution of a system in a viscous medium in which the momentum
parameter plays the role of mass. This evolution is not equivalent to that of
[37, 44] due to the effects of damping implied by the viscosity. Furthermore,
since our parameter updates are done at the level of the variables
$\mathbf{w}$ and not $\boldsymbol{\theta}$, the variables on which this
physical evolution is understood are not the same. We have found that these
differences can result in significant differences in solution quality for
large problems. For the case of more complex parameter updates, such as the
Adam [42] update that we use for our benchmarking, it is not clear if the
parameter evolution can be understood physically. In any case we note that, as
is typical with gradient descent methods, the solution quality can be quite
sensitive to the choice of step size, and it is the case that large step sizes
often outperform small step sizes, even for long optimisation times. Given the
considerations, our method could also be be viewed as a type of gradient-
descent based graduated optimisation [46, 47, 48, 49] (also called the
‘continuation method’).

The performance of the tested algorithms is similar throughout all benchmarks,
with no algorithm clearly outperforming another. This is perhaps not
surprising, since although they are all designed from quite different starting
points, the parameter updates all make use of the same matrix-vector product
between the coupling matrix $J$ and a parameter dependent vector that encodes
the solution, which becomes more dominant as the optimisation progresses. For
LQA, this is the product $J\mathbf{z}$ in equation (14). It is therefore
unclear whether one should expect any large difference in performance between
any of these algorithms since they may all feature basins of attraction to
similar solution qualities despite their differing parameter updates. We would
argue however, that LQA may be the most versatile of the options. Firstly, it
is pure gradient-descent based; it can therefore make use of myriad of tools
from machine learning for gradient based optimisation, as well as be used as a
subroutine in any other gradient-based algorithm. The use of Adam in our tests
makes it quite stable to initial step size variation, which to some extent
removes the burden of setting this hyperparameter. Although not done here,
second order derivatives of the cost could be calculated analytically via
(14). Thus, methods that make use of second order information could be
incorporated exactly and may give a performance boost.

We believe that the main value of the algorithm however is the form of the
time-dependent cost function, since this results in significantly better
solutions than using a static Hamiltonian. It would be interesting to
investigate if this approach could be improved. For example, although we use a
simple linear annealing schedule here, non-linear schedules may give better
results. On this note, there have been a number of works which use alternative
transverse Hamiltonians in order to improve the performance of quantum
annealing [50, 51]. It would be interesting to investigate if the different
cost landscapes implied by these Hamiltonians could lead to improvements.
Finally, it would be interesting to see if the product ansatz (8) could be
expanded to include a wider range of quantum states without significantly
sacrificing the speed of the algorithm. Here, previous works connecting QUBO
optimisation to neural network quantum states [34, 52] and matrix product
state [53] may be valuable.

_Related work_ —We note that while preparing this draft we became aware of an
independent work [54] that also suggests using a gradient approach that is
similar to ours.

_Acknowledgments_ — This work was supported by the ICFO-Quside Joint
Laboratory in Quantum Processing, Fundacio Cellex, Fundacio Mir-Puig,
Generalitat de Catalunya (SGR 1381, SGR 1341, QuantumCAT, and CERCA
Programme), ERC AdG NOQIA and CERQUTE, Spanish MINECO (Severo Ochoa
CEX2019-000910-S, Plan National FIDEUA and TRANQI, Retos QuSpin, FPI, QUANTERA
MAQS PCI2019-111828-2 / 10.13039/501100011033), EU Horizon 2020 FET-OPEN
OPTOLogic (Grant No 899794), and the National Science Centre, Poland (Symfonia
Grant No. 2016/20/W/ST4/00314), Marie Skłodowska-Curie grant STRETCH No
101029393, a fellowship granted by la Caixa Foundation (ID 100010434,
fellowship code LCF/BQ/PR20/11770012), and the AXA Chair in Quantum
Information Science.

## References

  * Boros and Hammer [1991] E. Boros and P. L. Hammer, The max-cut problem and quadratic 0–1 optimization; polyhedral aspects, relaxations and bounds, Annals of Operations Research 33, 151 (1991).
  * Glover  _et al._ [2019] F. Glover, G. Kochenberger, and Y. Du, Quantum bridge analytics i: a tutorial on formulating and using qubo models, 4OR 17, 335 (2019).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv preprint arXiv:1411.4028 (2014).
  * Date  _et al._ [2019] P. Date, R. Patton, C. Schuman, and T. Potok, Efficiently embedding qubo problems on adiabatic quantum computers, Quantum Information Processing 18, 1 (2019).
  * Jünger  _et al._ [2019] M. Jünger, E. Lobe, P. Mutzel, G. Reinelt, F. Rendl, G. Rinaldi, and T. Stollenwerk, Performance of a quantum annealer for ising ground state computations on chimera graphs, arXiv preprint arXiv:1904.11965 (2019).
  * Lee  _et al._ [2021] J. Lee, A. B. Magann, H. A. Rabitz, and C. Arenz, Towards favorable landscapes in quantum combinatorial optimization, arXiv preprint arXiv:2105.01114 (2021).
  * Venturelli  _et al._ [2015] D. Venturelli, S. Mandrà, S. Knysh, B. O’Gorman, R. Biswas, and V. Smelyanskiy, Quantum optimization of fully connected spin glasses, Physical Review X 5, 031040 (2015).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse ising model, Physical Review E 58, 5355 (1998).
  * Inagaki  _et al._ [2016] T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, _et al._ , A coherent ising machine for 2000-node optimization problems, Science 354, 603 (2016).
  * Boettcher [2019] S. Boettcher, Analysis of the relation between quadratic unconstrained binary optimization and the spin-glass ground-state problem, Physical Review Research 1, 033142 (2019).
  * [11] D-wave systems documentation, https://docs.dwavesys.com/docs/latest/index.html.
  * Boothby  _et al._ [2020] K. Boothby, P. Bunyk, J. Raymond, and A. Roy, Next-generation topology of d-wave quantum processors, arXiv preprint arXiv:2003.00133 (2020).
  * Goto  _et al._ [2019] H. Goto, K. Tatsumura, and A. R. Dixon, Combinatorial optimization by simulating adiabatic bifurcations in nonlinear hamiltonian systems, Science advances 5, eaav2372 (2019).
  * Goto  _et al._ [2021] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura, High-performance combinatorial optimization based on classical mechanics, Science Advances 7, eabe7953 (2021).
  * Okuyama  _et al._ [2019] T. Okuyama, T. Sonobe, K.-i. Kawarabayashi, and M. Yamaoka, Binary optimization by momentum annealing, Physical Review E 100, 012111 (2019).
  * Patel  _et al._ [2020] S. Patel, L. Chen, P. Canoza, and S. Salahuddin, Ising model optimization problems on a fpga accelerated restricted boltzmann machine, arXiv preprint arXiv:2008.04436 (2020).
  * Haribara  _et al._ [2017] Y. Haribara, H. Ishikawa, S. Utsunomiya, K. Aihara, and Y. Yamamoto, Performance evaluation of coherent ising machines against classical neural networks, Quantum Science and Technology 2, 044002 (2017).
  * Tiunov  _et al._ [2019] E. S. Tiunov, A. E. Ulanov, and A. Lvovsky, Annealing by simulating the coherent ising machine, Optics express 27, 10288 (2019).
  * Yamaoka  _et al._ [2015] M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, A 20k-spin ising chip to solve combinatorial optimization problems with cmos annealing, IEEE Journal of Solid-State Circuits 51, 303 (2015).
  * Yamamoto  _et al._ [2017] K. Yamamoto, W. Huang, S. Takamaeda-Yamazaki, M. Ikebe, T. Asai, and M. Motomura, A time-division multiplexing ising machine on fpgas, in  _Proceedings of the 8th International Symposium on Highly Efficient Accelerators and Reconfigurable Technologies_ (2017) pp. 1–6.
  * Kihara  _et al._ [2017] Y. Kihara, M. Ito, T. Saito, M. Shiomura, S. Sakai, and J. Shirakashi, A new computing architecture using ising spin model implemented on fpga for solving combinatorial optimization problems, in  _2017 IEEE 17th International Conference on Nanotechnology (IEEE-NANO)_ (IEEE, 2017) pp. 256–258.
  * Tsukamoto  _et al._ [2017] S. Tsukamoto, M. Takatsu, S. Matsubara, and H. Tamura, An accelerator architecture for combinatorial optimization problems, Fujitsu Sci. Tech. J 53, 8 (2017).
  * Yamamoto  _et al._ [2020] K. Yamamoto, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, H. Teramoto, A. Sakai, S. Takamaeda-Yamazaki, and M. Motomura, 7.3 statica: A 512-spin 0.25 m-weight full-digital annealing processor with a near-memory all-spin-updates-at-once architecture for combinatorial optimization with complete spin-spin interactions, in  _2020 IEEE International Solid-State Circuits Conference-(ISSCC)_ (IEEE, 2020) pp. 138–140.
  * Kalinin  _et al._ [2020] K. P. Kalinin, A. Amo, J. Bloch, and N. G. Berloff, Polaritonic xy-ising machine, Nanophotonics 9, 4127 (2020).
  * Bauckhage  _et al._ [2017] C. Bauckhage, E. Brito, K. Cvejoski, C. Ojeda, R. Sifa, and S. Wrobel, Ising models for binary clustering via adiabatic quantum computing, in  _International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition_ (Springer, 2017) pp. 3–17.
  * Date  _et al._ [2021] P. Date, D. Arthur, and L. Pusey-Nazzaro, Qubo formulations for training machine learning models, Scientific Reports 11, 1 (2021).
  * Calude  _et al._ [2017] C. S. Calude, M. J. Dinneen, and R. Hua, Qubo formulations for the graph isomorphism problem and related problems, Theoretical Computer Science 701, 54 (2017).
  * Bauckhage  _et al._ [2019] C. Bauckhage, N. Piatkowski, R. Sifa, D. Hecker, and S. Wrobel, A qubo formulation of the k-medoids problem., in  _LWDA_ (2019) pp. 54–63.
  * Papalitsas  _et al._ [2019] C. Papalitsas, T. Andronikos, K. Giannakis, G. Theocharopoulou, and S. Fanarioti, A qubo model for the traveling salesman problem with time windows, Algorithms 12, 224 (2019).
  * Wang  _et al._ [2015] W. Wang, J. Machta, and H. G. Katzgraber, Comparing monte carlo methods for finding ground states of ising spin glasses: Population annealing, simulated annealing, and parallel tempering, Physical Review E 92, 013303 (2015).
  * Aramon  _et al._ [2019] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. G. Katzgraber, Physics-inspired optimization for quadratic unconstrained problems using a digital annealer, Frontiers in Physics 7, 48 (2019).
  * Schuetz  _et al._ [2021] M. J. Schuetz, J. K. Brubaker, and H. G. Katzgraber, Combinatorial optimization with physics-inspired graph neural networks, arXiv preprint arXiv:2107.01188 (2021).
  * Hibat-Allah  _et al._ [2021] M. Hibat-Allah, E. M. Inack, R. Wiersema, R. G. Melko, and J. Carrasquilla, Variational neural annealing, arXiv preprint arXiv:2101.10154 (2021).
  * Gomes  _et al._ [2019] J. Gomes, K. A. McKiernan, P. Eastman, and V. S. Pande, Classical quantum optimization with neural network quantum states, arXiv preprint arXiv:1910.10675 (2019).
  * Zhao  _et al._ [2020] T. Zhao, G. Carleo, J. Stokes, and S. Veerapaneni, Natural evolution strategies and variational monte carlo, Machine Learning: Science and Technology 2, 02LT01 (2020).
  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem, Science 292, 472 (2001).
  * Smolin and Smith [2014] J. A. Smolin and G. Smith, Classical signature of quantum annealing, Frontiers in physics 2, 52 (2014).
  * Irie  _et al._ [2021] H. Irie, H. Liang, T. Doi, S. Gongyo, and T. Hatsuda, Hybrid quantum annealing via molecular dynamics, Scientific reports 11, 1 (2021).
  * Perera  _et al._ [2020a] D. Perera, F. Hamze, J. Raymond, M. Weigel, and H. G. Katzgraber, Computational hardness of spin-glass problems with tile-planted solutions, Physical Review E 101, 023316 (2020a).
  * Hamze  _et al._ [2020] F. Hamze, J. Raymond, C. A. Pattison, K. Biswas, and H. G. Katzgraber, Wishart planted ensemble: A tunably rugged pairwise ising model with a first-order phase transition, Physical Review E 101, 052102 (2020).
  * Perera  _et al._ [2020b] D. Perera, I. Akpabio, F. Hamze, S. Mandra, N. Rose, M. Aramon, and H. G. Katzgraber, Chook–a comprehensive suite for generating binary optimization problems with planted solutions, arXiv preprint arXiv:2005.14344 (2020b).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Lucas [2014] A. Lucas, Ising formulations of many np problems, Frontiers in physics 2, 5 (2014).
  * Hatomura and Mori [2018] T. Hatomura and T. Mori, Shortcuts to adiabatic classical spin dynamics mimicking quantum annealing, Physical Review E 98, 032136 (2018).
  * Qian [1999] N. Qian, On the momentum term in gradient descent learning algorithms, Neural networks 12, 145 (1999).
  * Hazan  _et al._ [2016] E. Hazan, K. Y. Levy, and S. Shalev-Shwartz, On graduated optimization for stochastic non-convex problems, in  _International conference on machine learning_ (PMLR, 2016) pp. 1833–1841.
  * Gargiani  _et al._ [2020] M. Gargiani, A. Zanelli, Q. Tran-Dinh, M. Diehl, and F. Hutter, Convergence analysis of homotopy-sgd for non-convex optimization, arXiv preprint arXiv:2011.10298 (2020).
  * Mobahi and Fisher III [2015] H. Mobahi and J. Fisher III, A theoretical analysis of optimization by gaussian continuation, in  _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 29 (2015).
  * Liu and Qiao [2013] Z.-Y. Liu and H. Qiao, Gnccp—graduated nonconvexityand concavity procedure, IEEE transactions on pattern analysis and machine intelligence 36, 1258 (2013).
  * Nishimori and Takada [2017] H. Nishimori and K. Takada, Exponential enhancement of the efficiency of quantum annealing by non-stoquastic hamiltonians, Frontiers in ICT 4, 2 (2017).
  * Susa  _et al._ [2018] Y. Susa, Y. Yamashiro, M. Yamamoto, and H. Nishimori, Exponential speedup of quantum annealing by inhomogeneous driving of the transverse field, Journal of the Physical Society of Japan 87, 023002 (2018).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Bauer  _et al._ [2015] B. Bauer, L. Wang, I. Pizorn, and M. Troyer, Entanglement as a resource in adiabatic quantum optimization, arXiv:1501.06914 (2015).
  * Veszeli and Vattay [2021] M. T. Veszeli and G. Vattay, Mean field approximation for solving qubo problems, arXiv preprint arXiv:2106.03238 (2021).

##  Appendix A parameter setting for Wishart planting benchmark

The table below displays the step sizes used for the Wishart planting
benchmarking of Figure 3 for each algorithm. Step sizes were set close to the
largest values possible before instabilities deteriorate the solution quality.
instance: 0 1 2 3 4 5 6 7 8 9 LQA 2 2 2 2 2 2 2 2 2 2 SB 0.5 1 1 1 1 1 1 1 1 1
SBB 1.0 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 SBD 0.1 0.1 1.25 1.25
1.25 1.25 1.25 1.25 1.25 1.25 SIMCIM 1 1 1 1 1 1 1 1 1 1


# PDBL: Improving Histopathological Tissue Classification with Plug-and-Play
Pyramidal Deep-Broad Learning

Jiatai Lin†, Guoqiang Han†, Xipeng Pan†, Hao Chen, Danyi Li, Xiping Jia,
Zhenwei Shi, Zhizhen Wang, Yanfen Cui, Haiming Li, Changhong Liang, Li Liang,
Zaiyi Liu, Chu Han This work was supported by the Key R&D Program of Guangdong
Province, China (No. 2021B0101420006), the National Science Fund for
Distinguished Young Scholars (No.81925023), the National Natural Science
Foundation of China (No. 62102103, 81771912, 82071892, 82001789, 81901704,
81702322, 81472712, 81772918, 81972277 and 62002082), High-level Hospital
Construction Project (No. DFJH201805 and DFJH201914), Natural Science
Foundation of Guangdong Province (No. 2017A030313896), Guanzhou Science &
Technology Project (No. 201804010319). (Corresponding author: Li Liang, Zaiyi
Liu, Chu Han.)Jiatai Lin and Guoqiang Han are with the School of Computer
Science and Engineering, South China University of Technology, Guangzhou,
China.Chu Han, Xipeng Pan, Zhenwei Shi, Changhong Liang and Zaiyi Liu are with
the Department of Radiology, Guangdong Provincial People’s Hospital, Guangdong
Academy of Medical Sciences, Guangzhou, Guangdong, 510080, China. E-mail:
<EMAIL_ADDRESS>Chen is with the Department of Computer Science and
Engineering, The Hong Kong University of Science and Technology, Clear Water
Bay, Hong Kong.Xiping Jia is with Guangdong Polytechnic Normal University,
School of Computer Science, Guangzhou, China, 510665.Li Liang and Danyi Li are
with Department of Pathology, Nanfang Hospital and Basic Medical College,
Southern Medical University, Guangzhou 510515, Guangdong, China. Guangdong
Province Key Laboratory of Molecular Tumor Pathology, Guangzhou 510515,
Guangdong, China.Haiming Li is with the Department of Radiology, Fudan
University Shanghai Cancer Center, Shanghai 200032, China and the Department
of Oncology, Shanghai Medical College, Fudan University, Shanghai 200032,
China.Yanfen Cui is with the Department of Radiology, Shanxi Province Cancer
Hospital, Taiyuan, Shanxi, 030013, China.†The first three authors contribute
equally.

###### Abstract

Histopathological tissue classification is a fundamental task in pathomics
cancer research. Precisely differentiating different tissue types is a benefit
for the downstream researches, like cancer diagnosis, prognosis and etc.
Existing works mostly leverage the popular classification backbones in
computer vision to achieve histopathological tissue classification. In this
paper, we proposed a super lightweight plug-and-play module, named Pyramidal
Deep-Broad Learning (PDBL), for any well-trained classification backbone to
further improve the classification performance without a re-training burden.
We mimic how pathologists observe pathology slides in different magnifications
and construct an image pyramid for the input image in order to obtain the
pyramidal contextual information. For each level in the pyramid, we extract
the multi-scale deep-broad features by our proposed Deep-Broad block (DB-
block). We equipped PDBL in three popular classification backbones,
ShuffLeNetV2, EfficientNetb0, and ResNet50 to evaluate the effectiveness and
efficiency of our proposed module on two datasets (Kather Multiclass Dataset
and the LC25000 Dataset). Experimental results demonstrate the proposed PDBL
can steadily improve the tissue-level classification performance for any CNN
backbones, especially for the lightweight models when given a small among of
training samples (less than 10%), which greatly saves the computational time
and annotation efforts.

###### Index Terms:

Pyramidal deep-broad learning, Histopathological tissue classification, Broad
learning system

##  I Introduction

Histopathological slides not only play a vital role in cancer diagnosis, but
also deliver valuable tumor microenvironment information for cancer research
[1, 2]. To analyze the whole slide images by computer algorithms is crucial
for precision medicine on cancers, such as diagnosis prediction [3, 4],
molecular status prediction [5, 6] and even the origins of the unknown primary
cancers prediction [7]. Segmenting and recognizing various tissue types is the
very first step of histopathological image analysis. Semantic segmentation [8]
is the best way to define tissue types for every single pixel. However, due to
the gigapixel resolution and the expertise requirement, obtaining pixel-level
annotations is extremely difficult and time-consuming [9]. Therefore, patch-
level classification now becomes an alternative solution [10, 11], which can
greatly save the annotation efforts.

Convolutional neural network (CNN) has demonstrated outstanding performance in
image classification problem [12], with a series of classification backbones,
i.e., ResNet [13], ShuffLeNet [14], Efficient-Net [15] and etc. They have been
rapidly extended to medical image classification [16, 17], including
histopathological image classification [18, 19]. Typically, Han et al. [20]
presented a CNN-based multi-classification model for histopathological tissue
classification of breast cancer, which achieves over 94% patch-level accuracy
under four magnification factors. Tsai et al. [21] tested five common
classification backbones on colorectal tissue classification using 100,000
training image patches. All the backbones can achieve over 95% accuracy. The
above researches prove that the current CNN classification backbones have
already demonstrated strong feature representation ability and achieved
promising results for histopathology tissue classification. In this paper, we
reconsider how to make good use of the features extracted from the existing
CNN backbones, to further improve the classification performance as well as to
increase the model generalizability, adaptability, and robustness.

In clinical practice, pathologists read histopathological slides by switching
the object lens to observe the slides under different magnifications.
Therefore, considering multi-scale contextual information is critical for
histopathological image analysis. Inspired by this observation, we proposed a
lightweight plug-and-play module for any CNN classification backbones, named
Pyramidal Deep-Broad Learning (PDBL). An image pyramid is constructed to
extract the pyramidal contextual information. For each level in the pyramid,
we proposed a Deep-Broad block (DB-block) to fully discover the multi-scale
deep-broad features extracted by the CNN backbones from low level to high
level. Our proposed PDBL can be plugged on any classification backbone and
effectively improves the classification performance with very few extra
computational resources.

We tested PDBL on three representative CNN backbones, ShuffLeNetV2,
EfficientNetb0, and ResNet50 on Kather Multiclass Dataset [11] and Lung Colon
Cancer Histopathological Image Dataset [22]. We conduct two main experiments,
one is the effectiveness of PDBL with different proportions of the training
samples, the other is the robustness of PDBL with only $1\%$ training samples.
Experimental results demonstrate that PDBL effectively improves the
classification performance on both datasets. When very limited training
samples are involved (1% only), PDBL can maintain a standout improvement
compared with the baseline models without PDBL. It can significantly reduce
the annotation efforts and computational resources. Experimental results also
show that PDBL improves domain adaptation abilities for CNN models. The
contributions of this paper are summarized as follows:

  * $\bullet$

We proposed a lightweight plug-and-play module (PDBL), which can be easily
applied on almost any common CNN-based classification backbone. It can
generally improve all the three CNN backbones we have tested for
histopathological tissue classification with no re-training burden.

  * $\bullet$

We proposed a Deep-Broad block to fully discover the multi-scale deep-broad
features from low level to high level.

  * $\bullet$

The proposed PDBL demonstrates outstanding improvement of the performance for
the lightweight models with very limited training samples (1% only).

  * $\bullet$

Models with PDBL can relieve the requirement of large scale training data and
be easily and efficiently adapted to a new domain with only a few training
samples, which greatly saves the computational resources and annotation
efforts.

##  II Related Works

###  II-A Histopathological Image Classification

Automated analysis of whole slide images (WSIs) plays a crucial role in
computer-assisted tumor diagnosis [3, 23]. Due to the giga-pixel resolution,
directly processing the entire WSI is not feasible. Moreover, obtaining pixel-
level annotation is extremely difficult. Hence, histopathological tissue
classification has been widely employed as an alternative solution for tissue
semantic segmentation of WSIs [11].

With the development of CNN models, most of the histopathological image
classification models [24, 25] are originated from the popular classification
backbones from the natural image classification. However, histopathological
image classification faces different challenges, such as extremely large image
resolution, deficiency of labels and multi-scale information integration [26].
WSI-Net [27] model was proposed to add an additional classification branch to
discard the normal tissue in order to save computational resources. Raczkowski
et al. [28] proposed a pathologist-in-the-loop model to solve the insufficient
labeling problem. Xue et al. [10] proposed to synthesis histopathological
patch images using GAN to enhance the feature representation and improve the
classification performance. Many studies attempted to extract multi-scale
features to better solve the classification problem of histopathological
images with end-to-end deep learning models, such as the Deep-Hipo model [29]
and multi-resolution model [30].

In this paper, we used a broad learning strategy to fully discover the deep
features extracted by deep learning models and leverage the multi-scale
contextual information to improve the performance of the CNN-based models
without excessive computational costs. Next, we will introduce some common
deep learning architectures and broad learning system approaches.

###  II-B Deep Learning Architectures

Deep learning models have already dominated the image classification problem
[31]. They usually go through several stages to reduce the feature dimension
and to extract higher-level semantic features, defined in Eq. 1.

$\mathcal{F}=\mathbf{stage}_{1}\odot\mathbf{stage}_{2}\odot...\odot\mathbf{stage}_{j}(X),j=1,2,...,h.$
(1)

where $\mathcal{F}$ denotes the CNN model and $X$ is the input image. Each
stage is composed of a series of cascaded convolutional blocks, such as Res-
block [13], Efficient-block [15], Shuffle-block [32] and Inception-block [33].
These blocks were designed to prevent the gradient vanishing problem and to
increase the capacity of CNN models by balancing the depth and width of the
deep architecture. Some skip connections were also introduced to transmit the
features between different convolutional layers to avoid information loss and
enhance the feature representation.

The current CNN classification backbones have already had strong capacity and
feature representation ability. So in this paper, we aim to discover the
potential of the multi-scale deep features extracted from different stages,
and further improve the performance for any well-trained CNN backbones for
histopathological tissue classification.

Figure 1: Overview of the proposed Pyramidal Deep-Broad Learning (PDBL). In
PDBL, we create an image pyramid for each image in the training set to obtain
the pyramidal contextual information. For each image in the pyramid, we
extract the multi-scale deep-broad features by a Deep-Broad block (DB-block).
Finally, histopathological image categories can be inferred by the broad
learning system. For simplification, we only show the notations of the left
input image.

###  II-C Broad Learning System

With the breakthrough of the GPU architectures over the past decade,
researches have kept increasing the depth of CNN models and achieved
outstanding performance in most of the computer vision and medical imaging
tasks. However, the deeper the network is, the more computational time of the
training we spend. Chen et al. [34] proposed an opposite direction of neural
network by expanding the width instead of increasing the depth, called Broad
Learning system (BL). BL tends to breadth-wise expand the feature space by
multi-group feature mapping, and uses a shallow fully-connection layer to
calculate output, which greatly saves the computational resources comparing
with deep learning (DL).

In the past years, a series of BL approaches [35, 36, 37, 38] have been
proposed. The motivation of these approaches is to provide an architecture to
breadth-wisely combine multiple groups of features by solving the following
optimization problems:

$W_{opt}=\mathop{\arg\min}_{W_{init}}\parallel
AW_{init}-Y\parallel_{2}^{2}+\gamma\parallel W_{init}\parallel_{2}^{2}$ (2)

where $Y$, $A$ represent target matrix(vector) and combined feature matrix
that concatenates with all groups of feature nodes and enhance nodes.
$W_{init}$ and $W_{opt}$ are pre-update and post-update weights of output
layer, which can be updated rapidly by pseudo-inverse method:

$W_{opt}=A^{+}Y=\lim_{\lambda\rightarrow 0}(AA^{T}+\lambda E)^{-1}A^{T}Y$ (3)

where $E$ and $\lambda$ represent an identity matrix and a constant parameter.
When $\lambda=0$, the updating method is ridge regression which requires $A$
to be a non-singular matrix.

In short, DL has stronger semantic feature extraction ability, while BL is
faster and more lightweight. So in this paper, we want to gather the strengths
of both DL and BL by using DL to extract multi-scale semantic features and
using BL for inference. Such deep-broad design is effective and will not
introduce extra training burden. Considering histopathological tissue
classification, we designed a pyramidal structure for the deep-broad learning
to consider the pyramidal contextual information of the histopathological
images.

TABLE I: Symbol annotations $\mathbf{Symbol}$ | $\mathbf{Meaning}$
---|---
$I$ | Input image (WSI patch)
$X$ | Sub-image of the image pyramid
$\mathcal{X}$ | Image pyramid
$f$ | Intermediate deep features
e | Channel-wise feature vector of $f$
z | Multi-scale deep-broad feature vector of $X$
b | Pyramidal deep-broad feature vector of $\mathcal{X}$
B | Feature matrix of the complete training data
C | Covariance matrix
U | Dimension reduction matrix
A | Feature matrix after PCA
$W_{\text{PBDL}}$ | Weights of PDBL
$n$ | Number of the training samples
$c$ | Number of the categories
$q$ | Dimension of z
$p$ | Dimension of b
$d$ | Target dimension in dimensionality reduction

##  III Pyramidal Deep-Broad Learning

Deep learning (DL) has the powerful feature extraction ability while Broad
Learning (BL) is good at combining multiple groups of features for fast
inference. Theoretically, associating DL with BL can effectively improve the
performance of existing CNN-based models. In this section, we proposed a novel
Pyramidal Deep-Broad Learning (PDBL) with a Deep-Broad block (DB-block) for
histopathological tissue classification. Fig. 1 demonstrates the overview of
the proposed PDBL. We first construct an image pyramid for the input image.
And then we extract the multi-scale deep-broad features by DB-block, shown in
Fig. 2. Finally, a broad learning system was introduced for the inference. The
annotations in this article are defined in Table I.

###  III-A Image Pyramid Construction

Typically, pathologists observe pathological sections under different
magnifications. Inspired by this, we want the CNN-based models to be able to
consider the multi-scale contextual information too. Given an input WSI patch
$I$, we construct an image pyramid by downscaling the input image for $s$
times, defined as follows:

$X_{i}=\xi(I,S_{i}),i=1,...,s\ $ (4)

where $\xi$ denotes the scale transformation of bilinear interpolation with
the scaling factor $S_{i}$.

Now we have an image pyramid $\mathcal{X}$ with $s$ sub-images, including the
input image $X_{1}=I$.

$\mathcal{X}=\\{X_{1},X_{2},...,X_{s}\\}$ (5)

And then, each sub-image is passed into our proposed DB-block for feature
extraction.

Figure 2: Illustration of DB-block. Given any CNN backbone, intermediate
feature maps of each stage are compressed into channel-wise features e by the
adaptive global average pooling. Then we can obtain the multi-scale deep-broad
features z by concatenating e.

###  III-B Multi-Scale Deep-Broad Feature Extraction via Deep-Broad Block

We proposed a Deep-Broad block (DB-block) to extract features for the image
pyramid $\mathcal{X}$, as shown in Fig. 2. The DB-block broadens the deep
features from each stage of the CNN backbone and forms the multi-scale deep-
broad features.

For the sub-image $X_{i}$, we first extract its intermediate deep features $f$
from the last layer of each stage of the CNN backbones. The intermediate deep
features $f_{k}$ at the $k$-th stage are squeezed into a channel-wise feature
vector by adaptive global average pooling in Eq. 6.

$\textbf{e}_{k}=\frac{1}{H_{k}\times
W_{k}}\sum_{i=1}^{H_{k}}\sum_{j=1}^{W_{k}}f_{k}(i,j)$ (6)

where $H_{k}$ and $W_{k}$ represent the height and width of the intermediate
feature maps at the $k$-th stage.

The multi-scale deep-broad features $\textbf{z}_{i}$ of the sub-image $X_{i}$
can be obtained by concatenating all the channel-wise feature vector e, as
follows:

$\textbf{z}_{i}=\textbf{e}_{1}\oplus\textbf{e}_{2}\oplus,...,\oplus\textbf{e}_{h}$
(7)

where $\textbf{z}_{i}\in\mathbb{R}^{1\times q}$, and $q$ represents the
dimensions of this feature group, which is the summation of channel numbers of
all the stages. $h$ denotes the number of stages of the CNN backbone. $\oplus$
represents the concatenation operation.

Now for the image pyramid $\mathcal{X}=\\{X_{1},X_{2},...,X_{s}\\}$, we can
obtain pyramidal deep-broad feature vector b with pyramidal contextual
information by concatenating all the multi-scale deep-broad features z of the
sub-images as follows:

$\textbf{b}=\textbf{z}_{1}\oplus\textbf{z}_{2}\oplus,...,\oplus\textbf{z}_{s}$
(8)

where $\textbf{b}\in\mathbb{R}^{1\times p}$ is the pyramidal deep-broad
feature vector of the image pyramid $\mathcal{X}$. And $p$ is the dimension of
b where $p=s\times q$.

In DB-block, we extract the deep learning features by the baseline model pre-
trained with only train the CNN backbone once using the original training set.
And we extract features of all the sub-images from the image pyramid by this
CNN backbone.

###  III-C Broad Learning Inference

With the pyramidal deep-broad feature vector b of the image pyramid
$\mathcal{X}$, we apply a broad learning system for inference. Let us denote
the complete training samples as
$\mathbb{X}=\\{{\mathcal{X}_{i}|i=1,2,...,n}\\}$, where $n$ is the number of
training samples. We can obtain a set of feature vectors
$\\{\textbf{b}_{i}|i=1,2,...,n\\}$. Then a broad feature matrix B is
constructed by combining all the feature vectors $\textbf{b}_{i}$ in Eq. 9.

$\textbf{B}=\delta([\textbf{b}_{1}^{\mathsf{T}},\textbf{b}_{2}^{\mathsf{T}},...,\textbf{b}_{n}^{\mathsf{T}}]^{\mathsf{T}}),\textbf{B}\in\mathbb{R}^{n\times
p}$ (9)

where $\delta$ denotes the matrix normalization transformation.

In order to reduce the feature dimension and redundancy, principal component
analysis (PCA) is employed for B. We first calculate the covariance matrix C
by Eq. 10:

$\textbf{C}=\frac{1}{n}\textbf{B}^{\mathsf{T}}\textbf{B},\textbf{C}\in\mathbb{R}^{p\times
p}.$ (10)

Furthermore, dimension reduction matrix is obtained by SVD algorithm:

$[\textbf{U},\Sigma,\textbf{V}^{\mathsf{T}}]=SVD(\textbf{C}).$ (11)

where $\textbf{U}\in\mathbb{R}^{p\times d}$, and $d$ represents the target
dimension in dimensionality reduction.

According to U, the dimension of feature matrix B can be reduced to the matrix
A by Eq. 12:

$\textbf{A}=\textbf{B}\times\textbf{U},\textbf{A}\in\mathbb{R}^{n\times d}.$
(12)

Finally, the probabilities $Y$ of all the categories can be calculated by:

$Y=\textbf{A}W_{\text{PDBL}},W_{\text{PDBL}}\in\mathbb{R}^{d\times c}$ (13)

where $W_{\text{PDBL}}$ is the weights of PDBL, which can be calculated by the
pseudo-inverse algorithm. $c$ denotes the number of categories.

In the training phase, we use $\textbf{A}_{train}$ to represent the feature
matrix for the complete training set after PCA. According to the ground truth
labels $Y_{train}$, the weights $W_{\text{PDBL}}$ of PDBL can be calculated by
the pseudo-inverse algorithm as follows:

$W_{\text{PDBL}}=\textbf{A}_{train}^{+}Y_{train}$ (14)

where $\textbf{A}_{train}^{+}$ can be calculated by:

$\textbf{A}_{train}^{+}=\mathbf{lim}_{\lambda\rightarrow 0}(\lambda
E+\textbf{A}_{train}\textbf{A}_{train}^{T})^{-1}\textbf{A}_{train}^{T}$ (15)

where $\lambda$ and $E$ represent a constant and unit matrix, respectively.
Since pseudo-inverse algorithm only updates the weights for once, it greatly
saves the computational resources.

In the testing phase, we can obtain the feature matrix $\textbf{A}_{test}$ by
the same steps of the training phase. Then we can infer the probabilities
$Y_{test}$ of tissue categories by:

$Y_{test}=\textbf{A}_{test}W_{\text{PDBL}}$ (16)

The final classification results is the tissue categories with the largest
probabilistic value.

Since the proposed PDBL is a plug-and-play module. It can be applied to any
CNN backbone and further improve the classification performance.

###  III-D Implementation and Training Details

In our experiments, all the CNN backbones were implemented in PyTorch. The
backbones were trained on an NVIDIA RTX 2080Ti with the cross entropy loss and
the SGD optimizer with $1e-3$ learning rate and $0.9$ momentum. The batchsize
was set to $20$. Patches were resized into $224\times 224$ and normalized in
both the training and test phase. All the backbones were fine-tuned for $50$
epochs with weights provided by PyTorch which were pretrained on ImageNet.

For our proposed PDBL, we created image pyramid under three different
resolution ($112\times 112$, $160\times 160$ and $224\times 224$). The deep-
broad features of the images of different scales were all extracted from the
network trained by the original image size ($224\times 224$). And the target
dimension $d$ of PCA was decided by the total number $n$ of the training
samples as follows:

$b=\left\\{\begin{array}[]{cc}0.9*n,&n<=2000\\\
2000,&n>2000\end{array}\right.$ (17)

where $n$ is the number of training samples.

##  IV Datasets

We evaluate our proposed PDBL in the following two datasets Kather Multiclass
Dataset [11] and LC25000 Dataset [22].

###  IV-A Kather Multiclass Dataset

Kather Multiclass Dataset is a multi-class colorectal dataset composed of H&E
stained histopathological tissue patches, which was published by J. N.
Kather222https://zenodo.org/record/1214456. Kather Multiclass Dataset is
composed of two subsets. Kather et al. [11] manually delineated tissue regions
in $86$ colorectal (CRC) tissue slides and they extracted $100,000$ H&E
histopathological tissue patches from these regions as Kather Multiclass
Internal (KMI) subset at $20\times$ magnification. They also extracted an
additional independent Kather Multiclass External (KME) subset include $7180$
H&E stained histopathological patches. Histopathological images in Kather
Multiclass Dataset are cropped to a square size of $224\times 224$ pixels from
primeval whole slide images (WSIs) at $20\times$ magnification.

As shown in Fig. 3, each histopathological image belongs to one category of
tissues and there are $9$ categories of tissues in Kather Multiclass Dataset
such as adipose (ADI), background (BACK), debris (DEB), lymphocytes (LYM),
mucus (MUC), smooth muscle (MUS), normal colon mucosa (NORM), cancer-
associated stroma (STR), and colorectal adenocarcinoma epithelium (TUM).

###  IV-B LC25000 Dataset

To advance computer-aided automated analysis of lung and colon carcinomas,
[22] released a lung and colon histopathological image dataset (LC25000
Dataset)333https://github.com/tampapath/lung_colon_image_set. In LC25000
Dataset, histopathological images are cropped to a square size of $768\times
768$ pixels from H&E stain WSIs of lung carcinoma and colon carcinoma. As
shown in Fig. 4, LC25000 Dataset has $5$ categories such as benign lung
tissues (LN), lung adenocarcinomas (LAC), lung squamous cell carcinomas
(LSCC), benign colonic tissues (CN), and colon adenocarcinomas (CAC). LC25000
Dataset is a balanced dataset that each class that has $5000$
histopathological images.

|  |  |  |
---|---|---|---|---
ADI | BACK | DEB | LYM | MUC

|  |  |
---|---|---|---
MUS | NORM | STR | TUM

Figure 3: Kather Multiclass Dataset includes Adipose (ADI), background (BACK), debris (DEB), lymphocytes (LYM), mucus (MUC), smooth muscle (MUS), normal colon mucosa (NORM), cancer-associated stroma (STR), colorectal adenocarcinoma epithelium (TUM). |  |  |  |
---|---|---|---|---
LN | LAC | LSCC | CN | CAC
Figure 4: LC25000 Dataset has 5 classes which includes benign lung tissues
(LN), lung adenocarcinomas (LAC), lung squamous cell carcinomas (LSCC), benign
colonic tissues (CN) and colon adenocarcinomas (CAC)

##  V Experimental Results

In this section, we conducted several studies to evaluate the proposed PDBL.
In the following experiments, PDBL was respectively plugged on three common
classification architectures, including EfficientNet-b0 [15], ResNet50 [13]
and a lightweight model ShuffleNetV2 [32]. In Section V-A, we evaluate the
effectiveness of PDBL with different proportions of training samples. In
Section V-B, we test the limit of the proposed PDBL by an extremely difficult
task by leveraging only $1\%$ training samples to inference the rest of them
(99%). In Section V-C, we conduct an ablation study to verify the
effectiveness and the necessity of the pyramidal design. Next, we demonstrate
the advantages of rapid domain adaptation on PDBL in Section V-D. Finally, we
also show the WSI-level semantic segmentation results by stitching the patch-
level classification results.

Accuracy and macro F1 score were used to evaluate the patch-level
classification performance of the proposed PDBL in all the experiments. Due to
the page limit, in Table II and Table III, we only demonstrate the accuracy.
F1 scores can be found in the supplementary materials.

TABLE II: Quantitative comparison (ACC) with different proportions of training samples in Kather Dataset (100% training samples: 100,000 patches) and LC25000 Dataset (100% training samples: 15,000 patches). Models | Dataset | Accuracy
---|---|---
1% | 5% | 10% | 25% | 35% | 50% | 70% | 100%
ShuffLeNetV2 | Baseline+PDBL | Kather | 0.86880 | 0.92242 | 0.94095 | 0.94666 | 0.94568 | 0.94791 | 0.95139 | 0.95097
Baseline* | 0.17187 | 0.73607 | 0.91448 | 0.94805 | 0.96058 | 0.96114 | 0.95097 | 0.95195
Baseline*+PDBL | 0.87869 | 0.93816 | 0.95696 | 0.96114 | 0.96365 | 0.96253 | 0.96198 | 0.96156
Baseline+PDBL | LC25000 | 0.93260 | 0.96730 | 0.97070 | 0.98670 | 0.98950 | 0.99300 | 0.99320 | 0.99400
Baseline* | 0.39280 | 0.62380 | 0.78980 | 0.96810 | 0.97310 | 0.98480 | 0.98890 | 0.99660
Baseline*+PDBL | 0.94140 | 0.95970 | 0.96540 | 0.98990 | 0.99140 | 0.99370 | 0.99520 | 0.99720
EfficientNetb0 | Baseline+PDBL | Kather | 0.86685 | 0.90933 | 0.92214 | 0.92841 | 0.93064 | 0.93036 | 0.93287 | 0.93398
Baseline* | 0.91267 | 0.93788 | 0.94109 | 0.94011 | 0.94415 | 0.93816 | 0.93524 | 0.94471
Baseline*+PDBL | 0.92256 | 0.94930 | 0.94847 | 0.94930 | 0.95432 | 0.95557 | 0.95641 | 0.96086
Baseline+PDBL | LC25000 | 0.94170 | 0.97290 | 0.98200 | 0.99000 | 0.99360 | 0.99450 | 0.99530 | 0.99500
Baseline* | 0.87220 | 0.96320 | 0.98100 | 0.98960 | 0.99230 | 0.99630 | 0.99780 | 0.99870
Baseline*+PDBL | 0.95300 | 0.97990 | 0.98750 | 0.99460 | 0.99610 | 0.99710 | 0.99860 | 0.99940
ResNet50 | Baseline+PDBL | Kather | 0.86797 | 0.92159 | 0.93315 | 0.93538 | 0.93788 | 0.93538 | 0.93858 | 0.93663
Baseline* | 0.93844 | 0.95682 | 0.94930 | 0.94582 | 0.95153 | 0.95850 | 0.95292 | 0.95432
Baseline*+PDBL | 0.93928 | 0.95933 | 0.95752 | 0.95543 | 0.96212 | 0.96253 | 0.96142 | 0.96421
Baseline+PDBL | LC25000 | 0.93400 | 0.96960 | 0.97750 | 0.98910 | 0.98800 | 0.99140 | 0.99280 | 0.99250
Baseline* | 0.95570 | 0.97940 | 0.98470 | 0.99380 | 0.99510 | 0.99770 | 0.99860 | 0.99890
Baseline*+PDBL | 0.95120 | 0.98080 | 0.98560 | 0.99470 | 0.99550 | 0.99810 | 0.99910 | 0.99950

###  V-A Effectiveness of PDBL with Different Proportions of Training Samples

In this experiment, we evaluate the effectiveness and efficiency of the
proposed PDBL with different proportions of the training set in both datasets.
For Kather dataset, Kather Multiclass Internal set ($100k$ patches) and Kather
Multiclass External set ($7k$ patches) are the complete training set and test
set. For LC25000 dataset, we let 60% and 40% samples as the training set and
the test set. And then, the training sets of two datasets were randomly split
into eight incremental subsets with the proportions of [1%, 5%, 10%, 25%, 35%,
50%, 70%, 100%], respectively. We conducted this experiment by comparing three
models for every CNN backbone. (1) PDBL directly plugged on the baseline
models pre-trained by ImageNet [12], denoted as Baseline+PDBL. (2) Baseline
models pre-trained by ImageNet fine-tuned for $50$ epochs without PDBL,
denoted as Baseline*. (3) Baseline models pre-trained by ImageNet fine-tuned
for $50$ epochs with PDBL, denoted as Baseline*+PDBL. The above notations are
used in all the experiments.

Quantitative results for three CNN backbones on two datasets are demonstrated
in Table II. First, let us compare the baseline models trained by different
training set proportions with and without PDBL (Baseline* and Baseline*+PDBL).
We can observe an overall improvement for nearly all the CNN backbones when
equipping with PDBL. With the increasing training samples, the improvements
became less significant for three CNN backbones. It can be observed that the
degree of the improvements actually depends on three factors, the complexity
of the models, the difficulties of the datasets, and the ratio of the training
samples. For example, for the deeper backbone with more parameters like
ResNet50, the overall improvement on LC25000 is more obvious than Kather.
Because LC25000 only has 5 classes but Kather contains 9 classes. When with
enough training samples in LC25000, say 35%, ResNet50 with PDBL gets only less
than 0.001 improvement (0.99510$\rightarrow$0.99550) but it can bring 0.018
improvement (0.97310$\rightarrow$0.99140) for a lightweight backbone
ShuffLeNetV2. It is because ResNet50 has a stronger feature representation
capability than ShuffLeNetV2, which is a trade-off between computational
resources and performance. Despite this, the proposed PDBL is able to steadily
improve the classification performance for all the CNN backbones.

Since the proposed PDBL is a model-agnostic plug-and-play module, we are also
curious about if semantic features extracted from the backbones trained with
ImageNet can be utilized for inference in histopathology images. So for each
CNN backbone, we directly applied PDBL on them and updated PDBL for once
(Baseline+PDBL). In this comparison, different proportions of training samples
were only used to update PDBL and did not involve in network training. As can
be seen in Table II, baseline models with PDBL can achieve stable
classification performances in both datasets comparing with other two
competitors (Baseline* and Baseline*+PDBL). This observation proves that the
proposed PDBL can effectively make the inference only relying on the visual
semantic features from natural images. For the easier dataset LC25000 with
only 5 classes, the performance gap between Baseline+PDBL and Baseline*+PDBL
is relatively small (less than 0.005). For the more difficult dataset Kather
with 9 classes, the performance gap between them is around 0.03 accuracy, but
the accuracy is still over 0.90. It is interesting that when with smaller
training samples (less than 10%), Baseline+PDBL occasionally outperforms
Baseline* even Baseline*+PDBL, especially for the lightweight model
ShuffLeNetV2. It is because lack of training samples may lead to overfitting
or unstable training problems. For extremely small training data, PDBL
provides an alternative solution that uses a well-trained CNN model on
ImageNet for feature extraction and to use PDBL for inference. This solution
does not need to re-train the network and can greatly save computational
resources and annotation efforts.

Based on this experiment, we further discuss the robustness of PDBL with an
extremely small training set (only 1%) in the next experiment.

TABLE III: Robustness of PDBL in an extreme experiment in Kather Dataset (training set: 1,000 patches; test set: 99,000 patches) and LC25000 Dataset (training set: 250 patches; test set: 24,750 patches). Models | Dataset | Accuracy
---|---|---
Fold1 | Fold2 | Fold3 | Fold4 | Fold5 | Mean$\pm$SD
ShuffLeNetV2 | Baseline+PDBL | Kather | 0.90529 | 0.91216 | 0.91456 | 0.90605 | 0.90610 | 0.90883$\pm$0.004
Baseline* | 0.14412 | 0.14964 | 0.15066 | 0.14828 | 0.14817 | 0.14817$\pm$0.002
Baseline*+PDBL | 0.91164 | 0.91679 | 0.91831 | 0.91517 | 0.91819 |  0.91602$\pm$0.003
Baseline+PDBL | LC25000 | 0.95733 | 0.95956 | 0.95095 | 0.95301 | 0.95051 |  0.95427$\pm$0.004
Baseline* | 0.42461 | 0.43818 | 0.54149 | 0.41709 | 0.62816 | 0.48991$\pm$0.092
Baseline*+PDBL | 0.95499 | 0.95414 | 0.95079 | 0.95123 | 0.95103 | 0.95244$\pm$0.002
EfficientNetb0 | Baseline+PDBL | Kather | 0.92166 | 0.92842 | 0.92677 | 0.92209 | 0.91803 | 0.92339$\pm$0.004
Baseline* | 0.90577 | 0.90841 | 0.91234 | 0.90911 | 0.91201 | 0.90953$\pm$0.003
Baseline*+PDBL | 0.94894 | 0.95266 | 0.95602 | 0.94843 | 0.94926 |  0.95106$\pm$0.003
Baseline+PDBL | LC25000 | 0.95564 | 0.95681 | 0.95487 | 0.95725 | 0.94978 | 0.95487$\pm$0.003
Baseline* | 0.91640 | 0.93689 | 0.93345 | 0.94101 | 0.93782 | 0.93311$\pm$0.010
Baseline*+PDBL | 0.96230 | 0.96194 | 0.95984 | 0.95960 | 0.96299 |  0.96133$\pm$0.002
ResNet50 | Baseline+PDBL | Kather | 0.92028 | 0.92207 | 0.92342 | 0.91214 | 0.91886 | 0.91935$\pm$0.004
Baseline* | 0.96063 | 0.95318 | 0.95742 | 0.95080 | 0.95435 | 0.95528$\pm$0.004
Baseline*+PDBL | 0.96078 | 0.95559 | 0.95957 | 0.95207 | 0.95637 |  0.95688$\pm$0.003
Baseline+PDBL | LC25000 | 0.95091 | 0.95147 | 0.94513 | 0.93996 | 0.94404 | 0.94630$\pm$0.005
Baseline* | 0.96465 | 0.96756 | 0.95891 | 0.96752 | 0.96311 | 0.96435$\pm$0.004
Baseline*+PDBL | 0.96537 | 0.96642 | 0.96040 | 0.96655 | 0.96570 |  0.96489$\pm$0.003

###  V-B Robustness of PDBL with an extreme proportions of training and test
sets.

In the previous experiment, we observed a surprising outstanding performance
of PDBL with a small training set. In this experiment, we conducted a cross-
validation-like experiment, but more difficult than cross-validation, to
further discuss the robustness of the proposed PDBL. We randomly sampled 1%
patches (1000 for Kather; 250 for LC25000) from the training data and used the
rest 99% patches (99000 for Kather; 24750 for LC25000) of the training data
for testing. To alleviate the sampling bias, we repeated this experiment five
times. Each repetition is regarded as a fold in this experiment. All
Baselines* were pre-trained by ImageNet and fine-tuned for 50 epochs.

According to the accuracy in Table III, we notice that the performances of the
fine-tuned lightweight models ShuffLeNetV2 and EfficientNetb0 are unstable
(Baseline*). However, even the quantitative results of ShuffLeNetV2 are nearly
random guesses, the proposed PDBL (Baseline*+PDBL) can still drastically boost
the accuracy from 0.14 to 0.91. Because CNN backbones make an inference by
only relying on the high-level semantic features from the last layer. PDBL
enriches the features from multiple stages of the CNN backbones, which is
beneficial to make an inference when the training samples are very limited.
Furthermore, results of Baseline*+PDBL are more stable with smaller standard
deviation.

In this experiment, we also directly plugged PDBL on the pre-trained models
(Baseline+PDBL). Quantitative results prove that the features extracted from
the model pre-trained by ImageNet are also stable and effective. For
ShuffLeNetV2, it is even higher than Baseline*+PDBL. The reason may be that
ShuffLeNetV2 failed to learn a stable model with such limited training
samples. Since 1% training samples of LC25000 is 250, it means that there are
only 50 patches for each class. Instead of training an unstable model for
histopathology images, PDBL coupling with a more stable model trained from
natural images may be an alternative solution for very few training samples.

###  V-C Ablation Study of Pyramidal Design

In this experiment, an ablation study was conducted to evaluate the
effectiveness of the pyramidal design of the proposed PDBL. We keep the same
experimental setting of Section V-B by using 1% samples for training and 99%
samples for testing five times. Here, we only show the mean and standard
deviation of the results. Note that +DBL and +PDBL represent our method
without and with the image pyramid, respectively.

Table IV and Table V demonstrate the accuracy and F1 score on both datasets
with only 1% training samples. We can observe a constantly stable improvement
with pyramidal design for all three backbones on both datasets. According to
the quantitative results, we found that the improvement of pyramidal design is
more significant on the lightweight model. In Kather, PDBL can introduce
around 0.09 improvement for ShuffLeNetV2, 0.03 for EfficientNetb0 and only
0.0006 for ResNet50 respectively. In fact, comparing ShuffLeNetV2 with
ResNet50, ResNet50 has stronger capacity and feature representation ability,
which are great advantages when very few training samples are involved.

In addition, the same finding is also verified in this experiment that
pyramidal design brings more improvement for more difficult datasets. For
Kather with 9 classes, pyramidal design can introduce multi-level contextual
features, which can support the feature representation for more classes.

TABLE IV: Ablation of Pyramidal Design (Accuracy) Models | Accuracy (Mean$\pm$SD)
---|---
Kather | LC25000
ShuffLeNetV2 | +DBL | 0.81173$\pm$0.013 | 0.94772$\pm$0.005
+PDBL |  0.90883$\pm$0.004 |  0.95427$\pm$0.004
EfficientNetb0 | +DBL | 0.89725$\pm$0.004 | 0.95175$\pm$0.004
+PDBL |  0.92339$\pm$0.004 |  0.95487$\pm$0.003
ResNet50 | +DBL | 0.91912$\pm$0.004 | 0.94462$\pm$0.003
+PDBL |  0.91935$\pm$0.004 |  0.94630$\pm$0.005
TABLE V: Ablation of Pyramidal Design (F1 score) Models | F1 (Mean$\pm$SD)
---|---
Kather | LC25000
ShuffLeNetV2 | +DBL | 0.81464$\pm$0.013 | 0.94793$\pm$0.005
+PDBL |  0.90947$\pm$0.004 |  0.95440$\pm$0.004
EfficientNetb0 | +DBL | 0.89798$\pm$0.004 | 0.95191$\pm$0.004
+PDBL |  0.92342$\pm$0.005 |  0.95496$\pm$0.003
ResNet50 | +DBL | 0.91927$\pm$0.004 | 0.94476$\pm$0.003
+PDBL |  0.91984$\pm$0.004 |  0.94638$\pm$0.005
TABLE VI: Domain adaptation from Kather Multiclass Dataset to Zhao et al. [39]. * means fine-tuning CNN backbones for one epoch. $\dagger$ means updating PDBL for once. | Models | ShuffLeNetV2 | EfficientNetb0 | ResNet50
---|---|---|---|---
| Acc | F1 | Acc | F1 | Acc | F1
No re-training | Baseline | 0.77931 | 0.69545 | 0.79387 | 0.74014 | 0.83666 | 0.76890
Baseline+PDBL | 0.74612 | 0.67446 | 0.79874 | 0.75906 | 0.85536 | 0.80056
1% | Baseline+PDBL$\dagger$ | 0.96336 | 0.95980 | 0.97064 | 0.96790 | 0.96848 | 0.96524
Baseline* | 0.94626 | 0.94032 | 0.93037 | 0.92082 | 0.95875 | 0.95416
Baseline*+PDBL$\dagger$ | 0.96275 | 0.95910 | 0.97122 | 0.96845 | 0.97183 | 0.96906
10% | Baseline+PDBL$\dagger$ | 0.97952 | 0.97775 | 0.98349 | 0.98203 | 0.98190 | 0.98031
Baseline* | 0.96624 | 0.96295 | 0.96451 | 0.96071 | 0.97678 | 0.97448
Baseline*+PDBL$\dagger$ | 0.97930 | 0.97767 | 0.98325 | 0.98175 | 0.98454 | 0.98454
100% | Baseline+PDBL$\dagger$ | 0.98482 | 0.98347 | 0.98587 | 0.98463 | 0.98444 | 0.98302
Baseline* | 0.98133 | 0.97952 | 0.98114 | 0.97941 | 0.98753 | 0.98641
Baseline*+PDBL$\dagger$ | 0.98590 | 0.98463 | 0.98658 | 0.98496 | 0.98987 | 0.98887

###  V-D Domain Adaptation Study

Domain adaptation is a crucial ability for a neural network model. In this
experiment, we tested if our proposed PDBL can be easily adapted from the
source domain to the target domain. The same with Kather Multiclass Dataset,
Zhao et al. [39] also released a histopathological tissue classification
dataset for colorectal cancer with night tissue types, which were collected
from four different centers, including TCGA, Kather, Guangdong Provincial
People’s Hospital and Yunnan Cancer Center. We use the data from Guangdong
Provincial People’s Hospital ($105k$ patches) as the target domain, which were
divided into a training set ($63k$ patches) and a test set ($42k$ patches) in
this experiment. Let CNN backbones with and without PDBL be trained by the
source domain (Kather), we first directly applied them to the target domain
[39]. Next, we fine-tuned the models with 100% ($63k$ patches), 10% ($6.3k$
patches) and 1% ($636$ patches) of the training set from the target domain,
respectively. (1) The models trained on the source domain Kather Dataset (100%
KMI) are denoted as Baseline and Baseline+PDBL in this study. (2) Baseline*
represents fine-tuning CNN backbones for only one epoch and +PDBL$\dagger$
means updating PDBL by the training set of target domain. 1%, 10%, and 100%
indicate the ratio of training samples we used for fine-tuning CNN backbones
and updating PDBL.

Table VI demonstrates the quantitative results of the domain adaptation study.
When we directly apply baseline models (with and without PDBL) trained by
Kather to Zhao et al. [39], the performance drastically decreases. It means
that there exists a domain shift between two datasets. Then we fine-tuned the
models by 1%, 10%, and 100% training samples respectively, with only one epoch
fine-tuning of baselines and weights updating of PDBL. Baseline*+PDBL$\dagger$
came back to relatively high performance with the complete training set
(100%). It is interesting that even we do not re-train the baseline models,
Baseline+PDBL$\dagger$ can still obtain an outstanding performance by updating
PDBL for once with only 1% training samples (636 patches). This observation
also supports our conclusion in Section V-A that using a more stable baseline
model for feature extraction and PDBL for inference is a good solution when
training samples are limited. It greatly saves computational resources and
annotation efforts.

###  V-E Timing Statistics

Table VII demonstrates the timing statistics of updating weights of PDBL
versus training CNN backbones on Kather Multiclass Dataset. The time of one
epoch training of CNN models is the average of 50 epochs training time. Since
our proposed PDBL only needs to calculate the weights once. According to the
timing statistics, the CNN backbones can get a performance boost by only
spending around one epoch training time for the calculation of PDBL.

TABLE VII: Average training time of baseline models and PDBL on Kather Multiclass Dataset (second). |  CNN (1 epoch) | PDBL (Total)
---|---|---
|  Feature extraction |  Subsequent calculation | Total
ResNet34 | 210 | 187 | 43 | 230
EfficientNetb0 | 345 | 394 | 85 | 479
ShuffLeNetV2 | 157 | 180 | 33 | 213

###  V-F Semantic Segmentation for Whole Slide Images

The intention of patch-level tissue classification is to achieve semantic
segmentation for whole slide images. So in this experiment, we show a
colorectal WSI example of semantic segmentation by our proposed method using
the WSI from the department of pathology, Guangdong Provincial People’s
Hospital. The model is Baseline*+PDBL of ResNet50 with 100% training set in
this experiment.

Given a WSI in Fig. 5 (a), we first cut it into $224\times 224$ patches under
20$\times$ magnification using sliding windows with the step size of $104$
pixels. For the overlapping region, we decide the tissue class by a voting
strategy. So the smaller the step size is, the more the semantic segmentation
precision will be, but the more inference time it will spend. Fig. 5 (b)
demonstrates the predicted semantic segmentation results. We overlaid the
result on the whole slide image for better visualization in Fig. 5 (c).

|  |
---|---|---
(a) The original colorectal WSI | (b) Semantic segmentation result | (c) Overlaid result
Figure 5: Semantic segmentation results of the colorectal WSI. The original
colorectal WSI (a) was cut into 224$\times$224 patches then combined together
and formed the semantic segmentation mask (b) and the original image and the
mask were overlaid into an overlay (c).

##  VI Conclusion

In this paper, we perform histopathological image classification in a new
perspective by reconsidering how to make use of the deep features in order to
further improve the performance of existing CNN classification backbones.
Thus, we proposed a lightweight plug-and-play module called Pyramidal Deep-
Broad Learning for any CNN backbone without re-training burden.

We have equipped this plug-and-play module on three representative CNN
backbones and achieved a steady improvement of the performance using different
proportions of training samples. Specifically, the proposed PDBL demonstrated
good feature representation capability and inference ability when very few
training samples were involved (less than 10%), especially for the lightweight
models. With PDBL, we even do not have to re-train the baseline models. Such
properties can greatly save computational time and annotation efforts. We also
look forward to applying this plug-and-play module to more excellent CNN
backbones on the other datasets from different tumors in the future.

## References

  * [1] G. Campanella, M. G. Hanna, L. Geneslaw, A. Miraflor, V. W. K. Silva, K. J. Busam, E. Brogi, V. E. Reuter, D. S. Klimstra, and T. J. Fuchs, “Clinical-grade computational pathology using weakly supervised deep learning on whole slide images,” _Nature medicine_ , vol. 25, no. 8, pp. 1301–1309, 2019.
  * [2] J. N. Kather, A. T. Pearson, N. Halama, D. Jäger, J. Krause, S. H. Loosen, A. Marx, P. Boor, F. Tacke, U. P. Neumann _et al._ , “Deep learning can predict microsatellite instability directly from histology in gastrointestinal cancer,” _Nature medicine_ , vol. 25, no. 7, pp. 1054–1056, 2019.
  * [3] N. Coudray, P. S. Ocampo, T. Sakellaropoulos, N. Narula, M. Snuderl, D. Fenyö, A. L. Moreira, N. Razavian, and A. Tsirigos, “Classification and mutation prediction from non–small cell lung cancer histopathology images using deep learning,” _Nature medicine_ , vol. 24, no. 10, pp. 1559–1567, 2018.
  * [4] M. Gehrung, M. Crispin-Ortuzar, A. G. Berman, M. O’Donovan, R. C. Fitzgerald, and F. Markowetz, “Triage-driven diagnosis of barrett’s esophagus for early detection of esophageal adenocarcinoma using deep learning,” _Nature Medicine_ , vol. 27, no. 5, pp. 833–841, 2021.
  * [5] A. Binder, M. Bockmayr, M. Hägele, S. Wienert, D. Heim, K. Hellweg, M. Ishii, A. Stenzinger, A. Hocke, C. Denkert _et al._ , “Morphological and molecular breast cancer profiling through explainable machine learning,” _Nature Machine Intelligence_ , vol. 3, no. 4, pp. 355–366, 2021.
  * [6] Y. Fu, A. W. Jung, R. V. Torne, S. Gonzalez, H. Vöhringer, A. Shmatko, L. R. Yates, M. Jimenez-Linan, L. Moore, and M. Gerstung, “Pan-cancer computational histopathology reveals mutations, tumor composition and prognosis,” _Nature Cancer_ , vol. 1, no. 8, pp. 800–810, 2020.
  * [7] M. Y. Lu, T. Y. Chen, D. F. Williamson, M. Zhao, M. Shady, J. Lipkova, and F. Mahmood, “Ai-based pathology predicts origins for cancers of unknown primary,” _Nature_ , vol. 594, no. 7861, pp. 106–110, 2021.
  * [8] M. van Rijthoven, M. Balkenhol, K. Siliņa, J. van der Laak, and F. Ciompi, “Hooknet: Multi-resolution convolutional neural networks for semantic segmentation in histopathology whole-slide images,” _Medical Image Analysis_ , vol. 68, p. 101890, 2021.
  * [9] C. L. Srinidhi, O. Ciga, and A. L. Martel, “Deep neural network models for computational histopathology: A survey,” _Medical Image Analysis_ , vol. 67, p. 101813, 2021.
  * [10] Y. Xue, J. Ye, Q. Zhou, L. R. Long, S. Antani, Z. Xue, C. Cornwell, R. Zaino, K. C. Cheng, and X. Huang, “Selective synthetic augmentation with histogan for improved histopathology image classification,” _Medical Image Analysis_ , vol. 67, p. 101816, 2021.
  * [11] J. N. Kather, J. Krisam, P. Charoentong, T. Luedde, E. Herpel, C.-A. Weis, T. Gaiser, A. Marx, N. A. Valous, D. Ferber _et al._ , “Predicting survival from colorectal cancer histology slides using deep learning: A retrospective multicenter study,” _PLoS medicine_ , vol. 16, no. 1, p. e1002730, 2019.
  * [12] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [13] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [14] X. Zhang, X. Zhou, M. Lin, and J. Sun, “Shufflenet: An extremely efficient convolutional neural network for mobile devices,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 6848–6856.
  * [15] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 6105–6114.
  * [16] J. Zhang, Y. Xie, Q. Wu, and Y. Xia, “Medical image classification using synergic deep learning,” _Medical image analysis_ , vol. 54, pp. 10–19, 2019\.
  * [17] L. Cai, J. Gao, and D. Zhao, “A review of the application of deep learning in medical image classification and segmentation,” _Annals of translational medicine_ , vol. 8, no. 11, 2020.
  * [18] B. Gecer, S. Aksoy, E. Mercan, L. G. Shapiro, D. L. Weaver, and J. G. Elmore, “Detection and classification of cancer in whole slide breast histopathology images using deep convolutional networks,” _Pattern recognition_ , vol. 84, pp. 345–356, 2018.
  * [19] S. Wang, Y. Zhu, L. Yu, H. Chen, H. Lin, X. Wan, X. Fan, and P.-A. Heng, “Rmdl: Recalibrated multi-instance deep learning for whole slide gastric image classification,” _Medical image analysis_ , vol. 58, p. 101549, 2019\.
  * [20] Z. Han, B. Wei, Y. Zheng, Y. Yin, K. Li, and S. Li, “Breast cancer multi-classification from histopathological images with structured deep learning model,” _Scientific Reports_ , vol. 7, no. 1, p. 4172, 2017.
  * [21] M.-J. Tsai and Y.-H. Tao, “Deep learning techniques for the classification of colorectal cancer tissue,” _Electronics_ , vol. 10, no. 14, p. 1662, 2021\.
  * [22] A. A. Borkowski, M. M. Bui, L. B. Thomas, C. P. Wilson, L. A. DeLand, and S. M. Mastorides, “Lung and colon cancer histopathological image dataset (lc25000),” _arXiv preprint arXiv:1912.12142_ , 2019.
  * [23] T. C. Hollon, B. Pandian, A. R. Adapa, E. Urias, A. V. Save, S. S. S. Khalsa, D. G. Eichberg, R. S. D’Amico, Z. U. Farooq, S. Lewis _et al._ , “Near real-time intraoperative brain tumor diagnosis using stimulated raman histology and deep neural networks,” _Nature medicine_ , vol. 26, no. 1, pp. 52–58, 2020.
  * [24] L. Hou, D. Samaras, T. M. Kurc, Y. Gao, J. E. Davis, and J. H. Saltz, “Patch-based convolutional neural network for whole slide tissue image classification,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2424–2433.
  * [25] Y. Xu, Z. Jia, L.-B. Wang, Y. Ai, F. Zhang, M. Lai, I. Eric, and C. Chang, “Large scale tissue histopathology image classification, segmentation, and visualization via deep convolutional activation features,” _BMC bioinformatics_ , vol. 18, no. 1, pp. 1–17, 2017.
  * [26] S. M. Ayyad, M. Shehata, A. Shalaby, M. A. El-Ghar, and A. El-Baz, “Role of ai and histopathological images in detecting prostate cancer: A survey,” _Sensors_ , vol. 21, no. 8, p. 2586, 2021.
  * [27] H. Ni, H. Liu, K. Wang, X. Wang, X. Zhou, and Y. Qian, “Wsi-net: branch-based and hierarchy-aware network for segmentation and classification of breast histopathological whole-slide images,” in _International Workshop on Machine Learning in Medical Imaging_. Springer, 2019, pp. 36–44.
  * [28] L. Rczkowski, M. Mozejko, J. Zambonelli, and E. Szczurek, “Ara: accurate, reliable and active histopathological image classification framework with bayesian deep learning,” _Scientific reports_ , vol. 9, no. 1, pp. 1–12, 2019.
  * [29] S. C. Kosaraju, J. Hao, H. M. Koh, and M. Kang, “Deep-hipo: Multi-scale receptive field deep learning for histopathological image analysis,” _Methods_ , vol. 179, pp. 3–13, 2020.
  * [30] J. Li, W. Li, A. Sisk, H. Ye, W. D. Wallace, W. Speier, and C. W. Arnold, “A multi-resolution model for histopathology image classification and localization with multiple instance learning,” _Computers in Biology and Medicine_ , vol. 131, p. 104253, 2021.
  * [31] I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio, _Deep learning_. MIT press Cambridge, 2016, vol. 1, no. 2.
  * [32] N. Ma, X. Zhang, H.-T. Zheng, and J. Sun, “Shufflenet v2: Practical guidelines for efficient cnn architecture design,” in _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 116–131.
  * [33] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. A. Alemi, “Inception-v4, inception-resnet and the impact of residual connections on learning,” in _Thirty-first AAAI conference on artificial intelligence_ , 2017.
  * [34] C. P. Chen and Z. Liu, “Broad learning system: An effective and efficient incremental learning system without the need for deep architecture,” _IEEE transactions on neural networks and learning systems_ , vol. 29, no. 1, pp. 10–24, 2017.
  * [35] C. L. P. Chen, Z. Liu, and S. Feng, “Universal approximation capability of broad learning system and its structural variations,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 30, no. 4, pp. 1191–1204, 2019.
  * [36] M. Han, S. Feng, C. P. Chen, M. Xu, and T. Qiu, “Structured manifold broad learning system: A manifold perspective for large-scale chaotic time series analysis and prediction,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 31, no. 9, pp. 1809–1821, 2018.
  * [37] M. Xu, M. Han, C. L. P. Chen, and T. Qiu, “Recurrent broad learning systems for time series prediction,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 4, pp. 1405–1417, 2020.
  * [38] S. Feng and C. P. Chen, “Fuzzy broad learning system: A novel neuro-fuzzy model for regression and classification,” _IEEE transactions on cybernetics_ , vol. 50, no. 2, pp. 414–424, 2018.
  * [39] K. Zhao, Z. Li, S. Yao, Y. Wang, X. Wu, Z. Xu, L. Wu, Y. Huang, C. Liang, and Z. Liu, “Artificial intelligence quantified tumour-stroma ratio is an independent predictor for overall survival in resectable colorectal cancer,” _EBioMedicine_ , vol. 61, p. 103054, 2020.