# Elaborating Higgs to dimuon decay from gluon fusion
by decorrelation and jet substructure

Subin Han  Hyung Do Kim  Department of Physics and Astronomy and Center for
Theoretical Physics, Seoul National University, Seoul, Korea 08826

###### Abstract

Discovery of the Higgs boson decay to dimuon is anticipated soon based on the
current evidence. Precise categorization of the events without affecting the
invariant mass shape is crucial in the analysis. Decorrelation of the
invariant mass and the output of discriminators (the score of discriminators)
is essential for consistent and precise analysis. In this paper we use
distance correlation as the additional loss function to achieve the
decorrelation for discriminators and examine various analysis methods. The
analyses with and without jet substructure variables are presented. Adding jet
substructure variables considerably improves the significance of the Higgs to
dimuon signal from gluon fusion.

##  I Introduction

After Higgs discovery Chatrchyan  _et al._ (2012)Aad  _et al._ (2012), various
properties of Higgs have been measured including Higgs couplings to Standard
Model (SM) particles. The SM predicts the Higgs coupling to a particle is
proportional to its mass. Accordingly, Higgs couplings to heavy particles,
such as gauge bosons and third generation fermions, are well measured Aad  _et
al._ (2020)Aaboud  _et al._ (2018)Sirunyan  _et al._ (2018)Aaboud  _et al._
(2019). However, couplings to relatively light particles have not been
satisfactorily measured so far. The Higgs to dimuon channel is the most
sensitive channel among the second generation fermions at the LHC due to the
relatively clean leptonic final states. The property of Higgs couplings to
light fermions can verify the Higgs mechanism of the SM more accurately.
Furthermore, the growing interest towards a muon collider makes this channel
more interesting.

In 2020, the decay of Higgs to dimuon was barely confirmed to be consistent
with the SM Sirunyan  _et al._ (2021). The observed (expected) significance
was 3.0$\sigma$ (2.5$\sigma$). The analysis was performed separately for each
Higgs production channel. However, the signal-to-background ratio for each
channel was not sufficient to extract the number of signal events ($\sim
0.3\%$ for ggH channel at Higgs mass). To increase the sensitivity, each
channel was further divided into several categories with different signal-to-
background ratios based on the output of discriminators. The dimuon invariant
mass $m_{\mu\mu}$ is an excellent variable for distinguishing the resonance
peak of the Higgs signal from the background events. Therefore extraction of
the number of signal and background events as well as the calculation of
discovery significance were conducted by fitting data to the expected
$m_{\mu\mu}$ distribution. Several uncertainties arising from the simulation
could be reduced through shape analysis during the distribution fitting.

A crucial requirement for the analyses using categorization above is that the
information about the invariant mass should not be used in the categorization
process. In other words, the output of discriminators, or the discriminator
score, should be statistically independent of the invariant mass.
Categorization by a score which is statistically dependent on the invariant
mass generally results in different invariant mass distributions for each
category, $f(m|\text{score})\neq f(m)$, where $f$ represents the probability
density function. If the invariant mass distribution is distorted, the number
of events extracted from each category can not be reliably trusted. In extreme
cases, a resonance peak near the Higgs mass may appear in the background
$m_{\mu\mu}$ distribution (See Fig.1 for example) adding a significant
uncertainties on the number of signal events near the Higgs bump region. In
addition, if we fit to the invariant mass after the categorization by a
discriminator correlated to the invariant mass, it results in double usage of
the invariant mass, compounding the error. Therefore the standard likelihood
fitting to the invariant mass distribution can be consistently combined with
machine discriminator results only when two results are
decorrelated/deassociated. 111Technically, the term
‘deassociate/deassociation’ is more appropriate Altman and Krzywinski (2015),
but we will use the term ‘decorrelate/decorrelation’ following other
literature Kasieczka and Shih (2020)

To avoid the distortions in the invariant mass distribution that arise when
events are selected by a discriminator correlated with the invariant mass,
researchers have not only securely excluded the invariant mass but also
carefully selected other variables to make the categories as insensitive to
the invariant mass as possible. However, as will be seen in the paper, machine
discriminators generally learn about $m_{\mu\mu}$ since it has strong
discrimination power. An ad hoc choice of variables which seemingly unrelated
to the $m_{\mu\mu}$, qualitative check of distortion in distributions after
categorization Bourilkov  _et al._ (2019) and a posteriori modification of
fitting functions are not satisfactory solutions. In this paper we provide a
quantitative and consistent framework to deal with this problem. Among the
various measures of dependence between two sets of feature variables, we
choose ‘distance correlation’ Kasieczka and Shih (2020) Székely  _et al._
(2007) since it is differentiable, so being compatible with deep learning
frameworks, well defined in both population and sample spaces, and capable of
capturing both linear and non-linear dependence. At the same time the range is
well defined such that it is 0 when there is no correlation and is 1 when it
is fully deterministic/dependent. By leveraging the flexiblility of deep
neural network (DNN), we use the distance correlation between $m_{\mu\mu}$ and
output of discriminator as an additional loss for training the discriminator.
As a result, it has been found that almost perfect statistical independence
(for both signal and background data classes simultaneously) between
$m_{\mu\mu}$ and the discriminator output can be systematically achieved,
enabling us to choose the required degrees of independence, which is crucial
for accurate analysis, despite being in a trade-off relation to the
performance of the discriminator. By carefully balancing this trade-off, this
approach ensures reliable categorization without compromising the integrity of
the invariant mass distribution.

Compared to the analysis using only the kinematic variables of reconstructed
particle objects, incorporating jeet substructure variables of the jets
associated with the dimuon system as additional input features for the
discriminator can provide extra power, especially in discriminating the gluon
dominant ggH signal from other background processes. Currently the vector
boson fusion (VBF) is the leading production channel for the Higgs to dimuon
discovery despite the fact that the ggH production cross section is more than
10 times larger than VBF one at $\sqrt{s}=13\ $TeV CERN (2017). In Cho  _et
al._ (2020), it has been demonstrated that ggH channel can be exploited to
improve the sensitivity of the Higgs invisible signal if jet substructure
variables are fully utilized to distinguish the flavors and kinematics of ISR
by DNN analysis. The difference in parton luminosities and QCD color factors
of vertices results in the signal ggH ISR to be gluon-rich, while the
background Drell-Yan (DY) ISR to be quark-rich. In addition, due to the
helicity conservation, ggH ISR has a deficit of quark jets at the central
rapidity region. 222There is no issue of decorrelation in Higgs to invisible
study, the treatment of systematic uncertainties originating from the
simulation remain as open question. In this study, the same method is used to
improve the sensitivity of Higgs to dimuon signal from ggH channel.

The paper is composed as following. In section.II, we first outline the event
simulation and ggH selection cuts. We also summarize the analysis strategy
including training variables, common training schemes, and the categorization
procedure. In section.III, we introduce three types of discriminators and
compare their result in detail. In particular, with the $m_{\mu\mu}$
decorrelated discriminator implemented using the distance correlation, the
degree of statistical independence and the discrimination performance will be
compared. In section. IV, we compute the discovery significance using the
categorized events and check the effect of including jet substructure
variables.

##  II Events and analysis strategy

When an event is classified being from ggH, it means that the event passed all
selection cuts to identify ggH, regardless of genuine origin of the event.
Therefore, all possible Higgs production channels have to be simulated.
However, the contribution from other channels to the ggH signals is much
smaller than genuine ggH channel due to their small cross sections. Next to
leading order contribution is from VBF channel which is less than 5% of ggH
channel. Thus we simulated only ggH and VBF channel. Similarly, there are many
channels which act as backgrounds against ggH. Among them, we included the
dominant background DY and t$\bar{\text{t}}$ which has non-negligible
contribution to multi-jet process (more than 20% of DY). All the events are
simulated by MadGraph5_aMC@NLO Alwall  _et al._ (2014) at $\sqrt{s}=13\ $TeV.
Each simulation was done inclusively for the number of jets, so that we can
properly include both soft and hard jet emissions.

After the simulation, the selection of events for ggH was conducted.
Basically, an event is required to have two oppositely charged muons with
$p_{T}>20\ $GeV and $|\eta|<2.4$. Any events including b-jet and/or additional
charged lepton are vetoed to avoid overlap with t$\bar{t}$H and VH channel.
Finally, events with two energetic and largely separated jets ($m_{jj}>400\
$GeV, $|\Delta\eta_{jj}|>2.5)$ are rejected as those are to be identified as
VBF.

To enhance the signal-to-background ratio, various information other than
$m_{\mu\mu}$ is used to divide the total data into some categories. Recently,
with huge development of machine learning techniques, output of machine
discriminators are frequently used as a criterion for categorization to use
various variables simultaneously. The set of training variables used by the
CMS collaboration is denoted by $S_{n_{\text{jet}}}$ for each number of jet
$n_{\text{jet}}$. $S_{n_{\text{jet}}}$ are carefully selected to exclude the
$m_{\mu\mu}$ information. Please note that $n_{\text{jet}}=2$ means events
with at least two jets and the kinematics of jets are included in terms of
pseudo rapidity $\eta$ of leading and next to leading jet. To exploit ISR
characteristic of ggH channel, we also included set of jet substructure
variables $JS\equiv\\{n_{\text{track}},\text{girth},$
$\text{broadening},C_{1}^{\beta=0.2},\text{RMS}-p_{T},\text{Pull-vector}\\}$
for each jet. Therefore total set of variables including jet substructure
variables are $J_{0}\equiv S_{0}$, $J_{1}\equiv S_{1}\cup JS_{1}$ and
$J_{2}\equiv S_{2}\cup JS_{1}\cup JS_{2}$. See Appendix A for more detail.

The discriminators are trained separately for events with different number of
jets $n_{\text{jet}}$ because we have different number of training variables
for each $n_{\text{jet}}$. The training sample consists of one million signals
and backgrounds respectively, 70% for training and 30% for validation. We will
train three types of machine discriminators to compare the results; BDT as a
benchmark, (ordinary) DNN and decorrelated DNN. After the training, total data
set is divided into categories according to the output of discriminators, or
equivalently, by the signal efficiency. The boundaries are recursively
selected to maximize the total expected significance with $S/\sqrt{B}$ as a
figure of merit, while minimizing the number of categories (until the gain
from additional boundary is less than 1%).

A typical next step is to extract the number of events from fitting of data to
$m_{\mu\mu}$ distribution. However we will skip the extraction procedure and
use the number of events by simulation directly because we do not have real
observed data to fit. Since we already used large data, fitting to another set
of simulated data is not worth much. Though our result depends on the
simulated data, there is no significant problem in comparison between
discrimination schemes and check the impact of including jet substructure
variables. Therefore after the categorization, we will compute the discovery
significance by profile likelihood fitting to the $m_{\mu\mu}$ distribution.

##  III  Discriminators

The astonishing success of DNN is based on the ability to build a hidden
representation from the training variables. However, due to the very feature,
DNN generally reconstruct the invariant mass even if it is not directly
included in training variables. As a result, the distributions of categorized
events are largely distorted. We used PyTorch Paszke  _et al._ (2019) to
implement DNN and used SCIKIT LEARN Pedregosa  _et al._ (2011) for
preprocessing. Selected models are given in Appendix.B.2. Fig.1 shows the
actual distortion by ordinary DNN discriminator. In categories with high
signal purity, even the background distribution resembles the resonance peak
of the signal. It can mimic signals and overestimate the number of signal
events. Therefore, the extraction of the number of events by fitting to
$m_{\mu\mu}$ is not valid any more under such distortion. In such cases, it
would be more appropriate to extract the number of events by fitting signal
strength modifier to the distribution of output directly. However, this method
suffers from large systematic uncertainties from simulation and training and
therefore fittint to the invariant mass is more preferred generally.

(a) DNN; without jet substructure variables($S_{n_{\text{jet}}}$)

(b) DNN; with jet substructure variables($J_{n_{\text{jet}}}$)

Figure 1: $m_{\mu\mu}$ distribution binned by optimal boundaries chosen by
output of (ordinary) DNN discriminators.

The CMS collaboration used boosted decision tree (BDT) to categorize their ggH
data. The tree-based models including BDT are essentially complex collections
of cuts on training variables. Therefore, generally BDT has good performance
for theoretically well-motivated training variables such as $p_{T},\eta$ and
so on. It is also easy to interpret the result. However, compared to network-
based models, it lacks the ability of capturing complex non-linear
relationships among training variables. Generally it is the reason why network
based discriminators such DNN are preferred for more complex data. Ironically,
this intrinsic weakness of tree-based model would be profitable for making the
resultant discriminator to be independent to $m_{\mu\mu}$.

To compare our result, we also conducted the BDT discrimination. We used the
gradient-boosted tree implemented in SCIKIT LEARN. The Friedman MSE is used as
a splitting criterion and a deviance is used as a loss. Selected models are
given in Appendix.B.1. Fig.2 shows $m_{\mu\mu}$ distribution categorized by
output of BDT discriminators. It is obvious that the $m_{\mu\mu}$
distributions are much more similar compared to the ordinary DNN case(Fig.1).
However, the distributions are not exactly the same. There is a difference in
slope of background distributions which also can be problematic for extracting
the number of signal events. Actually the CMS collaboration was aware of it
and use mass side band events to constrain the background distribution with
modified fitting function having common shape and category specific slope (The
‘core-pdf’ method).

(a) BDT; without jet substructure variables($S_{n_{\text{jet}}}$)

(b) BDT; with jet substructure variables($S_{0},J_{1}$ and $J_{2}$)

Figure 2: $m_{\mu\mu}$ distribution binned by optimal boundaries chosen by
output of BDT discriminators.

For both BDT and DNN, the discrimination result is not statistically
independent to $m_{\mu\mu}$. To achieve the independence, we exploit the
flexibility of DNN. The network in DNN is an analytic function and it can
handle wide variety of loss functions, which only requires to be
differentiable. Therefore we included ‘decorrelation loss’
$\mathcal{L}_{\text{decor}}$ with a parameter $\lambda$ as
$\mathcal{L}_{\text{tot}}\equiv\mathcal{L}_{\text{discrim}}+\lambda\mathcal{L}_{\text{decor}}$,
where the discrimination loss $\mathcal{L}_{\text{discrim}}$ is the binary
cross entropy. The $\mathcal{L}_{\text{discrim}}$ and
$\mathcal{L}_{\text{decor}}$ are in conflict since $m_{\mu\mu}$ has strong
discrimination power. However what we want to achieve is the decorrelation
even if it costs discrimination performance. It can be achieved by setting
$\lambda$ to be large enough. The distance correlation ($\equiv\mathcal{R}$)
is used as a measure of statistical dependence. We required the discrimination
score to be independent to $m_{\mu\mu}$ separately for signal and background
events,
$\mathcal{L}_{\text{decor}}\equiv\mathcal{R}(m_{\mu\mu},\text{score}|\text{signal})+\mathcal{R}(m_{\mu\mu},\text{score}|\text{background})$.
We will call this DNN model with decorrelation loss as ‘decorrleated DNN’.

The decorrelated DNN is implemented with PyTorch again. We added decorrelation
loss term every two epochs for calculation efficiency. Furthermore, since the
distance matrix requires large memory, 1% of every batch is sampled to
calculate $\mathcal{R}$. To test the degree of decorrelation at the end of
training, the distance correlation is calculated with sample of size 25,000
for signal and background respectively (which is actually the limitation of
our gpu memory). Selected models are given in Appendix.B.3. Fig.3 shows the
$m_{\mu\mu}$ distributions which are not distorted by the categorization.
Therefore, there will be no systematic error from dependence.

(a) decorrelated DNN; without jet substructure variables
($S_{n_{\text{jet}}}$)

(b) decorrelated DNN; with jet substructure variables ($J_{n_{\text{jet}}}$)

Figure 3: $m_{\mu\mu}$ distribution binned by optimal boundaries chosen by
output of decorrelated DNN discriminators.

###  III.1 Comparison

As mentioned before, the decorrelation is in conflict with discrimination.
Therefore the comparison among three machine discrimminators - BDT,
(ordinary)DNN and decorrelated DNN - can be summarized in terms of two
properties, the performance of discrimination and the degree of decorrelation.

Figure 4: Receiver operating characteristic(ROC) curve for three
discrimination schemes - BDT, (ordinary) DNN and decorrelated DNN - with two
different set of variables($S_{n_{\text{jet}}}$ and $JS_{n_{\text{jet}}})$.
The area under curve(AUC) is also given. Note that the true positive rate is
equivalent to signal efficiency and the false positive rate is equivalent to
(1-background rejection rate).

Fig. 4 shows the receiver operating characteristic (ROC) curve with area under
curve (AUC). It shows the performance of each discriminator. For all three
discriminators, including jet substructure variables (solid lines) improve the
discrimination performance of ggH channel. The ordinary DNN’s (black) have
large AUC since they exploited $m_{\mu\mu}$ internally. BDT (blue) and
decorrelated DNN (red) have similar AUCs, which is around 0.1 less than those
of ordinary DNN. The gap between ordinary DNN’s and other methods represents
the discrimination power of $m_{\mu\mu}$, or more conservatively, power of the
information of training variables correlated to $m_{\mu\mu}$. The small AUC’s
of decorrelated DNN reflects that the decorrelation intentionally discards the
information about $m_{\mu\mu}$ and those of BDT reflects that the BDT cannot
reconstruct $m_{\mu\mu}$ internally. The small difference of AUC between BDT
and decorrelated can be interpreted as the the uncontrolled contamination by
$m_{\mu\mu}$.

The degree of decorrelation between $m_{\mu\mu}$ and the output of
discriminators can be seen directly from the distortion of invariant mass
distribution categorized by the output of discriminators (Fig.1, 2, 3). It can
also be quantitatively checked by a test of independence using asymptotic
distribution of distance covariance $\mathcal{V}$ Székely  _et al._ (2007).
The test statistic is $V\equiv\frac{n\times\mathcal{V}^{2}(X,Y)}{S_{2}}$,
where $n$ is the size of smaple and
$S_{2}\equiv\frac{1}{n^{2}}\sum_{k,l}|X_{k}-X_{l}|\frac{1}{n^{2}}\sum_{k,l}|Y_{k}-Y_{l}|$.
The statistical independence would be rejected if
$\sqrt{V}>\Phi^{-1}(1-\alpha/2)$ for selected significance level $\alpha$,
where $\Phi$ is the Gaussian cumulative density function.

Fig.5 shows the statistic $\sqrt{V}$ between $m_{\mu\mu}$ and the output of
discriminator. Samples of size 10000 with 1000 iterations are used.

Figure 5: The statistic $\sqrt{V}$ of signal and background separately for
various discriminators. The error bar comes from 1,000 iterations and sample
size 10,000 was used. Please note that the scale of y-axis is arbitrary.

It is interesting that the result becomes more independent for all the cases
when the jet substructure variables are included for training. This is because
jet substructure variables are useful for discrimination of ggH signal from
background but relatively less correlated to $m_{\mu\mu}$ compared to
$S_{n_{\text{jet}}}$ variables including information about muon and dimuon
system. It is obvious that $m_{\mu\mu}$ and the output of discriminators are
enormously correlated for (ordinary) DNN. The independence hypothesis is
always rejected for DNN case even under $\alpha<0.05$. The situation is
improved in the case of BDT but the uncontrolled impact of $m_{\mu\mu}$ still
exist. The BDT trained with jet substructure variables barely pass the test up
to $\alpha=0.05$. On the other hand, the output of decorrelated DNN and
$m_{\mu\mu}$ is independent up to $\alpha=0.215$, which is the upper limit of
significance level given in Székely  _et al._ (2007).

###  III.2 Decorrelated DNN result in detail

After imposing distance correlation loss, the decorrelated DNN loses its
accuracy but becomes independent of $m_{\mu\mu}$. Therefore the $m_{\mu\mu}$
can be used without systematic error after categorization. Now, have a
detailed look at the categorization result of decorrelated DNN.

Fig.6 shows the number of events for each category and production channel. The
number of signal and background is stacked respectively. The vertical dotted
lines are the optimal boundaries.

(a) without jet substructure variables ($S_{n_{\text{jet}}}$)

(b) with jet substructure variables($J_{n_{\text{jet}}}$)

Figure 6: Number of events according to optimized bins according to signal
efficiency. (a)/(b): For a model trained without/with jet substructure
variables.

Table.1, 2 shows the number of signal, portion of each production channel,
number of background, $S/(S+B)$ and $S/\sqrt{B}$ for optimized categories
without and with jet substructure variables. In Fig.6a and Table.1, the
portion of VBF among signal is larger for higher purity category. This
tendency stems from the discrimination power of two characteristic jets of
VBF. On the other hand, the tendency of growing VBF portion is totally
reversed for the discriminators trained with jet substructure variables. In
Fig.6b, Table.2, the portion of VBF events are getting smaller for higher
signal purity category. This reversal of VBF portion represents the improved
discrimination power of ggH signals by exploiting jet substructure variables.
This property can be used to reduce errors from analysis done separately for
each production channel and improve the precision of ggH channel.

Signal | Sig | ggH | VBF | Bkg | $\frac{S}{S+B}$ | $\frac{S}{\sqrt{B}}$
---|---|---|---|---|---|---
efficiency |  | (%) | (%) |  | (%) |
0-25% | 212.45 | 97.5 | 2.5 | 545713 | 0.04 | 0.29
25-45% | 169.93 | 96.6 | 3.4 | 193690 | 0.09 | 0.39
45-75% | 170.17 | 95.6 | 4.4 | 130068 | 0.13 | 0.47
75-90% | 212.65 | 94.2 | 5.8 | 96033 | 0.22 | 0.69
90-95% | 42.65 | 92.5 | 7.5 | 11317 | 0.38 | 0.40
95-100% | 42.41 | 91.3 | 8.6 | 7254 | 0.58 | 0.50
Table 1: Categorization for decorrelated DNN without jet substructure variables. Signal | Sig | ggH | VBF | Bkg | $\frac{S}{S+B}$ | $\frac{S}{\sqrt{B}}$
---|---|---|---|---|---|---
efficiency |  | (%) | (%) |  | (%) |
0-10% | 85.09 | 92.6 | 7.4 | 389125 | 0.02 | 0.14
10-25% | 127.42 | 94.0 | 6.0 | 216714 | 0.06 | 0.27
25-60% | 297.62 | 94.7 | 5.3 | 260441 | 0.11 | 0.58
60-75% | 127.66 | 96.7 | 3.3 | 61487 | 0.21 | 0.51
75-90% | 127.53 | 97.9 | 2.1 | 40971 | 0.31 | 0.63
90-95% | 42.52 | 98.3 | 1.7 | 9348 | 0.45 | 0.44
95-100% | 42.41 | 98.6 | 1.4 | 5933 | 0.71 | 0.55
Table 2: Categorization for decorrelated DNN with jet substructure variables.

##  IV Discovery significance

The statistical significance is calculated by the profile likelihood ratio fit
to the $m_{\mu\mu}$. We used uniform binning of size $50$ MeV. To incorporate
systematic errors, nuisance parameters $\theta$ are included with Gaussian
suppression. Since the simulation was done for arbitrary luminosity, the
number of events are normalized to that of Sirunyan  _et al._ (2021), $137\
$fb-1. We used non-integers from the normalization directly to Poisson counts
to avoid the subtle risk of getting the over- or under-estimated results from
the rounding off. The problem of non-integers for Poisson distribution is not
a genuine problem since it will be canceled off by calculating likelihood
ratio Cowan  _et al._ (2011).

To compare the result, significance is computed for BDT and decorrelated DNN.
We do not compute for ordianry DNN case since it is totally distorted. Please
note that since we do not extract the number of events by fitting and use the
simulation data directly, the significance computed for BDT do not suffer from
the problem in extraction though there is distortion indeed. The true
significance will be different but cannot predict the tendency since it
largely depend on the prescriptions such as the ‘core-pdf’ method used in
Sirunyan  _et al._ (2021). Table.3 shows the AUC and significance. Clearly,
significance increases about 0.15 by exploiting ggH ISR characteristic. The
slightly difference of AUC and significance between BDT and decorrelated dNN
can be interpreted as uncontrolled impact of $m_{\mu\mu}$.

Machine | BDT | decorrelated DNN
---|---|---
Variables | $S_{n_{\text{jet}}}$ | $J_{n_{\text{jet}}}$ | $S_{n_{\text{jet}}}$ | $J_{n_{\text{jet}}}$
significance | (1.64) | (1.78) | 1.61 | 1.76
AUC | 0.72 | 0.77 | 0.71 | 0.76
Table 3: Discovery significance table. CMS benchmark was 1.56$\sigma$ for ggH
and 1.77 $\sigma$ for VBF. The numbers in the parenthesis are obtained without
sufficient decorrelation and should be interpreted cautiously.

##  V Summary

We applied decorrelation using distance correlation to the Higgs to dimuon
discovery from ggH channel. The decorrelation between output of discriminator
and $m_{\mu\mu}$ allows us to categorize the events without any systematic
error from $m_{\mu\mu}$ correlations. The decorrelated DNN achieved the
independence simultaneously for both signal and background events up to
significance level $\alpha=0.215$. The jet substructure variables are helpful
for improving discovery significance of ggH channel and the significance was
improved from $1.61$ to $1.76$.

Decorrelation can be helpful in general for both precision and discovery
researches as it can preserve the variable of interest-distribution of events
after the cut and get rid of systematic errors from distortion in the process
of categorization. In addition, all the variables can be used for
categorization without any concern of correlation.

##  VI Acknowledgements

We thank Wonsang Cho for collaboration at an initial stage of the work and for
discussion and comments. This work is supported by NRF of Korea grant, No.
0426-20230001.

##  Appendix A Training variables

The set of training variables for events with each number of jets
$n_{\text{jet}}$ will be summarized. We will denote the variables used by the
CMS collaboration without jet substructure variables as $S_{n_{\text{jet}}}$
and variables with jet substructure variables as $J_{n_{\text{jet}}}$. Please
note that $n_{\text{jet}}=2$ means events with at least two jets. Therefore
$S_{2}$ and $J_{2}$ mean training variables for events with multi-jets, not
only events with exactly two jets.

###  A.1 CMS used variables

The training variables were carefully selected to exclude $m_{\mu\mu}$
information. The set of training variables $S_{n_{\text{jet}}}$ used by CMS
collaboration are

$\displaystyle S_{0}\equiv$
$\displaystyle~{}\\{p_{T}^{\mu\mu},y_{\mu\mu},\phi_{CS},\cos\theta_{CS},\eta(\mu_{1}),\eta(\mu_{2})$
$\displaystyle~{},\frac{p_{T}^{\mu 1}}{m_{\mu\mu}},\frac{p_{T}^{\mu
2}}{m_{\mu\mu}}\\}$ $\displaystyle S_{1}\equiv$
$\displaystyle~{}S_{0}\cup\\{p_{T}(j_{1}),\eta(j_{1}),\Delta
R(\mu\mu,j_{1})\\}$ $\displaystyle S_{2}\equiv$
$\displaystyle~{}S_{1}\cup\\{p_{T}(j_{2}),\eta(j_{2}),m_{j_{1}j_{2}},\Delta\eta_{j_{1}j_{2}},\Delta\phi_{j_{1}j_{2}}$
$\displaystyle~{},\Delta\eta_{\mu\mu,j_{1}},\Delta\eta_{\mu\mu,j_{2}},\Delta\phi_{\mu\mu,j_{1}},\Delta\phi_{\mu\mu,j_{2}}$
$\displaystyle~{},z_{*},n\\}$

In the above notation, CS means those computed in the dimuon Collins-Soper
rest frame, $j_{1,2}$ respectively denotes the highest and second highest
$p_{T}$ jet,
$z*\equiv\frac{y_{\mu\mu}-(y_{j_{1}}+y_{j_{2}})/2}{|y_{j_{1}}-y_{j_{2}}|}$ is
the Zeppendfeld variable and $n$ is the number of jets with $p_{T}>25$GeV,
$|\eta|<4.7$.

###  A.2 Jet substructure variables

To exploit the ISR characteristic of ggH channel, we used jet substructure
variables in addition to the conventional variables. The jet substructure
variables are generally defined as weighted moment sum of constituent
transverse momenta normalized by transverse momentum of the jet to capture the
fatness of gluon jet. For each jet, we used the number of charged tracks in
the jet $n_{\text{track}}$Gallicchio and Schwartz (2011), girth (linear radial
moment) Gallicchio and Schwartz (2011)Gallicchio  _et al._ (2011),
broadeningCatani  _et al._ (1992), Energy-energy correlation(EEC)
$C^{\beta=0.2}_{1}$ with $\beta=0.2$ Larkoski  _et al._ (2013)Dasgupta  _et
al._ (2013), RMS-$p_{T}$Gallicchio and Schwartz (2011), Pull-vectorGallicchio
and Schwartz (2010).

girth $\displaystyle
G=\frac{1}{p^{\text{jet}}_{T}}\displaystyle\sum_{i\in\textrm{jet}}p^{i}_{T}|\Delta\vec{r}_{i}|$
broadening $\displaystyle B=\frac{1}{\sum_{i}|\vec{p\
}^{i}|}\displaystyle\sum|\vec{p\ }^{i}\times\hat{p}^{\textrm{jet}}|$
$\displaystyle~{}~{}=\frac{1}{\sum_{i}|\vec{p\
}^{i}|}\displaystyle\sum|\vec{k\ }^{i}_{T}|$ EEC $\displaystyle
C^{\beta}_{1}=\frac{1}{(\sum_{i}p^{i}_{T})^{2}}\displaystyle\sum_{i<j}p^{i}_{T}p^{j}_{T}(\Delta
R_{ij})^{\beta}$ RMS-$p_{T}$
$\displaystyle\sqrt{<p^{2}_{T}>}=\frac{1}{p^{\textrm{jet}}_{T}}\sqrt{\frac{1}{n_{\textrm{tk}}}\displaystyle\sum_{i}(p^{i}_{T})^{2}}$
Pull-vector
$\displaystyle\vec{v}_{p}=\frac{1}{p^{\textrm{jet}}_{T}}\displaystyle\sum_{i}p^{i}_{T}|\Delta\vec{r\
}^{i}|\Delta\vec{r\ }^{i}$

##  Appendix B Discrimination models in detail

###  B.1 BDT

A minimum sample size for a node to be split (min_samples_split) is 1% to
regularize the over-fitting. We used stochastic gradient boosting with $70\%$
of sample for each tree. Table 4 shows the resultant (hyper-)parameters and
accuracy for each machine. The minimum sample size after the splitting is once
more limited by min_samples_leaf (MSL).

| w/o jet sub. | w/ jet sub.
---|---|---
Variables | $S_{0}$ | $S_{1}$ | $S_{2}$ | $J_{1}$ | $J_{2}$
$n_{\text{tree}}$ | 80 | 1500 | 1000 | 1500 | 1000
max depth | 3 | 5 | 7 | 7 | 9
learning rate | 0.4 | 0.2 | 0.2 | 0.1 | 0.1
MSL | 0.001 | 0.01 | 0.02 | 0.001 | 0.001
accuracy | 0.66 | 0.65 | 0.66 | 0.72 | 0.73
Table 4: (hyper-)Parameters and accuracy of BDT. MSL stands for
min_samples_leaf.

###  B.2 (Ordinary) DNN

Table 5 shows the resultant (hyper-)parameters and accuracy for each ordinary
DNN machine. Through all the machines, 256 nodes per layer and 100 batches are
commonly used.

| w/o jet sub. | w/ jet sub.
---|---|---
Variables | $S_{0}$ | $S_{1}$ | $S_{2}$ | $J_{1}$ | $J_{2}$
$n_{\text{layers}}$ | 2 | 3
learning rate | 0.007 | 0.01 | 0.003 | 0.04 | 0.005
epochs | 100 | 300 | 400 | 300 | 300
accuracy | 0.76 | 0.75 | 0.75 | 0.79 | 0.79
$\mathcal{R}$(sig) | 0.35 | 0.39 | 0.33 | 0.32 | 0.24
$\mathcal{R}$(bkg) | 0.31 | 0.31 | 0.30 | 0.28 | 0.18
Table 5: (hyper-)Parameters and accuracy of DNN. The resultant distance
correlation $\mathcal{R}$ the distance correlation $\mathcal{R}$ is calculated
with 25,000 events each from the test sample.

###  B.3 Decorrelated DNN

Table 6 shows the resultant (hyper-)parameters and accuracy for each
decorrelated DNN machine. Through all the machines, 128 nodes per layers and
100 batches are commonly used. The relative multiplier $\lambda=8\times
10^{-3}$ is also common. We made the $\lambda$ to be magnified to 100 times of
its input value if $\mathcal{L}_{\text{decor}}>10^{-3}$ during the training.
Therefore the $\lambda=8\times 10^{-1}$ for almost all the time during the
training. Since the $\mathcal{L}_{\text{decor}}$ is the sum of distance
correlations of signal and background and it is computed with smaller sample,
it is larger than the value of distance correlation computed for the test
sample(**).

| w/o jet sub. | w/ jet sub.
---|---|---
Variables | $S_{0}$ | $S_{1}$ | $S_{2}$ | $J_{1}$ | $J_{2}$
$n_{\text{layers}}$ | 3 | 3 | 6 | 3 | 6
learning rate | 0.0004 | 0.001 | 0.01 | 0.001 | 0.001
epochs | 30 | 40 | 60 | 60 | 60
$\mathcal{L}_{\text{decor}}$* | 0.046 | 0.047 | 0.044 | 0.037 | 0.038
accuracy | 0.65 | 0.64 | 0.65 | 0.71 | 0.72
$\mathcal{R}$(sig)** | 0.006 | 0.009 | 0.011 | 0.008 | 0.008
$\mathcal{R}$(bkg)** | 0.019 | 0.018 | 0.022 | 0.010 | 0.010
Table 6: Resultant model (hyper-)parameters and accuracy. *:
$\mathcal{L}_{\text{decor}}$ is the sum of distance correlation loss between
score and $m_{\mu\mu}$ respectively for signal and background. 1% of every
batch is ampled to calculate $\mathcal{L}_{\text{decor}}$. **: Here the
distance correlation $\mathcal{R}$ is calculated with 25,000 events each from
the test sample.

##  Appendix C Independence test

We checked the performance of test using $V$ with the well-known distribution,
the two-dimensional Gaussian with various correlation. Fig.7 shows the
fraction of rejected samples among 500 iterations(equivalent to type 1 error
rate for zero correlation case and statistical power for non-zero correlation
cases) with different sample size and correlation. More than sample size 5000,
the samples with correlation $\geq$ 0.1 are all rejected for significance
level $\alpha\geq$ 0.1.333In Table 1 of Székely  _et al._ (2007), they
reported type 1 error rate about 0.1 for various independent distributions,
sample sizes under $\alpha=0.1$. In Fig. 2, the power of test for Gaussian
with correlation $\rho=0.1$ and under $\alpha=0.1$ is also given.

(a) $\alpha=0.1$

(b) $\alpha=0.215$

(c) $\alpha=0.215$

Figure 7: Fraction of rejected samples among 500 iterations with two-
dimensional Gaussian, for various sample size, correlation and significance
level $\alpha$

Thus we conducted the test of independence with size 5000 and significance
level $\alpha=0.215$ for 500 iterations. The rejected fraction is given in
Table.7

Machine | BDT | DNN | decor.DNN
---|---|---|---
Variable | $S_{n_{\text{jet}}}$ | $J_{n_{\text{jet}}}$ | $S_{n_{\text{jet}}}$ | $J_{n_{\text{jet}}}$ | $S_{n_{\text{jet}}}$ | $J_{n_{\text{jet}}}$
sig(%) | 100 | 76.6 | 100 | 100 | 30.2 | 17.4
bkg(%) | 99.8 | 64.6 | 100 | 100 | 12.0 | 16
Table 7: Rejection rate among 500 iterations for the test of dependence.
Significance level $\alpha=0.215$ and size = 5000 was used.

## References

  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS Collaboration), Physics Letters B 716, 30 (2012).
  * Aad  _et al._ (2012) G. Aad _et al._ (ATLAS Collaboration), Physics Letters B 716, 1 (2012).
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS Collaboration), Phys. Rev. D 101, 012002 (2020).
  * Aaboud  _et al._ (2018) Aaboud _et al._ (ATALS Collaboration), Physics Letters B 784, 173 (2018).
  * Sirunyan  _et al._ (2018) Sirunyan _et al._ (CMS Collaboration), Physical review letters 121, 121801 (2018).
  * Aaboud  _et al._ (2019) Aaboud _et al._ (ATLAS Collaboration), Physical Review D 99, 072001 (2019).
  * Sirunyan  _et al._ (2021) A. M. Sirunyan _et al._ (CMS Collaboration), JHEP 01, 148 (2021), arXiv:2009.04363 [hep-ex] .
  * Altman and Krzywinski (2015) N. Altman and M. Krzywinski, Nature Methods 12, 899 (2015).
  * Kasieczka and Shih (2020) G. Kasieczka and D. Shih, Physical Review Letters 125 (2020), 10.1103/physrevlett.125.122001.
  * Bourilkov  _et al._ (2019) D. Bourilkov, D. Acosta, P. Bortignon, A. Brinkerhoff, A. Carnes, S. Gleyzer, and B. Regnery (CMS Collaboration), EPJ Web Conf. 214, 06002 (2019).
  * Székely  _et al._ (2007) G. J. Székely, M. L. Rizzo, and N. K. Bakirov, The Annals of Statistics 35 (2007), 10.1214/009053607000000505.
  * CERN (2017) CERN, “Cern yellow reports: Monographs, vol 2 (2017): Handbook of lhc higgs cross sections: 4. deciphering the nature of the higgs sector,”  (2017).
  * Cho  _et al._ (2020) W. S. Cho, H. D. Kim, and D. Lee, Physical Review D 102 (2020), 10.1103/physrevd.102.115007.
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” (2019), arXiv:1912.01703 [cs.LG] .
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , the Journal of machine Learning research 12, 2825 (2011).
  * Cowan  _et al._ (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, The European Physical Journal C 71, 1 (2011).
  * Gallicchio and Schwartz (2011) J. Gallicchio and M. D. Schwartz, Physical Review Letters 107 (2011), 10.1103/physrevlett.107.172001.
  * Gallicchio  _et al._ (2011) J. Gallicchio, J. Huth, M. Kagan, M. D. Schwartz, K. Black, and B. Tweedie, Journal of High Energy Physics 2011, 1 (2011).
  * Catani  _et al._ (1992) S. Catani, G. Turnock, and B. R. Webber, Phys. Lett. B 295, 269 (1992).
  * Larkoski  _et al._ (2013) A. J. Larkoski, G. P. Salam, and J. Thaler, Journal of High Energy Physics 2013 (2013), 10.1007/jhep06(2013)108.
  * Dasgupta  _et al._ (2013) M. Dasgupta, A. Fregoso, S. Marzani, and G. P. Salam, Journal of High Energy Physics 2013 (2013), 10.1007/jhep09(2013)029.
  * Gallicchio and Schwartz (2010) J. Gallicchio and M. D. Schwartz, Phys. Rev. Lett. 105, 022001 (2010), arXiv:1001.5027 [hep-ph] .


# dynamicMF: A Matrix Factorization Approach to Monitor Resource Usage in High
Performance Computing Systems

Niyazi Sorkunlu
Computer Science and Engineering
University at Buffalo,
State University of New York
Buffalo, New York 14260
<EMAIL_ADDRESS>
Duc Thanh Anh Luong
Computer Science and Engineering
University at Buffalo,
State University of New York
Buffalo, New York 14260
<EMAIL_ADDRESS>Varun Chandola
Computer Science and Engineering
University at Buffalo,
State University of New York
Buffalo, New York 14260
<EMAIL_ADDRESS>
###### Abstract

High performance computing (HPC) facilities consist of a large number of
interconnected computing units (or nodes) that execute highly complex
scientific simulations to support scientific research. Monitoring such
facilities, in real-time, is essential to ensure that the system operates at
peak efficiency. Such systems are typically monitored using a variety of
measurement and log data which capture the state of the various components
within the system at regular intervals of time. As modern HPC systems grow in
capacity and complexity, the data produced by current resource monitoring
tools is at a scale that it is no longer feasible to be visually monitored by
analysts. We propose a method that transforms the multi-dimensional output of
resource monitoring tools to a low dimensional representation that facilitates
the understanding of the behavior of a High Performance Computing (HPC)
system. The proposed method automatically extracts the low-dimensional signal
in the data which can be used to track the system efficiency and identify
performance anomalies. The method models the resource usage data as a three
dimensional tensor (capturing resource usage of all compute nodes for
difference resources over time). A dynamic matrix factorization algorithm,
called dynamicMF, is proposed to extract a low-dimensional temporal signal for
each node, which is subsequently fed into an anomaly detector. Results on
resource usage data collected from the Lonestar 4 system at the Texas Advanced
Computing Center show that the identified anomalies are correlated with actual
anomalous events reported in the system log messages.

##  1 Introduction

High performance computing (HPC) is at the forefront of scientific discovery
and engineering innovation. To support the growing computing demands in these
areas, modern HPC systems have seen a rapid transformation in terms of compute
capacity and the underlying complexity. To ensure that an HPC system is
operating at its peak efficiency, analysts typically rely on system monitoring
tools to monitor the system performance and identify anomalies.

Several tools exist for collecting and visualizing resource usage data from
large scale HPC installations (e.g. Texas Advanced Computing Center TACC_Stats
[10], XSEDE Metrics on Demand or XDMoD [17], etc.). Such tools can produce
large amounts of high dimensional resource usage data at a high temporal
frequency for each computational node in the system. The data collected by
such tools provides a real-time view of the performance of the system which is
typically fed into a visual interface (e.g., XDMoD [17], OVIS [3], Bright
Cluster [4], etc.). However, such data is typically large and unwieldy, and
visual monitoring is often challenging and inefficient.

Existing automated methods typically monitor each system resource or state for
every compute node independently for potential deviations or anomalies [18],
relying on a pre-defined threshold. Such methods can identify only those
anomalous scenarios in which an individual node exhibits significant deviation
for individual resources. These methods often miss anomalies that are under
the threshold. Such anomalies, are weakly manifested across several nodes and
multiple system resources and can be potentially detected by understanding the
interactions between the different aspects of the system. A tensor
decomposition based method [21] jointly models the interactions across the
nodes, resources, and time, and produces a single time varying signal that can
be used for tracking the overall system performance. However, the output does
not allow for finer resolution analysis, for instance, identifying the
specific nodes that are performing anomalously.

In this paper, we propose a method for understanding the system behavior at
the node resolution. The key assumption here is that the observed system
behavior, captured as resource usage information by tools such as TACC_Stats,
can be decoupled into node and metric specific behaviors. Further, the node
behavior can be decoupled into time-invariant (or static) and time-dependent
(or dynamic) terms. This “decoupling” is achieved using a dynamic matrix
factorization method which operates on the sequence of node-usage matrices,
collected over time (See Figure 1). The output of the proposed algorithm is a
set of low-dimensional representations for the metrics and the nodes that
facilitate understanding of the system in multiple ways. In particular, the
dynamic node representation allows for tracking the performance of each node.
Additionally, we use the proposed algorithm to produce a node level anomaly
statistic and show that many performance anomalies, identified from message
logs, co-occur with the identified anomalies.

Figure 1: Typical resource usage data collected by tools such as TACC_Stats.
Each panel corresponds to a compute node and consists of time-varying usage
metrics for a variety of system resources.

The paper is organized as follows. The key problem addressed in this paper is
outlined in Section 2, along with the proposed model, dynamicMF. Results on a
week long TACC_Stats data for the Lonestar 4 system in the Texas Advanced
Computing Center are discussed in Section 3. Section 4 provides a brief
overview of related works in this area and conclusion is presented in Section
5.

##  2 Proposed Methodology

We describe the proposed dynamicMF algorithm in this section. The algorithm is
motivated by the fact that, while the performance of an HPC system is captured
for many nodes and in the context of several resources, the metrics and the
nodes exhibit a clustering pattern: several resources produce correlated usage
metrics and many nodes behave similarly, either due to their computing
specifications or the workload.

###  2.1 Problem Setting

We consider an HPC system consisting of $N$ compute nodes. A resource usage
monitoring program, such as TACC_Stats periodically reports usage data
corresponding to $M$ different resource usage metrics (CPU consumption, memory
usage, network usage are three typical examples). At a given time $t$, the
variable $z_{nmt}$ denotes the $m^{th}$ resource’s usage for node $n$ in the
time interval ending at time $t$. We collect all of these variables at time
$t$ in a matrix ${\bf Z}_{t}\in\mathbb{R}^{N\times M}$, such that
$Z_{t}[n,m]=z_{nmt}$.

###  2.2 Proposed Model

The proposed model assumes that the nodes and the metrics can be represented
as vectors in two distinct (and unobserved/latent) $K$-dimensional spaces
($K\ll N,M$), respectively. Thus each metric is represented as a
$K$-dimensional vector, ${\bf v}_{m}\in\mathbb{R}^{K}$, where each of the $K$
dimensions denote a canonical characteristic of the corresponding metric. For
the nodes, we assume two representations. The first representation is static
over time, denoted as ${\bf\bar{u}}_{n}\in\mathbb{R}^{K}$, and characterizes
the base behavior of the node. The second representation is time-dependent,
denoted as ${\bf\widehat{u}}_{nt}\in\mathbb{R}^{K}$, and captures the dynamic
behavior of the node.

Assuming that the above representations, including $\\{{\bf
v}_{m}\\}_{m=1}^{M}$, $\\{{\bf\bar{u}}_{n}\\}_{n=1}^{N}$,
$\\{{\bf\widehat{u}}_{nt}\\}_{n=1,t=1}^{N,T}$, are known, the observed
behavior for each node and metric combination, i.e., $z_{nmt}$, is generated
by a functional interaction between the three representations. In this paper,
we assume a linear interaction of the following form:

$z_{nmt}=({\bf\bar{u}}_{n}\odot\widehat{\bf u}_{nt})^{\top}{\bf v}_{m}$ (1)

where $\odot$ is an element-wise product. For the entire resource usage
matrix, the interaction can be written as:

${\bf Z}_{t}=({\bf\bar{U}}\odot{\bf\widehat{U}}_{t}){\bf V}^{\top}$ (2)

${\bf V}\in\mathbb{R}^{M\times K}$ consists of the latent representations of
the $M$ metrics, i.e., ${\bf V}=[{\bf v}_{1},{\bf v}_{2},\ldots{\bf
v}_{M}]^{\top}$. ${\bf\bar{U}}\in\mathbb{R}^{N\times K}$ consists of the
static latent representations of the $N$ nodes, i.e, ${\bf\bar{U}}=[\bar{\bf
u}_{1},\bar{\bf u}_{2},\ldots{\bf\bar{u}_{N}}]^{\top}$.
${\bf\widehat{U}}_{t}\in\mathbb{R}^{N\times K}$ consists of the dynamic latent
representations of the $N$ nodes at time $t$, i.e,
${\bf\widehat{U}_{t}}=[\widehat{\bf u}_{1t},\widehat{\bf
u}_{2t},\ldots{\bf\widehat{u}_{Nt}}]^{\top}$. Figure 2 illustrates the
interaction between components of the proposed model to generate the
observation for each metric and node pair in a given time window.

###  2.3 Inference Problem

The proposed generative model (See Figure 2) shows how the observations are
generated through interactions between the three low-dimensional
representations. However, given that these representations are unknown, the
actual problem here is that of inference, i.e., given the observed resource
usage data, $\\{{\bf Z}_{t}\\}_{t=1}^{T}$, estimate ${\bf V},\bar{\bf U}$, and
$\\{\widehat{\bf U}_{t}\\}_{t=1}^{T}$.

Each of the inferred entities provides a distinct insight into the system
performance. The latent metric characteristics (${\bf v}_{m}$) map the metrics
into a common space, where one can understand the similarity between the
various metrics as inferred from the data. For instance, one would expect all
metrics pertaining to the network resource appear similar in the latent space.
However, how would input/output related metrics compare with the network
metrics? Such relations are quantifiable using the latent space
representation. Similar insights can be derived for the nodes using the static
node characteristics ($\bar{\bf u}_{n}$).

${\bf z}_{nmt}$$\widehat{\bf u}_{nt}$$\bar{\bf u}_{n}$${\bf v}_{m}$$M$$T$$N$
Figure 2: Illustration of the assumed interactions between the observed data
(${\bf z}_{nmt}$) and the latent variables ($\widehat{\bf u}_{nt},\bar{\bf
u}_{n},{\bf v}_{m}$) at a given time $t$ for node $n$ and metric $m$ in the
proposed dynamicMF model.

However, the dynamic node characteristics ($\widehat{\bf u}_{nt}$) is the key
output of the proposed model that allows for tracking the behavior of each
node over time using a few latent characteristics. The information in matrices
$\\{\widehat{\bf U}_{t}\\}_{t=1}^{T}$ is obtained from the observed resource
usage data after “explaining away” the static node and metric behavior. Thus,
the residual $K$-dimensional signal is expected to better provide information
about the true dynamic behavior of each node.

###  2.4 Estimating Latent Representations

We pose the problem of estimating the latent entities, ${\bf V},\bar{\bf U}$,
and $\\{\widehat{\bf U}_{t}\\}_{t=1}^{T}$, as an optimization problem that
minimizes the error between the observed data ($\\{{\bf Z}_{t}\\}_{t=1}^{T}$)
and reconstruction obtained in (2). The optimization problem can be written
as:

$\displaystyle\min\sum_{t=1}^{T}\left\lVert{\bf Z}_{t}-(\bar{\bf
U}\odot\widehat{\bf U}_{t}){\bf V}^{T}\right\rVert_{F}^{2}$ (3) with respect
to $\bar{\bf U},{\bf V},\\{\widehat{\bf U}_{t}\\}_{t=1}^{T}$

where $\left\lVert\cdot\right\rVert_{F}$ represents the Frobenius norm.

###  2.5 Optimization

In order to solve the optimization problem as stated in (3), we use Adam
algorithm [15] as the main building block to jointly optimize over the entire
set of variables. The Adam algorithm performs first-order gradient-based
optimization of stochastic objective functions, based on adaptive estimates of
lower-order moments, and has been highly successful in the training of deep
neural networks. Typically, as a gradient-based method, one step of Adam
algorithm only requires one pass over the data. In addition, its optimization
with momentum allows the algorithm to work well for mini-batches of data
instead of the entire dataset. Beside its fast computation, this algorithm is
also efficient in terms of memory usage and well-suited for training model
with many parameters. Algorithm 1 provides the pseudo-code for our approach.

Initialization: initialize the values of $\bar{\bf U}$, ${\bf V}$ and
$\\{\widehat{\bf U}\\}_{t=1}^{T}$

for $iter\in\\{1,\cdots,max\\_iter\\}$ do

Update $\bar{\bf U},{\bf V},\\{\widehat{\bf U}_{t}\\}_{t=1}^{T}$
simultaneously using one step of Adam optimizer

end for

Algorithm 1 Adam Optimizer

For the experiments conducted in this paper, we implement the dynamicMF
algorithm using TensorFlow [1]. The hyper-parameters of the Adam optimizer are
set by following the suggestion in the original paper [15], i.e.
$\alpha=0.001$ (step size), $\beta_{1}=0.9$, $\beta_{2}=0.999$ (decay rate in
exponential moving average).

###  2.6 Anomaly Detection

By using the resulting latent values $\bar{\bf U}$, $\\{\hat{\bf
U}_{t}\\}_{t=1}^{T}$ and ${\bf V}$, we can derive an anomaly score for node
$n$ at time step $t$ which can be defined as a discrepancy between model and
observation, i.e.:

${\bf
a}_{nt}=\frac{1}{m}\sum_{m}^{M}\left|z_{nmt}-({\bf\bar{u}}_{n}\odot\widehat{\bf
u}_{nt})^{\top}{\bf v}_{m}\right|$ (4)

##  3 Experimental Results

In this section, we present the experimental results by applying the
dynamicMFalgorithm to analyze real resource usage logs for a large HPC system,
and identify potential system anomalies.

###  3.1 Dataset

The dataset used in the experiments is a subset of data obtained from Lonestar
4 cluster located at the Texas Advanced Computing Center
(TACC)111https://portal.tacc.utexas.edu/archives/lonestar4. This cluster
consists of $1888$ computing nodes from which the records of $1709$ nodes are
used in this study. The resource usage data is collected using the TACC_Stats
system monitor [10] which records various resources usage statistics at each
computational node every $10$ minutes. In our experiments, we use a set of
$86$ resource usage statistics with a resolution of $10$ minutes from 1:10:01
March $1^{st}$ 2013 to 23:40:01 March $7^{th}$ 2013\. Totally, we have data of
1000 timesteps that will be used for analysis in our experiments. See table 1
for the complete list of resource usage statistics used in this paper. As
different statistics have different scale and units, we have normalized them
to have zero mean and unit variance.

Component |  Resource usage metrics
---|---
CPU |  user, nice, system, idle, iowait, irq, softirq
I/O |  rd_ios, rd_merges, rd_sectors, rd_ticks, wr_ios, wr_merges, wr_sectors, wr_ticks
Lustre /scratch, /work |  read_bytes, write_bytes, dirty_pages_hits, dirty_pages_misses, ioctl, open, close, mmap, seek, fsync, setattr, truncate, getattr, statfs, alloc_inode, setxattr, getxattr, listxattr, inode_permission, readdir, lookup, link, unlink, rename
Lustre network |  tx_msgs, rx_msgs, tx_bytes, rx_bytes,
Virtual Memory |  pgpgin, pgpgout, pgalloc_normal, pgfree, pgactivate, pgdeactivate, pgfault, pgmajfault, pgrefill_normal, pgsteal_normal, pgscan_kswapd, pgscan_direct, pginodesteal, slabs_scanned, kswapd_steal, kswapd_inodesteal, pageoutrun, allocstall, pgrotated
Table 1: List of resource usage metrics used for experiments.

In addition to resource usage statistics, we also have information of
computing jobs that performed in Lonestar 4. In particular, for each computing
node that the job was assigned to, the system records its start time and end
time. There is also the set of system log messages that outputs log traces
from programs running on the HPC system, including the messages from the
applications running on the nodes.

###  3.2 Optimization performance and effect of different K

The performance of the optimization algorithm (See Algorithm 1) in estimating
the latent representations is shown in Figure 3. Results are shown for $K$ =
3, 5, and 10. For each $K$, the objective function decays with each
optimization step and eventually converges.

Figure 3: Objective function (See (3)) optimization for different number of
latent dimensions ($K$). Inset shows the average final objective function
value (or error) for $z_{nmt}$. (Best viewed in color.)

As expected, using a larger number of latent dimensions results in a better
fit. However, larger $K$ also results in a larger number of latent variables
to monitor. In this paper, we choose $K=10$ to report the subsequent findings.
Note that the average error between the model and observations is
approximately $0.57$ for $K=10$ (See Figure 3 inset).

###  3.3 Analyzing latent representations

In this section we illustrate how the latent representations obtained using
dynamicMF can be used to gain insights about the underlying system. The first
latent representation (${\bf V}$) provides a characterization of the metrics.
A two-dimensional visualization of the metric representation is shown in
Figure 4. We use Principal Component Analysis (PCA) [14] to reduce the
$K$-dimensional data to two dimensions for visualization.

Figure 4: 2-dimensional view of the latent representation for metrics
estimated by dynamicMFwith $K=10$. The $K$-dimensional data was reduced to two
dimensions using PCA. Metrics belonging to a single coarse grouping (llite,
vm, cpu, cpu, lnet, block_sda) are displayed in the same way. (Best viewed in
color.)

Figure 4 reveals several interesting insights regarding the characteristics of
various metrics. While metrics within most coarse categories appear to cluster
together in the latent space (e.g., llite, vm), there are some clear outliers
(e.g., llite_work_listxattr, vm_pgmajfault).

In the same way, the static latent representations for the nodes in the matrix
${\bf\bar{U}}$ can be used to visualize the nodes. Figure 5 shows the two-
dimensional plot for the static node behavior. An interesting pattern with two
distinct and equal sized clusters is revealed. While a thorough post-analysis
of the node specifications is needed to explain the patterns, we distinguish
between the regular compute nodes with 24 GB memory and the large memory nodes
with 1 TB memory. The large memory nodes are mapped in the periphery of two
clusters, indicating that dynamicMF is able to capture the specification
related characteristics of the nodes.

Figure 5: 2-dimensional view of the static latent representation for the
compute nodes estimated by dynamicMF with $K=10$, using PCA for reduction to
two dimensions. Large memory nodes are displayed in a different way. (Best
viewed in color.)

As mentioned earlier, the key output that is produced by dynamicMF is the
dynamic latent representation for each node. An important point here is that
each latent dimension contains unique information about the nodes. While we do
not explicitly enforce orthogonality in the optimization formulation, the
results show that the latent dimensions are not significantly correlated (See
Figure 6).

Figure 6: Pearson R correlation among the 10 latent node representations.
(Best viewed in color.)

Figure 7 shows the $K$ components of the time-varying ${\widehat{u}_{nt}}$ for
two different nodes in the system. For both nodes, the different latent
dimensions reveal a unique behavior of the node. Moreover, the output for the
two nodes is significantly different. This is due to the fact that the two
nodes had significantly different workloads during the target week (obtained
from the job information present in the TACC_Stats data). Node c325-312 (blue
plot) had a significantly light workload with 9 executed jobs, and
consequently the latent dynamic behavior has low variance. On the other hand,
node c336-203 (red plot) had a heavy workload with 60 executed jobs, which is
reflected in the high variance for the latent dimensions. Figure 7 shows 10
latent dimensions of the time-varying component ${\widehat{u}_{nt}}$ for two
different nodes with two different workloads in the system. In particular,
node c325-312 (blue plot) had a significantly light workload with seven jobs
receiving from users. On the other extreme, node c336-203 (red plot) had the
most active jobs during the target week with $60$ jobs in total. In Figure 7,
these two nodes exhibit similar patterns of having sporadic spikes at some
certain time steps. Although it looks counterintuitive to have similar
patterns for nodes with two different workloads, one may note that even for
node with zero job submitted from users in HPC system, it still has system-
level jobs running in the background which make the node have different
resource usage at different time, causing sporadic spikes in the temporal
behavior as seen in the figure. In addition, what we observe in Figure 7 is
the latent temporal representation of nodes after “explaining away” the
overall time-independent component ${\bf\bar{u}}_{n}$. For this reason, spikes
appear in the temporal representation as a result of unusual resource
consumptions in each node.

Figure 7: Temporal 10-dimensional representation estimated by dynamicMFfor two
selected nodes. Node c325-312 was inactive during the target week with seven
executed jobs, shown in blue plot. Node c336-203 was the most active during
the target week with 60 executed jobs, shown in red plot.. (Best viewed in
color.)

###  3.4 Finding anomalies using dynamicMF

In this section, we analyze the association of our resulting hidden values and
common errors obtained from system log messages. In particular, there are
three frequent error types that we extracted from system log messages: (1)
write error, (2) segmentation fault and (3) inode error. Write error is a type
of error that happens in the system log with the format “ost_write operation
failed with $-122$”. This is a linux error that happens when disk quota is
exceeded222https://github.com/torvalds/linux/blob/master/include/uapi/asm-
generic/errno.h. Segmentation fault is a common error type when a program
tries to write or read in an invalid memory location. Inode error happens when
ll_inode_revalidate_fini() fails and returns with error code $-43$. We
associate the error types from system log with our resource usage statistics
by assigning the error that occurs within 10 minutes of the occurence of the
present time window.

In our analysis, we use the anomaly score defined in section 2.6 and analyze
its correspondence to the frequent system log error types mentioned above. To
compare, we also show results using an adapted anomaly detection method that
constructs a three-way tensor using the resource usage data and then uses
tensor reconstruction error as the anomaly score for the entire system [21].
While the original method provides a single score for the whole system, we
adapt the method to produce a score for each node individually.

Figure 8 shows three typical examples of three computing nodes for each system
log error type. In this figure, for write error, we do not observe any
association between write error and anomaly score. As the write error happens
as a result of exceeding disk quota, this kind of erratic behavior does not
reflect in resource consumption and therefore it does not reveal any
association with our residual values, as can be seen in Figure 8a. On the
other hand, the segmentation faults in Figure 8b and inode errors in Figure 8c
occur very close and even identical to the peaks of the anomaly score. As
segmentation fault and inode errors are typical examples of buggy programs,
these errors are usually results of faulty programming practice. This can lead
to the programs with inefficient resource usage, which can be subsequently
observed by our resource usage statistics. Therefore, in Figure 8b and 8c, the
spikes in anomaly score values often appear roughly around the errors in terms
of time.

The tensor based method (red plot in Figure 8) produces a metric shows high
anomaly values that are temporally close to some of the errors. However, the
method also produces a large number of false alarms, which can make such a
method practically unusable.

(a) write error

(b) segmentation fault

(c) inode error

Figure 8: Association between anomaly scores computed using (4) and errors in
system log messages. Anomaly scores produced by dynamicMFare shown in blue and
the anomaly score produced by the state of art tensor reconstruction method is
shown in red. Vertical gray lines indicate where the errors happen. The node
name is shown above each figure. (Best viewed in color.)

##  4 Related Work

The dynamic matrix factorization algorithm proposed here falls under the
general purview of unsupervised representation learning [2], that includes a
variety of tasks such as dictionary learning, independent component analysis,
autoencoders, matrix factorization and various forms of clustering. In
particular, a relevant topic is sparse coding [7], where the objective is to
learn over-complete and sparse bases for a given data set. However, methods in
these categories primarily focus on a static representation (typically as a
matrix) and do not handle temporally evolving data.

The idea of dynamic matrix factorization has been used in designing the
recommendation system. In particular, dynamic Poisson factorization (dPF)
model works effectively when we assume each observation follows a Poisson
distribution [6]. In this model, because of their probabilistic setting,
variational inference has been used to perform inference. After that, they use
the block coordinate approach to solve the optimization problem yield as a
result of variational inference. Although our model is similar to dPF in the
design of shared components and local components, our model is different from
dPF in several ways: (1) we do not assume the temporal relationship between
consecutive timestep as in dPF, (2) because of our non-probabilistic
formulation, we can avoid the problem of non-conjugacy in dPF and use more
effective gradient-based optimization such as Adam algorithm instead of using
block coordinate optimization with L-BFGS (a more computationally intensive
optimization method) as a building block.

Solutions for detecting, diagnosing, and predicting faults and failures in
large high performance computing installations have typically relied on
message logs [22, 16, 11, 20, 19] or resource usage data [12, 5], or both [9,
13, 8]. Since this paper focuses on the detection task, we present a brief
overview of related methods that deal with detecting faults. Methods operating
on message logs typically aggregate message logs by an entity of interest,
e.g., a computational node/block [16, 22] or a job [11], and identify
anomalous entities using different data representations, such as message
arrival statistics [16], vector representation derived from the message
content [22] or a state machine that models the dynamic behavior of the entity
[11]. However, message logs are typically noisy and often incomplete, which
has led to methods that analyze alternate data sources such as resource usage
metrics or performance counter data [12, 5]. However, these solutions perform
a node-specific or job-specific analysis of resource usage to identify
anomalous nodes [12] or jobs [5]. Recently, solutions that combine resource
usage data and message logs to improve fault detection have been proposed [13,
8]. The Crude system uses resource usage data to improve the performance of a
PCA driven anomaly detection method [22] that operates at node and job level.
However, none of the existing methods model the temporal dimension to better
identify faults. In a recent work, a tensor based representation of the
resource usage data was used to identify anomalies at the system level [21].
The method relies on low rank approximation of a three way tensor that
captures interactions between nodes, metrics, and time, to establish an
anomaly score for the entire system at every time instance. However, the
method does not provide a node level estimate of system behavior and
anomalies.

##  5 Conclusions and Future Directions

In this study, we propose dynamicMF \- a dynamic matrix factorization method
which takes both global and local behaviors of the system into account to
analyze for variations. This method is applied to a list of time-varying
resource usage matrices. In our experiments with real data from Lonestar 4
system, dynamicMF helps reduce the dimension of resource usage statistics and
provide visualization at a node-specific level. The proposed anomaly statistic
allows for an easy identification of performance anomalies. Future direction
of our research will be on developing a fully-automated anomaly detection
method to detect failure in HPC system by using the low dimensional latent
space obtained from dynamicMF.

Future direction of our research will be to leverage job information for
analysis and incorporate it into our model. Successfully modeling the
interaction between jobs and nodes can help us better explain the data and
obtain a more accurate model to detect inefficiencies in resource usage.

## References

  * [1] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al.  Tensorflow: A system for large-scale machine learning.  In OSDI, volume 16, pages 265–283, 2016.
  * [2] Y. Bengio, A. C. Courville, and P. Vincent.  Unsupervised feature learning and deep learning: A review and new perspectives.  CoRR, abs/1206.5538, 2012.
  * [3] J. Brandt, A. Gentile, J. Mayo, P. Pebay, D. Roe, D. Thompson, and M. Wong.  Resource monitoring and management with ovis to enable hpc in cloud computing environments.  In 2009 IEEE International Symposium on Parallel Distributed Processing, pages 1–8, 2009.
  * [4] https://www.microway.com/download/datasheet/Bright_Cluster_Manager_Brochure.pdf, 2017\.
  * [5] G. Bronevetsky, I. Laguna, B. R. de Supinski, and S. Bagchi.  Automatic fault characterization via abnormality-enhanced classification.  In Proceedings of the 2012 42Nd Annual IEEE/IFIP International Conference on Dependable Systems and Networks (DSN), pages 1–12, 2012.
  * [6] L. Charlin, R. Ranganath, J. McInerney, and D. M. Blei.  Dynamic poisson factorization.  In Proceedings of the 9th ACM Conference on Recommender Systems, pages 155–162. ACM, 2015.
  * [7] S. S. Chen, D. L. Donoho, and M. A. Saunders.  Atomic decomposition by basis pursuit.  SIAM Rev., 43(1):129–159, 2001.
  * [8] E. Chuah, A. Jhumka, J. C. Browne, N. Gurumdimma, S. Narasimhamurthy, and B. Barth.  Using message logs and resource use data for cluster failure diagnosis.  In 23rd annual IEEE International Conference on High Performance Computing, Data, and Analytics (HiPC 2016), 2016.
  * [9] E. Chuah, A. Jhumka, S. Narasimhamurthy, J. Hammond, J. C. Browne, and B. Barth.  Linking resource usage anomalies with system failures from cluster log data.  In Proceedings of the 2013 IEEE 32Nd International Symposium on Reliable Distributed Systems, SRDS ’13, pages 111–120, Washington, DC, USA, 2013\. IEEE Computer Society.
  * [10] T. Evans, W. Barth, J. Browne, R. DeLeon, T. Furlani, S. Gallo, M. Jones, and A. Patra.  Comprehensive resource use monitoring for HPC systems with TACC stats.  In International Workshop on HPC User Support Tools, pages 13–21, 2014.
  * [11] Q. Fu, J.-G. Lou, Y. Wang, and J. Li.  Execution anomaly detection in distributed systems through unstructured log analysis.  In International conference on Data Mining, 2009.
  * [12] Q. Guan, D. Smith, and S. Fu.  Anomaly detection in large-scale coalition clusters for dependability assurance.  In 2010 International Conference on High Performance Computing, pages 1–10, 2010.
  * [13] N. Gurumdimma, A. Jhumka, M. Liakata, E. Chuah, and J. Browne.  Crude: Combining resource usage data and error logs for accurate error detection in large-scale distributed systems.  In 2016 IEEE 35th Symposium on Reliable Distributed Systems (SRDS), pages 51–60, 2016.
  * [14] I. Jolliffe.  Principal Component Analysis.  Springer Verlag, 1986.
  * [15] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [16] A. J. Oliner, A. Aiken, and J. Stearley.  Alert detection in system logs.  In 2008 Eighth IEEE International Conference on Data Mining, pages 959–964, Dec 2008.
  * [17] J. Palmer, S. Gallo, T. Furlani, M. Jones, R. DeLeon, J. White, N. Simakov, A. Patra, J. Sperhac, T. Yearke, R. Rathsam, M. Innus, C. Cornelius, J. Browne, W. Barth, and R. Evans.  Open XDMoD: A tool for the comprehensive management of high-performance computing resources.  Computing in Science Engineering, 17(4):52–62, 2015.
  * [18] M. Peiris, J. H. Hill, J. Thelin, S. Bykov, G. Kliot, and C. Konig.  Pad: Performance anomaly detection in multi-server distributed systems.  In Proceedings of the 2014 IEEE International Conference on Cloud Computing, pages 769–776, 2014.
  * [19] A. Pelaez, A. Quiroz, J. C. Browne, E. Chuah, and M. Parashar.  Online failure prediction for hpc resources using decentralized clustering.  In 2014 21st International Conference on High Performance Computing (HiPC), pages 1–9, 2014.
  * [20] T. Reidemeister, M. A. Munawar, M. Jiang, and P. A. S. Ward.  Diagnosis of recurrent faults using log files.  In Proceedings of the 2009 Conference of the Center for Advanced Studies on Collaborative Research, pages 12–23, 2009.
  * [21] N. Sorkunlu, V. Chandola, and A. K. Patra.  Tracking system behavior from resource usage data.  In 2017 IEEE International Conference on Cluster Computing, CLUSTER 2017, Honolulu, HI, USA, September 5-8, 2017, pages 410–418, 2017\.
  * [22] W. Xu, L. Huang, A. Fox, D. Patterson, and M. Jordan.  Mining console logs for large-scale system problem detection.  In Proceedings of the Third Conference on Tackling Computer Systems Problems with Machine Learning Techniques, SysML’08, pages 4–4, Berkeley, CA, USA, 2008. USENIX Association.