Active Inference-Based Optimization of Discriminative Neural Network
  Classifiers
Commonly used objective functions (losses) for a supervised optimization of
discriminative neural network classifiers were either distribution-based or
metric-based. The distribution-based losses could compromise the generalization
or cause classification biases towards the dominant classes of an imbalanced
class-sample distribution. The metric-based losses could make the network model
independent of any distribution and thus improve its generalization. However,
they could still be biased towards the dominant classes and could suffer from
discrepancies when a class was absent in both the reference (ground truth) and
the predicted labels. In this paper, we proposed a novel optimization process
which not only tackled the unbalancedness of the class-sample distribution of
the training samples but also provided a mechanism to tackle errors in the
reference labels of the training samples. This was achieved by proposing a
novel algorithm to find candidate classification labels of the training samples
from their prior probabilities and the currently estimated posteriors on the
network and a novel objective function for the optimizations. The algorithm was
the result of casting the generalized Kelly criterion for optimal betting into
a multiclass classification problem. The proposed objective function was the
expected free energy of a prospective active inference and could incorporate
the candidate labels, the original reference labels, and the priors of the
training samples while still being distribution-based. The incorporation of the
priors into the optimization not only helped to tackle errors in the reference
labels but also allowed to reduce classification biases towards the dominant
classes by focusing the attention of the neural network on important but
minority foreground classes.
Faezeh Fallah
# Active Inference-Based Optimization of Discriminative Neural Network Classifiers

###### Abstract

Commonly used objective functions (losses) for a supervised optimization of discriminative neural network classifiers were either distribution-based or metric-based. The distribution-based losses were mostly based on the cross entropy and fitted the network model to the distribution of the training samples. This could compromise the generalization (predictive performance on unseen samples) or cause classification biases towards the dominant classes of an imbalanced class-sample distribution. The metric-based losses could make the network model independent of any distribution and thus improve its generalization. However, the metrics involved in them were binary classification metrics. This implied to decompose a multiclass classification into a series of one-vs-all classifications and then form the overall loss from an average of the one-vs-all losses. This averaging could naturally lead to a bias towards the dominant classes. Moreover, the metric-based losses could suffer from discrepancies when a class was absent in both the reference (ground truth) labels and the predicted labels. To tackle these issues, recent works have used a combination of the distribution-based and metric-based losses. In this paper, we formulated the optimization of a discriminative neural network classifier within the framework of active inference and showed that the cross entropy-based losses were indeed the variational free energy of a retrospective active inference. Then, we proposed a novel optimization process which not only tackled the unbalancedness of the class-sample distribution of the training samples but also provided a mechanism to tackle errors in the reference (ground truth) labels of the training samples. This was achieved by proposing a novel algorithm to find candidate classification labels of the training samples during the network optimization and a novel objective function for the optimizations. The algorithm could find the candidate labels of the training samples from their prior probabilities and the currently estimated posteriors on the network. The proposed objective function incorporated these candidate labels along with the original reference labels and the priors of the training samples while still being distribution-based. The proposed algorithm was the result of casting the generalized Kelly criterion for optimal betting into a multiclass classification problem. To this end, we showed that the objective function of the generalized Kelly criterion was a tight upper bound of the expected complexity of the expected free energy of a prospective active inference. This in turn allowed us to derive our proposed objective function from such an expected free energy. The incorporation of the priors into the optimization not only helped to tackle errors in the reference labels but also allowed to reduce classification biases towards the dominant classes by focusing the attention of the neural network on important but minority foreground classes.
Background and Motivation

### Active Inference

Bayesian inference enabled perception, learning, and decision making in a passive or active perceptual task. This perception could be over a categorical (multinomial) distribution of independent and mutually exclusive states. This distribution assigned one probability to each state of each observation with the sum of these probabilities for each observation being one. That is, each observation could only be in one state at a time. In an active perception, an agent actively engaged with its environment to gather information, seek preferred observations, avoid unpreferred observations, and take actions which could reduce uncertainty and maximize reward. If the states, observations, and policies (actions) could be discretized, then the tasks could be formulated over categorical distributions of the states, observations, and policies. These formed a discrete state-space model in which the time could be discrete as well. An active perception ruled by the Bayesian inference was called an **active inference**. The Bayesian inference inferred joint/posterior distribution of a generative/discriminative model by using the Bayes' theorem. For the classification/segmentation tasks addressed in this dissertation, a discriminative model was sufficient. Thus, we restricted the use of the active inference to a discriminative model and only involved the posteriors in our formulations [Smith 2022].

According to the Bayes' theorem, for each observation (o), state (s), and policy (\(\pi\)), the posterior \(p(s|o,\pi)\) could be deduced from the likelihood \(p(o|s,\pi)\) as

\[p(s|o,\pi)=\frac{p(o|s,\pi)\cdot p(s|\pi)}{p(o|\pi)} \tag{1}\]

with \(p(o|\pi)=\sum_{s|\pi}p(o|s,\pi)\) being the model evidence or the marginal likelihood. This way, the Bayesian inference enabled perception, learning, and decision making by model inversion, i.e. deduction of the posterior \(p(s|o,\pi)\) from the likelihood \(p(o|s,\pi)\). This resulted in a maximum a posteriori estimation. In a simpler approach, a maximum likelihood estimation might be followed. However, the maximum likelihood estimation was prone to overfitting because the likelihoods only encoded the aleatoric uncertainty of the model caused by noise (disturbances) in its process. The epistemic (cognitive) uncertainty of the model was reflected by the states' priors \(\left\{p(s|\pi)\right\}_{s}\) and the model evidence \(p(o|\pi)\) included in the posteriors. The computation of the model evidence implied to sum the likelihoods of every observation over all possible states. For most of the categorical distributions this computation was intractable. Also, by increasing the number of the states the number of the summation terms increased exponentially. For continuous distributions this summation mostly turned into a nonconvex integration of no closed-form (analytical) solution. To enable a computationally tractable active inference, the Bayes' theorem got approximated by minimizing

* variational free energy (VFE)1 for perception and learning Footnote 1: The term **free energy** stemmed from connections between the Bayesian inference and the Bayesian mechanics ruling free energy in particular (quantum) physics elaborated by neuroscientists [Friston 2019].
* expected free energy (EFE) for optimal decision making, planning, and action selection.

Each of the aforementioned objective functions depended on the policies (actions). Accordingly, the minimization of each of them provided an estimate of the posteriors conditioned on the policies. However, the VFE resulted from a course of policies based on the observations in the past and present but the EFE resulted from a course of policies based on the observations in the future. Thus, the VFE and the EFE respectively enabled retrospective and prospective policy evaluations. This difference mattered in the cases where optimal policies for the past or present were not the optimal policies for the future or vice versa. To derive the aforementioned objectives, negative logarithm of both sides of the Bayes' formula was taken and \(-\ln\bigl{(}p(o|\pi)\bigr{)}\) was introduced to be the self-information or surprisal2 of the model evidence \(p(o|\pi)\). Then, the VFE got defined to be the upper bound of this quantity. This way, by minimizing the VFE,
the surprisal or deviation between observations and predictions of the model got minimized or the amount of evidence an observation could provide for the model got maximized, i.e. the model evidence got maximized.

As detailed in [Smith 2022], the objective function of the VFE was given by

\[\begin{split}\mathcal{L}_{\mathrm{VFE}}&=\mathrm{KL} \Big{[}p(s|\pi)||q(s|\pi)\Big{]}-\mathrm{E}_{p(s|\pi)}\Big{[}\mathrm{ln} \big{(}q(o|s)\big{)}\Big{]}\\ &=\underbrace{\mathrm{E}_{p(s|\pi)}\Big{[}\mathrm{ln}\big{(}p(s| \pi)\big{)}-\mathrm{ln}\big{(}q(s|\pi)\big{)}\Big{]}}_{\mathrm{complexity}}- \underbrace{\mathrm{E}_{p(s|\pi)}\Big{[}\mathrm{ln}\big{(}q(o|s)\big{)}\Big{]} }_{\mathrm{accuracy}}\\ &=\sum_{s|\pi}p(s|\pi)\cdot\mathrm{ln}\big{(}p(s|\pi)\big{)}- \sum_{s|\pi}p(s|\pi)\cdot\mathrm{ln}\big{(}q(s|\pi)\big{)}-\sum_{s|\pi}p(s|\pi )\cdot\mathrm{ln}\big{(}q(o|s)\big{)}\\ &=\underbrace{\sum_{s|\pi}p(s|\pi)\cdot\mathrm{ln}\big{(}p(s|\pi )\big{)}}_{-\mathrm{entropy}}+\underbrace{\sum_{s|\pi}-p(s|\pi)\cdot\mathrm{ln} \big{(}q(o|\pi)\big{)}}_{\mathrm{cross\ entropy}}\end{split} \tag{2}\]

with \(q(\cdot)\) being the distribution approximating the true distribution \(p(\cdot)\), \(\mathrm{KL}[p(\cdot)||q(\cdot)]\) being the Kullback-Leibler (KL) divergence (dissimilarity) between \(p(\cdot)\) and \(q(\cdot)\), and \(\mathrm{E}_{p(s|\pi)}[\cdot]\) being the expectation with respect to \(p(s|\pi)\). The KL divergence was derived from the Akaike information criterion (AIC) measuring the goodness of a model in terms of its underfitting (estimation bias on seen samples) and overfitting (predictive variance on unseen samples). The AIC measured the amount of information loss (relative entropy) resulted from representing a model with another model. Here, the cross entropy was not a distance metric because the cross entropy of two identical distributions equaled their entropy. However, after subtracting the entropy from the cross entropy, the KL divergence become a distance metric. That is, the KL divergence of two identical distributions was zero [Kullback 1951, McMillan 1956]. This way, the minimization of \(\mathcal{L}_{\mathrm{VFE}}\) amounted to finding the distribution \(q(\cdot)\) which best fitted \(p(\cdot)\). The best fit was the minimizer of the **complexity** (overfitting) and the maximizer of the **accuracy**. The minimization of \(\mathcal{L}_{\mathrm{VFE}}\) was independent of \(p(s|\pi)\). Thus by adding the **entropy** term to \(\mathcal{L}_{\mathrm{VFE}}\), an objective function called the cross entropy loss was obtained as

\[\mathcal{L}_{\mathrm{CE}}=-\sum_{s|\pi}p(s|\pi)\cdot\mathrm{ln}\big{(}q(o|\pi )\big{)}. \tag{3}\]

If \(q(\cdot)\) was Gaussian, then the cross entropy loss become a sum of squared errors.

The minimization of the EFE selected optimal policies (actions) by solving the explore-exploit dilemma [Friston 2019]. That is, when information about the states were not enough, it emphasized on exploration (maximization of information gain or minimization of uncertainty). When the information was enough, it emphasized on exploitation (maximization of reward or minimization of expected complexity). The choice of the exploratory or the exploitative optimization depended on the current uncertainty and the future (expected) reward. This way, the minimization of the EFE sought the policies which could lead to future observations optimizing the trade-off between the maximization of the information gain and the maximization of the reward. These self-evidencing observations were called to be **preferred**. The incidence probability of a preferred observation \(o\) was denoted by \(p(o)\). As detailed in [Smith 2022], the objective function of the EFE was given by

\[\mathcal{L}_{\mathrm{EFE}} =\mathrm{KL}\Big{[}p(o)||q(o|\pi)\Big{]}+\mathrm{E}_{p(s|\pi)} \Big{[}\mathrm{H}\big{[}q(o|\pi)\big{]}\Big{]} \tag{4}\] \[=\underbrace{\mathrm{E}_{p(o)}\Big{[}\mathrm{ln}\big{(}p(o))- \mathrm{ln}\big{(}q(o|\pi)\big{)}\Big{]}}_{\mathrm{expected\ complexity}}+ \underbrace{\mathrm{E}_{p(s|\pi)}\Big{[}\mathrm{H}\big{[}q(o|\pi)\big{]}\Big{]} }_{\mathrm{uncertainty}}\] \[=\underbrace{\sum_{o}p(o)\cdot\Big{[}\mathrm{ln}\big{(}p(o) \big{)}-\mathrm{ln}\big{(}q(o|\pi)\big{)}\Big{]}}_{\mathrm{expected\ complexity}}+ \underbrace{\sum_{s|\pi}-p(s|\pi)\cdot\sum_{o|\pi}q(o|\pi)\cdot\mathrm{ln}\big{(} q(o|\pi)\big{)}}_{\mathrm{uncertainty}}\]
with \(\mathrm{H}\big{[}q(o|\pi)\big{]}=-\sum_{o|\pi}q(o|\pi)\cdot\ln\big{(}q(o|\pi)\big{)}\) being the entropy of \(q(o|\pi)\). This way, active inference provided a unified mathematical framework to model interdependent aspects of perception, learning, and decision making. This framework could build highly flexible and generalizable generative models which could explain neuro-cognitive behavioral processes as well as partially observable Markov decision processes [17, 20].

### Optimization of Discriminative Neural Network Classifiers

A neural network was composed of several perceptrons (nodes) in multiple layers. The layers included an input layer, some hidden layers, and an output layer. A perceptron contained a nonlinear function called an activation and was connected to other perceptrons in neighboring layers via some weights and a bias. These weights, biases, and the nonlinear activations formed **main parameters** of the neural network. Besides, the neural network had some **hyperparameters** defining its architecture and its optimization process. Neural networks have demonstrated promising results in a wide range of applications. This was due to the **universal approximation theorem** stating that a feed-forward network with a hidden layer containing a finite number of neurons (perceptrons) could approximate any continuous function on a compact subset of \(\mathbb{R}^{d}\) if and only if the used activations (perceptrons' nonlinearities) were nonpolynomial. The number of the parameters of such an approximating model defined its capacity to represent and to predict patterns. For a fully connected neural network, this number was \(\mathcal{O}(n_{\mathrm{layer}}\cdot n_{\mathrm{width}}^{2})\) where \(n_{\mathrm{layer}}\) was the number of layers (depth of the network) and \(n_{\mathrm{width}}^{2}\) was the number of perceptrons per layer (width of the network). Thus, an increase in the width increased the number of the parameters faster than an increase in the number of layers. An increase in the number of parameters increased the chance of overfitting. Moreover, a wide shallow network could fit to the patterns in the seen (training) samples but could not predict the patterns in unseen (validation or test) samples. To enhance the generalization (predictive performance on unseen samples), the neural network should contain more layers (become deeper) [1, 13, 14].

In a fully connected neural network, every perceptron was connected to all the perceptrons in its neighboring layers. This network lacked the capability of capturing regional (intra-layer) neighborhood patterns and thus needed handcrafted features to accomplish its task. To have an end-to-end neural network, directly applicable to the input samples without any preprocessing or explicit feature extraction, the features should be extracted by the network itself. This implied to capture regional (intra-layer) neighborhood patterns through limited receptive fields. The receptive field of a perceptron defined the size and the shape of the region at the input of the network affecting the output of the perceptron. The receptive field was determined by the kernel and the dept of the perceptron in the neural network. The deeper the perceptron in the network was the larger its receptive field become.

The application of a perceptron's kernel to its inputs returned a number of feature maps. By increasing the receptive field of the perceptron, the number and the abstraction level of its feature maps got increased but the size of each map got decreased. Accordingly, by using different kernels and locating the perceptrons at different depths of the network, features of different resolutions and abstraction levels could be obtained. Besides capturing subtle features and patterns, a kernel-based network enabled **weight sharing** by applying the same kernel coefficients to various regions in space. This resulted in a significantly lower number of parameters than a fully connected network and thus reduced the chance of overfitting and improved the generalization (predictive performance on unseen samples). In addition, it reduced the number of samples needed to train (optimize) the network. An easy-to-implement kernel for estimating a categorical distribution in a classification problem or a continuous distribution in a regression task was **convolutional1**. This type of kernel formed a convolutional neural network (CNN) which could be end-to-end and deep as well.

Footnote 1: In practice, many machine learning libraries avoided the _sign flip_ action involved in the convolution and thus simply implemented a cross correlation between the inputs and the kernels of each layer.

As shown in Figure 1, a neural network could be **plain** or **Bayesian**. In the plain network, each parameter, i.e. each weight, bias, or activation, had a single value. In the Bayesian network, each parameter had a vector of values representing its distribution and uncertainty. The Bayesian network was formed from an ensemble of plain networks. That is, multiple plain networks got built and then the Bayesian network's parameters got derived from a weighted
average of the plain networks' parameters with the weight of each network being the posteriors estimated by it for the training samples. Accordingly, whatever derived or concluded for the plain networks could be extended to the Bayesian networks. In the following, we simply referred to the plain neural network as the neural network. Such a network demanded an objective function and a process to optimize its parameters as well as a regularization to mitigate overfitting. A commonly used objective function for such a network was the cross entropy loss introduced in (10). The commonly used optimization processes were based on the gradient (first derivative) descent of the objective function [Kingma 2015]. The regularization was mostly done by penalizing large perceptrons' weights or dropping perceptrons of low confident weights in a method called Dropout [Gal 2015, Jospin 2022].

The gradient descent optimization relied on the fact that the opposite direction of the gradient (first derivative) of the scalar field of the objective function pointed to the minimum of the function. Accordingly, in each iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\) of this optimization, a movement in the direction of the negative gradient of the objective function at the current point updated the network's parameters. This optimization had a linear complexity with regard to the number of network's parameters. The gradient at each iteration was the average gradient of the training samples passed through the network's layers. The samples could be passed one-by-one or all at once. The former led to a stochastic and the latter led to a batch-based optimization. A complete pass through all the training samples was called an _epoch_[Dean 2012, Ruder 2016, Goodfellow 2016].

The averaging of the gradients of the batch's samples resulted in a smooth variation of the cost versus the iterations. In addition, the batch-based optimization allowed to apply vectorized and parallelized operations. However, it was restricted to convex or relatively smooth error manifolds and could only find local minima. Moreover, feeding a large batch of samples become memory intensive. The stochastic gradient descent optimization updated the network's parameters by passing one sample through the network in each iteration. This could avoid memory issues, could address nonconvex optimizations, and could even find global minima. However, due to a more frequent update of the network's parameters it resulted in fluctuating cost versus the iterations. Depending on the samples' gradients the fluctuations might never reach a minimum but rather dance around it. Moreover, the stochastic optimization could not benefit from the vectorized or the parallelized operations.

An intermediate between the stochastic and the batch-based optimization was a mini-batch-based optimization. In this approach, the training samples got divided into \(n_{\mathrm{batch}}\) disjoint batches, i.e. \(\mathbb{T}_{\mathrm{train}}=\cup_{b=1}^{n_{\mathrm{batch}}}\mathbb{T}_{b}\). Then, in each iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\), the samples of one batch got passed through the network and the average gradient of these samples updated the network's parameters. The size or the number of the batches was a hyperparameter. This way, by adapting the size or the number of the batches, the mini-batch-based optimization could utilize the vectorized and the parallelizable operations to speed up its computations while fitting the fluctuations of the cost versus the iterations to the nonconvexity of the addressed problem. Accordingly, if \(n_{\mathrm{epoch}}\) was the number of epochs, then the network was optimized by \(n_{\mathrm{it}}=(|\mathbb{T}_{\mathrm{train}}|/|\mathbb{T}_{b}|)\times n_{ \mathrm{epoch}}\) iterations. In each epoch, the batches and the samples of each batch got randomly shuffled to avoid overfitting to some of the samples.

Figure 1: A neural network with plain (single-valued) weights and biases (a), plain activations (b), Bayesian (distributed) weights and biases (c), and Bayesian activations (d).
With \(\alpha_{\mathrm{lr}}\in(0,1)\) being the learning rate (step size), \(\mathbf{\eta}^{(i)}\) being the vector of the main parameters of the neural network in the iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\), and \(\nabla_{\mathbf{\eta}^{(i)}}(\mathcal{L})\) being the gradient of a generic objective function \(\mathcal{L}\) with regard to these parameters, we had

\[\mathbf{\eta}^{(i)}=\mathbf{\eta}^{(i-1)}-\alpha_{\mathrm{lr}}\cdot\mathbf{\delta}^{(i)}. \tag{5}\]

In the gradient descent optimization, \(\mathbf{\delta}^{(i)}=\nabla_{\mathbf{\eta}^{(i-1)}}(\mathcal{L})\). This resulted in a slow convergence and sensitivity to abrupt variations of the gradient due to noise and perturbations. To speed up the convergence, to propel out of local minima, and to smooth out the gradient variations, in the method of _momentum_, \(\mathbf{\delta}^{(i)}\) got defined to be an exponentially weighted moving average (first moment) of the current and past gradients. The averaging weight was a decay rate called first moment rate \(\beta_{\mathrm{fm}}\in[0,1)\). It emphasized the importance of recent gradients to the older ones. For \(\beta_{\mathrm{fm}}=0\), the momentum boiled down to the gradient descent. For \(\beta_{\mathrm{fm}}=1\) and \(\alpha_{\mathrm{lr}}\approx 0\) it resulted in endless fluctuations of the cost versus the iterations like the movements of a ball in a frictionless bowl. Two major bottlenecks of the gradient descent and the momentum were the possibility of being trapped into saddle points (i.e. points of zero gradients in all directions) and a slow update in the directions of sparse features of weak gradients. To tackle these, the adaptive gradient algorithm (AdaGrad) defined \(\mathbf{\delta}^{(i)}\) to be the instant (current) gradient divided (normalized) by the square root of the sum of the squared gradients. This scaling allowed to avoid saddle points and adapted the gradient and thus the optimization rate in each direction to its history of updates. That is, the more a feature (direction) was updated in the past the less it would be updated in the future.

Despite of these improves, the AdaGrad was slow since the sum of the squared gradients only grew but never shrank. This growth also resulted in a rapid decay of \(\mathbf{\delta}^{(i)}\) and thus a poor performance in dealing with nonconvex objective functions and dense features (directions of strong gradients). The root mean square propagation (RMSprop) fixed these issues by replacing the sum of the squared gradients with an exponentially weighted moving average of the squared gradients. This was called second moment of the gradient. The averaging weight was a decay rate called the second moment rate \(\beta_{\mathrm{sm}}\in[0,1)\). It emphasized the importance of recent gradients to the older ones. Moreover, in the formation of \(\mathbf{\delta}^{(i)}\), the division (normalization) of the instant gradient by the second moment balanced the step size. More specifically, it decreased the step size for large gradients to prevent their explosion and increased the step size for small gradients to prevent their vanishing. The exploding and the vanishing gradients were common issues of deep neural networks.

The adaptive moment estimation (Adam) combined the momentum (first moment) with the RMSprop (second moment) to take advantages of both. This was done by defining the \(\mathbf{\delta}^{(i)}\) to be the first moment divided (normalized) by the second moment. This way, the Adam got the convergence speed from the momentum and the ability to adapt the gradients in different directions from the RMSprop (Kingma, 2015). More specifically,

\[\mathbf{\delta}^{(i)}=\hat{\mathbf{m}}^{(i)}\odot\left(\sqrt[6]{\hat{\bm {v}}^{(i)}}\oplus 10^{-8}\right) \mathbf{g}^{(i)}=\nabla_{\mathbf{\eta}^{(i-1)}}(\mathcal{L}) \tag{6}\] \[\text{biased first moment:} \mathbf{m}^{(i)}=\beta_{\mathrm{fm}}\odot\mathbf{m}^{(i-1)}\oplus(1-\beta _{\mathrm{fm}})\odot\mathbf{g}^{(i)}\] bias-corrected first moment: \[\hat{\mathbf{m}}^{(i)}=\mathbf{m}^{(i)}\odot(1-\beta_{\mathrm{fm}}^{i})\] biased second moment: \[\mathbf{v}^{(i)}=\beta_{\mathrm{sm}}\odot\mathbf{v}^{(i-1)}\oplus(1-\beta _{\mathrm{sm}})\odot\mathbf{g}^{(i)}\odot\mathbf{g}^{(i)}\] bias-corrected second moment: \[\hat{\mathbf{v}}^{(i)}=\mathbf{v}^{(i)}\odot(1-\beta_{\mathrm{sm}}^{i}).\]

All the aforementioned techniques relied on the gradient (first derivative) of the scalar field of the objective function of the neural network. The second derivative of this scalar field was represented by a Hessian matrix. Commonly used optimization techniques based on the Hessian matrix were the Newton and the quasi-Newton method, the conjugate gradient method, and the Levenberg-Marquardt algorithm (Dean, 2012; Ruder, 2016). A common way to optimize a network's parameters by any one of the derivative-based techniques was a backpropagation. This method demanded the objective function to be expressed in terms of the network's outputs (goodness of the model) and to be differentiable with respect to the outputs of every layer. In case of using the gradient of the objective function with respect to the
network's parameters, this gradient got expressed as a product of the layerwise errors. Then, the backpropagation took the following steps:

* initialized the network's parameters with random numbers.
* passed a batch through all the layers and computed the outputs of every layer.
* computed the error at the last layer by comparing the predictions with the references.
* propagated the error from the last layer to the first layer to find the error of each layer.
* expressed the gradient of the objective function as a product of the layerwise errors.
* updated the network's parameters according to (5).

### Commonly Used Objective Functions

For a probabilistic estimate, the outputs of the neural network got converted to probabilities (posteriors) by using a softmax (normalized exponential) function. This function converted a vector to another vector whose elements summed up to one and each element of the output had a monotonic relationship with an element of the input. In our case, the input vector was the network's outputs for each sample and had a length of \(n_{\mathrm{clas}}=\left|\mathbb{L}\right|\). This way, the output of the softmax function could be interpreted as a categorical probability distribution of a multinomial classification over \(n_{\mathrm{clas}}\) mutually exclusive classes. That is, every sample could only have one reference classification label. A special case of the softmax function was the sigmoid function. This function assumed that the classes were independent but not mutually exclusive. Thus, every sample could have multiple reference labels. The sigmoid function cast a multinomial classification into a series of binary (one-vs-all) classifications. Accordingly, its outputs did not necessarily sum up to one. For a sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\), the network's outputs at the \(i^{\mathrm{th}}\) iteration of the optimization formed a vector \(\mathbf{z}_{b,j}^{(i)}=\left[z_{b,j,c}^{(i)}\right]_{c\in\mathbb{L}}\). Then, the posteriors \(\hat{\mathbf{p}}_{b,j}^{(i)}=\left[\hat{p}_{b,j,c}^{(i)}\right]_{c\in\mathbb{L}}\) produced by applying the softmax function to these outputs were

\[\hat{p}_{b,j,c}^{(i)}=\frac{\exp\bigl{(}z_{b,j,c}^{(i)}\bigr{)}}{\sum_{k\in \mathbb{L}}\exp\bigl{(}z_{b,j,k}^{(i)}\bigr{)}}\in(0,1)\quad\text{ with }\quad\sum_{c\in\mathbb{L}}\hat{p}_{b,j,c}^{(i)}=1. \tag{7}\]

Accordingly, if the training samples \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) were used to optimize the network's parameters in the iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\), then \(\mathbf{L}_{b}=\left[\mathbf{l}_{b,j}\right]_{j}=\left[\mathbf{l}_{b,c} \right]_{c}=\left[b_{b,j,c}\right]_{j,c}\) was the \(\left|\mathbb{T}_{b}\right|\times n_{\mathrm{clas}}\) matrix of vectorized reference labels of these samples, \(\mathbf{Z}_{b}^{(i)}=\left[\mathbf{z}_{b,j}^{(i)}\right]_{j}=\left[z_{b,j,c}^ {(i)}\right]_{j,c}\) was the \(\left|\mathbb{T}_{b}\right|\times n_{\mathrm{clas}}\) matrix of the network's outputs for these samples, and \(\hat{\mathbf{P}}_{b}^{(i)}=\left[\hat{\mathbf{p}}_{b,j}^{(i)}\right]_{j}= \left[\hat{p}_{b,j,c}^{(i)}\right]_{j,c}\) was the \(\left|\mathbb{T}_{b}\right|\times n_{\mathrm{clas}}\) matrix of their classification posteriors estimated by the network.

If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vector \(\mathbf{l}_{b,j}\) was a one-hot-encoding of its reference label \(l_{b,j}\in\mathbb{L}\) and was given by

\[\mathbf{l}_{b,j}=\left[l_{b,j,c}\right]_{c\in\mathbb{L}}\quad\text{ with }\quad l_{b,j,c}=\begin{cases}1&\text{if }c=l_{b,j}=\text{reference label of }v_{b,j}\in\mathbb{T}_{b}\\ 0&\text{otherwise}\end{cases}. \tag{8}\]

If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were not provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vector \(\mathbf{l}_{b,j}\) was

\[\mathbf{l}_{b,j}=\left[l_{b,j,c}\right]_{c\in\mathbb{L}}=\frac{1}{n_{\mathrm{ clas}}}\odot\mathbf{1}_{n_{\mathrm{clas}}=\left|\mathbb{L}\right|}. \tag{9}\]

For a discriminative neural network classifier acting on \(\left|\mathbb{L}\right|=n_{\mathrm{clas}}\) classes, a common way to evaluate the estimated posteriors against the reference labels was to use the cross entropy loss introduced in (3). In this application, the policies \(\pi\) incorporated in (3) represented the network's parameters. Each state \(s\) was a class \(c\in\mathbb{L}\) and each observation \(o\) was a sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\). Accordingly, \(p(s|\pi)=p(s)\) was the occurrence probability of a class (state) \(s\) which could be represented by the vectorized reference labels of the samples (observations). Also, \(q(o|\pi)\) was the classification
posterior estimated by the network's parameters \(\pi\) for the reference classification label of a sample (observation) \(o\). With these, the cross entropy loss of the discriminative neural network classifier become

\[\mathcal{L}_{\mathrm{CE}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{-1}{| \mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in\mathbb{L}} l_{b,j,c}\cdot\ln\bigl{(}\hat{p}_{b,j,c}^{(i)}\bigr{)}. \tag{10}\]

If the posteriors were generated by the softmax function, then this loss was called a **softmax cross entropy loss**. As detailed in (2), the cross entropy loss resulted from the minimization of the VFE through minimizing the KL divergence (dissimilarity) between the reference distribution \(p(\cdot)\) and the estimated distribution \(q(\cdot)\). In a categorical classification, the reference distribution \(p(\cdot)\) was the histogram of the class-sample distribution of the training samples. The estimated distribution \(q(\cdot)\) was a known function parametrized with the network's parameters. This way, the cross entropy loss and the objective functions of the active inference compared the distributions and thus were **distribution-based**. If the class-sample distribution of the training samples was imbalanced, then it had maxima at the dominant classes. These maxima formed minima of the cross entropy loss. Thus, any minimizer of the cross entropy loss could be trapped into those minima and could thus return classifications biased towards the dominant classes of the training samples.

To reduce the impacts of the dominant classes on the optimization of a neural network, the cross entropy loss got weighted and/or modulated. The resulting losses included

1. **weighted cross entropy loss** which weighted the contribution of each class \(c\in\mathbb{L}\) by the inverse of its frequency \(w_{b,c}\in(0,1)\) in the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) and (optionally) weighted the contribution of each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) by its distance \(d_{b,j,1}\in\mathbb{R}_{\geq 0}\) to the border of the nearest class and its distance \(d_{b,j,2}\in\mathbb{R}_{\geq 0}\) to the border of the second nearest class through the weight \(w_{b,j}\in(0,1)\)[Ronneberger, 2015, Badrinarayanan, 2016] \[\mathcal{L}_{\mathrm{WCE}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{-1 }{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in \mathbb{L}}w_{b,j,c}\cdot l_{b,j,c}\cdot\ln\bigl{(}\hat{p}_{b,j,c}^{(i)}\bigr{)}\] (11) \[w_{b,j,c}=w_{b,c}+w_{b,j}=\frac{\sum_{k\in\mathbb{L}}|\mathbb{T}_{b,k}|}{ \mathbb{T}_{b,c}|+10^{-8}}+\underbrace{w_{\mathrm{mo}}\cdot\exp\Bigl{(}- \frac{(d_{b,j,1}+d_{b,j,2})^{2}}{2\cdot\sigma_{\mathrm{mo}}^{2}}\Bigr{)}}_{w_{ b,j}\in(0,1)}\] (12) with \(w_{\mathrm{mo}}=10\), \(\sigma_{\mathrm{mo}}=5\), and \(|\mathbb{T}_{b,c}|=\mathrm{card}\Bigl{(}\{l_{b,j,c}=1\}\Bigr{)}\). The distances to the classification borders could be computed by applying morphological operators to the samples in the classification domain, e.g. the spatial domain in an image segmentation task.
2. **focal (modulated cross entropy) loss** which weighted the contribution of each class by the difficulty of classifying its samples with the difficulties being highlighted with a modulation factor \(\gamma_{\mathrm{mod}}\in\mathbb{R}_{+}\). That is, the higher the \(\gamma_{\mathrm{mod}}\in\mathbb{R}_{+}\) was, the more the easy samples got downweighted to emphasize the role of the difficult samples [Lin, 2018] \[\mathcal{L}_{\mathrm{FL}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{-1 }{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in \mathbb{L}}\bigl{(}1-\hat{p}_{b,j,c}^{(i)}\bigr{)}^{\gamma_{\mathrm{mod}}} \cdot l_{b,j,c}\cdot\ln\bigl{(}\hat{p}_{b,j,c}^{(i)}\bigr{)}.\] (13)
3. **weighted focal loss** which additionally weighted the contribution of each class \(c\in\mathbb{L}\) by the inverse of its frequency \(w_{b,c}\in(0,1)\) in the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\)[Lin, 2018] \[\mathcal{L}_{\mathrm{WFL}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{-1 }{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in \mathbb{L}}w_{b,c}\cdot\bigl{(}1-\hat{p}_{b,j,c}^{(i)}\bigr{)}^{\gamma_{\mathrm{ mod}}}\cdot l_{b,j,c}\cdot\ln\bigl{(}\hat{p}_{b,j,c}^{(i)}\bigr{)}.\] (14)

The weighted cross entropy and the weighted focal loss highlighted the role of the minority classes over the role of the majority classes by including the weight \(w_{b,c}\in(0,1)\) in their terms. This way, the more a class had training samples, the less its classification errors contributed to the overall loss. In a so-called class-balanced cross entropy loss [Cui, 2019], each weight \(w_{b,c}\in(0,1)\) got defined based on the _effective number_\(n_{b,c}\in(0,1)\) of the training samples of
the class \(c\in\mathbb{L}\) in the feature space as

\[w_{b,c}=\bigg{[}1-\frac{n_{b,c}-1}{n_{b,c}}\bigg{]}/\bigg{[}1-\Big{(}\frac{n_{b,c }-1}{n_{b,c}}\Big{)}^{|\mathbb{T}_{b,c}|}\bigg{]}. \tag{15}\]

This method assumed that each sample in the feature space covered a subspace and the overall samples' subspaces of each class formed its prototypical subspace. Then, the volume of this prototype defined the effective number of the class. However, in most of the applications, the feature space was hardly accessible. In a neural network, it was also variable across the network's layers. Moreover, the computation of the subspace coverages in the feature space was expensive and depending on the dimensionality and the geometry of the space. Accordingly, in [15], each number \(n_{b,c}\in(0,1)\) got handled as a hyperparameter.

The aforementioned weighting and modulation schemes could reduce the impacts of the dominant classes of the seen (training) samples on the network's optimization. However, they were still based on the cross entropy loss and thus fitted the network's model to the seen distribution. This could compromise the network's generalization (predictive performance on unseen samples) when the distribution of the unseen (validation or test) samples differed from the distribution of the seen (training) samples. An objective evaluation of a classifier on unseen samples could be done through several metrics. Among these metrics, the Dice coefficient (DICE) and its equivalent the Jaccard index (JI) provided perceptual clues, scale invariance, and counts of false positive and false negative mispredictions. The JI was also called the intersection over union (IoU) and the DICE was the F-\(\beta\) score with \(\beta=1\). These metrics could be computed with a low complexity. This enabled their integration into an iterative optimization of neural network classifiers in the form of **metric-based** losses. Then, the optimum network's parameters were the **maximizers** of the DICE [10] or the **minimizers** of the Jaccard distance (JD)=1\(-\)JI=1\(-\)IoU [1].

The DICE=F-1 score and the JD=1\(-\)JI=1\(-\)IoU directly compared the binary masks of the predicted and the reference labels of the training samples without considering their distribution. This made the network's model independent of any distribution and could thus tackle the differences of the seen and unseen distributions. However, the binary masks compared by these metrics got formed from discrete-valued labels. This hindered to integrate those metrics into a continuous optimizer with backpropagation. More specifically, the predicted labels were the results of applying an **arg max** operation to the classification posteriors \(\hat{\mathbf{p}}_{b,j}^{(i)}=[\hat{p}_{b,j,c}^{(i)}]_{c\in\mathbb{L}}\) estimated by the network. This operation was nonlinear, irreversible, and indifferentiable. Thus, to integrate the metrics into a continuous optimizer with backpropagation, the network's outputs \(\mathbf{z}_{b,j}^{(i)}=[z_{b,j,c}^{(i)}]_{c\in\mathbb{L}}\) should be stored in each iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\) and for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\). These storages got retrieved during the backpropagation and thus increased the memory footprint of the network and hindered to optimize a large network with a large number of samples per batch [1].

To integrate the aforementioned metrics into a continuous optimization framework, they should be replaced by their continuous relaxed (real-valued) surrogates. For the DICE, this surrogate compared the vectorized reference labels \(\mathbf{L}_{b}=\left[\mathbf{l}_{b,j}\right]_{j}=\left[\mathbf{l}_{b,c}\right] _{c}=\left[\mathbf{l}_{b,j,c}\right]_{j,c}\) against the classification posteriors \(\hat{\mathbf{P}}_{b}^{(i)}=\left[\hat{\mathbf{p}}_{b,j}^{(i)}\right]_{j}=[\hat {p}_{b,j,c}^{(i)}]_{j,c}\) estimated by the network as

\[\mathcal{L}_{\mathrm{DICE}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{ 2}{\left|\mathbb{L}\right|}\sum_{c\in\mathbb{L}}\frac{\sum_{j\in\mathbb{T}_{b} }l_{b,j,c}\cdot\hat{p}_{b,j,c}^{(i)}}{\sum_{j\in\mathbb{T}_{b}}\left[l_{b,j,c} ^{2}+\hat{p}_{b,j,c}^{(i)^{2}}\right]}. \tag{16}\]

The above DICE loss was reversible and differentiable and could thus be integrated into a gradient descent optimization with backpropagation [10]. However, its **nonconvexity** hindered its wide use in many applications. Other metrics such as the mean symmetric surface distance and the Hausdorff distance were also nonconvex besides being too complex for an iterative optimization process [1]. In addition, each discrete-valued metric was a set function mapping from a set of mispredictions to a set of real numbers. However, among them, only the set function of the JD was submodular. This allowed to find a convex closure of the JD in a polynomial time. This convex closure was a **convex continuous** relaxed (real-valued) surrogate taking nonnegative real-valued mispredictions as inputs. Another
metric of these properties was the Hamming distance. The convex closure of the JD got derived according to the smooth convex Lovasz extension of submodular set functions [11, 12]. The JD was defined as

\[\text{Jaccard distance (JD)=1-JII=}\frac{\left|\mathbb{V}_{\text{prd}}\cup\mathbb{V}_{ \text{ref}}\big{|}\big{\backslash}\big{|}\mathbb{V}_{\text{prd}}\cap\mathbb{V}_ {\text{ref}}\right|}{\left|\mathbb{V}_{\text{prd}}\cup\mathbb{V}_{\text{ref}} \right|}\ \ \ \ =\ \frac{\left|\mathbb{V}_{\text{prd}}\big{\backslash}\mathbb{V}_{\text{ ref}}\big{|}+\left|\mathbb{V}_{\text{ref}}\right.\right\rangle\mathbb{V}_{ \text{prd}}\big{|}}{\left|\mathbb{V}_{\text{prd}}\cup\mathbb{V}_{\text{ref}} \right|}\,. \tag{17}\]

Based on this definition, the set function of the JD for the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\text{train}}\) and the class \(c\in\mathbb{L}\) in the iteration \(i\in\{1,\cdots,n_{\text{it}}\}\) was

\[\text{JD}:\quad\mathbb{M}^{(i)}_{b,c}\in\{0,1\}^{|\mathbb{T}_{b}|} \longmapsto\frac{\operatorname{nnz}\!\left(\mathbb{M}^{(i)}_{b,c}\right)}{ \operatorname{nnz}\!\left(\{l_{b,j,c}=1\}\cup\{\tilde{l}^{(i)}_{b,j,c}=1\} \right)}\in\mathbb{R} \tag{18a}\] \[\text{with}\ \ \tilde{l}^{(i)}_{b,j,c}=\begin{cases}1&\text{if }c= \operatorname{arg\,max}_{k}\{\tilde{p}^{(i)}_{b,j,k}\}\\ 0&\text{otherwise}\end{cases}\quad\text{forming}\ \ \tilde{l}^{(i)}_{b,j}=[\tilde{l}^{(i)}_{b,j,c}]_{c \in\mathbb{L}}\] (18b) \[\text{and}\ \ \mathbb{M}^{(i)}_{b,c}=\left[\left\{l_{b,j,c}=1, \tilde{l}^{(i)}_{b,j,c}\neq 1\right\}\cup\left\{l_{b,j,c}\neq 1,\tilde{l}^{(i)}_{b,j,c}=1\right\}\right]\in\{0,1\}^{|\mathbb{T}_{b}|} \tag{18c}\]

being the set of mispredictions defined over the discrete hypercube \(\{0,1\}^{|\mathbb{T}_{b}|}\). Also, \(\operatorname{nnz}\!\left(\mathbb{M}^{(i)}_{b,c}\right)\) was the number of nonzero elements of the binary set \(\mathbb{M}^{(i)}_{b,c}\). To form the convex continuous surrogate of the JD, first \(\mathbb{M}^{(i)}_{b,c}\in\{0,1\}^{|\mathbb{T}_{b}|}\) should be replaced by a nonnegative real-valued misprediction vector \(\mathbf{m}^{(i)}_{b,c}=\left[m^{(i)}_{b,j,c}\right]_{j}\in\mathbb{R}^{| \mathbb{T}_{b}|}_{\geq 0}\). Then, the surrogate should be found in \(\mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\). This search was NP-hard unless the JD was submodular. According to Proposition 11 in [13], the set function \(\text{JD}:\{0,1\}^{|\mathbb{T}_{b}|}\longmapsto\mathbb{R}\) was submodular. That is,

\[\forall\mathbb{M}_{1},\mathbb{M}_{2}\in\{0,1\}^{|\mathbb{T}_{b}|}:\quad\text{ JD}\!\left(\mathbb{M}_{1}\right)+\text{JD}\!\left(\mathbb{M}_{2}\right)\geq \text{JD}\!\left(\mathbb{M}_{1}\cup\mathbb{M}_{2}\right)+\text{JD}\!\left( \mathbb{M}_{1}\cap\mathbb{M}_{2}\right). \tag{19}\]

Under this condition, the convex closure of \(\text{JD}:\{0,1\}^{|\mathbb{T}_{b}|}\longmapsto\mathbb{R}\) in \(\mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\) was tight and continuous and could be computed in a polynomial time. This convex closure was called the Lovasz extension and was given in [14, 15] as

\[\begin{split}\overline{\text{JD}}:\ \ \mathbf{m}^{(i)}_{b,c}\in \mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\longmapsto\left[\frac{1}{|\mathbb{T}_{b}|} \sum_{j\in\mathbb{T}_{b}}m^{(i)}_{b,j,c}\cdot g_{j}\big{(}\mathbf{m}^{(i)}_{b, c}\big{)}\right]\in\mathbb{R}\\ \text{with}\ \ g_{j}\big{(}\mathbf{m}^{(i)}_{b,c}\big{)}=\text{JD}\! \left(\{u_{1},\cdots,u_{j}\}\right)-\text{JD}\!\left(\{u_{1},\cdots,u_{j-1}\} \right)\end{split} \tag{20}\]

being the \(j^{\text{th}}\) element of the gradient \(\mathbf{g}\big{(}\mathbf{m}^{(i)}_{b,c}\big{)}\) and \(\{u_{1},\cdots,u_{|\mathbb{T}_{b}|}\}\) denoting a permutation of the elements of \(\mathbf{m}^{(i)}_{b,c}=\left[m^{(i)}_{b,j,c}\right]_{j}\) in descending order, i.e. \(\left[\mathbf{m}^{(i)}_{b,c}\right]_{u_{1}}\geq\cdots\geq\left[\mathbf{m}^{(i) }_{b,c}\right]_{u_{|\mathbb{T}_{b}|}}\). Thus, the \(\overline{\text{JD}}\!\left(\mathbf{m}^{(i)}_{b,c}\right)\) was a weighted average of the elements of the misprediction vector \(\mathbf{m}^{(i)}_{b,c}\in\mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\) with the weights being the elements of the first derivative (gradient) of \(\overline{\text{JD}}\) with respect to \(\mathbf{m}^{(i)}_{b,c}\in\mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\). This way, the Lovasz extension \(\overline{\text{JD}}\) interpolated JD in \(\mathbb{R}^{|\mathbb{T}_{b}|}_{\geq 0}\setminus\{0,1\}^{|\mathbb{T}_{b}|}\) while having the same values as JD on \(\{0,1\}^{|\mathbb{T}_{b}|}\)[11, 11].

For a binary classification, the misprediction vector \(\mathbf{m}^{(i)}_{b,c}=\left[m^{(i)}_{b,j,c}\right]_{j}\in\mathbb{R}^{| \mathbb{T}_{b}|}_{\geq 0}\) was given by \(m^{(i)}_{b,j,c}=\max\!\left[(1-z^{(i)}_{b,j,c}\cdot l_{b,j,c}),\ 0\right]\) with \(\mathbf{z}^{(i)}_{b,j}=\left[z^{(i)}_{b,j,c}\right]_{c\in\mathbb{L}}\) being the network's outputs (before the softmax function) at the \(i^{\text{th}}\) iteration for the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\text{train}}\). This misprediction vector resulted in a convex piecewise linear surrogate called the Lovasz hinge loss [13].

For a multiclass classification, the misprediction vector \(\mathbf{m}^{(i)}_{b,c}=\left[m^{(i)}_{b,j,c}\right]_{j}\in\mathbb{R}^{| \mathbb{T}_{b}|}_{\geq 0}\) was formed from the classification posteriors \(\hat{\mathbf{p}}^{(i)}_{b,j}=\left[\hat{p}^{(i)}_{b,j,c}\right]_{c\in\mathbb{L}}\) produced by the softmax function in (7). This misprediction vector resulted in a convex continuous surrogate with regard to the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\text{train}}\) and the class \(c\in\mathbb{L}\) in the iteration \(i\in\{1,\cdots,n_{\text{it}}\}\). Thus, for the classification over \(n_{\text{clas}}=|\mathbb{L}|\) classes, the overall loss was an average of these class-specific surrogates. This
overall loss was called the Lovasz-Softmax loss and was given in [Berman 2018] as

\[\mathcal{L}_{\mathrm{LS}}(\mathbf{\hat{P}}_{b}^{(i)},\mathbf{L}_{b})= \frac{1}{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{c\in\mathbb{L}}\sum_{j\in \mathbb{T}_{b}}m_{b,j,c}^{(i)}\cdot g_{j}\big{(}\mathbf{m}_{b,c}^{(i)}\big{)} \tag{21}\] \[\text{with}\ \ \mathbf{m}_{b,c}^{(i)}=[m_{b,j,c}^{(i)}]_{j}\in \mathbb{R}_{\geq 0}^{|\mathbb{T}_{b}|}\ \ \text{and}\ \ m_{b,j,c}^{(i)}=\begin{cases}1-\hat{p}_{b,j,c}^{(i)}&\text{if }c=l_{b,j,c}\\ \hat{p}_{b,j,c}^{(i)}&\text{otherwise}\end{cases}\in(0,1).\]

The computation of the Lovasz extension \(\overline{\mathrm{JD}}\) in (20) implied to sort the elements of \(\mathbf{m}_{b,c}^{(i)}=[m_{b,j,c}^{(i)}]_{j}\in\mathbb{R}_{\geq 0}^{|\mathbb{T}_{b }|}\) and to call the JD with the permutation order. The sort had a complexity of \(\mathcal{O}\big{(}|\mathbb{T}_{b}|\cdot\log(|\mathbb{T}_{b}|)\big{)}\) and the call had a complexity of \(\mathcal{O}(|\mathbb{T}_{b}|)\). However, by keeping a track of the cumulative number of false positive and false negative mispredictions, the complexity of the call could be amortized to \(\mathcal{O}(1)\). That is, in each iteration, instead of computing the gradient from scratch only the gradient got updated. In this case, the overall complexity of computing (20) become \(\mathcal{O}\big{(}|\mathbb{T}_{b}|\cdot\log(|\mathbb{T}_{b}|)\big{)}\). The procedure of computing the gradient of the Lovasz-Softmax loss in (21) was given by Algorithm 1 in [Berman 2018].

The convexity and the differentiability of the Lovasz-Softmax loss in (21) allowed to use it as an objective function for optimizing a discriminative neural network classifier by a gradient descent optimizer with backpropagation. Also, the operations involved in its computation were differentiable and implementable on graphics processing units (GPUs).

### Baseline Architecture

Each convolutional layer of a neural network could extract features of a certain resolution while being capable of downsampling or reducing the spatial resolution by using an appropriate stride. These allowed to learn hierarchical (multiresolution) features by cascading multiple convolutional layers. The opposite of a convolutional layer was a transposed convolutional or a deconvolutional layer of similar feature learning capability but an inherent upsampling or increase of the spatial resolution. By following the convolutional layers with the deconvolutional layers an **encoder-decoder** architecture was obtained. The encoder was a downsampler, a compressor, or a contractor performing **analysis**. The decoder was an upsampler, a decompressor, or an expander performing **synthesis**. Each encoder/decoder was composed of multiple stages. Each **stage** processed features of a certain resolution through one or more convolutional/deconvolutional layers and then downsampled/upsampled its newly computed features to the next resolution. To avoid loss of information due to the downsampling, in each encoder stage, the number of the newly computed features got multiplied by the downsampling rate. Conversely, in each decoder stage, the number of the newly computed features got divided by the upsampling rate.

A widely used neural network of such an encoder-decoder architecture was the U-net. As the inputs passed through its encoder stages, the progressively expanding receptive fields of its convolutional layers increased the abstraction and the context of its extracted features. Thus, at the end of the encoder or bottom of the \(\mathbf{U}\), features of minimum resolution but maximum abstraction and context were obtained. The spatial resolution of these features got reconstructed by passing

Figure 2: Downsampling (left) and upsampling (right) in the V-net.
them through the deconvolutional layers of the decoder stages and combining them with original higher resolution features. The original features were directly obtained from the corresponding encoder stage through a skip connection. That is, features extracted by each encoder stage got forwarded to the corresponding decoder stage to compensate information loss due to the downsampling. This feature forwarding could enhance the delineation of boundaries between different classes and sped up the convergence of the optimization. At the end of the decoder, the resulting feature maps had a resolution and size like the input of the network. A weighted average of these feature maps combined them into the desired number of classes. This was done by passing them through a convolutional layer of \(1\times 1\times 1\) kernel size, \(0\) padding, and stride of \(1\) in each dimension. As given by (7), the resulting network's outputs got then passed through a softmax function to produce the estimated classification posteriors for the samples [Ronneberger 2015].

The downsampling and the upsampling of the U-net made it a hierarchical architecture capable of capturing, analyzing, and synthesizing features at different spatial resolutions. This way, the U-net could automatically extract local and contextual patterns. The local patterns got captured by the shallower layers and the contextual patterns by the deeper layers of a larger receptive field. At the end, the decoder synthesized (gathered and assembled) the local (high resolution) and the contextual (low resolution) features into the final classification. These enabled a localization as well as an accurate classification in any domain of any size and thus made the U-net a breakthrough for end-to-end optimizations. Moreover, making all the operations of the U-net 3D allowed to apply it to 3D volumetric domains. The 3D U-net got enhanced by making its encoder stages **residual**. That is, the input of each encoder stage got added to its output. This could mitigate vanishing gradients and speed up the convergence of the optimization [He 2016a]. In addition, the 3D U-net could learn 3D volumetric structures out of sparsely annotated 2D slices. This allowed to use it in a semi-automated annotation process as well as a fully automated 3D detection [Cicek 2016, Rakhlin 2018].

In the 3D U-net, each downsampling/upsampling had a factor of 2 and was done through a max-pooling/unpoolig over a \(2\times 2\times 2\) kernel with a stride of \(2\) in each dimension. Also, each convolutional layer applied \(0\) padding. Thus, the valid part of each feature map at the output of each convolutional layer had a smaller size than its input feature map. In addition, the 3D U-net learned the residual functions only in its encoder stages. In a so-called V-net, the 3D U-net become **fully convolutional** by applying each downsampling/upsampling through a convolutional/deconvolutional layer of a kernel size of \(2\times 2\times 2\), a 0 padding, and a stride of \(2\) in each dimension. To avoid loss of information, each downsampling doubled the number of feature maps. Conversely, each upsampling halved the number of feature maps. Figure 2 shows the downsampling and the upsampling in the V-net.

In contrast to the max-pooling/unpoolig operations, the convolution/deconvolution-based downsampling/upsampling was reversible and differentiable. These allowed to backpropagate each downsampling/upsampling without needing to store its inputs per sample and iteration. This way, the memory footprint of the V-net become much less than the 3D U-net while the analysis and comprehension of its internal process got simplified. Moreover, each convolution of the V-net applied an appropriate padding to make the feature maps at its output of the same size as its input. Furthermore, the V-net learned the residual functions not only in the encoder stages but also in the decoder stages. This further boosted its performance and sped up its optimization [Milletari 2016]. This way, the 3D U-net or the V-net got widely used in many applications [Rakhlin 2018, Li 2022]. Accordingly, we resorted to an end-to-end optimization of the 3D fully convolutional and residual V-net for our implementations and evaluations. For this, we tailored the number and

\begin{table}
\begin{tabular}{c|c c|c}
**Stage** & \multicolumn{2}{c|}{**Receptive Field**} & **Size of Feature Maps** \\ \hline \(1\) & \(5\times 5\times 5\) & \(551\times 551\times 551\) & \(128\times 352\times 256\) \\ \(2\) & \(22\times 22\times 22\) & \(546\times 546\times 546\) & \(64\times 176\times 128\) \\ \(3\) & \(72\times 72\times 72\) & \(528\times 528\times 528\) & \(32\times 88\times 64\) \\ \(4\) & \(172\times 172\times 172\) & \(476\times 476\times 476\) & \(16\times 44\times 32\) \\ \(5\) & \(372\times 372\times 372\) & \(372\times 372\times 372\) & \(8\times 22\times 16\) \\ \end{tabular}
\end{table}
Table 1: The receptive fields and the sizes of the feature maps at different stages of the V-net.
the sizes of the feature maps and the kernels of the convolutional/deconvolutional layers to our volumetric fat-water images. Also, through the network, we processed the data in an N\(\times\)D\(\times\)H\(\times\)W\(\times\)C format with N=\(\|\mathbb{T}_{b}|\) being the number of the volumetric fat-water images in each batch, C being the number of the feature maps, D being the depth, H being the height, and W being the width of each feature map. We trained (optimized) the V-net by using a mini-batch-based gradient descent optimizer with backpropagation and a sufficiently large input volume to capture as much contextual information as possible. Due to the memory limitations of the used GPU, we could only include 2 volumetric fat-water images in each batch. Moreover, each volumetric fat-water image had 2 channels containing its voxelwise fat and water intensities. Accordingly, at the input of the network, N\(\times\)D\(\times\)H\(\times\)W\(\times\)C=\(2\times 128\times 352\times 256\times 2\).

Each encoder/decoder stage of the V-net extracted and learned features of a certain spatial resolution by using one to three 3D (volumetric) convolutional/deconvolutional layers. In our case, each of these layers had a kernel size of \(5\times 5\times 5\), a padding of \(2\), and a stride of \(1\) in each dimension. Also, regarding the size of our images and the sizes of the addressed objects (tissues) in our segmentations, we found **5 stages (resolution levels)** to be sufficient for our hierarchical feature learning. Table 1 shows the receptive fields and the sizes of the feature maps at different stages. As can be seen, the innermost (deepest) stage of the network could already capture the entire context of the input volume. This allowed to perceive the whole anatomy of interest and ensured access to enough contextual information for reliably classifying each voxel at the output of the neural network classifier.

Besides the convolutional/deconvolutional layers, each **residual** encoder/decoder stage normalized its feature maps and applied nonlinearities to them. Like the original V-net, we used a parametric rectified linear unit (PReLU) with a parameter \(a_{\mathrm{prelu}}\in\mathbb{R}_{\geq 0}\) for each nonlinear activation. The parameter \(a_{\mathrm{prelu}}\in\mathbb{R}_{\geq 0}\) controlled the outputs for negative inputs and thus was called the coefficient of leakage. It got optimized along with the main parameters (weights and biases) of the network. The normalization of the feature maps decoupled the lengths of the network's gradients from their directions. This could accelerate the convergence of the optimizations and thus allowed higher learning rates. It could also stabilize the optimizations by mitigating the internal covariate shift1, enhancing the robustness against the initializations, and smoothing the objective function. Moreover, it could penalize large network's weights and thereby reduce the overfitting or improve the generalization.

Footnote 1: changes of stochastic distributions of the inputs of each layer of the network due to the changes of the parameters of the previous layers

We modified the V-net by changing the type of the normalization from batch normalization [10] to instance (contrast) normalization [10]. The commonly used batch normalization was based on mini-batch statistics. That is, during the training, the mean and the variance of each feature map of each batch got learned across all the dimensions (D, H, W) and all the N members of the batch to normalize (remove bias and scale of) the corresponding feature map in the evaluation phase. The instance normalization took a similar approach. However, it computed the mean and the variance of each feature map of each batch only across the dimensions (D, H, W). In case of having a small

Figure 3: Different normalization techniques applied to a feature map of size N\(\times\)D\(\times\)H\(\times\)W\(\times\)C with N denoting the number of batches, C denoting the number of channels, and D\(\times\)H\(\times\)W denoting the spatial dimensions. In each case, the blue voxels got normalized by the same mean and variance aggregated across them.
batch size, like our case, the exponential moving averages of the mean and the variance of each feature map of each batch had strong fluctuations across the training iterations. This was due to the poor statistical power of the small batch and thereby made the batch normalization ineffective. In this case, the instance normalization was more effective and consistent [16]. Other varieties of the normalization were the layer and the group normalization [20]. Figure 3 shows their differences to the batch and the instance normalization.

We also modified the V-net by changing the order of operations in each **residual** encoder/decoder stage. Instead of the convention of applying the normalization between the convolution/deconvolution and the nonlinear activation, as suggested in [14], we applied a full preactivation normalization and removed after-addition activation. Figure 4 compares the new and the original orders of the operations of a residual encoder/decoder stage comprising 2 convolutional/deconvolutional layers. The advantage of the new order was that it made the overall nonlinear function of each stage a real identity mapping. This enabled a direct and clean propagation of signals from one stage to another stage in both forward and backward directions. Other kinds of skip connections which involved a sort of scaling (like the Dropout), gating, or convolution/deconvolution on the signal path could hamper a clean propagation of the information and thus lead to optimization problems. Moreover, the new order could improve the generalization of the network's model by reducing its overfitting. That is, it increased the error on seen (training) samples but reduced the error on unseen (validation or test) samples. Furthermore, in the original order, addition of the shortcut to the normalized signal made the overall signal at the input of the last nonlinear activation unnormalized. However, in the new order, the signal at the input of each nonlinear activation was normalized. Figure 5 shows the described V-net architecture.

To mitigate **overfitting** and the **imbalanced class-sample distribution** of the training samples, **attention mechanisms** got proposed. These methods aimed to focus the attention of the network's parameters on important (foreground) minority classes. This attention could reduce the training samples to an **effective subset** of a lower unbalancedness than the original set. It could also vanish the redundant or irrelevant network's parameters by suppressing feature activations in irrelevant regions of the classification domain. These in turn reduced the overfitting and sped up the convergence of the network's optimization. The attention could be stimulated by **incorporating priors** into the optimization process and/or **modifying the network's architecture**. Neither the cross entropy-based nor the metric-based losses, defined in subsection 1.3, could accommodate the priors of the samples. Consequently, the attention mechanisms were restricted to architectural modifications.

Trainable (optimizable) attention mechanisms were categorized as **hard** or **soft**. The hard attention mechanisms iteratively cropped a region of interest through a Monte Carlo sampling optimized by a reinforcement learning. These

Figure 4: A residual encoder/decoder stage comprising 2 convolutional/deconvolutional layers (a) with the new (b) and the original order (c) of the operations.
sampling-based updates were indifferentiable and thus hard to optimize. The soft attention mechanisms involved a differentiable model composed of real-valued parameters. Thus, they could be optimized through a gradient descent optimizer with backpropagation. The output of the soft attention model for each feature map was a probabilistic map

Figure 5: Schematic of the 3D fully convolutional and residual V-net with the encoder and the decoder stages on its left and right side, respectively.
called **attention map**. In an additive or a multiplicative attention mechanism this map got computed by adding or multiplying the filtered feature map(s) by a filtered gating map, respectively. If the attention map was commuted by a convolutional neural network (CNN), then each filter was a convolutional layer. The attention mechanism turned into a self-attention if the gating maps were produced internally. The elementwise multiplication or addition of each attention map with its corresponding feature map highlighted salient features for the classification. This enabled an attention-based **feature pooling or pruning**. If the gating maps brought contextual information, then the feature pooling was with regard to the contextual dependencies of the features. Besides mitigating the overfitting and the imbalanced class-sample distribution of the training samples, the attention-based feature pooling could enhance the sensitivity, the prediction accuracy, and the robustness of the neural network classifier. A commonly used architecture for soft attention was a region proposing feed-forward CNN. A bottleneck of this approach was its excessive and redundant use of the model's parameters and features. This could increase the overall optimization overhead and the overfitting before the convergence of the optimization could realize any attention for a possible reduction of the network's parameters [20].

As mentioned earlier, the U-net and the V-net were capable of extracting (analyzing) and reconstructing (synthesizing) multiresolution (multiscale) features. This was done by extracting coarser features through downsampling the feature maps across the encoder stages and then reconstructing finer (higher resolution) features across the decoder stages. To this end, the receptive field at the coarsest resolution was to be large enough to capture all the contextual information highlighting the overall category and location of the foreground classes. After the localization, the finer (higher resolution) features delineated boundaries between different classes more precisely. These altogether allowed to capture large shape and size variations in the classification domain and thus improved the classification accuracy.

The reconstruction of the finer (higher resolution) features in each decoder stage was with the help of the features extracted by the corresponding encoder stage at the same spatial resolution. This feature forwarding reduced redundant and repeated computation of the features and thus enhanced efficiency in the usage of the computational power and memory. The plain skip connection of the feature forwarding path could be replaced by an **attention gate** realizing an **attention-based feature pooling**. This pooling vanished redundant features right before the concatenation of the original features with the reconstructed features. This way, it could suppress irrelevant regions in the classification domain by vanishing redundant network's perceptrons. This in turn reduced the overfitting of the network and the unbalancedness of the samples' distribution seen at the time of its training (optimization). Furthermore, the computational overhead of such an attention gate was much lower than the region proposing CNN. This and the reduction of the network's parameters could reduce the computational complexity of the optimizations and speed up their convergence [20].

A promising self-attention mechanism for integration into each feature forwarding path of the U-net or the V-net was a grid-based gating module. In this approach, each gating map was not fixed across the elements of its corresponding feature maps for which the attention maps were to be computed. Instead, it was a feature map of a lower (coarser) resolution already generated by the network itself. This way, the resulting attention maps were grid-based (i.e. variable across the elements of the feature maps) and could thus highlight salient features with respect to local patterns. The gating based on the feature maps of a lower (coarser) resolution allowed to consider a bigger context in the feature pooling and thereby disambiguated irrelevant and noisy features. Moreover, the grid-based gating module eliminated the need to an external explicit region proposing CNN by implicitly proposing soft (probabilistic) map of the target structures on the fly. This attention mechanism could be trained from scratch to focus on the target structures of varying shapes and sizes without additional supervision. Its filters (linear transformations) downweighted the gradients from irrelevant regions and could thus be implemented through convolutional layers filtering the network's activations in both forward and backward passes [20].

In [20, 17, 16], to reduce the number of the parameters and the computational complexity of the attention gates, each filter was a convolutional layer of \(0\) padding and \(1\times 1\times 1\) kernel size, i.e. without any spatial support. To downsample the input feature maps of each attention gate to the resolution of its gating maps, the convolutional filters of the feature maps had a stride of \(2\) in each dimension. Moreover, each attention gate handled a
binary classification and thus computed a common attention map for all the feature maps at its input. To this end, the downsampling convolutional filters of the feature maps linearly transformed them to an intermediate number of feature maps denoted by C'. Also, the convolutional filters of the gating maps linearly transformed them to C' intermediate maps. The intermediate feature/gating maps were to be more semantically discriminative than the original feature/gating maps in localizing the target structures. Thus, the number C' was a resolution-specific hyperparameter and needed to be optimized for each attention gate separately. Then, according to an **additive attention** mechanism, the intermediate downsampled feature maps got added to the intermediate gating maps and then passed through a nonlinear rectified linear unit (ReLU), a \(1\times 1\times 1\) convolutional layer of \(0\) padding and a stride of \(1\), and a nonlinear Sigmoid layer to form the attention map for all the input feature maps. This attention map had a lower resolution than the input feature maps and thus was upsampled by a grid-based trilinear interpolation to the same resolution as the input feature maps. In comparison to a multiplicative attention, the additive attention was more computationally demanding but more effective in enhancing the classification accuracy.

To handle a multiclass classification over \(n_{\mathrm{class}}=|\mathbb{L}|\) classes, we modified the aforementioned gating module by replacing the nonlinear Sigmoid function with a nonlinear Softmax function. Also, after the ReLU operation, the \(1\times 1\times 1\) convolutional layer did not map the outputs of the ReLU to one channel rather to the number of feature maps at the input of the gating module. That is, instead of computing one common attention map for all the input feature maps, we computed an attention map for each feature map separately and independently from other feature maps. Furthermore, to simplify the network's optimization we eliminated the resolution-specific hyperparameter C' defining the number of the intermediate feature/gating maps. To this end, the \(1\times 1\times 1\) convolutional layer directly applied to the input feature maps transferred them to the number of channels already existing in the input gating maps. This in turn eliminated the \(1\times 1\times 1\) convolutional layer directly applied to the input gating maps and thus further simplified the architecture of the gating module. Figure 6 compares the original gating module with our proposed one and Figure 7 shows the V-net architecture with such a gating module in each of its feature forwarding paths.

To reduce the overfitting of the baseline architectures to the seen (training) samples and thereby improve the generalization (predictive performance on unseen samples), we applied Dropout to every perceptron (node) of these architectures. This technique had a common root with a Bayesian neural network which, as described in subsection 1.2, was an ensemble of plain neural networks. In the training (optimization) phase, the Dropout dropped some of the perceptrons (nodes) of the network by vanishing their incoming and outgoing weights. The keep (retention) probability of each perceptron (node) was the occurrence probability of a Bernoulli distributed random variable. This probabi

Figure 6: Schematic of the original (**upper row**) and the proposed (**lower row**) grid-based gating module with \(\underline{\mathbf{F}}\) denoting the tensor of the input feature maps, \(\underline{\mathbf{G}}\) denoting the tensor of the gating maps, \(\underline{\mathbf{A}}\) denoting the tensor of the attention maps, \(\underline{\mathbf{E}}^{\prime}\) denoting the tensor of the output feature maps, and each blue box depicting a convolutional layer.
tunable hyperparameter indicating the confidence (inverse of the variance) of the node's estimations. We considered a common retention probability for all the perceptrons (nodes) of each encoder/decoder stage of the baseline architectures. For the \(s^{\text{th}}\) encoder/decoder stage, this probability was denoted by \(p_{s}\in[0,1]\). In the test phase, all the perceptrons (nodes) of the network were kept. However, the outgoing weights of each node got multiplied by its retention probability (\(p_{s}\)). The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in training phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase. The weights of each node are multiplied by the number of neurons in training phase phase. The weights of each node are multiplied by the number of neurons in training phase phase. The weights of each node are multiplied by the number of neurons in the training phase phase phase. The weights of each node are multiplied by the number of neurons in training phase phase phase. The weights of each node are multiplied by the number of neurons in training phase phase phase. The weights of each node are multiplied by the number of neurons in training phase phase phase. The weights of each node are multiplied by the number of neurons in training phase phase phase phase. The weights of each node are multiplied by the number of neurons in training phase
optimized during the hyperparameter optimization. The Dropout was shown to be superior to other regularization techniques such as the weight decay which penalized the weights of large \(l_{2}\) norms. This superiority come at the cost of a higher number of iterations for convergence of the optimizations (Srivastava 2014; Gal 2015; Jospin 2022).

## 2 Outline of Contributions

All the **metric-based** losses introduced in subsection 1.3 were independent of the class-sample distribution of the training samples and could thus enhance the generalization (predictive performance on unseen samples) of a neural network trained (optimized) with them. However, the metrics involved in those losses were binary classification metrics. This implied to decompose a multiclass classification into a series of one-vs-all classifications and then form its overall loss from an average of the one-vs-all losses. This was observable in the definition of the DICE loss in (16) and the Lovasz-Softmax loss in (21).

The averaging across the classes could naturally lead to a bias towards the dominant classes, i.e. classes of more samples. This bias could not be mitigated by a weighting mechanism such as the ones incorporated in the distribution-based losses introduced in page 8 and page 8. The reason was that such a weighting could diminish the false positive mispredictions on dominant classes and could thus mislead the optimization. Moreover, if a class was absent in both the reference labels and the predicted labels, then \(\text{DICE}=\text{JI}=1\) and \(\text{JD}=0\).

All the **distribution-based** losses introduced in subsection 1.3 were based on the cross entropy and had a common root with the variational free energy (VFE) of a retrospective active inference. These losses fitted the network's model to the class-sample distribution of the training samples and could thus compromise the network's generalization when the distribution of unseen (validation or test) samples differed from the distribution of the seen (training) samples. However, as described in page 8 and page 8, these losses could reduce the classification biases towards the dominant classes by weighting each class's term with regard to its number of samples or importance. In spite of this capability, there existed no optimal weighting which could be incorporated into the cross entropy-based losses to make them equivalent to any of the metric-based losses. Thus, to benefit from the advantages of the cross entropy-based and the metric-based losses while mitigating their drawbacks, a combination of them was used. Alternatively, to reduce the overfitting and thus to improve the generalization of the cross entropy-based losses, additional co-training with augmented training samples got conducted. Also, to reduce the classification biases towards the dominant classes, the false positive mispredictions of the network trained with the metric-based losses got post-corrected by using morphological operations (Isensee 2018; Bertels 2019; Jadon 2020; Chen 2022).

Despite of some improves, all the aforementioned schemes imposed extra overheads to the training or predictions of the neural networks. In addition, the augmentation of the training samples obtained from images was mostly done on the fly by applying gamma (luminance) modifications, mirroring, random scaling, random rotation, and random elastic deformation1 to the original images. These techniques could not be easily applied to medical images where pathological alterations should be differentiated from the augmentations. Moreover, none of the aforementioned schemes could completely mitigate the overfitting of a large network to a limited number of the training samples or the classification biases towards the dominant classes. Furthermore, none of the described losses could **incorporate priors** or **handle errors or uncertainties in the reference labels of the training samples**(Lo 2021). Errors in the reference labels of the training samples could arise from human errors in the manual annotations of the training samples and images or the errors induced by noise and artifacts. Uncertainties and ambiguities in the reference labels of the training samples could stem from similar features and textures of different classes. These similarities not only confused the manual annotators but also the neural network relying on those features and textures for learning boundaries between different classes.

Footnote 1: The elastic deformations were obtained from a B-spline interpolation over a grid of control points on a dense deformation field.

To mitigate the aforementioned bottlenecks, we proposed
1. a novel algorithm, based on the generalized (multinomial) Kelly criterion for optimal betting, to recompute the reference labels of the training samples by using their priors and the currently estimated classification posteriors on the network;
2. a novel objective function, based on the expected free energy (EFE) of a prospective active inference, with the capability of * incorporating prior probabilities of the training samples to focus the attention of the neural network on important but minority foreground classes and thereby reshape the effectively seen distribution for a reduction of the class-sample unbalancedness, the overfitting, and the classification biases towards the dominant classes; * representing the _precision_ and _recall_ metrics by its terms to enhance the robustness of the network's optimization against the class-sample unbalancedness;
3. a process to integrate the proposed algorithm and the proposed objective function into a mini-batch-based gradient descent optimizer with backpropagation.

The proposed algorithm for recomputing the reference labels was listed in Algorithm 1. This algorithm calculated a **set of candidate labels** for each training sample from its prior and currently estimated posterior probabilities on the network. This algorithm resulted from our reformulation of the generalized (multinomial) Kelly criterion for optimal betting on multiple horses in a horse race. This reformulation cast the generalized Kelly criterion into a multiclass classification problem by interpreting each training sample as a bettor, each class as a horse, and each iteration of the network's optimization as a horse race. Then, the classification prior of the training sample with regard to each class become the win probability of the corresponding horse. The classification posterior currently estimated by the network for the training sample with regard to the same class become the belief probability of the corresponding horse. The proposed sets of candidate labels got then plugged into the proposed objective function to form the current loss for an update (optimization) of the network's parameters in the current iteration. Thus, instead of a reference label, a set of candidate labels got considered for each training sample in each iteration.

This consideration allowed to mitigate the aforementioned uncertainties and ambiguities in the labels generated from manual annotations in the presence of noise, artifacts, and similar features or textures of different classes. In other words, the sets of candidate labels could handle possible overlaps between different classes and thus enhanced the reliability and the flexibility of the neural network's optimization. More specifically, these sets could help a gradient descent optimizer to escape from local optimums caused by the original reference labels. Moreover, if the reference labels of some training samples were missing, then their candidate labels could still be computed from their priors and posteriors. This **semi-supervised optimization** was of particular importance in the applications where the manual annotations of the reference labels were costly and cumbersome.

Our proposed Algorithm 1 for finding the candidate labels aimed to minimize the objective function of the generalized Kelly criterion. This minimized function was given by (36) and was indeed the **expected complexity** term of the EFE of a prospective active inference. That is, the objective function of the generalized Kelly criterion was a tight upper bound of the expected complexity of the EFE. The EFE was given by (4) and was composed of an **expected complexity** term plus an **uncertainty** term. As described in subsection 1.1, the minimization of the expected complexity was equivalent to the maximization of the reward. The reward maximization was also a goal of the Kelly criterion and could thus be partially fulfilled by finding the candidate labels through the proposed Algorithm 1.

More specifically, from the prior (win) and the posterior (belief) probabilities of each training sample (bettor), the generalized Kelly criterion computed optimal allocation fractions of the bettor's asset for betting on the candidate classes (horses)1. These allocation fractions maximized the geometric average of the growth rate of the bettor's asset or the **reward**. To further maximize the reward, the expected complexity of the EFE should be minimized further. This was doable by having enough information or maximizing the information gain, i.e. minimizing the **uncertainty** of the
EFE. Accordingly, to optimize a discriminative neural network classifier, we proposed a novel objective function based on the EFE of a prospective active inference. This function was given by (39) and was reversible and differentiable with respect to the outputs of every layer of the neural network. Thus, as described in subsection 1.2, it could be minimized by a gradient descent optimizer with backpropagation.

As explained in subsection 1.3, all the cross entropy-based losses were **distribution-based** and stemmed from the VFE given by (2) for a retrospective active inference. The VFE was **complexity** minus **accuracy**. The complexity reflected the overfitting of the neural network's model to the distribution of seen (training) samples and thus the variance of the predictions on unseen (validation or test) samples. The accuracy was inversely proportional to the bias (difference) of the predictions from their true values. Thus, the minimization of the VFE implied to minimize the complexity or the overfitting while maximizing the classification accuracy by minimizing the classification bias. This way, the VFE and the cross entropy-based losses addressed the bias-variance tradeoff of the classification problems without considering the unbalancedness of the class-sample distribution of the seen samples.

In contrast, the EFE given by (4) for a prospective active inference and thus our proposed objective function in (39) addressed the unbalancedness of the class-sample distribution of the seen (training) samples by representing the _precision_ and _recall_ metrics in their terms. The _precision_ and the _recall_ metrics were independent of the correct classification of unimportant majority samples (designated by true negatives) and instead focused on the correct classification of important minority samples (designated by true positives). This made them less sensitive than the other metrics to the imbalanced class-sample distributions and the classification biases towards the dominant classes.

As mentioned earlier, the minimization of the EFE or our proposed objective function implied to minimize the **expected complexity** and the **uncertainty**. The minimization of the expected complexity implied to maximize the reward and the reward was equivalent to the _recall_ (completeness or diversity). The minimization of the uncertainty implied to maximize the information gain or the _precision_ (exactness or confidence). This way, the EFE and our proposed objective function aimed to maximize the _precision_ and the _recall_ metrics. This allowed them to handle an imbalanced class-sample distribution while still being **distribution-based**[11, 12, 13].

Moreover, our proposed objective function could incorporate the prior probabilities of the training samples directly and indirectly. The indirect incorporation was through using the candidate classification labels computed from the priors and the posteriors of the training samples by the proposed Algorithm 1. This incorporation resulted in a grouping of the terms of the proposed objective function with regards to the candidate and noncandidate labels. More specifically, the priors or the posteriors of the noncandidate labels got summed together to form a collective prior or posterior for the noncandidate classes. This way, the noncandidate classes formed a collective class together and the neural network got enforced to find the boundary between each candidate class and the collective class of the noncandidates. In comparison to computing the boundaries between each pair of the classes, this grouping reduced the effective number of the classes and the boundaries needed to be computed. This in turn reduced the network's complexity and its overfitting to the seen (training) distribution and could thus enhance its generalization (predictive performance on unseen samples).

The direct incorporation of the prior probabilities of the training samples into the objective function of the network's optimization could focus the attention of the neural network on important but minority foreground classes. This could reshape the distribution effectively seen by the network during its optimization and could thereby reduce the class-sample unbalancedness, the overfitting, and the classification biases towards the dominant classes [12]. Similar effects could result from the architecture-based attention mechanisms described in subsection 1.4. That is, if no prior probabilities were provided, then **stronger posteriors** resulted from an **architecture-based attention mechanism** should help. In the baseline architecture described in subsection 1.4, an attention gate could be incorporated into each feature forwarding path between an encoder stage and its corresponding decoder stage. Without such a gate, the feature forwarding path was a plain skip connection.

Our proposed algorithm for finding the candidate labels and our proposed objective function for optimizing a discriminative neural network classifier got integrated into a mini-batch-based gradient descent optimizer with backpropagation
by using the process proposed in section 4. This process got evaluated against a similar process incorporating a representative of the cross entropy-based losses or a representative of the metric-based losses introduced in subsection 1.3. The representative of the cross entropy-based losses was the weighted focal loss. This loss comprised of a modulating factor and a weighting mechanism to alleviate classification biases towards the dominant classes of the training samples. The representative of the metric-based losses was the Lovasz-Softmax loss. Besides being smooth and differentiable, to the best of our knowledge, this loss was the only convex loss among the metric-based losses.

Accordingly, the evaluated losses were

1. the proposed objective function given by (39)
2. the weighted focal loss given by (14)
3. the Lovasz-Softmax loss given by (21).

These evaluations were on an end-to-end optimization of the baseline architecture described in subsection 1.4. For each case, the baseline architecture was once used without attention gates as depicted in Figure 5 and once used with the attention gates as depicted in Figure 7. Also, for (2) and (3) each training sample was accompanied by its reference (ground truth) label to fulfill the supervised nature of these objective functions. However, our proposed algorithm for finding the candidate labels and our proposed objective function got evaluated according to a fully supervised, a semi-supervised, and an unsupervised approach. These resulted in the training samples being

1. accompanied by their reference labels and their priors \(\rightarrow\) fully supervised
2. only accompanied by their reference labels \(\rightarrow\) semi-supervised
3. only accompanied by their priors \(\rightarrow\) semi-supervised
4. accompanied by neither their reference labels nor their priors \(\rightarrow\) unsupervised.

The unsupervised case only relied on the posteriors estimated by the neural network during its optimization and could thus be considered as a self-supervised case as well.

For the cases with the priors, the prior probabilities of the training samples could be computed by a multiatlas registration. If no prior probabilities were provided at the time of optimization (training), then uniform priors got assumed. If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vectorized reference label \(\mathbf{l}_{b,j}\) was the one-hot-encoding of its reference label \(l_{b,j}\in\mathbb{L}\) and was given by (8). If the reference labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were not provided at the time of optimization, then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vector \(\mathbf{l}_{b,j}\) was uniform and given by (9).

For each evaluation case, the main parameters and the hyperparameters of the baseline architecture got trained (optimized) to automatically segment \(n_{\mathrm{class}}=|\mathbb{L}|=8\) classes of vertebral bodies (VBs), intervertebral disks (IVDs), psoas major (PM) and quadratus lumborum (QL) muscles, epicardial adipose tissues (EpAT), pericardial adipose tissues (PeAT), cardiac perivascular adipose tissues (PvAT), and background on each volumetric fat-water image. To this end, the volumetric fat-water images got divided into a training and a test set. The training set formed the samples set \(\mathbb{T}_{\mathrm{train}}\) and got used to optimize the main parameters and the hyperparameters of the baseline architecture by each method. The test set formed the samples set \(\mathbb{T}_{\mathrm{test}}\) and got used to evaluate the classification performance of the baseline architecture after being fully optimized by each method. The training set was composed of samples accompanied by their reference labels and priors. The test set was composed of samples accompanied by their reference labels. The reference labels of the test samples were not fed to the neural network. They were rather compared against the corresponding labels predicted by the network to evaluate the classification performance of the network. The predicted label of each sample was the index of its maximum classification posterior estimated by the network.

Finally, our proposed optimization process was based on the generalized Kelly criterion for optimal betting and a prospective active inference. It addressed optimization of discriminative neural network classifiers with a feed-forward
architecture. Active inference-based optimizations could foster building highly flexible and generalizable generative models with and without memory. An example of a model with the memory was the one which could explain a partially observable Markov decision process. This model could be implemented by a recurrent or a long short-term memory network [14, 15, 16]. Accordingly, our proposed optimization process could be easily extended to generative or recurrent neural networks such as the networks in [1, 16, 17].

## 3 Application of the Kelly Criterion to Classification

The generalized (multinomial) Kelly criterion proposed optimal allocation fractions of a bettor's asset in betting on multiple horses in a horse race. Each horse had a win and a belief probability. The win probability was the chance of the horse to win the race. The belief probability was the collective belief of other bettors about the chance of the horse to win the race. Thus, for a specific bettor, an optimum betting strategy was to invest as much as possible on a horse of maximum win probability and minimum belief probability (minimum number of other bettors investing on it). This was based on the assumption that all the bettors followed the same strategy and the gain of a horse win got divided between all the bettors who have invested on it. Therefore, the lesser the belief probability was, the higher the paid gain to the investing bettor would be [13, 16].

To optimize a discriminative neural network classifier in a multiclass classification over \(n_{\mathrm{clas}}=|\mathbb{L}|\) classes by using the generalized Kelly criterion, we assumed

* every training sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) to be a bettor
* every class \(c\in\mathbb{L}\) to be a horse
* every iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\) of the optimization to be a round of horse race with its gambling competitions among the bettors (training samples)
* the win probability of each horse (class) \(c\in\mathbb{L}\) for each bettor (training sample) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) to be the prior probability \(a_{b,j,c}\in(0,1)\) estimated by another classifier1 Footnote 1: If no prior probabilities were provided then uniform priors got assumed.
* the belief probability of each class \(c\in\mathbb{L}\) for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) to be the classification posterior \(\hat{p}_{b,j,c}^{(i)}\in(0,1)\) estimated by the network in the current iteration \(i\).

It should be noted that in the betting, the win probabilities of the horses were shared across the bettors, but, in the classification, each sample had its own win probability for each class. Moreover, the interpretation of the estimated posteriors of the network as the belief probabilities might look counterintuitive because each sample (bettor) had no _other_ samples (bettors) to compete with. Thus the overall belief about a class (horse) could not be collected from other samples (bettors). Moreover, it was more tempting to select a class (invest on a horse) of maximum belief probability as this probability could be an indicator of the chance of the class (horse) to win. Our definition of the win probability and our counterintuitive definition of the belief probability could be explained under an **attention mechanism**.

On one hand, the selection of the classes (horses) of maximum win probability encouraged the network to focus on classes of confident (high) prior probabilities. In an image segmentation task conducted in a spatial domain, this implied to focus on important (relevant) regions highlighted by high prior probabilities in the image. On the other hand, the selection of the classes (horses) of minimum belief probability encouraged the network to focus on inconfident (low) posteriors and thus to improve its classification by tackling difficult examples.

In each iteration (race) \(i\), for each training sample (bettor) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\), the Kelly criterion proposed allocation fractions \(\hat{\mathbf{g}}_{b,j}^{(i)}=\left[\hat{g}_{b,j,c}^{(i)}\in[0,1]\right]_{c\in \mathbb{L}}\) of its asset for betting on \(n_{\mathrm{clas}}=|\mathbb{L}|\) classes (horses). If in the iteration (race) \(i\) the class (horse) \(c\in\mathbb{L}\) won, then the asset of \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) would be multiplied by \(\left[1-\sum_{k\in\mathbb{I}}\hat{g}_{b,j,k}^{(i)}+\frac{\hat{g}_{b,j,c}^{(i )}}{\hat{p}_{b,j,c}^{(i)}}\right]^{-1}\). We assumed that the outcomes of the iterations (horse races) were independent identically distributed (i.i.d.) random
variables. Thus, after \(i\) iterations, the geometric average of the growth rate of the asset of \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) with \(n_{c}^{(i)}\in[0,i]\) number of wins for each class \(c\in\mathbb{L}\) become

\[\eta_{b,j}^{(i)}=\prod_{c\in\mathbb{L}}\left[1-\sum_{k\in\mathbb{L}}\hat{g}_{b, j,k}^{(i)}+\frac{\hat{g}_{b,j,c}^{(i)}}{\hat{p}_{b,j,c}^{(i)}}\right]^{-n_{c}^{(i )}/i}\qquad i=\sum_{c\in\mathbb{L}}n_{c}^{(i)}. \tag{22}\]

By taking the \(\ln(\cdot)\) of both sides of (22), one obtained

\[\begin{split}\ln(\eta_{b,j}^{(i)})=\sum_{c\in\mathbb{L}}\frac{- n_{c}^{(i)}}{i}\cdot\ln\biggl{[}1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}+ \frac{\hat{g}_{b,j,c}^{(i)}}{\hat{p}_{b,j,c}^{(i)}}\biggr{]}\\ \lim_{i\to\infty}\frac{n_{c}^{(i)}}{i}=a_{b,j,c}\implies\lim_{i \to\infty}\ln(\eta_{b,j}^{(i)})=\sum_{c\in\mathbb{L}}-a_{b,j,c}\cdot\ln\biggl{[} 1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}+\frac{\hat{g}_{b,j,c}^{(i)}}{ \hat{p}_{b,j,c}^{(i)}}\biggr{]}.\end{split} \tag{23}\]

If the allocation fractions \(\mathbf{g}_{b,j}^{(i)}=\bigl{[}g_{b,j,c}^{(i)}\in[0,1]\bigr{]}_{c\in\mathbb{L}}\) proposed by the Kelly criterion for each sample (better) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) were **asymptotically optimum** over a long run \((i\to\infty)\), then they maximized the geometric average in (22). Due to the monotonic increase of the \(\ln(\cdot)\) function, the maximization of (22) was equivalent to the maximization of (23). This way, the asymptotically optimum allocation fractions were the maximizers of the averaged logarithms of the growth rate in (23). That is, \(\mathbf{g}_{b,j}^{(i)}=\operatorname*{arg\,max}_{\hat{\mathbf{g}}_{b,j}^{(i)} }\ \Bigl{[}\ln(\eta_{b,j}^{(i)})\Bigr{]}\) or

\[\mathbf{g}_{b,j}^{(i)}=\operatorname*{arg\,min}_{\hat{\mathbf{g}}_{b,j}^{(i) }}\ \Bigl{[}-\ln(\eta_{b,j}^{(i)})\Bigr{]}=\operatorname*{arg\,min}_{\hat{ \mathbf{g}}_{b,j}^{(i)}}\ \biggl{[}\sum_{c\in\mathbb{L}}a_{b,j,c}\cdot\ln\biggl{[}1-\sum_{k\in \mathbb{L}}\hat{g}_{b,j,k}^{(i)}+\frac{\hat{g}_{b,j,c}^{(i)}}{\hat{p}_{b,j,c}^ {(i)}}\biggr{]}\biggr{]}. \tag{24}\]

As detailed in [10], \(\hat{\mathbf{g}}_{b,j}^{(i)}=\bigl{[}\hat{g}_{b,j,c}^{(i)}\bigr{]}_{c\in \mathbb{L}}\in[0,1]^{n_{\rm clas}=|\mathbb{L}|}\) formed a convex set

\[\mathbb{G}_{b,j}^{(i)}=\left\{\hat{\mathbf{g}}_{b,j}^{(i)}\in[0,1]^{n_{\rm clas }=|\mathbb{L}|}\ \bigg{|}\ \biggl{[}1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}+\frac{\hat{g}_{b,j,c}^{( i)}}{\hat{p}_{b,j,c}^{(i)}}\biggr{]}>0\right\}\subseteq[0,1]^{n_{\rm clas}=| \mathbb{L}|} \tag{25}\]

which was an intersection of half spaces. Each half space was a side of a hyperplane. In addition, in the above optimization, \(\bigl{[}1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}\bigr{]}\in[0,1]\implies \sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}\in[0,1]\). That is, it was allowed to back a horse to win but not to lay a horse to lose. This condition constrained every \(\hat{\mathbf{g}}_{b,j}^{(i)}\in\mathbb{G}_{b,j}^{(i)}\) to a stricter convex set given by

\[\mathbb{G}_{b,j}^{(i)}=\left\{\hat{\mathbf{g}}_{b,j}^{(i)}\in\mathbb{G}_{b,j}^ {(i)}\ \right|\ \sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}\leq 1\ \ \text{and}\ \ \forall c\in\mathbb{L}:\hat{g}_{c,j}^{(i)}\geq 0\right\}\subseteq \mathbb{G}_{b,j}^{(i)}. \tag{26}\]

The definition of \(\ln(\eta_{b,j}^{(i)})\) in (23) showed that it was a finite linear combination of strictly concave logarithms with the coefficients being the priors \(\mathbf{a}_{b,j}=\bigl{[}a_{b,j,c}\in(0,1)\bigr{]}_{c\in\mathbb{L}}\). This way, the \(\ln(\eta_{b,j}^{(i)})\) become differentiable, strictly concave downwards, and of a unique maximum on the boundary of every bounded subset of \(\mathbb{G}_{b,j}^{(i)}\). Accordingly, to find the maximizers of \(\ln(\eta_{b,j}^{(i)})\) or the optimum allocation fractions \(\mathbf{g}_{b,j}^{(i)}=\bigl{[}g_{b,j,c}^{(i)}\in[0,1]\bigr{]}_{c\in\mathbb{L}}\), it was enough to only explore the boundaries of \(\mathbb{G}_{b,j}^{(i)}\subseteq\mathbb{G}_{b,j}^{(i)}\)[10]. This exploration (maximization) could be done by using the method of Lagrange multipliers and the Karush-Kuhn-Tucker (KKT) theory [1]. That is, instead of maximizing \(\ln(\eta_{b,j}^{(i)})\), we maximized

\[\gamma_{b,j}^{(i)}=\ln(\eta_{b,j}^{(i)})+\Bigl{[}\sum_{k\in\mathbb{L}}\lambda_{b,j,k}^{(i)}\cdot\hat{g}_{b,j,k}^{(i)}\Bigr{]}+\lambda_{b,j,0}^{(i)}\cdot \Bigl{[}1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i)}\Bigr{]} \tag{27}\]

with \(\bigl{\{}\lambda_{b,j,k}^{(i)}\in\mathbb{R}_{\geq 0}\bigr{\}}_{k=0}^{|\mathbb{L}|}\) being the Lagrange multipliers.
The KKT theory stated that every constrained maximizer of \(\ln(\eta^{(i)}_{b,j})\) was an unconstrained maximizer of \(\gamma^{(i)}_{b,j}\). The unconstrained maximization of \(\gamma^{(i)}_{b,j}\) was done through vanishing its gradient (derivatives) with respect to \(\hat{\mathbf{\beta}}^{(i)}_{b,j}=\left[\hat{\mathbf{\beta}}^{(i)}_{b,j,c}\in[0,1]\right]_ {c\in\mathbb{L}}\). That is,

\[\frac{\partial\gamma^{(i)}_{b,j}}{\partial\hat{g}^{(i)}_{b,j,c}}=\frac{-a_{b,j, c}+a_{b,j,c}/\hat{p}^{(i)}_{b,j,c}}{1-\sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}+ \hat{g}^{(i)}_{b,j,c}/\hat{p}^{(i)}_{b,j,c}}+\lambda^{(i)}_{b,j,c}-\lambda^{(i )}_{b,j,0}=0. \tag{28}\]

This resulted in the following KKT optimality constraints:

\[\text{if}\ \ \lambda^{(i)}_{b,j,c}\cdot\hat{g}^{(i)}_{b,j,c}=0 \implies\lambda^{(i)}_{b,j,c}=0\ \ \text{if}\ \ \hat{g}^{(i)}_{b,j,c}>0\] \[\text{if}\ \ \lambda^{(i)}_{b,j,0}\cdot\left[1-\sum_{k\in \mathbb{L}}\hat{g}^{(i)}_{b,j,k}\right]=0\implies\lambda^{(i)}_{b,j,0}=0\ \ \text{if}\ \ \sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}<1. \tag{29}\]

The allocation fractions \(\hat{\mathbf{\beta}}^{(i)}_{b,j}=\left[\hat{g}^{(i)}_{b,j,c}\in[0,1]\right]_{c\in \mathbb{L}}\) and the Lagrange multipliers \(\left\{\lambda^{(i)}_{b,j,k}\in\mathbb{R}_{\geq 0}\right\}_{k=0}^{|\mathbb{L}|}\) should fulfill (29) on the convex set \(\mathbb{G}^{\prime}_{\ b,j}\subseteq\mathbb{G}^{(i)}_{b,j}\). According to [11], the maximum of \(\ln(\eta^{(i)}_{b,j})\) under \(\sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}=1\) was less than its maximum under \(\sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}<1\). Thus, in (26), we replaced \(\sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}\leq 1\) with \(\sum_{k\in\mathbb{L}}\hat{g}^{(i)}_{b,j,k}<1\) and obtained \(\lambda^{(i)}_{b,j,0}=0\) from (29). For each sample (bettor) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\text{train}}\), the classes (horses) whose allocation fractions were nonzero were deemed to be **candidate** and formed the set \(\mathbb{L}^{(i)}_{b,j}\) with

\[\begin{split}\forall c\in\mathbb{L}^{(i)}_{b,j}\subseteq\mathbb{ L}:\ \hat{g}^{(i)}_{b,j,c}>0\ \ \text{and}\ \ \lambda^{(i)}_{b,j,c}=0\\ \forall c\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}:\ \hat{g}^{(i)}_{b,j,c}=0\ \ \text{and}\ \ \lambda^{(i)}_{b,j,c}\geq 0.\end{split} \tag{30}\]

Then, solving (28) under the above conditions gave

\[\forall c\in\mathbb{L}^{(i)}_{b,j}\subseteq\mathbb{L}:\ g^{(i)}_{b,j,c}=a_{b, j,c}-\hat{p}^{(i)}_{b,j,c}\cdot\frac{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b, j,k}}{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,k}} \tag{31}\]

\[\implies s^{(i)}_{b,j}=1-\sum_{c\in\mathbb{L}}g^{(i)}_{b,j,c}=1-\sum_{c\in \mathbb{L}^{(i)}_{b,j}}g^{(i)}_{b,j,c}=\overbrace{1-\sum_{c\in\mathbb{L}^{(i)} _{b,j}}a_{b,j}}^{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b,j,k}}+\frac{ \sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b,j,k}}{\sum_{k\in\mathbb{L}- \mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,k}}\cdot\sum_{c\in\mathbb{L}^{(i)}_{ b,j}}\hat{p}^{(i)}_{b,j,c}\]

\[=\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b,j,k}\cdot\left[1+\frac{\sum_{ c\in\mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,c}}{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{ b,j}}\hat{p}^{(i)}_{b,j,k}}\right]=\frac{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b,j,k}}{ \sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,k}} \tag{32}\]

\[\forall c\in\mathbb{L}^{(i)}_{b,j}\subseteq\mathbb{L}:\ s^{(i)}_{b,j}+\frac{g^{( i)}_{b,j,c}}{\hat{p}^{(i)}_{b,j,c}}=\frac{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}a_{b,j,k}}{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,k}}+\frac{a_ {b,j,c}}{\hat{p}^{(i)}_{b,j,c}}-\frac{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j} }a_{b,j,k}}{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}}\hat{p}^{(i)}_{b,j,k}}= \frac{a_{b,j,c}}{\hat{p}^{(i)}_{b,j,c}} \tag{33}\]

\[\forall c\in\mathbb{L}^{(i)}_{b,j}\subseteq\mathbb{L}\ \text{and}\ \ \forall l\in\mathbb{L}- \mathbb{L}^{(i)}_{b,j}:\ \ \frac{a_{b,j,l}}{\hat{p}^{(i)}_{b,j,l}}\leq s^{(i)}_{b,j}=\frac{\sum_{k\in\mathbb{L}- \mathbb{L}^{(i)}_{b,j}}a_{b,j,k}}{\sum_{k\in\mathbb{L}-\mathbb{L}^{(i)}_{b,j}} \hat{p}^{(i)}_{b,j,k}}<\frac{a_{b,j,c}}{\hat{p}^{(i)}_{b,j,c}}. \tag{34}\]

## 4 Proposed Objective and Process of Optimization

By using our classification-based formulation of the Kelly criterion in section 3 we proposed an objective function and a process for optimizing discriminative neural network classifiers. To be generic, we formulated the objective and the process in such a way that they could accommodate a fully supervised, a semi-supervised, or an unsupervised
optimization. In the fully supervised optimization, both the reference (ground truth) labels and the prior (win) probabilities of the training samples were provided at the time of optimization (training). In the semi-supervised optimization, either the reference labels or the prior (win) probabilities of the training samples were not provided at the time of optimization (training). In the unsupervised optimization, neither the reference labels nor the prior (win) probabilities of the training samples were provided at the time of optimization (training). If no prior probabilities were provided at the time of optimization (training), then uniform priors got assumed. If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vectorized reference label \(\mathbf{l}_{b,j}\) was a one-hot-encoding of its reference (ground truth) label \(l_{b,j}\in\mathbb{L}\) and was given by (8). If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were not provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vector \(\mathbf{l}_{b,j}\) was uniform and given by (9).

We denoted the vectorized reference labels, the fixed prior (win) probabilities, and the estimated posterior (belief) probabilities of the samples in the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) with the \(\left|\mathbb{T}_{b}\right|\times n_{\mathrm{class}}\) matrices of \(\mathbf{L}_{b}=\left[\mathbf{l}_{b,j}\right]_{j}=\left[l_{b,j,c}\right]_{j,c}\), \(\mathbf{A}_{b}=\left[\mathbf{a}_{b,j}\right]_{j}=\left[a_{b,j,c}\right]_{j,c}\), and \(\hat{\mathbf{P}}_{b}^{(i)}=\left[\hat{\mathbf{p}}_{b,j}^{(i)}\right]_{j}= \left[\hat{p}_{b,j,c}^{(i)}\right]_{j,c}\), respectively. Also, the allocation fractions estimated by the Kelly criterion for these samples formed a \(\left|\mathbb{T}_{b}\right|\times n_{\mathrm{class}}\) matrix denoted by \(\hat{\mathbf{G}}_{b}^{(i)}=\left[\hat{\mathbf{g}}_{b,j}^{(i)}\right]_{j}= \left[\hat{g}_{b,j,c}^{(i)}\right]_{j,c}\).

In each iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\) of optimizing a discriminative neural network classifier, we first found the set of **candidate** classification labels \(\mathbb{L}_{b,j}^{(i)}\subseteq\mathbb{L}\) for each sample (better) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\). To this end, we proposed Algorithm 1 by using (31), (32), (33), and (34). Through this algorithm, the set of candidate labels \(\mathbb{L}_{b,j}^{(i)}\subseteq\mathbb{L}\) got computed from the estimated posterior (belief) probabilities \(\hat{\mathbf{p}}_{b,j}^{(i)}=\left[\hat{p}_{b,j,c}^{(i)}\in(0,1)\right]_{c\in \mathbb{L}}\) and the fixed prior (win) probabilities \(\mathbf{a}_{b,j}=\left[a_{b,j,c}\in(0,1)\right]_{c\in\mathbb{L}}\) of the sample (better) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\).

The set \(\mathbb{L}_{b,j}^{(i)}\subseteq\mathbb{L}\) could contain multiple class labels or be empty. An empty set implied that the current posterior (belief) and the fixed prior (win) probabilities found no class label, even the reference label \(l_{b,j}\in\mathbb{L}\), to be reliable enough for the optimization of the neural network classifier. This could result in no further update of the posterior
(belief) probabilities in the following iterations. To avoid this standstill, at the end of the Algorithm 1, if \(\mathbb{L}_{b,j}^{(i)}=\emptyset\), then the reference label \(l_{b,j}\in\mathbb{L}\) of the sample (better) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) got inserted into it.

By extending (24) to all the samples in the batch \(\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\), one obtained

\[\mathbf{G}_{b}^{(i)}=\operatorname*{arg\,min}_{\mathbf{G}_{b}^{(i)}}\;\underbrace{ \frac{1}{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in \mathbb{L}}a_{b,j,c}\cdot\ln\!\left[1-\sum_{k\in\mathbb{L}}\hat{g}_{b,j,k}^{(i )}+\frac{\hat{g}_{b,j,c}^{(i)}}{\hat{p}_{b,j,c}^{(i)}}\right]}_{\mathcal{L}_{ \rm Kelly}(\mathbf{G}_{b}^{(i)})}. \tag{35}\]

However, the optimum allocation fractions \(\mathbf{G}_{b}^{(i)}=\left[\mathbf{g}_{b,j}^{(i)}\right]_{j}=\left[g_{b,j,c}^{ (i)}\right]_{j,c}\) had a closed form solution given by (31). This solution resulted in (32) and (33) and allowed to express

\[\min_{\mathbf{G}_{b}^{(i)}}\;\mathcal{L}_{\rm Kelly}(\mathbf{ \hat{G}}_{b}^{(i)})=\mathcal{L}_{\rm Kelly}(\mathbf{G}_{b}^{(i)})=\frac{1}{| \mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{c\in\mathbb{L} }a_{b,j,c}\cdot\ln\!\left[s_{b,j}^{(i)}+\frac{g_{b,j,c}^{(i)}}{\hat{p}_{b,j,c} ^{(i)}}\right] \tag{36}\] \[=\frac{1}{|\mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_ {b}}\left[\sum_{c\in\mathbb{L}_{b,j}^{(i)}}a_{b,j,c}\cdot\ln\left[\frac{a_{b, j,c}}{\hat{p}_{b,j,c}^{(i)}}\right]+\left[\sum_{k\in\mathbb{L}-\mathbb{L}_{b,j}^{(i )}}a_{b,j,k}\right]\cdot\ln\left[\frac{\sum_{k\in\mathbb{L}-\mathbb{L}_{b,j}^{( i)}}a_{b,j,k}}{\sum_{k\in\mathbb{L}-\mathbb{L}_{b,j}^{(i)}}\hat{p}_{b,j,k}^{(i)}} \right]\right].\]

As given by (10), the cross entropy loss for optimizing discriminative neural network classifiers was the variational free energy (VFE) of a retrospective active inference. That is,

\[\mathcal{L}_{\rm CE}(\mathbf{\hat{P}}_{b}^{(i)},\mathbf{L}_{b})=\frac{-1}{| \mathbb{L}|\cdot|\mathbb{T}_{b}|}\sum_{j\in\mathbb{T}_{b}}\sum_{e\in\mathbb{L }}l_{b,j,c}\cdot\ln\!\left(\hat{p}_{b,j,c}^{(i)}\right) \equiv -\sum_{s|\pi}p(s|\pi)\cdot\ln\!\left(q(o|\pi)\right). \tag{37}\]

Also, the expected free energy (EFE) of a prospective active inference was given in (4) as

\[\mathcal{L}_{\rm EFE}=\underbrace{\sum_{o}p(o)\cdot\left[\ln\!\left(p(o) \right)-\ln\!\left(q(o|\pi)\right)\right]}_{\rm expected\ complexity}+\underbrace{ \sum_{s|\pi}-p(s|\pi)\cdot\sum_{o|\pi}q(o|\pi)\cdot\ln\!\left(q(o|\pi)\right) }_{\rm uncertainty}. \tag{38}\]

Our proposed Algorithm 1 for finding the candidate labels \(\mathbb{L}_{b,j}^{(i)}\) aimed to minimize the objective function of the generalized Kelly criterion. This minimized function was given by (36). A comparison of (38) and (36) with regard to (37) revealed that the minimized objective of the Kelly criterion was the **expected complexity** term of the EFE of a prospective active inference. That is, the objective function of the generalized Kelly criterion was a tight upper bound of the expected complexity of the EFE. This equivalence got summarized in Table 2 and implied that the **preferred observations** denoted by \(o\) were realized through dividing \(\mathbb{L}\) into candidate \(\mathbb{L}_{b,j}^{(i)}\) and noncandidate classes \(\mathbb{L}-\mathbb{L}_{b,j}^{(i)}\) and then handling the noncandidate classes altogether as one class. To this end, in (36), the prior (win) probabilities of the noncandidate classes got summed together to form their collective prior (win) probability. Similarly, the estimated

\begin{table}
\begin{tabular}{|c|l|} \hline \(p(s|\pi)\) & \(l_{b,j,c}:\;c^{\rm th}\) entry of the vectorized reference label of the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) \\ \(p(o)\) & \(a_{b,j,c}:\) prior (win) probability of the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) \\ \(q(o|\pi)\) & \(\hat{p}_{b,j,c}^{(i)}:\) estimated posterior (belief) probability of the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) \\ \hline \hline \(\sum_{k\in\mathbb{L}-\mathbb{L}_{b,j}^{(i)}}a_{b,j,k}\): collective prior of noncandidate classes of the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) \\ \(\sum_{k\in\mathbb{L}-\mathbb{L}_{b,j}^{(i)}}b_{b,j,k}^{(i)}\): collective posterior of noncandidate classes of the sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\rm train}\) \\ \hline \end{tabular}
\end{table}
Table 2: Equivalence of the notations used in the objective functions of the active inference (left column) and the neural network optimization (right column).
posterior (belief) probabilities of the noncandidate classes got summed together to form their collective posterior (belief) probability.

The EFE in (38) was composed of an **expected complexity** term plus an **uncertainty** term. As described in subsection 1.1, the minimization of the expected complexity was equivalent to the maximization of the reward. The reward

Figure 8: Sagittal slices of the feature maps at the spatial regions enclosing the vertebral bodies and the intervertebral disks at the outputs of different encoder/decoder stages of the baseline architecture depicted in Figure 5 after being optimized by the proposed objective function and its associated optimization process.
maximization was also a goal of the Kelly criterion and could thus be partially fulfilled by finding the candidate labels through the proposed Algorithm 1. To further maximize the reward, the expected complexity should be minimized further. This was doable by having enough information or maximizing the information gain, i.e. minimizing the **uncertainty**. Accordingly, to optimize a discriminative neural network classifier, we proposed a novel objective function based on the EFE of a prospective active inference. The proposed function was given by

\[\mathcal{L}_{\mathrm{EFE}}(\hat{\mathbf{P}}_{b}^{(i)},\mathbf{A}_{b},\mathbf{ L}_{b})=\underbrace{\frac{-1}{\left\|\mathbb{L}\right|\cdot\left|\mathbb{T}_{b} \right|}}_{\begin{subarray}{c}j\in\mathbb{T}_{b}\end{subarray}}\sum_{c\in \mathbb{L}}\sum_{b,j,c}\mathbb{L}_{b,j,c}\cdot\hat{p}_{b,j,c}^{(i)}\cdot\ln \left[\hat{p}_{b,j,c}^{(i)}\right]\\ \mathrm{uncertainty}\]

This function was reversible and differentiable with respect to the posteriors \(\hat{\mathbf{P}}_{b}^{(i)}\). As given by (7), these posteriors were generated by applying the Softmax function to the network's outputs \(\mathbf{Z}_{b}^{(i)}=\left[\mathbf{z}_{b,j}^{(i)}\right]_{j}=\left[z_{b,j,c} ^{(i)}\right]_{j,c}\). Thus, the proposed function was also differentiable with respect to the \(\mathbf{Z}_{b}^{(i)}\) and the outputs of every layer. As described in subsection 1.2, these allowed to minimize it by a gradient descent optimizer with backpropagation.

We preceded the minimization of (3) with a partial minimization of its expected complexity term by finding the candidate classification labels \(\mathbb{L}_{b,j}^{(i)}\) of each sample (bettor) \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) through the Algorithm 1 proposed based on the Kelly criterion.

Accordingly, in each iteration \(i\in\{1,\cdots,n_{\mathrm{it}}\}\) of our proposed optimization process, every sample \(v_{j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) got passed through the network to estimate its classification posteriors \(\hat{\mathbf{P}}_{b}^{(i)}=\left[\hat{\mathbf{p}}_{b,j}^{(i)}\in(0,1)\right]_ {j}=\left[\hat{p}_{b,j,c}^{(i)}\right]_{j,c}\). From these posteriors and the fixed priors \(\mathbf{a}_{b,j}=\left[a_{b,j,c}\in(0,1)\right]_{c\in\mathbb{L}}\) of the sample, its candidate classification labels \(\mathbb{L}_{b,j}^{(i)}\subseteq\mathbb{L}\) got computed by using the proposed Algorithm 1. Then, the loss at the last network's layer got obtained by inputting the posteriors, the priors, and the candidate labels of the samples into the proposed function in (3). By propagating this loss from the last layer to the first layer, the loss of every layer got obtained. Then, the gradient (first derivative) of each layer's loss got calculated with respect to its outputs. The product of these layerwise gradients got used by the gradient descent optimizer to update the network's parameters.

In an image segmentation task, each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) was an image patch processed by a network's layer. In our baseline architecture described in subsection 1.4, each network's layer processed samples (patches) of a certain spatial resolution. The multiresolution hierarchy of the network was the result of downsampling and upsampling each volumetric fat-water image through convolutional and deconvolutional layers, respectively. For sake of simplicity, we omitted the resolution specifying indices from the samples' notations.

Figure 8 shows sagittal slices of the feature maps at the spatial regions enclosing the vertebral bodies and the intervertebral disks at the outputs of different encoder/decoder stages of the baseline architecture depicted in Figure 5 after being optimized by the proposed objective function and its associated optimization process.

## 5 Network's Parameters and Their Optimization

Our proposed algorithm for finding the candidate labels and our proposed objective function for optimizing a discriminative neural network classifier got integrated into a mini-batch-based gradient descent optimizer with backpropagation by using the process proposed in section 4. This process got evaluated against a similar process incorporating a representative of the cross entropy-based losses or a representative of the metric-based losses introduced in subsection 1.3. The representative of the cross entropy-based losses was the weighted focal loss. This loss comprised of a modulating factor and a weighting mechanism to alleviate classification biases towards the dominant classes of the training samples.
The representative of the metric-based losses was the Lovasz-Softmax loss. Besides being smooth and differentiable, to the best of our knowledge, this loss was the only convex loss among the metric-based losses.

Accordingly, the evaluated losses were

1. the proposed objective function (Po) given by (39)
2. the weighted focal loss (Fo) given by (14)
3. the Lovasz-Softmax loss (Lo) given by (21).

These evaluations were on an end-to-end optimization of the baseline architecture described in subsection 1.4. For each case, the baseline architecture was once used without attention gates (Na) as depicted in Figure 5 and once used with the attention gates (At) as depicted in Figure 7. Also, for (2) and (3) each training sample was accompanied by its reference (ground truth) label to fulfill the supervised nature of these objective functions. However, our proposed algorithm for finding the candidate labels and our proposed objective function got evaluated according to a fully supervised, a semi-supervised, and an unsupervised approach. These resulted in the training samples being

1. accompanied by their reference labels and their priors (GrPr) \(\rightarrow\) fully supervised
2. only accompanied by their reference labels (GrNp) \(\rightarrow\) semi-supervised
3. only accompanied by their priors (NgPr) \(\rightarrow\) semi-supervised
4. accompanied by neither their reference labels nor their priors (NgNp) \(\rightarrow\) unsupervised.

For the cases with the priors, the prior probabilities of the training samples could be computed by a multiatlas registration. If no prior probabilities were provided at the time of optimization (training), then uniform priors got assumed. If the reference (ground truth) labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were provided at the time of optimization (training), then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vectorized reference label \(\mathbf{l}_{b,j}\) was the one-hot-encoding of its reference label \(l_{b,j}\in\mathbb{L}\) and was given by (8). If the reference labels of the training samples \(\mathbb{T}_{\mathrm{train}}\) were not provided at the time of optimization, then for each sample \(v_{b,j}\in\mathbb{T}_{b}\subseteq\mathbb{T}_{\mathrm{train}}\) the vector \(\mathbf{l}_{b,j}\) was uniform and given by (9).

For each evaluation case, the main parameters and the hyperparameters of the baseline architecture got trained (optimized) to automatically segment \(n_{\mathrm{clas}}=|\mathbb{L}|=8\) classes of vertebral bodies (VBs), intervertebral disks (IVDs), psoas major (PM) and quadratus lumborum (QL) muscles, epicardial adipose tissues (EpAT), pericardial adipose tissues (PeAT), cardiac perivascular adipose tissues (PvAT), and background on each volumetric fat-water image. To this end, the volumetric fat-water images got divided into a training and a test set. The training set formed the samples set \(\mathbb{T}_{\mathrm{train}}\) and got used to optimize the main parameters and the hyperparameters of the baseline architecture by each method. The test set formed the samples set \(\mathbb{T}_{\mathrm{test}}\) and got used to evaluate the classification performance of the baseline architecture after being fully optimized by each method. The training set was composed of samples accompanied by their reference labels and priors. The test set was composed of samples accompanied by their reference labels. The reference labels of the test samples were not fed to the neural network. They were rather compared against the corresponding labels predicted by the network to evaluate the classification performance of the network. The predicted label of each sample was the index of its maximum classification posterior estimated by the network.

The main parameters of the baseline architecture included the weights and the biases of the convolutional and deconvolutional layers, the leakage coefficient \(a_{\mathrm{prelu}}\in\mathbb{R}_{\geq 0}\) of every nonlinear PReLU activation, and the means and variances of the (instance) normalizers introduced in page 14. Prior to the optimization of the main parameters, they should be initialized. This initialization was extremely important for the weights of the convolutional and deconvolutional layers of a residual network of several layers and thus different paths of signal propagation. Without a proper weight initialization, some parts of the network might have excessive activations and thus produce stronger gradients while some other parts might produce weaker gradients and thus get optimized less. To avoid this, a random initialization of the weights with the aim of breaking symmetries and making each feature map of a unit variance
was suggested. For this, the weights were drawn from a certain distribution. In networks with nonlinear Sigmoid or hyperbolic tangent activations as well as linear activations, the proper initializations of the weights of every layer were random numbers drawn from a uniform distribution in the range \([-\sqrt{6/(n_{\mathrm{in}}+n_{\mathrm{out}})},\;\sqrt{6/(n_{\mathrm{in}}+n_{ \mathrm{out}})}]\) with \(n_{\mathrm{in}}\) being the number of incoming network connections (fan-in) and \(n_{\mathrm{out}}\) being the number of outgoing network connections (fan-out) of the layer. This type of initialization was called a _Glorot_ or a _Xavier_ initialization and was shown to be improper for networks involving nonlinear rectified linear units, including the PReLU, as their activations [Glorot 2010]. For these networks, like our baseline architecture, the proper initializations of the weights of every convolutional/deconvolutional layer were random numbers drawn from a Gaussian distribution with a mean of 0 and a standard deviation of \(\sqrt{2/n_{\mathrm{in}}}\)[He 2015, Ronneberger 2015]. For a convolutional layer of a kernel size of \(5\times 5\times 5\), \(16\) input feature maps, and \(32\) output feature maps, the number of incoming network connections (fan-in) was \(5\times 5\times 5\times 16=2000\) and the number of outgoing network connections (fan-out) was \(32\). The biases of every convolutional/deconvolutional layer were initialized to 0. The leakage coefficient of every nonlinear PReLU activation got initialized to \(0.15\) to allow a small leakage of negative inputs. The means and the variances of the (instance) normalizers got initialized to 0 and 1 respectively.

The hyperparameters of the baseline architecture and their discretized values were

* number of convolutional/deconvolutional layers \(n_{s}\in\{1,2,\cdots,5\}\) of the \(s^{\mathrm{th}}\) encoder/decoder stage of the V-net of the baseline architecture
* Dropout's retention probability \(p_{s}\in\{0.1,0.2,\cdots,0.9\}\) of the perceptrons (nodes) of the \(s^{\mathrm{th}}\) encoder/decoder stage of the V-net of the baseline architecture.

To optimize the main parameters and the hyperparameters of the baseline architecture by each method, a random search over the discretized hyperparameter values and a 5-fold cross validation were conducted. To this end, the training set got divided into 5 subsets. Then, for each method, in each optimization trial, a set of hyperparameter values got randomly selected. With these hyperparameter values, 5 times training and validation got performed according to the 5-fold cross validation. In each fold, the main parameters of the baseline architecture got optimized on 4 subsets by using a mini-batch-based gradient descent optimizer with backpropagation. The gradient descent optimizer was the Adam optimizer described in subsection 1.2. The resulting network model got then evaluated on the remaining (validation) subset by calculating the _precision_ and the _recall_ metrics for each of the \(n_{\mathrm{clas}}-1=8-1=7\) foreground classes against the rest of the classes. This way, for the selected hyperparameter values, at the end of the 5-fold cross validation, 5 network models and 7 _precision_ and 7 _recall_ values per network model were obtained. For each model, the 7 _precision_ and the 7 _recall_ values got averaged. Then, for the selected hyperparameter values, the model of maximum averaged _precision_ and _recall_ was the best performing model. The optimization trials continued by randomly selecting another set of hyperparameter values until the best performing model resulted from the current hyperparameter values could not exceed the averaged _precision_ and _recall_ values of any of the best models in the last 50 trials. The _precision_ and _recall_ metrics were selected due to their robustness against the imbalanced class-sample distributions. Moreover, the aforementioned cross validation aimed to reduce the impacts of the randomized initialization of the main parameters on the resulting network models. The 5 folds were selected with regard to the maximum size of the baseline architecture and the sufficiency of the number of training and validation samples for the optimization and evaluation in each fold, respectively. The above process was done by using the tools provided in the distributed asynchronous hyperparameter optimization (Hyperopt) library in Python [Bergstra 2015]. For the hyperparameter selection, in addition to the randomization, this library provided a tree of Parzen estimators (TPE) and its adaptive variant. The TPE was more appropriate for belief neural networks of undirected graph topology than the feed-forward networks like our baseline architecture [Bergstra 2011, Bergstra 2012].

The evaluated objective functions and the Adam-based gradient descent optimizer involved following fixed parameters:

* N=\(|\mathbb{T}_{b}|=2\): As explained in page 13, due to the memory limitations of the used GPU, only 2 volumetric fat-water images were included in each mini-batch.
Figure 9: Convergence patterns of different evaluation cases with each case optimizing its main parameters with the best performing hyperparameters.
* \(\gamma_{\rm mod}=2\): Modulating factor of the focal loss given by (13).
* \(\alpha_{\rm lr}=0.001\): Learning rate (step size) of the gradient descent optimizer defined in (5). This learning rate did not need to be adapted manually as the Adam optimizer automatically changed the effective learning rate by the ratio of the exponential moving average of the first moment to the exponential moving average of the second moment.
* \(\beta_{\rm fm}=0.90\): Decay rate of the estimated first moments.
* \(\beta_{\rm sm}=0.99\): Decay rate of the estimated second moments.
* \(\mathbf{m}^{(0)}=\mathbf{0}\): Initial first moments.
* \(\mathbf{v}^{(0)}=\mathbf{0}\): Initial second moments.

The number of iterations \(n_{\rm it}\in\{10,\cdots,15000\}\) was determined according to an early stopping criterion. That is, when the exponential moving average of the validation error (loss) was not improved within the last 100 iterations, then the optimization got stopped.

Figure 9 shows convergence patterns of different evaluation cases with each case optimizing its main parameters with the best performing hyperparameters.

The aforementioned optimizations were conducted on 4 NVIDIA TITAN X(r) GPUs of 12 GB memory each and by using a memory efficient cuDNN3 implementation of the convolutional/deconvolutional layers and the TensorFlow(tm) library of version 2.3 [1].

Table 3 shows the optimized hyperparameters and the overall time of optimizing the main parameters and the hyperparameters for each evaluation case. After the optimizations, an automatic segmentation of the \(n_{\rm clas}=8\) classes on an unseen volumetric fat-water image took around 3 seconds for each evaluation case on the GPUs used for the optimizations.

\begin{table}
\begin{tabular}{|c|c|c c c c c c|} \hline  & & \multicolumn{6}{c|}{**Evaluation Case**} \\ \cline{2-9}  & & PoNaGrPr & PoNaGrNp & PoNaNpGr & PoNaNpNp & PoAtGrPr & PoAtGrNp \\ \hline \multirow{6}{*}{**Fashion**} & \(n_{1}\) & 2 & 2 & 2 & 3 & 1 & 1 \\  & \(n_{2}\) & 2 & 3 & 2 & 4 & 2 & 2 \\  & \(n_{3}\) & 3 & 3 & 4 & 4 & 2 & 2 \\  & \(n_{4}\) & 3 & 4 & 4 & 5 & 3 & 3 \\  & \(n_{5}\) & 3 & 4 & 5 & 5 & 3 & 4 \\  & \(p_{1}\) & 0.9 & 0.8 & 0.8 & 0.7 & 0.9 & 0.9 \\  & \(p_{2}\) & 0.8 & 0.7 & 0.7 & 0.7 & 0.7 & 0.8 \\  & \(p_{3}\) & 0.7 & 0.7 & 0.7 & 0.7 & 0.7 & 0.7 \\  & \(p_{4}\) & 0.7 & 0.7 & 0.7 & 0.6 & 0.6 & 0.7 \\  & \(p_{5}\) & 0.6 & 0.6 & 0.6 & 0.6 & 0.5 & 0.6 \\  & time & 83 & 90 & 95 & 107 & 82 & 85 \\ \hline \hline \multirow{6}{*}{**Fashion**} & \multicolumn{6}{c|}{**Evaluation Case**} \\ \cline{2-9}  & PoAtNgPr & PoAtNgNp & PoNaGrNp & PoAtGrNp & LoNaGrNp & LoAtGrNp \\ \cline{2-9}  & \(n_{1}\) & 1 & 2 & 2 & 2 & 1 & 1 \\  & \(n_{2}\) & 2 & 3 & 3 & 2 & 2 & 2 \\  & \(n_{3}\) & 3 & 4 & 4 & 3 & 3 & 2 \\  & \(n_{4}\) & 3 & 4 & 4 & 3 & 4 & 3 \\  & \(n_{5}\) & 4 & 5 & 4 & 4 & 4 & 3 \\  & \(p_{1}\) & 0.9 & 0.9 & 0.9 & 0.9 & 0.9 & 0.9 \\  & \(p_{2}\) & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 \\  & \(p_{3}\) & 0.7 & 0.8 & 0.8 & 0.7 & 0.7 & 0.7 \\  & \(p_{4}\) & 0.6 & 0.7 & 0.7 & 0.6 & 0.7 & 0.6 \\  & \(p_{5}\) & 0.5 & 0.6 & 0.6 & 0.5 & 0.5 & 0.5 \\  & time & 88 & 97 & 84 & 82 & 79 & 77 \\ \hline \end{tabular}
\end{table}
Table 3: Optimized hyperparameters and the overall time [hours] of optimizing the main parameters and the hyperparameters for each evaluation case.