# Learning Discrete Weights and Activations Using the Local Reparameterization
Trick

Guy Berger
Bar-Ilan University, Israel
<EMAIL_ADDRESS>
Aviv Navon
Bar-Ilan University, Israel
<EMAIL_ADDRESS>
Ethan Fetaya
Bar-Ilan University, Israel
<EMAIL_ADDRESS>

###### Abstract

In computer vision and machine learning, a crucial challenge is to lower the
computation and memory demands for neural network inference. A commonplace
solution to address this challenge is through the use of binarization. By
binarizing the network weights and activations, one can significantly reduce
computational complexity by substituting the computationally expensive
floating operations with faster bitwise operations. This leads to a more
efficient neural network inference that can be deployed on low-resource
devices. In this work, we extend previous approaches that trained networks
with discrete weights using the local reparameterization trick to also allow
for discrete activations. The original approach optimized a distribution over
the discrete weights and uses the central limit theorem to approximate the
pre-activation with a continuous Gaussian distribution. Here we show that the
probabilistic modeling can also allow effective training of networks with
discrete activation as well. This further reduces runtime and memory footprint
at inference time with state-of-the-art results for networks with binary
activations.

##  1 Introduction

As neural networks become more powerful and their applications more
commonplace, there is an important and growing need to reduce computational
costs and memory requirements. This is especially important in edge devices,
e.g., smartphones, that have weaker processors and need to optimize energy
consumption. To address this challenge, one promising approach is to train
binary or ternary networks, where the weights are constrained to take on only
a small number of discrete values [3, 9, 33, 8, 23]. One can also binarize the
activations (e.g., using the sign function) to further improve the efficiency
of the network, albeit with a larger reduction in accuracy [26, 17, 32, 24,
25, 10, 6, 28, 29, 11].

Most previous works that have utilized sign activation directly apply it
during the forward pass and use heuristic or approximation methods to estimate
the gradients during the backward pass, e.g., straight-through gradient or
other biased gradient estimators [17, 3, 26, 28, 4]. In this paper, we propose
a novel method to compute gradients for networks with discrete activation by
employing a smooth approximation. We construct a fully differentiable
probabilistic model to approximate the discrete network during training. After
training, we sample from our trained probabilistic model to get our discrete
weights.

Our method is based on the observation by Kingma et al. [21] which states that
if the weights of a certain layer are sampled from independent Gaussian
distributions, then one can get better stochastic gradient estimation by
modeling and sampling the Gaussian distribution of the pre-activations instead
of weights. They called it the local reparameterization trick. O. Shayer and
Fetaya [27] extended the local reparameterization trick to train a network
with discrete weights. When the weights are discrete and stochastic, the pre-
activations can still be well approximated by a Gaussian distribution
according to the (Lyapunov) central limit theorem (CLT). Then, by using the
reparameterization trick, we can compute the derivatives of a smooth
distribution. In this work, we extend [27] and train a network with discrete
weights and activations. Based on the observation that the pre-activations are
well approximated by a Gaussian distribution, we propose sampling discrete
values from the induced distribution given our Gaussian approximation. We
construct a probabilistic model which is fully differentiable, using the
Gaussian CDF function to calculate the activation probabilities and Gumbel-
Softmax [1, 7] as a smooth approximation for categorical sampling.

We demonstrate the effectiveness of our approach on several benchmark
datasets. The experiments show that our proposed method, which we name LAR-
nets (Local Activation Reparameterization networks), outperforms previous SoTA
binarization approaches. To summarize, we make the following contributions:

  * •

Introduce a novel approach for learning with discrete activations based on the
local reparameterization trick.

  * •

Propose a new variant for batch normalization to extend the applicability of
the normalization layer to distribution over weights.

  * •

Demonstrate the effectiveness of the proposed approach against multiple
baselines for training discrete networks on several vision benchmarks.

(a)

(b)

(c)

Figure 1:  (a) The structure of our LAR-Net layer with discrete weights and
activation. The reparameterization allows the gradients to flow to $\theta$ in
the backward pass through all the differentiable nodes; (b) The structure of a
basic Conv-BN-Act block during the training phase; (c) The structure of basic
Conv-BN-Act layer on evaluation phase (after sampling weights from $\theta$);

##  2 Related Work

Binary networks have been an active area of research in deep learning over the
past few years. Courbariaux et al. [3] introduced BinaryConnect, a method for
training binary neural networks using binary stochastic weights. They sample
binary weights and compute gradients as if it was a deterministic full-
precision network. BNN [17] and XNOR-Net [26] suggested a different approach;
They discretize the weights and activations during the forward pass and back-
propagate through this non-continuous discretization using the _straight-
through_ estimator. XNOR-Net also introduces a continuous scaling factor to
the binarized weights, which became popular in later works. Other similar
works based on the straight-through estimator extend it to ternary networks
(TWN [9], TTQ [33]).

Motivated by the high potential to reduce model complexity and improve
computational efficiency by also discretizing activations, works that
discretize both weights and activations became more common. DoReFa-Net [32]
extended XNOR-Net to accelerate the training process using quantized
gradients. ABC-Net [24] suggested reducing the quantization error by linearly
combining multiple binary weight matrices and scaling factors to fit the full-
precision weights and activations. Other works suggested making modifications
in the network architectures, such as adding tailored-made layers, to make
them more fit to binary neural network training [29, 25]. Some recent works
tried to achieve good quantization functions in forward propagation, which can
reduce the gradient error as well. Gong et al. [11] presented a differential
soft quantization (DSQ) that has more accurate gradients in the backward pass.
Diffenderfer and Kailkhura [5] proposed what they called Multi Prize Ticket
(MPT). They incorporated a binarization scheme with weight pruning (which
results in a ternary network). For the weights and activations binarization
they suggested a modified gradient estimator for the sign function. For the
weights pruning, inspired by Frankle and Carbin [10], they proposed a scheme
to prune a randomly initialized network using a learnable mask by updating
pruning scores during the training. Ding et al. [6], suggested reshaping the
data distribution before binarization by adding a distribution loss for
learning the proper binarization (BNN-DL). Qin et al. [28] suggested IR-Net,
which aimed to minimize the information loss by maximizing the information
entropy of the quantized parameters and minimizing the quantization error.
Several other works discretize the activation, but to more than one bit [8,
16, 23].

##  3 Background

In this section, we describe the local reparametrization trick and the LR-net
approach.

#### Notation.

Let $\\{(x_{i},y_{i})\\}_{i=1,...,N}$ denote our training examples. We denote
by $W$ the set of all model parameters and let $W^{(t)}$ denote the weights
matrix for layer $t$. Furthermore, let $w_{ij}^{(t)}$ denote the element of
$W^{(t)}$, i.e., $[W^{(t)}]_{ij}=w^{(t)}_{ij}$. The pre-activations are
defined as $z^{(t+1)}=W^{(t)}h^{(t)}$, where $h^{(t)}$ is the activations of
layer $t$ (we omit the bias term for simplicity). We let $h^{(0)}=x$ and
define $h^{(t+1)}=\phi(z^{(t+1)})$, where $\phi(\cdot)$ is a non-linear
activation function.

We assume a stochastic model in which each weight $w_{ij}^{(t)}$ is sampled
independently from a multinomial distribution $\mathcal{W}_{ij}^{(t)}$. The
objective is minimizing the expected loss w.r.t. $\>\mathcal{W}$, where
$\mathcal{W}$ denotes the the distribution over $W$,

$L(\mathcal{W})=\operatorname{\mathbb{E}}_{W\sim\mathcal{W}}[\sum_{i=1}^{N}\ell(f(x_{i},W),y_{i})].$
(1)

The classic approach for minimizing Eq. 1 with a discrete distribution is the
log-derivative trick [31]:

$\nabla
L(\mathcal{W})=\mathbb{E}_{W\sim\mathcal{W}}\left[\sum_{i=1}^{N}\ell\left(f(x_{i},W),y_{i}\right)\nabla\log(P(W))\right].$
(2)

While this allows us to get an unbiased estimation of the gradient, it suffers
from high variance, which makes optimization using this method challenging.

A popular alternative for _continuous_ distributions is the
_reparameterization trick_ [20] \- instead of optimizing $\mathbb{E}_{x\sim
p_{\varphi}}[f(x)]$ for distribution parameters $\varphi$ we parametrize
$x=g(\epsilon;\theta)$ where $\epsilon$ is drawn from a known fixed
distribution $p(\epsilon)$ (usually Gaussian) and optimize
$\mathbb{E}_{p(\epsilon)}[f(g(\epsilon,\theta))]$ for $\theta$. To estimate
the gradient w.r.t $\theta$, we sample $\epsilon_{1},...,\epsilon_{m}$ and use
the Monte-Carlo approximation:

$\nabla_{\theta}\mathbb{E}_{p(\epsilon)}[f(g(\epsilon,\theta))]\approx\sum_{i=1}^{m}\nabla_{\theta}f(g(\epsilon_{i},\theta)).$
(3)

Kingma et al. [21] proposed an alternative for the task of variational
approximation for Bayesian neural networks. The authors make the key
observation that if a weight matrix $W^{(t)}$ is sampled from an independent
Gaussians $w_{ij}\sim\mathcal{N}(\mu_{ij},\sigma^{2}_{ij})$, then the pre-
activations $z^{(t+1)}=W^{(t)}h^{(t)}$ are distributed according to

$z_{i}^{(t+1)}\sim\mathcal{N}\left(\sum_{j}\mu^{(t)}_{ij}h^{(t)}_{j},\sum_{j}{\sigma_{ij}^{(t)}}^{2}{h_{j}^{(t)}}^{2}\right)$
(4)

This allows sampling the pre-activations instead of the model weights, which
results in lower-variance gradient estimations [21] and better optimization.
This approach is termed the _local reparameterization trick_.

The local reparameterization trick has been further extended to training
networks with _discrete_ weights. O. Shayer and Fetaya [27] introduced a
method to train networks with binary $\\{\pm 1\\}$ or ternary $\\{-1,0,1\\}$
weights (but it can be applied to a wider range of discrete values). The key
idea behind this method is that while the pre-activations
$z_{i}=\sum_{j}w_{ij}h_{j}$ 111We omit the layer index $(t)$ from that point
to simplify the notation. are discrete, from the (Lyapunov) central limit
theorem (CLT) they are well approximated by the Gaussian distribution
$z_{i}\sim\mathcal{N}(\sum_{j}\mu_{ij}h_{j},\,\sum_{j}\sigma_{ij}^{2}h_{j}^{2})$.
We note that $\mu_{ij}$ and $\sigma_{ij}$ are now the mean and variance of a
multinomial distribution, and not the mean-field Gaussian distribution. By
sampling $\epsilon_{i}\sim\mathcal{N}(0,1)$ we can represent the output as
$z_{i}=m_{i}+\epsilon_{i}\cdot v_{i}$, where
$m_{i}=\operatorname{\mathbb{E}}(z_{i})=\sum_{j}\mu_{ij}h_{j}$ and
$v_{i}^{2}=Var(z_{i})=\sum_{j}\sigma_{ij}^{2}h_{i}^{2}$. Let $\theta_{ij}$ be
the parameters of the multinomial distribution over $w_{ij}$, then, our goal
is to minimize the expected loss $\ell$ $w.r.t.$ $\theta$ (where $\theta$
denotes the set of all parameters $\theta_{ij}$),

$\nabla_{\theta}L(\mathcal{\theta})=\operatorname{\mathbb{E}}_{p(\epsilon)}[\sum_{i=1}^{N}\ell(f(h_{i},\epsilon,\theta),y_{i})].$
(5)

##  4 Our Method

In this section, we describe our LAR-net approach. We propose a novel
extension to the LR-net [27] that allows for learning discrete activations.
Furthermore, we describe a novel batch-normalization layer to complement our
method. Finally, we outline the inference procedure using our method.

###  4.1 Learning Discrete Activations

Here we propose a novel extension to [27] for training networks with discrete
weights and activations. Assume we have a network with a distribution over
discrete weights and the $sign$ function as its non-linearity. If we
approximate the pre-activation output $z_{i}$ by the Gaussian distribution, we
can directly calculate the distribution of the values of the discrete
activation $sign(z_{i})$ using the probability that the sampled Gaussian would
be positive

$sign(z_{i})=\begin{cases}+1,&\text{$with\;probabailty\>\>p_{i}=1-\Phi(-\frac{m_{i}}{v_{i}})$}.\\\
-1,&\text{$with\;probabailty\>\>1-p_{i}=\Phi(-\frac{m_{i}}{v_{i}})$}.\end{cases}$
(6)

where $\Phi$ is the Gaussian CDF function. We note that the distributions of
$sign(z_{i})$ are independent, as each one depends on a different row in the
matrix $W$. We also note that if we tried to compute the distribution on
subsequent layers, the outputs would not be independent anymore due to the
shared inputs. As this would impair our CLT approximation, we instead choose
to sample the discrete activation at each layer.

In order to differentiate through the discrete activation sampling, we use the
Gumbel-Softmax approximation [1, 7]. This method builds on the Gumble-Max
trick [12], presented in eq. (7), that showed one can get a categorical sample
by taking the maximum of the logits with additional Gumbel distributed noise,

$\hat{h}=one\\_hot(\arg\max_{k}[g_{k}+\log\pi_{k}]),\>\>\>\>g_{k}\sim
Gumbel(0,1)$ (7)

The approximation in [1, 7] replaces the discrete max with a smooth Softmax
function with a temprature parameter $\tau$. It can be shown that as
$\tau\rightarrow 0$ the Softmax converges to the max and we get a sample from
the desired categorical distribution.

$\hat{h}_{k}=\frac{exp((log(\pi_{k})+g_{k})/\tau)}{\sum^{K}_{j=1}exp((log(\pi_{j})+g_{j})/\tau)}\>\>\>\>for\>\>k=1...K$
(8)

We are interested in discrete variables with two classes ($\\{\pm 1\\}$). We
define the class probabilities vector as $\pi_{i}=[1-p_{i},p_{i}]$. Using eq.
(8) we can generate a samples vector $[\hat{h}_{i,-1},\hat{h}_{i,1}]$ (where
$\hat{h}_{i,-1}+\hat{h}_{i,1}=1$) and then multiply it with our discrete
values $[-1,1]$ (each discrete value is multiplied with the corresponding
element in $\hat{h}_{i}$), resulting in a very simple formula:

$h^{(t+1)}_{i}=\hat{h}^{(t+1)}_{i,1}-\hat{h}^{(t+1)}_{i,-1}\approx
sign(z_{i}^{(t+1)})$ (9)

When $\tau\rightarrow 0$ and the Softmax is replaced with the argmax, we get
the exact sign activation. We note that the Gumbel-Softmax has two variants:
“hard" and “soft", our reported results are with the “hard" variant, which
achieved slightly better results (see [1] for details).

This leads to a simple algorithm for the forward pass in the training phase.
Let $\theta_{ij}$ be the parameters of the multinomial distribution over
$w_{ij}$. At the forward pass we compute the weights mean
$\mu_{ij}=\operatorname{\mathbb{E}}_{\theta_{ij}}[w_{ij}]$ and variance
$\sigma_{ij}^{2}=Var_{\theta_{ij}}[w_{ij}]$. Then, we calculate the mean and
variance of $z_{i}$, $m_{i}=\sum_{j}\mu_{ij}h_{j}$ and
$v_{i}=\sum_{j}\sigma_{ij}^{2}h_{j}^{2}$. Using the mean and variance, we
calculate the discrete distribution of $sign(z_{i})$, and then we use the
Gumbel Softmax to (approximately) sample from that distribution. The forward
pass during the training pahse is summarized in Algorithm 1. The backward pass
is straightforward, since all the functions are differentiable (see Fig.
1(a)).

Algorithm 1 Discrete Weights and activation layer

Input: $h\in\mathbb{R}^{d_{in}}$
Parameters: Multinomial parameters $\theta_{ij}$ for each weight
Forward Pass:
1: Compute $\mu_{ij}=\mathbb{E}_{\theta_{ij}}[w_{ij}]$ and
$\sigma_{ij}^{2}=Var_{\theta_{ij}}[w_{ij}]$
2: Compute $m_{i}=\sum_{j}\mu_{ij}h_{j}$ and
$v_{i}=\sum_{j}\sigma_{ij}^{2}h_{i}^{2}$
3: Compute
$\pi_{i}:\>\>\>p_{i}=P(sign(z_{i})=1)=\frac{1}{2}(1+erf(-\frac{m_{i}}{\sqrt{2}v_{i}}))$
and $1-p_{i}$
4: Sample $g_{i}\sim Gumbel(0,1)$
5: Use Gumbel-Softmax to generate $\hat{h}_{i}$ (eq. (8))
Return: $h\in\mathbb{R}^{d_{out}}$ with $h_{i}=\hat{h}_{i,1}-\hat{h}_{i,-1}$

###  4.2 Batch Normalization

Batch normalization [18] is a common part of convolutional neural networks
that is known to accelerate training and improve performance in many cases.
Specifically for discrete networks, the authors in [5] showed a high gain just
by including batch normalization layers. In our model, before the sign
activation function, we have the pre-activation distributions instead of
deterministic values (see 1(b)), and we cannot perform the classic batch
normalization layer.

To address this issue, we propose a new batch normalization layer for our
training process, which can be applied to distributions. The pre-activation
$z_{bcij}$ ($b$ is the batch index and $c$ is the channel index) of each
convolution are approximated by Gaussian random variables and are distributed
according to $z_{bcij}\sim\mathcal{N}(\mu_{bcij},\sigma_{bcij}^{2})$. The mean
of the normalized variable can be easily calculated as,

$\mu_{c}=E_{bij}[Z_{bcij}]=\frac{1}{B\cdot H\cdot
W}\sum\limits_{b,i,j}\mu_{bcij}.$ (10)

Where $B,W,H$ are the batch size, the width, and the height of the image
respectively. Using the Law of total variance,
$Var(Y)=Var[E[Y|X]]+E[Var[Y|X]]$, we can also calculate the variance of the
normalized variable:

$\sigma_{c}^{2}=\frac{1}{B\cdot H\cdot
W}\sum\limits_{b,i,j}(\mu_{bcij}-\mu{{}_{c}})^{2}+\frac{1}{B\cdot H\cdot
W}\sum\limits_{b,i,j}\sigma_{bcij}^{2}$ (11)

Batch normalization usually also includes learnable affine parameters $\gamma$
and $\beta$, which we incorporate here as well. The pre-activation after batch
normalization layer $z_{BN}$ is distributed according to

$z_{BN_{bcij}}\sim\mathcal{N}\left(\gamma_{c}\cdot\frac{\mu_{bcij}-\mu_{c}}{\sigma_{c}}+\beta_{c},\frac{\gamma^{2}_{c}\sigma_{bcij}^{2}}{\sigma_{c}^{2}}\right).$
(12)

###  4.3 Inference

At inference or test time, we wish to apply our model with the actual discrete
weights and activations and not the approximation used for training. To do
that, we sample discrete weights based on the probabilities we obtained during
training, and perform a standard inference using the sampled weights. This
process can be repeated several times to identify the best-performing weights.
However, it is worth noting that our trained distribution typically has low
entropy, resulting in little variation in test accuracy between weight
samples. As a result, only slight improvements in overall performance are
observed.

##  5 Implementation Details

###  5.1 Network Architecture

We use the ResNet-18 [14] and VGG [30], as is standard for experiments on
discrete neural networks. In all our experiments, similar to prior works, the
first layer and last fully-connected layer remain in full precision.

###  5.2 Optimization Details

We observed that incorporating multiple Monte-Carlo samples for each input
element significantly enhances the overall results. In every batch, we run our
model on each datum multiple times, each with a different random sample
resulting in the following gradient estimation

$\nabla_{\theta}L(\mathcal{\theta})\approx\sum_{i=1}^{N}\sum_{j=1}^{S}\nabla\ell(f(h_{i},\epsilon_{j},\theta),y_{i}).$
(13)

#### Weights Distributions Entropy.

Similar to [27], we added $probability\>\>decay$ regularization. We used it as
$L_{2}$ regularization on the distribution parameters. It helps to increase
the weights entropy and prevent many weights distributions converged to a
deterministic value, which degrades the CLT approximation and leads to
suboptimal results.

#### Activations Distributions Entropy.

In some cases, we observed that there is a need to reduce the entropy of the
activations distributions of the last binary layer. We found that many of the
activations probabilities $p_{i}$ in that layer converge to values around
$0.5$, which means an undesirable extremely high amount of randomness. We
found that giving a lower learning rate to the fully connected layer (after
the last binary layer) helps to reduce the entropy. We show this phenomenon in
Fig. 2.

#### Initialization.

Similar to [27] we use a pre-trained network to initialize the weight
distribution parameters $\theta$. We found that using a network with _Tanh_ as
an activation function instead of _ReLU_ results in a better initialization
for $\theta$, as it is closer to our discrete sign activation. For further
details on the initialization scheme, please refer to the supplementary
material.

(a)

(b)

(c)

Figure 2:  (a) Activation entropy example from an intermediate layer of the
network. The entropy is relatively low, meaning low randomness, but the model
is not deterministic; (b) Activation entropy example from the last binary
layer, without learning rate decay; (c) Example of activation entropy from the
last binary layer after reducing the learning rate of the last FC layer.

##  6 Experiments

In this section, we detail the conducted experiments. We use two benchmark
datasets, CIFAR-10 and CIFAR-100 [22], to verify the effectiveness of our
proposed method and compare it with other state-of-the-art methods. Our code
is available in the supplementary material.

#### CIFAR-10.

We compare our results with prior works on CIFAR-10, including IR-Net[28], MPT
[5], BNN-DL[6], DSQ[11], BNN[17] and XNOR-Net[26] using VGG-small and
ResNet-18. We compared networks with 1-bit weights and 1-bit activations (1/1)
as well as networks with 2-bits weights and 1-bit activations (2/1). The extra
bit resulting in weights sparsity, which brings additional efficiency. The
performance comparison using the different methods is shown in Table 1 (we
presented the best performance for each method over VGG-small and ResNet-18).
The results show that our approach achieves better results than other existing
methods. The entire experimental details are provided in the supplementary
material.

#### CIFAR-100.

We also compared our model performance on CIFAR-100 with IR-Net[28], MPT[5],
and BNN-DL[6]. The results are shown in Table 2. Our approach achieves better
results than other existing methods, with the exception of BNN-DL which on
par. However, we note that the BNN-DL uses a slightly larger network. All the
experimental details regarding our training process and the evaluation for MPT
and IR-Net are included in the supplementary material.

Method | Topology | Bit-Width (W/A) | Top-1 Acc
---|---|---|---
Full-Precision | ResNet-18 | 32/32 | 95.01
LR-net | ResNet-18 | 2/32 | 94.88
IR-Net | ResNet-18 | 1/1 | 91.5
BNN-DL | ResNet-18 | 1/1 | 90.5
DSQ | VGG-small | 1/1 | 91.72
MPT | VGG-small | 2/1 | 91.9
BNN | VGG-small | 1/1 | 89.9
XNOR-Net | VGG-small | 1/1 | 89.8
LAR-net | ResNet-18 | 2/1 | 92.34 $\pm$ 0.05
Table 1: Performance comparison with Binary SOTA methods on CIFAR-10. Method | Topology | Bit-Width (W/A) | Top-1 Acc
---|---|---|---
Full-Precision | ResNet-18 | 32/32 | 75.61
LR-net | ResNet-18 | 2/32 | 73.13
IR-Net | ResNet-18 | 1/1 | 67.42 $\pm$ 0.06
IR-Net | VGG-small | 1/1 | 67.8 $\pm$ 0.04
BNN-DL | ResNet-20†††They used their variant of ResNet with 20 layers, we refer to this for convenience as resnet-20 | 1/1 | 68.17
MPT | VGG-small | 2/1 | 63.71 $\pm$ 0.16
LAR-net | ResNet-18 | 2/1 | 68.35 $\pm$ 0.02
Table 2: Performance comparison with Binary SOTA methods on CIFAR-100.

##  7 Conclusions

In this paper, we introduced a novel scheme for training a network with
discrete weights and activations, for efficient inference of vision and
learning applications. We demonstrate that using a probabilistic approach is
not limited to networks with discrete weights, but can also be successfully
applied to networks with discrete activations. Furthermore, we show how
standard DNN layers, such as batch normalization can also be applied in this
scenario. Finally, we evaluate our method on various image classification
datasets obtaining state-of-the-art results.

## References

  * C. J. Maddison and Teh [2017] Mnih C. J. Maddison, A and Y. Whye Teh.  The concrete distribution: A continuous relaxation of discrete random variables.  _International Conference on Learning Representations (ICLR)_ , 2017\.
  * chrundle [2022] chrundle.  biprop.  _https://github.com/chrundle/biprop_ , 2022.
  * Courbariaux et al. [2015] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David.  Binaryconnect: Training deep neural networks with binary weights during propagations.  _Advances in Neural Information Processing Systems (NIPS)_ , 2015\.
  * Darabi et al. [2018] Sajad Darabi, Mouloud Belbahri, Matthieu Courbariaux, and Vahid Partovi Nia.  BNN+: improved binary network training.  _CoRR_ , abs/1812.11800, 2018.  URL http://arxiv.org/abs/1812.11800.
  * Diffenderfer and Kailkhura [2021] J. Diffenderfer and B. Kailkhura.  Multi-prize lottery ticket hypothesis: Finding accurate binary neural networks by pruning A randomly weighted network.  _International Conference on Learning Representations (ICLR)_ , 2021\.
  * Ding et al. [2019] Ruizhou Ding, Ting-Wu Chin, Zeye Liu, and Diana Marculescu.  Regularizing activation distribution for training binarized deep networks.  _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * E. Jang and Poole [2017] S. Gu E. Jang and B. Poole.  Categorical reparameterization with gumbel-softmax.  _International Conference on Learning Representations (ICLR)_ , 2017\.
  * Faraone et al. [2018] Julian Faraone, Nicholas J. Fraser, Michaela Blott, and Philip H. W. Leong.  SYQ: learning symmetric quantization for efficient deep neural networks.  _2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018_ , pages 4300–4309, 2018.
  * Fengfu Li and Liu [2016] Bo Zhang Fengfu Li and Bin Liu.  Ternary weight networks.  _Advances in Neural Information Processing Systems (NIPS Workshop)_ , 2016.
  * Frankle and Carbin [2019] Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Training pruned neural networks.  _International Conference on Learning Representations (ICLR)_ , 2019\.
  * Gong et al. [2019] Ruihao Gong, Xianglong Liu, Shenghu Jiang, Tianxiang Li, Peng Hu, Jiazhen Lin, Fengwei Yu, and Junjie Yan.  Differentiable soft quantization: Bridging full-precision and low-bit neural networks.  _International Conference on Computer Vision (ICCV)_ , 2019.
  * Gumbel [1954] E. J. Gumbel.  Statistical theory of extreme values and some practical applications: a series of lectures.  _Number 33. US Govt. Print. Office_ , 1954.
  * He et al. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  _International Conference on Computer Vision (ICCV)_ , 2015.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016\.
  * htqin [2020] htqin.  Ir-net.  _https://github.com/htqin/IR-Net_ , 2020.
  * Hu et al. [2018] Qinghao Hu, Peisong Wang, and Jian Cheng.  From hashing to cnns: Training binary weight networks via hashing.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Hubara et al. [2016] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks.  _Advances in Neural Information Processing Systems (NIPS)_ , 2016\.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _International Conference on Machine Learning (ICML)_ , 2015.
  * Kingma and Ba. [2015] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations (ICLR)_ , 2015\.
  * Kingma and Welling [2014] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  _International Conference on Learning Representations (ICLR)_ , 2014\.
  * Kingma et al. [2015] Diederik P. Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick.  _Advances in Neural Information Processing Systems (NIPS)_ , 2015\.
  * Krizhevsky [2009] A. Krizhevsky.  Learning multiple layers of features from tiny images.  _Tech Report_ , 2009.
  * Lin et al. [2016] Darryl Dexu Lin, Sachin S. Talathi, and V. Sreekanth Annapureddy.  Fixed point quantization of deep convolutional networks.  _Proceedings of the 33nd International Conference on Machine Learning, ICML 2016, New York City, NY, USA, June 19-24, 2016_ , 48:2849–2858, 2016.
  * Lin et al. [2017] Xiaofan Lin, Cong Zhao, and Wei Pan.  Towards accurate binary convolutional neural network.  _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 345–353, 2017.
  * Liu et al. [2018] Zechun Liu, Baoyuan Wu, Wenhan Luo, Xin Yang, Wei Liu, and Kwang-Ting Cheng.  Bi-real net: Enhancing the performance of 1-bit cnns with improved representational capability and advanced training algorithm.  _European conference on computer vision (ECCV)_ , 2018.
  * Mohammad Rastegari and Farhadi [2016] Joseph Redmon Mohammad Rastegari, Vicente Ordonez and Ali Farhadi.  Xnor-net: Imagenet classification using binary convolutional neural networks.  _ECCV_ , 2016.
  * O. Shayer and Fetaya [2018] D. Levi O. Shayer and E. Fetaya.  Learning discrete weights using the local reparameterization trick.  _International Conference on Learning Representations (ICLR)_ , 2018\.
  * Qin et al. [2020] Haotong Qin, Ruihao Gong, Xianglong Liu, Mingzhu Shen, Ziran Wei, Fengwei Yu, and Jingkuan Song.  Forward and backward information retention for accurate binary neural networks.  _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_ , 2020.
  * Shen et al. [2020] Mingzhu Shen, Xianglong Liu, Ruihao Gong, and Kai Han.  Balanced binary neural networks with gated residual.  _International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020.
  * Simonyan and Zisserman [2015] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  _International Conference on Learning Representations (ICLR)_ , 2015\.
  * Williams [1992] Ronald J Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine Learning_ , 1992.
  * Zhou et al. [2016] Shuchang Zhou, Zekun Ni, Xinyu Zhou, He Wen, Yuxin Wu, and Yuheng Zou.  Dorefa-net: Training low bitwidth convolutional neural networks with low bitwidth gradients.  _CoRR_ , abs/1606.06160, 2016.  URL http://arxiv.org/abs/1606.06160.
  * Zhu et al. [2017] Chenzhuo Zhu, Song Han, Huizi Mao, and William J. Dally.  Trained ternary quantization.  _International Conference on Learning Representations (ICLR)_ , 2017\.

##  Appendix A LAR-Net Initialization

In full-precision networks, it is common to use a random initializer for the
weight initialization, e.g., Kaiming Normal [13]. We are interested in
initializing the distributions over discrete weights. O. Shayer and Fetaya
[27] suggested using pretrained continuous deterministic weights for the
distributions initialization. Let $\tilde{W}$ be the normalized pretrained
weights (by dividing the weights in each layer $t$ by the standard deviation
of the weights $\sigma^{(t)}$).
Then, $p(w_{ij}^{(t)}=0)$ is initialized by

$p(w_{ij}^{(t)}=0)=p_{max}-(p_{max}-p_{min})\cdot\lvert\tilde{w}_{ij}^{(t)}\rvert$
(14)

Where $p_{max}$ and $p_{min}$ are hyperparameters (set to 0.05 and 0.95,
respectively, in our experiments). Next, $p(w_{ij}^{(t)}=1\lvert
w_{ij}^{(t)}\neq 0)$ is initialized by

$p(w_{ij}^{(t)}=1\lvert w_{ij}^{(t)}\neq
0)=0.5\cdot(1+\frac{\tilde{w}_{ij}^{(t)}}{1-p(w_{ij}^{(t)}=0)})$ (15)

Then, all the values from equations 14,15 are clipped in range
$[p_{min},p_{max}]$. We found that while this initialization works well also
to train networks with discrete activations, we achieved better results when
we used LR-Net (with only discrete weights) as a baseline. We first trained a
network with discrete weights and then transfer the weights distributions as
initialization to a network with discrete activations as well.

##  Appendix B CIFAR-10 Experimental Details

###  B.1 Hyperparameters for LAR-Net

CIFAR-10 [22] is an image classification benchmark dataset. It consists of
50,000 training images and 10,000 test images distributed among ten different
classes. Each image in the dataset is 32 × 32 pixels in RGB space, and it is
preprocessed by subtracting its mean and dividing by its standard deviation.
During training, we apply a padding of four pixels to each side of the image,
and a random 32 × 32 crop is sampled from the padded image. Furthermore, we
flip images horizontally at random during training. At test time, we evaluate
the original 32 × 32 image without any padding or multiple cropping. The loss
is minimized with Adam [19]. The weight decay parameter is set to $1e-4$. We
noticed that by executing multiple Monte-Carlo samples for every input
element, we can improve the results. In each batch, our model is executed
multiple times for each data, using a random sample, which results in eq. 13.
We use a batch size of 64 and run 2 times each batch ($S=2$). The initial
learning rate is 0.01, and we use a cosine decay learning rate policy for
training. The probability decay parameter is set to $1e-12$, and the
temperature of the Gumbel Soft-Max is 1.2 (fixed for all the training). We
report the test error rate after 300 training epochs. Our discrete weights are
$\\{-1,0,1\\}$, which also brings additional savings. In our method, the
sparsity level depends on the learned $\theta$ and is not pre-defined. We
achieved a sparsity of $44\%\pm 1.1\%$.

##  Appendix C CIFAR-100 Experimental Details

###  C.1 Hyperparameters for LAR-Net

CIFAR-100 [22] is comprised of 60,000 32x32 color images that are grouped into
100 classes, with 600 images for each class. The dataset is divided into
50,000 training images and 10,000 testing images. The loss is minimized with
Adam [19]. The weight decay parameter is set to $1e-4$. We use a batch size of
32 and run 4 times every batch ($S=2$). The initial learning rate is 0.01, and
we use a cosine decay learning rate policy for training. Similar to CIFAR-10,
The probability decay parameter is set to $1e-12$, and the temperature of the
Gumbel Soft-Max is 1.2. We report the test error rate after 300 training
epochs. We achieved a weight sparsity of $43\%\pm 1.0\%$.

###  C.2 Hyperparameters for MPT and IR-Net

We searched over all the next hyper-parameters to find the best results on
CIFAR-100. For IR-Net [28], we used the code from [15], and for MPT [5], we
used the code from [2]. We tried to find the best results with ResNet-18 and
VGG-small. We used Adam optimizer and SGD, and searched the learning rate (LR)
in the range [0.1, 0.001]. We checked different batch sizes [64, 128, 256] and
also different weight decay configurations. For MPT, we used the version in
which the batch normalization parameters are also learned (which they called
MPT+BN). They achieved better results with this version consistently. The
weights were initialized using the Kaiming Normal [13]. Table 3 presents the
best configurations for the models we checked.

Method | Model | Optimizer | LR | Momentum | Weight Decay | Batch | Epochs
---|---|---|---|---|---|---|---
IR-Net | ResNet-18 | SGD | 0.035 | 0.9 | 1e-4 | 128 | 400
IR-Net | VGG-small | SGD | 0.01 | 0.9 | 1e-4 | 128 | 400
MPT + BN | ResNet-18 | Adam | 0.1 |  | 5e-4 | 128 | 400
MPT + BN | VGG-small | Adam | 0.01 |  | 1e-4 | 128 | 400

Table 3: Hyperparameter Configurations for CIFAR-10 Experiments

##  Appendix D Ablation Study

We investigate the performance of LAR-Net without a Batch Normalization layer.
We separate our experiment into the different parts in the BN layer: the
normalization (using $mean$ and $var$) and the affine transform (using
learnable parameters $\gamma$ and $\beta$). Table 4 shows the results using
the different settings. We use no BatchNorm to indicate no Batch Normalization
layer at all, and Affine Transform only to indicate only affine transform
using $\gamma$ and $\beta$, but without the normalization.

Dataset | Model | BatchNorm | Top-1 Acc
---|---|---|---
| ResNet-18 | with BatchNorm | 68.31
CIFAR-100 | ResNet-18 | no BatchNorm | 60.71
| ResNet-18 | Affine Transform only | 61.82
| ResNet-18 | with BatchNorm | 92.34
CIFAR-10 | ResNet-18 | no BatchNorm | 85.65
| ResNet-18 | Affine Transform only | 83.07
Table 4: Ablation Study for LAR-Net


The Open University of<EMAIL_ADDRESS>https://orcid.org/0000-0002-6629-3243 Zeev Nutov [100]Theory of computation
Design and analysis of algorithms

###### Acknowledgements.

# Improved approximation ratio for covering pliable set families

Zeev Nutov

###### Abstract

A classic result of Williamson, Goemans, Mihail, and Vazirani [STOC 1993:
708–717] states that the problem of covering an uncrossable set family by a
min-cost edge set admits approximation ratio $2$, by a primal-dual algorithm
with a reverse delete phase. Recently, Bansal, Cheriyan, Grout, and Ibrahimpur
[ICALP 2023: 15:1–15:19] showed that this algorithm achieves approximation
ratio $16$ for a larger class of set families, that have much weaker
uncrossing properties. In this paper we will refine their analysis and show an
approximation ratio of $10$. This also improves approximation ratios for
several variants of the Capacitated $k$-Edge Connected Spanning Subgraph
problem.

###### keywords:

primal dual algorithm, pliable set family, approximation algorithm

##  1 Introduction

Let $G=(V,E)$ be graph. For $J\subseteq E$ and $S\subseteq V$ let
$\delta_{J}(S)$ denote the set of edges in $J$ with one end in $S$ and the
other in $V\setminus S$, and let $d_{J}(S)=|\delta_{J}(S)|$ be their number.
An edge set $J$ covers $S$ if $d_{J}(S)\geq 1$. Consider the following
problem:

Set Family Edge Cover
Input: A graph $G=(V,E)$ with edge costs $\\{c_{e}:e\in E\\}$ and a set family
${\cal F}$ on $V$.
Output: A min-cost edge set $J\subseteq E$ such that $d_{J}(S)\geq 1$ for all
$S\in{\cal F}$.

In this problem the set family ${\cal F}$ may not be given explicitly, but we
will require that some queries related to ${\cal F}$ can be answered in time
polynomial in $n=|V|$. Specifically, following previous work, we will require
that for any edge set $I$, the inclusion minimal members of the residual
family ${\cal F}^{I}=\\{S\in{\cal F}:d_{I}(S)=0\\}$ of ${\cal F}$ (the family
of sets in ${\cal F}$ that are uncovered by $I$) can be computed in time
polynomial in $n=|V|$. We will also assume that $V,\emptyset\notin{\cal F}$,
as otherwise the problem has no feasible solution.

Agrawal, Klein and Ravi [2] designed and analyzed a primal-dual approximation
algorithm for the Steiner Forest problem, and showed that it achieves
approximation ratio $2$. A classic result of Goemans and Williamson [7] from
the early 90’s shows by an elegant proof that the same algorithm applies for
Set Family Edge Cover with an arbitrary proper family ${\cal F}$. In fact, one
of the main achievements of the Goemans and Williamson paper was defining a
generic class of set families that models a rich collection of combinatorial
optimization problems, for which the primal dual algorithm achieves
approximation ratio $2$. Slightly later, Williamson, Goemans, Mihail, and
Vazirani [11] further extended this result to the more general class of
uncrossable families, by adding to the algorithm a novel reverse-delete phase;
a set family ${\cal F}$ is uncrossable if $A\cap B,A\cup B\in{\cal F}$ or
$A\setminus B,B\setminus A\in{\cal F}$ whenever $A,B\in{\cal F}$. They posed
an open question of extending this algorithm to a larger class of set families
and combinatorial optimization problems. However, for 30 years, the class of
uncrossable set families remained the most general generic class of set
families for which the WGMV algorithm achieves a constant approximation ratio.

Recently, Bansal, Cheriyan, Grout, and Ibrahimpur [4] analyzed the performance
of the WGMV algorithm for a class of set families that arise in several
variants of capacitated network design problems. Recall that ${\cal F}^{I}$
denotes the residual family of ${\cal F}$ w.r.t. an edge set $I$, and that two
sets $A,B$ cross if each of the four sets $A\cap B,V\setminus(A\cup
B),A\setminus B,B\setminus A$ is non-empty.

######  Definition 1.1.

A set family ${\cal F}$ is pliable if for any $A,B\in{\cal F}$ at least two of
the sets $A\cap B,A\cup B,A\setminus B,B\setminus A$ belong to ${\cal F}$. We
say that ${\cal F}$ is $\gamma$-pliable if it has the following additional
property:
Property $(\gamma)$: For any edge set $I$ and sets $S_{1}\subset S_{2}$ in
${\cal F}^{I}$, if an inclusion minimal set $C$ of ${\cal F}^{I}$ crosses each
of $S_{1},S_{2}$, then the set $S_{2}\setminus(S_{1}\cup C)$ is either empty
or belongs to ${\cal F}^{I}$.

Bansal, Cheriyan, Grout, and Ibrahimpur [4] showed that the WGMV algorithm
achieves approximation ratio $16$ for $\gamma$-pliable families, and that
Property $(\gamma)$ is essential – without it the cost of the solution found
by the WGMV algorithm can be $\Omega(\sqrt{n})$ times the cost of an optimal
solution. They also considered applications of their result to several
variants of capacitated network design problems, as follows.

Near Min-Cuts Cover
Input: A graph $G_{0}=(V,E_{0})$, an edge set $E$ on $V$ with costs
$\\{c_{e}:e\in E\\}$, and an integer $k$.
Output: A min-cost edge set $J\subseteq E$ that covers the set family
$\\{\emptyset\neq S\subsetneq V:d_{E_{0}}(S)<k\\}$.

It is known that the set family in Near Min-Cuts Cover is pliable, and [4]
showed that it satisfies Property $(\gamma)$, thus obtaining a
$16$-approximation.

Capacitated $k$-Edge Connected Spanning Subgraph (Cap-$k$-ECSS)
Input: A graph $G=(V,E)$ with edge costs $\\{c_{e}:e\in E\\}$ and edge
capacities $\\{u_{e}:e\in E\\}$, and an integer $k$.
Output: A mini-cost edge set $J\subseteq E$ such that $u(\delta_{J}(S))\geq k$
for all $\emptyset\neq S\subset V$.

One can see that Near Min-Cuts Cover is a particular case of Cap-$k$-ECSS,
when all edges in $E_{0}$ have cost $0$ and capacities $1$, and other edges
have capacity $k$. On the other hand, approximation ratio $\alpha$ for Near
Min-Cuts Cover implies approximation ratio $\alpha\cdot\lceil
k/u_{\min}\rceil$ for Cap-$k$-ECSS, where $u_{\min}$ is the minimum capacity
of an edge; see [4]. This gives approximation ratio $16\cdot\lceil
k/u_{\min}\rceil$ for Cap-$k$-ECSS.

Adjiashvili, Hommelsheim and Mühlenthaler [1] introduced the following related
problem, called $(k,q)$-Flexible Graph Connectivity. Suppose that there is a
subset $U\subseteq E$ of “unsafe” edges, and we want to find the cheapest
spanning subgraph $H$ that will be $k$-connected even if up to $q$ unsafe edge
are removed. Let us say that a subgraph $H=(V,J)$ of $G$ is $(k,q)$-flex-
connected if any cut $\delta_{H}(S)$ of $H$ has at least $k$ safe edges or at
least $k+q$ (safe and unsafe) edges. Namely, we require that $d_{H\setminus
U}(S)\geq k$ or $d_{H}(S)\geq k+q$ for all $\emptyset\neq S\subset V$.
Summarizing, we get the following problem.

$(k,q)$-Flexible Graph Connectivity ($(k,q)$-FGC)
Input: A graph $G=(V,E)$ with edge costs $\\{c_{e}:e\in E\\}$, $U\subseteq E$,
and integers $k,q\geq 0$.
Output: A min-cost subgraph $H$ of $G$ such that $H$ is $(k,q)$-flex-
connected.

For various approximation algorithms for this problem see recent papers [5, 4,
9, 6, 3]. Specifically, Bansal [3] showed that if the problem of covering a
$\gamma$-pliable familily achieves approximation ratio $\alpha$, then
$(k,3)$-FGC admits approximation ratio $6+\alpha$, thus obtaining a
$22$-approximation for $(k,3)$-FGC.

Another generalization of uncrossable families is considered in [8]. A set
family ${\cal F}$ is semi-uncrossable if for any $A,B\in{\cal F}$ we have that
$A\cap B\in{\cal F}$ and one of $A\cup B,A\setminus B,B\setminus A$ is in
${\cal F}$, or $A\setminus B,B\setminus A\in{\cal F}$. One can verify that
semi-uncrossable families are sandwiched between uncrossable and
$\gamma$-pliable families. In [8] it is shown that the WGMV algorithm achieves
the same approximation ratio $2$ for semi-uncrossable families (the proof is
essentially identical to that of [11] for uncrossable families), and are also
given several examples of problems that can be modeled by semi-uncrossable
families that are not uncrossable.

We note that Bansal, Cheriyan, Grout, and Ibrahimpur [4] introduced several
novel analysis techniques of the primal dual method, and derived a relevant
Property $(\gamma)$, that enables to obtain a constant approximation ratio for
pliable families without excluding known applications. Here we refine their
analysis to improve their approximation ratio $16$ as follows.

######  Theorem 1.2.

The Set Family Edge Cover problem with a $\gamma$-pliable set family ${\cal
F}$ admits approximation ratio $10$.

Naturally, this also improves approximation ratios for several applications of
$\gamma$-pliable families discussed in [4, 3] – these improvements are
summarized in Table 1.

| previous | this paper
---|---|---
Near Min-Cuts Cover |  $16$ [4] | $10$
Cap-$k$-ECSS |  $16\cdot\lceil k/u_{\min}\rceil$ [4] | $10\cdot\lceil k/u_{\min}\rceil$
$(k,3)$-FGC |  $22$ for $k$ even, $11+\epsilon$ for $k$ odd [3] |  $16$ for $k$ even
Table 1: Approximation ratios for Near Min-Cuts Cover, Cap-$k$-ECSS, and
$(k,3)$-FGC.

##  2 Proof of Theorem 1.2

Here we prove Theorem 1.2. We start with stating some simple properties of
pliable families. One can see that if an edge $e$ covers one of the sets
$A\cap B,A\cup B,A\setminus B,B\setminus A$ then it also covers one of $A,B$.
This implies the following.

######  Lemma 2.1.

If ${\cal F}$ is pliable or is $\gamma$-pliable, then so is ${\cal F}^{I}$,
for any edge set $I$.

An ${\cal F}$-core is an inclusion minimal member of ${\cal F}$; let ${\cal
C}_{\cal F}$ denote the family of ${\cal F}$-cores.

######  Lemma 2.2.

Let ${\cal F}$ be a pliable set family and let $A\in{\cal F}$ and $C\in{\cal
C}_{\cal F}$. Then either $C\subseteq A$, or $C\cap A=\emptyset$, or
$A\setminus C,A\cup C\in{\cal F}$. Consequently, the members of ${\cal
C}_{\cal F}$ are pairwise disjoint.

######  Proof 2.3.

Since $C\in{\cal C}_{\cal F}$, we can have $C\cap A\in{\cal F}$ only if
$C\subseteq A$, and we can have $C\setminus A\in{\cal F}$ only if $C\cap
A=\emptyset$. In any other case, we must have $A\setminus C,A\cup C\in{\cal
F}$, since ${\cal F}$ is pliable.

We now describe the algorithm. Consider the following LP-relaxation (P) for
Set Family Edge Cover and its dual program (D):

$\displaystyle\begin{array}[]{lllllll}&\hphantom{\bf(P)}&\min&\
\displaystyle\sum_{e\in E}c_{e}x_{e}&\hphantom{\bf(P)}&\max&\
\displaystyle\sum_{S\in{\cal F}}y_{S}\\\ &\mbox{\bf(P)}&\
\mbox{s.t.}&\displaystyle\sum_{e\in\delta(S)}x_{e}\geq 1\ \ \ \ \ \forall
S\in{\cal F}&\mbox{\bf(D)}&\ \mbox{s.t.}&\displaystyle\sum_{\delta(S)\ni
e}y_{S}\leq c_{e}\ \ \ \ \ \forall e\in E\\\ &\hphantom{\bf(P)}&&\ \ x_{e}\geq
0\ \ \ \ \ \ \ \ \ \ \ \forall e\in E&\hphantom{\bf(P)}&&\ \ y_{S}\geq 0\ \ \
\ \ \ \ \ \ \ \ \ \forall S\in{\cal F}\end{array}$

Given a solution $y$ to (D), an edge $e\in E$ is tight if the inequality of
$e$ in (D) holds with equality. The algorithm has two phases.

Phase 1 starts with $I=\emptyset$ an applies a sequence of iterations. At the
beginning of an iteration, we compute the family ${\cal C}={\cal C}_{{\cal
F}^{I}}$ of ${\cal F}^{I}$-cores. Then we raise the dual variables of the
${\cal F}^{I}$-cores uniformly (possibly by zero), until some edge $e\in
E\setminus I$ becomes tight, and add $e$ to $I$. Phase I terminates when
${\cal C}_{{\cal F}^{I}}=\emptyset$, namely when $I$ covers ${\cal F}$.

Phase 2 applies on $I$ “reverse delete”, which means the following. Let
$I=\\{e_{1},\ldots,e_{j}\\}$, where $e_{i+1}$ was added after $e_{i}$. For
$i=j$ downto $1$, we delete $e_{i}$ from $I$ if $I\setminus\\{e_{i}\\}$ still
covers ${\cal F}$. At the end of the algorithm, $I$ is output.

It is easy to see that the produced dual solution is feasible, hence
$\sum_{S\in{\cal F}}y_{S}\leq{\sf opt}$, by the Weak Duality Theorem. We prove
that at the end of the algorithm

$\sum_{e\in I}c(e)\leq 10\sum_{S\in{\cal F}}y_{S}\ .$

As any edge in $I$ is tight, the last inequality is equivalent to

$\sum_{e\in I}\sum_{\delta_{l}(S)\ni e}y_{S}\leq 10\sum_{S\in{\cal F}}y_{S}\
.$

By changing the order of summation we get:

$\sum_{S\in{\cal F}}d_{I}(S)y_{S}\leq 10\sum_{S\in{\cal F}}y_{S}\ .$

It is sufficient to prove that at any iteration the increase at the left hand
side is at most the increase in the right hand side. Let us fix some
iteration, and let ${\cal C}={\cal C}_{{\cal F}^{I}}$ be the family of cores
among the members of ${\cal F}$ not yet covered. The increase in the left hand
side is $\varepsilon\cdot\sum_{C\in{\cal C}}d_{I}(C)$, where $\varepsilon$ is
the amount by which the dual variables were raised in the iteration, while the
increase in the right hand side is $\varepsilon\cdot 10|{\cal C}|$.
Consequently, it is sufficient to prove that $\sum_{C\in{\cal C}}d_{I}(C)\leq
10|{\cal C}|$. As the edges were deleted in reverse order, the set
$I^{\prime}$ of edges in $I$ that were added after the iteration (and
“survived” the reverse delete phase), form an inclusion minimal edge-cover of
the family ${\cal F}^{\prime}$ of members in ${\cal F}$ that are uncovered at
the beginning of the iteration. Note also that $\bigcup_{C\in{\cal
C}}\delta_{I}(C)\subseteq I^{\prime}$. Hence to prove approximation ratio
$10$, it is sufficient to prove the following purely combinatorial statement,
in which due to Lemma 2.1 we can revise our notation to ${\cal
F}\leftarrow{\cal F}^{\prime}$ and $I\leftarrow I^{\prime}$.

######  Lemma 2.4.

Let $I$ be an inclusion minimal cover of a $\gamma$-pliable set family ${\cal
F}$ such that every edge in $I$ covers some $C\in{\cal C}$. Then:

$\sum_{C\in{\cal C}}d_{I}(C)\leq 5\cdot(2|{\cal C}|-1)\ .$ (1)

Furthermore, if ${\cal F}$ is symmetric ($V\setminus S\in{\cal F}$ whenever
$S\in{\cal F}$) then $\sum_{C\in{\cal C}}d_{I}(C)\leq 10\cdot(|{\cal C}|-1)$.

In the rest of this section we prove Lemma 2.4. A set family ${\cal L}$ is a
laminar if any two sets in ${\cal L}$ are disjoint or one of them contains the
other. Let $I$ be an inclusion minimal edge cover of a set family ${\cal F}$.
We say that a set $S_{e}\in{\cal F}$ is a witness set for an edge $e\in I$ if
$e$ is the unique edge in $I$ that covers $S_{e}$, namely, if
$\delta_{I}(S_{e})=\\{e\\}$. We say that ${\cal L}\subseteq{\cal F}$ is a
witness family for $I$ if $|{\cal L}|=|I|$ and for every $e\in I$ there is a
witness set $S_{e}\in{\cal L}$. By the minimality of $I$, there exists a
witness family ${\cal L}\subseteq{\cal F}$. The following lemma was proved in
[4]; we provide a proof for completeness of exposition.

######  Lemma 2.5 ([4]).

Let $I$ be an inclusion minimal cover of a pliable set family ${\cal F}$. Then
there exists a witness family ${\cal L}\subseteq{\cal F}$ for $I$ that is
laminar.

######  Proof 2.6.

Let $A,B\in{\cal F}$ be witness sets of edges $e,f\in I$, respectively. Note
that no edge in $I\setminus\\{e,f\\}$ covers a set from $A\cap B,A\cup
B,A\setminus B,B\setminus A$, as such an edge covers one of $A,B$,
contradicting that $A,B$ are witness sets. Thus all possible locations of such
$e,f$ are as depicted in Fig. 1. We claim that one of the sets $A\cap B,A\cup
B,A\setminus B,B\setminus A$ is a witness set for one of $e,f$. This follows
from the following observations.

  * •

$A\setminus B\notin{\cal F}$ in (a) and $B\setminus A\notin{\cal F}$ in (b),
since these sets are not covered by $e,f$; in both cases, at least one of
$A\cap B,A\cup B$ is in ${\cal F}$, and it is a witness set for one of $e,f$.

  * •

$A\cup B\notin{\cal F}$ in (c), and $A\cap B\notin{\cal F}$ in (d) since these
sets are not covered by $e,f$; in both cases, at least one of $A\setminus
B,B\setminus A$ is in ${\cal F}$, and it is a witness set for one of $e,f$.

This implies that if ${\cal L}$ is a witness family for $I$, then for any
$A,B\in{\cal L}$ there are $A^{\prime}\in\\{A\cap B,A\cup B,A\setminus
B,B\setminus A\\}$ and $B^{\prime}\in\\{A,B\\}$ such that ${\cal
L}^{\prime}=({\cal L}\setminus\\{A,B\\})\cup\\{A^{\prime},B^{\prime}\\}$ is
also a witness family for $I$. W.l.o.g. we may assume that $B^{\prime}=B$, so
${\cal L}^{\prime}$ is obtained from ${\cal L}$ by replacing just one set $A$
by some $A^{\prime}\in\\{A\cap B,A\cup B,A\setminus B,B\setminus A\\}$.

Let us say that two sets $A,B$ overlap if each of the sets $A\cap B,A\setminus
B,B\setminus A$ is non-empty. For $A\subseteq V$ let the overlapping number of
$A$ w.r.t. ${\cal L}$ be $\beta(A,{\cal L})=|\\{B\in{\cal L}:A,B\mbox{
overlap}\\}|$. Let $\beta({\cal L})=\sum_{A\in{\cal L}}\beta(A,{\cal L})$. By
the minimality of $I$ there exists a witness family for $I$. Let ${\cal L}$ be
a witness family for $I$ with $\beta({\cal L})$ minimal. We claim that ${\cal
L}$ is laminar. Suppose to the contrary that there are sets $A,B\in{\cal L}$
that overlap. Let $A^{\prime},B^{\prime},{\cal L}^{\prime}$ be as above, and
w.l.o.g. assume that $B^{\prime}=B$, so ${\cal L}^{\prime}=({\cal
L}\setminus\\{A\\})\cup\\{A^{\prime}\\}$ where $A^{\prime}\in\\{A\cap B,A\cup
B,A\setminus B,B\setminus A\\}$. By [10, Lemma 23.15], if $A$ overlaps some
$B\in{\cal L}$, then for any $A^{\prime}\in\\{A\cap B,A\setminus B,B\setminus
A,A\cup B\\}$, any $S\in{\cal L}$ overlapped by $A^{\prime}$ is also
overlapped by $A$. Consequently, $\beta(A^{\prime},{\cal
L}^{\prime})<\beta(A,{\cal L})$ (since $B$ is not overlapped by $A^{\prime}$),
while $\beta(S,{\cal L}^{\prime})\leq\beta(S,{\cal L})$ for any other set in
${\cal L}^{\prime}$. This implies $\beta({\cal L}^{\prime})<\beta({\cal L})$,
contradicting the choice of ${\cal L}$.

Figure 1: Illustration to the proof of Lemma 2.5.

Augment ${\cal L}$ by the set $V$. Then ${\cal L}$ can be represented by a
rooted tree ${\cal T}$ with node set ${\cal L}$ and root $V$, where the parent
of $S$ in ${\cal T}$ is the smallest set in ${\cal L}$ that properly contains
$S$. The (unique) edge in $I$ that covers $S$ corresponds to the edge in
${\cal T}$ from $S$ to its parent; a node $S$ of ${\cal T}$ corresponds to the
set of nodes in the set $S$ that do not belong to any child of $S$.

######  Definition 2.7.

A set $S\in{\cal L}$ owns a core $C\in{\cal C}$ if $S$ is the inclusion-
minimal set in ${\cal L}$ that contains $C$. We say that $S$ is hollow if it
owns no core. A sequence $(S_{0},S_{1},\ldots,S_{\ell})$ of sets in ${\cal L}$
is called a hollow chain (of length $\ell$) if $S_{1},\ldots,S_{\ell}$ are
hollow, $S_{\ell}\neq V$, and $S_{i-1}$ is the unique child of $S_{i}$ in
${\cal L}$, $i=1,\ldots,\ell$. For each $S_{i}$ let $a_{i}b_{i}$ be the unique
edge in $I$ that covers $S_{i}$, where $a_{i}\in S_{i}$ (possibly
$a_{i}=b_{i-1}$).

We will use for nodes of ${\cal T}$ the the same terminology as for sets in
${\cal L}$; specifically, a node of ${\cal T}$ is hollow if it represents a
hollow set, and a hollow chain in ${\cal T}$ is a path such that all internal
nodes in the path correspond to hollow sets (note that each of these nodes has
degree exactly $2$ in ${\cal T}$). We will need the following well known
statement.

######  Lemma 2.8.

Let $T=(V_{T},E_{T})$ be a tree rooted at $r$ and let $R,W$ be a partition of
$V_{T}$ such that every $w\in W\setminus\\{r\\}$ has at least $2$ children.
Then $|E_{T}|\leq 2|R|-1$. Furthermore, if $r$ has at least $2$ children then
$|E_{T}|\leq 2(|R|-1)$.

######  Proof 2.9.

No leaf of $T$ is in $W$ (note that $r$ is not a leaf), thus $|W|\leq|R|$; the
worse case is when $R$ is the set of leaves of $T$, $r$ has $1$ child, and
every node in $W\setminus\\{r\\}$ has exactly $2$ children. Consequently,
$|E_{T}|\leq|W|+|R|-1\leq 2|R|-1$. If $R$ has at least $2$ children then
$|W|\leq|R|-1$ and we get that $|E_{T}|\leq|W|+|R|-1\leq 2|R|-2$.

If ${\cal T}$ has no hollow chain, then Lemma 2.8 implies that $|{\cal L}|\leq
2|{\cal C}|-1$ and $|{\cal L}|\leq 2(|{\cal C}|-1)$ if ${\cal F}$ is
symmetric. To see this, just let $W$ to be the set of hollow nodes in ${\cal
T}$ in Lemma 2.8, and note that if ${\cal F}$ is symmetric then the root of
${\cal T}$ has at least $2$ children. Since the cores in ${\cal C}$ are
pairwise disjoint (by Lemmas 2.1 and 2.2), every edge contributes at most $2$
to $\sum_{C\in{\cal C}}d_{I}(C)$, hence (since $|I|=|{\cal L}|$)
$\sum_{C\in{\cal C}}d_{I}(C)\leq 2(2|C|-1)$, and $\sum_{C\in{\cal
C}}d_{I}(C)\leq 4(|C|-1)$ if ${\cal F}$ is symmetric.

Suppose now that every hollow chain has length $\ell$. Then the contribution
of the edges $a_{0}b_{0},\ldots,a_{\ell}b_{\ell}$ of each chain is at most
$2(\ell+1)$. If we “shortcut” every maximal hollow chain in ${\cal T}$, we
obtain a tree with at most $2|{\cal C}|-1$ edges. However, every such edge
might be a shortcut of a hollow chain, and thus may contribute $2(\ell+1)$ to
$\sum_{C\in{\cal C}}d_{I}(C)$. As the number of edges after the shortcuts is
at most $2|{\cal C}|-1$, we get $\sum_{C\in{\cal C}}d_{I}(C)\leq
2(\ell+1)(2|{\cal C}|-1)$. Bansal et. al. [4] showed that the maximum possible
length of a hollow chain is $3$, which gives the bound $\sum_{C\in{\cal
C}}d_{I}(C)\leq 8\cdot(2|{\cal C}|-1)<16\cdot|{\cal C}|$.

Let $U=\cup_{C\in{\cal C}}C$. To improve this bound of [4], we will show that
for any hollow chain with $\ell=2,3$, among the relevant nodes
$a_{0},b_{0},\ldots,a_{\ell},b_{\ell}$ at most $3$ belong to $U$. This reduces
the bound on the contribution of every hollow chain from $8$ to $5$, and gives
the bounds $\sum_{C\in{\cal C}}d_{I}(C)\leq 5\cdot(2|{\cal C}|-1)$ and
$\sum_{C\in{\cal C}}d_{I}(C)\leq 5\cdot 2(|{\cal C}|-1)$ if ${\cal F}$ is
symmetric. Specifically, in the next section we will prove the following.

Figure 2: Illustration to Lemma 2.10. Black nodes are in $U$, white nodes are
not in $U$, while gray nodes may or may not be in $U$.

######  Lemma 2.10.

Let $(S_{0},S_{1},\ldots,S_{\ell})$ be a hollow chain and let $a_{i}b_{i}$ be
as in Definition 2.7. Then $\ell\leq 3$ and the following holds, see Fig. 2:

  1. (a)

If $a_{1}\in U$ then $\ell=1$.

  2. (b)

If $b_{0}\in U$ then $\ell\leq 2$; if $\ell=2$ then $a_{1}\notin U$ and
$b_{0},b_{1},a_{2}$ belong to the same core.

  3. (c)

If $b_{1}\in U$ then $\ell\leq 3$; if $\ell=3$ then $a_{1},b_{0},a_{2}\notin
U$ and $b_{1},b_{2},a_{3}$ belong to the same core.

Let us show that Lemma 2.10 implies Lemma 2.4. First, note that at least one
of $a_{1},b_{1}$ is in $U$, since every edge covers some core $C\in{\cal C}$,
by the assumption in Lemma 2.4. Thus $\ell\leq 3$. In the tree representation
${\cal T}$ of ${\cal L}$ let us “shortcut” all maximal hollow chains, see Fig.
3. This means that we replace the chain – the edges of the chain (that
correspond to edges $a_{0}b_{0},\ldots,a_{\ell}b_{\ell}$) and the nodes that
correspond to sets $S_{1},\ldots,S_{\ell}$, by a new “shortcut edge” between
$S_{0}$ and $S_{\ell}$ (for illustration, in Fig. 3 we assume that this new
edge is $a_{0}b_{\ell}$). This operation also can be viewed as identifying
$S_{0}$ with $S_{\ell}$, by removing the nodes in $S_{\ell}\setminus S_{0}$
(and adding to $I$ the edge $a_{0}b_{\ell}$). Every such shortcut edge
contributes at most $5$ (at most $4$ if $\ell=1$) to $\sum_{C\in{\cal
C}}d_{I}(C)$, by Lemma 2.10. Furthermore, the obtained tree has no hollow
chains. By Lemma 2.8, the total number of edges in the obtained tree is at
most $2|{\cal C}|-1$, and at most $2(|{\cal C}|-1)$ if ${\cal F}$ is
symmetric. Each edge of ${\cal T}$ contributes to $\sum_{C\in{\cal
C}}d_{I}(C)$ at most $2$ if it is an ordinary edge and at most $5$ if it is a
shortcut edge. Consequently, $\sum_{C\in{\cal C}}d_{I}(C)\leq 5\cdot(2|{\cal
C}|-1)$, and $\sum_{C\in{\cal C}}d_{I}(C)\leq 10\cdot(|{\cal C}|-1)$ if ${\cal
F}$ is symmetric.

This concludes the proof of Theorem 1.2, provided that we will prove Lemma
2.10, which we will do in the next section.

Figure 3: Illustration of a shortcut of a hollow chain with $\ell=2$. Black
nodes are in $U$, white nodes are not in $U$, while gray nodes may or may not
be in $U$.

##  3 Proof of Lemma 2.10

To prove Lemma 2.10 we will apply Property $(\gamma)$ on the family ${\cal F}$
in Lemma 2.4 (that is in fact a residual family of the original family that we
want to cover), so we restate here Property $(\gamma)$ only for the current
family ${\cal F}$.

Property $(\gamma)$: For any sets $S_{1}\subset S_{2}$ in ${\cal F}$, if an
inclusion minimal set $C$ of ${\cal F}$ crosses each of $S_{1},S_{2}$, then
the set $S_{2}\setminus(S_{1}\cup C)$ is either empty or belongs to ${\cal
F}$.

In what follows, note that every $C\in{\cal C}$ is owned by some some set in
${\cal L}$, since $V\in{\cal L}$. The next two Lemmas is a preparation for
using Property $(\gamma)$.

######  Lemma 3.1.

Let $A$ be the child of a hollow set $S$. Then $S\setminus A$ contains no set
in ${\cal F}$.

######  Proof 3.2.

Suppose to the contrary that $S\setminus A$ contains a set $B\in{\cal F}$.
Then $B$ contains some core $C\in{\cal C}$. Let $S_{C}\in{\cal L}$ be the set
that owns $C$. Then we must have $S_{C}\subseteq S$ and $S_{C}\cap
A=\emptyset$. Consequently, $S$ has at least $2$ children in ${\cal L}$,
contradicting that $S$ is hollow.

######  Lemma 3.3.

Let $C\in{\cal C}$ and let $1\leq i\leq\ell$. If $(C\cap S_{i})\setminus
S_{0}\neq\emptyset$ then $C$ and $S_{i}$ cross.

######  Proof 3.4.

We need to show that $C\cap S_{i},S_{i}\setminus C,C\setminus S_{i}$ are
nonempty and that $S_{i}\cup C\neq V$.

  * •

$C\cap S_{i}\neq\emptyset$ by the assumption.

  * •

$S_{i}\setminus C\neq\emptyset$ since $C$ does not contain $S_{0}$, hence
there is $v\in S_{0}\setminus C\subset S_{i}\setminus C$.

  * •

$C\cup S_{i}\neq V$ since $C\cup S_{i}\in{\cal F}$ by Lemma 2.2 and since
$V\notin{\cal F}$.

It remains to show that $C\setminus S_{i}\neq\emptyset$. Suppose to the
contrary that $C\subseteq S_{i}$. Let $S_{C}$ be the set that owns $C$. Since
$C\subseteq S_{i}$, $S_{C}$ is a descendant of $S_{i}$. Since $C\setminus
S_{0}\neq\emptyset$, $S_{C}$ is not a descendant of $S_{0}$. Thus $S_{C}$ must
be one of the sets $S_{1},\ldots,S_{i}$, which is impossible, since all these
sets are hollow.

Now we will use Property $(\gamma)$. Note that since the cores are pairwise
disjoint (Lemma 2.2), then by the assumption that every edge covers some
$C\in{\cal C}$ (Lemma 2.4), for any edge $uv\in I$ the following holds:
$|\\{u,v\\}\cap U|\geq 1$ and $|\\{u,v\\}\cap C|\leq 1$ for any $C\in{\cal
C}$.

######  Lemma 3.5.

Let $(S_{0},S_{1},\ldots,S_{\ell})$ be a hollow chain and let $a_{i}b_{i}$ be
as in Definition 2.7.

  1. (i)

If $a_{i}\in U$ for some $i\geq 1$ then $\ell=i$.

  2. (ii)

If $b_{i-1}\in U$ for some $i\geq 0$ then $\ell\leq i+1$ and if $\ell=i+1$
then $b_{i-1},b_{i},a_{i+1}$ belong to the same core.

######  Proof 3.6.

For part (i), suppose to the contrary that $S_{i+1}$ exists, see Fig. 4(i).
Let $C\in{\cal C}$ be such that $a_{i}\in C$. By Lemma 3.3, $C$ crosses each
of $S_{i},S_{i+1}$. Note that $b_{i}\notin C$, hence the set
$S_{i+1}\setminus(S_{i}\cup C)$ is non-empty, and thus by Property $(\gamma)$
is in ${\cal F}$. This contradicts Lemma 3.1.

For part (ii), suppose that $C_{i+1}$ exists and let $C\in{\cal C}$ be such
that $b_{i-1}\in{\cal C}$, see Fig. 4(ii). By Lemma 3.3, $C$ crosses each of
$S_{i},S_{i+1}$. If $b_{i}\notin C$ or if $a_{i+1}\notin C$ then the set
$S_{i+1}\setminus(S_{i}\cup C)$ is non-empty, and thus by Property $(\gamma)$
is in ${\cal F}$, contradicting Lemma 3.1. Thus $b_{i},a_{i+1}\in C$. Since
$a_{i+1}\in U$, we get by part (i) that $\ell=i+1$.

Figure 4: Illustration to the proof of Lemma 3.5. (i) $a_{i}\in C$ and
$b_{i}\notin C$. (ii) $b_{i-1}\in C$ and one of $b_{i},a_{i+1}$ not in $C$.

Lemma 2.10 easily follows from Lemma 3.5.

  1. (a)

If $a_{1}\in U$ then by Lemma 3.5(i) we have case (a) of Lemma 2.10, so assume
that $a_{1}\notin U$.

  2. (b)

If $b_{0}\in U$ (and $a_{1}\notin U$), then by Lemma 3.5(ii) we have case (b)
of Lemma 2.10.

  3. (c)

If $a_{1},b_{0}\notin U$ then $b_{1}\in U$. Then either $\ell=2$, or $\ell=3$
and $b_{2},b_{2},a_{3}$ belong to the same core, so we have case (c) of Lemma
2.10.

This concludes the proof of Lemma 2.10, and thus also the proof of Lemma 2.4
and Theorem 1.2 is complete.

## References

  * [1] D. Adjiashvili, F. Hommelsheim, and M. Mühlenthaler.  Flexible graph connectivity.  Mathematical Programming, pages 1–33, 2021.
  * [2] A. Agrawal, P. Klein, and R. Ravi.  When trees collide: An approximation algorithm for the generalized Steiner problem on networks.  SIAM J. on Computing, 24(3):440–456, 1995.
  * [3] I. Bansal.  A constant factor approximation for the $(p,3)$-flexible graph connectivity problem.  CoRR, abs/2308.15714, August 2023.  URL: https://doi.org/10.48550/arXiv.2308.15714.
  * [4] I. Bansal, J. Cheriyan, L. Grout, and S. Ibrahimpur.  Improved approximation algorithms by generalizing the primal-dual method beyond uncrossable functions.  In ICALP, pages 15:1–15:19, 2023.
  * [5] S. C. Boyd, J. Cheriyan, A. Haddadan, and S. Ibrahimpur.  Approximation algorithms for flexible graph connectivity.  In FSTTCS, page 9:1–9:14, 2021.
  * [6] C. Chekuri and R. Jain.  Approximation algorithms for network design in non-uniform fault models.  In ICALP, volume 261, pages 36:1–36:20, 2023.
  * [7] M. X. Goemans and D. P. Williamson.  A general approximation technique for constrained forest problems.  SIAM J. Comput., 24(2):296–317, 1995.
  * [8] Z. Nutov.  Extending the primal-dual 2-approximation algorithm beyond uncrossable set families.  CoRR, abs/2307.08270, 2023.  URL: https://doi.org/10.48550/arXiv.2307.08270, arXiv:2307.08270.
  * [9] Zeev Nutov.  Improved approximation algorithms for some capacitated $k$ edge connectivity problems.  CoRR, abs/2307.01650, July 2023.
  * [10] V. V. Vazirani.  Approximation Algorithms.  Springer, 2001.
  * [11] D. P. Williamson, M. X. Goemans, M. Mihail, and V. V. Vazirani.  A primal-dual approximation algorithm for generalized Steiner network problems.  Combinatorica, 15(3):435–454, 1995.


∎

11institutetext: Research Program in Systems Oncology, 22institutetext:
Faculty of Medicine, University of Helsinki, Finland

Corresponding author. 22email<EMAIL_ADDRESS>
# Bayes in Wonderland! Predictive supervised classification inference hits
unpredictability ††thanks: AA devised the study and developed the theory, AA
and VK advanced the models and implemented the methods, JT provided the
funding, VK, JT, and AA wrote the paper.

Ali Amiryousefi  Ville Kinnula  Jing Tang

(Received: date / Accepted: date)

###### Abstract

The marginal Bayesian predictive classifiers (mBpc) as opposed to the
simultaneous Bayesian predictive classifiers (sBpc), handle each data
separately and hence tacitly assumes the independence of the observations.
However, due to saturation in learning of generative model parameters, the
adverse effect of this false assumption on the accuracy of mBpc tends to wear
out in face of increasing amount of training data; guaranteeing the
convergence of these two classifiers under de Finetti type of exchangeability.
This result however, is far from trivial for the sequences generated under
Partition exchangeability (PE), where even umpteen amount of training data is
not ruling out the possibility of an unobserved outcome (Wonderland!). We
provide a computational scheme that allows the generation of the sequences
under PE. Based on that, with controlled increase of the training data, we
show the convergence of the sBpc and mBpc. This underlies the use of simpler
yet computationally more efficient marginal classifiers instead of
simultaneous. We also provide a parameter estimation of the generative model
giving rise to the partition exchangeable sequence as well as a testing
paradigm for the equality of this parameter across different samples. The
package for Bayesian predictive supervised classifications, parameter
estimation and hypothesis testing of the Ewens Sampling Formula generative
model is deposited on CRAN as PEkit package and
https://github.com/AmiryousefiLab/PEkit.

###### Keywords:

Partition exchangeability Supervised classification Hypothesis testing

###### MSC:

MSC code1 MSC code2 more

††journal: Comp Stat

##  1 Introduction

Under the broad realm of inductive inference, the goal of the supervised
classification is to assign the test objects into $a\,priori$ defined number
of classes learned from the training data Solomonoff . One of the most
applicable machinery that can optimally handle these scenarios is Bayesian
which with a given prior information and accruing observed data, gradually
enhances the precision of the inferred population’s parameters Hand . We
consider here the general supervised classification case where the sets of
species observed for features are not closed $a\,priori$, leaving the
probability of observing new species at any stage non-negative. The de Finetti
type of exchangeability Berlow , seems intractable in these cases.
Nevertheless, one solution is to adhere to a form of partition exchangeability
due to Kingman Kingman2 . Assuming this type of exchangeability for each
competing class, the derivation here shows that given an infinite amount of
data, the simultaneous and marginal predictive classifiers will converge
asymptotically. This is congruent with the similar study under the de Finetti
exchangeability by Corander Corander . Due to the existence of marginal
dependency between the data points, the simultaneous and marginal classifiers
are not necessarily equal. On the other hand, their convergence is not
intuitive due to the complication posed with a priori unfixed set of
observable species. Upon availability of umpteen amount of data however, the
proof presented here justifies the replacement of the marginal classifiers
with the computationally expensive simultaneous ones.

The following section introduces partition exchangeability and a probability
distribution related to it. Hypothesis tests for its parameter are also
introduced. Finally, the predictive classifiers under partition
exchangeability are derived and the algorithms used to implement the
classifiers are presented along with classification results on simulated data
sets from partition exchangeable distributions.

##  2 Partition exchangeability

Assume that number of species related to our feature is unfixed a priori. Upon
availability of the vector of test labels $S$, under the ${\it partition}$
${\it exchangeability}$ framework, we can deduce the sufficient statistic for
each subset of data. To define this statistic, consider the assignment of
arbitrary permutation of integers 1, $\ldots$, —$s_{c}$— to the items in
$s_{c}$ where $n_{c}=|s_{c}|$ is the size of a given class $c$. Introducing
the $I(\ldots)$ indicator function and $n_{cl}=\Sigma_{i\in s_{c}}I(x_{i}=l)$
as the frequency of items in class $c$ having value $l\in{\bf\chi}$ , then in
terms of count in ${\bf x}^{(c)}$ one can write the sufficient statistic as

$\rho_{ct}=\sum_{l=1}^{\infty}I(n_{cl}=t),$ (1)

The vector of sufficient statistic $\rho_{c}=(\rho_{ct})_{t=1}^{n_{c}}$
indicates a partition of the integer $n_{c}$ such that $\rho_{ct}$ is the
frequency of specific feature values that have been observed only $t$ times in
class $c$ of test data. Given the above formulation Kingman , the random
partition is exchangeable if and only if two different sequences having the
same vector of sufficient statistics have the same probabilities of occurring.
According to Kingman’s representation theorem Kingman2 , the probability
distribution of the vector of sufficient statistics under partition
exchangeability will follow the Poisson-Dirichlet($\psi$) (PD) distribution
known also as the Ewens sampling formula Ewens ,

$\displaystyle
p(\boldsymbol{\rho}|\psi)=\frac{n!}{\psi(\psi+1)\ldots(\psi+n-1)}\prod_{t=1}^{n}\Bigg{\\{}{(\frac{\psi}{t})^{\rho_{t}}\frac{1}{\rho_{t}!}}\Bigg{\\}},$
(2)
$\displaystyle\;\forall\>\psi\in\mathbb{R}^{+},\;\;\boldsymbol{\rho}\in\mathfrak{S}_{\rho},$

where,

$\begin{split}\mathfrak{S}_{\rho}=\Bigg{\\{}(\rho_{1},\rho_{2},\ldots\rho_{n})\>|\>\sum_{i=1}^{n}i\rho_{i}=n,&\>\rho_{i}\in\mathbb{N}_{0},\>\\\
&i=1,2,\ldots n\Bigg{\\}}.\end{split}$

For a comprehensive review of the Ewens sampling formula and its history and
further applications, see Crane .

###  2.1 Parameter estimation

The parameter $\psi$ in (2) is called a dispersal parameter, as the higher its
value is the more distinct species will likely be observed in a sample from
this distribution. A Maximum Likelihood Estimate (MLE) for the parameter
$\psi$ can be derived from a sample as demonstrated in Ewens . The estimate
turns out to be the root of the equation

$\sum_{j=1}^{n}\frac{\psi}{\psi+j-1}=\sum_{t=1}^{n}\rho_{t}.$ (3)

The proof of this is provided in Appendix. The sum of sufficient statistics on
the right side of the equation equals the observed number of distinct species
in the sample. The left side of the equation equals, as shown in by Ewens in
Ewens , the expected number of distinct species observed given parameter
$\psi$. Thus, the MLE of $\psi$ is that number for which the observed number
of distinct species equals the expected number of species observed in a sample
of size $n$. There is no closed form solution for the equation for arbitrary
$n$, so the MLE $\hat{\psi}$ has to be numerically searched for. As the right
side of the equation is a strictly increasing function when $\psi>0$, a binary
search algorithm can find the root. In supervised classification, an estimate
for $\psi$ can be calculated for each class in the training data.

###  2.2 Hypothesis testing

A Lagrange multiplier test as defined in Rao  can be used for statistical
testing of a hypothesized parameter $\psi_{0}$ under the null hypothesis
$H_{0}:\>\psi=\psi_{0}$ for a single sample. The test statistic is constructed
as follows:

$S(\psi_{0})=\frac{U(\psi_{0})^{2}}{I(\psi_{0})},$ (4)

where $U$ is the gradient of the log-likelihood $L(\psi)$, and $I$ is the
Fisher information of the distribution. Under the null hypothesis the test
statistic $S$ follows the ${\bf\chi_{1}^{2}}$-distribution. In the case of the
PD distribution, these quantities become

$U(\psi_{0})=\sum_{i=1}^{n}\Bigg{(}\frac{\rho_{i}}{\psi_{0}}-\frac{1}{\psi_{0}+i-1}\Bigg{)}$
(5)

$I(\psi_{0})=\sum_{i=1}^{n}\Bigg{(}\frac{1}{\psi_{0}(\psi_{0}+i-1)}-\frac{1}{(\psi_{0}+i-1)^{2}}\Bigg{)}$
(6)

The proof of this is provided in Appendix.

For a multiple-sample test to infer whether there is a statistically
significant difference in the $\psi$ of each sample, we have a devised a
Likelihood Ratio Test (LRT) Neyman . The null hypothesis of the test is that
there is no difference in the $\psi$ of each of $s$ samples,
$H_{0}:\>\psi_{1}=\psi_{2}=\ldots=\psi_{s}$, and consequently the alternative
hypothesis is the inequality of at least two $\psi$. The test statistic
$\Lambda$ is constructed as follows:

$\Lambda=-2\>log\frac{sup\>\mathcal{L}(\theta_{0})}{sup\>\mathcal{L}(\hat{\theta})}\xrightarrow[]{\mathcal{D}}\chi_{d}^{2},$
(7)

where $\mathcal{L}(\theta_{0})$ is the likelihood function of the data given
the model under the null-hypothesis, and $\mathcal{L(\hat{\theta})}$ is the
unrestricted likelihood of the model. The $sup$ refers to supremum, so the
likelihood is evaluated at the MLE of the parameters. $\Lambda$ asymptotically
converges in distribution to the $\chi_{d}^{2}$, where $d$ equals the
difference in the amount of parameters between the models. When testing the
$\psi$ of $s$ different samples from the PD-distribution with possibly
different sample sizes $n_{s}$, the model under the null hypothesis has one
shared $\psi$, while the unrestricted model has $s$ different dispersal
parameters, ($\psi_{1},\psi_{2},\ldots,\psi_{t}$), so $d=t-1$. The likelihood
$\mathcal{L}$ of multiple independent samples from the PD-distribution is a
product of the density functions for the partitions $\rho$ of those samples,
and under $H_{1}$ the MLE of $\psi$ for each of the samples is evaluated as in
(3) from each sample independently, as the other samples have no effect on the
$\psi$ of a single sample. However, under the null-hypothesis, the samples
share a common $\psi$, the MLE of which according to 2.1 would be estimated by
solving,

$\sum_{i=1}^{s}\sum_{j=1}^{n_{s}}\frac{\psi}{\psi+j-1}=\sum_{i=1}^{s}\sum_{t=1}^{n_{s}}\rho_{st}.$
(8)

As the likelihood of multiple independent samples with identical $\psi$ under
the null-hypothesis is just a product of their likelihoods, the derivative of
the log-likelihood is just a sum of the derivatives of each sample’s
likelihood. Again, the $\psi$ has to be obtained as the root of the above
equation, which can be found with a small modification to the same binary
algorithm as the one we use to determine the MLE of a single PD-distribution.
Having found these MLEs of the $\psi$ under both hypotheses, the likelihood
ratio in the test statistic in equation (7) can then be expressed as,

$\prod_{j=1}^{s}\mathcal{L}(\boldsymbol{\rho}_{j}|\hat{\psi}_{j})\>/\>\prod_{j=1}^{s}\mathcal{L}(\boldsymbol{\rho}_{j}|\hat{\psi}_{j}),$
(9)

where the likelihood function
$\mathcal{L}(\boldsymbol{\rho}_{j}|\hat{\psi}_{j})$ is the likelihood function
of the PD($\boldsymbol{\rho}_{j}|\psi$) distribution for the partition
$\boldsymbol{\rho}$ of the $j$-th sample. As the restricted model in the
numerator can never have a larger likelihood than the unrestricted likelihood
in the denominator, and both are positive real numbers, this ratio is bounded
between 0 and 1.

###  2.3 A note on two-parameter PD

A two-parameter formulation of the distribution presented in Pitman . The
added parameter in PD($\alpha,\psi$) is called the discount parameter. The
role of the parameters is discussed in length in Crane . In short, the
$\alpha$ is defined to fall in the interval $[-1,1]$. When it is positive, it
increases the probability of observing new species in the future proportional
to the amount of already discovered species, while decreasing the probability
of seeing the newly observed species again. When $\alpha$ is negative, the
opposite is true and the number of new species to be discovered is bounded.
The single parameter $PD(\psi)$ is the special case of the two parameter
distribution with $\alpha$ set to 0, $PD(0,\psi)$. The two parameter
distribution is not considered further here, as the estimation of the
parameters becomes a daunting task compared to the simple one parameter
formalization.

###  2.4 Supervised classifiers

Consider the set of $m$ available training items by $M$ and correspondingly
the set of $n$ test items by $N$. For each item, we observe only one feature 1
that can take value from species set ${\bf\chi}=\\{1,2,...,r\\}$. Note that
each number in $\chi$ is represented with one species such that the first
species observed is represented with integer 1, the second species is
represented with integer 2, and so on. On the other hand, $r$ is not known
${\it priori}$ denoting the fact that we are uninformative about all of the
species possible in our population. A training item $i\in M$ is characterized
by a feature $z_{i}$ such that, $z_{i}\in\chi$. Similarly, we have for a test
item $i\in N$ the feature $x_{i}$ such that, $x_{i}\in\chi$. Collections of
the training and test data features are denoted by vectors ${\bf z}$ and ${\bf
x}$, respectively. Furthermore consider that the training data are allocated
into $k$ distinct classes and $T$ is a joint labeling of all the training
items into these classes. Simultaneous supervised classification will assign
labels to all the test data in $N$ in a joint manner. We can consider
partitioning of $N$ test elements into $k$ different classes similar to $T$
such that $S=(s_{1},\ldots,s_{k}),s_{c}\subseteq N,c=1,...,k$ be the joint
labeling of this partition. The $T$ and $S$ structures indicate a partition of
the training and test feature vectors, such that ${\bf z}^{(c)}$ and ${\bf
x}^{(c)}$ represent the subset of training and test items in class
$c=1,...,k,$ respectively. The $S$ denote the space of possible simultaneous
classifications for a given $N$ and so $S\in\mathbb{S}$.

The predictive probability of observing a new feature value of species $j$
given a set of prior observations from the PD distribution is

$p(x_{n+1}=j|{\bf n})=\frac{n_{j}}{N+\psi}$ (10)

where $n_{j}$ is the frequency of species $j$ in the observed set ${\bf n}$.
However, if the value $j$ is of a previously unobserved species, the
predictive probability is

$p(x_{n+1}=j|{\bf n})=\frac{\psi}{N+\psi}.$ (11)

The proof of this arises both from the mechanics of an urn model that Hoppe
showed to generate the Ewens sampling formula, as well as from Karlin . The
predictive probability of a previously unseen feature value is thus higher for
a population with a larger parameter $\psi$ than it is for a population of
equal size with a smaller $\psi$.

Now the product predictive distribution for all test data that is assumes to
be ${\it i.i.d}$ for the marginal classifier under the partition
exchangeability framework becomes:

$\begin{split}p_{M}&=\prod_{c=1}^{k}\prod_{i:S_{i}\in c}p(x_{i}=l_{i}|{\bf
z}^{(c)},T^{(c)},S_{i}=c)\\\ &=\prod_{c=1}^{k}\prod_{i:S_{i}\in
c}\Bigg{(}{\frac{m_{i;cl}}{m_{c}+\hat{\psi}_{c}}}\Bigg{)}^{I(m_{i;cl}\neq
0)}\times\\\
&\quad\Bigg{(}{\frac{\hat{\psi}_{c}}{m_{c}+\hat{\psi}_{c}}}\Bigg{)}^{I(m_{i;cl}=0)}.\end{split}$
(12)

Note that under a maximum $a\>posteriori$ -rule with training data of equal
size for each class, a newly observed feature value that has not been
previously seen in the training data for any class, will always be classified
to the class with the highest estimated dispersal parameter $\hat{\psi}$.
However, each class still has a positive probability of including previously
unseen values based on the observed variety of distinct values previously
observed within the class.

Analogously to the product predictive probability for the case of general
exchangeability in Corander , the product predictive distribution of the
simultaneous classifier under partition exchangeability, fitting the learning
algorithm that is described in a later chapter, is then defined as

$\begin{split}p_{S}&=\prod_{c=1}^{k}\prod_{i:S_{i}\in c}p(x_{i}=l_{i}|{\bf
z}^{(c)},T^{(c)},S_{i}=c)\\\ &=\prod_{c=1}^{k}\prod_{i:S_{i}\in
c}\Bigg{(}{\frac{m_{i;cl}+n_{i;cl}}{m_{c}+n_{i;cl}+\hat{\psi}_{c}}}\Bigg{)}^{I(m_{i;cl}\neq
0)}\times\\\
&\quad\Bigg{(}{\frac{\hat{\psi}_{c}}{m_{c}+n_{i;cl}+\hat{\psi}_{c}}}\Bigg{)}^{I(m_{i;cl}=0)}.\end{split}$
(13)

An asymptotic relationship between these classifiers is immediately apparent
in these predictive probabilities. As the amount of training data in each
class $m_{c}$ increases, the impact of class-wise test data $n_{c}$ becomes
negligible in comparison, and the difference in the predictive probabilities
approaches zero asymptotically. As the classifiers are searching for
classification structures $S$ that optimize the test data predictive
probability given the training data, and the predictive probabilities converge
asymptotically, the classifiers are searching for the same optimal labeling.

However, the classifiers handle values unseen in the training data
differently. This is the situation where $m_{i;cjl}=0$ in the predictive
probabilities $p_{M}$ and $p_{S}$ The marginal classifier’s predictive
probability for the test data is always maximized by assigning such a value
into the class with the highest $\hat{\psi}$. The simultaneous classifier,
however, considers the assignment of other instances of an unseen value as
well. This can lead to optimal classification structures where different
instances of an unseen value are classified into different classes. Thus, the
convergence of the test data predictive probabilities of the marginal and
simultaneous classifiers is not certain in the presence of unseen values. In
practice, though, as the amount of training data $m$ tends to infinity, the
probability of observing new values from the PD($\psi$)-distribution presented
in equation (11) tends to 0: $\psi/(\psi+m)\>\rightarrow\>0$ as
$m\rightarrow\infty$. Unexpected values would be very rare and would have a
minimal effect on the classifications made by the different classifiers.

###  2.5 Algorithms for the predictive classifiers

Table 1: The item-wise 0-1 classification error for the marginal and simultaneous classifiers, as well as the 0-1 difference between the predicted labels between the two classifiers. | $\hat{p}_{M}$ | $\hat{p}_{S}$ | $|\hat{p}_{M}-\hat{p}_{S}|$
---|---|---|---
$m=2*10^{6},n=2000,k=3,$ ${\bf\psi}=1,2$ | 0.491 | 0.491 | 0
$m=1000,n=2000,k=3,$ ${\bf\psi}=1,10,50$ | 0.3408 | 0.2823 | 0.0626
$m=1*10^{5},n=2000,k=3,$ ${\bf\psi}=1,10,50$ | 0.2768 | 0.2758 | 0.0010
$m=1000,n=2000,k=5,$ ${\bf\psi}=1,100,1000,5000,10000$ | 0.7535 | 0.5865 | 0.946
$m=2*10^{6},n=2000,k=5,$ ${\bf\psi}=1,100,1000,5000,10000$ | 0.434 | 0.364 | 0.115

In this chapter the learning algorithms defined in Corander  are described.
The predictive probabilities used in these algorithms are defined above up to
an unknown normalizing constant that, however, can be omitted as it doesn’t
affect the probability maximization step in the classification.

The marginal classifier is computationally attractive, as each test data point
is individually classified according to a maximum ${\it a\;posteriori}$ rule:

$\hat{p}_{M}:\underset{c=1,\ldots,k}{argmax}\;p(S_{i}=c|{\bf
z}^{(M)},x_{i},T),$ (14)

For the simultaneous classification algorithm, define classification structure
$S_{c}^{(i)}$ to be an identical classification structure to $S$ with the item
$i$ reclassified to class $c$. The greedy deterministic algorithm is then
defined as

$\hat{p}_{S}:\underset{S\in\mathbb{S}}{argmax}\;p(S|{\bf x}^{(N)},{\bf
z}^{(M)},T),$ (15)

  1. 1.

Set an initial $S_{0}$ with the marginal classifier algorithm $\hat{p}_{M}$.

  2. 2.

Until $S$ remains unchanged between iteration, do for each test item $i\in N$:

$\hat{p}_{S}:\underset{S\in\\{S_{1}^{(i)},\ldots,S_{k}^{(i)}\\}}{argmax}\;p(S|{\bf
x}^{(N)},{\bf z}^{(M)},T)$ (16)

The simultaneous classifier thus works the same as the marginal classifier,
but for each test item the potential labeling of other test items is also
taken into account.

##  3 Numerical illustrations of classifier performance and convergence

Figure 1: The classification error of the marginal and simultaneous
classifiers with data sets from PD-distribution with $\psi\in(1,10,50)$, as
well as the convergence of their labelings. Rows 2 and 3 in 1 are included in
the figure.

To study the classification performance and convergence of the two different
classifiers, we simulated both training and test data sets from the PD
distribution with a generative urn model described in Zabell  called the De
Morgan process. We varied the amount of training and test data, the values of
the distribution parameter $\psi$, as well as the amount of classes $k$. A
large training data set of 4 millino data points was created for each $\psi$.
A small sample of size 1000 was used to train the model first, and more of it
was subsequently added until finally half of the whole data set was used as
training data. A test data set created with the same parameters was kept
constant for classification with all training data samples. An excerpt of the
results can be found in Table 1, as well as Figure 1.

The table and the figures as well as the results not presented here show that
the simultaneous classifier performs better than the marginal classifier
especially as the amount of classes is higher and training data is small. As
the amount of training data used increases, the predictions of the classifiers
converge, as the extra information in the test data used by the simultaneous
classifier becomes negligible. This behaviour is illustrated in Figure 1.
Additionally, the first row of the table shows that the classifier fails on
the binary classification task with data sets created with $\psi$ of 1 and 2
respectively. This is because the classes created with $\psi$ of such similar
magnitude will be distributed very similarly. Thus the classification results
are no better than random guessing. The data is in fact so homogeneous in this
case, that the labelings of the two different classifiers already converge
with test data of size 1000.

In conclusion, the results support the hypothesis that the marginal and
simultaneous classifiers converge in their labelings with enough training
data. This justifies the use of the marginal classifier in place of the more
accurate simultaneous one with large data. Still, the simultaneous classifiers
is more accurate with smaller training data, as it benefits from the
information in the test data.

##  4 Discussion

Previously unseen or unanticipated feature values are challenging scenarios in
standard Bayesian inductive inference. Under general exchangeability as
formalized by de Finetti, upon such an observation the entire alphabet of
anticipated feature values must be retrospectively changed. The classifiers
introduced in this article are, however, equipped to update their predictions
autonomously based on a fully probabilistic framework. The classification of
previously unseen feature values is handled through the use of a parameter
learned from the training data, instead of for example using an uninformative
Dirichlet-prior.

The superiority in classification accuracy of labeling the test data
simultaneously instead of one by one was also shown under partition
exchangeability. The assumption that test data is $i.i.d.$ is obviously an
unrealistic one. The two classifiers considered here only converge in
prediction as the amount of training data approaches infinity. However, The
computational cost of the simultaneous classifier increases exponentially as
the amount of test data increases, making it unfeasible to use for large-scale
prediction. Additionally, the algorithm presented here is only capable of
arriving at local optima. Further research could be directed at the Gibbs
sampler -assisted algorithm presented in Corander , although the convergence
of such algorithms with large test data sets is also uncertain. An
implementation for the supervised classifiers considered here could also be
devised for the two parameter distribution Poisson-Dirichlet($\alpha,\psi$).

## Acknowledgements

We thank ERC grant No. 716063 for the financial support.

## References

  * (1) Barlow, R. E. (1992), Introduction to de Finetti (1937) foresight: Its logical laws, its subjective sources. Spr. Nat., 127–133.
  * (2) Corander, J., Cui, Y., Koski, T., and Siren, J. (2011). Have I seen you before? Principles of Bayesian predictive classification revisited. Springer, Stat. Comput. 23, 59–73.
  * (3) Crane, Harry. (2016). The Ubiquitous Ewens Sampling Formula. Statistical Science. 31. 10.1214/15-STS529.
  * (4) Ewens, W., (1972) The Sampling Theory of Selectively Neutral Alleles, Theoretical Population Biology 3, 87-112.
  * (5) W.J. Ewens,(1979),Testing the generalized neutrality hypothesis, Theoretical Population Biology, 15, 205-216, https://doi.org/10.1016/0040-5809(79)90035-2.
  * (6) Hand, D. J. and Yu, K., (2001), Idiot’s Bayes: Not so stupid after all?, Int. Stat. Rev. 69, (3), 385.
  * (7) Hoppe, F. M. (1984). Polya-like urns and the Ewens sampling formula, Journal of Mathematical Biology, 20, 91-94.
  * (8) Karlin, S., and McGregor, J. (1972). Addendum to a paper of Ewens, Theoretical Population Biology 3, 113-116.
  * (9) Kingman, J. F. C., (1977). The population structure associated with the Ewens sampling formula. Theoretical Population Biology 11, 274–283.
  * (10) Kingman, J. F. C., (1978), Random partitions in population genetics, Proc. of the Royal Society A: Math, Phys. and Eng. Sci. 361, 1704, 1–20.
  * (11) Neyman, J., Pearson, E. S. (1933). On the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical Or Physical Character, 231(694-706), 289-337. doi:10.1098/rsta.1933.0009.
  * (12) Pitman, J., Yor, M., (1997) ”The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator.” Ann. Probab. 25 (2) 855 - 900.
  * (13) Radhakrishna Rao, C, (1948), Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation. Mathematical Proceedings of the Cambridge Philosophical Society, 44(1), 50-57. doi:10.1017/S0305004100023987
  * (14) Solomonoff, R., (1964), A formal theory of inductive inference. Inf. and Ctrl. 7, 1–22.
  * (15) Zabell, S. L., (1992). Predicting the unpredictable. Springer, Synthese, 90, 205–232.
  * (16) Zhou, X., Huang, J., & Wu, X. (2017). Estimation of Poisson-Dirichlet parameters with monotone missing data. Mathematical Problems in Engineering, 2017, 1-10. doi:10.1155/2017/7892507.

## Appendix

### Maximum Likelihood Estimate

Here we provide the technical details for deriving the MLE in equation (3) as
well as the components $U$ and $I$ in the Lagrange Multiplier test described
in equations (5) and (6). The task is to find the first derivative of the
logarithm of the Ewens sampling formula and its root in equation (2), as well
as the second derivative needed for the Fisher information.

$\displaystyle\begin{split}L(\psi)&=log(n!)+\sum_{i=1}^{n}\\{-log(\psi+i-1)\,+\\\
&\rho_{i}\,log\,\psi-\rho_{i}\,log\,\rho_{i}-log(\rho_{i}!)\\}\\\
=>l(\psi)&=\sum_{i=1}^{n}\\{-log(\psi+i-1)\,+\rho_{i}\,log\,\psi\\}\end{split}$

$\displaystyle\begin{split}l^{\prime}(\psi)&=\sum_{i=1}^{n}\Bigg{(}-\frac{1}{\psi+i-1}\,+\frac{\rho_{i}}{\psi}\Bigg{)}\\\
&=U(\psi)\end{split}$

The MLE is found by finding the root of the equation:

$\displaystyle\begin{split}\sum_{i=1}^{n}\Bigg{(}-\frac{1}{\hat{\psi}+i-1}\,+\frac{\rho_{i}}{\hat{\psi}}\Bigg{)}=0\\\
\sum_{j=1}^{n}\frac{\hat{\psi}}{\hat{\psi}+j+1}=\sum_{t=1}^{n}\rho_{t}\end{split}$

### Lagrange Multiplier Test

According to Ewens in Ewens , the left side of the above equation equals the
expected number of unique values observed with this $\psi$ and this sample
size $n$, while the right side is the observed number of unique values:

$\displaystyle\begin{split}E[\sum_{t=1}^{n}\rho_{t}|\psi,n]=\sum_{i=1}^{n}\frac{\hat{\psi}}{\hat{\psi}+i-1}\end{split}$

This is needed for the Fisher information $I$, along with the second
derivative of $l(\psi)$:

$\displaystyle\begin{split}l^{\prime\prime}(\psi)&=\sum_{i=1}^{n}\Bigg{(}\frac{1}{(\psi+i-1)^{2}}\,-\frac{\rho_{i}}{\psi^{2}}\Bigg{)}\\\
&=\sum_{i=1}^{n}\frac{1}{(\psi+i-1)^{2}}-\frac{k}{\sum_{i=1}^{n}\psi^{2}},\end{split}$

where $k$ denotes the observed amount of unique values in the sample. The
Fisher information then becomes:

$\displaystyle\begin{split}I(\psi)&=-E[l^{\prime\prime}(k;\psi)|\psi]\\\
&=\frac{E[k|\psi]}{\sum_{i=1}^{n}\psi^{2}}-\sum_{i=1}^{n}\frac{1}{(\psi+i-1)^{2}}\\\
&=\sum_{i=1}^{n}\Bigg{(}\frac{1}{\psi(\psi+i-1)}-\frac{1}{(\psi+i-1)^{2}}\Bigg{)},\end{split}$

where the expectation of $k$ is as described from above.