# Self-Supervised RF Signal Representation Learning for NextG Signal
Classification with Deep Learning

Kemal Davaslioglu  University Technical Services, Inc. (UTS), Greenbelt, MD,
USA.  Serdar Boztaş  University Technical Services, Inc. (UTS), Greenbelt, MD,
USA.  Mehmet Can Ertem  University Technical Services, Inc. (UTS), Greenbelt,
MD, USA.  Yalin E. Sagduyu  Virginia Tech, National Security Institute,
Arlington, VA, USA.  Ender Ayanoglu  Center for Pervasive Communications and
Computing (CPCC), University of California, Irvine, Irvine CA, USA.

###### Abstract

Deep learning (DL) finds rich applications in the wireless domain to improve
spectrum awareness. Typically, DL models are either randomly initialized
following a statistical distribution or pretrained on tasks from other domains
in the form of transfer learning without accounting for the unique
characteristics of wireless signals. Self-supervised learning (SSL) enables
the learning of useful representations from Radio Frequency (RF) signals
themselves even when only limited training data samples with labels are
available. We present a self-supervised RF signal representation learning
method and apply it to the automatic modulation recognition (AMR) task by
specifically formulating a set of transformations to capture the wireless
signal characteristics. We show that the sample efficiency (the number of
labeled samples needed to achieve a certain performance) of AMR can be
significantly increased (almost an order of magnitude) by learning signal
representations with SSL. This translates to substantial time and cost
savings. Furthermore, SSL increases the model accuracy compared to the state-
of-the-art DL methods and maintains high accuracy when limited training data
is available.

###### Index Terms:

Automatic modulation recognition, wireless signal classification, contrastive
learning, deep learning, self-supervised learning, spectrum awareness.

##  I Introduction

_Deep learning_ (DL) finds diverse applications in wireless communications
such as in wireless signal classification, spectrum sensing, spectrum sharing,
resource allocation and provisioning, and attack detection and mitigation [1].
In particular, DL has emerged as a key enabler for _next-generation (NextG)
communications_ by providing novel means to solve complex tasks [2]. One
particular task that benefits from DL is _wireless signal classification_ that
serves a variety of applications ranging from user equipment (UE)
identification and PHY-layer authentication to jammer detection. To that end,
_automatic modulation recognition_ (AMR) with DL has drawn major attention to
showcase the use of deep neural networks (DNNs) in capturing the intrinsic
characteristics of wireless signals especially in the low signal-to-noise
ratio (SNR) regime, where statistical and conventional machine learning
methods with hand-crafted features (such as high-order cumulants as in [3])
often fail [4, 5, 6, 3].

Despite the unique performance promises, DL in Radio Frequency (RF)
applications faces several challenges. First, the received signal is subject
to noise, channel, traffic, and interference effects. Modeling or estimating
these effects individually or determining which signal features to use and how
to preprocess them is a very difficult task. For example, a recent study in
[3] transformed the complex-valued signal waveforms into images based on the
density of data features to express its statistical information. Second,
collecting representative RF data is generally an expensive and time consuming
endeavor. An example is dynamic spectrum sharing (DSA) as in the Citizens
Broadband Radio Service (CBRS) band, where the 5G communications need to share
the spectrum with incumbent users (such as radar) by classifying the received
signals and opportunistically accessing the spectrum when idle. Thus, it is
essential for wireless systems to train classifiers even when limited training
data samples are available.

Currently, a popular practice in using the DL models for spectrum
classification tasks (e.g., spectrum sensing and waveform identification) is
to employ _transfer learning_ where the DNN weights are initialized based on a
computer vision task trained using a vision dataset [6] or they are randomly
initialized [7, 8, 9]. However, these DL models do not exploit the unique
characteristics of RF spectrum data (e.g., channel distortions, phase and
frequency offsets, and noise effects) as they do not use the signal
representations during pretraining.

_Self-supervised representation learning_ from unlabeled data has recently
attracted major attention in the natural language processing (NLP) and
computer vision areas [10, 11, 12, 13, 14, 15]. The representations learned
using self-supervised learning (SSL) have been shown to outperform their
supervised counterparts [13]. One of the main goals in SSL is to learn
representations and pretrain models that can be transferred to _downstream
tasks_ by fine-tuning. In computer vision, examples of downstream tasks are
detection and segmentation. In wireless communication applications, the
downstream tasks may include determining the presence of a signal and
identifying its emitter, protocol, and modulation. In radar applications, the
downstream tasks may include resolving Pulse Descriptor Words (PDWs) such as
pulse repetition patterns, pulse width, SNR, duty cycle, and time, frequency
and direction of arrival.

SSL methods typically require _a pretext task_ to learn from and _a loss
function_ to optimize [13]. The pretext task is solved to learn a good
representation of the data, but typically is not the true purpose of learning.
For example, predicting the rotation of an image [10] or predicting the center
word given the surrounding words [11] are some commonly used pretexts. To the
best of our knowledge, there are only two prior studies that have applied SSL
to AMR task [16, 17]. Both of these references are based on the SimCLR method
[18], in which [17] uses adding Gaussian noise and carrier frequency offset
(CFO) for augmentation, while [16] uses rotations. In comparison to these
studies, we present an SSL method that is based on the MoCo-v3 framework. In
addition, we include five AMR-related transformations for the data
augmentation step. In Section IV, we compare the performance of our approach
with [16, 17] and demonstrate substantial gains.

In this letter, we study AMR as an example of wireless signal classification
using DL. Our contributions are twofold: (i) We propose an efficient self-
supervised signal representation method for AMR to learn the signal
representations in RF applications and reduce the need for labeled data. (ii)
We propose a set of data augmentation transformations that do not alter the
semantic information of the data. These transformations may significantly
change the amplitude or phase of the signal, but they preserve the semantic
meaning (e.g., modulation type, radar signal type, and emitter information).
The DL model is pretrained with contrastive learning using the proposed
transformations _without the need for any labels_. Then, the AMR task uses the
backbone of the pretrained model and fine-tunes its weights using labeled
data. The self-supervised pretraining significantly improves the sample
efficiency of the learning process, which is defined as the number of labeled
samples required to achieve a certain performance (in our case, classification
accuracy). Furthermore, SSL increases the accuracy of the DL model and
surpasses the performance of the existing supervised learning benchmarks.

The remainder of the letter is organized as follows: Section II defines the
AMR problem and data augmentations. Section III describes the SSL approach.
Section IV evaluates the performance and Section V concludes the letter.

##  II System Model

Suppose that a receiver receives the following signal with a single antenna
and over a single channel:

$\displaystyle y(i)=A(i)e^{j(\omega i+\phi)}x(i)+n(i)$ (1)

for $i=1,\cdots,N$, where $x(i)$, $y(i)$, and $n(i)$ are the $i$th transmitted
signal, received signal, and noise samples, respectively, $A(i)$ is the
channel gain, $\omega$ is the frequency offset, and $\phi$ is the phase
offset. The received signal can be expressed in a vector form
$\textbf{y}=[y_{1},\cdots,y_{N}]$, where all $N$ signal samples are modulated
using the same modulation.

The goal of the AMR classifier is to train a model that learns to map the
received signal sample y to a modulation class, that is,
$f_{\bm{\theta}}(\textbf{y}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{c}$, where
parameters (weights and biases) of the neural network are represented by
$\bm{\theta}$. The dimension of the received signal is $d$ and there are $c$
classes (labels) such that when a digitized RF waveform y is received and
input to the modulation recognition model, a label is returned, i.e.,
$f_{\bm{\theta}}(\textbf{y})\in\mathbb{R}^{c}$. This problem can be formulated
as a multiclass classification problem which can be solved by minimizing the
cross-entropy loss that is given by

$\displaystyle\ell_{CE}(\textbf{y})=-\sum_{i=0}^{c}\beta_{i}\log\left(p_{i}\right),$
(2)

where $\beta_{i}$ is a binary indicator that specifies if the vector y belongs
to class $i$ or not, that is, $\beta_{i}=1$ if y belongs to class $i$ and is
$0$ otherwise, and $p_{i}$ is the output of the neural network that denotes
the predicted probability. The cross-entropy loss is minimized to train the
DNN with supervision.

Figure 1: RF waveform augmentations. TABLE I: Parameter ranges used to generate waveform transformations.  Transformation | Min | Max
---|---|---
DC shift | 0 | 0.0001
Time shift (samples) | -40 | 40
Amplitude scale | 0.8 | 1.2
Zero-masking (samples) | 0 | 25
Additive White Gaussian noise | $\mathcal{N}(0,10^{-5})$

###  II-A Data augmentation methods

A key component of _contrastive learning_ is a set of _data augmentation
methods_ that do not alter the semantic information of the data. Most
contrastive learning algorithms learn representations that are maximally
similar for augmented instances of the same sample and minimally similar for
those from different samples. In computer vision, commonly used
transformations include random resized cropping, horizontal flipping, color
jittering, grayscale conversion, blurring, and solarization [15]. However,
most of these transformations cannot be applied to the RF signals. In this
letter, we carefully select five transformations to augment the RF signals
during the self-supervised training stage while maintaining the semantic
information of the signals. These transformations are adding a DC shift,
shifting the signal in time domain (time shift), amplitude scaling
(multiplying with a constant factor), zero-masking (nulling a set of
consecutive samples), and adding additive white Gaussian noise (AWGN).
Although these transformations can significantly alter the numerical values,
they preserve the semantic information of the data (e.g., constellation map
regions and frame structure). Fig. 1 illustrates an example of each
transformation where the original signal and its augmented version are shown
together. The x-axis denotes the time and y-axis denotes the amplitude of
signal. Table I presents the parameter ranges of these transformations. Each
transformation is applied with a random value within the specified range,
which depend on the data. For example, as the values in the RML2016.10a [7]
dataset are very small, we add a zero-mean noise with a small variance.

Figure 2: Contrastive learning using MoCo-v3.

##  III Self-Supervised Learning

Contrastive methods learn feature representations by enforcing similar
features to be close to each other (_positive samples_) and dissimilar
features to be far away from each other (_negative samples_). Fig. 2 shows the
MoCo-v3 framework [15] that uses two encoders and the pseudocode of the
MoCo-v3 algorithm is presented in Algorithm 1. The first encoder
$\textbf{f}_{q}$ includes a backbone (e.g., ResNet50 model [19]), a projection
head, and a prediction head. The backbone consists of a convolutional neural
network (CNN), whereas the projection and prediction heads are made up of
multi-layer perceptrons (MLPs), i.e., feedforward neural networks. The
projection and prediction heads consist of 3-layer and 2-layer MLPs,
respectively, where each linear layer (except the last one) is followed by a
batch normalization and a Rectifying Linear Unit (ReLU) activation. The second
encoder $\textbf{f}_{k}$ only has a backbone and a projection head. We denote
the output vectors of the encoders $\textbf{f}_{q}$ and $\textbf{f}_{k}$ as q
and k, respectively, both have 256 dimensions. MoCo-v3 uses the Information
Noise-Contrastive Estimation (InfoNCE) loss [20] that is given by

$\displaystyle\ell_{\textbf{q}}=-\log\left(\frac{\exp\left(\textbf{q}\textbf{k}_{+}/\tau\right)}{\exp\left(\textbf{q}\textbf{k}_{+}/\tau\right)+\sum_{\textbf{k}_{-}}\exp\left(\textbf{q}\textbf{k}_{-}/\tau\right)}\right),$
(3)

where $\tau$ is a temperature hyper-parameter. Each RF waveform in a batch is
augmented twice using the RF transformation function $A(\cdot)$ that applies
all the data transformations discussed in Table I consecutively. The output
vectors of the same waveform form a positive pair, denoted by
$\textbf{k}_{+}$, and others form the negative pairs, denoted by
$\textbf{k}_{-}$. In each batch, we calculate the InfoNCE loss in (3) and use
the queries q to retrieve the corresponding keys k. The backpropagation loss
is calculated and the query encoder $\textbf{f}_{q}$ is updated. The second
encoder $\textbf{f}_{k}$ is updated by the moving average of the query encoder
$\textbf{f}_{q}$ to maintain consistency.

Algorithm 1 MoCo-v3 algorithm for SSL.

1:procedure MoCo-v3($x$)

2: Input: Initial momentum encoder $\textbf{f}_{k}$, initial query encoder
$\textbf{f}_{q}$, RF transformation function $A(\cdot)$, momentum update
coefficient $\alpha$, and temperature coefficient $\tau$.

3: for $\textbf{x}\in\mathcal{X}$ do $\triangleright$ Load a batch of $N$
samples

4: $\textbf{x}_{1}=A(\textbf{x})$ and $\textbf{x}_{2}=A(\textbf{x})$
$\triangleright$ Augment x twice.

5: $\textbf{q}_{1}=f_{q}(\textbf{x}_{1})$ and
$\textbf{q}_{2}=f_{q}(\textbf{x}_{2})$ $\triangleright$ Forward pass
$\textbf{x}_{1}$ and $\textbf{x}_{2}$ through the query encoder.

6: $\textbf{k}_{1}=f_{k}(\textbf{x}_{1})$ and
$\textbf{k}_{2}=f_{k}(\textbf{x}_{2})$ $\triangleright$ Forward pass
$\textbf{x}_{1}$ and $\textbf{x}_{2}$ through the momentum encoder.

7: Calculate the InfoNCE loss in (3).

8: Back-propagate the loss.

9: Update the query encoder $\textbf{f}_{q}$.

10: Update the momentum encoder $\textbf{f}_{k}$ using the momentum by
$\textbf{f}_{k}=\alpha\textbf{f}_{k}+(1-\alpha)\textbf{f}_{q}$.

11: end for

12:end procedure

In comparison with MoCo-v1 [13] and MoCo-v2 [14], MoCo-v3 algorithm [15]
includes slight changes to provide simplicity, accuracy, and stability
improvements. It does not require the use of an external memory. Instead, it
uses the encoded features in the same batch. MoCo-v3 has an architectural
change where the query encoder $\textbf{f}_{q}$ uses an extra prediction head.
The MoCo-v3 algorithm was originally proposed for the computer vision domain
[15]. Our proposed algorithm builds on the MoCo-v3 framework and extends it to
RF applications with two subtle differences. First, the set of transformations
presented in this letter for data augmentation are different. Second, the
neural network architectures of the query and key encoders in this letter are
slightly different (see Section IV for details).

Overall, AMR model training consists of two steps. In the first step, the
signal representations are learned using the MoCo-v3 algorithm. This step is
called as _pretraining_. In the second step, the backbone of the encoder is
used to initialize the model. Its last layer is changed to have consistent
size with the AMR problem. The model is trained end-to-end until convergence,
where the weights and biases of the backbone are also updated. This step is
called as _end-to-end fine-tuning_. Note that there is also an alternate
approach called _linear probing_ [15], where the parameters of the backbone
are frozen and only the last layers of the classifier are updated. In this
letter, results with linear probing are omitted since its performance is
consistently lower than end-to-end fine-tuning.

##  IV Performance Evaluation for AMR

Performance is evaluated using the RML2016.10a dataset [7]. This dataset
includes 11 modulations (classes). Signals in the dataset cover an SNR range
from -20 dB to 18 dB in 2 dB increments. At each SNR, there are 1000 samples
per modulation type. The dataset is not readily split into training and test
subsets and we observe that the performance significantly depends on the way
data is split. In all simulations, we use the same training, validation, and
test datasets for repeatability. We leave 20% of the data for test and the
model is not evaluated on the test dataset during training. For the SSL stage,
the model is trained on the remaining 80% of the data without any labels. For
the supervised learning stage, we vary the training set portion of the
training and validation sets (constituting 80% data) from 0.1% to 90%, while
the remaining data is used for validation. After the training stage, all
models are evaluated on the test dataset using the model with the best
validation loss. The algorithms are implemented in Python using the PyTorch
framework. The experiments are conducted using a server with an NVIDIA GeForce
RTX 3090 Ti Graphics Processing Unit (GPU). For the SSL code, we modified the
code released by Facebook [15]. We used the temperature parameter as
$\tau=1.0$ as suggested in [15]. The models are trained for 100 epochs during
the SSL stage.

Our AMR model modifies the ResNet50 architecture [19] to be compatible with
the RF data. The first layer of convolutions is changed from 3 channels to 1
channel. We kept the kernel size as (7, 7) and stride size as (2, 2) as in the
original model. The output layer is changed to have 11 neurons representing
one for each label. The RF signals consist of 128 symbols in the RML2016.10a
dataset [7]. We reshape the in-phase and quadrature components (I/Q) of the RF
signal to $1\times 2\times 128$ using the channel-first notation. AdamW [21]
optimizer is used with a learning rate of 0.01. The learning rate in the
supervised learning is reduced to its half if the validation loss does not
improve for five consecutive epochs. All simulations are run for 500 epochs.
If the validation loss does not improve for 20 epochs, the simulations are
early terminated.

For comparison, we consider two supervised benchmarks using the same ResNet50
model. The first model uses transfer learning where the model is pretrained on
the ImageNet dataset [22]. The second one uses the Xavier initialization that
draws from a uniform distribution in $[-1,1]$ and then scales the weights by
$\sqrt{6/(m_{i}+m_{i+1})}$, where $m_{i}$ is the number of neurons at layer
$i$. We also include the Multi-Channel Convolutional Long short-term Deep
Neural Network (MCLDNN) [8] and Parameter Estimation and Transformation based
CNN-GRU Deep Neural Network (PET-CGDNN) models [9] in our comparison. These
models currently achieve the best performance reported in the literature for
this dataset. In addition, we implemented two SSL methods based on the SimCLR
method for comparison, namely, SimCLR-AWGN [17] that uses the AWGN
transformation and SimCLR-5TX that uses the five transformations considered in
this letter.

TABLE II: The AMR accuracy for different backbone initializations.

| Percent data used in training
---|---
Method | 0.5% | 1% | 5% | 10% | 50% | 75% | 90%
Xavier init. | 9.1 | 11.4 | 49.2 | 53.6 | 58.6 | - | 61.2
Pretrained | 14.9 | 34.5 | 53.1 | 54.4 | 60.0 | 61.0 | 61.1
MoCo-v3-512 (ours) | 50.4 | 53.1 | 55.2 | 54.6 | 61.4 | 61.8 | 62.4
MoCo-v3-1024 (ours) | 49.5 | 52.2 | 52.9 | 55.3 | 61.3 | 61.9 | 62.4
MoCo-v3-4096 (ours) | 45.9 | 49.0 | 53.1 | 55.2 | 60.9 | 62.2 | 62.6
SimCLR-5TX | 39.1 | 42.4 | 52.3 | 53.3 | 57.9 | 60.3 | 60.9
SimCLR-AWGN [17] | 26.9 | 34.8 | 50.2 | 53.1 | 57.7 | 60.5 | 60.4
MCLDNN [8] | - | - | - | - | - | 62.1 | -
PET-CGDNN [9] | - | - | - | - | - | 60.4 | -
# Training samples | 880 | 1760 | 8.8K | 17.6K | 88K | 132K | 158.4K

Table II studies the effect of different initializations. We observe that the
Xavier initialized model performs the worst across all training and validation
splits. The proposed MoCo-v3 pretrained model achieves the best performance in
all splits. When only 0.5% of training data is used over all SNR values, the
Xavier initialized model and the pretrained model achieve only 9.1% and 14.9%
accuracy (note that since there are 11 labels, purely random classification
would achieve $\frac{100}{11}$ $\approx$ 9.1% accuracy), whereas the proposed
model with projection size of 512 (MoCo-v3-512) achieves 50.4% accuracy
demonstrating a 35-41% improvement. The pretrained model achieves 53.1% only
when 5% of the data used. This indicates a 10$\times$ improvement in sample
efficiency which means it takes 10 times more labeled samples for the
pretrained model to achieve the same performance of MoCo-v3 pretrained model
in the low-labeled data regime. As the number of labeled samples increases,
the accuracy gap between the pretrained and self-supervised models decreases.
The accuracy difference between the pretrained model and MoCo-v3-512 starts at
35% at 0.5% data and drops to 1% improvement when 90% of the data is used,
since the importance of pretraining diminishes as the model is trained with
more labels. When we compare the performance of different SSL methods, MoCo-v3
significantly outperforms the SimCLR variants. For example, SimCLR-AWGN and
SimCLR-5TX achieve 26.9% and 39.1% accuracy at 0.5% data, respectively, while
MoCo-v3-512 achieves 50.4% indicating a significant difference between the
SimCLR and proposed algorithms. This result also indicates the importance of
transformations such that the SimCLR-5TX model consistently performs better
than the SimCLR-AWGN model.

Table III shows the effects of projection head size for different training and
validation splits. Projection head uses the features extracted by the backbone
CNN and improves its representation quality [12]. The projection head size
determines the width of this MLP. The models are pretrained using different
projection head sizes ranging from 128 to 4096. After pretraining, the models
are end-to-end trained (with labels) using different train and validation
splits. We observe that smaller projection head sizes generally perform better
in the low-data regime, while larger projection heads have better performance
for high-data regime. This is probably due to the fact that a large projection
head has more trainable parameters, which requires more samples to be trained
on and in the low-data regime, there are not simply enough samples.

Figs. 3(a)-(f) show the AMR accuracy (averaged over all labels) versus SNR for
different splits of training and validation data. For all splits, we observe
that as the SNR decreases below -4 dB, all models suffer in making accurate
classifications. This is an expected result since the noise power overwhelms
the signal strength and the accuracy decreases. When 0.5% data is used, the
self-supervised classifier achieves more than 60% accuracy for SNR$>$$-2$ dB
whereas the pretrained model achieves at most 20%. Similarly, when 1% of the
labeled data is used, the Xavier initialized and pretrained models achieve at
least 9% and 50% accuracy at SNR$>$$2$ dB, respectively, while all MoCo-v3
pretrained models achieve more than 70% accuracy. As the number of labeled
data increases, all models perform very similar to each other. However, self-
supervised representation learning still provides advantages. As an example,
when 90% of the data is used, the pretrained model achieves 91.3% accuracy at
18 dB SNR, whereas the MoCo-v3-512 model achieves 92.4%.

TABLE III: Effect of projection head size in MoCo-v3 on the AMR accuracy.

Projection Size | 0.5% | 1% | 5% | 10% | 50% | 75% | 90%
---|---|---|---|---|---|---|---
256 | 49.5 | 53.5 | 54.2 | 54.6 | 61.2 | 61.3 | 62.3
512 | 50.4 | 53.1 | 55.2 | 54.6 | 61.4 | 61.8 | 62.4
1024 | 49.5 | 52.2 | 52.9 | 55.3 | 61.3 | 61.9 | 62.4
2048 | 47.6 | 49.8 | 54.7 | 54.2 | 60.7 | 61.7 | 62.4
4096 | 45.9 | 49.0 | 53.1 | 55.2 | 60.9 | 62.2 | 62.6

Figure 3: Accuracy vs. SNR for different training and validation splits.

As a final note, we compare the size and computational complexity of the
models. The ResNet-50 model has 23.52M parameters, whereas MCLDNN and PETCGDNN
have 406.2K and 71.9K, respectively. The computational complexity of the
ResNet-50 model is 293.4M floating-point operations (flops), which is
substantially higher than the MCLDNN with 0.7M flops and PETCGDNN with 1.3M
flops. Despite the additional complexity of ResNet-50, we note that our main
contribution is to use the SSL method to improve the AMR performance rather
than developing a new architecture.

##  V Conclusion

In this letter, we studied the AMR problem using the DL models. The AMR models
using SSL can learn signal representations without any annotation. This
property can be leveraged to increase the model performance and sample
efficiency. For this goal, we proposed a contrastive learning algorithm using
the MoCo-v3 approach. In addition, we proposed five data augmentation methods
to learn the signal representations. The learned representations increase the
sample efficiency, which means that the model can be fine-tuned with less
labeled data to achieve high performance. On the other hand, when there is a
lot of labeled data, the learned representations are still able to improve the
classification accuracy. In future work, we plan to evaluate the performance
of other SSL algorithms and additional signal transformations. In addition, we
plan to extend the self-supervised signal representation learning framework to
other RF tasks such as emitter identification and radar waveform recognition.

## References

  * [1] T. Erpek, T. J. O’Shea, Y. E. Sagduyu, Y. Shi, and T. C. Clancy, “Deep learning for wireless communications,” in _Development and Analysis of Deep Learning Architectures_. Springer, 2020, pp. 223–266.
  * [2] E. Ayanoglu, K. Davaslioglu, and Y. E. Sagduyu, “Machine learning in NextG networks via generative adversarial networks,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 8, no. 2, pp. 480–501, 2022.
  * [3] Y. Lin, Y. Tu, Z. Dou, L. Chen, and S. Mao, “Contour stella image and deep learning for signal recognition in the physical layer,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 7, no. 1, pp. 34–46, 2021.
  * [4] Y. Shi, K. Davaslioglu, Y. E. Sagduyu, W. C. Headley, M. Fowler, and G. Green, “Deep learning for RF signal classification in unknown and dynamic spectrum environments,” in _IEEE International Symposium on Dynamic Spectrum Access Networks (DySPAN)_ , 2019.
  * [5] S. Soltani, Y. E. Sagduyu, R. Hasan, K. Davaslioglu, H. Deng, and T. Erpek, “Real-time and embedded deep learning on FPGA for RF signal classification,” in _IEEE Military Communications Conf._ , 2019.
  * [6] S.-H. Kim, J.-W. Kim, W.-P. Nwadiugwu, and D.-S. Kim, “Deep learning-based robust automatic modulation classification for cognitive radio networks,” _IEEE Access_ , vol. 9, pp. 92 386–92 393, 2021.
  * [7] T. J. O’Shea and N. West, “Radio machine learning dataset generation with GNU radio,” _Proceedings of the 6th GNU Radio Conference_ , 2016.
  * [8] J. Xu, C. Luo, G. Parr, and Y. Luo, “A spatiotemporal multi-channel learning framework for automatic modulation recognition,” _IEEE Wireless Communications Letters_ , vol. 9, no. 10, pp. 1629–1632, 2020.
  * [9] F. Zhang, C. Luo, J. Xu, and Y. Luo, “An efficient deep learning model for automatic modulation recognition based on parameter estimation and transformation,” _IEEE Communications Letters_ , vol. 25, no. 10, pp. 3287–3290, 2021.
  * [10] S. Gidaris, P. Singh, and N. Komodakis, “Unsupervised representation learning by predicting image rotations,” in _International Conference on Learning Representations (ICLR)_ , 2018.
  * [11] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” _arXiv:1301.3781_ , 2013.
  * [12] T. Chen, S. Kornblith, K. Swersky, M. Norouzi, and G. Hinton, “Big self-supervised models are strong semi-supervised learners,” _arXiv preprint arXiv:2006.10029_ , 2020.
  * [13] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 9726–9735.
  * [14] X. Chen, H. Fan, R. Girshick, and K. He, “Improved baselines with momentum contrastive learning,” _arXiv:2003.04297_ , 2020.
  * [15] X. Chen, S. Xie, and K. He, “An empirical study of training self-supervised vision transformers,” in _IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2021. [Online]. Available: https://github.com/facebookresearch/moco-v3
  * [16] D. Liu, P. Wang, T. Wang, and T. Abdelzaher, “Self-contrastive learning based semi-supervised radio modulation classification,” in _IEEE Military Communications Conf._ , 2021, pp. 777–782.
  * [17] X. Yun and X. Zhou, “Exploring self-supervised learning for radio signal recognition,” in _IEEE Int. Conf. on High Performance Computing & Communications_, 2021, pp. 2425–2430.
  * [18] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” _arXiv preprint arXiv:2002.05709_ , 2020.
  * [19] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _arXiv preprint arXiv:1512.03385_ , 2015.
  * [20] A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _arXiv preprint arXiv:1807.03748_ , 2018\.
  * [21] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” in _International Conference on Learning Representations (ICLR)_ , 2019.
  * [22] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “ImageNet: A large-scale hierarchical image database,” in _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2009, pp. 248–255.


# Incremental cutting-plane method and its application

Nagisa Sugishita, Andreas Grothey††footnotemark: , Ken McKinnon††footnotemark:
School of Mathematics, University of Edinburgh, James Clerk Maxwell Building,
Edinburgh (UK), EH9 3FD, Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider regularized cutting-plane methods to minimize a convex function
that is the sum of a large number of component functions. One important
example is the dual problem obtained from Lagrangian relaxation on a
decomposable problem. In this paper, we focus on an incremental variant of the
regularized cutting-plane methods, which only evaluates a subset of the
component functions in each iteration. We first consider a limited-memory
setup where the method deletes cuts after a finite number of iterations. The
convergence properties of the limited-memory methods are studied under various
conditions on regularization. We then provide numerical experiments where the
incremental method is applied to the dual problems derived from large-scale
unit commitment problems. In many settings, the incremental method is able to
find a solution of high precision in a shorter time than the non-incremental
method.

##  1 Introduction

The focus of this paper is a convex optimization problem whose objective
function is the sum of a large number of convex component functions. One
notable example is the dual problem derived from a decomposable optimization
problem.

Cutting-plane methods are popular solution methods to solve such a convex
optimization problem [13, 10]. Given the values and subgradients of the
component functions evaluated at the previous points, the methods construct a
piecewise affine model which minorizes the objective function. The minimizer
of the model is then used as the next point. Namely, for each component
function, its value and subgradient are evaluated at this new point and the
methods compute an affine function which minorizes the component. This affine
function is referred to as a cut and used to update the aforementioned
piecewise affine model. The modified model is then minimized and the above
procedure is repeated. The plain cutting-plane method suffers instability and
often regularization is added to the model to mitigate the issue [9].

In this paper, we consider an incremental variant of the regularized cutting-
plane methods in a sense that it only evaluates a subset of the components in
each iteration. The study on the incremental variant of the cutting-plane
methods is motivated by the success of other incremental first-order methods
such as incremental gradient descent methods and incremental subgradient
methods [14, 5]. One of the most notable examples is the training of neural
networks with stochastic gradient methods [6], which is closely related to
incremental gradient descent methods. Empirically it has been observed that
incremental methods often converge much faster than non-incremental methods,
especially when the iterate is far from the solution. It is also expected to
be advantageous for the cutting-plane method to evaluate only a subset of the
components and update the models and the iterate more frequently if the number
of components is large.

In a typical approach, the piecewise affine model described above grows
unlimitedly as the algorithm proceeds. That is, the model uses data obtained
in all of the previous iterations and its size increases every time a new
point is obtained. When the model gets large, it is of interest to delete some
data (i.e. delete some cuts) and keep the model size moderate. Heuristics may
be used to select cuts to be deleted, but they are not supported by any
convergence results [7, Remark 9.8]. We study the behaviour of the regularized
cutting-plane method in the limited-memory setup under various conditions on
the regularization.

Note that similar ideas have been explored in the context of the bundle
methods [12]. The bundle methods also construct a piecewise affine model of
the objective function and solve the model with regularization. Typically the
bundle methods support the limited-memory setting thanks to the technique
called cut aggregation. Recently, extensions of the methods to an incremental
setup have been studied [1, 11]. However, the current approaches often require
stronger assumptions (e.g., knowledge on the Lipschitz constance of the
objective) and further development is to be seen. See [1] for further
discussion on recent advances in the incremental bundle methods.

The rest of the paper is structured as follows. Section 2 introduces the
incremental cutting-plane method and studies its convergence property under a
limited-memory assumption. In section 3 we present some computational results.
Finally, in Section 4 conclusions of this paper are discussed.

##  2 Cutting-plane methods

In this section, we study the following convex optimization problem:

$\displaystyle\min\ $ $\displaystyle f(x)=\sum_{i=1}^{m}f_{i}(x)$ (2.1) s.t.
$\displaystyle x\in X,$

where we assume the following.

######  Assumption 2.1.

  1. 1.

$X\subset\mathbb{R}^{n}$ is a nonempty, closed and convex set.

  2. 2.

$f_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a convex function and is finite
and subdifferentiable over $X$ for $i=1,2,\ldots,m$.

  3. 3.

There exists $C>0$ for which

$\|g_{i}(x)\|\leq C,\qquad\forall i,x\in X,g_{i}(x)\in\partial f_{i}(x),$

where $\partial f_{i}(x)$ is the subdifferential of $f_{i}$ at $x$ for each
$i$.

This is a standard assumption to study incremental subgradient methods [14].
We note that for each $i$, given the convexity of $f_{i}$, the boundedness of
subgradients is equivalent to the $C$-Lipschitz continuity of $f_{i}$ [16,
Lemma 2.6]. We write $f^{*}=\inf_{x\in X}f(x)\in\mathbb{R}\cup\\{-\infty\\}$
and $X^{*}=\\{x\in X\mid f(x)=f^{*}\\}$, which may be empty.

###  2.1 Algorithm

The cutting-plane method is an iterative solution method which computes
candidate points $\\{x_{k}\\}$ to solve (2.1) approximately. In the $k$th
iteration, given the previous points $x_{0},x_{1},\ldots,x_{k}$, the method
computes the next point $x_{k+1}$ by minimizing a regularized model

$\displaystyle
x_{k+1}=\arg\min\sum_{i=1}^{m}\hat{f}_{k,i}(x)+\frac{1}{2t_{k}}\|x-x_{k}\|^{2},$

where $\hat{f}_{k,i}$ is a piecewise affine model of the $i$th component

$\displaystyle\hat{f}_{k,i}(x)=\max_{l=0,1,\ldots,k}f_{i}(x_{k-l})+g_{i}(x_{k-l})^{T}(x-x_{k-l}),$
(2.2)

$f_{i}(x_{l})$ and $g_{i}(x_{l})$ are the value and any subgradient of $f_{i}$
at $x_{l}$ ($l=0,1,\ldots,k$) and $t_{k}$ is a parameter to adjust the
strength of the regularization. We refer to $t_{k}$ as a step size in a reason
described later. The values and subgradients of the components are evaluated
at the new point $x_{k+1}$ to update the model (2.2) and the above process is
repeated.

In this paper, we consider an incremental variant of the method in a sense
that it only evaluates a subset of the components $f_{i}$ in each iteration.
We denote the index set of components evaluated in the $k$th iteration by
$I_{k}$. Since there are $m$ subproblems,

$I_{k}\subset\\{1,2,\ldots,m\\},\qquad\forall k.$

We note that the standard cutting-plane method is a special case where
$I_{k}=\\{1,2,\ldots,m\\}$ for all $k$. We refer to this non-incremental
method as the full-step method.

Furthermore, it is of interest to bound the memory requirement of the method.
To this end, we consider a method which only uses information obtained in the
last $W$ iterations, where $W$ is a prescribed number. That is, the piecewise
affine models only use cuts obtained within the last $W$ iterations.

The algorithm is formally described below.

######  Algorithm 2.2.

  1. 1.

Select $x_{0}\in X$, memory size $W\in\mathbb{N}$, step size $\\{t_{k}\\}$,
evaluation schedule $\\{I_{k}\\}$.

  2. 2.

For $k=0,1,\ldots$, repeat the following:

    1. (a)

Evaluate $f_{i}(x_{k})$ and $g_{i}(x_{k})\in\partial f_{i}(x_{k})$ for $i\in
I_{k}$.

    2. (b)

Let

$x_{k+1}=\arg\min{\breve{f}}_{k}(x),$ (2.3)

where

$\displaystyle{\breve{f}}_{k}(x)$
$\displaystyle=\sum_{i=1}^{m}{\hat{f}}_{k,i}(x)+\frac{1}{2t_{k}}\|x-x_{k}\|^{2},$
$\displaystyle{\hat{f}}_{k,i}(x)$
$\displaystyle=\left\\{\begin{array}[]{ll}\max_{l\in
B_{k,i}}f_{i}(x_{k-l})+g_{i}(x_{k-l})^{T}(x-x_{k-l}),&\text{if
}B_{k,i}\not=\emptyset\\\ 0,&\text{otherwise}\end{array}\right.$
$\displaystyle B_{k,i}$ $\displaystyle=\\{l\in\mathbb{Z}\mid 0\leq
l\leq\min\\{k,W-1\\},i\in I_{k-l}\\}.$

Note that the above algorithm is well-defined: For each $k$, $\breve{f}_{k}$
is a proper closed $(1/2t_{k})$-strongly convex function and a minimizer of
$\breve{f}_{k}$ exists and is unique [3, Theorem 5.25].

We assume the following on some analysis below.

######  Assumption 2.3.

  1. 1.

All components are evaluated at least once in every $W$ iterations:

$i\in\bigcup_{l=0}^{W-1}I_{k+l},\qquad\forall k,i.$

Under this assumption, each model has at least one cut after $W$ iterations
(that is, $B_{k,i}\not=\emptyset$ for any $k\geq W-1$ and $i$).

###  2.2 Convergence analysis

In this section, we study the behaviour of the algorithm.

The first lemma allows us to write the next point $x_{k+1}$ in terms of the
current point $x_{k}$ and the subgradients at the previous points.

######  Lemma 2.4.

Let Assumption 2.1 hold. For $k=0,1,\ldots$, there exists
$\alpha_{k}\in\mathbb{R}^{m\times W}_{+}$ such that

$x_{k+1}=P_{X}(u_{k}),$ $u_{k}=x_{k}-t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}\,g_{i}(x_{k-l}),$ $\sum_{l\in
B_{k,i}}\alpha_{k,i,l}=\left\\{\begin{array}[]{ll}1&\qquad\text{if
}B_{k,i}\not=\emptyset,\\\ 0&\qquad\text{otherwise},\\\
\end{array}\right.\quad\forall i,$ $\alpha_{k,i,l}\geq 0,\quad\forall i,l,$

where $P_{X}$ denotes the (Euclidean) projection onto the set $X$ and we adopt
the convention $\sum_{\emptyset}=0$.

We observe that the iteration resembles that of the projected (non-
incremental) subgradient methods. In our method, for each component, we
compute a convex combination of the subgradients evaluated in the last $W$
iterations (the weight is determined via an optimization problem in the proof
below), sum over all of the components, multiply the “step size” $t_{k}$ and
make a “subgradient step” followed by the projection.

###### Proof.

Let $J_{k}$ be the index set of the components with at least one cut

$J_{k}=\\{i\in\mathbb{Z}^{+}\mid B_{k,i}\not=\emptyset\\}.$

The minimization problem of ${\breve{f}}_{k}(x)$ can be written as

$\displaystyle\min\ $ $\displaystyle\sum_{i\in
J_{k}}r_{i}+\frac{1}{2t_{k}}\|x-x_{k}\|^{2}$ s.t. $\displaystyle r_{i}\geq
f_{i}(x_{k-l})+g_{i}(x_{k-l})^{T}(x-x_{k-l}),\quad\forall i\in J_{k},\ l\in
B_{k,i}$ $\displaystyle x\in X,r\in\mathbb{R}^{\lvert J_{k}\rvert}.$

This optimization problem satisfies the Slater condition [4, Proposition
5.3.1] and thus the (unique) optimal solution $(x^{*},r^{*})$ satisfies

$(x^{*},r^{*})=\arg\min_{x\in X,r\in\mathbb{R}^{\lvert
J_{k}\rvert}}L(x,r,\alpha)$

where

$\displaystyle L(x,r,\alpha)=\sum_{i\in
J_{k}}r_{i}+\frac{1}{2t_{k}}\|x-x_{k}\|^{2}-\sum_{i\in J_{k}}\sum_{l\in
B_{k,i}}\alpha_{i,l}\left(r_{i}-f_{i}(x_{k-l})-g_{i}(x_{k-l})^{T}(x-x_{k-l})\right)$

is the Lagrangian and $\alpha\geq 0$ is an optimal dual solution. The
optimality condition of $L$ gives

$\displaystyle 0\in-\partial_{x}L(x^{*},r^{*},\alpha)+N_{X}(x^{*}),\quad
0\in-\partial_{r}L(x^{*},r^{*},\alpha),$

where $N_{X}(x)$ is the normal cone to $X$ at $x$. From the second relation,
we obtain

$\displaystyle\sum_{l\in B_{k,i}}\alpha_{i,l}=1\qquad\forall i\in J_{k}$

while the first relation reads

$\displaystyle 0\in-\frac{1}{t_{k}}(x^{*}-x_{k})-\sum_{i\in J_{k}}\sum_{l\in
B_{k,i}}\alpha_{i,l}g_{i}(x_{k-l})+N_{X}(x^{*}).$

This implies

$\displaystyle x^{*}=P_{X}\left(x_{k}-t_{k}\sum_{i\in J_{k}}\sum_{l\in
B_{k,i}}\alpha_{i,l}g_{i}(x_{k-l})\right).$

By relabelling $\alpha$, we obtain the desired relations. ∎

With the preceding lemma, one can estimate the distance between iterates.

######  Lemma 2.5.

Let Assumption 2.1 hold. For any $k$ and $s$ with $k>s$, we have

$\|x_{k}-x_{s}\|\leq\overline{t}_{s,k-1}(k-s)mC,$

where $\overline{t}_{s,k-1}=\max\\{t_{s},t_{s+1},\ldots,t_{k-1}\\}$.

###### Proof.

Using Lemma 2.4 and the nonexpansion property of the projection, for all
$k\geq 0$, it follows

$\displaystyle\|x_{k+1}-x_{k}\|$ $\displaystyle=\|P_{X}(u_{k})-x_{k}\|$
$\displaystyle\leq\|u_{k}-x_{k}\|$
$\displaystyle=\left\|t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})\right\|$ $\displaystyle\leq
t_{k}\sum_{i=1}^{m}\sum_{l\in B_{k,i}}\alpha_{k,i,l}\|g_{i}(x_{k-l})\|$
$\displaystyle\leq t_{k}\sum_{i=1}^{m}\sum_{l\in B_{k,i}}\alpha_{k,i,l}C$
$\displaystyle\leq t_{k}mC,$

where we have used $\alpha_{k,i,l}\leq 1$ in the last inequality. Combining
the preceding inequality with the triangle inequality, we obtain the desired
inequality

$\displaystyle\|x_{k}-x_{s}\|$
$\displaystyle=\left\|\sum_{l=0}^{k-s-1}(x_{s+l+1}-x_{s+l})\right\|$
$\displaystyle\leq\sum_{l=0}^{k-s-1}\|x_{s+l+1}-x_{s+l}\|$
$\displaystyle\leq\sum_{l=0}^{k-s-1}t_{s+l}mC$
$\displaystyle\leq\overline{t}_{s,k-1}(k-s)mC.$

∎

Now we prove a fundamental inequality which plays a key role in the
convergence analysis of the method.

######  Lemma 2.6.

Let Assumption 2.1 hold. For any $y\in X$ and any $k\geq 0$ such that
$B_{k,i}\not=\emptyset$ for $i=1,2,\ldots,m$, we have

$\displaystyle\|x_{k+1}-y\|^{2}\leq\|x_{k}-y\|^{2}-2t_{k}(f(x_{k})-f(y))+4t_{k}\overline{t}_{k-W,k-1}m^{2}C^{2}W+t_{k}^{2}m^{2}C^{2},$

where $\overline{t}_{k-W,k-1}=\max\\{t_{k-W},t_{k-W+1},\ldots,t_{k-1}\\}$ and
we define $t_{i}=t_{0}$ for $i<0$ for notational convenience.

The condition $B_{k,i}\not=\emptyset$ for all $i$ is satisfied for any $k\geq
W-1$ if each component function is evaluated at least once in every $W$
iterations (that is, if Assumption 2.3 is satisfied).

###### Proof.

Using Lemma 2.4 and the nonexpansion property of the projection, for all
$k\geq 0$, it follows

$\displaystyle\|x_{k+1}-y\|^{2}-\|x_{k}-y\|^{2}$
$\displaystyle\quad=\|P_{X}(u_{k})-y\|^{2}-\|x_{k}-y\|^{2}$
$\displaystyle\quad\leq\|u_{k}-y\|^{2}-\|x_{k}-y\|^{2}$
$\displaystyle\quad=\left\|x_{k}-t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})-y\right\|^{2}-\|x_{k}-y\|^{2}$
$\displaystyle\quad=-2t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})^{T}(x_{k}-y)+t_{k}^{2}\left\|\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})\right\|^{2}.$

We shall bound the two terms in the last line. Using the property of
subgradient, for all $k\geq 0$, $i$ and $l\in B_{k,i}$, we obtain

$\displaystyle-g_{i}(x_{k-l})^{T}(x_{k}-y)$
$\displaystyle=-g_{i}(x_{k-l})^{T}(x_{k-l}-y)-g_{i}(x_{k-l})^{T}(x_{k}-x_{k-l})$
$\displaystyle\leq-(f_{i}(x_{k-l})-f_{i}(y))+\|g_{i}(x_{k-l})\|\cdot\|x_{k}-x_{k-l}\|$
$\displaystyle\leq-(f_{i}(x_{k})-f_{i}(y))+(f_{i}(x_{k})-f_{i}(x_{k-l}))+C\|x_{k}-x_{k-l}\|.$

The convexity of $f_{i}$ yields

$f_{i}(x_{k})-f_{i}(x_{k-l})\leq
g_{i}(x_{k})^{T}(x_{k}-x_{k-l})\leq\|g_{i}(x_{k})\|\|x_{k}-x_{k-l}\|\leq
C\|x_{k}-x_{k-l}\|$

where $g_{i}(x_{k})$ is a subgradient of $f_{i}$ at $x_{k}$. Therefore,

$\displaystyle-g_{i}(x_{k-l})^{T}(x_{k}-y)$
$\displaystyle\leq-(f_{i}(x_{k})-f_{i}(y))+C\|x_{k-l}-x_{k}\|+C\|x_{k}-x_{k-l}\|$
$\displaystyle\leq-(f_{i}(x_{k})-f_{i}(y))+2\overline{t}_{k-W,k-1}mC^{2}W,$

where the last inequality follows from Lemma 2.5. Thus, for any $k\geq 0$ such
that $B_{k,i}\not=\emptyset$ for $i=1,2,\ldots,m$,

$\displaystyle-2t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})^{T}(x_{k}-y)$ $\displaystyle\quad\leq
2t_{k}\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}\left(-(f_{i}(x_{k})-f_{i}(y))+2\overline{t}_{k-W,k-1}mC^{2}W\right)$
$\displaystyle\quad=-2t_{k}(f(x_{k})-f(y))+4t_{k}\overline{t}_{k-W,k-1}m^{2}C^{2}W,$

where the last equality holds since $\sum_{l\in B_{k,i}}\alpha_{k,i,l}=1$ for
all $i$ and $f(x)=\sum_{i=1}^{m}f_{i}(x)$ for any $x\in X$.

Now we shall bound the other term. With the triangle inequality, for all
$k\geq 0$,

$\displaystyle\left\|\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}g_{i}(x_{k-l})\right\|^{2}$
$\displaystyle\leq\left(\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}\|g_{i}(x_{k-l})\|\right)^{2}$
$\displaystyle\leq\left(\sum_{i=1}^{m}\sum_{l\in
B_{k,i}}\alpha_{k,i,l}C\right)^{2}$ $\displaystyle\leq m^{2}C^{2}.$

Combining the two bounds, we obtain the desired relation. ∎

The preceding lemma is an analogue of Lemma 2.1 from Nedic [14], which is a
basis of the convergence analysis of incremental subgradient methods.
Following their analysis of incremental subgradient methods, now we can obtain
various convergence properties. The three propositions below correspond to
Proposition 2.1, 2.4 and 2.6 from Nedic [14], respectively. Proofs are similar
to those in the reference and omitted.

######  Proposition 2.7.

Let Assumption 2.1 and 2.3 hold. Assume that step size $\\{t_{k}\\}$ is fixed
to a positive constant $t$.

  1. 1.

If $f^{*}=-\infty$, then

$\liminf_{k\rightarrow\infty}f(x_{k})=-\infty.$

  2. 2.

If $f^{*}>-\infty$, then

$\liminf_{k\rightarrow\infty}f(x_{k})\leq
f^{*}+2m^{2}tC^{2}W+\frac{1}{2}m^{2}tC^{2}.$

The above lemma estimates the worst-case suboptimality the method may suffer
when the step size is kept constant. It is counter-intuitive that the worst-
case bound gets larger as the memory size $W$ increases. We speculate that
tightening the analysis in Lemma 2.6 would yield a better bound, but it is
left as future work.

With a diminishing step size rule, the method is capable of solving the
problem asymptotically.

######  Proposition 2.8.

Let Assumption 2.1 and 2.3 hold. Assume that step size $\\{t_{k}\\}$ is such
that

$t_{k}>0\ (\forall k),\ t_{k}\geq t_{k+1}\ (\forall k),\
\lim_{k\rightarrow\infty}t_{k}=0,\ \sum_{k=0}^{\infty}t_{k}=\infty.$

Then, we have

$\liminf_{k\rightarrow\infty}f(x_{k})=f^{*}.$

######  Proposition 2.9.

Let Assumption 2.1 and 2.3 hold. Assume that step size $\\{t_{k}\\}$ is such
that

$t_{k}>0\ (\forall k),\ t_{k}\geq t_{k+1}\ (\forall k),\
\sum_{k=0}^{\infty}t_{k}=\infty,\ \sum_{k=0}^{\infty}t_{k}^{2}<\infty,$

and assume that the set of optimal solutions $X^{*}$ is nonempty. Then, the
sequence $\\{x_{k}\\}$ converges to some optimal solution.

##  3 Numerical experiments

In this section, the incremental cutting-plane method is applied to the dual
problems of unit commitment (UC) problems and the performance is compared with
that of the full-step method. The UC problem is an important optimization
problem in the energy industry to compute the optimal operating schedules of
power plants given demand over a fixed period. The formulation of the problem
is given in the appendix. One of the standard approaches is to apply the
Lagrangian relaxation to decompose the problem by generators and solve the
dual problem with the cutting-plane method [2].

We consider UC problems with 200 generators whose data is based on [8]. The
length of the planning horizon is 48 hours with a time resolution of 1 hour.
The demand data is based on the historical demand data in the UK published by
National Grid ESO.111https://www.nationalgrideso.com/ We sample 8 different
demand data and create test instances based on these demand data.

The incremental cutting-plane method used in this section is the one described
in Algorithm 2.2. For simplicity, we consider the constant step size rule
($t_{k}=t$ for all $k$) and we do not delete any cuts ($W=\infty$). When the
method is full-step, which evaluates all the component values in each
iteration, the objective value at each point is readily available. In
practice, this can be used to monitor the progress. In our numerical
experiments, we only update the iterate $x_{k}$ when the objective value gets
improved. Namely, if the objective value of the point computed by the model is
worse than that of the current point, only the model gets updated and the
iterate (i.e. the regularization centre) is kept to the same point. This
modification improves the performance of the full-step method and is often
used in practice (e.g. [15]).

The evaluation schedule $\\{I_{k}\\}$ is computed based on permutation. In our
experiments, the number of evaluated components per iteration, denoted by $p$,
is fixed to be the same in all iterations. Namely, $\lvert I_{k}\rvert=p$ for
all $k$. We consider a sequence of component functions where every $m$
elements are a permutation of $f_{1},f_{2},\ldots,f_{m}$

$\underbrace{f_{h_{1}^{1}},f_{h_{2}^{1}},\ldots,f_{h_{m}^{1}}}_{\text{1st
permutation}},\underbrace{f_{h_{1}^{2}},f_{h_{2}^{2}},\ldots,f_{h_{m}^{2}}}_{\text{2nd
permutation}},\ldots.$

In every iteration, the first $p$ components are removed from the sequence and
evaluated. This loosely follows the idea considered in the study on
incremental subgradient methods [14], which corresponds to $p=1$ where a
single component is evaluated per iteration. Setting $p=m$ we obtain the full-
step method which evaluates all of the components in each step.

For later use, we introduce an adjusted step size $\widetilde{t}=mt/p$ where
$t$ is the actual step size. We note that the method with a small value of $p$
makes more steps (i.e. solves model (2.3) more often) than the method with a
large value of $p$. If all of the component functions are affine, the method
makes the same progress per component evaluation no matter which $p$ value is
used, as long as the adjusted step size $\widetilde{t}$ is the same. Thus, the
adjusted step size $\widetilde{t}$ facilitates the comparison of experiments
with different values of $p$.

Two approaches to initialize the cutting-plane method are considered:
warmstart and coldstart. Coldstart is used when there is no prior knowledge of
the problem to be solved. However, if available, warmstart may significantly
reduce the computational time. It is of interest to verify if the incremental
update helps in this practical setup or not. In our numerical experiments, to
warmstart the method we first solve the linear programming relaxation (LPR) of
the original problem and use the optimal dual value to the LPR as the initial
value of the cutting-plane method. This initialization method was used in [15]
to warmstart the cutting-plane method. In the coldstart case, the initial dual
value is set to 0.

To observe the effect of the incremental update, the cutting-plane method is
applied to the test instances with various numbers of component evaluations
per iteration $\lvert I_{k}\rvert=p$. Table 1 shows the average computational
time, the average number of component evaluations and the average number of
iterations to find a solution within 0.1% and 0.05% tolerance. All times in
this section are wall-clock times. If the method is warmstarted, the
computational time includes the time spent to solve the LPR, which is
approximately 5 seconds. For each value of $p$, the adjusted step size
$\widetilde{t}$ is set to $0.01$ (i.e. the step size is $t=0.01p/m$). The
incremental method does not compute the function value $f(x_{k})$ for
$k=0,1,\ldots$ and thus they are evaluated after the experiments.

Table 1: Performance of the cutting-plane method with various evaluation schedule sizes $p$ |  | tol: 0.1% | 0.05%
---|---|---|---
initializer | $p/m$ | time | comp. eval. | iter. | time | comp. eval. | iter.
warmstart | 0.05 | 14.4 | 345.0 | 34.5 | 19.8 | 531.2 | 53.1
| 0.10 | 13.8 | 347.5 | 17.4 | 18.9 | 537.5 | 26.9
| 0.20 | 14.4 | 375.0 | 9.4 | 20.1 | 585.0 | 14.6
| 0.30 | 15.5 | 420.0 | 7.0 | 20.0 | 607.5 | 10.1
| 0.40 | 16.7 | 480.0 | 6.0 | 23.3 | 740.0 | 9.2
| 0.50 | 19.2 | 587.5 | 5.9 | 25.3 | 825.0 | 8.2
| 1.00 | 22.2 | 700.0 | 3.5 | 30.7 | 1025.0 | 5.1
coldstart | 0.05 | 23.5 | 911.2 | 91.1 | 39.0 | 1226.2 | 122.6
| 0.10 | 20.0 | 910.0 | 45.5 | 32.7 | 1232.5 | 61.6
| 0.20 | 19.4 | 950.0 | 23.8 | 29.2 | 1235.0 | 30.9
| 0.30 | 19.9 | 1005.0 | 16.8 | 29.5 | 1305.0 | 21.8
| 0.40 | 22.3 | 1110.0 | 13.9 | 31.8 | 1400.0 | 17.5
| 0.50 | 25.6 | 1237.5 | 12.4 | 36.1 | 1537.5 | 15.4
| 1.00 | 51.2 | 2075.0 | 10.4 | 66.9 | 2450.0 | 12.2

As shown in Table 1, for any value of $p$, the warmstarted method solves the
problem quicker than the coldstarted method. Namely, the LPR gives a dual
value close to the optimal one, which helps the method to find a solution in a
shorter time.

Interestingly, whichever initialization method is used, the incremental method
($p/m<1$) performs better than the full-step method ($p/m=1$). In general, an
incremental first-order method is expected to be more efficient than the full-
step counterpart when the iterate is far from the solution, but the
convergence to the solution of high precision is slower [5]. The incremental
cutting-plane method, however, successfully finds a solution of high precision
(e.g. 0.05%) in a short time. We note that as the number of component
evaluations per iteration $p$ gets smaller, the method tends to find a
solution with fewer component evaluations in total. However, at the same time,
the number of iterations and hence the number of solutions of model (2.3) gets
larger. If $p$ is too small, the time spent to solve (2.3) becomes significant
and the overall computational time increases. The best value of $p$ balances
the number of component evaluations in total and the growth of model solution
time.

Table 2 shows the performances of the method with various step sizes. As shown
in the table, the performance of the method is subject to the choice of the
step size. When the step size is very small, the numbers of component
evaluations required by the incremental method and the full-step method tend
to be close. With a small step size, the regularization in the model (2.3)
becomes strong and the methods cannot make much progress per iteration. The
behaviours of the two methods are similar in this situation. However, in the
setup in Table 2, the incremental method has to solve the model ten times more
often than the full-step method per component evaluation. As a result, it
requires longer computational times. On the other hand, an excessively large
step size slows down the method as well. Although this trend holds both for
the incremental method and for the full-step method, the growth of the
computational time is more significant in the full-step method. In other
words, the incremental method is more robust and tends to work relatively well
even if the step size is set larger than the best value.

Table 2: Performance of the cutting-plane method with various step sizes |  |  | tol: 0.1% | 0.05%
---|---|---|---|---
initializer | $p/m$ | $\widetilde{t}$ | time | comp. eval. | iter. | time | comp. eval. | iter.
warmstart | 0.1 | 0.100 | 26.1 | 677.5 | 33.9 | 34.0 | 907.5 | 45.4
|  | 0.050 | 20.4 | 530.0 | 26.5 | 28.1 | 772.5 | 38.6
|  | 0.010 | 13.8 | 347.5 | 17.4 | 18.9 | 537.5 | 26.9
|  | 0.005 | 13.8 | 360.0 | 18.0 | 22.4 | 690.0 | 34.5
|  | 0.001 | 27.9 | 940.0 | 47.0 | 77.1 | 2735.0 | 136.8
| 1.0 | 0.100 | 65.0 | 1775.0 | 8.9 | 101.1 | 2650.0 | 13.2
|  | 0.050 | 46.2 | 1350.0 | 6.8 | 66.0 | 1950.0 | 9.8
|  | 0.010 | 22.2 | 700.0 | 3.5 | 30.7 | 1025.0 | 5.1
|  | 0.005 | 18.7 | 600.0 | 3.0 | 27.7 | 975.0 | 4.9
|  | 0.001 | 26.6 | 950.0 | 4.8 | 70.5 | 2775.0 | 13.9
coldstart | 0.1 | 0.100 | 28.9 | 975.0 | 48.8 | 38.9 | 1195.0 | 59.8
|  | 0.050 | 24.4 | 900.0 | 45.0 | 32.7 | 1095.0 | 54.8
|  | 0.010 | 20.0 | 910.0 | 45.5 | 32.7 | 1232.5 | 61.6
|  | 0.005 | 25.2 | 1145.0 | 57.2 | 56.9 | 1902.5 | 95.1
|  | 0.001 | 143.3 | 4507.5 | 225.4 | 416.3 | 8847.5 | 442.4
| 1.0 | 0.100 | 267.9 | 4300.0 | 21.5 | 347.9 | 4850.0 | 24.2
|  | 0.050 | 186.9 | 3750.0 | 18.8 | 218.8 | 4050.0 | 20.2
|  | 0.010 | 51.2 | 2075.0 | 10.4 | 66.9 | 2450.0 | 12.2
|  | 0.005 | 35.3 | 1725.0 | 8.6 | 56.5 | 2350.0 | 11.8
|  | 0.001 | 98.5 | 4500.0 | 22.5 | 235.2 | 8875.0 | 44.4

Figure 1 shows the objective values at the points computed by the method on
one of the test instances. The left and right figures correspond to the
warmstart and the coldstart respectively. In each plot, the full-step method
and incremental method are shown in solid and dashed lines, respectively. The
warmstarted method first solves the LPR, which approximately take 5 seconds,
while the coldstarted method starts outputting the iterates immediately. In
both cases, the full-step method with step size 0.1 often outputs points whose
objective values are out of the range and do not appear in the figures.

$0$$5$$10$$15$$20$0.9900.9940.9981.002elapse (second)obj. valuewarmstart

$0$$5$$10$$15$$20$$0.8$$0.85$$0.9$$0.95$$1$elapse (second)obj.
valuecoldstartstepsize

Figure 1: Objective values of the iterates by the warmstarted (left) and
coldstarted (right) method. Solid and dashed lines correspond to the full-step
method and the incremental method respectively.

The figure on the left shows that the strength of the regularization (i.e. the
reciprocal of the step size) has a significant effect on the behaviour of the
full-step method, which is plotted in the solid lines. We note that all of the
cases use the same initial point which is close to the optimal solution,
computed by the LPR. However, the full-step method with the weakest
regularization (step size 0.1) is quite unstable and it outputs dual values
far away from the optimal one. This is expected behaviour since the model does
not have enough cuts in the first few iterations and without suitable
regularization the model tends to output points far from the solution, even if
it is initialized at a near-optimal point. When the step size is 0.01, the
full-step method still shows instability at first but makes constant progress
afterwards. If the step size decreased to the further smaller value 0.001, the
full-step method does not show such unstable behaviour at all but the entire
progress becomes significantly slow. As a result, the full-step method with
step size 0.001 takes longer than that of step size 0.01.

The same discussion applies to the warmstarted incremental method. However,
the instability of the incremental method observed with step size 0.1 and 0.01
is drastically improved compared with the full-step method. The incremental
method also shows the unstable behaviour at first but becomes stable more
quickly than the full-step method, which leads to the reduction of the overall
computational time. However, the incremental update does not fix the slow
progress caused by the too small step size 0.001 and the full-step method and
incremental method behave quite similarly.

Similar trends are observed with the coldstart case which is plotted on the
right. We still observe the advantage of the incremental update with the large
step sizes (0.1 and 0.01). Namely the incremental update mitigate the
instability of the cutting-plane method. However, the incremental method with
the small step size (0.001) is rather slower than the full-step method.

##  4 Conclusion

In this paper, we considered an incremental variant of regularized cutting-
plane methods. First, the convergence property of the method was studied under
the limited-memory assumption. That is, we considered a setting where cuts
were deleted after a finite number of iterations. By showing the similarity of
the incremental cutting-plane method and incremental subgradient methods, the
convergence of the method was established. Then, the method was applied to the
dual problems of large-scale unit commitment problems. In many settings, the
incremental method found a solution of high precision in a shorter time than
the full-step method. The incremental method mitigated the instability of the
cutting-plane method observed in the first few iterations, which seemed to be
one of the factors to account for the superior performance of the incremental
method.

## References

  * [1] W. Ackooij and A. Frangioni  “Incremental Bundle Methods using Upper Models”  In _SIAM Journal on Optimization_ 28.1, 2018, pp. 379–410  DOI: 10.1137/16M1089897
  * [2] W. Ackooij et al.  “Large-scale unit commitment under uncertainty: an updated literature survey”  In _Annals of Operations Research_ 271.1 New York: Springer US, 2018, pp. 11–85
  * [3] A. Beck  “First-Order Methods in Optimization”  Philadelphia, PA: Society for IndustrialApplied Mathematics, 2017  DOI: 10.1137/1.9781611974997
  * [4] D.. Bertsekas  “Convex optimization theorem”  Belmont, Massachusetts: Athena Scientific, 2009
  * [5] D.. Bertsekas  “Incremental Gradient, Subgradient, and Proximal Methods for Convex Optimization: A Survey”  In _Optimization for Machine Learning_ Cambridge, Massachusetts: The MIT Press, 2011, pp. 85–119
  * [6] D.. Bertsekas and J.. Tsitsiklis  “Neuro-dynamic programming”  In _Neuro-dynamic programming_ , Optimization and neural computation series ; 3  Belmont, Massachusetts: Athena Scientific, 1996
  * [7] J.. Bonnans, J.. Gilbert, C. Lemaréchal and C.. Sagastizábal  “Numerical optimization theoretical and practical aspects”  In _Numerical optimization theoretical and practical aspects_ , Universitext  Berlin: Springer, 2006
  * [8] A. Borghetti, A. Frangioni, F. Lacalandra and C.. Nucci  “Lagrangian Heuristics Based on Disaggregated Bundle Methods for Hydrothermal Unit Commitment”  In _IEEE Power Engineering Review_ 22.12 IEEE, 2002, pp. 60–60
  * [9] O. Briant et al.  “Comparison of bundle and classical column generation”  In _Mathematical programming_ 113.2 Berlin, Heidelberg: Springer-Verlag, 2008, pp. 299–344
  * [10] E.. Cheney and A.. Goldstein  “Newton’s method for convex programming and Tchebycheff approximation”  In _Numerische Mathematik_ 1.1, 1959, pp. 253–268  DOI: https://doi.org/10.1007/BF01386389
  * [11] G. Emiel and C. Sagastizábal  “Incremental-like bundle methods with application to energy planning”  In _Computational optimization and applications_ 46.2 Boston, Massachusetts: Springer US, 2010, pp. 305–332  DOI: 10.1007/s10589-009-9288-8
  * [12] J. Hiriart-Urruty and C. Lemaréchal  In _Convex analysis and minimization algorithms II_ , Grundlehren der mathematischen Wissenschaften 306  Berlin: Springer, 1993  DOI: 10.1007/978-3-662-06409-2
  * [13] J.. Kelley  “The Cutting-Plane Method for Solving Convex Programs”  In _Journal of the Society for Industrial & Applied Mathematics_ 8.4 Society for IndustrialApplied Mathematics, 1960, pp. 703–712
  * [14] A. Nedic and D.. Bertsekas  “Incremental Subgradient Methods for Nondifferentiable Optimization”  In _SIAM Journal on Optimization_ 12.1, 2001, pp. 109–138  DOI: 10.1137/S1052623499362111
  * [15] T. Schulze, A. Grothey and K. McKinnon  “A Stabilised Scenario Decomposition Algorithm Applied to Stochastic Unit Commitment Problems”  In _European Journal of Operational Research_ 261.1 Elsevier B.V., 2017, pp. 247–259  DOI: 10.1016/j.ejor.2017.02.005
  * [16] S. Shalev-Shwartz  “Online Learning and Online Convex Optimization”  In _Foundations and Trends in Machine Learning_ 4.2, 2012, pp. 107–194  DOI: 10.1561/2200000018

##  Appendix A Problem formulation

We follow one of the standard formulations in literature and formulate the
following constraints:

  * •

Load balance: Generators have to meet the all demand in each time period.

  * •

Reserve: To deal with contingencies, it is required to keep a sufficient
amount of back up in each time period, which can be activated quickly.

  * •

Power output bounds: Each generator’s power output has to be within its limit.

  * •

Ramp rate: Generators can only change their outputs within the ramp rates.

  * •

Minimum up/downtime: If switched on (off), each generator has to stay on (off)
for a given minimum period. This is to avoid thermal stress in the generators
which may cause wear and tear of the turbines.

To formulate the model, we use the following notation.

  * •

Sets

    * –

$G=\\{1,2,\ldots,n_{G}\\}$: set of generators

    * –

$T=\\{1,2,\ldots,n_{T}\\}$: set of time indices where decisions are taken

  * •

Parameters

    * –

$C^{\mathrm{nl}}_{g}$: no-load cost of generator $g$

    * –

$C^{\mathrm{mr}}_{g}$: marginal cost of generator $g$

    * –

$C^{\mathrm{up}}_{g}$: startup cost of generator $g$

    * –

$P^{\max/\min}_{g}$: maximum/minimum generation limit of generator $g$

    * –

$P^{\mathrm{ru}/\mathrm{rd}}_{g}$: operating ramp up/down limits of generator
$g$

    * –

$P^{\mathrm{su}/\mathrm{sd}}_{g}$: startup/shutdown ramp limits of generator
$g$

    * –

$T^{\mathrm{u}/\mathrm{d}}_{g}$: minimum uptime/downtime of generator $g$

    * –

$P^{\mathrm{d}}_{t}$: power demand at time $t$

    * –

$P^{\mathrm{r}}_{t}$: reserve requirement at time $t$

  * •

Variables

    * –

$\alpha_{gt}\in\\{0,1\\}$: 1 if generator $g$ is on in period $t$, and 0
otherwise

    * –

$\gamma_{gt}\in\\{0,1\\}$: 1 if generator $g$ starts up in period $t$, and 0
otherwise

    * –

$\eta_{gt}\in\\{0,1\\}$: 1 if generator $g$ shut down in period $t$, and 0
otherwise

    * –

$p_{gt}\geq 0$: power output of generator $g$ in period $t$

The objective is the total cost

$\min\sum_{t\in T}\sum_{g\in
G}\left(C^{\mathrm{nl}}_{g}\alpha_{gt}+C^{\mathrm{mr}}_{g}p_{gt}+C^{\mathrm{up}}_{g}\gamma_{gt}\right).$

This is to be minimized subject to the following constraints.

  * •

Load balance

$\sum_{g\in G}p_{gt}\geq P^{\mathrm{d}}_{t}\qquad t\in T$

  * •

Reserve

$\sum_{g\in G}(P^{\max}_{g}\alpha_{gt}-p_{gt})\geq P^{\mathrm{r}}_{t}\qquad
t\in T$

  * •

Power output bounds

$P^{\min}_{g}\alpha_{gt}\leq p_{gt}\leq P^{\max}_{g}\alpha_{gt}\qquad g\in
G,t\in T$

  * •

Ramp rate

$p_{gt}-p_{g\,t-1}\leq
P^{\mathrm{ru}}_{g}\alpha_{g\,t-1}+P^{\mathrm{su}}_{g}\gamma_{gt}\qquad g\in
G,t\in T\backslash\\{1\\}$ $p_{g\,t-1}-p_{gt}\leq
P^{\mathrm{rd}}_{g}\alpha_{gt}+P^{\mathrm{sd}}_{g}\eta_{gt}\qquad g\in G,t\in
T\backslash\\{1\\}$

  * •

Minimum up/downtime

$\sum_{i=\max\\{t-T^{\mathrm{u}}_{g}+1,1\\}}^{t}\gamma_{gi}\leq\alpha_{gt}\qquad
g\in G,t\in T$ $\sum_{i=\max\\{t-T^{\mathrm{u}}_{g}+1,1\\}}^{t}\eta_{gi}\leq
1-\alpha_{gt}\qquad g\in G,t\in T$

  * •

Polyhedral/Switching constraints (to enforce binaries to work as we expect)

$\alpha_{gt}-\alpha_{g\,t-1}=\gamma_{gt}-\eta_{gt}\qquad g\in G,t\in T$
$1\geq\gamma_{gt}+\eta_{gt}\qquad g\in G,t\in T$