Probing hadronic interactions with measurements from the Pierre Auger
Observatory

B. Andrada1 on behalf of the Pierre Auger Collaboration2$\star$

1 Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), Buenos Aires, Argentina

2 Observatorio Pierre Auger, Av. San Martín Norte 304, 5613 Malargüe,
Argentina

(Full author list: https://www.auger.org/archive/authors_2021_07.html)

*<EMAIL_ADDRESS>
50th International Symposium on Multiparticle Dynamics
(ISMD2021)
12-16 July 2021 10.21468/SciPostPhysProc.?

## Abstract

The Pierre Auger Observatory is the largest facility in the world to study
ultra-high-energy cosmic rays. It has a hybrid detection technique that
combines the observation of the longitudinal development of extensive air
showers and the measurement of their particles at the ground. This capability
has opened the possibility to probe hadronic interactions taking place at
energies well beyond those accessible by human-made accelerators. In this
report, we present a selection of the latest results on hadronic interactions
with measurements from the Pierre Auger Observatory. These data span over
three decades in energy, showing the tension between data from the muon
component of air showers and predictions based on the most updated hadronic
interaction models.

##  1 Introduction

When ultra-high-energy cosmic rays reach the upper atmosphere, they interact
with air nuclei at a center of mass energy up to $\sqrt{s}\sim 400$ TeV.
Meanwhile, human-made accelerators are able to produce p-p collisions at
$\sqrt{s}=13$ TeV. Thus, ultra-high-energy cosmic rays provide access to
hadronic interactions that are not accessible with accelerators. However, they
cannot be measured directly with good statistics because of their extremely
low flux. All their properties have to be inferred from the extensive air
showers of secondary particles they produce when colliding with atmospheric
nuclei. These showers have two components. The electromagnetic part is formed
predominantly by the decay of neutral pions and grows through pair production
and bremsstrahlung. The hadronic part, on the other hand, is formed mainly by
chain reactions of secondary mesons and baryons which interact until their
energies are low enough such that the probability for decay becomes greater
than that of further interactions. Muons are one of the products of these
decays.

The Pierre Auger Observatory [1], located near the city of Malargüe, in
Argentina, is the largest facility in the world designed to observe ultra-
high-energy cosmic rays and it combines two independent measurement techniques
providing a hybrid detection mechanism. The Surface Detector (SD) studies the
lateral distribution of air shower particles that reach ground level with more
than 1600 water-Cherenkov detectors (WCD) that have a duty cycle close to
100%. WCDs are distributed in three hexagonal grids with different spacing
between detectors and covering different areas. The SD-1500 covers an area of
about $3000$ km2 in a hexagonal grid of $1500$ m spacing, while the SD-750 and
SD-433 cover areas of $27$ km2 and $1.9$ km2 with hexagonal grids of $750$ m
and $433$ m spacing respectively. Each of these arrays aims at covering
different primary energies: the larger the spacing and area covered, the
higher the observed energy. The Fluorescence Detector (FD) views the
atmosphere above the SD with $27$ telescopes distributed in four sites
surrounding the SD. It measures the nitrogen fluorescence light produced by
the passage of relativistic charged particles through the atmosphere and
operates only on clear and moonless nights, having a $\sim 15\%$ duty cycle. A
major upgrade is taking place in the Auger Observatory [2], AugerPrime, that
will enable the disentanglement of the electromagnetic and muonic components
of the shower, improving the sensitivity to hadronic interactions and mass
composition. This upgrade includes the installation of an Underground Muon
Detector (UMD) to allow for the direct detection of the muon content of
extensive air showers. This detector consists of 30 m2 scintillators located
next to each WCD of the SD-750 and buried 2.3 m underground to shield all
particles of the shower except for muons with energies greater than 1 GeV.

One of the main purposes of the Pierre Auger Observatory is to understand the
mass composition of ultra-high-energy cosmic rays and its astrophysical
implications [3]. An important observable involved in these studies is the
atmospheric depth of maximum development of the shower, $X_{\rm{max}}$, which
is sensitive to mass composition because showers from heavier primaries
develop higher in the atmosphere (lower $X_{\rm{max}}$) and their profiles
fluctuate less than those of showers produced by lighter primaries.
$X_{\rm{max}}$ and its fluctuations can be directly measured with the FD.
Further information relevant to mass composition studies can be obtained by
looking at the muon content of air showers, due to heavier primaries being
expected to produce larger number of muons than lighter primaries. This is
done with Auger data using different observables and techniques, coming both
from SD and UMD measurements.

The measurements of extensive air shower properties, however, need to be
interpreted using the most updated hadronic interaction models which are tuned
to the latest LHC data. When doing so, a deficit in the number of muons
predicted by the models has been observed by a number of different
experiments, using diverse detection techniques and covering a wide energy
range (above $10^{16}$ eV) [4]. The Auger Observatory, being able to study
different hadronically-sensitive observables, has contributed to determine and
characterize the muon deficit in simulations with several results, a selection
of which will be given here.

##  2 Top-Down analysis

A novel _top-down_ analysis was performed on a set of high-quality hybrid
events, with energies between $10^{18.8}$ eV and $10^{19.2}$ eV. This study is
based on the fact that when an attempt is made to reproduce measured events
with simulations, the longitudinal profiles, as seen by the FD, can be well
reproduced but the lateral distribution of particles at the ground, as
measured with the SD, produce systematically smaller signals in the
simulations than the ones found in the data. For each event, batches of
simulations were generated using the measured energy and geometry as inputs
and considering different primaries. In those cases where a matching
longitudinal profile was found, the discrepancies in the SD signals were
scrutinized.

The measured signals in the WCDs derive from both muons and electromagnetic
particles and these contributions cannot be easily decoupled. However, it is
straightforward to obtain both the electromagnetic and hadronic parts of the
signal when using simulations. Two scaling factors were defined to allow for
shifts in each component of the simulated signals: $R_{E}$ rescales the
primary energy affecting both components of the signal and $R_{\rm{had}}$
rescales only the contribution to the ground signal of inherently hadronic
origin, which consists mostly of muons. Given an appropriate set of events,
these two factors can be separately determined because of the different
behaviour of the components with zenith angle, namely that the electromagnetic
part of the shower is more strongly attenuated in the atmosphere than the
muonic and the slant depth increases with zenith angle. A maximum likelihood
fit was performed to find the parameters that minimize the difference between
the rescaled simulated signals and the measured ones. Results are shown in
Figure 1 (left) for two hadronic interaction models and two composition
scenarios: pure proton and a mixed proton-iron composition. For the mixed
composition scenario, no energy rescaling of the primary was necessary, but
the hadronic signal in measurements is between $30\%$ and $60\%$ larger than
the one predicted by hadronic interaction models. The details of this analysis
can be found in [5].

Figure 1: _Left:_ Rescaling factors for the simulated hadronic signal
$R_{\rm{had}}$ and the primary energy $R_{E}$ for two hadronic interaction
models considering alternatively a pure proton composition and a mixed proton-
iron composition. The grey rectangles indicate the systematic uncertainties.
[5]. _Right:_ Energy evolution of the average muon content as measured with
the UMD (white circles) and with the SD (black circles). Predictions for
proton and iron primaries for two hadronic interaction models are also shown.
[6]

##  3 Highly inclined events

When the incoming cosmic ray has a large zenith angle, $\theta>60^{\circ}$,
the electromagnetic component of the corresponding extensive air shower is
strongly attenuated in the atmosphere. By the time the shower front reaches
ground level, the particles that produce a signal in the WCDs are mostly muons
or electrons arising from the decay of muons. Thus, a data set of hybrid, high
energy, highly inclined events allows for the simultaneous and independent
measurement of the two shower components: the electromagnetic part of the
shower is seen by the FD while the SD provides information on the muon
content.

In Figure 1 (right) the average muon content $R_{\mu}$ (black circles) is
shown as a function of the primary energy reconstructed by the FD. Predictions
from hadronic interaction models for proton and iron primaries are also shown.
A higher number of muons is found in data than in simulations. Even though the
data is compatible within the uncertainties with the predictions for iron
primaries, when independent $X_{\rm{max}}$ measurements (directly related with
the primary composition) are contemplated in the analysis, as can be seen in
Figure 2 (left) for the events in the energy bin corresponding to $10^{19}$
eV, the measured number of muons clearly exceeds all model predictions. A
thorough description of this analysis can be found in [7] and an update of it
is available in [8].

##  4 Direct muon measurement

The direct measurement of the muon content with the engineering array of the
UMD provides insight into the muon deficit issue on a lower energy region,
between $2\times 10^{17}$ eV and $2\times 10^{18}$ eV, for the first time with
Auger data. The muon content obtained with UMD data is shown in Figure 1
(right) as open circles in the lower energy region of the plot, together with
predictions from two hadronic interaction models for proton and iron
primaries. Once again, predictions from the models fall short with respect to
Auger data points, even though they are compatible within uncertainties with
heavy primary predictions. However, when independently obtained $X_{\rm{max}}$
information is considered, as was done for the analysis of the highly inclined
events, Auger data is found to lie completely outside the phase space
predicted by the models, as can be seen in Figure 2 (right). So, even for
cosmic rays with a primary energy much closer to the maximum center of mass
energy achieved in the latest LHC run, a deficit in the production of muons in
the hadronic models is found. The complete analysis of UMD data can be found
in [9].

Figure 2: _Left:_ Auger SD measurement and model predictions of the average
muon number and the depth of shower maximum for the events in the energy bin
corresponding to $10^{19}$ eV from Figure 1 (right). [8] _Right:_ Auger UMD
measurement and model predictions of the average muon number and the depth of
shower maximum for the events in the energy bin corresponding to $10^{17.5}$
eV from Figure 1 (right). [9]

##  5 Fluctuations in the number of muons

The latest Auger contribution to the muon deficit issue is an analysis of high
energy, highly inclined events, in which the fluctuations in the number of
muons were considered for the first time. The fluctuations relative to the
average muon number are shown as a function of energy in Figure 3 (left).
Auger data is well contained between the model predictions and, furthermore,
it is compatible with predictions contemplating the mass composition
information based on independent $X_{\rm{max}}$ measurements, shown in the
plot as a grey band. Figure 3 (right) shows the measured relative fluctuations
and the muon number (black marker) for the events in the energy bin
corresponding to $10^{19}$ eV. The colored contours determine the regions
containing the expected values for any combination of the four indicated
primaries (p, He, N and Fe) while the star symbols indicate the expected
values for the mixture of primaries that results from $X_{\rm{max}}$
measurements and the shaded areas represent their uncertainties. This plot
shows simultaneously how the fluctuations in the number of muons are
adequately predicted while the average muon number falls short for all models.
This result favours a small muon deficit present at every stage of the shower
over a large discrepancy only in the first interactions. The details of this
analysis and further discussion of its results can be found in [8].

Figure 3: _Left:_ Relative fluctuations in the number of muons as a function
of energy. Predictions from hadronic interaction models for proton and iron
primaries are also shown. The grey band shows the expected values from mass
composition measurements. _Right:_ Auger measurement and model predictions for
the fluctuations and average muon number for the events in the energy bin
corresponding to $10^{19}$ eV. [8]

##  6 Conclusion

The hybrid detection techniques that the Pierre Auger Observatory uses to
study ultra-high-energy cosmic rays provide insights into the hadronic
interactions taking place in the development of extensive air showers.
Different analysis have shown that the most updated hadronic interaction
models fail to predict the muon content at the ground that is seen with Auger
data, when mass composition information is considered. The fluctuations in the
number of muons have also been studied and shown to be adequately described by
the models. It all seems to indicate that the muon production mechanism of the
models has a small deficit at every stage of shower development, and several
theoretical models are currently under discussion to explain and, eventually,
fix this issue. The upgrade of the Pierre Auger Observatory will provide
further insight into both the primary mass composition and hadronic
interactions of ultra-high-energy cosmic rays.

## References

  * [1] A. Aab _et al._ ,  _The Pierre Auger Cosmic Ray Observatory_ ,  Nucl. Instrum. Meth. A 798, 172 (2015),  10.1016/j.nima.2015.06.058.
  * [2] A. Aab _et al._ ,  _The Pierre Auger Observatory Upgrade - Preliminary Design Report_ ,  http://arxiv.org/abs/1604.03637.
  * [3] A. Yushkov,  _Mass Composition of Cosmic Rays with Energies above $10^{17.2}\,\mathrm{eV}$ from the Hybrid Data of the Pierre Auger Observatory_,  PoS ICRC2019, 482 (2020),  10.22323/1.358.0482.
  * [4] L. Cazon,  _Working Group Report on the Combined Analysis of Muon Density Measurements from Eight Air Shower Experiments_ ,  PoS ICRC2019, 214 (2020),  10.22323/1.358.0214.
  * [5] A. Aab _et al._ ,  _Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory_ ,  Phys. Rev. Lett 177, 192001 (2016),  10.1103/PhysRevLett.177.192001.
  * [6] L. Cazon,  _Measurements and tests of hadronic interactions at ultra-high-energies with the Pierre Auger Observatory_ ,  EPJ Web Conf. 210, 02002 (2019),  10.1051/epjconf/201921002002.
  * [7] A. Aab _et al._ ,  _Muons in Air Showers at the Pierre Auger Observatory: Mean Number in Highly Inclined Events_ ,  Phys. Rev. D 91, 032003 (2015),  10.1103/PhysRevD.91.032003,  [Erratum: Phys.Rev.D 91, 059901 (2015)].
  * [8] A. Aab _et al._ ,  _Measurement of the fluctuations in the number of muons in extensive air showers with the Pierre Auger Observatory_ ,  Phys. Rev. Lett. 126, 152002 (2021),  10.1103/PhysRevLett.126.152002.
  * [9] A. Aab _et al._ ,  _Direct measurement of the muonic content of extensive air showers between $2\times 10^{17}$ and $2\times 10^{18}\,\mathrm{eV}$ at the Pierre Auger Observatory_,  Eur. Phys. J. C 80, 751 (2020),  10.1140/epjc/s10052-020-8055-y.


# Neural Networks with Divisive normalization
for image segmentation with application in cityscapes dataset

Pablo Hernández-Cámara, Valero Laparra and Jesús Malo

Image Processing Lab., Universitat de València
46980 Paterna, Spain
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

One of the key problems in computer vision is adaptation: models are too rigid
to follow the variability of the inputs. The canonical computation that
explains adaptation in sensory neuroscience is _divisive normalization_ , and
it has appealing effects on image manifolds. In this work we show that
including _divisive normalization_ in current deep networks makes them more
invariant to non-informative changes in the images. In particular, we focus on
U-Net architectures for image segmentation. Experiments show that the
inclusion of _divisive normalization_ in the U-Net architecture leads to
better segmentation results with respect to conventional U-Net. The gain
increases steadily when dealing with images acquired in bad weather
conditions. In addition to the results on the Cityscapes and Foggy Cityscapes
datasets, we explain these advantages through visualization of the responses:
the equalization induced by the _divisive normalization_ leads to more
invariant features to local changes in contrast and illumination111An extended
version is under consideration at Patter Recognition Letters..

Keywords— Divisive normalization, Adaptation, Manifold alignment, Image
segmentation, Neural networks.

##  1 Introduction

A fundamental problem in image analysis is the variability of the image
manifold depending on acquisition conditions [1]. Usually sources are not
stationary: the texture and color of an object can be different depending on
the acquisition conditions (see Fig. 1). This has an obvious negative impact
on applications such as segmentation. Adaptation to (and compensation of) non-
informative factors of variation is key for optimal segmentation. However it
is not clear how deep-learning, e.g. the popular U-Nets [2], actually achieve
this.

In this work we propose the use of the canonical computation that accounts for
adaptation in biological neurons, namely _divisive normalization_ [3], to
improve adaptation of artificial networks in an explainable way. We illustrate
this adaptation ability in image segmentation under different atmospheric
conditions.

Figure 1: The problem: changes along the image due to non-informative factors.
In this example, texture and color locally change due to the shadows,
illumination, atmospheric scattering and depth (see highlighted regions).
Associated problems for statistical learning are shown in the scatter plots
(where samples are colored according to their region). The top-right plot
shows changes in the spatial texture defined by the luminance of neighbor
pixels: the contrast (distance from the diagonal) is smaller for lighter
regions (away from the origin), and the energy in the different directions of
the image space (spatial frequency) locally changes along the image.
Atmospheric scattering implies local changes in the hue: see the shift in
dominant wavelengths in the CIExy color diagram.

##  2 Why Using Divisive normalization?

The linear+nonlinear structure of conventional artificial neurons [4] comes
from the seminal computation proposed for sensory neurons [5]:
$\mathbf{x}\overset{\mathcal{L}}{\longrightarrow}\mathbf{z}\overset{\mathcal{N}}{\longrightarrow}\mathbf{y}$,
where the intermediate linear response,
$\mathbf{z}=\mathcal{L}\cdot\mathbf{x}$, is given by the matrix $\mathcal{L}$,
that contains the so called (linear) receptive fields, and original versions
of $\mathcal{N}(\cdot)$ were simple point-wise threshold or saturating
functions [6] eventually rectified. These were the inspiration for current
sigmoids/ReLU in deep-learning [4].

In sensory neuroscience, the _divisive normalization_ model of
$\mathcal{N}(\cdot)$ was a way to account for the inhibitory effect of
neighbor neurons within a layer [3]:

$\displaystyle y_{ip}$ $\displaystyle=$
$\displaystyle\frac{z_{ip}}{\left(\beta_{i}+\sum_{jp^{\prime}}{\gamma_{ijpp^{\prime}}|z_{jp^{\prime}}|^{\alpha_{j}}}\right)^{\epsilon_{i}}}$
(1)

where the linear response of a neuron tuned to the $i$-th feature and to the
$p$-th spatial location, $z_{ip}$, is inhibited (normalized) by a pool of the
activity of neighbor neurons tuned to other features and spatial locations
$z_{jp^{\prime}}$. The constant $\beta_{i}$ determines the level in which the
pool generates effective inhibition, and the exponents control the norm of the
pool. While the interaction kernel in the denominator, $\gamma$, can have
whatever dense structure [7], in visual neuroscience the spatial interaction
is usually assumed to be convolutional in the fovea [8].

However, why using this biological transform to improve artificial networks
devoted to image segmentation? There are examples of the effectiveness of
_divisive normalization_ in many applications as image coding [9, 10, 11],
restoration [12, 13], distortion metrics [14, 15, 16], classification [17, 18,
19, 20] or even referring to its appealing statistical properties [21, 22, 23,
24, 25]. Here we illustrate the effect of _divisive normalization_ with an
explicit example. Fig. 2 considers a patch from the forest (from Fig. 1) that
displays a space varying contrast due to atmospheric effects, and shows how
_divisive normalization_ modifies the features and the image manifold in a
positive way to identify this object class.

The spatial texture problem is enough to stress the point of local
normalization so in this motivating example we restrict ourselves to the
luminance channel222Note that local color compensation could also be done
using modern divisive normalization formulations of classical Von Kries ideas
[26, 27].. In this achromatic context, a reasonable physiological model of
texture perception includes a linear layer of wavelet-like filters (in this
illustration $\mathcal{L}$ is a steerable wavelet transform as in [21, 7]),
and a _divisive normalization_ (in this case with the psychophysically tuned
parameters in [28]).

The neighborhood in $\gamma$ that defines the pooling region (illustrated in
Fig. 2 by the circles in the feature vectors or wavelet bands) has a
qualitatively relevant effect: the division by the local pool moderates the
response in regions where contrast is high, in orange, while boosts the
response in regions of low contrast, in white. See the difference between the
features before and after the local normalization, and note the equalization
of the amplitude of the edges along the field of view. The local normalization
compensates for the contrast variation along the image due to atmospheric
conditions. As a result, while the distribution of the samples from the
original luminance channel display separate clusters corresponding to spatial
variations of the texture, the different samples have been compacted due to
the contrast equalization effect.

In the discussion we will explicitly check if the artificial networks for
segmentation equiped with divisive normalization also develop this contrast
enhancement and equalization effects that have been found in biological
neurons.

Figure 2: The idea for the solution: manifold equalization though Divisive
Normalization. At the top-left we have the luminance input to a
linear+nonlinear model of the V1 cortex based on wavelets and divisive
normalization [21, 7, 28]. Below we have the wavelet-like linear features,
$\mathbf{z}=\mathcal{L}\cdot\mathbf{x}$, and the corresponding divisive
normalized features, $\mathbf{y}=\mathcal{N}(\mathbf{z})$. The scatter plots
display the distributions of image patches taken from $\mathbf{x}$ and
$\mathbf{x^{\prime}}$.

##  3 Models and Experiments

In order to confirm the intuition obtained from the 1-layer, purely
biological, non-optimized model presented in the previous section, here we
implemented the generic _divisive normalization_ layer in Eq. 1 in an
automatic differentiation context. Therefore, one may include our layer at any
place of any architecture and optimize the network for the desired task.

Qualitatively critical points are (a) the nature of the interaction kernel,
$\gamma$, and (b) the specific locations to include the normalization in the
U-Nets for segmentation.

Regarding the nature of $\gamma$, the above example points out that the
consideration of spatial neighborhoods is convenient to get the local contrast
equalization along the visual field required to overcome shadow, scattering or
fog. The recent literature that exploits automatic differentiation shows a
range of kernel structures: some do not consider spatial interactions (either
in a dense [24, 11, 16] or convolutional [20] combinations of features); while
others do with some restrictions (either uniform weights [29], a ring of
locations [18, 19], or special symmetries in the space [30]). Following
biology [8, 21, 3, 28] and the intuition pointed out in the previous section.
In our experiments we use the version in Eq. 1 which includes a general
spatial surround, as in [29], but allowing variable weights over space.
Additionally a general (dense) interaction over channels was considered.
Besides the $\gamma$ parameter, the $\beta$ parameter is also trained, but we
set $\alpha$ and $\epsilon$ to $1$ which we found it makes the training more
stable.

Regarding where to include the proposed normalization in U-Nets, we recall the
_color_ and _texture_ problems pointed out in Fig. 1: following color
appearance models [26, 27], constancy maybe addressed by a single
normalization of the photoreceptors (right before the first convolution),
while texture equalization over space may require normalization at deeper
stages where filters tuned to specific patterns have emerged. This intuition
lead us to the specific scheme in Fig. 3 where _green_ layers stand for the
proposed _divisive normalization_. Note that we propose modifications only in
the _encoding/compressive_ part, where the features that will lead to the
segmentation are computed. The color/texture distinction shapes the
experiments because we will compare the behavior of the improved network (4 DN
layers) with regard to the No-DN network, but we will also consider a 1 DN
layer network eventually addressing only the color problem, but not the
texture problem.

Figure 3: U-Net for segmentation including 4 DN layers (in green). Black
arrows represent the skip unions. The model with only one DN layer had only
the first green layer. The model without DN layers did not have any of the
green layers. The consideration of the DN layers only increases the number of
parameters by an 1.8$\%$.

In the experiments we used the Cityscapes dataset [31], that includes scenes
with annotated segmentation ground truth and the Foggy Cityscapes dataset
[32], that includes degraded versions of the scenes simulating poor weather
conditions of controlled (low, medium, high) severity. In fog, spatial varying
degradation is a major problem. The bottom pixels of each image were cropped
from $1024\times 2048$ to $768\times 2048$ and then normalized and resized to
$96\times 256$ to maintain the aspect ratio. We used 2675, 300 and 500 images
for train, validate and test. We used MAE loss function, a batch-size of 64
and Adam optimiser with an initial learning rate of $0.001$ during 500 epochs
with the original scenes (in regular weather conditions). We kept the models
with higher Intersection over Union (IoU) in validation (in regular weather)
and used these for the test (in poor weather). Testing in more general
datasets where visibility is strongly reduced will confirm or refute the
intuition in Fig. 2.

##  4 Results and Discussion

Dataset | No DN | 1 DN | 4 DN
---|---|---|---
Original | 0.73 | 0.73 (0.0%) | 0.78 (6.8%)
Low fog | 0.65 | 0.63 (-3.1%) | 0.70 (7.7%)
Middle fog | 0.55 | 0.52 (-5.5%) | 0.61 (10.9%)
High fog | 0.39 | 0.42 (7.7%) | 0.46 (17.9%)

IoU reductions (due to fog)
---
Dataset change | No DN | 1 DN | 4 DN
Original-low fog | -11.0% | -14.3% | -9.8%
Original-middle fog | -24.9% | -28.8% | -21.9%
Original-high fog | -46.9% | -42.6% | -40.9%

Table 1: Top: IoU performance in test. Improvements of the use of DN layers
with regard no using DN in parenthesis. Bottom: Reductions in test IoU when
comparing each fog level with the original.

The top panel in Table 1 shows the test results of each model in the four
weather conditions for models trained only using the original (good weather
conditions) images. The bottom panel shows how the test IoU changes in poor
conditions with regard to regular weather.

The first observation is that using 4-DN layers gives substantial increase in
IoU in all cases. Using just 1-DN only helps in really heavy fog. This is
consistent with the previous explanation since introducing fog into the images
mainly changes the textures, which as previously mentioned requires several
normalizations in deep stages. Second, as expected, for progressively heavier
fog, IoU gets reduced in all cases. However the conventional architecture is
more sensitive to the decrease in visibility (bigger reductions in
performance) than the architecture with 4-DN layers. And third, the use of
4-DN always leads to improvements with regard the no-DN case, but it is
important to see that the gains get progressively bigger when the acquisition
conditions are poor.

Figure 4: Segmentation results for the four datasets (original, low, middle
and high fog) and the three models.

Figure 5: Effect of the 2nd-DN layer in two channels, one tuned to vertical
and one tuned to horizontal features: activity before (input) and after
(output) the divisive normalization layer. Results are shown for two different
fog situations.

Figure 4 shows illustrative examples of the predictions, which are consistent
with average IoUs in the tables. Note that only the 4-DN model is able to get
the border in the shadow and preserve the detection of the car in heavy fog.

In order to understand the above advantages we explored the effect of the two
initial DN layers. The first DN, applied to the input RGB signal performs a
Weber-like saturating non-linearity (not shown), very much like a square root
equalization in the separate channels: it boosts low values and moderates high
values. More interestingly, Fig. 5 shows the effect of the second DN layer on
two illustrative features tuned to horizontal and vertical edges (key for
object recognition). We see that in the foggy image the deepest parts of the
image were impossible to distinguish. However, after applying the DN
transformation that the model had learnt with the original (no foggy) images,
the output is much more like the case without fog, now being able to
distinguish details that are further away.

##  5 Conclusion

We presented results in image segmentation when including _divisive
normalization_ in the U-Net architecture. We found that introducing more than
one _divisive normalization_ layer leads us to an improvement in IoU in the
Cistyscapes dataset. Moreover, a model trained only with images acquired in
good weather conditions obtains an increase in IoU over conventional U-Net
when tested with images in bad weather conditions. In fact, the higher the
level of fog introduced in the images, the higher the advantage, which shows
that this transformation helps the models to be more invariant under the non-
informative variations introduced by contrast and luminance changes.

## References

  * [1] Y. Bengio,  “Learning deep architectures for AI,”  Fdns and Trends in ML, vol. 2, no. 1, pp. 1–127, 2009.
  * [2] O. Ronneberger et al.,  “U-net: Convolutional networks for biomedical image segmentation,”  in 2015 MICCAI, 2015, pp. 234–241.
  * [3] M. Carandini et al.,  “Normalization as a canonical neural computation,”  Nature Reviews Neuroscience, vol. 13, pp. 51–62, 2012.
  * [4] S.S. Haykin,  Neural networks and learning machines,  Pearson Education, Upper Saddle River, NJ, third edition, 2009.
  * [5] D. H. Hubel et al.,  “Receptive fields of single neurons in the cat’s striate cortex,”  The Journal of Physiology, vol. 148, no. 3, pp. 574–591, 1959.
  * [6] K. I. Naka et al.,  “S‐potentials from luminosity units in the retina of fish (cyprinidae),”  The Journal of Physiology, vol. 185, no. 3, pp. 587–599, 1966.
  * [7] M. Martinez-Garcia et al.,  “Derivatives and inverse of cascaded linear+ nonlinear neural models,”  PloS One, vol. 13, no. 10, pp. e0201326, 2018.
  * [8] A. Watson et al.,  “Model of visual contrast gain control and pattern masking,”  JOSA, vol. 14, no. 9, pp. 2379–2391, 1997.
  * [9] I. Epifanio et al.,  “Linear transform for simultaneous diagonalization of covariance and perceptual metric matrix in image coding,”  Pattern Recognition, vol. 36, pp. 1799–1811, 08 2003.
  * [10] J. Malo et al.,  “Nonlinear image representation for efficient perceptual coding,”  TIP, vol. 15, pp. 68–80, 2006.
  * [11] J. Ballé et al.,  “End-to-end optimized image compression,”  in 2017 ICLR. 2017, OpenReview.net.
  * [12] J. Gutiérrez et al.,  “Regularization operators for natural images based on nonlinear perception models,”  TIP, vol. 15, no. 1, pp. 189–200, 2005.
  * [13] V. Laparra et al.,  “Perceptually optimized image rendering,”  JOSA, vol. 34, no. 9, pp. 1511–1525, Sep 2017.
  * [14] A.M. Pons et al.,  “Image quality metric based on multidimensional contrast perception models,”  Displays, vol. 20, no. 2, pp. 93–110, 1999.
  * [15] V. Laparra et al.,  “Divisive normalization image quality metric revisited,”  JOSA, vol. 27, no. 4, pp. 852–864, Apr 2010.
  * [16] A. Hepburn et al.,  “Perceptnet: A human visual system inspired neural network for estimating perceptual distance,”  in 2020 ICIP. IEEE, 2020, pp. 121–125.
  * [17] R. Coen-Cagli et al.,  “The impact on midlevel vision of statistically optimal divisive normalization in V1,”  Journal of vision, vol. 13, no. 8, pp. 13–13, 2013.
  * [18] L.G. Giraldo et al.,  “Integrating flexible normalization into midlevel representations of deep convolutional neural networks,”  Neural Computation, vol. 31, pp. 2138–2176, 2019.
  * [19] X. Pan et al.,  “Brain-inspired weighted normalization for CNN image classification,”  bioRxiv, 2021.
  * [20] M. Michelle et al.,  “Divisive feature normalization improves image recognition performance in alexnet,”  in 2022 ICLR, 2022.
  * [21] O. Schwartz et al.,  “Natural signal statistics and sensory gain control,”  Nature Neuroscience, vol. 4, no. 8, pp. 819–825, 2001.
  * [22] J. Malo et al.,  “V1 non-linear properties emerge from local-to-global non-linear ICA,”  Network: Computation in Neural Systems, vol. 17, no. 1, pp. 85–102, 2006.
  * [23] J. Malo et al.,  “Psychophysically tuned divisive normalization approximately factorizes the pdf of natural images,”  Neural Computation, vol. 22, no. 12, pp. 3179–3206, 2010.
  * [24] J. Ballé et al.,  “Density modeling of images using a generalized normalization transformation,”  in 2016 ICLR, 2016.
  * [25] J. Malo,  “Spatio-chromatic information available from different neural layers via gaussianization,”  J. Math. Neurosci., vol. 10, no. 1, pp. 1–40, 2020.
  * [26] J.M. Hillis et al.,  “Do common mechanisms of adaptation mediate color discrimination and appearance? uniform backgrounds,”  JOSA, vol. 22, no. 10, pp. 2090–2106, 2005.
  * [27] M.D. Fairchild,  Color Appearance Models,  IS and T Wiley Series. Wiley, Sussex, UK, 2013.
  * [28] M. Martinez-Garcia et al.,  “In praise of artifice reloaded: Caution with natural image databases in modeling vision,”  Frontiers in Neuroscience, vol. 13, 2019.
  * [29] M. Ren et al.,  “Normalizing the normalizers: Comparing and extending network normalization schemes,” 2017.
  * [30] M.F. Burg et al.,  “Learning divisive normalization in primary visual cortex,”  PLOS Computational Biology, vol. 17, no. 6, pp. e1009028, 2021.
  * [31] M. Cordts et al.,  “The cityscapes dataset for semantic urban scene understanding,”  in 2016 CVPR, 2016.
  * [32] C. Sakaridis et al.,  “Semantic foggy scene understanding with synthetic data,”  Int. J. Comput. Vis., vol. 126, no. 9, pp. 973–992, 2018.


# Maximum-likelihood Soft-decision Decoding for Binary Linear Block Codes
Based on Their Supercodes

Yunghsiang S. Han1, Hung-Ta Pai2, Po-Ning Chen3 and Ting-Yi Wu3  1Dep. of
Electrical Eng. National Taiwan University of Science and Technology, Taipei,
Taiwan
Email<EMAIL_ADDRESS>2Dep. of Communication Eng. National Taipei
University, Taipei, Taiwan  3Dep. of Electrical and Computer Eng., National
Chiao Tung University, Hsinchu, Taiwan

###### Abstract

Based on the notion of supercodes, we propose a two-phase maximum-likelihood
soft-decision decoding (tpMLSD) algorithm for binary linear block codes in
this work. The first phase applies the Viterbi algorithm backwardly to a
trellis derived from the parity-check matrix of the supercode of the linear
block code. Using the information retained from the first phase, the second
phase employs the priority-first search algorithm to the trellis corresponding
to the linear block code itself, which guarantees finding the ML decision.
Simulations on Reed-Muller codes show that the proposed two-phase scheme is an
order of magnitude more efficient in average decoding complexity than the
recursive maximum-likelihood decoding (RMLD) [1] when the signal-to-noise
ratio per information bit is $4.5$ dB.

##  I Introduction

Linear block codes have been deployed for error control in communication
systems for many years, while algebraic structures of such codes are generally
used for their decoding [2]. Since the inputs of algebraic decoders are
commonly required to be quantized into two levels, they are classified as
hard-decision decoding technique. In comparison with soft-decision decoding
technique, a loss of information is induced due to quantization and hence the
decoding performance is restricted.

By contrast, the soft-decision decoding is developed to eliminate the
performance loss due to quantization. The input of soft-decision decoding is
thus unquantized (or practically quantized into more than two levels). In the
literature, many maximum-likelihood (ML) soft-decision decoding algorithms for
linear block codes have been proposed [3, 4, 5, 6, 7, 8, 9, 1, 10, 11], and
the priority-first search algorithm (PFSA) is one of them [6]. It has been
shown in [6] that the PFSA can provide the optimal ML decoding performance
within practically acceptable decoding complexity.

In this paper, a novel two-phase maximum-likelihood soft-decision decoding
(tpMLSD) scheme based on supercodes of linear block codes is proposed.
Specifically, in the first phase, the Viterbi algorithm (VA) is applied to a
trellis, derived from the parity-check matrix of the supercode of the linear
block code to be decoded, in a backward fashion (i.e., operated from the last
trellis level to the first trellis level). Upon the completion of the first
phase, each state will retain a path metric that is used later in the second
phase. Because the trellis derived from the parity-check matrix of the
supercode of the linear block code has fewer states than that derived from the
parity-check matrix of the linear block code itself, the computational
complexity is considerably reduced.

In the second phase, the priority-first search algorithm is applied to the
trellis corresponding to the parity-check matrix of the linear block code.
With a properly designed evaluation function, the optimal ML decision is
guaranteed to be located. Notably, the path metric information obtained from
the first phase is incorporated into the evaluation function for priority-
first search, by which the decoding procedure can be significantly sped up.
Simulations on Reed-Muller codes are then performed to confirm the efficiency
of the proposed two-phase ML soft-decision decoding scheme.

It should be pointed out that the idea of decoding linear block codes based on
their super codes is not new in the literature. It has been used in the hard-
decision decoding in [12], where super codes are designed based upon covering
sets and split syndromes. In addition, a suboptimal hard-decision list
decoding of linear block codes based on trellises of supercodes was presented
in [13]. Further generalization of [13] to ML soft-decision list decoding and
to soft-output decoding can be found in [14] and [15], respectively.

The rest of this paper is organized as follows. Notions of supercodes and ML
soft-decision decoding of linear block codes are introduced in Section II. The
proposed two-phase ML soft-decision decoding algorithm for binary linear block
codes based on their supercodes is presented in Section III. The optimality of
the proposed algorithm is proved in Section IV. Section V evaluates the
complexity of the proposed algorithm for practical linear block codes, and
Section VI concludes the paper.

##  II Notions of Supercodes and ML Soft-decision Decoding of Linear Block
Codes

Let ${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$ be an $(n,k)$ binary
linear block code with parity-check matrix $\mathbb{H}$ of size $(n-k)\times
n$. Denote by $\overline{}\mbox{${\textstyle\,{\cal
C}\\!\\!\\!\\!\\!\\!\\!\sim}$}$ an $(n,\bar{k})$ supercode of
${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$ with parity-check matrix
$\overline{\mathbb{H}}$ of size $(n-\bar{k})\times n$, satisfying that

$\mathbb{H}=\begin{bmatrix}\ \overline{\mathbb{H}}\ \ \\\
\mathbb{P}\end{bmatrix}$ (1)

for some matrix $\mathbb{P}$ of size $(\bar{k}-k)\times n$, where $\bar{k}>k$.

A trellis corresponding to linear block code ${\textstyle\,{\cal
C}\\!\\!\\!\\!\\!\\!\\!\sim}$ can then be constructed below. Denote by
$\boldsymbol{h}_{j}$, $0\leq j\leq n-1$, the $(j+1)$th column of $\mathbb{H}$.
Let $\boldsymbol{v}=(v_{0},v_{1},\ldots,v_{n-1})$ denote a codeword of
${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$. By defining recursively a
sequence of states $\\{\boldsymbol{s}_{\ell}\\}_{\ell=-1}^{n-1}$ as:

$\boldsymbol{s}_{\ell}=\begin{cases}{\bf 0},&\ell=-1\\\
\boldsymbol{s}_{\ell-1}+v_{\ell}\boldsymbol{h}_{\ell},&\ell=0,1,\ldots,n-1,\end{cases}$

a path corresponding to codeword $\boldsymbol{v}$ on a trellis $\mathscr{T}$
of $(n+1)$ levels can be identified, where ${\bf 0}$ is the all-zero vector of
proper size. Obviously,

$\boldsymbol{s}_{\ell}=\sum_{j=0}^{\ell}v_{j}\boldsymbol{h}_{j}\ \ \text{ for
}\ell=0,1,\ldots,n-1$

and

$\boldsymbol{s}_{n-1}={\bf 0}\ \ \text{ for all codewords of
\mbox{${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$}}.$

The trellis $\mathscr{T}$ derived from $\mathbb{H}$ is then formed by picking
up all paths corresponding to codewords of ${\textstyle\,{\cal
C}\\!\\!\\!\\!\\!\\!\\!\sim}$.

By convention, state $\boldsymbol{s}_{\ell}$ identifies a node on trellis
$\mathscr{T}$ at level $\ell$. In particular, $\boldsymbol{s}_{-1}$ and
$\boldsymbol{s}_{n-1}$ identify the initial node and the final node on trellis
$\mathscr{T}$ at levels $-1$ and $n-1$, respectively. In addition, the branch
connecting state $\boldsymbol{s}_{\ell-1}$ and state $\boldsymbol{s}_{\ell}$
is labeled with code bit $v_{\ell}$. As such, the one-to-one mapping between
codewords of ${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$ and paths over
$\mathscr{T}$ is built. This completes the construction of trellis
$\mathscr{T}$ based on parity-check matrix $\mathbb{H}$. The super-trellis
$\overline{\mathscr{T}}$ corresponding to supercode
$\overline{}\mbox{${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$}$ and its
parity-check matrix $\overline{\mathbb{H}}$ can be similarly constructed, of
which its state at level $\ell$ is denoted by $\bar{\boldsymbol{s}}_{\ell}$.

We next introduce the ML soft-decision decoding for codes with trellis
representation. Denote again by
$\boldsymbol{v}\triangleq(v_{0},v_{1},\ldots,v_{n-1})$ a binary zero-one
codeword of ${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$. Define the
hard-decision sequence $\boldsymbol{y}=(y_{0},y_{1},\ldots,y_{n-1})$
corresponding to the received vector
$\boldsymbol{r}=(r_{0},r_{1},\ldots,r_{n-1})$ as

$y_{j}\triangleq\begin{cases}1,&\text{if}\ \phi_{j}<0;\\\
0,&\text{otherwise,}\end{cases}$

where

$\phi_{j}\triangleq\log{{\Pr(r_{j}|0})\over\Pr(r_{j}|1)}$

is the log-likelihood ratio, and $\Pr(r_{j}|0)$ and $\Pr(r_{j}|1)$ are the
conditional probabilities of receiving $r_{j}$ given $0$ and $1$ were
transmitted, respectively. Here, $\Pr(r_{j}|0)$ can be either a probability
density function (pdf) for continuous (unquantized) $r_{j}$ or a probability
mass function (pmf) for discrete (softly quantized) $r_{j}$.

The syndrome of $\boldsymbol{y}$ is given by
$\boldsymbol{y}\mathbb{H}^{\textsf{T}}$, where superscript “T” denotes the
matrix transpose operation. Let $E(\boldsymbol{a})$ be the collection of all
error patterns whose syndrome is $\boldsymbol{a}$. Then, the maximum-
likelihood (ML) decoding output $\hat{\boldsymbol{v}}$ for received vector
$\boldsymbol{r}$ satisfies:

$\hat{\boldsymbol{v}}=\boldsymbol{y}\oplus\boldsymbol{e}^{\ast},$

where
$\boldsymbol{e}^{\ast}=(e_{0}^{\ast},e_{1}^{\ast},\ldots,e_{n-1}^{\ast})\in
E(\boldsymbol{y}\mathbb{H}^{\textsf{T}})$ is the error pattern satisfying

$\sum^{n-1}_{j=0}\ e^{\ast}_{j}|\phi_{j}|\leq\sum^{n-1}_{j=0}\
e_{j}|\phi_{j}|$

for all $\boldsymbol{e}=(e_{0},e_{1},\ldots,e_{n-1})\in
E(\boldsymbol{y}\mathbb{H}^{\textsf{T}})$, and “$\oplus$” denotes component-
wise modulo-two addition. We thereby define a new metric for paths in a
trellis as follows.

######  Definition 1 (ML path metric)

For a path with labels
$\boldsymbol{x}_{(\ell)}=(x_{0},x_{1},\ldots,x_{\ell})$, which ends at level
$\ell$ on trellis $\mathscr{T}$, define the metric associated with it as

$M\left(\boldsymbol{x}_{(\ell)}\right)\triangleq\sum^{\ell}_{j=0}M(x_{j}),$

where $M(x_{j})\triangleq(y_{j}\oplus x_{j})|\phi_{j}|$ is the bit metric.
Similarly, for a backward path with labels
$\bar{\boldsymbol{x}}_{[\ell]}=(\bar{x}_{\ell},\bar{x}_{\ell+1},\ldots,\bar{x}_{n-1})$
on super-trellis $\overline{\mathscr{T}}$, define the metric associated with
it as

$M\left(\bar{\boldsymbol{x}}_{[\ell]}\right)\triangleq\sum^{n-1}_{j=\ell}M(\bar{x}_{j}).$
(2)

After giving the notions of supercode and super-trellis as well as path
metrics, we proceed to present the proposed two-phase decoding scheme in the
next section.

##  III Two-phase ML Soft-Decision Decoding Algorithm for Binary Linear Block
Codes

As mentioned in the introduction section, the proposed decoding algorithm has
two phases.

The first phase applies the Viterbi algorithm backwardly to the supe-trellis
derived from the parity-check matrix $\overline{\mathbb{H}}$ of supercode
$\overline{}\mbox{${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$}$ using
the path metric defined in (2), during which the path metric of the backward
survivor starting from the final node at level $n-1$ and ending at a node
corresponding to state $\bar{\boldsymbol{s}}_{\ell}$ at level $\ell$ is
retained for use in the second phase. For convenience of referring it, we
denote this path metric by $c(\bar{\boldsymbol{s}}_{\ell})$. At the end of the
first phase, a backward survivor path ending at the initial node at level $-1$
is resulted. The backward Viterbi algorithm in the first phase is summarized
below.

$\langle$Phase 1: The backward Viterbi Algorithm$\rangle$

  * Step 1.

Associate zero initial metric with the backward path111It is clear that a path
on a trellis can not only be identified by its labels, but also be determined
by the states it traverses. Accordingly, path $\bar{\boldsymbol{x}}_{[\ell]}$
can be equivalently designated by
$(\bar{\boldsymbol{s}}_{\ell},\bar{\boldsymbol{s}}_{\ell+1},\ldots,\bar{\boldsymbol{s}}_{n-1})$.
containing only the final state $\bar{\boldsymbol{s}}_{n-1}$ on super-trellis
$\overline{\mathscr{T}}$, and let $c(\bar{\boldsymbol{s}}_{n-1})=0$. Set
$\ell=n-1$.

  * Step 2.

Decrease $\ell$ by one. Compute the metrics for all backward paths extending
from the backward survivors ending at level $\ell+1$ (and hence entering a
state at level $\ell$). For each state $\bar{\boldsymbol{s}}_{\ell}$ at level
$\ell$, keep the entering path with the least metric as its survisor, and
delete the remaining. Let $c(\bar{\boldsymbol{s}}_{\ell})$ be this least
metric.

  * Step 3.

If $\ell=0$, stop the algorithm; otherwise, go to Step 2.

In the second phase, the priority-first search algorithm is operated on
trellis $\mathscr{T}$ in the usual forward fashion (i.e., from level $0$ to
level $n-1$); hence, the second phase always outputs a codeword in
${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$.

Now for each path with labels $\boldsymbol{x}_{(\ell)}$ on trellis
$\mathscr{T}$, an evaluation function $f$ associated with it is defined as:

$f\left(\boldsymbol{x}_{(\ell)}\right)=g\left(\boldsymbol{x}_{(\ell)}\right)+h\left(\boldsymbol{x}_{(\ell)}\right),$

where the value of $g$-function is assigned according to:

$g\left(\boldsymbol{x}_{(\ell)}\right)=\begin{cases}0,&\ell=-1;\\\
g\left(\boldsymbol{x}_{(\ell-1)}\right)+M(x_{(\ell)}),&\ell=0,1,\ldots,n-1\end{cases}$
(3)

and the value of $h$-function is given by:

$h\left(\boldsymbol{x}_{(\ell)}\right)=c\left(\beta(\boldsymbol{s}_{\ell})\right).$
(4)

In (4), $\boldsymbol{s}_{\ell}$ is the ending state of the path with label
$\boldsymbol{x}_{(\ell)}$, and $\beta(\boldsymbol{s}_{\ell})$ is the state
$\bar{\boldsymbol{s}}_{\ell}$ on super-trellis $\overline{\mathscr{T}}$ that
has the same first $(n-\bar{k})$ components as $\boldsymbol{s}_{\ell}$. Note
that $\beta(\boldsymbol{s}_{\ell})$ exists and is well-defined for every
$\boldsymbol{s}_{\ell}$ on trellis $\mathscr{T}$ since the parity-check
matrices of ${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$ and
$\overline{}\mbox{${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$}$ satisfy
(1).

It can be verified that $f(\boldsymbol{x}_{(n-1)})=g(\boldsymbol{x}_{(n-1)})$
since $\beta(\boldsymbol{s}_{n-1})={\bf 0}$ and
$h\left(\boldsymbol{x}_{(n-1)}\right)=c(\beta(\boldsymbol{s}_{n-1}))=0$. This
implies that the path with the minimum $f$-function value on trellis
$\mathscr{T}$ is exactly the one with the minimum ML path metric.

Two storage spaces are necessary for the priority-first search over trellis
$\mathscr{T}$. The Open Stack records the paths visited thus far by the
priority-first search, while the Close Table keeps the starting and ending
states and ending levels of the paths that have ever been on top of the Open
Stack. They are so named because the paths in the Open Stack can be further
extended and hence remain open, but the paths with information in the Closed
Table are closed for further extension.

We summarize the priority-first search algorithm over trellis $\mathscr{T}$ in
the following.

$\langle$Phase 2: The Priority-First Search Algorithm$\rangle$

  * Step 1.

Let $\rho=\infty$, and assign $\boldsymbol{x}=\emptyset$.

  * Step 2.

Load into the Open Stack the path containing only the initial state
$\boldsymbol{s}_{-1}$ at level $-1$.

  * Step 3.

If the Open Stack is empty, output $\boldsymbol{x}$ as the final ML decision,
and stop the algorithm.

  * Step 4.

If the starting and ending states and ending level of the top path in the Open
Stack have been recorded in the Close Table, discard the top path from the
Open Stack, and go to Step 3; otherwise, record the starting and ending states
and ending level of this top path in the Close Table.

  * Step 5.

Compute the $f$-function values of the successors of the top path in the Open
Stack, and delete the top path from the Open Stack. If the $f$-function value
of any successor is equal to or greater than $\rho$, just delete it.

  * Step 6.

For all remaining successor paths that reach level $n-1$, set $\rho$ to be the
least path metric among them, and update $\boldsymbol{x}$ as the successor
path corresponding to this least path metric and discard all the others.

  * Step 7.

Insert the remaining successor paths (from Steps 5 and 6) into the Open Stack,
and re-order the paths in the Open Stack according to ascending $f$-function
values. Go to Step 3.

##  IV Optimality of the Proposed Algorithm

This section proves the optimality of the proposed two-phase decoding
algorithm. We begin with two essential lemmas required for the optimality
proof.

###### Lemma 1

Let path $\boldsymbol{x}_{(\ell+1)}$ be an immediate successor of path
$\boldsymbol{x}_{(\ell)}$ on trellis $\mathscr{T}$. Denote the ending states
of $\boldsymbol{x}_{(\ell+1)}$ and $\boldsymbol{x}_{(\ell)}$ by
$\boldsymbol{s}_{\ell+1}$ and $\boldsymbol{s}_{\ell}$, respectively. Then,

$\beta(\boldsymbol{s}_{\ell+1})=\beta(\boldsymbol{s}_{\ell})+x_{\ell+1}\bar{\boldsymbol{h}}_{\ell+1},$

where $\bar{\boldsymbol{h}}_{\ell+1}$ is the $(\ell+2)$th column of parity-
check matrix $\overline{\mathbb{H}}$. In other words, there exists a branch
between $\beta(\boldsymbol{s}_{\ell})$ and $\beta(\boldsymbol{s}_{\ell+1})$
with label $x_{\ell+1}$ over super-trellis $\overline{\mathscr{T}}$.

###### Proof:

Recall that

$\boldsymbol{s}_{\ell+1}=\boldsymbol{s}_{\ell}+x_{\ell+1}\boldsymbol{h}_{\ell+1}$

and

$\boldsymbol{h}_{\ell+1}=\begin{bmatrix}\ \bar{\boldsymbol{h}}_{\ell+1}\ \\\ \
\boldsymbol{p}_{\ell+1}\ \end{bmatrix}$

for some $\boldsymbol{p}_{\ell+1}$ according to (1). It is thus obvious that

$\beta(\boldsymbol{s}_{\ell+1})=\beta(\boldsymbol{s}_{\ell})+x_{\ell+1}\bar{\boldsymbol{h}}_{\ell+1}$

since $\bar{\boldsymbol{h}}_{\ell+1}$ contains the first $(n-\bar{k})$
components of $\boldsymbol{h}_{\ell+1}$. ∎

###### Lemma 2

$f$ is a non-decreasing function along any path on trellis $\mathscr{T}$,
i.e.,

$f\left(\boldsymbol{x}_{(\ell)}\right)\leq
f\left(\boldsymbol{x}_{(\ell+1)}\right),$

where path $\boldsymbol{x}_{(\ell+1)}$ is an immediate successor of path
$\boldsymbol{x}_{(\ell)}$ over trellis $\mathscr{T}$.

###### Proof:

The fundamental attribute of the backward Viterbi algorithm in the first phase
gives that $c(\beta(\boldsymbol{s}_{\ell}))$ is the minimum metric among all
backward paths that end at state $\beta(\boldsymbol{s}_{\ell})$ at level
$\ell$. By Lemma 1, we have:

$c(\beta(\boldsymbol{s}_{\ell}))\leq
c(\beta(\boldsymbol{s}_{\ell+1}))+M(x_{\ell+1}),$

where $\boldsymbol{s}_{\ell+1}$ and $\boldsymbol{s}_{\ell}$ are respectively
the states that paths $\boldsymbol{x}_{(\ell+1)}$ and
$\boldsymbol{x}_{(\ell)}$ end at. Hence, we derive:

$\displaystyle f\left(\boldsymbol{x}_{(\ell+1)}\right)$ $\displaystyle=$
$\displaystyle
g\left(\boldsymbol{x}_{(\ell+1)}\right)+h\left(\boldsymbol{x}_{(\ell+1)}\right)$
$\displaystyle=$ $\displaystyle
g\left(\boldsymbol{x}_{(\ell)}\right)+M(x_{\ell+1})+c(\beta(\boldsymbol{s}_{\ell+1}))$
$\displaystyle\geq$ $\displaystyle
g\left(\boldsymbol{x}_{(\ell)}\right)+c(\beta(\boldsymbol{s}_{\ell}))$
$\displaystyle=$ $\displaystyle f\left(\boldsymbol{x}_{(\ell)}\right).$

∎

Based on these two lemmas, the next theorem proves the optimality of the
proposed two-phase algorithm.

###### Theorem 1

In the second phase, the priority-first search algorithm always output an ML
path.

###### Proof:

It suffices to prove that if the Open Stack is empty, the algorithm will
output an ML path as claimed in Step 3. This can be confirmed by showing that
Steps 4 and 5 never delete any ML path.

Suppose that in Step 4, the starting and ending states and ending level of the
new top path $\boldsymbol{x}_{(\ell)}$ have been recorded in the Close Table
at some previous time due to path $\hat{\boldsymbol{x}}_{(\ell)}$. Since path
$\boldsymbol{x}_{(\ell)}$ must be an offspring of some path
$\boldsymbol{x}_{(j)}$ that once coexisted with path
$\hat{\boldsymbol{x}}_{(\ell)}$ in the Open Stack at the time path
$\hat{\boldsymbol{x}}_{(\ell)}$ was on top of the Open Stack, where $j<\ell$,
we have

$f\left(\boldsymbol{x}_{(\ell)}\right)\geq
f\left(\boldsymbol{x}_{(j)}\right)\geq
f\left(\hat{\boldsymbol{x}}_{(\ell)}\right).$ (5)

Notably, the first inequality in (5) follows from Lemma 2, and the second
inequality in (5) is valid because the top path in the Open Stack always
carries the minimum $f$-function value among all coexisting paths. As a
result, the offsprings of path $\boldsymbol{x}_{(\ell)}$ ending at level $n-1$
cannot yield smaller metrics than those length-$n$ offsprings of path
$\hat{\boldsymbol{x}}_{(\ell)}$, and hence deletion of path
$\boldsymbol{x}_{(\ell)}$ will not compromise the optimality of the decoding
algorithm.

For Step 5, we argue that $\rho$ is either a trivial upper bound of the final
ML path (cf. Step 1) or the metric of a valid path that reaches level $n-1$
(cf. Step 6), so deletion of any successor paths whose $f$-function values are
no less than $\rho$ will never eliminate any ML path. This completes the proof
of optimality of the proposed algorithm. ∎

##  V Evaluation of Computational Efforts

In this section, we investigate by simulations the computational effort of the
proposed decoding algorithm over the additive white Gaussian noise (AWGN)
channels. We assume that the codeword is antipodally modulated, and hence the
received vector is given by

$r_{j}=(-1)^{v_{j}}\sqrt{\cal E}+\lambda_{j},$

for $0\leq j\leq n-1$, where ${\cal E}$ is the signal energy per channel bit,
and $\\{\lambda_{j}\\}_{j=0}^{n-1}$ are independent noise samples of a white
Gaussian process with single-sided noise power per hertz $N_{0}$. The signal-
to-noise ratio (SNR) for the channel is therefore given by
$\text{SNR}\triangleq{\cal E}/N_{0}$. In order to account for the code
redundancy for different code rates, we will use the SNR per information bit
in the following discussion, which is defined as

$\text{SNR}_{\text{b}}=\frac{n{\cal E}/k}{N_{0}}=\frac{n}{k}\left(\frac{\cal
E}{N_{0}}\right).$

It can be easily verified that for antipodal-input AWGN channels, the log-
likelihood ratio $\phi_{j}$ is a fixed multiple of the received scalar
$r_{j}$; thus, the metric associated with a path $\boldsymbol{x}_{(\ell)}$ can
be equivalently simplified to

$M\left(\boldsymbol{x}_{(\ell)}\right)\triangleq\sum^{\ell}_{j=0}(y_{j}\oplus
x_{j})|r_{j}|,$

where

$y_{j}\triangleq\begin{cases}1,&\text{if}\ r_{j}<0;\\\
0,&\text{otherwise.}\end{cases}$

The decoding complexity in the first phase is clearly determined by the number
of bit metric computations performed. We emphasize that the decoding
complexity in the second phase can also be regulated by the number of
$f$-function evaluations (equivalently, the number of bit metric computations
as indicated in (3)) during the priority-first search. This is due to that the
cost of searching and re-ordering of stack elements can be made a constant
multiple of the computational complexity by adopting a priority-queue data
structure in stack implementation [16]. One can even employ a hardware-based
stack structure [17] and attain constant complexity in stack maintenance.
Therefore, to use the number of overall metric computations as the key
determinant of algorithmic complexity for our proposed two-phase decoding
algorithm is justified.

We now turn to empirical examination of the average decoding complexity of the
proposed tpMLSD algorithm. The linear block code considered is the $r$th order
binary Reed-Muller code, RM$(r,m)$, which is an $(n,k)$ linear block code with
$n=2^{m}$ and $k=1+\sum_{i=1}^{r}{m\choose i}$. It is known [18] that
RM$(r+i,m)$ is a supercode of RM$(r,m)$ for $i\geq 1$. In our simulations,
${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$ is RM$(2,6)$ and
$\overline{}\mbox{${\textstyle\,{\cal C}\\!\\!\\!\\!\\!\\!\\!\sim}$}$ is
RM$(4,6)$; hence, $n=64$, $k=22$ and $\bar{k}=57$. Under the same optimal ML
performance, we compare the proposed two-phase ML soft-decision decoding
(tpMLSD) algorithm with the recursive ML decoding (RMLD) algorithm [1] and the
list ML decoding (LMLD) algorithm [14] in average decoding complexity, and
summarize the results in Table I.

Note that instead of listing the decoding complexity of the LMLD, lower bounds
obtained from decoding its supercode counterpart using the marking algorithm
are given [14]. Apparently, the real decoding complexity of the LMLD is higher
than this lower bound. The table then shows that the LMLD is much more complex
than the other two algorithms, and our two-phase decoding algorithm consumes
only 1/13 of the computational effort of the RMLD at
$\text{SNR}_{\text{b}}=4.5$ dB, in which circumstance the bit error rate (BER)
is around $10^{-5}$. Further, when $\text{SNR}_{\text{b}}$ is reduced to $3$
dB, the average computational complexity of the proposed two-phase decoding
scheme can still reach $1/8$ of that of the RMLD.

TABLE I: Average computational complexities (I.e., average number of metrics
evaluated) of the RMLD, the LMLD, and the tpMLSD. The linear block code
considered is RM$(2,6)$, while the suppercode used in the tpMLSD is RM$(4,6)$.

$\text{SNR}_{\text{b}}$ | 3 dB | 3.5 dB | 4 dB | 4.5 dB | 5 dB
---|---|---|---|---|---
RMLD [1] | 78209 | 78209 | 78209 | 78209 | 78209
*LMLD [14] | 2097152 | 2097152 | 2097152 | 2097152 | 2097152
tpMLSD | 10078 | 7863 | 6602 | 6010 | 5695

*What are listed here are lower bounds to the decoding complexities of the LMLD.

##  VI Conclusion

In this work, we proposed a two-phase scheme for ML soft-decision decoding of
linear block codes. This novel decoding algorithm has two phases, where the
backward Viterbi algorithm is employed on a supercode of the linear block code
in the first phase, while the priority-first search algorithm is performed on
the trellis of the linear block code in the second phase. Simulations showed
that the computational complexity of the proposed two-phase scheme is one
order of magnitude better than that of the RMLD when
$\text{SNR}_{\text{b}}=4.5$ dB. Since such a new approach can be extended to
decoding any linear block codes when their supercodes are obtainable, a
possible future work is to extend this two-phase decoding scheme to codes like
Reed-Solomon, for which maximum-likelihood soft-decision decoding is generally
considered a challenging task.

## References

  * [1] T. Fujiwara, H. Yamamoto, T. Kasami, and S. Lin,  “A trellis-based recursive maximum-likelihood decoding algorithm for binary linear block codes,”  IEEE Trans. Inform. Theory, pp. 714–729, March 1998.
  * [2] S. Lin and D. J. Costello,  Error Control Coding,  Prentice-Hall, Upper Saddle River, NJ, second edition, 2004.
  * [3] J. K. Wolf,  “Efficient maximum likelihood decoding of linear block codes using a trellis,”  IEEE Trans. Inform. Theory, pp. 76–80, January 1978.
  * [4] J. Snyders,  “Reduced lists of error patterns for maximum likelihood soft decoding,”  IEEE Trans. Inform. Theory, pp. 1194–1200, July 1991.
  * [5] N. J. C. Lous, P. A. H. Bours, and H. C. A. van Tilborg,  “On maximum likelihood soft-decision decoding of binary linear codes,”  IEEE Trans. Inform. Theory, pp. 197–203, January 1993.
  * [6] Y. S. Han, C. R. P. Hartmann, and C.-C. Chen,  “Efficient priority-first search maximum-likelihood soft-decision decoding of linear block codes,”  IEEE Trans. Inform. Theory, vol. 39, no. 5, pp. 1514–1523, September 1993.
  * [7] T. Kaneko, T. Nishijima, H. Inazumi, and S. Hirasawa,  “An efficient maximum-likelihood decoding algorithm for linear block codes with algebraic decoder,”  IEEE Trans. Inform. Theory, pp. 320–327, March 1994.
  * [8] T. Kaneko, T. Nishijima, and S. Hirasawa,  “An improvement of soft-decision maximum-likelihood decoding algorithm using hard-decision bounded-distance decoding,”  IEEE Trans. Inform. Theory, pp. 1314–1319, July 1997.
  * [9] D. Gazelle and J. Snyders,  “Reliability-based code-search algorithm for maximum-likelihood decoding of block codes,”  IEEE Trans. Inform. Theory, pp. 239–249, January 1997.
  * [10] L. E. Aguado and P. G. Farrell,  “On hybrid stack decoding algorithms for block codes,”  IEEE Trans. Inform. Theory, pp. 398–409, January 1998.
  * [11] Y. S. Han,  “A new treatment of priority-first search maximum-likelihood soft-decision decoding of linear block codes,”  IEEE Trans. Inform. Theory, vol. 44, no. 7, pp. 3091–3096, November 1998.
  * [12] A. Barg, E. Krouk, and H. C. A. van Tilborg,  “On the complexity of minimum distance decoding of long linear codes,”  IEEE Trans. Inform. Theory, pp. 1392–1405, July 1999.
  * [13] J. Freudenberger,  “On bounded distance list decoding based on supercodes,”  in International Symposium on Communication Theory and Applications, Ambleside, Lake District, UK, 2003.
  * [14] J. Freudenberger and M. Bossert,  “Maximum-likelihood decoding based on supercodes,”  in International ITG Conference Source and Channel Coding, Erlangen, Germany, 2004.
  * [15] J. Freudenberger and U. Kaiser,  “Symbol-by-symbol APP decoding based on supercode decoding,”  in International Zurich Seminar on Communications, Zurich, Switzerland, 2012.
  * [16] T. H. Cormen, C. E. Leiserson, and R. L. Rivest,  Introduction to Algorithms,  MIT Press, Cambridge, MA, 1991.
  * [17] P. Lavoie, D. Haccoun, and Y. Savaria,  “A systolic architecture for fast stack sequential decoders,”  IEEE Trans. Commun., vol. 42, no. 5, pp. 324–335, May 1994.
  * [18] F. J. MacWilliams and N. J. A. Sloane,  The Theory of Error-Correcting Codes,  New York, NY: Elsevier Science Publishing Company, Inc., 1977.