# Stochastic Subnetwork Annealing: A Regularization Technique for Fine Tuning
Pruned Subnetworks

Tim Whitaker  Department of Computer Science
Colorado State University
Fort Collins, CO, USA
<EMAIL_ADDRESS>
Darrell Whitley  Department of Computer Science
Colorado State University
Fort Collins, CO, USA
<EMAIL_ADDRESS>

###### Abstract

Pruning methods have recently grown in popularity as an effective way to
reduce the size and computational complexity of deep neural networks. Large
numbers of parameters can be removed from trained models with little
discernible loss in accuracy after a small number of continued training
epochs. However, pruning too many parameters at once often causes an initial
steep drop in accuracy which can undermine convergence quality. Iterative
pruning approaches mitigate this by gradually removing a small number of
parameters over multiple epochs. However, this can still lead to subnetworks
that overfit local regions of the loss landscape. We introduce a novel and
effective approach to tuning subnetworks through a regularization technique we
call Stochastic Subnetwork Annealing. Instead of removing parameters in a
discrete manner, we instead represent subnetworks with stochastic masks where
each parameter has a probabilistic chance of being included or excluded on any
given forward pass. We anneal these probabilities over time such that
subnetwork structure slowly evolves as mask values become more deterministic,
allowing for a smoother and more robust optimization of subnetworks at high
levels of sparsity.

##  I Introduction

Deep neural networks have seen a steady increase in size as large amounts of
compute and high quality data have become more accessible. Models with
billions of parameters are becoming commonplace and demonstrating incredible
performance across a wide range of difficult machine learning tasks. However,
these models also bring important challenges related to computational needs,
storage cost, and training efficiency. The resource-intensive nature of these
large networks have spurred a growing interest in techniques that can reduce
the size and computational complexity associated with training and deploying
these models.

One of the most popular methods used to compress large networks is weight
pruning. It’s long been known that you can remove a significant number of
parameters from these trained models, and after a small number of epochs of
continued training, they can maintain or even exceed the performance of the
full size network [1, 5]. While saliency metrics are often used to prune the
most unnecessary weights, even random sampling has been show to produce
accurate models at moderate levels of sparsity [27].

Several researchers have investigated ways tune these subnetworks more
efficiently and with better accuracy. One of the most effective techniques
used in state of the art pruning methods involves iterative pruning/tuning
cycles [1]. The key insight being that pruning too many parameters at once can
lead to drastic performance collapse as subnetworks get stuck in lower
performing regions of the optimization landscape. Pruning a small number of
parameters over several epochs allows the model to better adapt to the
changing network structure. However, these iterative methods can still result
in subnetworks that are prone to overfitting local optimization regions.

We introduce a novel regularization approach for fine tuning these subnetworks
by leveraging stochasticity for the network structure. Rather than using fixed
subnetworks and discrete pruning operations, we instead represent subnetworks
with probability matrices that determine how likely it is that a parameter is
retained on any given forward pass. The probability matrices are then adjusted
during the tuning phase through the use of dynamic annealing schedules which
allows for a subnetwork to be slowly revealed over several epochs. The
probabilistic inclusion of extra parameters early in the tuning process allows
for gradient information to bleed through into the target subnetwork,
encouraging robust adaptation and avoiding the drastic performance collapse
observed with one-shot pruning methods.

We conduct a large scale ablation study with ResNets and Vision Transformers
to explore the efficacy and dynamics of several hyperparameters and their
effects on convergence, including: degree of stochasticity, number of
annealing epochs, amount of sparsity, constant/phasic learning rate schedules,
and random/saliency based parameter selection strategies. Our experiments
demonstrate significant improvements over both one-shot and iterative pruning
methods across several tuning configurations with especially large
improvements in highly sparse convolutional subnetworks (95-98%).

We additionally explore how Stochastic Subnetwork Annealing can be leveraged
in a benchmark low-cost ensemble learning method that generates diverse child
networks through random pruning of a trained parent network. Implementing our
technique to tune the child networks results in better ensemble generalization
on benchmark image classification tasks, improving on the training
compute/accuracy Pareto Frontier for efficient ensemble methods.

##  II Background

Subnetworks are represented with binary bit mask matrices with the same
dimensions of the weight matrices for each layer in the parent network.
Consider a weight matrix $W\in\mathbb{R}^{m\times n}$ representing the weights
of a particular layer in a neural network. We introduce a matrix
$M\in\\{0,1\\}^{m\times n}$ with the same dimensions as W. The elements of M
are binary values, where $M_{ij}=1$ if the corresponding weight $W_{ij}$ is
retained, and $M_{ij}=0$ if the weight is masked. The mask is generated with
an arbitrary discrete stochastic process $\phi$. The subnetwork weights
$\hat{W}$ can then be computed as the element-wise (Hadamard) product of W and
M.

$\displaystyle M$ $\displaystyle\sim\phi^{m\times n}$ (1)
$\displaystyle\hat{W}$ $\displaystyle=W\circ M$ (2)

The topology of the subnetwork can be further described by the granularity in
which parameters are masked. This granularity refers to the unstructured or
structured distribution of masked weights, where unstructured methods refer to
weight-level or connection-level masking and structured methods refer to
neuron-level, channel-level, or layer-level masking. Removing entire rows,
columns, or blocks from a layer’s weight matrix can be effective at reducing
computational complexity as the reduced size of the weight matrix can be
leveraged for hardware optimizations. This is more difficult to achieve with
the sparse matrices resulting from unstructured masking. However, unstructured
masking tends to result in networks with better generalization as the number
of masked parameters increases [1]. We choose to use unstructured masking for
our experiments with Stochastic Annealing due to the higher sparsity
potential, however our work is fully compatible with structured methods.

It’s important to also consider the distribution of masked weights throughout
non-homogenous networks. It’s common for layer configurations in deep neural
networks to vary significantly in size and shape. Severe pruning of small
layers may result in bottlenecks that restrict gradient flow. The distribution
of masked weights can be controlled with global and local masking methods
where global methods are applied uniformly across the entire network while
local methods are applied independently within each layer or sub-region.
Global methods tend to result in higher compression rates as more parameters
from the larger layers are removed whereas local methods offer more fine-
grained control and reduced variance [1].

Iterative pruning has been shown to be highly effective at improving accuracy
of subnetworks compared to one-shot pruning [23, 5, 10, 6]. Instead of pruning
all of the weights at once, instead an iterative cycle is implemented where a
small number of weights are pruned, the network is tuned, and this repeats
until a target sparsity level is reached. Iterative pruning results in a less
destructive effect on network performance which can allow for greater levels
of sparsity at improved accuracy.

###  II-A Optimizing Pruned Subnetworks

Continued training of the subnetwork is crucial in order to recover the lost
accuracy from pruning. Only a small number of epochs are generally needed as
subnetworks inherited from trained networks tend to converge quickly. The
standard practice for fine-tuning subnetworks involves training with a small
constant learning rate consistent with the final phase of the parent network’s
training procedure [23, 28, 10].

The Lottery Ticket Hypothesis introduced the concept of weight rewinding,
where both the network weights and the learning rate is rewound to a previous
state $t$ epochs ago. Training continues from this previous state but with the
new subnetwork structure fixed [5]. Subsequent research has shown that
rewinding both the weights and the learning rate is less effective than
rewinding only the learning rate itself [35].

The efficacy of learning rate rewinding has been tangentially observed in many
optimization papers where a decaying learning rate schedule improves the
training efficiency and convergence quality in deep networks [44]. The one-
cycle policy is one such schedule, which consists of a warm-up phase that
anneals from a small learning rate to a large learning rate, followed by a
cosine annealed cool-down to a very small value. This schedule can lead to a
phenomenon called super-convergence, where network training is greatly
accelerated on some datasets [36]. This kind of phasic learning rate schedule
is effective when the number of tuning epochs is relatively small, as seen in
several efficient ensemble methods, as large learning rates help to encourage
more movement in parameter space before using small learning rates to converge
to local optima [36, 15, 7, 42].

###  II-B Probabilistic Subnetworks

Several researchers have explored the connection between stochasticity and
sparse network structure. The most prominent example being Dropout, which
randomly disables neurons on every forward pass during training according to a
layerwise constant probability value [37]. Variational and Adaptive Dropout
methods introduced a probabilistic approach to determining the importance of
neurons by learning optimal dropout rates during training [32, 17, 21].
Sparsity inducing regularization techniques have been used to prune networks
during training by incorporating additional loss terms and learning objectives
[23, 40, 24]. Probabilistic gating techniques have also been used to route
signals through sparse network paths and is effective in multi-task or
continual learning environments [14, 4]. While these methods leverage
stochasticity to find sparse network structures, they often require special
training procedures, network architectures, or are limited in that they learn
only a single subnetwork throughout the entirety of training. Stochastic
Subnetwork Annealing has a completely distinct methodology and purpose by
instead focusing on improving the optimization of any randomly sampled
subnetwork from any fully trained model in a small number of epochs.

##  III Stochastic Subnetwork Annealing

In pruning literature, sparse network structures are generally static and
represented with binary bit masks. We propose a model of representing neural
subnetworks with probabilistic masks, where each parameter is assigned a score
that determines how likely it is that the parameter will be retained on any
given forward pass. This introduces stochasticity into the subnetwork sampling
process, which can act as a form of implicit regularization analagous to a
reverse dropout, where parameters that would have been pruned have a chance to
activate. This technique encourages exploration of a larger space of
subnetworks during the fine-tuning phase, preventing it from becoming overly
reliant on a single fixed topological configuration, resulting in more robust
and generalized subnetworks.

Consider a weight matrix $W\in\mathbb{R}^{m\times n}$ representing the weights
of a particular layer in a neural network. We introduce a probability matrix
$P\in\mathbb{R}^{m\times n}$ containing scores that represent the probability
that a parameter $W_{ij}$ will be will be masked on any given forward pass.
The subnetwork mask $M\in\\{0,1\\}^{m\times n}$ is determined with a Bernoulli
realization of the probability matrix $P$.

We then anneal these probability values over some number of epochs such that
subnetworks are slowly ”revealed” throughout fine-tuning. This is done by
introducing an annealing schedule, where the probability values for each
parameter slowly move towards $0$ or $1$ depending on the target subnetwork
sparsity. At the beginning of the training process, a high level of
stochasticity is desirable as it encourages exploration of the weight space
which may help to prevent overfitting and search for more effective subnetwork
configurations. As training progresses, the stochasticity should gradually
decrease via some decay function (linear, cosine, exponential, etc.) to allow
for stable convergence over several epochs.

$\displaystyle y_{lin}(t)$
$\displaystyle=\tau_{max}-(\tau_{max}-\tau_{min})\left(\frac{t}{T_{max}}\right)$
(3) $\displaystyle y_{cos}(t)$
$\displaystyle=\tau_{min}+\frac{1}{2}(\tau_{max}-\tau_{min})\left(1+cos\left(\frac{t}{T_{max}}\right)\right)$
(4)

where $(\tau_{max},\tau_{min})$ are the initial and final leaning rates, $t$
is the current time step, and $T_{max}$ is the total number of annealing
steps.

###  III-A Random Annealing

The probability matrix can be generated in arbitrary ways. For example, a
uniform distribution $P\sim U([0,1])^{m\times n}$ can be used to randomly
assign probabilities to each parameter. All parameters with a value less than
the sparsity target will then anneal towards 0 while parameters with a
probability value greater than the sparsity target will anneal towards 1.

0$S$1

With this implementation, the mean activation of parameters at the beginning
of tuning will be 50% regardless of the target subnetwork sparsity. Other
distributions can offer more fine grained control over network structure. For
example, assume that a binary matrix $X\in\\{0,1\\}^{m\times n}$ is randomly
generated and used to index into a probability matrix $P\in\mathbb{R}^{m\times
n}$. Using this index matrix $X$, we can sample from Gaussian distributions
with different means and variances.

$P=\begin{cases}P_{ij}\sim\mathcal{N}(\mu_{1},\sigma_{1}^{2}),\ \text{if}\
X_{ij}=0\\\ P_{ij}\sim\mathcal{N}(\mu_{2},\sigma_{2}^{2}),\ \text{if}\
X_{ij}=1\end{cases}$

$0$$1$$x$$y$

This approach to constructing multi-modal distributions can allow for
interesting formulations of stochastic subnetworks. Future work may find
natural applications to multi-task learning, where certain groups of
parameters can be strongly correlated with each other while allowing for
overlap with other task specific subnetworks. This stochastic overlap may
encourage shared portions of the network to learn generalized features,
avoiding the problem in typical network-splitting approaches where task
specific subnetworks become increasingly narrow with degraded performance.

###  III-B Temperature Annealing

Temperature scaling may be applied to binary matrices in order to allow for an
even application of stochasticity, reducing some variance relative to random
annealing. Assume that some binary matrix $X\in\\{0,1\\}^{m\times n}$ is
generated to represent a target subnetwork. A temperature scaling constant
$\tau$ is introduced such that the values in $X$ with a 1 are decayed by
$\tau$ and the values with 0 are increased by $\tau$. The mask is then
determined on every forward pass with an altered probability according to
$\tau$. Altering the initial value of $\tau$ allows for a much more controlled
approach to noise injection during the early phases of tuning, which can be
highly desirable when the number of total tuning epochs is limited.

A variation of temperature scaling, where tau is only applied to parameters
that are not a part of the target subnetwork, will be shown to be a highly
effective form of regularization analagous to a reverse dropout. That is, the
subnetwork is always active, but parameters that would have been pruned away
now have a chance to pop back in during tuning. This formulation is less
destructive than full temperature scaling as the target subnetwork will be
optimized on every training step. Allowing other parameters to become active
allows for gradient information to contribute to optimization of the target
subnetwork which can help to encourage avoidance of local minima.

$P=\begin{bmatrix}1&1&0+\tau&1\\\ 0+\tau&1&0+\tau&0+\tau\\\
1&0+\tau&0+\tau&1\\\ 0+\tau&0+\tau&1&1\end{bmatrix}$

##  IV Experiments

###  IV-A Ablations

Figure 1: Ablations exploring how the initial temperature and the number of annealing epochs affect convergence behavior. The rightmost graph displays the best performing annealing configuration for each sparsity and plots it against the one-shot pruning baseline. Stochastic annealing outperformed the other methods for each sparsity, with more drastic improvements appearing in extremely sparse networks. Table I: Comparison between various baseline approaches when used with random/magnitude pruning on CIFAR-100. We report the best accuracy for each method at various levels of target sparsity when tuned with both a constant learning rate and a one-cycle learning rate schedule. The parent network has a baseline accuracy of 71.5%. Our annealing strategies consistently outperform established one-shot and iterative pruning methods. Resnet-18/CIFAR-100 |  Method | 50% | 70% | 90% | 95% | 98%
---|---|---|---|---|---|---
Random/Constant |  One-Shot Baseline | $64.0\pm 0.07$ | $63.4\pm 0.13$ | $54.9\pm 0.12$ | $48.4\pm 0.24$ | $41.1\pm 0.47$
|  Iterative Pruning | $65.1\pm 0.10$ | $63.4\pm 0.05$ | $57.8\pm 0.06$ | $49.8\pm 0.22$ | $42.8\pm 0.26$
|  Random Annealing | $65.3\pm 0.17$ | $63.6\pm 0.15$ | $58.9\pm 0.17$ | $54.5\pm 0.16$ | $45.7\pm 0.22$
|  Temperature Annealing | $65.6\pm 0.09$ | $64.3\pm 0.12$ | $59.5\pm 0.02$ | $55.4\pm 0.19$ | $46.5\pm 0.15$
Random/One-Cycle |  One-Shot Baseline | $70.1\pm 0.07$ | $68.9\pm 0.19$ | $61.8\pm 0.07$ | $56.8\pm 0.21$ | $49.3\pm 0.27$
|  Iterative Pruning | $70.1\pm 0.03$ | $69.2\pm 0.11$ | $63.9\pm 0.05$ | $57.4\pm 0.12$ | $49.8\pm 0.51$
|  Random Annealing | $70.2\pm 0.09$ | $69.2\pm 0.12$ | $64.5\pm 0.07$ | $59.7\pm 0.16$ | $52.7\pm 0.30$
|  Temperature Annealing | $70.6\pm 0.06$ | $69.6\pm 0.09$ | $64.9\pm 0.08$ | $60.2\pm 0.08$ | $53.3\pm 0.29$
Magnitude/Constant |  One-Shot Baseline | $67.1\pm 0.16$ | $67.3\pm 0.37$ | $65.7\pm 0.24$ | $63.3\pm 0.14$ | $54.4\pm 0.24$
|  Iterative Pruning | $68.2\pm 0.09$ | $67.9\pm 0.04$ | $65.9\pm 0.07$ | $63.3\pm 0.02$ | $56.2\pm 0.11$
|  Temperature Annealing | $68.5\pm 0.03$ | $68.3\pm 0.03$ | $66.3\pm 0.02$ | $63.4\pm 0.11$ | $56.6\pm 0.19$
Magnitude/One-Cycle |  One-Shot Baseline | $70.9\pm 0.12$ | $70.1\pm 0.26$ | $69.7\pm 0.06$ | $67.5\pm 0.16$ | $60.2\pm 0.10$
|  Iterative Pruning | $71.1\pm 0.04$ | $71.1\pm 0.06$ | $69.9\pm 0.16$ | $67.6\pm 0.08$ | $60.8\pm 0.19$
|  Temperature Annealing | $71.5\pm 0.06$ | $71.4\pm 0.02$ | $70.3\pm 0.10$ | $67.8\pm 0.02$ | $61.5\pm 0.02$

We begin with an exploration of several stochastic subnetwork annealing
configurations with the goal of investigating how different hyperparameters
impact the efficiency and quality of subnetwork convergence, compared to
established one-shot and iterative pruning techniques.

We use the benchmark CIFAR-10/CIFAR-100 datasets to conduct our explorations
[18]. CIFAR consists of 60,000 small natural colored images that are 32x32
pixels in size. Each dataset is split into a training set containing 50,000
images and a test set containing 10,000 images. CIFAR-10 samples images from
10 different classes, or target labels, while CIFAR-100 samples from 100
different classes. Thus, CIFAR-10 contains 5,000 images for each class while
CIFAR-100 is comparatively more difficult containing only 500 images for each
class.

All ablations use the same ResNet-18 trained for 100 epochs using a
standardized optimization configuration [13]. We use PyTorch’s Stochastic
Gradient Descent optimizer with an initial learning rate of $0.1$ and Nesterov
momentum of $0.9$. After 50 epochs, the learning rate is decayed to $0.01$ and
again to $0.001$ for the final 10 epochs. We use standard data augmentations
including random crop, random horizontal flip, and mean standard
normalization.

Pruning is done in a layerwise unstructured fashion, after which each
subnetwork is tuned for an additional 20 epochs. We experiment with both a
constant learning rate of $0.01$ and a one-cycle policy with a max learning
rate of $0.1$. We explore the results for each configuration with different
levels of final subnetwork sparsity $\rho\in[0.5,0.7,0.9,0.95,0.98]$. We
additionally include results for subnetworks created through L1 unstructured
pruning. In this case, parameters with the smallest magnitudes at each layer
are pruned.

The one-shot baseline prunes the parent network to the target sparsity before
we start tuning. Iterative pruning calculates the number of parameters to
remove at the beginning of each epoch over some number of pruning epochs
$\varphi\in[5,10,15,20]$ such that the final target sparsity is hit. Random
Annealing uses a probability matrix that is generated according to a random
uniform distribution, where each parameter is assigned a value between 0 and
1. Parameters with a value less than the target sparsity value are linearly
annealed to 0 and values greater than the target sparsity are linearly
annealed to 1, over some number of annealing epochs $\varphi\in[5,10,15,20]$.
Temperature Annealing uses a randomly generated binary bitmask according to
the target sparsity. All parameters with a value of 0 are modified to an
initial temperature value of $\tau\in[0.2,0.4,0.6,0.8,1.0]$. Those values are
then cosine annealed to 0 over some number of annealing epochs
$\varphi\in[5,10,15,20]$.

Table 1 includes the mean accuracies for the best configurations of each
method on CIFAR-100. We see consistent improvement with both of our stochastic
annealing methods over the baseline one-shot and iterative pruning techniques
across all sparsities and with both a constant learning rate and a one-cycle
policy. As networks become more sparse, the benefits from our annealing
approaches become more significant, with a 6% and 4% improvement at 98%
sparsity over the one-shot baseline with a constant and one-cycle rate
schedule respectively. We also observed improved performance with magnitude
pruning, however the differences were smaller at 2% and 1% improvement
respectively.

Figure 1 includes an exploration of the initial temperature, the number of
annealing epochs and the test trajectories for the best models tuned with a
constant learning rate schedule. The hyperparameter with the most significant
impact on all pruning methodologies was the number of epochs that were pruned
or annealed over. We saw best results for all subnetwork sparsities with $e=5$
or $e=10$ annealing epochs. It’s important for the final subnetwork topology
to be established for a sufficient number of epochs to allow for optimal
convergence behavior. This pattern holds for both constant and one-cycle
learning rate schedules. The initial temperature has a smaller impact on
performance than the number of annealing epochs. A small value of $\tau$ means
that the target subnetwork is always active and other parameters have a small
chance to turn on, while a high value of $\tau$ gives other parameters a high
chance to turn on. We saw best results when $\tau$ was in the $0.4$ to $0.6$
range. This corresponds to a higher state of entropy regarding network
structure which results in a stronger regularization effect for those initial
training examples.

Figure 2: Ablations over the number of annealing epochs for random annealing
vs temperature annealing and temperature annealing vs iterative pruning. The
rightmost graph displays the best performing annealing configuration for each
sparsity and plots it against the baseline. Stochastic annealing outperformed
the other methods for each sparsity when used with a one-cycle policy.

Figure 2 includes an exploration of random annealing vs temperature annealing,
iterative pruning vs temperature annealing, and the test trajectories for the
best models with a one-cycle learning rate schedule. Despite the additional
variance associated with random annealing, we found that the performance was
very good and nearly approached that of temperature annealing. Temperature
annealing consistently outperformed iterative pruning. The graph displays
results for annealing epochs $e\in[5,10]$ and a target subnetwork sparsity of
$0.9$. While the early accuracy of temperature annealing appears worse, it
quickly surpasses iterative pruning as the annealing phase ends. The One-Cycle
schedule does reduce the accuracy gap between our method and the baselines due
to the much improved generalization compared to the constant learning rate
models. However, temperature annealing still consistently outperformed the
baseline across all sparsities.

###  IV-B Vision Transformers

Table II: Comparison between one-shot and temperature annealing with Vision Transformers on Tiny Imagenet. We report the best accuracy for each method at various levels of target sparsity when tuned with both SGD and Adam using a constant learning rate and a one-cycle learning rate schedule. Temperature annealing consistently outperforms one-shot pruning methods for various levels of sparsity. ViT/Tiny-ImageNet |  Method | 50% | 60% | 70% | 80% | 90%
---|---|---|---|---|---|---
SGD (Constant) |  One-Shot Baseline | $38.6\pm 0.25$ | $37.8\pm 0.31$ | $29.3\pm 0.38$ | $7.1\pm 0.68$ | $3.0\pm 0.01$
|  Temperature Annealing | $40.6\pm 0.48$ | $38.3\pm 0.19$ | $30.7\pm 0.26$ | $7.1\pm 0.24$ | $3.2\pm 0.09$
SGD (One-Cycle) |  One-Shot Baseline | $42.5\pm 0.11$ | $39.6\pm 0.21$ | $33.4\pm 0.05$ | $11.4\pm 1.46$ | $5.1\pm 0.04$
|  Temperature Annealing | $42.8\pm 0.10$ | $39.8\pm 0.08$ | $33.9\pm 0.31$ | $11.7\pm 0.61$ | $6.4\pm 0.13$
Adam (Constant) |  One-Shot Baseline | $42.0\pm 0.06$ | $42.4\pm 0.08$ | $42.3\pm 0.05$ | $40.0\pm 0.23$ | $35.6\pm 0.25$
|  Temperature Annealing | $42.6\pm 0.20$ | $43.1\pm 0.13$ | $42.3\pm 0.27$ | $40.6\pm 0.12$ | $36.4\pm 0.16$
Adam (One-Cycle) |  One-Shot Baseline | $47.0\pm 0.05$ | $46.8\pm 0.03$ | $45.7\pm 0.18$ | $41.9\pm 0.02$ | $36.9\pm 0.21$
|  Temperature Annealing | $47.6\pm 0.20$ | $47.3\pm 0.20$ | $46.4\pm 0.12$ | $42.0\pm 0.10$ | $37.6\pm 0.19$

We also include a comparison of temperature annealing and one-shot pruning
with modern vision transformer networks on the Tiny Imagenet benchmark [20,
3]. Borrowing from their natural language counterparts, Vision Transformers
make use of attention mechanisms and multilayer percpetrons. These models are
quickly growing in popularity as they continue to produce state-of-the-art
results on large scale and difficult image classification tasks [46]. For our
ablation experiments, we use an open source implementation of the Vision
Transformer with hyperparameters selected to be appropriate for small scale
datasets [34]. These include twelve transformer encoder layers and twelve
attention heads where each encoder contains a self attention block with a
hidden size of 384 followed by a two layer MLP with an expansion factor of 4
for a hidden size of 1536.

The parent models are trained for 200 epochs using Adam with a learning rate
of 1e-4, weight decay of 5e-5, and a batch size of 128. We incorporate a
warmup phase for 5 epochs where the learning rate is linearly increased from 0
to 1e-4 which is known to reduce divergent variance between runs [29].
Standard data augmentation schemes are used including random crop, random
horizontal flip, and autoaugment [2].

We then implement the one-shot pruning and temperature annealing approaches
for target subnetworks of various levels of sparsity
$s\in[0.5,0.6,0.7,0.8,0.9]$. Subnetworks are trained for 20 additional epochs
with both the stochastic gradient descent and adam optimizers [16]. We
additionally explore the efficacy of the one-cycle learning rate schedule with
both of these optimizers. With SGD, the constant configuration uses a learning
rate of 1e-2 and the one-cycle learning rate uses a maximum value of 1e-1,
both with a momentum value of 0.9. With Adam, the constant configuration uses
a learning rate of 1e-3 and the one-cycle learning rate uses a maximum value
of 1e-2.

Table II reports the mean results for each technique over three runs where we
observe small yet consistent improvement with temperature annealing over one-
shot pruning. This holds for both SGD and Adam with both constant and one-
cycle learning rate schedules. There is a stark difference between SGD and
Adam as the target networks become more sparse as Vision Transformers appear
to be much more sensitive to pruning over convolutional networks.

##  V Stochastic Prune and Tune Ensembles

Input: N, the number of ensemble members

Input: ($X_{train},X_{test}$), the training and test data

Input: ($t_{parent},t_{child},t_{anneal}$), the number of parent, child, and
annealing epochs

Input: ($\eta_{parent},\eta_{child}$), the parent and child optimizer
hyperparameters

Input: $\phi(\rho)$, the subnetwork sampling process with target sparsity
$\rho$

Input: $\Theta(M,\tau,e,t)$, the temperature scaling function applied to a
mask $M$, with initial temperature $\tau$, for current epoch $e$, over $t$
total annealing epochs

$O=optimizer(\eta_{parent})$

$F=initialize()$

$F.w=train(F,O,X_{train},t_{parent})$

$ensemble=[]$

for  _i in 1 to $N$_ do

$f_{i}=initialize()$

$f_{i}.w=F.w$

for  _$j$ in $f_{i}.layers$_ do

$M_{j}\sim\phi(\rho)^{f_{i,j}.w}$

end for

$O=optimizer(\eta_{child})$

for  _$e$ in $t_{child}$_ do

$P\leftarrow\Theta(M,\tau,e,t_{anneal})$

for  _$(x,y)$ in $X_{train}$_ do

$\hat{f}_{i}.w\leftarrow f_{i}.w\circ Bernoulli(P)$

$f_{i}.w=train(\hat{f}_{i},O,(x,y))$

end for

end for

$ensemble\leftarrow f_{i}$

end for

for  _$(x,y)$ in $X_{test}$_ do

$outputs=[F(x),\ c(x)\text{ for c in ensemble}$]

$predictions=softmax(mean(outputs))$

end for

Algorithm 1 Stochastic Prune and Tune Ensembles

Ensemble learning is a powerful technique for improving the generalization of
machine learning systems [11, 19]. These algorithms train multiple models
which are then evaluated independently on test data. The combination of
several predictions allows for bias and variance reduction that results in
reliable performance improvement. However, as datasets and neural networks
have grown larger, traditional ensemble methods have become prohibitively
expensive to implement. Recent research has shown that low-cost ensemble
methods can achieve performance that rivals full size ensembles at a
significantly reduced cost. Several of the most powerful low-cost methods do
this by leveraging sparse subnetworks within large parent networks. [42, 38,
25, 12].

Prune and Tune Ensembling (PAT) is one such technique that demonstrates
excellent efficiency for the training compute/accuracy tradeoff. These work by
first training a single parent network for the majority of the training
budget. Child networks are then spawned by cloning and dramatically pruning
the parent using random or anti-random sampling strategies. Each of the child
networks are then fine tuned with a cyclic learning rate schedule for a small
number of epochs.

As child networks are all derived from an identical parent network, anti-
random pruning and one-cycle tuning are used to encourage diversity and reduce
correlation among ensemble members by ensuring that topological structures are
distant and that parameters move far apart in optimization space before
converging [31, 43, 36, 42].

Anti-random pruning creates mirrored pairs of child networks, such that
whenever we randomly prune the parent to create a child, a sibling is created
that inherits the opposite set of parameters. Consider a binary bit string
$M=\\{x_{0},...,x_{n}:x\in\\{0,1\\}\\}$, that is randomly generated where $1$
represents parameters that are kept and $0$ represents parameters that are
pruned. The anti-random network is created by reversing the polarity of all
the bits in the mask $M$, such that:

$\displaystyle\hat{\theta}_{1}$ $\displaystyle=\theta\circ M$ (5)
$\displaystyle\hat{\theta}_{2}$ $\displaystyle=\theta\circ(1-M)$ (6)

where $\hat{\theta}_{i}$ are the parameters of the child network, $\theta$ are
the parameters of the parent network and $\circ$ denotes the Hadamard product.

Stochastic Subnetwork Annealing can be naturally implemented to tune child
networks within the context of this low-cost ensemble algorithm. Random
annealing also extends to the ideas of anti-random pruning through anti-
probability matrices. When a probability matrix is generated, a mirrored
probability matrix can be generated such that the subnetworks anneal to
opposite topological structures $P^{\prime}=1-P$.

Table III: Results for ensembles of WideResNet-28-10 models on both CIFAR-10
and CIFAR-100. Methods with ∗ denote results obtained from [12, 26]. Best low-
cost ensemble results are bold. cAcc, cNLL, and cECE correspond to corrupted
test sets. M corresponds to the number of models in the ensemble. We report
the mean values over 10 runs for Prune and Tune and Temperature Annealing.

Methods (CIFAR-10/WRN-28-10) |  Acc $\uparrow$ |  NLL $\downarrow$ |  ECE $\downarrow$ |  cAcc $\uparrow$ |  cNLL $\downarrow$ |  cECE $\downarrow$ |  Epochs $\downarrow$
---|---|---|---|---|---|---|---
Independent Model∗ | 96.0 | 0.159 | 0.023 | 76.1 | 1.050 | 0.153 | 200
Monte Carlo Dropout∗ | 95.9 | 0.160 | 0.024 | 68.8 | 1.270 | 0.166 | 200
Snapshot (M=5) | 96.3 | 0.131 | 0.015 | 76.0 | 1.060 | 0.121 | 200
Fast Geometric (M=12) | 96.3 | 0.126 | 0.015 | 75.4 | 1.157 | 0.122 | 200
Prune and Tune (M=6) | 96.5 | 0.113 | 0.005 | 76.2 | 0.972 | 0.081 | 200
Temperature Annealing (M=6, $\tau$=0.5, e=3) | 96.7 | 0.110 | 0.005 | 76.3 | 0.968 | 0.079 | 200
TreeNet (M=3)∗ | 95.9 | 0.258 | 0.018 | 75.5 | 0.969 | 0.137 | 250
BatchEnsemble (M=4)∗ | 96.2 | 0.143 | 0.021 | 77.5 | 1.020 | 0.129 | 250
Multi-Input Multi-Output (M=3)∗ | 96.4 | 0.123 | 0.010 | 76.6 | 0.927 | 0.112 | 250
FreeTickets (EDST) (M=7)∗ | 96.4 | 0.127 | 0.012 | 76.7 | 0.880 | 0.100 | 850
FreeTickets (DST) (M=3)∗ | 96.4 | 0.124 | 0.011 | 77.6 | 0.840 | 0.090 | 750
Dense Ensemble (M=4)∗ | 96.6 | 0.114 | 0.010 | 77.9 | 0.810 | 0.087 | 800

Methods (CIFAR-100/WRN-28-10) |  Acc $\uparrow$ |  NLL $\downarrow$ |  ECE $\downarrow$ |  cAcc $\uparrow$ |  cNLL $\downarrow$ |  cECE $\downarrow$ |  Epochs $\downarrow$
---|---|---|---|---|---|---|---
Independent Model∗ | 79.8 | 0.875 | 0.086 | 51.4 | 2.700 | 0.239 | 200
Monte Carlo Dropout∗ | 79.6 | 0.830 | 0.050 | 42.6 | 2.900 | 0.202 | 200
Snapshot (M=5) | 82.1 | 0.661 | 0.040 | 52.2 | 2.595 | 0.145 | 200
Fast Geometric (M=12) | 82.3 | 0.653 | 0.038 | 51.7 | 2.638 | 0.137 | 200
Prune and Tune (M=6) | 82.7 | 0.634 | 0.013 | 52.7 | 2.487 | 0.131 | 200
Temperature Annealing (M=6, $\tau$=0.5, e=3) | 83.1 | 0.633 | 0.010 | 52.8 | 2.440 | 0.131 | 200
TreeNet (M=3)∗ | 80.8 | 0.777 | 0.047 | 53.5 | 2.295 | 0.176 | 250
BatchEnsemble (M=4)∗ | 81.5 | 0.740 | 0.056 | 54.1 | 2.490 | 0.191 | 250
Multi-Input Multi-Output (M=3)∗ | 82.0 | 0.690 | 0.022 | 53.7 | 2.284 | 0.129 | 250
FreeTickets (EDST) (M=7)∗ | 82.6 | 0.653 | 0.036 | 52.7 | 2.410 | 0.170 | 850
FreeTickets (DST) (M=3)∗ | 82.8 | 0.633 | 0.026 | 54.3 | 2.280 | 0.140 | 750
Dense Ensemble (M=4)∗ | 82.7 | 0.666 | 0.021 | 54.1 | 2.270 | 0.138 | 800

We evaluate the efficacy of Stochastic Subnetwork Annealing in this context
via a benchmark low-cost ensemble comparison on CIFAR-10 and CIFAR-100 with
WideResnet-28x10. We include results for these ensembles on corrupted versions
of the CIFAR datasets as well, in order to test robustness on out-of-
distribution images [33]. These corrupted datasets are generated by adding 20
different kinds of image corruptions (gaussian noise, snow, blur, pixelation,
etc.) at five different levels of severity to the original test sets. The
total number of images in each of these additional sets is 1,000,000.

We take the training configuration, ensemble size and parameter settings
directly from studies of three state-of-the-art benchmark low-cost ensemble
methods: MotherNets [39], Snapshot Ensembles [15], and Fast Geometric
Ensembles [7]. We also compare our results with published results of several
recent low-cost ensemble methods including: TreeNets [22], BatchEnsemble [41],
FreeTickets [25], and MIMO [12].

All methods compared use the same WideResnet model and are optimized using
Stochastic Gradient Descent with Nesterov momentum $\mu=0.9$ and weight decay
$\gamma=0.0005$. The ensemble size and training schedule is as used in
previous comparisons [39, 7]. We use a batch size of 128 for training and use
random crop, random horizontal flip, and mean standard normalization data
augmentations for all approaches [7, 12, 26, 15]. The parent learning rate
uses a step-wise decay schedule. An initial learning rate of $\eta_{1}=0.1$ is
used for 50% of the training budget which decays linearly to $\eta_{2}=0.001$
at 90% of the training budget. The learning rate is kept constant at
$\eta_{2}=0.001$ for the final 10% of training.

We train a single parent network for 140 epochs. Six children are then created
by randomly pruning 50% of the connections in the parent network. Neural
partitioning is implemented where pairs of children are generated with
opposite sets of inherited parameters. Each child is tuned with a one-cycle
learning rate for 10 epochs. The tuning schedule starts at $\eta_{1}=0.001$,
increases to $\eta_{2}=0.1$ at 1 epoch and then decays to $\eta_{3}=1e-7$
using cosine annealing for the final 9 epochs.

We implement Stochastic Annealing when tuning the child networks by using the
procedure illustrated above for Temperature Annealing. We initialize the
binary subnetwork masks with a target sparsity of $0.5$. We then modify the
masks such that the 0 parameters are initialized with a temperature of
$\tau=0.5$. That value is decayed to $\tau=0$ over 3 annealing epochs using a
cosine decay.

The Independent Model is a single WideResNet-28-10 model trained for 200
Epochs. The Dropout Model includes dropout layers between convolutional layers
in the residual blocks at a rate of 30% [45]. Snapshot Ensembles use a cosine
annealing learning rate with an initial learning rate $\eta=0.1$ for a cycle
length of 40 epochs [15]. Fast Geometric Ensembles use a pre-training routine
for 156 epochs. A curve finding algorithm then runs for 22 epochs with a cycle
length of 4, each starting from checkpoints at epoch 120 and 156. TreeNets,
[22], BatchEnsemble [41] and Multiple-Input Multiple-Output [12] are all
trained for 250 epochs. FreeTickets introduces several configurations for
building ensembles. We include their two best configurations for Dynamic
Sparse Training (DST, M=3, S=0.8) and Efficient Dynamic Sparse Training (EDST,
M=7, S=0.9).

We compare the results of our approach to a wide variety of competitive
benchmarks which include both low-cost and full-size ensemble approaches. The
details for these implementations are taken from baseline results reported in
[42, 12, 25], and is informed by each original implementation in [15, 7, 22,
41, 12]. We observe improved performance with temperature annealing over these
established methods on both CIFAR-10 and CIFAR-100.

##  VI Conclusions

Stochastic Subnetwork Annealing offers a novel approach to tuning pruned
models by representing subnetworks with probabilistic masks. Rather than
discretely removing parameters, we instead create probability matrices that
alter the chance for parameters to be retained on any given forward pass. The
probability values are then annealed towards a deterministic binary mask over
several epochs such that the subnetwork is slowly revealed. We introduce
several variations of Subnetwork Annealing, including Random Annealing, which
uses random probabilities for every parameter, and Temperature Annealing,
which applies an even amount of stochasticity governed by a hyperparameter
value $\tau$.

The efficacy of Stochastic Subnetwork Annealing is built upon the principles
behind iterative pruning. Instead of pruning a large number of parameters at
once, which can lead to a destructive effect on convergence quality,
Stochastic Annealing slowly evolves subnetworks over time. Recent insights in
optimization research have revealed that early epochs are critical during
optimization as they set the foundational trajectory for the rest of the
training process. [8, 30, 9]. Stochastic Subnetwork Annealing provides
effective regularization during the early epochs of subnetwork tuning to
promote training stability and encourage robust adaptation.

Stochastic Subnetwork Annealing is flexible with a variety of hyperparameter
choices to fit a wide range of learning tasks, network architectures, and
optimization goals. Annealing towards a deteministic binary structure allows
for probability matrices to be tossed out after tuning, resulting in small and
sparse networks that are able to take advantage of hardware optimization
during inference. Our extensive ablation studies explore the dynamics of
Stochastic Annealing with regard to different network architectures, benchmark
datasets, annealing configurations, pruning methodologies, and optimization
hyperparameters.

We observe marked improvement over established one-shot and iterative pruning
benchmarks in a wide variety of cases. This technique is especially effective
for ResNets at high levels of sparsity (95-98%), and we see a similar
performance gain with Vision Transformers on the more difficult Tiny Imagenet
dataset. We additionally demonstrate the utility of this method in the context
of a low-cost ensemble learning algorithm called Prune and Tune Ensembles,
where we used Stochastic Annealing to fine tune the child networks to achieve
increased generalization and robustness performance over benchmark methods,
resulting in better accuracy than methods that require 4-5x more training
compute. Stochastic Annealing is an elegant and flexible technique that can be
used to better tune sparse subnetworks obtained from trained parent models.

## References

  * [1] Davis Blalock, Jose Javier Gonzalez Ortiz, Jonathan Frankle, and John Guttag.  What is the state of neural network pruning?, 2020.
  * [2] Ekin D. Cubuk, Barret Zoph, Dandelion Mane, Vijay Vasudevan, and Quoc V. Le.  Autoaugment: Learning augmentation policies from data, 2019.
  * [3] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale, 2021.
  * [4] Sara Elkerdawy, Mostafa Elhoushi, Hong Zhang, and Nilanjan Ray.  Fire together wire together: A dynamic pruning approach with self-supervised mask prediction, 2022.
  * [5] Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks, 2019.
  * [6] Trevor Gale, Erich Elsen, and Sara Hooker.  The state of sparsity in deep neural networks, 2019.
  * [7] Timur Garipov, Pavel Izmailov, Dmitrii Podoprikhin, Dmitry Vetrov, and Andrew Gordon Wilson.  Loss surfaces, mode connectivity, and fast ensembling of dnns, 2018.
  * [8] Justin Gilmer, Behrooz Ghorbani, Ankush Garg, Sneha Kudugunta, Behnam Neyshabur, David Cardoze, George Dahl, Zachary Nado, and Orhan Firat.  A loss curvature perspective on training instability in deep learning, 2021.
  * [9] Akhilesh Gotmare, Nitish Shirish Keskar, Caiming Xiong, and Richard Socher.  A closer look at deep learning heuristics: Learning rate restarts, warmup and distillation, 2018.
  * [10] Song Han, Jeff Pool, John Tran, and William J. Dally.  Learning both weights and connections for efficient neural networks, 2015\.
  * [11] L.K. Hansen and P. Salamon.  Neural network ensembles.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(10):993–1001, 1990.
  * [12] Marton Havasi, Rodolphe Jenatton, Stanislav Fort, Jeremiah Zhe Liu, Jasper Snoek, Balaji Lakshminarayanan, Andrew M. Dai, and Dustin Tran.  Training independent subnetworks for robust prediction, 2021.
  * [13] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In IEEE Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [14] Weizhe Hua, Yuan Zhou, Christopher De Sa, Zhiru Zhang, and G. Edward Suh.  Channel gating neural networks, 2019.
  * [15] Gao Huang, Yixuan Li, Geoff Pleiss, Zhuang Liu, John E. Hopcroft, and Kilian Q. Weinberger.  Snapshot ensembles: Train 1, get m for free, 2017.
  * [16] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [17] Diederik P. Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick, 2015.
  * [18] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical Report, University of Toronto, May 2012.
  * [19] Anders Krogh and Jesper Vedelsby.  Neural network ensembles, cross validation, and active learning.  In G. Tesauro, D. Touretzky, and T. Leen, editors, Advances in Neural Information Processing Systems, volume 7. MIT Press, 1994.
  * [20] Ya Le and Xuan S. Yang.  Tiny imagenet visual recognition challenge.  2015\.
  * [21] Jaeho Lee, Sejun Park, Sangwoo Mo, Sungsoo Ahn, and Jinwoo Shin.  Layer-adaptive sparsity for the magnitude-based pruning, 2021.
  * [22] Stefan Lee, Senthil Purushwalkam, Michael Cogswell, David Crandall, and Dhruv Batra.  Why m heads are better than one: Training a diverse ensemble of deep networks.  arXiv preprint:1511.06314, 2015.
  * [23] Hao Li, Asim Kadav, Igor Durdanovic, Hanan Samet, and Hans Peter Graf.  Pruning filters for efficient convnets, 2017.
  * [24] Ji Lin, Yongming Rao, Jiwen Lu, and Jie Zhou.  Runtime neural pruning.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [25] Shiwei Liu, Tianlong Chen, Zahra Atashgahi, Xiaohan Chen, Ghada Sokar, Elena Mocanu, Mykola Pechenizkiy, Zhangyang Wang, and Decebal Constantin Mocanu.  Deep ensembling with no overhead for either training or testing: The all-round blessings of dynamic sparsity, 2022.
  * [26] Shiwei Liu, Tianlong Chen, Zahra Atashgahi, Xiaohan Chen, Ghada Sokar, Elena Mocanu, Mykola Pechenizkiy, Zhangyang Wang, and Decebal Constantin Mocanu.  Deep ensembling with no overhead for either training or testing: The all-round blessings of dynamic sparsity, 2022.
  * [27] Shiwei Liu, Tianlong Chen, Xiaohan Chen, Li Shen, Decebal Constantin Mocanu, Zhangyang Wang, and Mykola Pechenizkiy.  The unreasonable effectiveness of random pruning: Return of the most naive baseline for sparse training, 2022.
  * [28] Zhuang Liu, Mingjie Sun, Tinghui Zhou, Gao Huang, and Trevor Darrell.  Rethinking the value of network pruning, 2019.
  * [29] Jerry Ma and Denis Yarats.  On the adequacy of untuned warmup for adaptive optimization, 2019.
  * [30] Jerry Ma and Denis Yarats.  On the adequacy of untuned warmup for adaptive optimization, 2021.
  * [31] Yashwant Malaiya.  Antirandom testing: getting the most out of black-box testing.  In Proceedings of Sixth International Symposium on Software Reliability Engineering. ISSRE’95, pages 86 – 95, 1995.
  * [32] Dmitry Molchanov, Arsenii Ashukha, and Dmitry Vetrov.  Variational dropout sparsifies deep neural networks, 2017.
  * [33] Zachary Nado, Neil Band, Mark Collier, Josip Djolonga, Michael Dusenberry, Sebastian Farquhar, Angelos Filos, Marton Havasi, Rodolphe Jenatton, Ghassen Jerfel, Jeremiah Liu, Zelda Mariet, Jeremy Nixon, Shreyas Padhy, Jie Ren, Tim Rudner, Yeming Wen, Florian Wenzel, Kevin Murphy, D. Sculley, Balaji Lakshminarayanan, Jasper Snoek, Yarin Gal, and Dustin Tran.  Uncertainty Baselines: Benchmarks for uncertainty & robustness in deep learning.  arXiv preprint arXiv:2106.04015, 2021.
  * [34] Omiita.  Vit-cifar.  https://github.com/omihub777/ViT-CIFAR, 2022.
  * [35] Alex Renda, Jonathan Frankle, and Michael Carbin.  Comparing rewinding and fine-tuning in neural network pruning, 2020.
  * [36] Leslie N. Smith and Nicholay Topin.  Super-convergence: Very fast training of neural networks using large learning rates, 2018.
  * [37] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(56):1929–1958, 2014.
  * [38] Johannes von Oswald, Seijin Kobayashi, Alexander Meulemans, Christian Henning, Benjamin F. Grewe, and João Sacramento.  Neural networks with late-phase weights, 2022.
  * [39] Abdul Wasay, Brian Hentschel, Yuze Liao, Sanyuan Chen, and Stratos Idreos.  Mothernets: Rapid deep ensemble learning.  arXiv preprint:1809.04270, 2018.
  * [40] Wei Wen, Chunpeng Wu, Yandan Wang, Yiran Chen, and Hai Li.  Learning structured sparsity in deep neural networks, 2016.
  * [41] Yeming Wen, Dustin Tran, and Jimmy Ba.  Batchensemble: An alternative approach to efficient ensemble and lifelong learning, 2020.
  * [42] Tim Whitaker and Darrell Whitley.  Prune and tune ensembles: Low-cost ensemble learning with sparse independent subnetworks, 2022.
  * [43] Shen Wu, Jandhyala Sridhar, Yashwant Malaiya, and Anura Jayasumana.  Antirandom testing: A distance-based approach.  Journal of VLSI Design, 2008.
  * [44] Kaichao You, Mingsheng Long, Jianmin Wang, and Michael I. Jordan.  How does learning rate decay help modern neural networks?, 2019.
  * [45] Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks, 2017.
  * [46] Xiaohua Zhai, Alexander Kolesnikov, Neil Houlsby, and Lucas Beyer.  Scaling vision transformers, 2022.


# Generating GHZ state in $2m$-qubit spin network

M. A. Jafarizadeha,d , R. Sufiania , F. Eghbalifama, M. Azimia,
S. F. Taghavic and E. Baratie,
aDepartment of Theoretical Physics and Astrophysics, University of Tabriz,
Tabriz 51664, Iran.
bCenter of excellence for photonic, University of Tabriz, Tabriz 51664, Iran.
cInstitute for Studies in Theoretical Physics and Mathematics, Tehran
19395-1795, Iran.
dResearch Institute for Fundamental Sciences, Tabriz 51664, Iran.
eInstitute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52,
01-224 Warszawa, Poland<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider a pure $2m$-qubit initial state to evolve under a particular
quantum mechanical spin Hamiltonian, which can be written in terms of the
adjacency matrix of the Johnson network $J(2m,m)$. Then, by using some
techniques such as spectral distribution and stratification associated with
the graphs, employed in [1, 2], a maximally entangled $GHZ$ state is generated
between the antipodes of the network. In fact, an explicit formula is given
for the suitable coupling strengths of the hamiltonian, so that a maximally
entangled state can be generated between antipodes of the network. By using
some known multipartite entanglement measures, the amount of the entanglement
of the final evolved state is calculated, and finally two examples of four
qubit and six qubit states are considered in details.

Keywords: maximal entanglement , GHZ states, Johnson network, Stratification,
Spectral distribution

PACs Index: 01.55.+b, 02.10.Yn

##  1 Introduction

The idea to use quantum spin chains for short distance quantum communication
was put forward by Bose [3]. After the work of Bose, the use of spin chains
[4]-[18] and harmonic chains [19] as quantum wires have been proposed. In the
previous work [1], the so called distance-regular graphs have been considered
as spin networks (in the sense that with each vertex of a distance-regular
graph, a qubit or a spin was associated) and perfect state transfer (PST) of a
single qubit state over antipodes of these networks has been investigated. In
that work, a procedure for finding suitable coupling constants in some
particular spin Hamiltonians has been given so that perfect and optimal
transfer of a quantum state between antipodes of the corresponding networks
can be achieved, respectively.

Entanglement is one of the other important tasks in quantum communication.
Quantum entanglement in spin systems is an extensively-studied field in recent
years [20,21,22,23], in the advent of growing realization that entanglement
can be a resource for quantum information processing. Within this general
field, entanglement of spin$1/2$ degrees of freedom, qubits, has been in focus
for an obvious reason of their paramount importance for quantum computers, not
to mention their well-known applicability in various condensed-matter systems,
optics and other branches of physics. In [24], authors attempted to generate a
Bell state between distant vertices in a permanently coupled spin network
interacting via invariant stratification graphs (ISGs). At the first step,
they established an upper bound over achievable entanglement between the
reference site and the other vertices. Due to this upper bound they found that
creation of a Bell state between the reference site and a vertex is possible
if the stratum of that vertex is a single element, e.g. antipodal ISGs. The
present work focuses on the provision of $GHZ$ state, by using a $2m$-qubit
initial product state. To this end, we will consider the Johnson networks
$J(2m,m)$ (which are distance-regular) as spin networks. Then, we use the
algebraic properties of these networks in order to find suitable coupling
constants in some particular spin Hamiltonians so that $2m$-qubit $GHZ$ state
can be achieved.

The organization of the paper is as follows: In section 2, we review some
preliminary facts about graphs and their adjacency matrices, spectral
distribution associated with them; In particular, some properties of the
networks derived from symmetric group $S_{n}$ called also Johnson networks are
reviewed. Section 3 is devoted to $2m$-qubit $GHZ$ state provision by using
algebraic properties of Johnson network $J(2m,m)$, where a method for finding
suitable coupling constants in particular spin Hamiltonians is given so that
maximal entanglement in final state is possible. The paper is ended with a
brief conclusion and two appendices.

##  2 preliminaries

###  2.1 Graphs and their adjacency matrices

A graph is a pair $\Gamma=(V,E)$, where $V$ is a non-empty set called the
vertex set and $E$ is a subset of $\\{(x,y):x,y\in V,x\neq y\\}$ called the
edge set of the graph. Two vertices $x,y\in V$ are called adjacent if
$(x,y)\in E$, and in that case we write $x\sim y$. For a graph $\Gamma=(V,E)$,
the adjacency matrix $A$ is defined as

$\bigl{(}A)_{\alpha,\beta}\;=\;\cases{1&if $\;\alpha\sim\beta$\cr
0&\mbox{otherwise}\cr}.$ (2-1)

Conversely, for a non-empty set $V$, a graph structure is uniquely determined
by such a matrix indexed by $V$. The degree or valency of a vertex $x\in V$ is
defined by

$\kappa(x)=|\\{y\in V:y\sim x\\}|$ (2-2)

where, $|\cdot|$ denotes the cardinality. The graph is called regular if the
degree of all of the vertices be the same. In this paper, we will assume that
graphs under discussion are regular. A finite sequence
$x_{0},x_{1},...,x_{n}\in V$ is called a walk of length $n$ (or of $n$ steps)
if $x_{i-1}\sim x_{i}$ for all $i=1,2,...,n$. Let $l^{2}(V)$ denote the
Hilbert space of $C$-valued square-summable functions on $V$. With each
$\beta\in V$ we associate a vector $|\beta\rangle$ such that the $\beta$-th
entry of it is $1$ and all of the other entries of it are zero. Then
$\\{|\beta\rangle:\beta\in V\\}$ becomes a complete orthonormal basis of
$l^{2}(V)$. The adjacency matrix is considered as an operator acting in
$l^{2}(V)$ in such a way that

$A|\beta\rangle=\sum_{\alpha\sim\beta}|\alpha\rangle.$ (2-3)

###  2.2 Spectral distribution associated with the graphs

Now, we recall some preliminary facts about spectral techniques used in the
paper, where more details have been given in Refs. [26,27,28,29]

Actually the spectral analysis of operators is an important issue in quantum
mechanics, operator theory and mathematical physics [30,31]. As an example
$\mu(dx)=|\psi(x)|^{2}dx$ ($\mu(dp)=|\widetilde{\psi}(p)|^{2}dp$) is a
spectral distribution which is assigned to the position (momentum) operator
$\hat{X}(\hat{P})$. The mathematical techniques such as Hilbert space of the
stratification and spectral techniques have been employed in [32,33] for
investigating continuous time quantum walk on graphs. Moreover, in general
quasi-distributions are the assigned spectral distributions of two hermitian
non-commuting operators with a prescribed ordering. For example the Wigner
distribution in phase space is the assigned spectral distribution for two non-
commuting operators $\hat{X}$ (shift operator) and $\hat{P}$ (momentum
operator) with Wyle-ordering among them [34, 35]. It is well known that, for
any pair $(A,|\phi_{0}\rangle)$ of a matrix $A$ and a vector
$|\phi_{0}\rangle$, one can assign a measure $\mu$ as follows

$\mu(x)=\langle\phi_{0}|E(x)|\phi_{0}\rangle,$ (2-4)

where $E(x)=\sum_{i}|u_{i}\rangle\langle u_{i}|$ is the operator of projection
onto the eigenspace of $A$ corresponding to eigenvalue $x$, i.e.,

$A=\int xE(x)dx.$ (2-5)

Then, for any polynomial $P(A)$ we have

$P(A)=\int P(x)E(x)dx,$ (2-6)

where for discrete spectrum the above integrals are replaced by summation.
Therefore, using the relations (2-4) and (2-6), the expectation value of
powers of adjacency matrix $A$ over reference vector $|\phi_{0}\rangle$ can be
written as

$\langle\phi_{0}|A^{m}|\phi_{0}\rangle=\int_{R}x^{m}\mu(dx),\;\;\;\;\
m=0,1,2,....$ (2-7)

Obviously, the relation (2-7) implies an isomorphism from the Hilbert space of
the stratification onto the closed linear span of the orthogonal polynomials
with respect to the measure $\mu$.

###  2.3 Underlying networks derived from symmetric group $S_{n}$

Let $\lambda=(\lambda_{1},...,\lambda_{m})$ be a partition of $n$, i.e.,
$\lambda_{1}+...+\lambda_{m}=n$. We consider the subgroup $S_{m}\otimes
S_{n-m}$ of $S_{n}$ with $m\leq[\frac{n}{2}]$. Then we assume the finite set
$M^{\lambda}=\frac{S_{n}}{S_{m}\otimes S_{n-m}}$ with
$|M^{\lambda}|=\frac{n!}{m!(n-m)!}$ as vertex set. In fact, $M^{\lambda}$ is
the set of $(m-1)$-faces of $(n-1)$-simplex (recall that, the graph of an
$(n-1)$-simplex is the complete graph with $n$ vertices denoted by $K_{n}$).
If we denote the vertex $i$ by $m$-tuple $(i_{1},i_{2},...,i_{m})$, then the
adjacency matrices $A_{k}$, $k=0,1,...,m$ are defined as

$\bigl{(}A_{k})_{i,j}\;=\left\\{\begin{array}[]{c}1\quad\mathrm{if}\;\;\
\partial(i,j)=k,\\\ 0\quad\quad\mathrm{otherwise}\quad\quad\quad(i,j\in
M^{\lambda})\\\ \end{array}\right.,\;\;\ k=0,1,...,m.$ (2-8)

where, we mean by $\partial(i,j)$ the number of components that
$i=(i_{1},...,i_{m})$ and $j=(j_{1},...,j_{m})$ are different (this is the
same as Hamming distance which is defined in coding theory). The network with
adjacency matrices defined by (2-8) is known also as the Johnson network
$J(n,m)$ and has $m+1$ strata such that

$\kappa_{0}=1,\;\ \kappa_{l}=\left(\begin{array}[]{c}m\\\
m-l\end{array}\right)\left(\begin{array}[]{c}n-m\\\ l\end{array}\right),\;\;\
l=1,2,...,m.$ (2-9)

One should notice that for the purpose of maximal entanglement provision, we
must have $\kappa_{m}=1$ which is fulfilled if $n=2m$, so we will consider the
network $J(2m,m)$ (hereafter we will take $n=2m$ so that we have
$\kappa_{m}=1$). If we stratify the network $J(2m,m)$ with respect to a given
reference node $|\phi_{0}\rangle=|i_{1},i_{2},...,i_{m}\rangle$, where
$|i_{1},i_{2},...,i_{m}\rangle\equiv|0...0\underbrace{1}_{i_{1}}0...0\underbrace{1}_{i_{2}}0...0\underbrace{1}_{i_{m}}0\rangle$
and $i_{1}\neq i_{2}\neq...\neq i_{m}$. The unit vectors $|\phi_{i}\rangle$,
$i=1,...,m$ are defined as

$|\phi_{1}\rangle=\frac{1}{\sqrt{\kappa_{1}}}(\sum_{i^{\prime}_{1}\neq
i_{1}}|i^{\prime}_{1},i_{2},...,i_{m}\rangle+\sum_{i^{\prime}_{2}\neq
i_{2}}|i_{1},i^{\prime}_{2},i_{3},...,i_{m}\rangle+...+\sum_{i^{\prime}_{m}\neq
i_{m}}|i_{1},...,i_{m-1},i^{\prime}_{m}\rangle),$
$|\phi_{2}\rangle=\frac{1}{\sqrt{\kappa_{2}}}\sum_{k\neq
l=1}^{m}\sum_{i^{\prime}_{l}\neq i_{l},i^{\prime}_{k}\neq
i_{k}}|i_{1},...i_{l-1},i^{\prime}_{l},i_{l+1},...,i_{k-1},i^{\prime}_{k},i_{k+1}...,i_{m}\rangle,$
$\vdots$ $\hskip
14.22636pt|\phi_{j}\rangle=\frac{1}{\sqrt{\kappa_{j}}}\sum_{k_{1}\neq
k_{2}\neq...\neq k_{j}=1}^{m}\sum_{i^{\prime}_{k_{1}}\neq
i_{k_{1}},...,i^{\prime}_{k_{j}}\neq
i_{k_{j}}}|i_{1},...,i_{k_{1}-1},i^{\prime}_{k_{1}},i_{k_{1}+1},...,i_{k-1},i^{\prime}_{k_{j}},i_{k_{j}+1}...,i_{m}\rangle,$
$\vdots$ $|\phi_{m}\rangle=\frac{1}{\sqrt{\kappa_{m}}}\sum_{i^{\prime}_{1}\neq
i_{1},...,i^{\prime}_{m}\neq
i_{m}}|i^{\prime}_{1},i^{\prime}_{2},...,i^{\prime}_{m}\rangle.$ (2-10)

Since the network $J(2m,m)$ is distance-regular, the above stratification is
independent of the choice of reference node. The intersection array of the
network is given by

$b_{l}=(m-l)^{2}\;\ ;\;\;\ c_{l}=l^{2}.$ (2-11)

Then, by using the Eq. (B-49), the QD parameters $\alpha_{i}$ and $\omega_{i}$
are obtained as follows

$\alpha_{l}=2l(m-l),\;\ l=0,1,...,m;\;\;\ \omega_{l}=l^{2}(m-l+1)^{2},\;\
l=1,2,...,m.$ (2-12)

Then, one can show that [27]

$A|\phi_{l}\rangle=(l+1)(m-l)|\phi_{l+1}\rangle+2l(m-l)|\phi_{l}\rangle+l(m-l+1)|\phi_{l-1}\rangle.$
(2-13)

##  3 $GHZ$ state generation by using quantum mechanical Hamiltonian in the
network $J(2m,m)$

The model we consider is the distance-regular Johnson network $J(2m,m)$
consisting of $N=C^{2m}_{m}=\frac{(2m)!}{m!m!}$ sites labeled by
$\\{1,2,...,N\\}$ and diameter $m$. Then, we stratify the network with respect
to a chosen reference site, say 1 (the discussion about stratification has
been given in appendix A; In these particular networks, the first and the last
strata possess only one node, i.e., $|\phi_{0}\rangle=|1\rangle$ and
$|\phi_{m}\rangle=|N\rangle$). At time $t=0$, a $2m$-qubit state is prepared
in the first (reference) site of the network. We wish to provide a maximal
quantum entanglement between the state of this site and the state of the
$N$-th site after a well-defined period of time, in which the corresponding
network is evolved under a particular Hamiltonian.

If the network be assumed as a spin network, in which a spin-$1/2$ particle is
attached to each vertex (node) of the network, the Hilbert space associated
with the network is given by ${\mathcal{H}}=(C^{2})^{\otimes 2m}$. The
standard basis for an individual qubit is chosen to be
${|0\rangle=|\downarrow\rangle,|1\rangle=|\uparrow\rangle}$. Then we consider
the Hamiltonian

$H_{s}=\frac{1}{2}\sum\limits_{1\leq i<j\leq 2m}H_{ij}$ (3-14)

where, $H_{ij}=\sigma_{i}\cdot\sigma_{j}$ and $\sigma_{i}$ is a vector with
familiar Pauli matrices $\sigma_{i}^{x},\sigma_{i}^{y}$ and $\sigma_{i}^{z}$.
One can easily see that, the Hamiltoniaan (3.14) commutes with the total $z$
component of the spin, i.e., $[\sigma^{z}_{total},H_{s}]=0$, hence the Hilbert
space ${\mathcal{H}}$ decompose into invariant subspaces, each of which is a
distinct eigenspace of the operator $\sigma^{z}_{total}$. So the total number
of up and down spins are invariant under action of Hamiltonian or time
evolution operator. Now, we recall that the kets
$|i_{1},i_{2},\ldots,i_{2m}\rangle$ with
$i_{1},\ldots,i_{2m}\in\\{\uparrow,\downarrow\\}$ form an orthonormal basis
for Hilbert space ${\mathcal{H}}$. Then, one can easily obtain

$H_{ij}|...\underbrace{\uparrow}_{i}...\underbrace{\uparrow}_{j}...\rangle=|...\underbrace{\uparrow}_{i}...\underbrace{\uparrow}_{j}...\rangle$

and

$H_{ij}|...\underbrace{\uparrow}_{i}...\underbrace{\downarrow}_{j}...\rangle=-|...\underbrace{\uparrow}_{i}...\underbrace{\downarrow}_{j}...\rangle+2|...\underbrace{\downarrow}_{i}...\underbrace{\uparrow}_{j}...\rangle.$
(3-15)

Equation (3.15) implies that the action of $H_{ij}$ on the basis vectors is
equivalent to the action of the operator $2P_{ij}-I$, i.e. we have

$H_{ij}=2P_{ij}-I$ (3-16)

where $P_{ij}$ is the permutation operator acting on sites $i$ and $j$. So

$\frac{1}{2}\sum\limits_{1\leq i<j\leq
2m}\sigma_{i}\cdot\sigma_{j}=\sum\limits_{1\leq i<j\leq
2m}P_{ij}-\frac{1}{2}\left(\begin{array}[]{c}2m\\\ 2\\\ \end{array}\right)I,$
(3-17)

In fact restriction of the operator $\sum_{1\leq i<j\leq 2m}P_{ij}$ on the
$m$-particle subspace (subspace spanned by the states with $m$ spin up) which
has dimension $C^{2m}_{m}$, is written as the adjacency matrix $A$ of the
Johnson network $J(2m,m)$, as

$\sum_{1\leq i<j\leq 2m}P_{ij}=A+m(m-1)I.$ (3-18)

For more details see Ref.[1]. Then we stratify the network with respect to a
chosen reference site, say $|\phi_{0}\rangle$. At time $t=0$, the state is
prepared in the $2m$-qubit state $|\psi(t=0)\rangle=|\underbrace{11\ldots
1}_{m}\rangle|\underbrace{00\ldots 0}_{m}\rangle$. Now, we consider the
dynamics of the system to be governed by the Hamiltonian

$H=\sum_{k=0}^{m}J_{k}P_{k}(1/2\sum_{{}_{1\leq i<j\leq
2m}}{\mathbf{\sigma}}_{i}\cdot{\mathbf{\sigma}}_{j}+\frac{m}{2}I),$ (3-19)

Then, by using (3.17)-(3.20), the Hamiltonian can be written as

$H=\sum_{k=0}^{m}J_{k}P_{k}(A)$ (3-20)

$J_{k}$ is the coupling strength between the reference site $|\phi_{0}\rangle$
and all of the sites belonging to the $k$-th stratum with respect to
$|\phi_{0}\rangle$, and $P_{k}(A)$ are polynomials in terms of adjacency
matrix of the Johnson network. Then, the total system is evolved under unitary
evolution operator $U(t)=e^{-iHt}$ for a fixed time interval, say $t$. The
final state becomes

$|\psi(t)\rangle=\sum_{j=1}^{N}f_{jA}(t)|j\rangle$ (3-21)

where, $N$ is the number of vertices, $|j\rangle$s have $2m$ entries inclusive
$m$ entries equal to 1 and the other entries are 0 and
$|A\rangle=|\underbrace{11\ldots 1}_{m}\underbrace{00\ldots 0}_{m}\rangle$ so
that $f_{jA}(t):=\langle j|e^{-iHt}|A\rangle$.

The evolution with the adjacency matrix $H=A\equiv A_{1}$ for distance-regular
networks (see Appendix B) starting in $|\phi_{0}\rangle$, always remains in
the stratification space. For distance-regular network $J(2m,m)$ for which the
last stratum, i.e., $|\phi_{m}\rangle$ contains only one site, then maximal
entanglement between the starting site $|\phi_{0}\rangle\equiv|A\rangle$ and
the last stratum $|\phi_{m}\rangle$ (the corresponding antipodal node) is
generated, by choosing suitable coupling constants $J_{k}$. In fact, for the
purpose of a maximally entangled $GHZ$ state generation between the first and
the last stratum of the network, we impose the constraints that the amplitudes
$\langle\phi_{i}|e^{-iHt}|\phi_{0}\rangle$ be zero for all $i=1,...,m-1$ and
$\langle\phi_{0}|e^{-iHt}|\phi_{0}\rangle=f$,
$\langle\phi_{m}|e^{-iHt}|\phi_{0}\rangle=f^{\prime}$. Therefore, these
amplitudes must be evaluated. To do so, we use the stratification and spectral
distribution associated with the network $J(2m,m)$ to write

$\langle\phi_{i}|e^{-iHt}|\phi_{0}\rangle=\langle\phi_{i}|e^{-it\sum_{l=0}^{m}J_{l}P_{l}(A)}|\phi_{0}\rangle=\frac{1}{\sqrt{\kappa_{i}}}\langle\phi_{0}|A_{i}e^{-it\sum_{l=0}^{m}J_{l}P_{l}(A)}|\phi_{0}\rangle$

Let the spectral distribution of the graph is
$\mu(x)=\sum_{k=0}^{m}\gamma_{k}\delta(x-x_{k})$ (see Eq. (B-53)). The Johnson
network is a kind of network with a highly regular structure that has a nice
algebraic description; For example, the eigenvalues of this network can be
computed exactly (see for example the notes by Chris Godsil on association
schemes [39] for the details of this calculation). Indeed, the eigenvalues of
the adjacency matrix of the network $J(2m,m)$ (that is $x_{k}$’s in $\mu(x)$)
are given by

$x_{k}=m^{2}-k(2m+1-k),\;\;\ k=0,1,\ldots,m.$ (3-22)

Now, from the fact that for distance-regular graphs we have
$A_{i}=\sqrt{\kappa_{i}}P_{i}(A)$ [27],
$\langle\phi_{i}|e^{-iHt}|\phi_{0}\rangle=0$ implies that

$\sum_{k=0}^{m}\gamma_{k}P_{i}(x_{k})e^{-it\sum_{l=0}^{m}J_{l}P_{l}(x_{k})}=0,\;\;\
i=1,...,m-1$

Denoting $e^{-it\sum_{l=0}^{m}J_{l}P_{l}(x_{k})}$ by $\eta_{k}$, the above
constraints are rewritten as follows

$\sum_{k=0}^{m}P_{i}(x_{k})\eta_{k}\gamma_{k}=0,\;\;\ i=1,...,m-1,$
$\sum_{k=0}^{m}P_{0}(x_{k})\eta_{k}\gamma_{k}=f$
$\sum_{k=0}^{m}P_{m}(x_{k})\eta_{k}\gamma_{k}=f^{\prime}.$ (3-23)

From invertibility of the matrix ${\mathrm{P}}_{ik}=P_{i}(x_{k})$ (see Ref.
[2]) one can rewrite the Eq. (3-23) as

$\left(\begin{array}[]{c}\eta_{0}\gamma_{0}\\\ \eta_{1}\gamma_{1}\\\ \vdots\\\
\eta_{d-1}\gamma_{d-1}\\\
\eta_{d}\gamma_{d}\end{array}\right)=P^{-1}\left(\begin{array}[]{c}f\\\ 0\\\
\vdots\\\ 0\\\ f^{\prime}\end{array}\right).$ (3-24)

The above equation implies that $\eta_{k}\gamma_{k}$ for $k=0,1,...,m$ are the
same as the entries in the first column of the matrix
${\mathrm{P}}^{-1}=WP^{t}$ multiplied with $f$ and the entries in the last
column multiplied with $f^{\prime}$, i.e., the following equations must be
satisfied

$\eta_{k}\gamma_{k}=\gamma_{k}e^{-it\sum_{l=0}^{m}J_{l}P_{l}(x_{k})}={(W{\mathrm{P}}^{t})}_{k0}f+{(W{\mathrm{P}}^{t})}_{km}f^{\prime}\;\
,\;\;\ \mbox{for}\;\ k=0,1,...,m,$ (3-25)

with $W:=diag(\gamma_{0},\gamma_{1},\ldots,\gamma_{m})$. By using the fact
that $\gamma_{k}$ and ${(W{\mathrm{P}}^{t})}_{km}$ are real for
$k=0,1,\ldots,m$, and so we have $\gamma_{k}=|{(W{\mathrm{P}}^{t})}_{km}|$ and
$\gamma_{k}={(W{\mathrm{P}}^{t})}_{k0}$. The Eq. (3.26) can be rewritten as

$\eta_{k}=e^{-it\sum_{l=0}^{m}J_{l}P_{l}(x_{k})}=f+\sigma(k)f^{\prime}$ (3-26)

where $\sigma(k)$ is defined as

$\sigma(k)=\;\cases{-1&for odd $k$\cr 1&\mbox{otherwise}\cr}.$ (3-27)

Assuming $f=|f|e^{i\theta}$ and
$f^{\prime}=|f^{\prime}|e^{i{\theta}^{\prime}}$, it should be considered
${\theta}^{\prime}=\theta\pm\frac{\pi}{2}$ then

$e^{-it\sum_{l=0}^{m}J_{l}P_{l}(x_{k})}=e^{i\theta}(|f|\pm
i\sigma(k)|f^{\prime}|)=e^{i(\theta\pm\arctan{(\frac{\sigma(k)|f^{\prime}|}{|f|})}+2c_{k}\pi)};\;\;\;c_{k}\in\mathcal{Z}$
(3-28)

One should notice that, the Eq. (3.29) can be rewritten as

$(J_{0},J_{1},\ldots,J_{m})=-\frac{1}{t}[\theta+2c_{0}\pi\pm\arctan{(\frac{\sigma(k)|f^{\prime}|}{|f|})},\theta+2c_{1}\pi\pm\arctan{(\frac{\sigma(k)|f^{\prime}|}{|f|})},$
$,\ldots,\theta+2c_{m}\pi\pm\arctan{(\frac{\sigma(k)|f^{\prime}|}{|f|})}](W{\mathrm{P}}^{t})$
(3-29)

or

$J_{k}=-\frac{1}{t}\sum_{j=0}^{m}[\theta+2c_{j}\pi\pm\arctan{(\frac{\sigma(k)|f^{\prime}|}{|f|})}](W{\mathrm{P}}^{t})_{jk}$
(3-30)

where $c_{j}$ for $j=0,1,\ldots,m$ are integers. The result (3-30) gives an
explicit formula for suitable coupling constants so that GHZ state in the
final state can be achieved. The final state is as the form

$|\psi(t)\rangle=f|11\ldots 100\ldots 0\rangle+f^{\prime}|00\ldots 011\ldots
1\rangle$ (3-31)

One attempt to provide a computationally feasible and scalable quantification
of entanglement in multipartite systems was made in Refs. [40,41,42]. For a
pure $n$-qubit state $|\psi\rangle$, the so-called global entanglement is
defined as

$Q(|\psi\rangle)=2(1-\frac{1}{N}\sum\limits_{i=0}^{N-1}Tr[\rho_{i}^{2}])$
(3-32)

where $\rho_{i}$ represents the density matrix of $i$th qubit after tracing
out all other qubits. As seen from this definition, the global entanglement
can be interpreted as the average over the (bipartite) entanglements of each
qubit with the rest of the system. The global entanglement for state in Eq.
(3.32) will be

$Q(|\psi\rangle)=4|f|^{2}|f^{\prime}|^{2}$ (3-33)

Also we introduce a simple multiqubit entanglement quantifier based on the
idea of bipartition and the measure negativity (which is two times the
absolute value of the sum of the negative eigenvalues of the corresponding
partially transposed matrix of a state $\rho$) [43]. For an arbitrary
$N$-qubit state $\rho_{s_{1}s_{2}...s_{N}}$ , a multiqubit entanglement
measure can be formulated as [44]

$\overline{\varrho}=\frac{N}{2}\sum\limits_{1}^{\frac{N}{2}}\varrho_{k|N-k}(\rho_{s_{1}s_{2}...s_{N}})$
(3-34)

where $N$ is assumed even, otherwise $\frac{N}{2}$ should be replaced by
$\frac{N-1}{2}$, and $\varrho_{k|N-k}(\rho_{s_{1}s_{2}...s_{N}})$ is the
entanglement in terms of negativity between two blocks of a bipartition
$k|N-k$ of the state $\rho_{s_{1}s_{2}...s_{N}}$. We can define the following
partition-dependent residual entanglements (PREs)

$\Pi_{q_{1}...q_{m}q_{m+1}...q_{k}|q_{k+1}...q_{n}q_{n+1}...q_{N}}={\varrho}^{2}_{q_{1}...q_{m}q_{m+1}...q_{k}|q_{k+1}...q_{n}q_{n+1}...q_{N}}$
$-\varrho^{2}_{q_{1}...q_{m}|q_{k+1}...q_{n}}-\varrho^{2}_{q_{1}...q_{m}|q_{n+1}...q_{N}}-\varrho^{2}_{q_{m+1}...q_{k}|q_{k+1}...q_{n}}-\varrho^{2}_{q_{m+1}...q_{k}|q_{n+1}...q_{N}}$
(3-35)

and

$\Pi^{\prime}_{q_{1}...q_{k}|q_{k+1}...q_{N}}=\varrho^{2}_{q_{1}...q_{k}|q_{k+1}...q_{N}}-\sum\limits_{i-1}^{k}\sum\limits_{j=k+1}^{N}\varrho^{2}_{q_{i}q_{j}}$
(3-36)

For the state in Eq.(3.32), we have

$\Pi_{q_{1}...q_{m}q_{m+1}...q_{k}|q_{k+1}...q_{n}q_{n+1}...q_{N}}=\Pi^{\prime}_{q_{1}...q_{k}|q_{k+1}...q_{N}}=\varrho^{2}_{q_{1}...q_{k}|q_{k+1}...q_{N}}=4|f|^{2}|f^{\prime}|^{2}$
(3-37)

Another useful entanglement measure was introduced in Refs.[45,46] for
$n$-qubit state $|\psi\rangle=\sum_{i=0}^{2^{n}-1}a_{i}|i\rangle$ with even
$n$, as

$\tau(\psi)=2|{\chi}^{*}(a,n)|$ (3-38)

where

${\chi}^{*}(a,n)=\sum\limits_{i=0}^{2^{n-2}-1}sgn^{*}(n,i)(a_{2i}a_{(2^{n-1}-1)-2i}-a_{2i+1}a_{(2^{n-2}-2)-2i}),$
(3-39) $sgn^{*}(n,i)=\;\cases{{(-1)^{N(i)}}&$0\leq
i\leq{2^{n-3}-1}$\cr{(-1)^{N(i)+n}}&$2^{n-3}\leq i\leq{2^{n-2}-1}$\cr}$ (3-40)

where, $N(i)$ is the number of the occurrences of 1 in the $n$-bit binary
representation of $i$ as $i_{n-1}...i_{1}i_{0}$ ( in binary representation,
$i$ is written as $i=i_{n-1}2^{n-1}+...+i_{1}2^{1}+i_{0}2^{0}$). For the state
Eq.(3.32), one can see that

$\tau(\psi)=2|{\chi}^{*}(a,n)|=2|a_{2^{m}-1}a_{2^{2m}-2^{m}}|=2|ff^{\prime}|.$
(3-41)

In order to achieve maximal entanglement ($GHZ$ state), we should have

$|f|=|f^{\prime}|=\frac{1}{\sqrt{2}}$ (3-42)

Then
$Q(|\psi\rangle)=\Pi_{q_{1}...q_{m}q_{m+1}...q_{k}|q_{k+1}...q_{n}q_{n+1}...q_{N}}=\Pi^{\prime}_{q_{1}...q_{k}|q_{k+1}...q_{N}}=\tau(\psi)=1$.

In the following we consider the four qubit state (the case $m=2$)
$|\psi(t=0)\rangle=|1100\rangle$ and the six qubit state (the case $m=3$)
$|\psi(t=0)\rangle=|111000\rangle$ in details: From Eq. (2-12), for $m=2$, the
QD parameters are given by

$\alpha_{1}=2,\;\ \alpha_{2}=0;\;\;\ \omega_{1}=\omega_{2}=4,$

Then by using the recursion relations (B-48) and (B-51), we obtain

$Q_{2}^{(1)}(x)=x^{2}-2x-4,\;\;\ Q_{3}(x)=x(x-4)(x+2),$

so that the stieltjes function is given by

$G_{\mu}(x)=\frac{Q_{2}^{(1)}(x)}{Q_{3}(x)}=\frac{x^{2}-2x-4}{x(x-4)(x+2)}.$

Then the corresponding spectral distribution is given by

$\mu(x)=\sum_{l=0}^{2}\gamma_{l}\delta(x-x_{l})=\frac{1}{6}\\{3\delta(x)+\delta(x-4)+2\delta(x+2)\\},$

which indicates that

$W=\left(\begin{array}[]{ccc}\gamma_{0}&0&0\\\ 0&\gamma_{1}&0\\\
0&0&\gamma_{2}\\\
\end{array}\right)=\frac{1}{6}\left(\begin{array}[]{ccc}1&0&0\\\ 0&3&0\\\
0&0&2\\\ \end{array}\right).$

In order to obtain the suitable coupling constants, we need also the
eigenvalue matrix $P$ with entries
$P_{ij}=P_{i}(x_{j})=\frac{1}{\sqrt{\omega_{1}\ldots\omega_{i}}}Q_{i}(x_{j})$.
By using the recursion relations (B-48), one can obtain $P_{0}(x)=1,\;\
P_{1}(x)=\frac{x}{2}$ and $P_{2}(x)=\frac{1}{4}(x^{2}-2x-4)$, so that

$P=\left(\begin{array}[]{ccc}1&1&1\\\ 2&0&-1\\\ 1&-1&1\\\ \end{array}\right).$

Then, Eq. (3-30) leads to

$-t(J_{0}+2J_{1}+J_{2})=\theta\pm\frac{\pi}{4}\pm 2c_{0}\pi,$
$-t(J_{0}-J_{2})=\theta\mp\frac{\pi}{4}\pm 2c_{1}\pi,$
$-t(J_{0}-J_{1}+J_{2})=\theta\pm\frac{\pi}{4}\pm 2c_{2}\pi.$

Now, by considering $c_{0}=c_{1}=c_{2}=0$ we obtain

$J_{0}=-\frac{\theta}{t},\;\ J_{1}=0,\;\ J_{2}=\mp\frac{\pi}{4t}.$

Also by considering $c_{0}=0,c_{1}=c_{2}=1$ the coupling constants will be

$J_{0}=-\frac{3\theta\pm 5\pi}{3t},\;\ J_{1}=\pm\frac{2\pi}{3t},\;\
J_{2}=\pm\frac{\pi}{12t}$

and by considering $c_{0}=1,c_{1}=c_{2}=0$

$J_{0}=-\frac{3\theta\pm\pi}{3t},\;\ J_{1}=\mp\frac{2\pi}{3t},\;\
J_{2}=\mp\frac{7\pi}{12t}$

From Eq. (2-12), for $m=3$, the QD parameters are given by

$\alpha_{1}=4,\;\ \alpha_{2}=4,\;\ \alpha_{3}=0;\;\;\
\omega_{1}=\omega_{3}=9,\;\omega_{2}=16$

Then by using the recursion relations (B-48) and (B-51), we obtain

$Q_{3}^{(1)}(x)=x^{3}-8x^{2}-9x+36,\;\;\ Q_{4}(x)=(x^{2}-9)(x-9)(x+1),$

so that the stieltjes function is given by

$G_{\mu}(x)=\frac{Q_{3}^{(1)}(x)}{Q_{4}(x)}=\frac{x^{3}-8x^{2}-9x+36}{(x^{2}-9)(x-9)(x+1)}.$

Then the corresponding spectral distribution is given by

$\mu(x)=\sum_{l=0}^{3}\gamma_{l}\delta(x-x_{l})=\frac{1}{20}\\{\delta(x-9)+5\delta(x-3)+9\delta(x+1)+5\delta(x+3)\\},$

which indicates that

$W=\left(\begin{array}[]{cccc}\gamma_{0}&0&0&0\\\ 0&\gamma_{1}&0&0\\\
0&0&\gamma_{2}&0\\\ 0&0&0&\gamma_{3}\\\
\end{array}\right)=\frac{1}{20}\left(\begin{array}[]{cccc}1&0&0&0\\\
0&5&0&0\\\ 0&0&9&0\\\ 0&0&0&5\\\ \end{array}\right).$

By using the recursion relations (B-48), one can obtain $P_{0}(x)=1,\;\
P_{1}(x)=\frac{x}{3}$, $P_{2}(x)=\frac{1}{12}(x^{2}-4x-9)$ and
$P_{3}(x)=\frac{1}{36}(x^{3}-8x^{2}-9x+36)$, so that

$P=\left(\begin{array}[]{cccc}1&1&1&1\\\ 3&1&-\frac{1}{3}&-1\\\
3&-1&-\frac{1}{3}&1\\\ 1&-1&1&-1\\\ \end{array}\right).$

Then, Eq. (3-30) gives

$-t(J_{0}+3J_{1}+3J_{2}+J_{3})=\theta\pm\frac{\pi}{4}\pm 2c_{0}\pi,$
$-t(J_{0}+J_{1}-J_{2}-J_{3})=\theta\mp\frac{\pi}{4}\mp 2c_{1}\pi,$
$-t(J_{0}-\frac{1}{3}J_{1}-\frac{1}{3}J_{2}+J_{3})=\theta\pm\frac{\pi}{4}\pm
2c_{2}\pi,$ $-t(J_{0}-J_{1}+J_{2}-J_{3})=\theta\mp\frac{\pi}{4}\pm 2c_{3}\pi.$

Again, by considering $c_{0}=c_{1}=c_{2}=c_{3}=0$ we obtain

$J_{0}=-\frac{\theta}{t},\;\ J_{1}=J_{2}=0,\;\ J_{3}=\frac{\mp\pi}{4t}.$

##  4 Conclusion

A $2m$-qubit initial state was prepared to evolve under a particular spin
Hamiltonian, which could be written in terms of the adjacency matrix of the
Johnson graph $J(2m,m)$. By using spectral analysis methods and employing
algebraic structures of the Johnson networks, such as distance-regularity and
stratification, a method for finding a suitable set of coupling constants in
the Hamiltonians associated with the networks was given so that in the final
state, the maximal entanglement of the form $GHZ$ state, could be generated.
In this work we imposed a constraint so that all amplitudes in the final state
were equal to zero except to two amplitudes corresponding to the first and the
final strata (any pair of antinodes of the network), where for $J(2m,m)$ these
strata contain only one vertex, then $GHZ$ state was generated. We hope to
generalize this method to arbitrary Johnson networks $J(n,m)$ and other
various graphs, in order to investigate the entanglement of such systems by
using some multipartite entanglement measures.

## Appendix

##  Appendix A Stratification technique

In this section, we recall the notion of stratification for a given graph
$\Gamma$. To this end, let $\partial(x,y)$ be the length of the shortest walk
connecting $x$ and $y$ for $x\neq y$. By definition $\partial(x,x)=0$ for all
$x\in V$. The graph becomes a metric space with the distance function
$\partial$. Note that $\partial(x,y)=1$ if and only if $x\sim y$. We fix a
vertex $o\in V$ as an origin of the graph, called the reference vertex. Then,
the graph $\Gamma$ is stratified into a disjoint union of strata (with respect
to the reference vertex $o$) as

$V=\bigcup_{i=0}^{\infty}\Gamma_{i}(o),\;\ \Gamma_{i}(o):=\\{\alpha\in
V:\partial(\alpha,o)=i\\}$ (A-43)

Note that $\Gamma_{i}(o)=\emptyset$ may occur for some $i\geq 1$. In that case
we have $\Gamma_{i}(o)=\Gamma_{i+1}(o)=...=\emptyset$. With each stratum
$\Gamma_{i}(o)$ we associate a unit vector in $l^{2}(V)$ defined by

$|\phi_{i}\rangle=\frac{1}{\sqrt{\kappa_{i}}}\sum_{\alpha\in\Gamma_{i}(o)}|\alpha\rangle,$
(A-44)

where, $\kappa_{i}=|\Gamma_{i}(o)|$ is called the $i$-th valency of the graph
($\kappa_{i}:=|\\{\gamma:\partial(o,\gamma)=i\\}|=|\Gamma_{i}(o)|$).

One should notice that, for distance regular graphs, the above stratification
is independent of the choice of reference vertex and the vectors
$|\phi_{i}\rangle,i=0,1,...,d-1$ form an orthonormal basis for the so called
Krylov subspace $K_{d}(|\phi_{0}\rangle,A)$ defined as

$K_{d}(|\phi_{0}\rangle,A)=\mathrm{span}\\{|\phi_{0}\rangle,A|\phi_{0}\rangle,\cdots,A^{d-1}|\phi_{0}\rangle\\}.$
(A-45)

Then it can be shown that [25], the orthonormal basis $|\phi_{i}\rangle$ are
written as

$|\phi_{i}\rangle=P_{i}(A)|\phi_{0}\rangle,$ (A-46)

where $P_{i}=a_{0}+a_{1}A+...+a_{i}A^{i}$ is a polynomial of degree $i$ in
indeterminate $A$ (for more details see for example [25,26]).

##  Appendix B Spectral distribution associated with the graphs

In this section we recall some facts about spectral techniques used in the
paper. From orthonormality of the unit vectors $|\phi_{i}\rangle$ given in
Eq.(A-44) (with $|\phi_{0}\rangle$ as unit vector assigned to the reference
node) we have

$\delta_{ij}=\langle\phi_{i}|\phi_{j}\rangle=\int_{R}P_{i}(x)P_{j}(x)\mu(dx).$
(B-47)

By rescaling $P_{k}$ as $Q_{k}=\sqrt{\omega_{1}\ldots\omega_{k}}P_{k}$, the
spectral distribution $\mu$ under question will be characterized by the
property of orthonormal polynomials $\\{Q_{k}\\}$ defined recurrently by

$Q_{0}(x)=1,\;\;\;\;\;\ Q_{1}(x)=x,$
$xQ_{k}(x)=Q_{k+1}(x)+\alpha_{k}Q_{k}(x)+\omega_{k}Q_{k-1}(x),\;\;\ k\geq 1.$
(B-48)

The parameters $\alpha_{k}$ and $\omega_{k}$ appearing in (B-48) are defined
by

$\alpha_{0}=0,\;\;\ \alpha_{k}=\kappa-b_{k}-c_{k},\;\;\;\;\
\omega_{k}\equiv\beta^{2}_{k}=b_{k-1}c_{k},\;\;\ k=1,...,d,$ (B-49)

where, $\kappa\equiv\kappa_{1}$ is the degree of the networks and $b_{i}$’s
and $c_{i}$’s are the corresponding intersection numbers. Following Refs.
[34], we will refer to the parameters $\alpha_{k}$ and $\omega_{k}$ as $QD$
(Quantum Decomposition) parameters (see Refs. [26,27,28,34] for more details).
If such a spectral distribution is unique, the spectral distribution $\mu$ is
determined by the identity

$G_{\mu}(x)=\int_{R}\frac{\mu(dy)}{x-y}=\frac{1}{x-\alpha_{0}-\frac{\omega_{1}}{x-\alpha_{1}-\frac{\omega_{2}}{x-\alpha_{2}-\frac{\omega_{3}}{x-\alpha_{3}-\cdots}}}}=\frac{Q_{d}^{(1)}(x)}{Q_{d+1}(x)}=\sum_{l=0}^{d}\frac{\gamma_{l}}{x-x_{l}},$
(B-50)

where, $x_{l}$ are the roots of the polynomial $Q_{d+1}(x)$. $G_{\mu}(x)$ is
called the Stieltjes/Hilbert transform of spectral distribution $\mu$ and
polynomials $\\{Q_{k}^{(1)}\\}$ are defined recurrently as

$Q_{0}^{(1)}(x)=1,\;\;\;\;\;\ Q_{1}^{(1)}(x)=x-\alpha_{1},$
$xQ_{k}^{(1)}(x)=Q_{k+1}^{(1)}(x)+\alpha_{k+1}Q_{k}^{(1)}(x)+\omega_{k+1}Q_{k-1}^{(1)}(x),\;\;\
k\geq 1,$ (B-51)

respectively. The coefficients $\gamma_{l}$ appearing in (B-50) are calculated
as

$\gamma_{l}:=\lim_{x\rightarrow x_{l}}[(x-x_{l})G_{\mu}(x)]$ (B-52)

Now let $G_{\mu}(z)$ is known, then the spectral distribution $\mu$ can be
determined in terms of $x_{l},l=1,2,...$ and Gauss quadrature constants
$\gamma_{l},l=1,2,...$ as

$\mu=\sum_{l=0}^{d}\gamma_{l}\delta(x-x_{l})$ (B-53)

(for more details see Refs. [35,36,37,38]).

## References

  * [1] M. A. Jafarizadeh and R. Sufiani, (2008), Phys. Rev. A 77, 022315\.
  * [2] M. A. Jafarizadeh, R. Sufiani, S. F. Taghavi and E. Barati, (2008), J. Phys. A: Math. Theor. 41, 475302.
  * [3] Sougato Bose, Phys. Rev. Lett. 91, 207901 (2003).
  * [4] V. Subrahmanyam, Phys. Rev. A 69, 034304 (2004).
  * [5] M. Christandl, N. Datta, A. Ekert and A. J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
  * [6] M. Christandl, N. Datta, T. C. Dorlas, A. Ekert, A. Kay and A. J. Landahl, (2005), Phys. Rev. A 71, 032312.
  * [7] C. Albanese, M. Christandl, N. Datta and A. Ekert, Phys. Rev. Lett. 93, 230502 (2004).
  * [8] T. J. Osborne and N. Linden, Phys. Rev. A 69, 052315 (2004).
  * [9] F. Verstraete, M. A. Martin-Delgado and J. I. Cirac, Phys. Rev. Lett. 92, 087201 (2004).
  * [10] F. Verstraete, M. Popp and J. I. Cirac, Phys. Rev. Lett. 92, 027901.
  * [11] B.Q. Jin and V. E. Korepin, Phys. Rev. A 69, 062314 (2004).
  * [12] M. H Yung, D. W Leung and S. Bose, Quant. Inf. and Comp. 4, 174 (2004).
  * [13] L. Amico, A. Osterloh, F. Plastina, R. Fazio and G. M. Palma, Phys. Rev. A 69, 022304 (2004).
  * [14] V. Giovannetti and R. Fazio, Phys. Rev. A 71, 032314 (2005).
  * [15] D. Burgarth and S. Bose, (2005), Phys. Rev. A 71, 052315.
  * [16] D. Burgarth and S. Bose, (2005), New J. Phys. 7, 135.
  * [17] A. Chan and R. Hosoya, (2004), J. Algebraic Combinatorics 20, 341-351.
  * [18] A. Chan and C. D. Godsil, (2004), J. Combin. Th. Ser. A 106, 165-191.
  * [19] M.B. Plenio, J. Hartley and J. Eisert, New J. Phys. 6, 36 (2004).
  * [20] C. H. Bennett, G. Brassard, C. Crepeau, R. Jozsa, A. Peres, and W. K. Wooters, Phys. Rev. Lett. 70, 1895(1993).
  * [21] J. I. Cirac, A. K. Ekert, S. F. Huelga, and C. Macchiavello, Phys. Rev. A 59, 4249 (1999).
  * [22] C. H. Bennett and S. J. Wiesner, Phys. Rev. Lett. 69, 2881 (1992).
  * [23] M. Ghojavand, M. A. Jafarizadeh, and S. Rouhani, J. Stat. Mech. (2010) P03001.
  * [24] M. Ghojavand, M. A. Jafarizadeh, and S. Rouhani, J. Stat. Mech. (2010) P12023.
  * [25] M. A. Jafarizadeh and S. Salimi, (2006), J. Phys. A : Math. Gen. 39, 1-29.
  * [26] M. A. Jafarizadeh, R. Sufiani and S. Jafarizadeh, (2007), J. Phys. A: Math. Theor. 40, 4949-4972.
  * [27] M. A. Jafarizadeh, R. Sufiani and S. Jafarizadeh, (2008), Journal of Mathematical Physics 49, 073303.
  * [28] M. A. Jafarizadeh, S. Salimi, (2007), Annals of physics 322, 1005-1033.
  * [29] M. A. Jafarizadeh, R. Sufiani, S. Salimi and S. Jafarizadeh, Eur. Phys. J. B 59, 199-216.
  * [30] M. A. Jafarizadeh, R. Sufiani, (2007), Physica A, 381, 116-142.
  * [31] M. A. Jafarizadeh and R. Sufiani, (2007), International Journal of Quantum Information Vol. 5, No. 4, 575-596.
  * [32] H. Cycon, R. Forese, W. Kirsch and B. Simon Schrodinger operators (Springer-Verlag, 1987).
  * [33] P. D. Hislop and I. M. Sigal, Introduction to spectral theory: With applications to schrodinger operators (1995).
  * [34] Y. S. Kim, Phase space picture of quantum mechanics:group theoretical approach, (Science, 1991).
  * [35] H. W. Lee, Physics. Report, 259, 147 (1995).
  * [36] N. Obata, (2004), Quantum Probabilistic Approach to Spectral Analysis of Star Graphs, Interdisciplinary Information Sciences, Vol. 10, 41-52.
  * [37] A. Hora, and N. Obata, (2003), Fundamental Problems in Quantum Physics, World Scientific, 284.
  * [38] J. A. Shohat, and J. D. Tamarkin, (1943), The Problem of Moments, American Mathematical Society, Providence, RI.
  * [39] C. Godsil, (2005), Association schemes, Combinatorics and Optimization, University of Waterloo.
  * [40] D. A. Meyer, and N. R. Wallach, J. Math. Phys. 43 (2002) 4273.
  * [41] G. K. Brennen, Quantum Inf. Comput. 3 (2003) 619.
  * [42] T. Radtke, S. Fritzsche, Computer Physics Communications. 175 (2006) 145 166.
  * [43] G. Vidal, and R. F. Werner, (2002), Phys. Rev. A. 65 032314.
  * [44] Z. X. Man, Y. J. Xia and N. B. An, New Journal of Physics. 12 (2010) 033020\.
  * [45] D. Li, X. Li, H. Huang and X. Li, Phys. Rev. A. 77 (2008) 056302.
  * [46] D. Li, X. Li, H. Huang and X. Li, J. Math. Phys. 50 (2009) 012104.