# CTMQ: Cyclic Training of Convolutional Neural Networks with Multiple
Quantization Steps

HyunJin Kim
School of Electronics and Electrical Engineering,
Dankook University
152, Jukjeon-ro, Suji-gu, Yongin-si, Gyeonggi-do, Republic of Korea
<EMAIL_ADDRESS>
&Jungwoo Shin
School of Electronics and Electrical Engineering,
Dankook University
152, Jukjeon-ro, Suji-gu, Yongin-si, Gyeonggi-do, Republic of Korea
<EMAIL_ADDRESS>
Alberto A. Del Barrio
Department of Computer Architecture and System Engineering, Complutense
University of Madrid, Madrid, ES
<EMAIL_ADDRESS>
Corresponding Author, homepapge: https://www.empaslab.com

###### Abstract

This paper proposes a training method having multiple cyclic training for
achieving enhanced performance in low-bit quantized convolutional neural
networks (CNNs). Quantization is a popular method for obtaining lightweight
CNNs, where the initialization with a pretrained model is widely used to
overcome degraded performance in low-resolution quantization. However, large
quantization errors between real values and their low-bit quantized ones cause
difficulties in achieving acceptable performance for complex networks and
large datasets. The proposed training method softly delivers the knowledge of
pretrained models to low-bit quantized models in multiple quantization steps.
In each quantization step, the trained weights of a model are used to
initialize the weights of the next model with the quantization bit depth
reduced by one. With small change of the quantization bit depth, the
performance gap can be bridged, thus providing better weight initialization.
In cyclic training, after training a low-bit quantized model, its trained
weights are used in the initialization of its accurate model to be trained. By
using better training ability of the accurate model in an iterative manner,
the proposed method can produce enhanced trained weights for the low-bit
quantized model in each cycle. Notably, the training method can advance Top-1
and Top-5 accuracies of the binarized ResNet-18 on the ImageNet dataset by
5.80% and 6.85%, respectively.

##  1 Introduction

In recent years, quantized CNNs are widely adopted for reducing the hardware
complexity (Gupta et al. (2015); Han et al. (2015); Wu et al. (2016); Hubara
et al. (2017); Wang et al. (2018); Wu et al. (2020); Yu et al. (2020)).
Moreover, it has been proved that low-bit quantized CNNs using binary (Hubara
et al. (2016); Courbariaux et al. (2016); Rastegari et al. (2016); Liu et al.
(2018)), ternary (Zhou et al. (2016); Deng et al. (2017); Wan et al. (2018))
quantizations can be trainable, producing acceptable performance. However,
there is a significant performance gap between real-valued CNNs and low-bit
quantized CNNs. Quantization noise disturbs weight updating in the
quantization-aware training (Wang et al. (2018)), thus degrading the
performance of trained models in low-bit quantized CNNs. The retraining method
uses pretrained weights to initialize a target model, and the initialized
model is retrained. Because the retraining method automatically searches for
suitable parameter values, it can produce better classification results
compared with naive post-training quantization (Gholami et al. (2021)). It is
known that the initial weights from accurate models are helpful to produce
better training results. However, low-bit quantized CNNs show significant
performance degradation immediately after the initialization.

In each batch of training, after updating real-valued weights in back-
propagation, quantized weights are used in forward paths, making quantization
noise. The quantization noise hinders the training step to find its optimized
weights (Helwegen et al. (2019); Xu et al. (2021a)). Besides, the scaling and
biasing parameters used in normalization have disturbances due to the gap
between real-valued parameters and their quantized values. There have been
several works to overcome the degraded performance in low-bit quantized models
(Helwegen et al. (2019); Kim et al. (2020); Martinez et al. (2020); Xu et al.
(2021a, b)). However, the existing studies for low-bit quantized models do not
consider the training method for achieving better weight initialization from
accurate models.

This paper proposes a method consisting of multiple cyclic training for
achieving better performance in low-bit quantized CNNs. Firstly, the proposed
method softly transfers the learned knowledge based on accurate models to low-
bit quantized models with multiple quantization steps. In a quantization step,
the pretrained weights are used to initialize the weights of the next model
with the quantization bit depth reduced by one. With small change of the
quantization bit depth, the performance gap can be bridged, thus producing
better weight initialization. In the multiple quantization steps, the
pretrained weights in a step are used in the initialization of weights to be
trained in the next quantization step. As proceeding with each quantization
step, the quantization bit depth of trained models decreases gradually.
Secondly, the proposed cyclic training uses the trained weights of a low-bit
quantized model in the weight initialization of its accurate model to be
trained. Then, after training the accurate model, its trained weights are used
in the training of the low-bit quantized model in an iterative manner. The
proposed method uses better training ability of the accurate model in each
cycle, producing better weight initialization for the low-bit quantized model
in an iterative manner. The proposed training method is applied to the
training of the binarized ResNet-18 (He et al. (2016)). Quantized CNNs based
on DoReFa-Net (Zhou et al. (2016)) and XNOR-Net (Rastegari et al. (2016)) are
trained with the proposed method. Without any structural changes in ResNet-18,
the proposed training method outperforms the training results from scratch and
using the initialization of pretrained real-valued models. In the experiments
on the ImageNet dataset, Top-1 and Top-5 accuracies of the binarized ResNet-18
are enhanced by 5.80% and 6.85%, respectively.

This paper is organized as follows: in the preliminaries, the difficulties in
training highly approximate CNNs are described in detail. Then, the proposed
cyclic training having multiple quantization steps is explained along with a
detailed process description. Next, we explain our experimental environments
and experimental results for training the binarized ResNet-18.

##  2 Preliminaries

###  2.1 Quantized CNNs

As the easiest way to reduce the computational resources and storage
requirements of CNNs, quantized CNNs have been studied recently. Normalized
activations and weights used in the forward path of CNNs are expressed using a
fixed-point format. Quantized CNNs quantize activations and weights to have
8-bit, 4-bit, 2-bit, ternary, and binary formats. The simplified format and
its lightweight operations can dramatically reduce the hardware complexity of
convolutional operations in CNNs, compared with real-valued CNNs. Notably, the
binarized CNNs utilize 1-bit activations and weights in binarized
convolutions.

It has been proven in several existing studies that quantized CNNs can achieve
a certain level of inference performance (Gholami et al. (2021)). However, due
to quantization errors, their performance can be degraded compared with real-
valued CNNs. Notably, there may be a significant performance degradation in
the quantized CNN using weights based on a format of less than 2-bit (Zhou et
al. (2016)). Therefore, it is necessary to improve the performance of
quantized CNNs.

When a pretrained CNN model is given, the pretrained model can be used to
initialize its quantized CNN model. Quantization of CNNs has been categorized
depending on whether retrained is performed or not. Post-training quantization
reuses the parameters of the pretrained model without retraining, where
quantization levels are uniformly spaced or unequally determined to minimize
the effects of quantization errors (Gholami et al. (2021)). However, due to
the non-linearity in each layer of the CNN model, it is difficult to guarantee
the best value by reducing the amount of quantization errors. On the other
hand, the retraining scheme (Han et al. (2015)) uses the existing pretrained
model for the weight initialization of a model and performs training on the
quantized CNN model. In quantization-aware training (Wang et al. (2018)), the
error of the forward path in the quantized model is considered as the amount
of loss. Besides, the parameters for scaling and normalizing activations can
be optimized during retraining. Therefore, it is known that the performance
after retraining can be higher than that without retraining in quantized CNN
models.

###  2.2 Training of Quantized CNNs and its Difficulties

In general, the objective function of a CNN is a non-convex loss function
denoted as $\mathcal{L}$, where training tries to find a local optimum using
the gradient descent method (Ruder (2016)). The gradient descent is described
in Algorithm 1.

Input : initial weights $w^{(0)}$, number of iterations $T$, learning rate
$\gamma$

Output : final weights $w^{(T)}$

1 for  _$t\leftarrow 0$ to $T-1$_ do

2 $w^{(t+1)}:=w^{(t)}+\gamma\nabla$ $\mathcal{L}$(_$w^{(t)}$_)

3 end for

Return(_$w^{(T)}$_).

Algorithm 1 Gradient Descent

Although other optimizers such as ADAM (Kingma and Ba (2014)) have their own
formulations, the gradient descent was explained using (1) for its simplicity.
The local optimum is repeatedly searched by updating the gradient of
$\mathcal{L}$. The updated weights $w^{(t+1)}$ are considered when obtaining
the next $\mathcal{L}$ in an iteration.

On the other hand, in quantized CNNs, the quantization error could nullify or
suppress the effects of the updating with $\gamma\nabla$
$\mathcal{L}$(_$w^{(t)}$_) in Algorithm 1. Weights $w^{(t)}$ at the $t$-th
iteration can be rewritten using $k$-bit quantized weights $w^{(t)}_{q_{k}}$
and their following error $\epsilon(w^{(t)}_{q_{k}})$ as:

$w^{(j)}_{q_{k}}:=w^{(j)}-\epsilon(w^{(j)}_{q_{k}}),0\leq t\leq T.$ (1)

From (1), $\epsilon(w^{(t)}_{q_{k}})$ and $\epsilon(w^{(t+1)}_{q_{k}})$ affect
the quantized weights used in the next forward path. The equation in line 2 of
Algorithm 1 can be rewritten as:

$w^{(t+1)}_{q_{k}}:=w^{(t)}_{q_{k}}+\epsilon(w^{(t)}_{q_{k}})+\gamma\nabla\mathcal{L}(w^{(t)}_{q_{k}})-\epsilon(w^{(t+1)}_{q_{k}}).$
(2)

From (2), we can conclude two facts: firstly, terms
$\epsilon(w^{(t)}_{q_{k}})$ and $-\epsilon(w^{(t+1)}_{q_{k}})$ have different
signs, which could somewhat cancel out the quantization error of the updated
weights $\epsilon(w^{(t+1)}_{q_{k}})$.

Secondly, entire values of $\gamma\nabla\mathcal{L}(w{(t)}^{q_{k}})$ in (2),
cannot be considered in the weight updating when $\epsilon(w^{(t)}_{q_{k}})$
and $-\epsilon(w^{(t+1)}_{q_{k}})$ are not the same. When $k$ becomes small in
low-bit quantized CNNs, the cancellation can be hard, which can mislead
searches for the local optimum in training steps. It is assured that when a
model has high quantization bit depth with small
$\epsilon(w^{(t+1)}_{q_{k}})$, it can produce good training results.

Several training methods for quantized CNNs empirically adopt a straight-
through estimator (STE) (Bengio et al. (2013); Yin et al. (2019)) for back-
propagating gradients instead of using derivatives of the quantizer. For
example, binarized CNNs can use the $sign$ function to binarize activations
(Zhou et al. (2016); Rastegari et al. (2016); Liu et al. (2018)). The
derivative of the $sign$ function can be the Dirac Delta function (Hassani
(2009)), which is hard to implement. Using a STE, the approximation of the
derivative of the $sign$ function can be shown by (3).

$\nabla sign(x)\approx\left\\{\begin{matrix}&+1&\textrm{if }|x|\leq 1,\\\
&0&\textrm{else.}\end{matrix}\right.$ (3)

Instead of using the Dirac Delta function, low-bit quantized CNNs adopt the
STE, making the incoming gradients equal to their output gradients in the
back-propagation.

In spite of long training time, a STE in low-bit quantized and binarized CNNs
can produce acceptable performance in many several works (Zhou et al. (2016);
Rastegari et al. (2016); Liu et al. (2018); Kim (2021); Shin and Kim (2022)).
However, term $\epsilon(w^{(t)}_{q_{k}})$ in (1) still exists, producing
instability in training steps.

##  3 CTMQ: cyclic training with multiple quantization steps

###  3.1 Proposed training method

As discussed in the preliminaries section, accurate models with small
quantization errors can provide better training performance along with more
stable weight convergence. In existing works of binarized CNNs (Liu et al.
(2018); Bulat and Tzimiropoulos (2019); Kim et al. (2020), pretrained weights
are used in the weights initialization of the binarized models, which can
transfer the training results of accurate models into low-bit quantized
models. Let $w(i)^{(T)}$ be the pretrained weights after training an $i$-bit
quantized CNN model. Let us assume that the pretrained weights are used in
weight initialization for $k$-bit quantized model. The weight updating in (2)
can be rewritten by $w^{(0)}=w(i)^{(T)}$ as:

$w^{(t+1)}_{q_{k}}:=w(i)^{(T)}_{q_{k}}+\epsilon(w(i)^{(T)}_{q_{k}})+\gamma\nabla\mathcal{L}(w(i)^{(t)}_{q_{k}})-\epsilon(w^{(t+1)}_{q_{k}}).$
(4)

Depending on the difference between $i$ and $k$, the distribution of quantized
values can be different. For example, small weights in a Gaussian distribution
can be quantized into non-zero values in accurate models. On the other hand,
when quantization bit depth becomes low in (2), many small weights are
quantized into zeros. The significantly different quantization bit depths
could prevent the knowledge of pretrained models from being well transferred
into low-bit quantized models.

To overcome the problem, we proposed a new method called CTMQ, which is an
abbreviation for cyclic training with multiple quantization steps. In the
proposed method, the knowledge of a pretrained model is softly transferred
into its low-bit quantized model using weight initialization. Then, the
knowledge transfer using initial weights is cyclically repeated between $k$
and $(k+1)$-bit quantized models. In the cyclic training, the pretrained
weights of the $k$-bit quantized model can be used in the initialization for
training the $(k+1)$-bit quantized model in an iterative manner. Then, the
$(k+1)$-bit quantized model is trained with the initial weights. After
obtaining the pretrained weights of the $(k+1)$-bit quantized model, the
pretrained weights can be used the training of the $k$-bit quantized model.
The cyclic training is repeated several times. In each cycle, based on the
initial weights from the $k$-bit quantized model, the training of the
$(k+1)$-bit quantized model can produce better training weights for the
initialization of the $k$-bit quantized model. The proposed CTMQ for training
a $k$-bit quantized model is described in Algorithm 2.

Input : initial weights $w^{(0)}$, number of cycles $C$, numbers of iterations
$T_{s},T_{c},T_{f}$, learning rate $\gamma$, $n$-bit quantizations $q_{n}$

Output : final weights $w^{(T_{l})}$

1 for  _$n\leftarrow 8$ to $k+2$_ do // 1st part: soft knowledge transfer

2 for  _$t\leftarrow 0$ to $T_{s}-1$_ do

3 $w^{(t+1)}=w^{(t)}+\gamma\nabla$ $\mathcal{L}_{q_{n}}$(_$w^{(t)}_{q_{n}}$_);

4

5 end for

6 $w^{(0)}=w^{(T_{s})}$;

7

8 end for

9for  _$c\leftarrow 0$ to $C-1$_ do // 2nd part: cyclic training with
$q_{k+1}$ and $q_{k}$

10 for  _$n\leftarrow k+1$ to $k$_ do

11 for  _$t\leftarrow 0$ to $T_{c}-1$_ do

12 $w^{(t+1)}=w^{(t)}+\gamma\nabla$
$\mathcal{L}_{q_{n}}$(_$w^{(t)}_{q_{n}}$_);

13

14 end for

15 $w^{(0)}=w^{(T_{c})}$;

16

17 end for

18

19 end for

20for  _$t\leftarrow 0$ to $T_{c}-1$_ do

21 $w^{(t+1)}=w^{(t)}+\gamma\nabla$
$\mathcal{L}_{q_{n}}$(_$w^{(t)}_{q_{k+1}}$_);

22

23 end for

24$w^{(0)}=w^{(T_{c})}$;

25 for  _$t\leftarrow 0$ to $T_{f}-1$_ do // 3rd part: full training with
$q_{k}$

26 $w^{(t+1)}=w^{(t)}+\gamma\nabla$
$\mathcal{L}_{q_{n}}$(_$w^{(t)}_{q_{k}}$_);

27

28 end for

29Return(_$w^{(T_{l})}=w^{(T_{f})}$_).;

Algorithm 2 CTMQ: cyclic training with multiple quantization steps for $k$-bit
quantized CNNs

In Algorithm 2, there are three parts: In the first part, after training an
$n$-bit quantized model with $T_{s}$ iterations, the pretrained model is used
in the initialization of an $(n-1)$-bit quantized model. The initial weights
from the pretrained real-valued model can be used as $w^{(0)}$. In $n$-bit
quantization, both weights and activations are quantized into $n$ bits. This
iteration commented as soft knowledge transfer is repeated until $n=k+2$
(lines 1-6). By letting $w^{(0)}$ be $w^{(T_{s})}$ (line 5), the training
results of an accurate model are delivered into an inaccurate model. In the
second part, cyclic training with $q_{k+1}$ and $q_{k}$ quantized models are
performed during $C$ cycles (lines 7-18). Alternatively, $(k+1)$-bit and
$k$-bit quantized models are trained during $T_{c}$ iterations. The trained
weights are used to initialize the next quantized model by
$w^{(0)}=w^{(T_{c})}$ (line 12). Then, a $(k+1)$-bit quantized model is
trained after finishing the cyclic training (lines 15-17). In the third part,
the trained weights with the final $(k+1)$-bit quantized model is used to
initialize the final $k$-bit quantized model. Finally, the $k$-bit quantized
model is trained during $T_{f}$ iterations.

###  3.2 Target quantizations

In the multiple quantization steps of the proposed method, we adopts the
quantization scheme shown in DoReFa-Net (Zhou et al. (2016)) for quantizing
weights and activations, which is summarized in the following. It is noted
that the proposed training cannot depend on any specific quantization method.

Except for binarized weights, weights $w$ are normalized into $w_{norm}$ for
obtaining $k$-bit quantized weights as:

$w_{norm}=\frac{w}{2*max(abs(tanh(w)))}+0.5.$ (5)

Terms $tanh$ and $abs$ mean the hyperbolic tangent and absolute functions,
respectively. Notably, $tanh$ function is used to limit the value range. Then,
normalized weights $w_{norm}$ are quantized into $w_{q_{k}}$ as:

$w_{q_{k}}=2\times\frac{round(w_{norm}\times(2^{k}-1))}{2^{k}-1}-1.$ (6)

Besides, binarized weights can be obtained using the $sign$ function as:

$w_{q_{1}}=\left\\{\begin{matrix}&sign(w/mean(abs(w)))*mean(abs(w))\textrm{ if
}w\neq 0,\\\ &1\textrm{ else.}\end{matrix}\right.$ (7)

In (7), when $w=0$, let $sign(w/mean(w))$ be $1$ to avoid the cases with
$w_{q_{1}}=0$.

On the other hand, $k$-bit quantized activations denoted as $x_{q_{k}}$ are
obtained as:

$x_{q_{k}}=\frac{round(clamp(x,0,1)\times(2^{k}-1)))}{2^{k}-1}.$ (8)

In (8), term $clamp(x,0,1)$ denotes the clamp function to limit $x$ to the
range $[0,1]$. Unlike the $1$-bit quantized weights, the proposed method does
not use the $sign$ function for $1$-bit quantized activations. Therefore,
$x_{q_{1}}$ can be $0$ or $1$.

A value $2^{k}-1$ is multiplied by the denominator and the numerator in (6)
and (8). In (7), $mean(abs(w))$ is multiplied by both the denominator and
numerator. The range of values in $k$-bit quantized weights are limited by
normalizing weights in (5). In $k$-bit quantized activations, the clamp
function in (8) limits the range of values. Regardless of $k$, the range of
values is fixed so that two quantized models with different $k$ have the same
range of weights and activations.

##  4 Experimental results and analysis

In experiments, existing binarized CNN models were trained to evaluate the
proposed method on the training of low-bit quantized CNNs. Binarized CNNs
adopted $1$-bit weights and activations in the convolutions so that
significant computational and storage resources can be saved. However, the
performance degradation of binarized CNNs has been the most critical issue
(Qin et al. (2020)). It is expected that the proposed method can solve this
problem by achieving high classification accuracy. The pyramid structure with
highly stacked convolution layers and residual shortcuts have been adopted in
many studies of binarized CNNs (Rastegari et al. (2016); Lin et al. (2017);
Liu et al. (2018); Bulat and Tzimiropoulos (2019); Phan et al. (2020); Liu et
al. (2020); Kim (2021); Xu et al. (2021a); Shin and Kim (2022)). Considering
the popularity of the residual binarized networks, our experiments adopted a
residual CNN model called ResNet-18 He et al. (2016). Standardized datasets
(CIFAR-100 (Krizhevsky et al. (2009))) and ImageNet (Russakovsky et al.
(2015))) were in the evaluations. ResNet-18 was quantized for achieving the
final trained model with binarized (1-bit) weights and activations. For fair
comparison, we adopted the quantization scheme in Zhou et al. (2016). Besides,
we maintained the structure of basic blocks in DoReFa-Net (Zhou et al. (2016))
and XNOR-Net (Rastegari et al. (2016) in the experiments.

###  4.1 Experimental frameworks

We evaluated the proposed training method using the Pytorch deep learning
frameworks (Paszke et al. (2019)). We coded our training script that supported
multiple GPU-based data parallelism. After finishing training, its real-valued
trained model was stored in a file with pth extension. As shown in Algorithm
2, the file having pretrained weights was loaded for the initialization used
in the next training. The above process was repeated in the proposed cyclic
training using shell scripts. It must be noted that quantizations were not
implemented in silicon-based hardware circuits. Instead, quantizers were coded
as functions to emulate the quantizations of weights and activations in the
forward path.

###  4.2 Target dataset and model structures

Figure 1:  ResNet-18 for ImageNet dataset.  Figure 2:  Basic blocks for
quantized ResNet-18: (a) TypeI; (b) TypeII.

The summary of target datasets is as follows: the CIFAR-100 dataset consists
of 60K $32\times 32$ colour images with 100 classes, where each class contains
500 training and 100 test images. For the data augmentation during training,
$32\times 32$ input images were randomly cropped and flipped from $40\times
40$ padded images. The ImageNet dataset contains 1.3M training and 50K
validation images with 1K classes. During training on the ImageNet dataset,
images were randomly resized between 256 and 480 pixels, and then randomly
cropped and flipped into $224\times 224$ images. In inference, we adopted
$224\times 224$ center-cropped images without data augmentations from the
validation dataset.

Figure 1 describes the structure of ResNet-18 for the ImageNet dataset. The
layers have 5 groups from Conv-1 to Conv-5, where the convolutional layer in a
group has the same number of output channels. Firstly, 3 channels from RGB
coloured image are used as $7\times 7$ convolutional layer with $stride=2$
(Conv-1). With $stride=2$, the width and height of output features are reduced
to half. After the max pooling with $stride=2$, the convolutional layers of
Conv-2 are performed with 64 input channels.For the CIFAR dataset, a $3\times
3$ convolutional layer with $stride=1$ is used in Conv-1. A shortcut with
direct mapping are described as a thin round arrow. Per two convolutional
layers, one shortcut is placed so that a basic block contains two
convolutional layers and one shortcut. After finishing two basic blocks, the
downsampling is performed by doubling the number of channels and having
$stride=2$. The thick round arrow denotes the downsampling using an $1\times
1$ convolution. This process is repeated until finishing convolutions in
Conv-5. Then, the average pooling is applied to 512 channels with $7\times 7$
features. After performing fully-connected and softmax layers, the final
classification results can be obtained. If the quantization is adopted in the
first convolution layer of Conv-1, the original $8$-bit image pixel data could
be degraded. Because the width and height of features in Conv-5 are small,
large quantization errors can be critical in the fully-connected layer. In our
experiments, the convolution in Conv-1 and fully-connected layers adopted
real-valued weights and activations like Zhou et al. (2016); Rastegari et al.
(2016).

Figure 2 illustrates two basic blocks denoted as TypeI and TypeII. Our
experments designed the TypeI and TypeII basic blocks based on DoReFa-Net
(Zhou et al. (2016)) and XNOR-Net (Rastegari et al. (2016)). Term BN denotes
the batch normalization (Ioffe and Szegedy (2015)). The learnable mean and
variance were used to scale and shift the normalized activations. In the two
convolutional layers and downsampling, $k$-bit weights and activations were
used except for the first convolution in Conv-2. Besides, whereas TypeI
adopted $k$-bit quantization in the downsampling, real-valued weights and
activations were used in the downsampling of TypeII. In TypeII, after the
scaling of the batch normalization layer, $k$-bit quantized activations were
obtained, where the quantization in (8) was used instead of the $sign$
function in the original XNOR-Net.

###  4.3 Training on CIFAR-100 and ImageNet datasets

Training on the CIFAR-100 dataset used ADAM optimizer (Kingma and Ba (2014))
having $\beta_{1}=0.9\&\beta_{2}=0.999$ without weight decay. When training a
$k$-bit quantized model during $E_{epochs}$ epochs, the base learning rate
$lr_{base}$ was set as 0.001. The learning rate $lr$ was decayed based on poly
policy in the $e_{epochs}$-th epoch, limiting the maximum learning rate of the
ADAM optimizer by $lr_{base}\times(1-e_{epochs}/E_{epochs})$. After finishing
the training of an $n$-bit quantized model, $lr_{base}$ of the next training
was reset as 0.001.

The inputs in (2) were given as follows: pretrained real-valued weights were
used as the initial weights $w^{(0)}$. Let $C$ and $q_{n}$ be 9 and 8,
respectively. Therefore, an $8$-bit quantized model was firstly trained. When
finishing the first part in (2), a trained $3$-bit quantized model was
obtained. In the second part, $2$-bit and $1$-bit quantized models were
alternatively trained during 9 cycles. Its mini-batch size was set as 512. 20
training epochs were adopted for $T_{s}$ and $T_{c}$ so that
$T_{s}=T_{c}=1,960$. The third part peformed 200 training epochs for
$T_{f}=19,600$.

Several hyperparameters used in the training on the ImageNet dataset were the
same as those on the CIFAR-100 dataset. Besides, $C$ and $q_{n}$ were set as 9
and 8 in the training on the ImageNet dataset. Considering a large number of
images on the ImageNet dataset, let $T_{s}=T_{c}=12,515$ for 5 training epochs
in the first and second parts. In the third part, 50 training epochs were
performed for $T_{f}=125,150$.

Figure 3:  Inference accuracies of quantized ResNet-18 models on ImageNet
dataset.

Figure 3 illustrates inference accuracies of quantized ResNet-18 models on the
ImageNet dataset during training. At each epoch, the inference accuracies of
two binarized ResNet-18 models based on TypeI and TypeII were shown. In Torch
(2022), it must be noted that 1-crop Top-1 and Top-5 inference accuracies of
the original ResNet-18 were 69.76% and 89.08%, respectively. At the starting
point of the first part, the quantization error lowered the inference
accuracy. As training continued, the accuracies rapidly increased. In the
first part, the accuracies from TypeI and TypeII were nearly close, which
means that the effects of the quantized downsampling were not significant in
the first part. In the second part, compared with those of TypeII, the models
using TypeI had low inference accuracies immediately after $1$-bit quantized
models were initialized. Therefore, it was concluded that the $1$-bit
quantization of the downsampling can degrade performance in binarized CNNs.
The initialization of binarized CNNs with pretrained $2$-bit quantized models
was performed at the 36-th, 46-th, and 56-th epochs. The inference accuracies
at the epochs increased as the cyclic training proceeded. For example, the
accuracies at the 46-th epoch were better than those at the 36-th epoch so
that better initial weights can be obtained in the cyclic training. After the
third cycle, the increase slowed down, as shown in Figure 3. In the third
part, compared with the earlier epochs of the second part, the $1$-bit
quantized CNNs at the third part can start their training with higher
inference accuracies after weight initializations. As training proceeded, the
inference accuracies gradually increased in the third part.

Table 1: Comparison of $k$-bit quantized models in terms of inference accuracies Dataset | Model | CTMQ | $k$ | Initialization | Top-1 | Top-5
---|---|---|---|---|---|---
CIFAR-100 | TypeI | N | 32 | N | 73.91% | -
N | 8 | N | 68.28% | -
N | 1 | N | 64.13% | -
N | 1 | Y | 65.41% | -
Y | 1 | - | 67.93% | -
TypeII | N | 32 | N | 72.61% | -
N | 8 | N | 67.21% | -
N | 1 | N | 67.69% | -
N | 1 | Y | 68.89% | -
Y | 1 | Y | 68.84% | -
ImageNet | ResNet-18 | N | 32 | N | 69.76% | 89.08%
TypeI | N | 32 | N | 67.40% | 87.71%
N | 1 | N | 47.51% | 72.71%
N | 1 | Y | 42.38% | 67.91%
Y | 1 | Y | 51.54% | 76.15%
TypeII | N | 32 | N | 67.03% | 87.27%
XNOR-Net | N | 1 | N | 51.2% | 73.20%
TypeII | Y | 1 | N | 57.0% | 80.05%

Table 1 lists the comparison in terms of inference accuracies. The proposed
training method was evaluated based on the average of five runs. The
accuracies of ResNet-18 and XNOR-Net in Table 1 were referenced from Torch
(2022) and Rastegari et al. (2016). There were no published results on the
quantized ResNet-18 on the CIFAR-100 dataset in Zhou et al. (2016) and
Rastegari et al. (2016) so that we trained the quantized models during 200
epochs on the CIFAR-100 dataset and 50 epochs on the ImageNet dataset for the
counterparts without applying the proposed method. In the training of the
counterparts, most of the hyperparameters were identical with those using the
proposed method. However, weight decay was set as 0.0001 for the
regularization. In Table 1, when the initialization was Y, real-valued
pretrained models were used in weight initialization. Besides, when item CTMQ
was N, the soft knowledge transfer and cyclic training were not applied.

Due to the small number of training images and small image size in the
CIFAR-100 dataset, overfitting was shown in the experiments of ResNet-18.
Therefore, the enhancements with the proposed method were limited in the
experiments on the CIFAR-100 dataset. On the other hand, the proposed method
significantly enhanced the final inference accuracies on the ImageNet dataset
compared with other cases. On the ImageNet dataset, when not using the
proposed method, TypeI with initialization showed significant performance drop
only having 42.38% Top-1 and 67.91% Top-5 inference accuracies. TypeI using
the proposed CTMQ achieved 51.54% Top-1 and 76.15% Top-5 accuracies, which
produced better performance over the original XNOR-Net Rastegari et al.
(2016). The binarized ResNet-18 with TypeII had Top-1 and Top-5 inference
accuracies up to 57.0% and 80.05%, respectively, which showed 5.80% and 6.85%
enhancements over the training results of the original XNOR-Net.

Considering the experimental results, it is concluded that the proposed
training method was effective for training complex residual networks on large
dataset.

##  5 Conclusion

We focus on a new training method performing cyclic training with different
$k$-bit quantizations. The training method performs multiple cyclic training
for softly delivering the knowledge of pretrained accurate models and using
better training ability of accurate networks. Although the number of training
epochs increases during three training parts, it has been proven that the
final inference accuracies for complex residual networks on large dataset are
significantly enhanced in experimental results. It is noted that there are no
structural changes in CNNs for enhancing the final inference accuracy.
Although several binarized CNNs are trained in experiments, the proposed
method can be applied to other quantized CNNs based on the existing
quantization and optimizer. Notably, the training method can advance Top-1 and
Top-5 accuracies of the binarized ResNet-18 on the ImageNet dataset by 5.80%
and 6.85%, respectively. Considering the enhanced classification outputs, it
is concluded that the proposed cyclic training method is useful for producing
high-performance systems in quantized CNNs.

## References

## References

  * Bengio et al. (2013) Yoshua Bengio, Nicholas Léonard, and Aaron Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation.  _arXiv preprint arXiv:1308.3432_ , 2013.
  * Bulat and Tzimiropoulos (2019) Adrian Bulat and Georgios Tzimiropoulos.  Xnor-net++: Improved binary neural networks.  _arXiv preprint arXiv:1909.13863_ , 2019.
  * Courbariaux et al. (2016) Matthieu Courbariaux, Itay Hubara, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks: Training deep neural networks with weights and activations constrained to+ 1 or-1.  _arXiv preprint arXiv:1602.02830_ , 2016.
  * Deng et al. (2017) Lei Deng, Peng Jiao, Jing Pei, Zhenzhi Wu, and Guoqi Li.  Gated xnor networks: deep neural networks with ternary weights and activations under a unified discretization framework.  _arXiv preprint arXiv:1705.09283_ , 2, 2017.
  * Gholami et al. (2021) Amir Gholami, Sehoon Kim, Zhen Dong, Zhewei Yao, Michael W Mahoney, and Kurt Keutzer.  A survey of quantization methods for efficient neural network inference.  _arXiv preprint arXiv:2103.13630_ , 2021.
  * Gupta et al. (2015) Suyog Gupta, Ankur Agrawal, Kailash Gopalakrishnan, and Pritish Narayanan.  Deep learning with limited numerical precision.  In _International Conference on Machine Learning_ , pages 1737–1746, 2015.
  * Han et al. (2015) Song Han, Huizi Mao, and William J Dally.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  _arXiv preprint arXiv:1510.00149_ , 2015.
  * Hassani (2009) Sadri Hassani.  Dirac delta function.  In _Mathematical methods_ , pages 139–170. Springer, 2009.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Helwegen et al. (2019) Koen Helwegen, James Widdicombe, Lukas Geiger, Zechun Liu, Kwang-Ting Cheng, and Roeland Nusselder.  Latent weights do not exist: Rethinking binarized neural network optimization.  In _Advances in neural information processing systems_ , pages 7531–7542, 2019.
  * Hubara et al. (2016) Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Binarized neural networks.  In _Advances in neural information processing systems_ , pages 4107–4115, 2016.
  * Hubara et al. (2017) Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Quantized neural networks: Training neural networks with low precision weights and activations.  _The Journal of Machine Learning Research_ , 18(1):6869–6898, 2017.
  * Ioffe and Szegedy (2015) Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _arXiv preprint arXiv:1502.03167_ , 2015.
  * Kim et al. (2020) Hyungjun Kim, Kyungsu Kim, Jinseok Kim, and Jae-Joon Kim.  Binaryduo: Reducing gradient mismatch in binary activation network by coupling binary activations.  _arXiv preprint arXiv:2002.06517_ , 2020.
  * Kim (2021) HyunJin Kim.  Aresb-net: accurate residual binarized neural networks using shortcut concatenation and shuffled grouped convolution.  _PeerJ Computer Science_ , 7:e454, 2021.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Lin et al. (2017) Xiaofan Lin, Cong Zhao, and Wei Pan.  Towards accurate binary convolutional neural network.  In _Advances in Neural Information Processing Systems_ , pages 345–353, 2017.
  * Liu et al. (2018) Zechun Liu, Baoyuan Wu, Wenhan Luo, Xin Yang, Wei Liu, and Kwang-Ting Cheng.  Bi-real net: Enhancing the performance of 1-bit cnns with improved representational capability and advanced training algorithm.  In _Proceedings of the European conference on computer vision (ECCV)_ , pages 722–737, 2018.
  * Liu et al. (2020) Zechun Liu, Zhiqiang Shen, Marios Savvides, and Kwang-Ting Cheng.  Reactnet: Towards precise binary neural network with generalized activation functions.  _arXiv preprint arXiv:2003.03488_ , 2020.
  * Martinez et al. (2020) Brais Martinez, Jing Yang, Adrian Bulat, and Georgios Tzimiropoulos.  Training binary neural networks with real-to-binary convolutions.  _arXiv preprint arXiv:2003.11535_ , 2020.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , pages 8024–8035, 2019.
  * Phan et al. (2020) Hai Phan, Yihui He, Marios Savvides, Zhiqiang Shen, et al.  Mobinet: A mobile binary network for image classification.  In _The IEEE Winter Conference on Applications of Computer Vision_ , pages 3453–3462, 2020.
  * Qin et al. (2020) Haotong Qin, Ruihao Gong, Xianglong Liu, Xiao Bai, Jingkuan Song, and Nicu Sebe.  Binary neural networks: A survey.  _Pattern Recognition_ , 105:107281, 2020.
  * Rastegari et al. (2016) Mohammad Rastegari, Vicente Ordonez, Joseph Redmon, and Ali Farhadi.  Xnor-net: Imagenet classification using binary convolutional neural networks.  In _European conference on computer vision_ , pages 525–542. Springer, 2016.
  * Ruder (2016) Sebastian Ruder.  An overview of gradient descent optimization algorithms.  _arXiv preprint arXiv:1609.04747_ , 2016.
  * Russakovsky et al. (2015) Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei.  ImageNet Large Scale Visual Recognition Challenge.  _International Journal of Computer Vision (IJCV)_ , 115(3):211–252, 2015.  doi: 10.1007/s11263-015-0816-y.
  * Shin and Kim (2022) Jungwoo Shin and HyunJin Kim.  Presb-net: parametric binarized neural network with learnable activations and shuffled grouped convolution.  _PeerJ Computer Science_ , 8:e842, 2022.
  * Torch (2022) Contributors Torch.  Models and Pre-trained weights.  https://pytorch.org/vision/stable/models.html, 2022.  [Online; accessed 17-May-2022].
  * Wan et al. (2018) Diwen Wan, Fumin Shen, Li Liu, Fan Zhu, Jie Qin, Ling Shao, and Heng Tao Shen.  Tbn: Convolutional neural network with ternary inputs and binary weights.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 315–332, 2018.
  * Wang et al. (2018) Naigang Wang, Jungwook Choi, Daniel Brand, Chia-Yu Chen, and Kailash Gopalakrishnan.  Training deep neural networks with 8-bit floating point numbers.  In _Advances in neural information processing systems_ , pages 7675–7684, 2018.
  * Wu et al. (2020) Hao Wu, Patrick Judd, Xiaojie Zhang, Mikhail Isaev, and Paulius Micikevicius.  Integer quantization for deep learning inference: Principles and empirical evaluation.  _arXiv preprint arXiv:2004.09602_ , 2020.
  * Wu et al. (2016) Jiaxiang Wu, Cong Leng, Yuhang Wang, Qinghao Hu, and Jian Cheng.  Quantized convolutional neural networks for mobile devices.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 4820–4828, 2016.
  * Xu et al. (2021a) Weixiang Xu, Qiang Chen, Xiangyu He, Peisong Wang, and Jian Cheng.  Improving binary neural networks through fully utilizing latent weights.  _arXiv preprint arXiv:2110.05850_ , 2021a.
  * Xu et al. (2021b) Zihan Xu, Mingbao Lin, Jianzhuang Liu, Jie Chen, Ling Shao, Yue Gao, Yonghong Tian, and Rongrong Ji.  Recu: Reviving the dead weights in binary neural networks.  _arXiv preprint arXiv:2103.12369_ , 2021b.
  * Yin et al. (2019) Penghang Yin, Jiancheng Lyu, Shuai Zhang, Stanley Osher, Yingyong Qi, and Jack Xin.  Understanding straight-through estimator in training activation quantized neural nets.  _arXiv preprint arXiv:1903.05662_ , 2019.
  * Yu et al. (2020) Haibao Yu, Tuopu Wen, Guangliang Cheng, Jiankai Sun, Qi Han, and Jianping Shi.  Low-bit quantization needs good distribution.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_ , pages 680–681, 2020.
  * Zhou et al. (2016) Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, He Wen, and Yuheng Zou.  Dorefa-net: Training low bitwidth convolutional neural networks with low bitwidth gradients.  _arXiv preprint arXiv:1606.06160_ , 2016.


# Statistical Properties of the log-cosh Loss Function Used in Machine
Learning

Resve A. Saleh  Dept. of Electrical and Computer Engineering
University of British Columbia
Vancouver, Canada
<EMAIL_ADDRESS>A.K.Md. Ehsanes Saleh  School of Mathematics and Statistics
Carleton University
Ottawa, Canada
<EMAIL_ADDRESS>
###### Abstract

This paper analyzes a popular loss function used in machine learning called
the log-cosh loss function. A number of papers have been published using this
loss function but, to date, no statistical analysis has been presented in the
literature. In this paper, we present the distribution function from which the
log-cosh loss arises. We compare it to a similar distribution, called the
Cauchy distribution, and carry out various statistical procedures that
characterize its properties. In particular, we examine its associated pdf,
cdf, likelihood function and Fisher information. Side-by-side we consider the
Cauchy and Cosh distributions as well as the MLE of the location parameter
with asymptotic bias, asymptotic variance, and confidence intervals. We also
provide a comparison of robust estimators from several other loss functions,
including the Huber loss function and the rank dispersion function. Further,
we examine the use of the log-cosh function for quantile regression. In
particular, we identify a quantile distribution function from which a maximum
likelihood estimator for quantile regression can be derived. Finally, we
compare a quantile M-estimator based on log-cosh with robust monotonicity
against another approach to quantile regression based on convolutional
smoothing.

###### Index Terms:

log-cosh function, machine learning, distribution function, quantile
regression

##  I Introduction

According to several authors [1][2][3][4], the log-cosh loss function is one
of the most important loss functions in machine learning today but very little
has been published about its statistical characteristics. It has been
implemented in many software environments for machine learning and used in a
number of different research efforts. From a programming perspective, it may
be found in R in the limma package [5] and can be implemented in a single line
of python using the numpy library [6]. It is also available in TensorFlow2 in
the Keras library [7], and likewise in PyTorch [8] for deep learning.

Research areas of application include variational autoencoders [9][10] and
cancer detection [11]. It has also been used in tree-based learning algorithms
such as XGBoost[3], and most recently for a solution to the crossing problem
in quantile regression [12]. One caveat is that a function like
$\log(\cosh(x))$ can overflow if not implemented correctly. However, this
issue has been addressed through built-in library functions and the problem is
rarely encountered. As a result, it has been routinely used in machine
learning for a number of years now.

The log-cosh loss function belongs to the class of robust estimators that tend
to prefer solutions in the vicinity of the median rather than the mean.
Another view is that robust estimators are more tolerant to outliers in the
data set and this is perhaps one of the key reasons to select the log-cosh
loss function over others. It is also continously differentiable, unlike the
Huber loss function which does not have a continuous second derivative. These
aspects are well-known in the machine learning community. The gap is the
literature on log-cosh is relative to its origin and statistical properties.
In this paper, we seek to remedy this gap by deriving the log-cosh loss
function from first principles, starting with its distribution, and studying
properties such as bias, variance, confidence intervals, Fisher information
and standard errors of estimation.

##  II M-estimators

The log-cosh loss function falls in the category of M-estimators having the
general form:

$Q(x,\theta)=\frac{1}{n}\sum_{i=1}^{n}\rho(x_{i},\theta)$ (1)

where $\rho(x_{i},\theta)$ is a term in a given loss function. Here $X$ is an
iid random variable associated with the residuals and $\theta$ is the
parameter to be estimated. Any desired function may be used for
$\rho(x_{i},\theta)$ as long as it satisfies some minimal set of conditions.
The estimate of $\theta$ is produced as follows:

$\hat{\theta}=\text{argmin}_{\theta\in\mathbb{R}}\frac{1}{n}\sum_{i=1}^{n}\rho(x_{i},\theta).$
(2)

This is typically carried out using numerical optimization employing some form
of gradient descent in the context of machine learning. Of particular interest
here is the log-cosh loss function given by

$\rho_{L}(x_{i},\theta)=\text{log}(\text{cosh}(x_{i}-\theta)).$ (3)

In the sections to follow, we show how to derive this loss function starting
with a probability distribution and then provide a more intuitive view of its
origin.

##  III Background

###  III-A Hyperbolic Functions

Recall the basic hyperbolic functions as follows:

$\text{sinh}(x)=\frac{e^{x}-e^{-x}}{2}$

$\text{cosh}(x)=\frac{e^{x}+e^{-x}}{2}$

and

$\text{tanh}(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}.$

We also have that

$\text{sech}(x)=\frac{1}{\text{sinh}(x)}=\frac{2}{e^{x}-e^{-x}}$

$\text{csch}(x)=\frac{1}{\text{cosh}(x)}=\frac{2}{e^{x}+e^{-x}}$

and

$\text{coth}(x)=\frac{1}{\text{tanh}(x)}=\frac{e^{x}+e^{-x}}{e^{x}-e^{-x}}.$

###  III-B The Cosh Distribution

Of primary interest here is the cosh(x) function which is plotted in Fig.
1(a). The reciprocal of this function, i.e., $\text{csch}(x)$, is plotted in
Fig. 1(b).

Figure 1: (a) $\text{cosh}(x)$ (b) $1/\text{cosh}(x)$

This function resembles a probability density function (pdf) but in order to
satisfy the necessary conditions, it must be non-negative and integrate to 1.
The first requirement is already satisfied by inspection. For the second
requirement, if we take the integral of the function, we obtain the
normalizing constant as

$\int_{-\infty}^{\infty}\frac{1}{\cosh(x)}dx=\pi.$

Therefore, the pdf of the Cosh distribution is given by

$f(x)=\frac{1}{\pi\,\cosh(x)}$ (4)

The corresponding cumulative distribution function (cdf) is given by:

$F(x)=\frac{1}{2}+\frac{1}{\pi}{\text{tan}^{-1}\big{[}\sinh\big{(}{x}\big{)}}\big{]}$
(5)

which is invertible. The cdf and pdf are plotted in Fig. 2.

Figure 2: Cosh distribution (a) cdf (b) pdf

More generally, we have a location-scale family of distributions of the form

$f(x;\theta,\sigma)=\frac{1}{\pi\sigma\,\text{cosh}(\frac{x-\theta}{\sigma})}$
(6)

where $\theta$ is the location parameter and $\sigma$ is the scale parameter.
We note here that

$\mathbb{E}[X]=\theta\,\,\,\,\,\text{and}\,\,\,\,\,\mathbb{E}[X^{2}]=\frac{\pi^{2}\sigma^{2}}{4}+\theta^{2}$

so the asymptotic variance is

$\text{Var}(X)=\mathbb{E}[X^{2}]-\mathbb{E}[X]^{2}=\frac{\pi^{2}\sigma^{2}}{4}.$

The cdf given by:

$F(x;\theta,\sigma)=\frac{1}{2}+\frac{1}{\pi}\tan^{-1}\bigg{[}\sinh\bigg{(}\frac{x-\theta}{\sigma}\bigg{)}\bigg{]}.$
(7)

###  III-C Maximum Likelihood Estimation

We next develop the maximum likelihood estimator (MLE) for the case of
$\sigma=1$, as follows. If $x_{1},x_{2},...,x_{n}$ are i.i.d random variables,
then the likelihood function is given by:

$L(x_{1},x_{2},...,x_{n};\theta)=\prod_{i=1}^{n}\frac{1}{\pi\,\text{cosh}(x_{i}-\theta)}.$
(8)

The negative log-likelihood expression is given by:

$-\ell(x_{1},x_{2},...,x_{n};\theta)=\text{log}\pi+\sum_{i=1}^{n}\text{log}(\text{cosh}(x_{i}-\theta)).$
(9)

The estimator $\hat{\theta}$ is the solution to the equation:

$\sum_{i=1}^{n}\tanh(x_{i}-\theta)=0$ (10)

or equivalently

$\hat{\theta}=\text{argmin}_{\theta\in\mathbb{R}}\sum_{i=1}^{n}\text{log(cosh}(x_{i}-\theta)).$
(11)

We see that this estimator is, in fact, the MLE of the Cosh distribution. The
equation can be solved in a straight-forward manner using a convex
optimization procedure due to the fact that it is globally convex. That is,
the second derivative can be shown to be non-negative. This can be
demonstrated by noting that

$\psi_{L}(x)=\rho_{L}^{\prime}(x)=\tanh(x)$ (12)

and

$\psi_{L}^{\prime}(x)=\rho_{L}^{\prime\prime}(x)=\text{sech}^{2}(x).$ (13)

Therefore, both the gradient and Hessian terms can be easily obtained and this
is important for machine learning applications.

##  IV An Intuitive View of log-cosh

It is useful to consider the log-cosh loss function from an intuitive
standpoint to understand its relative importance. We first examine the loss
function in terms of its behavior relative to the L1 and L2 functions.
Consider M-estimators of the form

$\hat{\theta}=\text{argmin}_{\theta\in\mathbb{R}}\frac{1}{n}\sum_{i=1}^{n}\rho(x_{i},\theta).$
(14)

For example, we can define $\rho(x)$ in any manner we choose such as

$\rho_{L2}(x)=x^{2}$ (15)

for least squares estimate (LSE), also called the L2 loss function, and

$\rho_{L1}(x,\theta)=|x|$ (16)

for least absolute deviation (LAD), also called the L1 loss function. For log-
cosh, we already derived the MLE to be

$\rho_{L}(x)=\text{log}(\text{cosh}(x))$ (17)

Next, we compare the log-cosh loss function relative to the L1 and L2 cases
above. As $x\rightarrow+\infty$, the log-cosh loss tends towards

$\text{log(cosh}(x))\approx x-\text{log}(2)$ (18)

and as $x\rightarrow-\infty$, the function tends towards

$\text{log(cosh}(x))\approx-x-\text{log}(2).$ (19)

Therefore, far from $x=0$, the function behaves as

$\text{log(cosh}(x))\approx|x|-\text{log}(2)$ (20)

which is like L1 except for the constant term. Now, as $x\rightarrow 0$, we
can take a Taylor series expansion at $x=0$ as follows

$\text{log(cosh}(x))=\frac{x^{2}}{2}-\frac{x^{4}}{12}+\frac{x^{6}}{45}+O(x^{8})$
(21)

and therefore it is, to first-order, like L2. That is,

$\text{log(cosh}(x))\approx\frac{x^{2}}{2}$ (22)

for small $x$. Since it behaves like L2 close to the origin and L1 far from
the origin, one can view it as a smoothed out L1 using L2 around the origin.
Hence, both the Jacobian and Hessian matrix terms exist for this loss
function. For robust regression, this function can be used as an alternative
to L1.

Huber [13] proposed a method of combining the best of L1 and L2 by explicitly
using L2 in the vicinity of the origin where the discontinuity lies, and then
switching to L1 a certain distance, $\delta$, away from the origin. The small
interval around the origin is defined by $[-\delta,+\delta]$. Inside this
region, the L2 loss function is used since it is continuous. Outside this
region, the L1 loss is used but great care is taken to match the derivatives
at the interface between the two regions. The resulting Huber loss function is
given by:

$\hat{\theta}_{n}^{\rm
H}=\text{argmin}_{\theta}\sum_{i=1}^{n}\rho_{H}(x_{i}-\theta)$ (23)

where

$\rho_{H}(x)=\begin{cases}x^{2}/2&\quad\text{if }|x|\leq\delta\\\
\delta(|x|-\frac{\delta}{2})&\quad\text{if }|x|>\delta\end{cases}$ (24)

The derivative of the Huber function is given by:

$\psi_{H}(x)=\begin{cases}x&\quad\text{if }|x|\leq\delta\\\
+\delta&\quad\text{if }x>+\delta\\\ -\delta&\quad\text{if
}x<-\delta\end{cases}$ (25)

There are a total of 3 regions defined by the Huber function: $x<-\delta$,
$x>+\delta$, and $|x|<\delta$. As shown in Figure 3, the Huber function and
its first derivative are both continuous in all regions. In this case,
$\delta=1$ for illustrative purposes. However, the second derivative (not
shown) is discontinuous. It is clear that $\rho_{H}^{\prime\prime}(x)$ will be
discontinuous at both $-\delta$ and $+\delta$. Hence, the Hessian is not
defined at those points.

Figure 3: Huber loss and derivative as a function of $x$ for $\delta=1$.

Another intuitive view of log-cosh can be obtained by starting with the L1
loss function. This loss function produces the median as the estimate but it
has a discontinuous first derivative. We seek to find a way to construct a
continuous function that has continuous first and second derivatives to
replace the L1 function. Consider how one might do this intuitively starting
with the L1 loss function. This is shown in Fig. 4. The absolute value
function, $|x|$, has a “V” shape. When we take its derivative, we obtain
sgn($x$) which is discontinuous at $x=0$. We can replace the discontinuous
derivative function with a continuous version using $\tanh(x)$. They look
about the same but we know that one is discontinuous and the other is
continuous. Then, by integrating the $\tanh(x)$ function, we obtain
$\log(\cosh(x))$ which does not have a kink (although it may appear to have
one in the figure).

Figure 4: Developing a continuous L1 function.

##  V Statistical Properties

Consider the general form of the Cosh distribution of Eqn. (6). The MLE of
$\theta$ can be derived as given in Eqn. (11). The estimator is represented as
$\hat{\theta}$. Then, as $n\rightarrow\infty$, we find that

$\sqrt{n}(\hat{\theta}-\theta)\overset{\mathcal{D}}{\to}\mathcal{N}(0,\mathcal{I}^{-1}(\theta))$
(26)

where $\mathcal{I}(\theta)$ is the Fisher information which can be derived as

$\mathcal{I}(\theta)=\frac{1}{2\sigma^{2}}.$ (27)

Furthermore, the MLE estimator is consistent such that

$\hat{\theta}\overset{\mathcal{P}}{\to}\theta$ (28)

with variance

$\text{Var}(\hat{\theta})=\frac{1}{n\mathcal{I}(\theta)}.$

Hence, we find the asymptotic variance of $\hat{\theta}$ to be

$\text{Var}(\hat{\theta})=\frac{2\sigma^{2}}{n}.$ (29)

Next, given that $\hat{\theta}\overset{\mathcal{P}}{\to}\theta$, the
asymptotic bias is

$\text{bias}=\mathbb{E}[\hat{\theta}]-\theta=0.$

The log-cosh estimator is therefore asymptotically unbiased. Further, the
$(1-\alpha)$ confidence interval is given by

$CI=\bigg{[}\hat{\theta}-\frac{z_{\alpha/2}}{\sqrt{n\mathcal{I}(\theta)}},\,\,\hat{\theta}+\frac{z_{\alpha/2}}{\sqrt{n\mathcal{I}(\theta)}}\bigg{]}.$

We can use the bootstrap method to validate the results. Specifically, we
select a sample of size $n$ from a uniform distribution $\mathcal{U}[0,1]$,
say $(u_{1},u_{2},...,u_{n})$, to produce the set $x_{i}=F^{-1}(u_{i})$, where
$F(x)$ is defined in Eqn. (7). Its inverse is given by

$F^{-1}(u)=\sigma\,\sinh^{-1}(\text{tan}\big{[}(u-\frac{1}{2})\pi\big{]})+\theta.$
(30)

The bootstrapped distributions generated in this manner are shown in Fig. 5. A
number of different combinations of $\theta$ and $\sigma$ were selected. For
each case, estimates of the location parameter, $\hat{\theta}$, and the
variance term, $n\widehat{\text{Var}}(\hat{\theta})$, are obtained and listed
in the figure. The same results are also provided in Table I along with the
estimate of $\hat{\sigma}$. The estimates in all cases validate the equations
derived earlier (see Eqn. (29)).

Figure 5: Histogram of estimates of $\hat{\theta}$ and $n\widehat{\text{Var}}(\hat{\theta})$ from 10000 samples with $n=100$. TABLE I: Locationand Variance estimates using bootstrapping. | Location/Variance
---|---
Plot |  |  |  |  |
label | $\theta$ | $\sigma$ | $\hat{\theta}$ | $n\widehat{\text{Var}}(\hat{\theta})$ | $\hat{\sigma}$
(a) | 0.0 | 1.0 | 0.00 | 2.02 | 1.00
(b) | 0.0 | 3.0 | 0.00 | 19.11 | 3.09
(c) | 5.0 | 2.0 | 4.99 | 8.32 | 2.04
(d) | 5.0 | 3.0 | 4.99 | 18.96 | 3.08

##  VI Comparison with Other Distributions

Our next step is to compare the Cosh distribution with the Normal and Cauchy
distributions. We begin with a comparison with the standard Normal case. We
see in Fig. 6 that the Cosh distribution has heavier tails than the Gaussian
$\mathcal{N}(0,1)$ distribution. This is also illustrated in the Q-Q plot of
Fig. 7 for a random set of points from the Cosh distribution.

Figure 6: Comparison of standard Normal distribution with Cosh distribution.
Figure 7: Q-Q plot of Cosh against Normal distribution.

The second comparison is with the Cauchy distribution. One may recognize the
form of the pdf of the Cosh distribution as being similar to the Cauchy
distribution (Student t distribution with df=1) given by:

$f(x;\theta,\sigma)=\frac{1}{\pi\sigma\,(1+(\frac{x-\theta)}{\sigma})^{2})}$
(31)

with corresponding cdf given by

$F(x;\theta,\sigma)=\frac{1}{2}+\frac{1}{\pi}{\tan^{-1}\bigg{(}\frac{x-\theta}{\sigma}\bigg{)}}.$
(32)

We note here that for the Cauchy distribution

$\mathbb{E}[X]=\text{undefined}$

and

$\text{Var}(X)=\text{undefined}.$

This may be one of the reasons why the Cauchy distribution is not used in
practice. A graphical comparison of the Cosh and Cauchy distributions is
provided in Fig. 8.

Figure 8: Comparison of Cauchy distribution with Cosh distribution.

We can still develop an estimator based on the Cauchy distribution. Applying
MLE produces the following estimator:

$\hat{\theta}=\text{argmin}_{\theta\in\mathbb{R}}\sum_{i=1}^{n}\text{log}(1+(x_{i}-\theta)^{2})$
(33)

The associated loss function can be shown to be non-convex so it is not
suitable for regression. However, there are simple cases for which it will
produce an acceptable result.

Figure 9: $\rho(x)$ function. Figure 10: $\psi(x)$ function.

Figs. 9 and 10 show the characteristics of the $\rho(x)$ and $\psi(x)$
functions derived from the double exponential (L1), Gaussian (L2), Cosh and
Cauchy distributions. Comparing the different graphs, it is clear that the
Cosh characteristics are very similar to the L1 characteristics. It inherits
some of the good properties of L1 without the discontinuities in the first
derivative. This is the main reason for its recent rise in popularity in the
machine learning community.

##  VII Regression Examples

###  VII-A Location Problem

Consider a location problem of the form:

$y_{i}=\theta+\epsilon_{i},\quad i=1,\ldots,n.$

The one-dimensional data set under consideration, ranging from -2.8 to 20.5,
is given in Table II. For this set, the LS, Cauchy and log-cosh estimates for
$\theta$ were computed.

TABLE II: Location problem data set.

(1) | -2.80 | -1.98 | -1.70 | -1.20 | -1.10 | -0.82
---|---|---|---|---|---|---
(7) | -0.79 | -0.73 | -0.66 | -0.51 | -0.41 | -0.35
(13) | -0.23 | 0.10 | 0.22 | 0.25 | 0.37 | 0.52
(19) | 0.93 | 0.95 | 1.36 | 1.52 | 1.76 | 3.07
(25) | 20.50 | - | - | - | - | -

The results are provided in Table III and illustrated graphically in Fig. 11.
We see that the estimates of Cauchy and log-cosh are quite close in value
whereas the LSE is much further away due to the outlier on the right-hand
side. Asymptotically, the Cauchy estimator should produce the median, which in
this case is -0.23 while the mean is 0.73. We also note the standard errors
(s.e.’s), obtained using bootstrapping and given in Table III, indicate a 3
times larger value for LSE compared to the relatively smaller values for
Cauchy and log-cosh.

TABLE III: Comparison of LSE, Cosh and Cauchy for location problem $n=25$. Statistic | Technique Used
---|---
Used | LSE | Cosh | Cauchy
$\hat{\theta}$ | 0.73 | -0.06 | -0.19
s.e.($\hat{\theta}$) | 0.84 | 0.28 | 0.28
Figure 11: Estimates for Cauchy, Cosh and LSE on location problem.

###  VII-B Simple Linear Regression Example

We will now compare a number of different methods on a simple linear example
from the telephone data set in the Rfit package [14] in R. Rfit uses another
robust method called rank-based regression [15] [16] to be described shortly.
For our purposes, we will use the Belgium telephone data set [14][15]. This
data set is provided in Table IV and contains the number of calls (in units of
10’s of millions) made in Belgium in the years between 1950 - 1973. It has
$24$ data points, with $6$ outliers. The outliers are due to a change in
measurement technique without re-calibration for $6$ years, as is often cited.

TABLE IV: Belgium telephone data set.

x | 1950 | 1951 | 1952 | 1953 | 1954 | 1955 | 1956 | 1957
---|---|---|---|---|---|---|---|---
y | 0.44 | 0.47 | 0.47 | 0.59 | 0.66 | 9.73 | 0.81 | 0.88
x | 1958 | 1959 | 1960 | 1961 | 1962 | 1963 | 1964 | 1965
y | 1.06 | 1.2 | 1.35 | 1.49 | 1.61 | 2.12 | 11.9 | 12.4
x | 1966 | 1967 | 1968 | 1969 | 1970 | 1971 | 1972 | 1973
y | 14.2 | 15.9 | 18.2 | 21.2 | 4.3 | 2.4 | 2.7 | 2.9

Rank-based methods [15][16] also offer robust regression so it is useful to
compare log-cosh against this approach. The rank-based loss function is
governed by:

$\text{loss}_{R}(\beta)=\sum_{i=1}^{n}(y_{i}-x_{i}\beta)a_{n}(R_{n_{i}}(\beta))$
(34)

where $R_{n_{i}}(\beta)$ are the ranks of the residuals $(y_{i}-x_{i}\beta)$
and

$a_{n}(R_{n_{i}}(\beta))=\phi\left(\frac{R_{n_{i}}(\beta)}{n+1}\right)$ (35)

with $\phi(u)=2u-1$. Further details of these quantities may be found in
[15][16].

The results of the regression for LS, Huber, log-cosh, and rank are shown in
Fig. 12. The estimates are given in Table V. We note that log-cosh, Huber and
rank are all aligned whereas the LSE model is greatly affected by the
outliers. This example clearly illustrates the robustness property of the log-
cosh estimator. In addition, the Huber and rank results are very similar in
this case, while the log-cosh method produces a slightly higher slope. This
will vary from dataset to dataset.

Figure 12: Results for LSE, log-cosh, rank and Huber on a simple linear
problem. TABLE V: Slope and intercept for 4 regression methods on telephone
data set.

| $\beta_{1}$ | $\beta_{0}$
---|---|---
least squares | 0.504 | -983.9
log-cosh | 0.173 | -338.1
rank-based | 0.146 | -284.3
Huber ($\delta=0.1$) | 0.143 | -280.0

Figure 13: Komlogorov-Smirnov Goodness-of-Fit Test using residuals of
Telephone data set comparing p-values of Gaussian, Cauchy and Cosh
distributions.

A Goodness-of-fit test can be performed using the Komlogorov-Smirov (K-S)
method to examine what distribution best fits the residuals of the telephone
data set following least squares estimation (LSE). This test uses the
empirical cdf of the distributions. The results are shown graphically in Fig.
13 along with their p-values. We note that the null hypothesis is rejected in
the case of the Gaussian distribution at the $\alpha=0.05$ level whereas it
cannot be rejected for the Cauchy or Cosh distributions. Hence, one should use
a robust method. If the null hypothesis fails, one should consider a robust
method. On the other hand, the log-cosh loss provides equivalent estimates in
cases when outliers are not present in the data set. Therefore, it may be used
whether or not outliers exist.

##  VIII Multiple Linear Regression

The different approaches for robust regression can be compared in the context
of multiple linear regression to study standard errors of the estimates. We
use a well-known Swiss Fertility data set with 5 variables and 47
observations. This data set is part of the R environment. The data itself may
also be found in [15]. The variable abbreviations are as follows: A =
Agriculture, Ex = Examination, Ed = Education, C = Catholic, IM = Infant
Mortality. The goal is to build a model that predicts Fertility based on these
five explanatory variables. In Table VI, we provide the estimates for Huber,
Rank and log-cosh, along with their associated standard errors (s.e.). They
all produce similar s.e. values, although the results may vary slightly from
dataset to dataset.

TABLE VI: Estimates and standard errors $(s.e.)$ for Huber (H), Rank (R) and log-cosh (L) on Swiss data set. | $\hat{\beta}_{n}^{H}$ | s.e.($\hat{\beta}_{n}^{H}$) | $\hat{\beta}_{n}^{R}$ | s.e.($\hat{\beta}_{n}^{R}$) | $\hat{\beta}_{n}^{L}$ | s.e.($\hat{\beta}_{n}^{L}$)
---|---|---|---|---|---|---
A | -0.19 | 0.071 | -0.20 | 0.069 | -0.20 | 0.075
Ex | -0.28 | 0.258 | -0.25 | 0.249 | -0.26 | 0.264
Ed | -0.84 | 0.186 | -0.88 | 0.179 | -0.89 | 0.190
C | 0.10 | 0.035 | 0.10 | 0.034 | 0.10 | 0.035
IM | 1.21 | 0.388 | 1.19 | 0.375 | 1.40 | 0.395

##  IX Quantile Regression

Quantile regression is a robust method for studying the effect of explanatory
variables on the entire conditional distribution of the response variable
rather than just on the median. It has been used extensively since the initial
concepts were developed in the late 1970’s and early 1980’s [17]. It is based
on the so-called check function whereby the quantile of interest is set by a
parameter, $\tau$. The original check function exhibits a non-monotonic
behavior which has been the subject of a number of research papers over the
years [19][20][21][22][23]. In [12], we used the log-cosh function to develop
a new M-estimator to overcome the crossing problem. In this section, we derive
an MLE equivalent so that a statistical analysis can be carried out.

The original check function can be written concisely as follows:

$\rho_{\tau}(x_{i},\tau)=\begin{cases}-(1-\tau)x_{i}&\quad\text{if }x_{i}<0\\\
\tau x_{i}&\quad\text{if }x_{i}\geq 0\end{cases}$ (36)

where the parameter $\tau\in(0,1)$ is the quantile and $x_{i}$ is the $i$th
residual. Different regression quantiles represented by lines or hyperplanes
are obtained by selecting $\tau$ and minimizing the conditional quantile
function

$Q_{\tau}(Y|x,\tau)=\sum_{i=1}^{n}\rho_{\tau}(x_{i},\tau).$ (37)

A plot of $\rho_{\tau}(x_{i},\tau)$ for different values of $\tau$ is provided
in Fig. 14. Suppose we are interested in the median, which is also the 50th
percentile and the 2nd quartile.

Figure 14: Check function.

We would set $\tau=0.5$ and use the following check function,

$\rho_{\tau}(x_{i},0.5)=\begin{cases}-(0.5)x_{i}&\quad\text{if }x_{i}<0\\\
0.5x_{i}&\quad\text{if }x_{i}\geq 0.\end{cases}$ (38)

This is equivalent to the least absolute deviation (LAD) function or the L1
loss function. We see from the figure (black line, $\tau=0.5$) that it is
symmetric about $x_{i}=0$ but has a kink at 0. This is the same characteristic
that causes problems for the L1 function. Because of this kink, the derivative
of this function is discontinuous at 0. Other cases shown in the figure for
$\tau=0.1,0.25,0.75,0.9$ all exhibit a kink except that the function is
asymmetric in one direction or the other. The red lines associated with
$\tau=0.75$ are reminiscent of a check mark, and hence the name check
function. In any case, the kinks in this function imply discontinuous
derivatives leading to a number of mathematical and numerical problems when
solving for quantiles.

###  IX-A Continuous Check Function

To circumvent this problem, one could employ a continuous check function using
log-cosh as follows,

$\rho_{S}(x,\tau)=\text{log}(\text{cosh}(x))+(\tau-\frac{1}{2})x$ (39)

This function does not have any kinks as will be demonstrated shortly, but we
can already anticipate that it will be smoother simply due to the
characteristics of log-cosh.

The quantile regression problem involves minimizing an associated convex loss
function which is the conditional quantile function for each $\tau$ given by

$Q_{S}(Y|x,\tau)=\sum_{i=1}^{n}\rho_{S}(r_{i},\tau).$ (40)

To derive this continuous check function, we first postulate a pdf given by

$f(x)=\frac{e^{-(\tau-\frac{1}{2})x}}{\kappa\,\text{cosh}(x)}$ (41)

where $\kappa$ is a normalizing constant term such that the pdf integrates to
1. The value of $\kappa$ will vary as a function of the selected $\tau$. A
table of values for $\kappa$ for selected $\tau$ values is given in Table VII.
Note that for $\tau=0.5$, we obtain that $\kappa=\pi$ as expected, since it is
the original log-cosh pdf.

TABLE VII: Selected values of $\kappa$

$\tau$ | $\kappa$
---|---
0.0 | $\pi\sqrt{2}$
0.25 | $\pi\sqrt{4-2\sqrt{2}}$
0.5 | $\pi$
0.75 | $\pi\sqrt{4-2\sqrt{2}}$
1.0 | $\pi\sqrt{2}$

The MLE for this distribution can be derived as follows:

$L(x_{1},...,x_{n})=\bigg{(}\frac{e^{-(\tau-\frac{1}{2})x_{1}}}{\kappa\,\text{cosh}(x_{1})}\bigg{)}\dots\bigg{(}\frac{e^{-(\tau-\frac{1}{2})x_{n}}}{\kappa\,\text{cosh}(x_{n})}\bigg{)}.$
(42)

Then,

$-\ell(x_{1},...,x_{n})=\sum_{i=1}^{n}[\text{log(cosh}(x_{i}))+(\tau-\frac{1}{2})x_{i}+\text{log}(\kappa)]$
(43)

Therefore, after removing the constant term, we obtain Eqns. (39) and (40).
The Fisher information can be derived to be

$\mathcal{I}(\theta)=\frac{1}{2}.$ (44)

Figure 15: pdf of quantile function.

The distributions and loss functions can now be plotted. Starting with Eqn.
(41) and the normalizing constants $\kappa$ given in Table VII, different
values of $\tau$ were selected and their associated distributions plotted in
Fig. 15. The cases shown are for $\tau=0,0.25,0.5,0.75,1.0$. The distributions
shown on the right-hand side of the figure can be seen to skew from one side
to the other as $\tau$ decreases in value from 1 to 0. On the left-hand side,
we have plotted the loss term of the MLE in each case. We see that the use of
the logcosh function smooths out the kinks of the original loss function (see
Fig. 14).

A more general M-estimator based on log-cosh for quantile regression is given
by:

$\rho_{SMRQ}(x,\tau)=\frac{1}{2c}\text{log}(\text{cosh}(c(x-h)))+(\tau-s)x+v$
(45)

where $c$ is the desired curvature of the function (i.e. the severity of the
kink), $v$ controls vertical shift, $h$ controls the horizontal shift and
$\tau-s$ is used to produce the asymptotic slopes on the two sides of the
check function itself.

Quantile regression using the above equation is referred to as SMRQ for
‘smoother regression quantiles’. Fig. 16 shows five different cases of Eqn.
(45) by varying parameters $c$ and $v$ (while holding $s=0.5$ and $h=0$ fixed)
for $\tau=0.5$ and $\tau=0.7$. One can observe the different levels of
smoothness offered by the flexible check function, which can be varied easily
as the need arises.

We propose the use of $c=1/2$ and $v=1/2$ to produce:

$\rho_{SMRQ}(x,\tau)=\text{log}(\text{cosh}\big{(}\frac{x}{2}\big{)})+(\tau-\frac{1}{2})x+\frac{1}{2}.$
(46)

Figure 16: Loss functions for SMRQ for $\tau=0.5$ and $\tau=0.7$ cases.

The variance of the estimates due to this loss function cannot be obtained in
closed form. Therefore, a bootstrapping method can be used to extract standard
errors for estimates, which are naturally related to the variance.

###  IX-B Comparison with Convolutional Smoothing

Recently, a technique was published [24] that utilizes convolutional smoothing
[25] applied to the loss function and is implemented in the conquer library in
R. Convolutional smoothing provides a different approach to reducing the
crossing problem. Therefore, it is appropriate to compare it with SMRQ in
terms of solving the crossing problem.

One natural way to monitor the monotonicity of a quantile regression method is
to simply count the number of data points below a given regression line (or
hyperplane in the multivariate case). The median line (50th percentile) should
result in half the points below it. The 75th percentile should results in 75%
of the data below that line, and so on. One rule that should be enforced is
that a higher percentile will have more points below it than a lower
percentile. That is, we should not have more points below the 75th percentile
than we do at the 80th percentile. This would violate the basic definition of
percentiles, which is not desirable. However, in the original formulation of
quantile regression [18], this was often the case and was called the quantile
crossing problem.

To investigate the monotonicity of conquer and SMRQ, we selected the well-
known swiss data set in R, which has 5 variables. We show the monotonicity
results for conquer and SMRQ in Fig. 17. Note that the plot should be non-
decreasing if the crossing problem is avoided. This is true for SMRQ but
conquer exhibits 1 non-monotonic event.

Figure 17: Conquer vs. SMRQ on the swiss data set.

We now compare the standard errors (s.e.) of conquer and SMRQ in Table VIII.
Rather than listing each s.e. for all 5 variables in the 5 cases of $\tau$
(total of 25 values each), we compute the $L_{2}$-norm of the s.e.’s and list
those values in the table. The numbers represent an aggregate standard error
for each $\tau$ which is sufficient for this comparison.

We see that the errors are on average 10% smaller for SMRQ compared to
conquer. A similar result was also obtained on the diabetes dataset in R. In
addition, SMRQ possesses the monotonocity property for quartiles, deciles and
percentiles whereas conquer may exhibit some non-monotonic behavior. However,
conquer was intended to handle large-scale problems and requires less runtime
for the same size problem. This is a relatively small example for conquer, but
it is expected to perform much better with large data sets.

TABLE VIII: s.e. for conquest and SMRQ as a function of $\tau$ on swiss data set. Quantile | conquer (C) | SMRQ (S)
---|---|---
$\tau$ | $||(s.e._{C})||_{2}$ | $||(s.e._{S})||_{2}$
0.01 | 0.851 | 0.781
0.25 | 0.781 | 0.792
0.5 | 0.755 | 0.710
0.75 | 0.921 | 0.756
0.99 | 0.825 | 0.753

##  X Conclusions

In this paper, we identified the Cosh distribution from which the MLE for the
log-cosh function can be derived. We then derived the asymptotic variance,
asymptotic bias and confidence intervals. The log-cosh loss function was
compared to the Huber, rank and LSE loss functions in the case of simple
linear regression. The estimates and standard errors were found to be similar
for Huber, rank and log-cosh for multiple linear regression. Next, the use of
log-cosh in quantile regression was described to resolve the crossing problem,
the details of which can be found in [12]. The M-estimator for quantile
regression was compared and found to have a smaller standard error than
conquer, which uses convolutional smoothing. From the analysis provided
herein, it is clear that the log-cosh is an important loss function for
machine learning and is expected to increase in use over the coming years.

## References

  * [1] Grover, P. (2019, September 25). “5 Regression Loss Functions All Machine Learners Should Know”. Retrieved from https://heartbeat.fritz.ai/5-regression-loss-functions-all-machine-learners-should-know-4fb140e9d4b0.
  * [2] Gupta,S. (2022, April 14). “The 7 Most Common Machine Learning Loss Functions”. Retrieved from https://builtin.com/machine-learning/common-loss-functions.
  * [3] Wang, Q. et. al, “A Comprehensive Survey of Loss Functions in Machine Learning”, Annals of Data Science 9, April 2020.
  * [4] Jadon, S. “A survey of loss functions for semantic segmentation”, arXiv:2006.14822v4, Sept. 2020.
  * [5] R software, limma package.
  * [6] numpy documentation, https://numpy.org/doc.
  * [7] TensorFlow documentation regarding log-cosh loss, Retrieved from https://www.tensorflow.org/api-docs/python/tf/keras/losses/log-cosh.
  * [8] PyTorch documentation, https://pytorch.org/doc.
  * [9] Chen P., Chen G., Zhang S. “Log Hyperbolic Cosine Loss Improves Variational Auto-Encoder”, submitted to Sixth International Conference on Learning Representations, openreview.net, 2018.
  * [10] X. Xu, J. Li, Y. Yang and F. Shen, ”Toward Effective Intrusion Detection Using Log-Cosh Conditional Variational Autoencoder,” in IEEE Internet of Things Journal, vol. 8, no. 8, pp. 6187-6196, April 2021.
  * [11] Kalifi, E. Y. et. al, “Classification of Breast Cancer Lesions in Ultrasound Images by Using Attention Layer and Loss Ensemble in Deep Convolutional Neural Networks,” National Library of Medicine, published online Oct. 2021.
  * [12] Saleh, R., and Saleh, A.K.Md. E., “Solution to the Non-Monotonicity and Crossing Problems in Quantile Regression”, arXiv:2111.04805v2, Nov. 2021.
  * [13] Huber, P. J., Ronchetti, E. M., “Robust Statistics”, Second Edition, 2009.
  * [14] Kloke, J.D., McKean, J.W., “Rfit:Rank-based estimation for linear models”, The R Journal, 4(2), 57-64, 2012.
  * [15] Saleh, A.K.Md.E. et. al, “Rank-based Shrinkage and Selection with Application to Machine Learning”, John Wiley and Sons Publishers, Feb. 2022.
  * [16] Hettmansperger, T. P., McKean, J. W., “Robust Nonparametric Statistical Methods”, 2nd Ed. Chapman Hall, New York, 2011.
  * [17] Koenker, R., and Bassett, “Regression Quantiles,” Econometrica 46, 33-50 (1978).
  * [18] Koenker, R., Quantile Regression, Cambridge University Press, 2005.
  * [19] He, X., “Quantile Curves without Crossing”, The American Statistician 51, 186-192 (1997).
  * [20] Bondell, H.D., B.J. Reich and X. Wang, “Noncrossing quantile regression curve estimation”, Biometrika, 97: 825-838, 2010.
  * [21] Amerise, I., “Quantile Regression Estimation Using Non-Crossing Constraints”, Journal of Mathematics and Statistics, Volume 14: 107-118, 2018.
  * [22] Neocleous, T. and Portnoy, S. “On Monotonicity of Regression Quantile Functions”, Statistics and Probability Letters, 78: 1226-1229, 2007.
  * [23] Chernozhukov, V. and Fernandez-Val, I. and Galichon, A., “ Quantile and probability curves without crossing”, Econometrica 78, 1093-1125, 2010.
  * [24] He, X., Pan, X., Tan, K. M., and Zhou, W., “Smoothed Quantile Regression with Large-Scale Inference”, arXiv:202.05187v1, 2020.
  * [25] Fernandes, M. and Guerre, E. and Horta, E., “Smoothing quantile regressions”, arXiv:1905.08535v3, Aug. 2019.