# DBNet: A Dual-branch Network Architecture Processing on Spectrum and
Waveform for Single-channel Speech Enhancement

###### Abstract

In real acoustic environment, speech enhancement is an arduous task to improve
the quality and intelligibility of speech interfered by background noise and
reverberation. Over the past years, deep learning has shown great potential on
speech enhancement. In this paper, we propose a novel real-time framework
called DBNet which is a dual-branch structure with alternate interconnection.
Each branch incorporates an encoder-decoder architecture with skip
connections. The two branches are responsible for spectrum and waveform
modeling, respectively. A bridge layer is adopted to exchange information
between the two branches. Systematic evaluation and comparison show that the
proposed system substantially outperforms related algorithms under very
challenging environments. And in INTERSPEECH 2021 Deep Noise Suppression (DNS)
challenge, the proposed system ranks the top 8 in real-time track 1 in terms
of the Mean Opinion Score (MOS) of the ITU-T P.835 framework.

Index Terms: Speech enhancement, time domain, frequency domain, real-time

##  1 Introduction

With the COVID-19 pandemic, lots people are working online, and the demand for
reliable real-time speech enhancement algorithms has increased sharply. During
these times, we need to ensure clear call quality and effective communication
and cooperation with others without delay. Our communication are usually
disturbed by a lot of background noise, such as washing machines, passing
trucks, and the strong reverberation. These will affect the efficiency of our
work and communication. Recently, many researchers from academia and industry
have made significant contributions to monaural speech enhancement. However,
due to the diversity of noise types in reality, even the the state-of-the-art
algorithms cannot handle challenging environments well.

With the development of deep learning, many research regard speech enhancement
as a supervised learning problem [1] [2] [3], and obtained excellent
performance. Usually the neural network input is time domain signal or STFT
domain signal. In [2] [4], the authors study the noise reduction problem in
the short-time Fourier transform (STFT) domain. In [1] [5], the time-domain
signal after framing is directly feed to the neural network. Both frequency
domain and time domain has its own advantages. The STFT method is more in line
with human hearing perception, and the characteristics of speech is more
explicit. Time domain method does not damage the signal by STFT, avoids the
well-known invalid STFT problem.

Lim et al. [6] introduced a time-frequency network to jointly optimize the
time and frequency domains of a signal for the task of audio super resolution.
They show that combine these two domains can boost the audio super resolution
performance and obtain the state-of-the-art in both quantitatively and
qualitatively.

(a) Time domain noisy with an impulse-like noise

(b) Frequency domain noisy with an impulse-like noise

(c) Time domain noisy with a narrowband-like noise

(d) Frequency domain noisy with a narrowband-like noise

Figure 1: Time- and frequency-domain noisy signal with two kinds of noises.

As we know, for an impulse-like noise as shown in Fig. 1(a) and 1(b), it is
easy to eliminate in the time domain. Only a few samples need to be removed.
In the frequency domain, the impulse-like nose pollutes the entire frequency
band, and it is difficult to be eliminated by a speech enhancement method
based on the frequency domain. In contrast, for a narrowband-like noise as
shown in Figure 1(c) and 1(d), the noise is distributed on the narrowband and
the frequency domain-based method covers well. In the time-domain, the noise
and the speech are coupled together, and it is hard to decouple by an
enhancement method based on time-domain. In this paper, in order to reduce
noise better, we proposed a novel speech enhancement algorithm called DBNet,
which combines the time domain and frequency domain together. The DBNet is a
dual-branch structure with alternate interconnection. Each branch incorporates
an encoder decoder architecture with skip connections. The two branches are
responsible for spectrum and waveform modeling, respectively. A bridge layer
is adopted to exchange information between the two branches. Experiments show
that the proposed method has excellent results on the WSJ0 SI-84 [7] and DNS
Challenge [8].

The rest of the paper is organized as follows. The structure of the proposed
DBNet is described in Section 2. Section 3 describes the experimental setup.
The experimental results are revealed in Section 4. We conclude this paper in
Section 5.

Figure 2: Diagram of the proposed DBNet model

##  2 Model Description

The overall structure of the proposed method is shown in Fig. 2. In the
following sections, three important parts: SRS module, GLSTM, bridge layer,
will be introduced one by one.

###  2.1 Shift Real Spectral (SRS) Module

Recently, Liu et al. [9] proposed a new separation target based on the time-
frequency representation of SRS in their work, which proved the superiority of
SRS over STFT. First, SRS takes the phase into consideration, which improves
the speech intelligibility and quality. Second, SRS is a spectral
representation method in real field instead of the complex field, all of the
elements of input are real numbers. So it reduces the modeling difficulty and
provide convenience for the information interaction module of our model.
Therefore in this paper, SRS is adopted as our frequency domain input based on
the two advantages above.

###  2.2 Gated Convolution and Group LSTM

Dauphin et al. [10] improved the masked convolution for image convolution
modeling and proposed gated convolution (GCNN) which is described as:

$\displaystyle y$ $\displaystyle=(x\ast W_{1}+b_{1})\odot\sigma{(x\ast
W_{2}+b_{2})}$ (1) $\displaystyle=v_{1}\odot\sigma{(v_{2})}$

where $W$ and $b$ represent kernel and bias, respectively. $\ast$ and $\odot$
denote operation of convolution and element-wise multiplication, respectively.
$\sigma$ represents a nonlinear activation function. The GCNN can reduce the
vanishing gradient problem for deep architectures by providing a linear path
for the gradients, so we replaced the convolution in the original crn with it.
A diagram of gated convolution is shown in Fig. 3.

Model efficiency is pretty important, and many application scenarios have
higher requirements for processing speed and memory usage. However, due to the
introduction of the dual-branch structure, the amount of calculation and
memory occupied by the model are much higher than that of the standard
convolutional recurrent network. Gao et al. [11] proposed a grouped recursive
neural network (RNN) strategy, which reduces the complexity of the model while
ensuring the performance of the model. The process of a group RNN is shown in
Fig. 4.

In this paper, the group LSTM contains two layers of RNN, and each layer has
two LSTMs to learn the features within each group. Between the two layers, a
frame-level rearrangement is used to establish the inter-group relationship of
features, which guarantees the utilization of inter-group correlations to a
certain extent.

(a) Gated Convolution

(b) Gated Deconvolution

Figure 3: Gated Convolution and Gated Deconvolution, where $\alpha$ denotes a
sigmoid funtion Figure 4: The process of group LSTM

###  2.3 Bridge Layer

Bridge Layer is a linear unit, which is responsible for converting information
from one branch to another. Actually, the bridge layer consists of vectors
with the same length as the frame length. We take the real part of the fast
Fourier transform (FFT) variable as the initialization parameter of the
vectors to fit the situation that SRS is used as the frequency domain
representation.

###  2.4 Parameters Setting

The encoder has six layers, each consist of a bridge layer, a gated
convolution followed by batch normalization and ELU nolinearity. Note that a
module is added before the input of frequency branch to calculate the
frequency-domain representation of the real number field. The input size of
the model is [batch_size, 1, seq_len, features]. The first layer of the
encoder expands the input channel from 1 to 64. After that,the rest layer of
the encoder do the following operations. The bridge layer first translates the
features from the other branch and then concatenate it with the output of the
previous layer along the channel axis, finally pass the features to a gated
convolution. The kernel and the stride of convolution are set to (1, 3) and
(1, 2), respectively. The final output size of the encoder is [batch_size, 64,
seq_len, features / 64]. The output is fed into group LSTM layers to extract
the long-term relationship of the features in time domain and frequency
domain, respectively.

Similarly, the decoder also has six layers. In addition to the features from
another domain, each layer also gets the skip connection of the corresponding
layer of the encoder, and concatenate them with the output of the previous
layer along the channel axis. The decoder uses gated deconvolution to double
the feature dimensions layer by layer to reconstruct the signal to the
original size. The last layer of decoder enhance the signal to one channel.
Finally, the output is converted into speech through the overlap-and-add
operation. Note that the speech from frequency branch is regarded as final
result.

###  2.5 Loss Function

In the early experiments, we used a loss based on STFT magnitude which was
proposed in [1] and can be described as:

$\displaystyle L_{Mag}(s,\hat{s})=\frac{1}{T\cdot F}\sum_{t=0}^{T-1}$
$\displaystyle\sum_{f=0}^{F-1}[(|S_{r}(t,f)|+|S_{i}(t,f)|)$ (2)
$\displaystyle-(|\hat{S}_{r}(t,f)|+|\hat{S}_{i}(t,f)|)]$

where $T$ and $F$ represent the number of time frames and frequency
dimensions, $S$ and $\hat{S}$ denote STFTs of s and $\hat{s}$, respectively.
$S_{r}$ and $\hat{S}_{i}$ represent the real and imaginary parts of $S$,
respectively. Note that the output of the network contains two enhanced
utterances, one of which is from the time branch and the other from the
frequency branch and they are optimized independently. So the total loss is
defined as:

$L=L_{Mag}(s,\hat{s}_{t})+L_{Mag}(s,\hat{s}_{f})$ (3)

However, we found that magnitude loss introduced a large number of unknown
artifacts. Although it does not affect the objective evaluation score, it
would bring terrible auditory perception. Therefore, in DNS Challenge, the
magnitude loss is replaced with the phase constrained magnitude loss proposed
in [12] and achieved good subjective evaluation scores in the competition.

##  3 Experiment

###  3.1 Datasets

In this study, we evaluated the performance of our proposed model on the WSJ0
SI-84 dataset [7] which includes 7138 utterances from 83 speakers (42 males
and 41 females). We used the utterances of 77 speakers for training and the
rest for test. We used 10000 non-speech sounds from a sound effect library
(available at www.sound-ideas.com) [13] and generated 320000 and 3000
utterances at the SNRs uniformly sampled from {-5dB, -4dB, -3dB, -2dB, -1dB,
-0dB} for training and validation, respectively. For the test set, two noises
(babble and cafeteria) from Auditec CD (available at http://www.auditec.com)
are used to generate 300 mixtures at each SNR of -5dB, 0dB, and 5dB.

###  3.2 Baselines

In this study, we compared the proposed dual-branch network with another 3
baselines, namely CRN [14], GCRN [2] and AECNN [1], which are given as
follows:

  * •

CRN: it is a casual convolutional recurrent network in T-F domain. The network
uses 5 convolution layers as the encoder and 5 deconvolution layers as the
decoder. Two LSTM layers are used for sequence modeling. This network
receiving magnitude as input. The number of channels is decreased and the
number of parameters is 4.5M.

  * •

GCRN: it is a causal gated convolutional recurrent network for complex
spectral mapping. The structure is similar to CRN except that GCRN has two
decoders to model real and imaginary, respectively. The input of network is
complex spectral. We kept the best configuration in [2] and the number of
parameters is 9.76M.

  * •

AECNN: it is a autoencoder based fully convolutional neural network in the
time domain. Raw waveform is chunked into frames with a large tiem frame size
(1.024s). We kept the best configuration in [1]. The number of parameters is
18M.

  * •

DBNet: the structure of two branches are same. 6 (de) convolution block are
set for encoder and decoder. The number of channels is 64 for each layer. The
kernel size of (1, 3) and stride of (1, 2) are set for time and frequency
axis. The input is time frames and SRS for time branch and frequency branch,
respectively. The numebr of parameters is 2.9M.

###  3.3 System Settings

All utterances are sample at 16kHz. The frames are extracted using a
rectangular window and a hamming window of size 20ms for time domain and
frequency domain, respectively. The overlap is 10ms. The models are trained
using the Adam optimizer [15] with a learning rate of 0.001. And the batch
size is set to 32 at the utterance level. Note that a random segment of 7
seconds is intercepted from an utterance if it is larger than 7 seconds. The
smaller utterances are zero-padded to match the size of the largest utterance
in the batch.

###  3.4 Evaluate Metrics

The performance is evaluated with two objective metrics: short-time objective
intelligibility (STOI) [16] and perceptual evaluation of speech quality (PESQ)
[17]. STOI values typically range from 0 to 1, which can be roughly
interpreted as percent correct. PESQ values range from -0.5 to 4.5. For both
of the STOI and PESQ metrics, the higher number indicates better performance.

##  4 Result and Analysis

Table 1: STOI and PESQ comparisons between DBNet and the baseline models

Metrix | STOI | PESQ | Casual?
---|---|---|---
Test Noise | Babble | Cafeteria | Babble | Cafeteria
Test SNR | -5db | 0db | 5db | AVG | -5db | 0db | 5db | AVG | -5db | 0db | 5db | AVG | -5db | 0db | 5db | AVG
Mixture | 59.1 | 71.6 | 83.3 | 71.3 | 57.5 | 70.1 | 82.4 | 70.0 | 1.52 | 1.76 | 2.07 | 1.78 | 1.40 | 1.69 | 2.07 | 1.72
CRN | 77.9 | 88.0 | 93.2 | 86.4 | 75.7 | 86.6 | 92.7 | 85.0 | 1.99 | 2.50 | 2.91 | 2.47 | 2.01 | 2.47 | 2.89 | 2.46 | ✓
GCRN | 81.1 | 90.7 | 94.7 | 88.8 | 78.0 | 88.8 | 94.1 | 87.0 | 2.03 | 2.65 | 3.07 | 2.58 | 1.97 | 2.56 | 3.03 | 2.52 | ✓
AECNN | 80.5 | 90.6 | 94.2 | 88.4 | 80.0 | 89.4 | 94.0 | 87.8 | 2.00 | 2.57 | 2.88 | 2.48 | 2.03 | 2.55 | 2.93 | 2.50 | ✓
DBNet | 83.1 | 91.8 | 95.3 | 90.0 | 79.6 | 89.9 | 94.8 | 88.1 | 2.15 | 2.78 | 3.16 | 2.70 | 2.08 | 2.69 | 3.10 | 2.62 | ✓
GCRN-NC | 84.1 | 92.1 | 95.5 | 90.5 | 81.3 | 90.5 | 95.0 | 89.0 | 2.22 | 2.81 | 3.17 | 2.73 | 2.08 | 2.72 | 3.07 | 2.62 | ×
DBNet-NC | 87.3 | 93.5 | 96.2 | 92.3 | 84.3 | 91.9 | 95.7 | 90.6 | 2.56 | 3.06 | 3.34 | 2.99 | 2.43 | 2.95 | 3.31 | 2.90 | ×

###  4.1 Objective Comparisons

First, we compare our model with baselines in terms of both STOI and PESQ. The
results are given in Table 1 and the best results are marked in bold. For
STOI, DBNet is better than all baselines at all SNRs and noises except AECNN.
However, AECNN is a model based on large frame size so it is not suitble for
real-time scenes. Compared with GCRN, an average improvement of 1.20% and
1.10% is obtained for babble and cafeteria, respectively. For PESQ, GCRN is
the best baseline and an average improvement of 0.12 and 0.10 is obtained for
babble and cefeteria, respectively. As for non causal system, DBNet-NC
outperforms GCRN-NC and obtained an average improvement of 1.7 for STOI and
0.27 for PESQ.

In conclusion, the proposed dual-branch model outperforms both AECNN which is
a time-domain based model and GCRN which is a frequency-domain based model for
complex spectrogram mapping, indicating that alternate interconnection of the
information of two domains can significantly improve the performance of the
model and improve parameter utilization.

###  4.2 Deep Noise Suppression Challenge

The Deep Noise Suppression (DNS) challenge is designed to foster innovation in
the area of noise suppression to achieve superior perceptual speech quality.
To evaluate the performance of the proposed method on more complicated and
real acoustic scenarios, the proposed model was trained with the DNS-Challenge
wide band dataset which contains more complex acoustic scenarios including
reverberation, singing, emotions, and non-English speech. The settings of
generating training set is described as follow. The SNRs of training mixtures
vary from -5dB to 25dB. Around 30% of utterances are convolved with provided
synthetic and real room impulse responses (RIRs) before mixed with different
noise signals, and we process speech, noise and reverberation by using the
spectral augmentation filters proposed in [18]. Moreover, there is a 5% chance
that there may be multiple compound noises in one utterance. To meet the
requirement of DNS-Challenge, the number of channels is appropriately
decreased. In addition, kernel size of convolution blocks is set to (2, 3).

We use DNSMOS [19] which is a reliable non-intrusive objective speech quality
metric as our evaluation metrics at training stag and take Mean Opinion Score
(MOS) of the ITU-T P.835 framework as result. The results of the evaluation
using the ITU-T P.835 criterion [20] which is provided by the organizer are
shown in Table 2. It is obvious that the proposed model outperforms the
baseline (NSnet2) by overall 0.12 DMOS score. Then we calculated the
processing latency of our algorithm according to the competition requirements.
In this model, the frame size T = 20ms, and the overlap between consecutive
frames Ts = 10ms, so the latency T =30ms, which meets the requirements. We
also evaluated the memory access cost (MAC), and the result is 2.847G per
second.

Table 2: Subjective evaluation with P.835 criterion on the DNS Challenge

| Stationary | Emotional | Tonal | Non-English | Musical | English | Overall
---|---|---|---|---|---|---|---
Noisy | 3.03 | 2.28 | 3 | 3.04 | 2.57 | 2.52 | 2.77
NSnet2 | 3.28 | 2.75 | 3.31 | 3.25 | 2.78 | 2.93 | 3.07
DBNet | 3.44 | 3.07 | 3.37 | 3.43 | 2.78 | 2.93 | 3.19

##  5 Conclusions

In this study, we propose a novel single-channel speech enhancement system,
which consists of two denoising branches on time domain and frequency domain.
The results turns out that the proposed model outperforms other advanced
models in terms of objective intelligibility and quality scores.

We explain our work has excellent performance because the information in the
time domain and the frequency domain is not exactly the same. According to the
principle of STFT, convolution in the time domain is equivalent to direct
product in the frequency domain. Operations in the time domain tend to focus
more on local information, and operations in the frequency domain focus more
on the relationship between frames. A reasonable combination of the two can
achieve better performance. And the proposed model has fewer parameters which
indicate the two branch structure improves parameter utilization. Subjective
results showed that the proposed system ranks the top 8 in the Mean Opinion
Score (MOS) of the ITU-T P.835 for real-time track 1 of the INTERSPEECH 2021
Deep Noise Suppression (DNS) challenge.

##  6 Acknowledgements

The author would like to thank Yongjie Yan, Tailong Zhang and Pengjie Shen for
their valuable comments. This research was partly supported by the China
National Nature Science Foundation (No. 61876214).

## References

  * [1] A. Pandey and D. Wang, “A new framework for cnn based speech enhancement in the time domain,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. PP, pp. 1–1, 04 2019.
  * [2] K. Tan and D. Wang, “Learning complex spectral mapping with gated convolutional recurrent networks for monaural speech enhancement,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 380–390, 2019.
  * [3] D. Wang and J. Chen, “Supervised speech separation based on deep learning: An overview,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 10, pp. 1702–1726, 2018.
  * [4] Q. Li, F. Gao, H. Guan, and K. Ma, “Real-time monaural speech enhancement with short-time discrete cosine transform,” _arXiv preprint arXiv:2102.04629_ , 2021.
  * [5] A. Pandey and D. Wang, “Densely connected neural network with dilated convolutions for real-time speech enhancement in the time domain,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6629–6633.
  * [6] T. Y. Lim, R. A. Yeh, Y. Xu, M. N. Do, and M. Hasegawa-Johnson, “Time-frequency networks for audio super-resolution,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 646–650.
  * [7] D. B. Paul and J. M. Baker, “The design for the wall street journal-based csr corpus,” in _The Second International Conference on Spoken Language Processing, ICSLP 1992, Banff, Alberta, Canada, October 13-16, 1992_ , 1992.
  * [8] C. K. Reddy, H. Dubey, K. Koishida, A. Nair, V. Gopal, R. Cutler, S. Braun, H. Gamper, R. Aichner, and S. Srinivasan, “Interspeech 2021 deep noise suppression challenge,” _arXiv preprint arXiv:2101.01902_.
  * [9] Y. Liu, H. Zhang, and X. Zhang, “Using shifted real spectrum mask as training target for supervised speech separation,” 09 2018, pp. 1151–1155.
  * [10] Y. N. Dauphin, A. Fan, M. Auli, and D. Grangier, “Language modeling with gated convolutional networks,” 2016.
  * [11] F. Gao, L. Wu, L. Zhao, T. Qin, X. Cheng, and T.-Y. Liu, “Efficient sequence learning with group recurrent networks,” _NAACL-HLT_ , pp. 799–808, 2018\.
  * [12] A. Pandey and D. Wang, “Dense CNN with Self-Attention for Time-Domain Speech Enhancement,” _arXiv e-prints_ , p. arXiv:2009.01941, Sep. 2020.
  * [13] J. Chen, Y. Wang, S. E. Yoho, D. L. Wang, and E. W. Healy, “Large-scale training to increase speech intelligibility for hearing-impaired listeners in novel noises,” _Journal of the Acoustical Society of America_ , vol. 139, no. 5, p. 2604, 2016.
  * [14] K. Tan and D. L. Wang, “A convolutional recurrent neural network for real-time speech enhancement,” in _Proc. Interspeech 2018_ , 2018, pp. 3229–3233. [Online]. Available: http://dx.doi.org/10.21437/Interspeech.2018-1405
  * [15] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” _arXiv e-prints_ , p. arXiv:1412.6980, Dec. 2014.
  * [16] C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen, “An algorithm for intelligibility prediction of time–frequency weighted noisy speech,” _IEEE Transactions on Audio Speech & Language Processing_, vol. 19, no. 7, pp. 2125–2136, 2011.
  * [17] A. W. Rix, M. P. Hollier, A. P. Hekstra, and J. G. Beerends, “Perceptual evaluation of speech quality (pesq) the new itu standard for end-to-end speech quality assessment: Part i: Time-delay compensation,” _Journal of the Audio Engineering Society_ , vol. 50, no. 10, pp. 755–764, 2002.
  * [18] J. M. Valin, “A hybrid dsp/deep learning approach to real-time full-band speech enhancement,” 2017.
  * [19] C. K. Reddy, V. Gopal, and R. Cutler, “Dnsmos: A non-intrusive perceptual objective speech quality metric to evaluate noise suppressors,” _arXiv e-prints_ , pp. arXiv–2010, 2020.
  * [20] B. Naderi and R. Cutler, “A crowdsourcing extension of the itu-t recommendation p. 835 with validation,” _arXiv preprint arXiv:2010.13200_ , 2020.


# An updated view and perspectives on high-energy gamma-ray emission from SGR
J1935+2154 and its environment

Jaziel G. Coelho11footnotetext: Corresponding author. Luana N. Padilha  Rita
C. dos Anjos  Cynthia V. Ventura  and Geanderson A. Carvalho

###### Abstract

SGR J1935+2154 was discovered in 2016 and is currently one of the most burst-
active Soft Gamma-ray Repeaters (SGR), having emitted many X-ray bursts in
recent years. In one of our previous articles, we investigated the
contribution to high-energy and very high-energy gamma-ray emission (VHE,
$E>100$ GeV) due to cosmic-ray acceleration of SNR G57.2+0.8 hosting SGR
J1935+2154 using the GALPROP propagation code. However, follow-up observations
of SGR 1935+2154 were made for 2 hours on April 28, 2020, using the High
Energy Stereoscopic System (H.E.S.S.). The observations coincide with X-ray
bursts detected by INTEGRAL and Fermi/Gamma-ray Burst Monitor (GBM). These are
the first high-energy gamma-ray observations of an SGR in a flaring state, and
upper limits on sustained and transient emission have been derived. Now that
new H.E.S.S. observations have been made, it is interesting to update our
model with respect to these new upper limits. We extend our previous results
to a more general situation using the new version of GALPROP. We obtain a
hadronic model that confirms the results discussed by H.E.S.S.. This leads to
an optimistic prospect that cosmic ray gamma rays from SGR J1935+2154 can
contribute to the overall gamma energy density distribution and in particular
to the diffusion gamma rays from the Galactic center.

##  1 Introduction

The Soft Gamma-ray Repeater (SGR) J1935+2154 was initially detected by the BAT
(Burst Alert Telescope) instrument aboard the Swift satellite as an X-ray
burst. Subsequent observations of this source allowed to classify it as a
magnetar and they found that the source became active again in April 2020,
whilst it exhibited multiple and severe X-ray burst activity (see [1, 2, 3]
and references therein). SGRs are a diverse set of sources with huge magnetic
fields of the order of $10^{14}-10^{15}$ G, rotational periods $P\sim(2-12)$
s, slowing down rates $\dot{P}\sim(10^{-15}-10^{-10})$ s/s, persistent X-ray
luminosity as large as $10^{35}$ erg/s, transient activity in the form of
outbursts of energy around $10^{41}-10^{43}$ erg, and, for some sources, the
presence of giant flares, whose typical luminosities are $10^{44}-10^{47}$ erg
(see, e.g., [4] and references therein). The emission nature of SGRs remains a
cause for debate, and several scenarios were proposed to explain their
observed spectra and properties [see 5, 6, 7, for detailed reviews]. Examples
thereof are magnetars [8, 9, 10, 11, 12, 13], accreting neutron stars [14,
15], rotation-powered pulsars [16], quark stars [17], or massive, fast-
spinning and highly magnetized white dwarfs [18, 19, 20, 21, 22, 23].

On April 28, 2020, a very bright radio outburst of SGR J1935+2154 was
identified that turned out to be brighter than any radio burst seen from any
galactic source to date. An extremely bright millisecond-duration radio burst
was emitted by SGR J1935+2154 and detected with the Canadian Hydrogen
Intensity Mapping Experiment (CHIME) and STARE2 telescopes [24, 25]. The radio
burst turned out to be temporally coincident with a bright hard X-ray burst
detected with the NICER [26], INTEGRAL [2], Konus-Wind [27], Insight Hard
X-ray Modulation Telescope - HXMT [28] and AGILE satellites [29].

Fast radio bursts are extremely intense millisecond-long radio pulses of
primarily extragalactic origin and showcasing the characteristic dispersion
sweep of radio pulsars. The physical nature of these bursts is unknown, and a
variety of explanations have been considered, including synchrotron maser
emission from young SGRs in supernova remnants (SNRs) [see 30, and references
therein] and magnetospheric emission [see 31, for detais]. SGRs have been
proposed as sources of FRBs [see 32, 33, 34, 35]. In particular, [36] suggest
that repeating FRBs are generated not far from the surface of the source, as a
result of ultrarelativistic internal shocks in the magnetar wind, which are
launched by the magnetospheric flares.

In one of our previous articles we have obtained the contribution to the high
energy and very high-energy gamma-ray (VHE, $E>100$ GeV) emission due to
cosmic-ray acceleration from SNR G57.2+0.8 hosting SGR J1935+2154 with the use
of the GALPROP propagation code222http://galprop.stanford.edu [see 37]. We
consider the SGR J1935+2154 and its environment as a single source near to the
Galactic center. Then, we have proposed that the above setting can provide a
more comprehensive scenario for the generation of GeV-TeV gamma-rays. In fact,
the diffuse gamma-ray emission from the Galactic plane was first detected by
EGRET [38] and then followed by Fermi LAT measurements [39]. Additionally, a
TeV diffuse emission originated in the central part of our Galaxy was detected
by ground-based imaging atmospheric Cherenkov telescopes [40, 41, 42]. The
observations by the High Energy Stereoscopic System (H.E.S.S.) of a gradient
in the cosmic-ray (CR) profile derived from the diffuse VHE gamma-ray
emission, with a peak at the inner regions, imply an injection by a steady
source located at the center of the Milk Way [43].

The main goal of the present paper is to extend our previous study where we
have investigated particle acceleration models for SNR G57.2+0.8 hosting SGR
J1935+2154 [see 37, for details] motivated by H.E.S.S., Fermi-GBM and INTEGRAL
VHE gamma-ray observations [see 44, 45]. Here, we have used GALPROP
propagation code [46, 47, 48] in its enhanced version [49] performing 3D
diffusion/re-acceleration model simulations. Further, we use the new flux from
the pulsar spin-down timescale for SGR J1935+2154. It is assumed that the SGR
J1935+2154 is injecting accelerated particles into the interstellar medium
(ISM) in equal numbers with a fraction of its spin-down power converted to
pairs. After injection, the particles propagate via a diffusive process [see
e.g., 50]. The mechanism we demonstrate here is an indication that improves
the understanding of the leptonic and hadronic origins of the gamma-ray
emission. Also, such analysis would show whether the gamma-ray emission around
this region is due to leptonic or hadronic processes. H.E.S.S. measurements at
VHE - TeV energies of SGR J1935+2154 are crucial for constraining SGR particle
acceleration scenarios up to the knee. As a result, they may provide
independent properties of the magnetospheric regions of SGRs and hence unveil
some of the physics taking place in their outermost regions. To do so this
paper is organized as follows. Section 2 is devoted to an update and brief
summary of H.E.S.S. observations. In Section 3 and Section 4 we present how
the calculations are done and discuss the obtained results.

##  2 High-energy gamma-ray emission from SGR J1935+2154

Follow-up observations of SGR 1935+2154 were carried by H.E.S.S. for 2 hr on
2020 April 28. The observations are coincident with X-ray bursts detected by
INTEGRAL and Fermi-GBM, thus providing the first very high energy gamma-ray
observations of a SGR in a flaring state. High-quality data acquired during
these follow-up observations allow us to perform a search for short-time
transients. No significant signal at energies $E>0.6$ TeV is found, and upper
limits on the persistent and transient emission were derived [see 45, for
details].

The integral upper limits derived from the 2 hr of H.E.S.S. observation,
assuming a spectral index of $E^{-2.5}$, can be translated into upper limits
on the flux $F(E>600$ GeV) $<2.4\times 10^{-12}$ $\rm erg$ $\rm
cm^{-2}s^{-1}$. SGR 1935+2154 is associated with the middle-aged galactic SNR
G57.2+0.8 at a distance of about 6.6 kpc (see [37, 51] and references
therein). Assuming this distance, [45] derive a luminosity upper limit
$L(E>600$ GeV) $<1.3\times 10^{34}$ erg/s. This places constraints on
persistent VHE emission from SGR 1935+2154 during the H.E.S.S. observations.

Moreover, the simultaneous X-ray emission and radio bursts offered the missing
link to correlate FRBs with SGRs. The X-ray emission can be an indication that
protons and electrons are being accelerated at TeV scales via interactions
with matter in the region or through inverse Compton scattering, respectively,
and H.E.S.S. observations can help better understand this scenario. As well
discussed in [37], if SGRs are cosmic-ray sources at the TeV-PeV energy range,
they could contribute to the diffuse gamma-ray emission in the Galaxy given
that it is expected to originate from the interactions of cosmic rays with
Galactic interstellar gas and radiation fields. It is also worth noting that
H.E.S.S. has shown a strong correlation between the Central Molecular Zone
(CMZ) and the brightness distribution of diffuse emission, possibly indicating
that these gamma rays originate from VHE protons interacting with matter in
these regions. It has been argued that the origin of these CRs via Inverse
Compton scattering is likely prevented by the strong radiative losses of TeV
electrons that limit their propagation in the CMZ [52]. Moreover, H.E.S.S. has
observed a high-energy CR density in the CMZ that is an order of magnitude
larger than in the rest of the Galaxy and that could be interpreted by the
existence of more than one accelerating source in the Galactic center [52].
This clearly motivates studies considering G57.2+0.8 and/or SGR J1935+2154 as
such sources. Recently, the Tibet AS$\gamma$ collaboration reported the first
detection of sub-PeV diffuse gamma rays in the Galactic disk. The lack of
correlation between the observed events and known TeV sources virtually rules
out a leptonic origin for this diffuse emission, and argues for an origin by
hadrons accelerated in this region at PeV energies [see 53, for details].

##  3 Simulations

In this section, we briefly review the main features of the model to highlight
the relevant improvements of our simulations compared with previous results
reported in [37]. In order to assess the distribution of gamma rays from the
propagation of cosmic rays (CRs) for SNR G57.2+0.8, which hosts SGR
J1935+2154, we used GALPROP code [46, 47, 48] in its latest version v57 [49]
performing 3D diffusion/re-acceleration model simulations. The corresponding
spectra of gamma rays reaching the Earth taking into account the propagation
effects and interactions with background radiation were determined. For a
given source distribution and boundary conditions, GALPROP solves the
transport equation. Energy loss, fragmentation and decay, convection,
diffusion, and re-acceleration processes are included in the simulations.
Based on the input CR source abundances, GALPROP computes a comprehensive
network of isotope production. A second-order implicit Crank-Nicholson
approach is used for the numerical solution. The spatial boundary conditions
are based on the assumption that free particles escape. The diffusion
coefficient for a given halo size is defined by the B/C nuclear ratio as a
function of rigidity and dispersion parameter [46, 47, 48]. To a better
comprehension, in this paper we improve the model described in [37] with its
environment and GALPROP parameters, used to better describe the SNR G57.2+0.8
hosting SGR J1935+215.

However, the calculation of the spectral energy distribution of the total
gamma-ray emission uses parameters of the source SGR J1935+2154. The modelling
assumes that accelerated particles are injected into the ISM and the spectral
model for the injected particles is obtained with a power law $dq(p)/dp\propto
p^{-\alpha}$, where $\alpha$ is the spectral index and $q(p)$ the injection
energy spectrum and, the total spectrum is normalized so that the total
injected power is given by the expression [54, 50, 49]

$L=L(t)+L^{\prime}\ \ \ \mathrm{and}\ \ \ L(t)=\eta
L_{0}\Biggl{(}1+\frac{t}{\tau_{0}}\Biggl{)}^{-2},$ (3.1)

where ${L_{0}}$ is the initial spin-down power of the SGR, $\eta=1.0$ is the
efficiency factor, $\tau_{0}=3.37$ $\mathrm{yr}$ is the pulsar time scale
defined as the ratio of the initial rotational energy to the initial spin-down
luminosity [see 54] and $L^{\prime}$ is the gamma-rays luminosity of the burst
state from H.E.S.S. upper limits [45]. The initial spin-down power is
calculated using the current spin-down power of $1.7\times 10^{34}$
$\mathrm{erg\;s^{-1}}$ assuming the characteristic age of $3.6\times 10^{3}$
$\mathrm{yr}$333In the magnetic dipole approximation, we can use
$\Omega/\dot{2\Omega}=-(t+\tau)$ [see Eq. A5 of 54], the characteristic age
$t$, and the current values of $\Omega$ and $\dot{\Omega}$ to calculate that
its pulsar time scale $\tau_{0}=3.37$ yr. Assuming NSs canonical parameters
(moment of inertia $I\sim 1.4\times 10^{45}$ $\rm gcm^{2}$) we can derive the
initial rotational energy. The similar approach was developed by [54] to
estimate the initial energy by using the current spin-down luminosity [see
also 50, 49].. The luminosity $L$ represents the total emission of particles
from SGR J1935+2154. We assume that this emission is caused by the rotational
energy loss of the pulsar [see e.g., 54], and by the luminosity in a flaring
state from the magnetar. Then, we have considered the magnetic energy
injection rate (for both quiescent and flaring/outburst emissions). The
distance of the SGR J1935+2154 is $6.6\pm 0.7$ kpc [55, 56, 51].

Calculations are performed in a Cartesian grid, with the Galactic plane
assumed to be the X-Y plane at whose origin the GC is located. The Galactic
volume spans up to 18 kpc in the X and Y directions, with a halo height (h) of
8 kpc in the X-Y direction. We use the treatment of inhomogeneous diffusion
properties in the space around CR sources, modelling two-zone diffusion
scenarios, as described in [37].

The generation of $\gamma$-rays is determined by the propagating CR
distributions, including primary, secondary, and inelastically scattered
protons. Inverse Compton scattering is calculated by the interactions between
anisotropic background photon distribution with a Galactic interstellar
radiation field. As well as different contributions of the pion-decay and
bremsstrahlung are calculated from emission levels using the column densities
of $H$ gas. Models of the Galactic magnetic field are used to determine the
synchrotron emission. The $\gamma$-ray and synchrotron sky maps result from
integrating the relevant emissivities with gas distributions, interstellar
radiation, and magnetic fields [see 46]. The gamma-ray flux is calculated from
G57.2+0.8 hosting SGR J1935+2154 and compared with those published in [37].

##  4 Discussions

Gamma-ray emission can be produced by hadronic and leptonic processes
triggered by the interaction of CRs accelerated in the molecular cloud region
and ISM. We simulate the CR flux as a result of SGR J1935+2154 injecting
particles into the surrounding space and account the contribution of gamma-ray
emission by the propagation of CRs for SNR G57.2+0.8.

In order to fix the same environment as in [37] we choose diffusion constant
$D_{0}=2.5\times 10^{28}$ $\mathrm{cm^{2}\;s^{-1}}$, slope of the diffusion
coefficient $\delta=0.6$, Alfvén velocity $v_{a}=28.0$ $\mathrm{km\;s^{-1}}$
and spectral slope $\alpha=2.2,2.4,2.6$. The results are described in Figure
1. The Figure displays the behavior of the $S0$ model described in [37], in
which the VHE gamma-ray emission from SNR G57.2+0.8 and SGR J1935+2154 was
normalized as the sum of pion decay and inverse Compton spectra, with an upper
limit on the integral flux of TeV gamma rays from H.E.S.S. Observatory up to
99.5% CL [52]. The model presented in [37], although consistent, shows a lower
flux than the current model. This is because the flux calculation in previous
work is associated with the upper limits of the gamma integral measured by
H.E.S.S. [52], which was updated in [45].

Figure 1 shows the data of diffuse TeV energy gamma-ray emission from the
Galactic Center of H.E.S.S. [43, 52] and MAGIC models [42]. According to our
findings, gamma-rays of CRs from SGR J1935+2154 can make contribution to the
total gamma energy density distribution and, in special, to the diffusion
gamma ray from Galactic Center. Up to $10^{4}$ $\mathrm{MeV}$, the pion
emission component dominates. The differences in the spectral indices show
that smaller indices make a larger contribution to the Galactic diffuse gamma,
see Fig. 1 \- (a,b,c). The figure 1 also describes the production of the VHE
gamma-ray emission from particles following the method discussed in [37]. The
method makes a conservative connection between the upper limit of the integral
flux of TeV gamma-rays detected by H.E.S.S. [45] and the cosmic-ray luminosity
at the source. The red line on figure 1 shows this result. We suggest that the
flare was the result of the interaction of particles accelerated, producing
$p-p$ collisions and generating gamma-ray emission and, therefore, dominating
the VHE spectrum, figure 1-(b,c). This emission is below than the quiescent
one, described by the orange line, figure 1-(a), because the spectral index is
harder, consequently generating more pions at higher energies.

(a) $\alpha=2.2$ (b) $\alpha=2.4$

(c) $\alpha=2.6$

Figure 1: Spectral energy distribution of the gamma-ray emission. The total
gamma-ray spectrum is the sum of pion decay, inverse Compton, and
bremsstrahlung from the spin-down model, $L(t)$, and the H.E.S.S. 2021 model,
described from the flaring state $L^{\prime}$, of the SGR J1935+2154. See text
for details. The models use the 2D gas distribution [48].

This simple model presented is an interesting window that improves the
understanding about the hadronic origins of gamma-ray emission at high
energies, as discussed by H.E.S.S. [45]. The approach presented here becomes
more predictive considering the energy losses during the CRs propagation,
magnetic field, distance and the spin-down age of SGR J1935+2154.

## Acknowledgments

The authors thank the referee for comments which helped to improve the quality
of the manuscript. J.G.C. is grateful for the support of CNPq (311758/2021-5)
and FAPESP Project No. 2021/01089-1. The research of R.C.A. is supported by
Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant
number 310448/2021-2. R.C.A. and L.N.P thank for the support of L’Oreal
Brazil, with partnership of ABC and UNESCO in Brazil. We acknowledge the
National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing
HPC resources of the SDumont supercomputer, which have contributed to the
research results reported within this paper. URL: http://sdumont.lncc.br.
J.G.C and R.C.A acknowledge the financial support of FAPESP Project No.
2021/01089-1. The authors acknowledge the financial support of NAPI “Fenômenos
Extremos do Universo” of Fundação Araucária.

## References

  * [1] G. L. Israel, P. Esposito, N. Rea, F. Coti Zelati, A. Tiengo, S. Campana, S. Mereghetti, G. A. Rodriguez Castillo, D. Götz, M. Burgay, A. Possenti, S. Zane, R. Turolla, R. Perna, G. Cannizzaro, and J. Pons, “The discovery, monitoring and environment of SGR J1935+2154,” MNRAS, vol. 457, pp. 3448–3456, Apr. 2016.
  * [2] S. Mereghetti and et al., “INTEGRAL Discovery of a Burst with Associated Radio Emission from the Magnetar SGR 1935+2154,” ApJ Letters, vol. 898, p. L29, Aug. 2020.
  * [3] A. Borghese, F. Coti Zelati, N. Rea, P. Esposito, G. L. Israel, S. Mereghetti, and A. Tiengo, “The X-Ray Reactivation of the Radio Bursting Magnetar SGR J1935+2154,” ApJ Letters, vol. 902, p. L2, Oct. 2020.
  * [4] S. A. Olausen and V. M. Kaspi, “The McGill Magnetar Catalog,” ApJS, vol. 212, p. 6, May 2014.
  * [5] S. Mereghetti, “The strongest cosmic magnets: soft gamma-ray repeaters and anomalous X-ray pulsars,” Astronomy and Astrophysics Reviews, vol. 15, pp. 225–287, July 2008.
  * [6] R. Turolla, S. Zane, and A. L. Watts, “Magnetars: the physics behind observations. A review,” Reports on Progress in Physics, vol. 78, p. 116901, Nov. 2015.
  * [7] V. M. Kaspi and A. M. Beloborodov, “Magnetars,” Annual Review of Astronomy and Astrophysics, vol. 55, pp. 261–301, Aug. 2017.
  * [8] R. C. Duncan and C. Thompson, “Formation of very strongly magnetized neutron stars - Implications for gamma-ray bursts,” ApJ Letters, vol. 392, pp. L9–L13, June 1992.
  * [9] C. Thompson and R. C. Duncan, “The soft gamma repeaters as very strongly magnetized neutron stars - I. Radiative mechanism for outbursts,” MNRAS, vol. 275, pp. 255–300, 7 1995.
  * [10] A. M. Beloborodov and C. Thompson, “Corona of Magnetars,” ApJ, vol. 657, pp. 967–993, Mar. 2007.
  * [11] N. Rea, J. A. Pons, D. F. Torres, and R. Turolla, “The Fundamental Plane for Radio Magnetars,” ApJ Letters, vol. 748, p. L12, Mar. 2012.
  * [12] A. M. Beloborodov, “On the Mechanism of Hard X-Ray Emission from Magnetars,” ApJ, vol. 762, p. 13, Jan. 2013.
  * [13] R. C. R. de Lima, J. G. Coelho, J. P. Pereira, C. V. Rodrigues, and J. A. Rueda, “Evidence for a Multipolar Magnetic Field in SGR J1745-2900 from X-Ray Light-curve Analysis,” ApJ, vol. 889, p. 165, Feb. 2020.
  * [14] J. van Paradijs, R. E. Taam, and E. P. J. van den Heuvel, “On the nature of the ‘anomalous’ 6-s X-ray pulsars,” A&A, vol. 299, p. L41, July 1995.
  * [15] M. A. Alpar, “On Young Neutron Stars as Propellers and Accretors with Conventional Magnetic Fields,” ApJ, vol. 554, pp. 1245–1254, June 2001\.
  * [16] J. G. Coelho, D. L. Cáceres, R. C. R. de Lima, M. Malheiro, J. A. Rueda, and R. Ruffini, “The rotation-powered nature of some soft gamma-ray repeaters and anomalous X-ray pulsars,” A&A, vol. 599, p. A87, Mar. 2017.
  * [17] R. Ouyed, D. Leahy, and B. Niebergal, “SGR 0418+5729 as an evolved Quark-Nova compact remnant,” MNRAS, vol. 415, pp. 1590–1596, Aug. 2011\.
  * [18] M. Malheiro, J. A. Rueda, and R. Ruffini, “SGRs and AXPs as Rotation-Powered Massive White Dwarfs,” Publications of the Astronomical Society of Japan, vol. 64, p. 56, 6 2012.
  * [19] J. G. Coelho and M. Malheiro, “Magnetic Dipole Moment of SGRs and AXPs Described as Massive and Magnetic White Dwarfs,” Publications of the Astronomical Society of Japan, vol. 66, pp. 14–14, 2 2014.
  * [20] R. V. Lobato, M. Malheiro, and J. G. Coelho, “Magnetars and white dwarf pulsars,” International Journal of Modern Physics D, vol. 25, no. 09, p. 1641025, 2016.
  * [21] D. L. Cáceres, S. M. de Carvalho, J. G. Coelho, R. C. R. de Lima, and J. A. Rueda, “Thermal X-ray emission from massive, fast rotating, highly magnetized white dwarfs,” MNRAS, vol. 465, pp. 4434–4440, Mar. 2017.
  * [22] S. V. Borges, C. V. Rodrigues, J. G. Coelho, M. Malheiro, and M. Castro, “A Magnetic White Dwarf Accretion Model for the Anomalous X-Ray Pulsar 4U 0142+61,” ApJ, vol. 895, p. 26, May 2020.
  * [23] M. F. Sousa, J. G. Coelho, and J. C. N. de Araujo, “Gravitational waves from SGRs and AXPs as fast-spinning white dwarfs,” MNRAS, vol. 498, pp. 4426–4432, Nov. 2020.
  * [24] CHIME/FRB Collaboration, B. C. Andersen, K. M. Bandura, M. Bhardwaj, A. Bij, M. M. Boyce, P. J. Boyle, C. Brar, T. Cassanelli, P. Chawla, T. Chen, J. F. Cliche, A. Cook, D. Cubranic, A. P. Curtin, N. T. Denman, M. Dobbs, F. Q. Dong, M. Fandino, E. Fonseca, B. M. Gaensler, U. Giri, D. C. Good, M. Halpern, A. S. Hill, G. F. Hinshaw, C. Höfer, A. Josephy, J. W. Kania, V. M. Kaspi, T. L. Landecker, C. Leung, D. Z. Li, H. H. Lin, K. W. Masui, R. McKinven, J. Mena-Parra, M. Merryfield, B. W. Meyers, D. Michilli, N. Milutinovic, A. Mirhosseini, M. Münchmeyer, A. Naidu, L. B. Newburgh, C. Ng, C. Patel, U. L. Pen, T. Pinsonneault-Marotte, Z. Pleunis, B. M. Quine, M. Rafiei-Ravandi, M. Rahman, S. M. Ransom, A. Renard, P. Sanghavi, P. Scholz, J. R. Shaw, K. Shin, S. R. Siegel, S. Singh, R. J. Smegal, K. M. Smith, I. H. Stairs, C. M. Tan, S. P. Tendulkar, I. Tretyakov, K. Vanderlinde, H. Wang, D. Wulf, and A. V. Zwaniga, “A bright millisecond-duration radio burst from a Galactic magnetar,” Nature, vol. 587, pp. 54–58, Nov. 2020.
  * [25] C. Bochenek, S. Kulkarni, V. Ravi, D. McKenna, G. Hallinan, and K. Belov, “Independent detection of the radio burst reported in ATel #13681 with STARE2,” The Astronomer’s Telegram, vol. 13684, p. 1, Apr. 2020.
  * [26] G. Younes, T. Güver, C. Kouveliotou, M. G. Baring, C.-P. Hu, Z. Wadiasingh, B. Begiçarslan, T. Enoto, E. Göğüş, L. Lin, A. K. Harding, A. J. van der Horst, W. A. Majid, S. Guillot, and C. Malacaria, “NICER View of the 2020 Burst Storm and Persistent Emission of SGR 1935+2154,” ApJ Letters, vol. 904, p. L21, Dec. 2020.
  * [27] A. Ridnaia, D. Svinkin, D. Frederiks, A. Bykov, S. Popov, R. Aptekar, S. Golenetskii, A. Lysenko, A. Tsvetkova, M. Ulanov, and T. L. Cline, “A peculiar hard X-ray counterpart of a Galactic fast radio burst,” Nature Astronomy, vol. 5, pp. 372–377, Apr. 2021.
  * [28] C. K. Li, L. Lin, S. L. Xiong, M. Y. Ge, X. B. Li, T. P. Li, F. J. Lu, S. N. Zhang, Y. L. Tuo, Y. Nang, B. Zhang, S. Xiao, Y. Chen, L. M. Song, Y. P. Xu, C. Z. Liu, S. M. Jia, X. L. Cao, J. L. Qu, S. Zhang, Y. D. Gu, J. Y. Liao, X. F. Zhao, Y. Tan, J. Y. Nie, H. S. Zhao, S. J. Zheng, Y. G. Zheng, Q. Luo, C. Cai, B. Li, W. C. Xue, Q. C. Bu, Z. Chang, G. Chen, L. Chen, T. X. Chen, Y. B. Chen, Y. P. Chen, W. Cui, W. W. Cui, J. K. Deng, Y. W. Dong, Y. Y. Du, M. X. Fu, G. H. Gao, H. Gao, M. Gao, Y. D. Gu, J. Guan, C. C. Guo, D. W. Han, Y. Huang, J. Huo, L. H. Jiang, W. C. Jiang, J. Jin, Y. J. Jin, L. D. Kong, G. Li, M. S. Li, W. Li, X. Li, X. F. Li, Y. G. Li, Z. W. Li, X. H. Liang, B. S. Liu, G. Q. Liu, H. W. Liu, X. J. Liu, Y. N. Liu, B. Lu, X. F. Lu, T. Luo, X. Ma, B. Meng, G. Ou, N. Sai, R. C. Shang, X. Y. Song, L. Sun, L. Tao, C. Wang, G. F. Wang, J. Wang, W. S. Wang, Y. S. Wang, X. Y. Wen, B. B. Wu, B. Y. Wu, M. Wu, G. C. Xiao, H. Xu, J. W. Yang, S. Yang, Y. J. Yang, Y.-J. Yang, Q. B. Yi, Q. Q. Yin, Y. You, A. M. Zhang, C. M. Zhang, F. Zhang, H. M. Zhang, J. Zhang, T. Zhang, W. Zhang, W. C. Zhang, W. Z. Zhang, Y. Zhang, Y. Zhang, Y. F. Zhang, Y. J. Zhang, Z. Zhang, Z. Zhang, Z. L. Zhang, D. K. Zhou, J. F. Zhou, Y. Zhu, Y. X. Zhu, and R. L. Zhuang, “HXMT identification of a non-thermal X-ray burst from SGR J1935+2154 and with FRB 200428,” Nature Astronomy, vol. 5, pp. 378–384, Apr. 2021.
  * [29] M. Tavani, C. Casentini, A. Ursi, F. Verrecchia, A. Addis, L. A. Antonelli, A. Argan, G. Barbiellini, L. Baroncelli, G. Bernardi, G. Bianchi, A. Bulgarelli, P. Caraveo, M. Cardillo, P. W. Cattaneo, A. W. Chen, E. Costa, E. Del Monte, G. Di Cocco, G. Di Persio, I. Donnarumma, Y. Evangelista, M. Feroci, A. Ferrari, V. Fioretti, F. Fuschino, M. Galli, F. Gianotti, A. Giuliani, C. Labanti, F. Lazzarotto, P. Lipari, F. Longo, F. Lucarelli, A. Magro, M. Marisaldi, S. Mereghetti, E. Morelli, A. Morselli, G. Naldi, L. Pacciani, N. Parmiggiani, F. Paoletti, A. Pellizzoni, M. Perri, F. Perotti, G. Piano, P. Picozza, M. Pilia, C. Pittori, S. Puccetti, G. Pupillo, M. Rapisarda, A. Rappoldi, A. Rubini, G. Setti, P. Soffitta, M. Trifoglio, A. Trois, S. Vercellone, V. Vittorini, P. Giommi, and F. D’Amico, “An X-ray burst from a magnetar enlightening the mechanism of fast radio bursts,” Nature Astronomy, vol. 5, pp. 401–407, Apr. 2021.
  * [30] B. Zhang, “The physical mechanisms of fast radio bursts,” Nature, vol. 587, pp. 45–53, Nov. 2020.
  * [31] W. Lu, P. Kumar, and B. Zhang, “A unified picture of Galactic and cosmological fast radio bursts,” MNRAS, vol. 498, pp. 1397–1405, Oct. 2020.
  * [32] S. B. Popov and K. A. Postnov, “Hyperflares of SGRs as an engine for millisecond extragalactic radio bursts,” in Evolution of Cosmic Objects through their Physical Activity (H. A. Harutyunian, A. M. Mickaelian, and Y. Terzian, eds.), pp. 129–132, Nov. 2010.
  * [33] Y. Lyubarsky, “A model for fast extragalactic radio bursts.,” MNRAS, vol. 442, pp. L9–L13, July 2014.
  * [34] B. D. Metzger, B. Margalit, and L. Sironi, “Fast radio bursts as synchrotron maser emission from decelerating relativistic blast waves,” MNRAS, vol. 485, pp. 4091–4106, May 2019.
  * [35] A. M. Beloborodov, “Blast Waves from Magnetar Flares and Fast Radio Bursts,” ApJ, vol. 896, p. 142, June 2020.
  * [36] A. M. Beloborodov, “A Flaring Magnetar in FRB 121102?,” ApJ Letters, vol. 843, p. L26, July 2017.
  * [37] R. C. dos Anjos, J. G. Coelho, J. P. Pereira, and F. Catalani, “High-energy gamma-ray emission from SNR G57.2+0.8 hosting SGR J1935+2154,” JCAP, vol. 2021, p. 023, Oct. 2021.
  * [38] S. D. Hunter, D. L. Bertsch, J. R. Catelli, T. M. Dame, S. W. Digel, B. L. Dingus, J. A. Esposito, C. E. Fichtel, R. C. Hartman, G. Kanbach, D. A. Kniffen, Y. C. Lin, H. A. Mayer-Hasselwander, P. F. Michelson, C. von Montigny, R. Mukherjee, P. L. Nolan, E. Schneid, P. Sreekumar, P. Thaddeus, and D. J. Thompson, “EGRET observations of the diffuse gamma-ray emission from the galactic plane,” The Astrophysical Journal, vol. 481, pp. 205–240, may 1997.
  * [39] M. Ackermann and et al., “Fermi-LAT Observations of the Diffuse $\gamma$-Ray Emission: Implications for Cosmic Rays and the Interstellar Medium,” ApJ, vol. 750, p. 3, May 2012.
  * [40] F. Aharonian and et al., “Discovery of very-high-energy $\gamma$-rays from the Galactic Centre ridge,” Nature, vol. 439, pp. 695–698, Feb. 2006.
  * [41] A. Archer and et al., “TeV Gamma-Ray Observations of the Galactic Center Ridge by VERITAS,” ApJ, vol. 821, p. 129, Apr. 2016.
  * [42] MAGIC Collaboration, “MAGIC observations of the diffuse $\gamma$-ray emission in the vicinity of the Galactic center,” A&A, vol. 642, p. A190, Oct. 2020.
  * [43] HESS Collaboration, “Acceleration of petaelectronvolt protons in the Galactic Centre,” Nature, vol. 531, pp. 476–479, Mar. 2016.
  * [44] L. Lin, E. Göğüş, O. J. Roberts, M. G. Baring, C. Kouveliotou, Y. Kaneko, A. J. van der Horst, and G. Younes, “Fermi/GBM View of the 2019 and 2020 Burst Active Episodes of SGR J1935+2154,” ApJ Letters, vol. 902, p. L43, Oct. 2020.
  * [45] H. E. S. S. Collaboration, “Searching for TeV Gamma-Ray Emission from SGR 1935+2154 during Its 2020 X-Ray and Radio Bursting Phase,” ApJ, vol. 919, p. 106, Oct. 2021.
  * [46] A. W. Strong and I. V. Moskalenko, “Propagation of cosmic-ray nucleons in the galaxy,” The Astrophysical Journal, vol. 509, pp. 212–228, dec 1998.
  * [47] T. A. Porter, G. Jóhannesson, and I. V. Moskalenko, “High-energy gamma rays from the milky way: Three-dimensional spatial models for the cosmic-ray and radiation field densities in the interstellar medium,” The Astrophysical Journal, vol. 846, p. 67, aug 2017.
  * [48] G. Jóhannesson, T. A. Porter, and I. V. Moskalenko, “The three-dimensional spatial distribution of interstellar gas in the milky way: Implications for cosmic rays and high-energy gamma-ray emissions,” vol. 856, p. 45, mar 2018.
  * [49] T. A. Porter, G. Johannesson, and I. V. Moskalenko, “The GALPROP Cosmic-ray Propagation and Non-thermal Emissions Framework: Release v57,” arXiv e-prints, p. arXiv:2112.12745, Dec. 2021.
  * [50] G. Jóhannesson, T. A. Porter, and I. V. Moskalenko, “Cosmic-Ray Propagation in Light of the Recent Observation of Geminga,” ApJ, vol. 879, p. 91, July 2019.
  * [51] P. Zhou, X. Zhou, Y. Chen, J.-S. Wang, J. Vink, and Y. Wang, “Revisiting the Distance, Environment, and Supernova Properties of SNR G57.2+0.8 that Hosts SGR 1935+2154,” ApJ, vol. 905, p. 99, Dec. 2020\.
  * [52] H. E. S. S. Collaboration, “Characterising the VHE diffuse emission in the central 200 parsecs of our Galaxy with H.E.S.S.,” A&A, vol. 612, p. A9, Apr. 2018.
  * [53] M. Amenomori and et al., “First Detection of sub-PeV Diffuse Gamma Rays from the Galactic Disk: Evidence for Ubiquitous Galactic Cosmic Rays beyond PeV Energies,” Physical Review Letters, vol. 126, p. 141101, Apr. 2021\.
  * [54] D. Malyshev, I. Cholis, and J. Gelfand, “Pulsars versus dark matter interpretation of ATIC/PAMELA,” PRD, vol. 80, p. 063005, Sept. 2009.
  * [55] X. H. Sun, P. Reich, W. Reich, L. Xiao, X. Y. Gao, and J. L. Han, “A Sino-German $\lambda$6 cm polarization survey of the Galactic plane. VII. Small supernova remnants,” A&A, vol. 536, p. A83, Dec. 2011.
  * [56] M. Z. Pavlović, D. Urošević, B. Vukotić, B. Arbutina, and Ü. D. Göker, “The Radio Surface-brightness-to-Diameter Relation for Galactic Supernova Remnants: Sample Selection and Robust Analysis with Various Fitting Offsets,” Astrophysical Journal Supplement Series, vol. 204, p. 4, Jan. 2013.


# Supervised learning of short and high-dimensional temporal sequences for
life science measurements

F.-M. Schleif 1 A. Gisbrecht 1 B. Hammer 1

(1Univ. of Bielefeld, CITEC Center of Excellence,
Universitätsstrasse 21-23, 33615 Bielefeld, Germany
01\. September 2011
Technical report follow of the Dagstuhl Seminar 11341
Learning in the context of very high dimensional data
21.08.11 - 26.08.11
Organizer: Michael Biehl (Univ. of Groningen, NL), Barbara Hammer (Univ.
Bielefeld, DE), Erzsébet Merényi (Rice Univ., US),
Alessandro Sperduti (Univ. of Padova, IT), Thomas Villmann (Univ. of Applied
Sc. Mittweida, DE))

###### Abstract

Motivation: The analysis of physiological processes over time is becoming
increasingly important. The measurements are often given by spectrometric or
gene expression profiles over time with only few time points but a large
number of measured variables. The analysis of such temporal sequences is
challenging and only few methods have been proposed. The information can be
encoded time independent, by means of classical expression differences for a
single time point or in expression profiles over time. Available methods are
limited to unsupervised and semi-supervised settings. The predictive variables
can be identified only by means of wrapper or post-processing techniques. This
is complicated due to the small number of samples for such studies. Here, we
present a supervised learning approach, termed _Supervised Topographic Mapping
Through Time_ (SGTM-TT). It learns a supervised mapping of the temporal
sequences onto a low dimensional grid. We utilize a hidden markov model (HMM)
to account for the time domain and relevance learning to identify the relevant
feature dimensions most predictive over time. The learned mapping can be used
to visualize the temporal sequences and to predict the class of a new
sequence. The relevance learning permits the identification of discriminating
masses or gen expressions and prunes dimensions which are unnecessary for the
classification task or encode mainly noise. In this way we obtain a very
efficient learning system for temporal sequences.

Results: The results indicate that using simultaneous supervised learning and
metric adaptation significantly improves the prediction accuracy for
synthetically and real life data in comparison to the standard techniques. The
discriminating features, identified by relevance learning, compare favorably
with the results of alternative methods. Our method permits the visualization
of the data on a low dimensional grid, highlighting the observed temporal
structure.

Contact<EMAIL_ADDRESS>
Keywords: high-dimensional time series, short time series, prototype learning,
relevance learning, topographic mapping

##  1 Introduction

The analysis of high-dimensional, short time series, or temporal sequences is
a challenging task. On the one hand side the data are not any longer identical
and independent distributed (i.i.d) due to the time constraint, on the other
hand the dimensionality of the data is large, complicating the learning of a
predictive model. Standard time series methods like auto-regressive moving
average (ARMA) or extensions thereof (see e.g. [9]) are in general not
applicable due to the limited number of time points and the large
dimensionality of the data. Some methods have been proposed to model this type
of data. In [20] an unsupervised projection techniques was proposed employing
a so called temporal context. The temporal data have been processed by a kind
of Self Organizing Map (SOM) [11] but the learning was modified such that it
depends on the the current temporal context. A further unsupervised proposal
has been made in [14] using the Generative Topographic Mapping Through Time
(GTM-TT) ([3]). Some new hidden variables were introduced to account for the
relevance of the different feature dimensions, to accounts, in a non-
discriminative manner, for the explained variance in the data over time. A
supervised two-class method solely based on hidden markov models was proposed
in [13]. It models the two different data distribution by independent HMMs and
evaluates the generated models to obtain a ranking of the input dimensions.
Subsequently the model was improved by selecting a set of features using a
wrapper strategy. In [6] a similar approach was proposed but in a semi-
supervised scenario, introducing classwise constraints in the hidden markov
model. The importance of the individual features was determined using a
complex post processing procedure. Another supervised method using all
features, based on Support Vector Machine (SVM) and a Kalman filter was
proposed in [5].

While the first two approaches have been found to be very effective for
unsupervised analysis, the last mentioned methods focus on supervised and
semi-supervised analysis. The results in [13] are very promising, with $85\%$
prediction accuracy on a real life multiple sclerosis data (MS) set, but make
strong pre-assumptions about the underlying HMM structure. Also, it is
proposed for two class scenarios, only. The approach in [5] improved this
result by $2-5\%$ but in a black box scenario, without additional feature
selection. The approach in [6] is evaluated also with respect to the results
of [13] achieving improved performance for the same MS data sets. There is
still ongoing work of research in this field, reflecting the high demand for
effective methods dealing with this type of data. The application field is not
limited to the bio-medical domain as considered in [13, 6, 8] but covers a
broader field of applications also in industry and geo-science as reflected in
[14, 20].

The identification of the relevant input dimensions of a temporal sequence is
very important as outlined in [14, 13] to obtain better understanding of the
data, to reduce the processing complexity and to improve the overall
prediction accuracy. As already motivated by some of the prior references,
prototype methods (see e.g. [11]) have been found to be very effective for the
analysis of high dimensional data also to analyze temporal sequences. In [3],
the _Generative Topographic Mapping - through time_ (GTM-TT), an unsupervised
prototype based method for the topographic projection of high-dimensional,
temporal sequences was proposed. GTM-TT learns a hidden markov model (HMM) of
a data generating process and represents the data by a prototype based
representation in time and space. Like in ordinary prototype methods the GTM-
TT approximates the data distribution by a vector quantization of the data
space. The temporal dependence between the prototype is modeled by an
appropriate HMM. Additionally the prototypes are assigned to a fixed grid
representation or lattice, which permits, provided the topology is preserved
(see [22]), the easy visualization and interpretation of the data trajectory
in a low dimensional space. In this contribution we extend the GTM-TT to a
supervised method and integrate relevance learning to identify the relevant
dimensions over time. Then we will briefly review Generative Topographic
Mapping (GTM) and Generative Topographic Mapping Through Time. Subsequently,
we outline our method and apply and discuss it for different experimental
data. The paper is closed with links to further extensions and open questions.

##  2 Approach and Methods

###  2.1 Generative Topographic Mapping

The Generative Topographic Mapping (GTM) as introduced in [4] models data
$\mathbf{x}\in\mathbb{R}^{D}$ by means of a mixture of Gaussians which is
induced by a lattice of points $\mathbf{w}$ in a low dimensional latent space
which can be used for visualization.

Figure 1: GTM-TT consisting of a HMM in which the hidden states are given by
the latent points of the GTM model. The emission probabilities are governed by
the GTM mixture distribution [3]. The different data distributions,
exemplified in 3D (bottom) and indicated by the color/shading are mapped to
the 2D grid (top). Here we consider $9$ hidden states on a $3\times 3$ grid.
The data distribution may change over time and hence also the mapping of the
GTM is effected over time, assuming smooth transitions.

The lattice points are mapped via
$\mathbf{w}\mapsto\mathbf{t}=y(\mathbf{w},\mathbf{W})$ to the data space,
where the function is parametrized by $\mathbf{W}$; one can, for example, pick
a generalized linear regression model based on Gaussian base functions

$y:\mathbf{w}\mapsto\Phi(\mathbf{w})\cdot\mathbf{W}$ (1)

where the base functions $\Phi$ are equally spaced Gaussians.The high-
dimensional points $y$ are so called prototypes of the original data space,
representing a larger set of points, they are responsible for, as measured by
Eq. (5). They can be directly inspected and permit to summarize the data.
Every latent point induces a Gaussian

$p(\mathbf{x}|\mathbf{w},\mathbf{W},\beta)=\left(\frac{\beta}{2\pi}\right)^{D/2}\exp\left(-\frac{\beta}{2}\|\mathbf{x}-y(\mathbf{w},\mathbf{W})\|^{2}\right)$
(2)

with variance $\beta^{-1}$, which gives the data distribution as mixture of
$K$ modes

$p(\mathbf{x}|\mathbf{W},\beta)=\sum_{k=1}^{K}p(\mathbf{w}^{k})p(\mathbf{x}|\mathbf{w}^{k},\mathbf{W},\beta)$
(3)

where, usually, $p(\mathbf{w}^{k})$ is taken as Dirac distributions of the
prototypes. Training of GTM optimizes the data log-likelihood

$\ln\left(\prod_{n=1}^{N}\left(\sum_{k=1}^{K}p(\mathbf{w}^{k})p(\mathbf{x}^{n}|\mathbf{w}^{k},\mathbf{W},\beta)\right)\right)$
(4)

by means of an expectation maximization (EM) approach with respect to the
parameters $\mathbf{W}$ and $\beta$. In the E step, the responsibility of
mixture component $k$ for point $n$ is determined as

$r^{kn}=p(\mathbf{w}^{k}|\mathbf{x}^{n},\mathbf{W},\beta)=\frac{p(\mathbf{x}^{n}|\mathbf{w}^{k},\mathbf{W},\beta)p(\mathbf{w}^{k})}{\sum_{k^{\prime}}p(\mathbf{x}^{n}|\mathbf{w}^{k^{\prime}},\mathbf{W},\beta)p(\mathbf{w}^{k^{\prime}})}$
(5)

In the M step, the weights $\mathbf{W}$ are determined solving the equality

$\mathbf{\Phi}^{T}\mathbf{G}_{\mathrm{old}}\mathbf{\Phi}\mathbf{W}_{\mathrm{new}}^{T}=\mathbf{\Phi}^{T}\mathbf{R}_{\mathrm{old}}\mathbf{X}$
(6)

where $\mathbf{\Phi}$ refers to the matrix of base functions $\Phi$ evaluated
at points $\mathbf{w}^{k}$, $\mathbf{X}$ to the data points, $\mathbf{R}$ to
the responsibilities, and $\mathbf{G}$ is a diagonal matrix with accumulated
responsibilities $G_{nn}=\sum_{k}r^{kn}(\mathbf{W},\beta)$. The variance can
be computed by solving

$\frac{1}{\beta_{\mathrm{new}}}=\frac{1}{ND}\sum_{k,n}r^{kn}(\mathbf{W}_{\mathrm{old}},\beta_{\mathrm{old}})\|{\Phi}(\mathbf{w}^{k})\mathbf{W}_{\mathrm{new}}-\mathbf{x}^{n}\|^{2}$
(7)

where $D$ is the data dimensionality and $N$ the number of data points.

###  2.2 Relevance learning

The principle of relevance learning has been introduced in [10] as a
particularly simple and efficient method to adapt the metric of prototype
based classifiers according to the given situation at hand. It takes into
account a relevance scheme of the data dimensions by substituting the squared
Euclidean metric by the weighted form

$d_{\boldsymbol{\lambda}}(\mathbf{x},\mathbf{t})=\sum_{d=1}^{D}\lambda_{d}^{2}(x_{d}-t_{d})^{2}\,.$
(8)

The principle is extended in [18, 17] to the more general metric form

$d_{\boldsymbol{\Omega}}(\mathbf{x},\mathbf{t})=(\mathbf{x}-\mathbf{t})^{T}\boldsymbol{\Omega}^{T}\boldsymbol{\Omega}(\mathbf{x}-\mathbf{t})$
(9)

Using a square matrix $\boldsymbol{\Omega}$, a positive semi-definite matrix
which gives rise to a valid pseudo-metric is achieved this way. In [18, 17],
these metrics are considered in local and global form, i.e. the adaptive
metric parameters can be identical for the full model, or they can be attached
to every prototype present in the model. Relevance learning for GTM has been
already introduced in [7] for i.i.d. data. In case of temporal sequences some
modification of the original principle are necessary and also the supervision
will be handled differently as pointed out subsequently. First however we
review the GTM through time as described in [3, 15] which is the basic method
to process i.i.d. data in our approach.

###  2.3 Generative Topographic Mapping Through-Time

The GTM through time (GTM-TT) has been introduced in [3]. For data vectors
$\mathbf{x}$ which have the form of a time series the vectors are no longer
independent. Nearby timepoints are likely to be correlated. As pointed out in
[3] such effects can be captured using Hidden Markov Models (HMM). Accordingly
in [3] the GTM is equipped by a HMM, constructing a kind of a topology-
constrained HMM

The structure of the GTM-TT is shown in Figure 1. Assuming a sequence length
$T$, of hidden states $Z=\\{z_{1},\ldots,z_{n},\ldots z_{T}\\}$ and the
observed multidimensional time series $X=\\{x^{1},x^{2},\ldots,x^{n},\ldots
x^{T}\\}$, the probability of the observations is given by

$p(X)=\sum_{\text{all sequences }Z}p(Z,X)$ (10)

where $p(Z,X)$ defines the complete data likelihood as in HMM models [4]
taking the following form:

$p(Z,X)=p(z_{1})\prod_{n=2}^{T}p(z_{n}|z_{n-1})\prod_{n=1}^{T}p(x^{n}|z_{n})$
(11)

So it consists of the initial state probability, the transition probability
between two hidden states, capturing the temporal dependence, and the
probability to observe a specific sequence in a given state also known as
emission probability (covered by Eq. (2)). The model parameters are
$\Theta=(\pi,A,W,\beta)$ where $\pi=\\{\pi_{j}\\}:\pi_{j}:=p(z_{1}=j)$ are the
initial state probabilities. $A=\\{a_{ij}\\}:a_{ij}=p(z_{n}=j|z_{n-1}=i)$ are
the transition state probabilities, and $\\{W,\beta\\}$ are given by Eq. (6).
Again we control the gaussians by a common invariance $\beta$. As in HMM the
above likelihood can be efficiently calculated using the _forward backward
procedure_ [23]. The probability being in state $\mathbf{w}_{k}$ at time $n$,
given the observation sequence and the model, also known as responsibility
$r^{kn}$ is calculated as:

$r^{kn}=p(z_{n}=\mathbf{w}^{k}|X,\Theta)=\frac{A_{kn}B_{kn}}{p(X|\Theta)}$
(12)

The forward variable $A_{kn}$ is the joint probability of the past sequences
$\\{x^{1},\ldots,x^{n}\\}$ and the state $z_{n}=\mathbf{w}^{k}$, i.e.
$A_{kn}=p(\\{x^{1},\ldots,x^{n}\\},z_{n}=\mathbf{w}^{k}|\Theta)$, given by the
following recursive equation:

$A_{kn}=\left(\sum_{i=1}^{K}A_{i,n-1}p_{i,k}\right)p_{k}(x^{n})$ (13)

where $A_{k,1}=\pi_{k}p_{k}(x^{1})$. The backward variable $B_{kn}$ which is
the probability of the future sequence $x^{n+1},x^{n+2},\ldots,x^{N}$ given
the hidden state $z_{n}=\mathbf{w}^{k}$, i.e.
$B_{kn}=p(\\{x^{n+1},x^{n+2},\ldots,x^{N}\\}|z_{n}=\mathbf{w}^{k},\Theta)$ is
calculated using the following recursive equation:

$B_{kn}=\sum_{i=1}^{K}p_{i,k}p_{i}(x^{n+1})B_{i,n+1}$ (14)

where $B_{k,T}=1$. The whole parameter estimation can be accomplished by a
maximum likelihood optimization using the EM algorithm as sketched above.
Details can be found in [19].

###  2.4 Supervised GTM-TT

Assume that data point $X$ is equipped with label information $l$ which is
element of a finite set of different labels $L$, e.g. $L=\\{0,1\\}$. Lets
assume we have only two labels 111An extension to multiple labels is straight
forward.. The data are divided into two groups, according to the labeling and
we train one separate GTM-TT per group. To keep the models comparable, the
$\beta$ update for the models is linked to each other and optimized in the
outer loop. The parameters $\mathbf{W}$ are determined for each model
individually leading to $\mathbf{W}_{0}$ and $\mathbf{W}_{1}$. We will further
assume that the grid structure is common for both models. The learning
procedure is thus similar to GTM-TT and depicted in Figure 1.

1:function Supervised GTM-TT($X$,$L$,$K$)

2: [Xn,Pars ] = normalize(X)

3: [X1,X2,L1,L2] = splitdata(Xn,L)

4: [$M_{0}$, $M_{1}$] = init both GTM-TT models

5: repeat

6: call train_single_step for $M_{0}$, $M_{1}$

7: call convergence_check for $M_{0}$, $M_{1}$

8: call optimize_beta for $M_{0}$, $M_{1}$

9: $\bar{\beta}=$ calculate mean of the  $\beta$

10: call update_beta($M_{0}$, $M_{1}$,$\bar{\beta}$)

11: until convergence is true for both models

12:end function

Algorithm 1 Pseudocode of supervised GTM-TT

We denote the obtained model as Supervised GTM-TT (SGTM-TT) and the submodels
as $M_{0}$ and $M_{1}$. The concept of the SGTM-TT is depicted schematically
in Figure 2.

Figure 2: Illustration of the SGTM-TT. It consists of multiple GTM-TT models.
It behaves similar to the regular GTM-TT but the training is classwise and the
$\beta$ parameter is common over the different models. The different classwise
models (top) are used to represent the data distribution (bottom) over time
(from left to right).

###  2.5 Classification using SGTM-TT

To classify new data points with the SGTM-TT model different approaches can be
taken. The simplest one is to make direct use of the samplewise likelihoods
considering the class wise models. In that case a new point is assigned to the
model with maximal likelihood considering one model against the rest. A more
interesting approach is to combine the performance of the generative SGTM-TT
model with a discriminative approach like the SVM [21]. Again we use the
likelihood values from the forward procedure (13) of the SGTM-TT and define a
kernel as follows:

$\displaystyle Lik_{j}^{l}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{K}A_{i,j}\text{ for a series $j$ and a sub-model }l$
(15) $\displaystyle K(X_{j},X_{k})$ $\displaystyle=$
$\displaystyle\sum_{l=1:\\#L}Lik_{j}^{l}\cdot Lik_{k}^{l}\;\text{ with equal
prior}$ (16)

Hence the kernel $K$ is based on a kernel function of inner-products in a one
dimensional feature space of the likelihood-values. In the following we will
make use of this approach employing a standard SVM implementation.

###  2.6 Relevance learning for SGTM-TT

Relevance learning for GTM has been introduced in [7], as the Relevance GTM
(R-GTM). The basic idea for Relevance GTM is to introduce an adaptive metric
for the GTM. The original Euclidean metric is replaced by a parametric
distance like the weighted Euclidean metric (8). After each GTM training step
the prototypes are post-labeled according to its responsibilities, employing
the labeling $L$ of the datapoints. Subsequently the metric parameters of the
distance are adapted according to an optimization criterion. In the article of
[7] different cost functions E where suggested.

The data of the GTM-TT are not any longer i.i.d. and, as mentioned before we
observe a sequence of states $Z$ for a given time series $X$. In the SGTM-TT
we know the labeling of the prototypes, assuming constant labels over time,
due to the split of the learning problem according to the data labeling.
Further, using a common metric and common $\beta$ parameters the prototypes
$Y$ exist still in the same common dataspace. Relevance learning can now be
done in the same way as for R-GTM. This however is often not useful because
the original relevance learning ignores the time domain. If data separation is
observed over time and not for a single time point the R-GTM approach will
fail. For temporal sequences we may also be interested on two views of
relevance, namely relevant, or separating input _dimensions_ $x_{i}$ but also
relevant _time points_ in a temporal sequence $x$. Taking this problem into
account we consider two distance measures, one for the time domain, denoted as
$d^{t}$ and one for the time-independent data space $d^{\lambda}$. A
parametrization of $d^{t}$ can be used to account for the relevance of
specific time points, e.g to prune out time points which are irrelevant for
the representation of the data in a discriminative manner. Parameters on
$d^{\lambda}$ can be used to identify discriminating feature dimensions, e.g.
to prune out noisy dimensions. Subsequently, we provide a distance measure
which can be used for $d^{t}$ and a specific form for $d^{\lambda}$. For
simplicity we will use a simple _global_ , _diagonal_ metric learning scheme
in the experiments.

SGTM-TT provides a probabilistic prediction of the internal representation
$\hat{\mathbf{x}}$ of a time series $\mathbf{x}$ considering the two GTM-TT
models, we obtain one reconstruction each:

$\displaystyle\hat{\mathbf{x}^{n}}^{l}_{i}$ $\displaystyle=$ $\displaystyle
Y_{l}(\underset{k}{\arg\max}\left(r^{kn}\right),i)\;\forall i\in[1,D]$
$\displaystyle\text{with }l\in\\{0,1\\}$

Now, two distances are calculated over time for each point and each dimension
$i$: $d^{t}(\hat{\mathbf{x}^{n}}^{0}_{i},\mathbf{x}^{n}_{i})$,
$d^{t}(\hat{\mathbf{x}^{n}}^{1}_{i},\mathbf{x}^{n}_{i})$. Using one of the
suggested cost functions in the paper of [7] we can calculate the relevance of
the individual dimensions for the separation between the two reconstructions
per point and hence between the different models.

Like for R-GTM the metric adaptation is done by an appropriate optimization
scheme on the cost functions, here we will use stochastic gradient descend,
with a fixed learning rate $\epsilon=0.1$. To avoid convergence to trivial
optima such as zero we pose constraints on the metric parameters of the form
$\|\boldsymbol{\lambda}\|=1$ or
$\mathrm{trace}(\boldsymbol{\Omega}^{T}\boldsymbol{\Omega})^{2}=1$, for matrix
learning. This is achieved by normalization of the values, i.e. after every
gradient step, $\boldsymbol{\lambda}$ is divided by its length, and
$\boldsymbol{\Omega}$ is divided by the square root of
$\mathrm{trace}(\boldsymbol{\Omega}^{T}\boldsymbol{\Omega})$.

A pseudo code of the SGTM-TT with relevance learning is depicted in 2.

1:function SGTM-TT-R($X$,$L$,$K$)

2: [Xn,Pars ] = normalize(X)

3: [X1,X2,L1,L2] = splitdata(Xn,L)

4: initialize the common metric

5: [$M_{0}$,$M_{1}$] = init both GTM-TT models

6: repeat

7: call train_single_step for each GTM-TT model

8: call convergence_check for each GTM-TT model

9: if $cycle>10$ then

10: $\forall X,\forall i=1:D$ call reconstruct($X_{i},M_{0},M_{1}$)

11: $\forall X,\forall i=1:D$ call $d^{t}(M_{0},\hat{x_{i}}_{0},x_{i})$

12: $\forall X,\forall i=1:D$ call $d^{t}(M_{1},\hat{x_{i}}_{1},x_{i})$

13: $\forall X$ call calculate_metric_update

14: average the metric updates and normalize

15: update the metric parameter annealed by $\epsilon$

16: end if

17: call optimize_beta for each GTM-TT model

18: $\bar{\beta}=$ calculate mean of the  $\beta$

19: call update_beta(M1,M2,$\bar{\beta}$)

20: until convergence is true for both models

21:end function

Algorithm 2 Pseudocode of supervised GTM-TT with relevance learning

Usually, we alternate between one EM step, one epoch of gradient descent, and
normalization in our experiments and start the metric learning after $10$
epochs of EM learning to allow a reasonable pre-positioning of the GTM-TT in
the dataspace. The metric learning is annealed by $\epsilon$. Since EM
optimization is much faster than gradient descent, this way, we can enforce
that the metric parameters are adapted on a slower time scale. Hence we can
assume an approximately constant metric for the EM optimization, i.e. the EM
scheme optimizes the likelihood as before. Metric adaptation takes place
considering quasi stationary states of the GTM solution due to the slower time
scale. The call of train_single_step is a regular EM optimization step of the
GTM-TT but without the adaptation of the parameter $\beta$ which is postponed
to allow a linking between the two GTM-TT models included in the SGTM-TT.

Now, we briefly review a concrete cost function $E$ of the relevance GTM for
the metric adaptation as already introduced in [7] but account for the
alternative distance calculations mentioned before.

#### Cost function - Generalized Relevance GTM (GRGTM)

Metric parameters have the form $\boldsymbol{\lambda}$ or
$\boldsymbol{\lambda}^{k}$ for a diagonal metric (8) and $\boldsymbol{\Omega}$
or $\boldsymbol{\Omega}^{k}$ for a full matrix (9), depending on whether a
local or global scheme is considered. In the following, we define the general
parameter $\Theta^{k}$ which can be chosen as one of these four possibilities
depending on the given setting. Thereby, we can assume that $\Theta^{k}$ can
be realized by a matrix which has diagonal form (for relevance learning) or
full matrix form (for matrix updates).

The cost function of generalized relevance GTM is taken from generalized
relevance learning vector quantization (GRLVQ), which can be interpreted as
maximizing the hypothesis margin of a prototype based classification scheme
[10, 18]. The cost function has the form

$E(\Theta)=\sum_{n}E_{n}(\Theta)=\sum_{n}\operatorname{sgd}\left(\frac{d_{\Theta^{+}}(\mathbf{x}^{n},\hat{\mathbf{x}^{n}}^{+})-d_{\Theta^{-}}(\mathbf{x}^{n},\hat{\mathbf{x}^{n}}^{-})}{d_{\Theta^{+}}(\mathbf{x}^{n},\hat{\mathbf{x}^{n}}^{+})+d_{\Theta^{-}}(\mathbf{x}^{n},\hat{\mathbf{x}^{n}}^{-})}\right)$
(17)

where $\operatorname{sgd}(x)=(1+\exp(-x))^{-1}$, $\hat{\mathbf{x}^{n}}^{\pm}$
is the reconstruction of $\mathbf{x}^{n}$ over time using the model $M_{0}$ or
$M_{1}$ depending on the label of $\mathbf{x}$, $+$ indicates the model with
the same level $-$ the model with a different label or the model for the
remaining data.

The adaptation formulas can be derived thereof by taking the derivatives with
respect to the metric parameter. Depending on the form of the metric, the
derivative of the metric is simple

$\frac{\partial
d_{\boldsymbol{\lambda}}(\mathbf{x},\hat{\mathbf{x}^{n}})}{\partial\lambda_{i}}=2\lambda_{i}d^{t}(x_{i},\hat{x^{n}}_{i})^{2}$
(18)

for a diagonal metric and

$\frac{\partial
d_{\boldsymbol{\Omega}}(\mathbf{x},\hat{\mathbf{x}^{n}})}{\partial\Omega_{ij}}=2d^{t}(x_{j},\hat{x^{n}}_{j})\sum_{d}\Omega_{id}d^{t}(x_{d},\hat{x^{n}}_{d})$
(19)

for a full matrix.

For simplicity, we denote the respective squared distances to the closest
correct and wrong model, respectively, by
$d^{+}=d_{\Theta^{+}}(\mathbf{x}^{n},\hat{\mathbf{x}}^{+})$ and
$d^{-}=d_{\Theta^{-}}(\mathbf{x}^{n},\hat{\mathbf{x}}^{-})$. The term
$\operatorname{sgd}^{\prime}$ is a shorthand notation for
$\operatorname{sgd}^{\prime}((d^{+}-d^{-})/(d^{+}+d^{-}))$. Given a data point
$\mathbf{x}^{n}$ the derivative of the corresponding summand of cost function
$E$ with respect to metric parameters yields

$\frac{\partial{E_{n}}}{\partial\Theta^{+}}=2\operatorname{sgd}^{\prime}\cdot\frac{d^{-}}{(d^{+}+d^{-})^{2}}\cdot\frac{\partial
d^{+}}{\partial\Theta^{+}}$ (20)

for the parameters of the closest correct prototype and

$\frac{\partial{E_{n}}}{\partial\Theta^{-}}=-2\operatorname{sgd}^{\prime}\cdot\frac{d^{+}}{(d^{+}+d^{-})^{2}}\cdot\frac{\partial
d^{-}}{\partial\Theta^{-}}$ (21)

for the parameters attached to the closest wrong model. All other parameters
are not affected. As pointed out before we choose only a global metric such
that the update corresponds to the sum of these two derivatives.

#### Distance measure for functional data

(a) Two functions: Euc = $L^{p}$-norm (b) Two functions: Euc $\neq$
$L^{p}$-norm

Figure 3: Illustration of the $L^{p}$-norm. Plot (a) indicates the case in
which the distance between two functions is equal, both for Euclidean or
$L^{p}$-norm. In plot (b) parts of the functions are interchanging (crossing).
The distance using Euc is still the same as in plot (a) but for the
$L^{p}$-norm the distance is changed, giving a more realistic measure of the
distance of the two functions.

Here we consider a functional distance measure as an extension of the $L^{p}$
norm proposed in ([12]) subsequently denoted as (FUNC). The functional
distance measure has the advantage of taking the functional nature of the data
into account, or in our case the dependence over time, which also constitutes
a function $f(t)$, with potentially discrete arguments $t$. It has been
already successfully used for the analysis of biomedical data as shown in
[16]. The standard Euclidean distance considers the individual features of a
signal independent, so that a change in the order of the positions does not
affect the calculated distance. However, the measurement points over time are
not independent, so that a distance taking this aspect into account can be
considered to be more appropriate for this type of data. Lee proposed a
distance measure taking the functional structure into account by involving the
previous and next values of a signal $v_{i}$ in the $i$-th term of the sum,
instead of $v_{i}$ alone. Assuming a constant sampling period $\tau$, the
proposed norm (FUNC) is:

$\mathcal{L}_{p}^{fc}\left(\mathbf{v}\right)=\left(\sum_{k=1}^{D}\left(A_{k}\left(\mathbf{v}\right)+B_{k}\left(\mathbf{v}\right)\right)^{p}\right)^{\frac{1}{p}}$
(22)

with

$\displaystyle
A_{k}\left(\mathbf{v}\right)=\begin{cases}\frac{\tau}{2}|v_{k}|&\text{if
}0\leq v_{k}v_{k-1}\\\
\frac{\tau}{2}\frac{v_{k}^{2}}{|v_{k}|+|v_{k-1}|}&\text{if
}0>v_{k}v_{k-1}\end{cases}$ (23) $\displaystyle
B_{k}\left(\mathbf{v}\right)=\begin{cases}\frac{\tau}{2}|v_{k}|&\text{if
}0\leq v_{k}v_{k+1}\\\
\frac{\tau}{2}\frac{v_{k}^{2}}{|v_{k}|+|v_{k+1}|}&\text{if
}0>v_{k}v_{k+1}\end{cases}$ (24)

representing the triangles on the left and right sides of $v_{i}$ and $D$
being the data dimensionality. For the data considered in this paper $v$ is a
time series or a prototype reconstruction. As for $L_{p}$, the value of $p$ is
assumed to be a positive integer. At the left and right extremes of the
sequence, $v_{0}$ and $v_{D}$ are assumed to be equal to zero. The concept of
the $L^{p}$-norm is shown in Figure 3. The calculation of this norm is
slightly more complex than that of the standard Euclidean.

###  2.7 Data set description

Subsequently we consider two data sets to evaluate our approach.

####  2.7.1 Simulated data sets

The first one is a simulated two class scenario, proposed in the paper of
[13]. It consists of $100$ samples divided into two classes of $50$ samples
each. For each sample $100$ features have been generated with $8$ time points.
Out of the $100$ features, only $10$ where substantially differentiating
between the classes. The generation mechanism behind the simulated data is to
sample the time series from a piecewise linear function. At a later step,
sample-specific variation is included by shrinking and expanding the curves.

####  2.7.2 Multiple sclerosis data

The second data set is taken from [2] (IBIS) in the prepared form, given in
[6]. The data are taken from a clinical study analyzing the response of
multiple sclerosis (MS) patients to the treatment. Blood sample entrenched
with mono-nuclear cells from $52$ relapsing-remitting MS patients were
obtained $0,3,6,7,12,18$ and $24$ months after initiation of IFN$\beta$
therapy. This resulted on an average of $7$ measurements across the $2$ years.
Expression profiles were obtained using one-step kinetic reverse-transcription
PCR over $70$ genes selected by the specialists to be potentially related to
IFN$\beta$ treatment. Overall, $8\%$ of the measurements were missing due to
patients missing the appointments. After the two year endpoint, patients were
classified as either good or bad responders, depending on strict clinical
criteria. Bad responders were defined as having suffered two or more relapses
or having a confirmed increase of at least one point on the expanded
disability status scale (EDSS). A good responder was to have a suppression of
relapses and not allowed to have an increase on the EDSS. From the $52$
patients, $33$ were classified as good and $19$ as bad responders. A more
detailed description of the data set is available in the paper of [2] and the
supplemented material, therein.

##  3 Results and Discussion

For the simulated and the MS data set, we reanalyzed the classification
accuracy of the SGTM-TT with $9$ hidden states and $4$ basis functions. The
analysis was done within a $4$ fold cross-validation with $5$ repetitions. We
compared it with the general HMM classifier (HMM-Lin) and the discriminative
HMM classifier (HHM-Disc-Lin) proposed in [13]. We also included the results
of [2] who originally proposed the MS study, the analysis of [1], employing a
Kalman Filter combined with an SVM approach and [6] proposing a semi-
supervised analysis coupled with a wrapper and cut-off technique to identify
discriminating features.

Figure 4:  Relevance profile as obtained using SGTM-TT with relevance
learning. The plot shows the average relevance (blue/dark), minimal relevance
(green/bright) and the standard deviation of the relevance, flipped to the
negative part of the relevance axis. We observe that the standard-deviation is
relatively small, hence the relevance profiles of different runs are very
stable. The most discriminative features (high-relevance), can in parts also
be found in [6] but some additional features appear to be relevant, and our
proposed set consists of $7$ genes rather $17$ like in [6]

###  3.1 Simulated data

We applied SGTM-TT with relevance learning for the simulated data set of [13].
We observed an overall prediction accuracy of $94\pm 4$. The relevance profile
identified all known $10$ features and effectively pruned out the remaining
irrelevant data dimensions. Our results are slightly better than those
reported in [13] $(90\%)$ and by [6] $(92\%)$.

###  3.2 Multiple sclerosis experiment

Method | Number of genes | Test accuracy (%)
---|---|---
SGTM-TT | $70$ | $85.66\pm 8.3$
SGTM-TT-R | $7$ | $93.43\pm 5.8$
IBIS | $3$ | $74.20$
Kalman-SVM | - | $87.80$
Lin-Best | $7$ | $85.00$
Costa-Best | $17$ | $92.70\pm 6.1$
Table 1:  Prediction accuracies on the test data for different models using
the MS data set. We observe improved predicition accuracy employing feature
selection. This is also true for SGTM-TT which improved by $\approx 6\%$ using
relevance learning and the SVM classifier. Interestingly also the prediction
accuracy on the full data set, including all features and without relevance
learning is quite good with nearly $84\%$ and hence close to the best result
proposed in [13].

In Table 1 we have summarized the prediction (test-set) results for the
classification of the MS data set in comparison to the results given in [2].
The obtained mappings of the SGTM-TT are topology preserving222In our
observations the topographic error was reasonable small. and we analyzed the
mapping of the points to its prototypes and the neighborhoods. The map for the
first class is depicted for two temporal sequences in Figure 5.

Figure 5:  Illustration of the $3\times 3$ SGTM-TT mapping for the responder
class. Plots in the first row show a typical state sequences. Also if the
state sequences $Z$ are not identical we can expect that the underlying
signals $X$ are similar due to its close neighborhood on the map. This is
reflected by such clustered signals at the bottom. The start of a sequence is
indicated by $\square$ and the termination state by a $\circ$.

As expected, results improved by integration of feature selection or relevance
learning compared to the full feature set. Overall the SGTM-TT with relevance
learning performed very well and achieved good results of $92.5\%$ with
respect to the best reported model and also a smaller number of necessary
features. 333We would like to stress that due to the small sample size and the
$4$ fold cross-validation a missclassification of $1$ point, accounts an error
of $8\%$.. Further the integrated relevance learning avoids multiple, time
consuming runs within a wrapper approach like for the techniques used in [13,
6]. The obtained relevance profile is depicted in Figure 4 and provides direct
access to an interpretation of the relevant features, or marker-candidates,
pruning irrelevant or noise dimensions. The values of the relevance profile
are roughly gaussian distributed with $\mu=0.1$. We calculate a threshold
$\zeta$ for the most relevant features using $\zeta=\mu+\sigma$ and obtain $7$
most relevant features, summarized in Table 2.

Genes | Relevance | found by Lin (7) | found by Costa (17)
---|---|---|---
MAP3K1 | 0.3014 | X | X
NFkBIB | 0.2609 | - | -
IRF8 | 0.2584 | - | X
Caspase 10 | 0.2471 | X | X
Jak2 | 0.1869 | X | X
FLIP | 0.1842 | - | -
RIP | 0.1647 | - | -
Table 2:  Most relevant genes using SGTM-TT with relevance learning.

The SGTM-TT also inherently models different subgroups by the probabilistic
regularizing model of the GTM and GTM-TT [4, 19]. Hence the model complexity
is not so critical provided the map is reasonable large. This is a plus with
respect to the approach presented in [6] which has the number of groups as an
additional meta parameter.

##  4 Conclusion

We have presented a theoretically sound approach for the analysis of short
temporal sequences. It is based on the novel idea to introduce supervision and
relevance learning into Generalized Topographic Mapping through time. Our
results show that we are able to achieve improved or similar performance to
alternative methods for the simulated and the MS data set. Further the
prototype concept of the underlying method permits a better understanding of
the model and extended visualization performance. We also obtain a direct
ranking of the individual features employing the relevance profile, rather by
use of wrapper techniques. In future work we will explore more advance metric
adaptation schemes and alternative functional distance measures. Further we
would like to apply our approach to non-clinical data and make it more
flexible with respect to missing values.

## Acknowledgment

The authors thank: Peter Tino, University of Birmingham for interesting
discussions about probabilistic modeling and support during the early stage of
this project and Falk Altheide, University of Bielefeld and Tien-ho Lin,
Carnegie Mellon University, USA for support with the simulation data. We would
also give extra thanks to Ivan Olier, University of Manchaster, UK; Iain
Strachan, AEA Technology, Harwell, UK and Markus Svensen, Microsoft Research,
Cambridge, UK for providing code and support with the GTM and GTM-TT.

##### Funding:

This work was supported by the German Res. Fund. (DFG), HA2719/4-1 (Relevance
Learning for Temporal Neural Maps) and by the Cluster of Excellence 277
Cognitive Interaction Technology funded in the framework of the German
Excellence Initiative.

## References

  * [1] Russ B. Altman, Tiffany Murray, Teri E. Klein, A. Keith Dunker, and Lawrence Hunter, editors.  Biocomputing 2006, Proceedings of the Pacific Symposium, Maui, Hawaii, USA, 3-7 January 2006. World Scientific, 2006.
  * [2] Sergio E Baranzini, Parvin Mousavi, Jordi Rio, Stacy J Caillier, Althea Stillman, Pablo Villoslada, Matthew M Wyatt, Manuel Comabella, Larry D Greller, Roland Somogyi, Xavier Montalban, and Jorge R Oksenberg.  Transcription-based prediction of response to ifnβ using supervised computational methods.  PLoS Biol, 3(1):e2, 12 2004.
  * [3] Christopher M. Bishop.  Gtm through time.  In In IEE Fifth International Conference on Artificial Neural Networks, pages 111–116, 1997.
  * [4] Christopher M. Bishop, Markus Svensén, and Christopher K. I. Williams.  Gtm: The generative topographic mapping.  Neural Computation, 10(1):215–234, 1998.
  * [5] Karsten M. Borgwardt, S. V. N. Vishwanathan, and Hans-Peter Kriegel.  Class prediction from time series gene expression profiles using dynamical systems kernels.  In Altman et al. [1], pages 547–558.
  * [6] Ivan G. Costa, Alexander Schönhuth, Christoph Hafemeister, and Alexander Schliep.  Constrained mixture estimation for analysis and robust classification of clinical time series.  Bioinformatics, 25(12), 2009.
  * [7] A. Gisbrecht and B. Hammer.  Relevance learning in generative topographic mapping.  Neurocomputing, 74(9):1359–1371, 2011.
  * [8] Christoph Hafemeister, Ivan G. Costa, Alexander Schönhuth, and Alexander Schliep.  Classifying short gene expression time-courses with bayesian estimation of piecewise constant functions.  Bioinformatics, in press, 2011.
  * [9] J. D. Hamilton.  Time Series Analysis.  Princeton University Press, 1994.
  * [10] B. Hammer and Th. Villmann.  Generalized relevance learning vector quantization.  Neural Networks, 15(8-9):1059–1068, 2002.
  * [11] Teuvo Kohonen.  Self-Organizing Maps, volume 30 of Springer Series in Information Sciences.  Springer, Berlin, Heidelberg, 1995.  (2nd Ed. 1997).
  * [12] J. Lee and M. Verleysen.  Generalizations of the lp norm for time series and its application to self-organizing maps.  In Marie Cottrell, editor, 5th Workshop on Self-Organizing Maps, volume 1, pages 733–740, 2005.
  * [13] Tien-ho Lin, Naftali Kaminski, and Ziv Bar-Joseph.  Alignment and classification of time series gene expression in clinical studies.  In ISMB, pages 147–155, 2008.
  * [14] Iván Olier and Alfredo Vellido.  Advances in clustering and visualization of time series using gtm through time.  Neural Networks, 21(7):904–913, 2008.
  * [15] Iván Olier and Alfredo Vellido.  A variational formulation for gtm through time.  In IJCNN, pages 516–521. IEEE, 2008.
  * [16] F.-M. Schleif, T. Riemer, U. Börner, and L. Schnapka-Hille M. Cross.  Genetic algorithm for shift-uncertainty correction in 1-D NMR based metabolite identifications and quantifications.  Bioinformatics, 27(4):524–533, 2011.
  * [17] P. Schneider, M. Biehl, and B. Hammer.  Distance learning in discriminative vector quantization.  Neural Computation, 21:2942–2969, 2009.
  * [18] P. Schneider, K. Bunte, H. Stiekema, B. Hammer, T. Villmann, and M. Biehl.  Regularization in matrix relevance learning.  IEEE Transactions on Neural Networks, 21:831–840, 2010.
  * [19] I. G. D. Strachan.  Latent Variable Methods for Visualization Through Time.  PhD thesis, University of Edinburgh, Edinburgh, UK, 2002.
  * [20] M. Strickert and B. Hammer.  Merge SOM for temporal data.  Neurocomputing, 64:39–72, 2005.
  * [21] Vladimir N. Vapnik.  The nature of statistical learning theory.  Springer New York, Inc., New York, NY, USA, 1995.
  * [22] Thomas Villmann, Ralf Der, Michael Herrmann, and Thomas M Martinetz.  Topology preservation in self-organizing feature maps: exact definition and measurement.  IEEE Transactions on Neural Networks, 8(2):256–266, 1997.
  * [23] Lloyd R. Welch.  Hidden Markov Models and the Baum-Welch Algorithm.  IEEE Information Theory Society Newsletter, 53(4), December 2003\.