# Commissioning and installation of the new small-Diameter Muon Drift Tube
(sMDT) detectors for the Phase-I upgrade of the ATLAS Muon Spectrometer

G. H. Eberwein, O. Kortner, H. Kroha, M. Rendel, P. Rieck, D. Soyk,
E. Voevodina1 and V. Walbrecht on behalf of the ATLAS Muon Group Manuscript
received December 1, 2021.1 Corresponding author. (e-mail:
voevodin@mpp.mpg.de).G. H. Eberwein was with Max-Planck-Institute for Physics,
Foehringer Ring 6, 80805 Munich, Germany. He is now with Department of
Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford
OX1 3RH, UK.O. Kortner, H. Kroha, E. Voevodina, M. Rendel, P. Rieck, D. Soyk
and V. Walbrecht are with Max-Planck-Institute for Physics, Foehringer Ring 6,
80805 Munich, Germany.

###### Abstract

The Monitored Drift Tubes, as a part of the ATLAS muon spectrometer, are
precision drift chambers designed to provide excellent spatial resolution and
high tracking efficiency independent of the track angle. Through the life of
the LHC and ATLAS experiment, this detector has already demonstrated that they
provide precise tracking over large areas. The aim of the ATLAS muon
spectrometer upgrade is to increase the muon trigger efficiency, precise muon
momentum measurement and to improve the rate capability of the muon system in
the high-background regions during the High-Luminosity LHC runs. To meet these
requirements, the proposed solution is based on the small (15 mm) diameter
Muon Drift Tube chamber (sMDT) technology. The new detector provides about an
order of magnitude higher rate capability and allows for the installation of
additional new triplet Resistive Plate Chambers (RPCs) trigger detectors in
the barrel inner layer of the muon system. A pilot project for the barrel
inner layer upgrade is underway during the 2019/21 LHC shutdown. For this
reason, the Max-Planck-Institut für Physik in Munich has built 16 sMDT
chambers, each will cover an area of about 2.5 $m^{2}$. To ensure their proper
operation in the experiment, the sMDT detectors have to pass a set of
stringent tests both at the production site and after their delivery at CERN.
After their installation in the ATLAS muon spectrometer, the muon stations are
further tested and commissioned with cosmic rays. The author will describe the
detector design, the quality assurance and certification path, as well as will
present the experience with the chamber tests, the integration procedure and
installation of the muon stations in the ATLAS experiment.

##  I Introduction

The success of the muon spectrometer is based on the high muon trigger
efficiency and muon momentum resolution up to the TeV scale. However, in the
innermost transition region between the ATLAS barrel and the endcap
($1.0<|\eta|<1.3$), a non-negligible rate of the fake muon trigger has been
observed due to the high rate based on secondary charged tracks generated by
beam halo protons and a lack of detector instrumentation. To reinforce the
fake rejection and to improve the selectivity of the muon trigger in this
problematic region, the two new generations of the gaseous detector
technology, such as the triplet thin gap Resistive Plate Chambers (RPC) and
small (15 mm) diameter Muon Drift Tube (sMDT) detectors, will be installed
during the LHC Long Shutdown 2 (LS2, 2020-2021). That Phase-I upgrade serves
as a pilot project, known as BIS78 (shown in Fig. 1), for the complete
replacement of the existing 96 BIS1-6 MDT chambers in the LS3 period
(2025-2027) by the combined sMDT+RPC muon modules that will have a very
similar design[1]. The first 16 new BIS78 sMDT detectors are currently under
installation and commissioning at CERN: this smaller size project mainly aims
to validate the final mechanical installation procedure, the correctness of
the services previously installed in the experimental cavern, and to obtain
early feedback about the front-end electronics response and detection
performance in the ATLAS environment with the detectors powered and read out
through the final services.

Figure 1: Schematic of the BIS78 project for the Phase-1 upgrade of the ATLAS
muon spectrometer.

##  II sMDT Detector Design

Each individual BIS78 muon station serves as the integrated tracking and
triggering detector and therefore composes of one sMDT chamber with two
different sizes of RPC triplets (BIS7 and BIS8). Furthermore, the sMDT
detectors consist of two multilayers of the high pressurized drift tubes, each
containing four layers of densely packed tubes mounted on an aluminum support
frame: the chamber sizes vary from about 2.67 $m^{2}$ to 3.04 $m^{2}$ while
tube lengths range from 1009 mm to 1669 mm. In order to uniformly fill the
space by taking into account the features of the toroid magnet coil at the
ends of the barrel inner layer and to reduce the number of dead areas, 12
different sMDT detectors versions have been properly designed in MPI Munich
and approved by ATLAS Muon Community. For detector control data, the BIS78
sMDT chambers are equipped with 12 temperature sensors distributed on the
multilayer and the support structure. On the top surface of the detector are
typically mounted eight platforms: four Axial-Praxial, two CCC and two B-field
platforms. In addition, each detector is instrumented with the in-plane
alignment system, consisting of two image sensors, four lenses, and two or
four masks, in order to monitor the deformations such as sag and torsion. For
the 16 A/C-side BIS78 sMDT detectors, the chamber-to-chamber alignment system
and B-field components will be re-used from the previous MDT detectors[2].

##  III Validation Tests and Integration at CERN

Upon sMDT chambers arrival at CERN BB5 construction facility, each detector is
subjected to a strict validation test program: gas tightness, connectivity,
noise rate measurement of the front-end electronics and test with cosmic rays,
in order to detect any damages which might have occurred during the
uncontrolled conditions of the transport from the MPI Munich production site
and to ensure their proper operation in the ATLAS muon spectrometer. A picture
of the sMDT detector during the commissioning test campaign at CERN BB5
facility is shown in Fig. 2.

Figure 2: Picture of the new sMDT chamber during the commissioning campaign at
CERN BB5 facility.

The leak rate of each multilayer is determined from the temperature corrected
pressure drop over a time period of at least 24 hr. In order to avoid the
possible gas contamination in the chamber, since it significantly can change
the maximum drift time, detectors are evacuated and then filled with the
operating $Ar/CO_{2}$ (93:7) gas mixture at the pressure 3 bar. The maximum
acceptable gas leak rate is $2N\times 10^{-8}$ bar$\times$l/s for N tube in a
detector. The measurement accuracy achievable after 24 hr obtained during the
commissioning campaign is 0.063 mbar/h. The high voltage stability at the
nominal operating high voltage (HV) point 2.73 kV is tested by measuring the
dark current of each tube layer. A noise rate test is performed without and
with HV to easily distinguish between electronic noise and tubes with internal
discharges. The average intrinsic noise rate of all tested chambers has been
found of $(0.10\pm 0.02)$ kHz/tube at the operating threshold value of -39 mV
and at the working point of 2.73 kV, well below the upper noise limit of the
sMDT detectors installed in the ATLAS muon system, that is above 5 kHz/tube.
Moreover, the control interface of the front-end electronics and communication
stability are checked as well. A complete system test of the sMDT detector is
performed with cosmic muons. The maximum drift time spectra, spatial
resolution and muon detection efficiency for each tube of the sMDT chamber are
computed and recorded as well as the detector hit profile is measured to allow
a check of the correct channel assignment.

Figure 3: BIS78 sMDT performance study during the validation tests at CERN BB5
facility.

The MDT chamber under the commissioning test reaches an efficiency of 98.1$\%$
for muons when operated with $Ar/CO_{2}$ (93:7) gas mixture, HV = 2.73 kV and
has a read-out with all mezzanines, configured with -39 mV threshold. The
average chamber efficiency is in an agreement both with the previous QA/QC
test campaign a the MPI production site and with the Monte Carlo (MC)
simulation prediction. The spatial resolution for the BIS78 sMDT detector
under validation test is measured to be (122 $\pm$ 11) $\mu$m with the
multiple scattering correction obtained from Monte Carlo simulation,
$\sigma_{MC}$ = (120 $\pm$ 10) $\mu$m [3]. This performance is fully in
agreement with the previous measurement at the muon beam facility,
$\sigma_{\mu}$ = (125 $\pm$ 4) $\mu$m [4]. Examples of typical results can be
seen in Fig 3.

Figure 4: Photograph of the sMDT and RPC detectors integration process.

After the commissioning campaign, the sMDT chamber is then integrated with its
two triplet RPC trigger detectors to form a so-called final muon station. The
challenges for the mechanical design of the BIS sMDT are the small available
space that needs to share with the new RPC detectors and minimizing the
gravitational deformation of the thin chambers with minimal space available
for the support structures. For this reason, the BIS78 muon stations require a
special common support frame to carry the RPC and sMDT chambers due to their
position on the ATLAS toroid magnet coil. The RPC is inserted in the support
frame and the sMDT chamber is mounted on kinematical bearings. The mechanical
adjustment accuracy is less than 1 mm, fulfilling the requirements on the
chamber alignment in the detector [1].

Figure 5: Final BIS78 muon station with installed B-field and Axial-Praxial,
CCC alignment sensors.

In order to study the effect of the integrated RPC detectors and its
electronics on the sMDT system, the interference test is performed. Then the
sMDT detector is equipped with B-field and Axial-Praxial, CCC alignment
sensors whose proper operation is verified. The photographs of the detectors
integration process and of the final muon station are shown in Fig. 4 and 5.

##  IV Installation and Commissioning of new sMDT stations in the ATLAS
experiment

At the surface above the ATLAS experimental cavern, the mechanical integrity
of the station is checked to ensure that no damages occurred during transport
from the validation and integration site. The new muon station is placed into
an installation frame and lowered into the underground experimental hall.
There, the installation frame is supported by two cranes, rotated to the
appropriate angle, and docked to the rail system of the ATLAS Muon
spectrometer. The muon station is slid onto the rails by using the two
winches, positioned with an accuracy of about a few mm in all three spatial
directions and fixed at one bearing on the rail. A photograph of the BIS78
muon station installation is shown in Fig. 6. Installation of the new BIS78
muon stations for the A-side of the ATLAS muon spectrometer began in mid-
September 2020 and lasted until the end of January 2021.

Figure 6: Photograph of the installation process of the new BIS78 muon station
in the ATLAS muon spectrometer.

The underground commissioning phases are required to bring the new sMDT sub-
system into an operational state suitable for future ATLAS physics data-taking
campaigns. Any potential damages and performance losses are immediately
identified and repaired after installation. The correct sag adjustment is
verified with the on-chamber alignment system. ATLAS global alignment systems
are checked as well. No failures have been observed after installation.
Moreover, the detector’s basic parameters, such as no leaks in the gas, high-
voltage stability of the chamber with the operating gas mixture, low
background noise, appropriate connectivity of the electronic components, and
the operational characteristics of the front-end/back-end electronics, are set
to optimal settings.

The sMDT sub-detectors have been fully added and integrated into the ATLAS
operational Data Acquisition/Detector Control/Data Quality Monitoring systems.

##  V Conclusion

In the coming years, the ATLAS muon spectrometer will go through a series of
modernizations in order to cope with the foreseen increase in LHC performance.
New BIS78 detector stations (the triplet thin gap Resistive Plate Chamber and
small (15 mm) diameter Muon Drift Tube detector) have been constructed for the
Phase-1 upgrade as a pilot project for the innermost barrel region of the
ATLAS muon system to increase the trigger acceptance and to prepare for the
high background rate requirements of the HL-LHC. The Max Planck Institute for
Physics in Munich has built 16 new sMDT chambers. All sMDT detectors passed
the acceptance tests at CERN and were integrated with their RPC trigger
chambers. A total of 8 out of 16 muon stations were successfully installed in
the ATLAS experiment before the start of the Run-3, while the muon system
modernization campaign with another BIS78 modules is continued during the next
third Long Shutdown (LS3) of the LHC.

## References

  * [1] ATLAS Collaboration, _”ATLAS Muon Spectrometer Phase-II Upgrade”_. In: CERN-LHCC-2017-017, ATLAS-TDR-026.
  * [2] H. Kroha, _”Design and construction of integrated small-diameter drift tube and thin-gap resistive plate chambers for the phase-1 upgrade of the ATLAS muon spectrometer”_. In: Nucl. Instrum. Methods A 936, 445 (2019), doi:10.1016/j.nima.2018.10.139.
  * [3] O. Kortner, _”Test-Stand Simulation Programme”_. In: ATL-MUON-2000-021.
  * [4] O. Bittner, _”Construction and Test of a Prototype Chamber for the Upgrade of the ATLAS Muon Spectrometer”_. In: Phys.Procedia 37 (2012) 515-521; doi: 10.1016/j.phpro.2012.02.398.


# Deciphering solar turbulence from sunspots records

F. Plunian1, G. R. Sarson2, R. Stepanov3
1 LGIT, UJF, CNRS, B.P. 53, 38041 Grenoble Cedex 9, France
2 School of Mathematics and Statistics, Newcastle University, Newcastle NE1
7RU, UK
3 Institute of Continuous Media Mechanics, Korolyov 1, 614013 Perm, Russia

(Accepted 2009 …. Received 2009 ….; in original form 2009)

###### Abstract

It is generally believed that sunspots are the emergent part of magnetic flux
tubes in the solar interior. These tubes are created at the base of the
convection zone and rise to the surface due to their magnetic buoyancy. The
motion of plasma in the convection zone being highly turbulent, the surface
manifestation of sunspots may retain the signature of this turbulence,
including its intermittency. From direct observations of sunspots, and
indirect observations of the concentration of cosmogenic isotopes 14C in tree
rings or 10Be in polar ice, power spectral densities in frequency are plotted.
Two different frequency scalings emerge, depending on whether the Sun is
quiescent or active. From direct observations we can also calculate scaling
exponents. These testify to a strong intermittency, comparable with that
observed in the solar wind.

###### keywords:

MHD, turbulence, statistics, sunspots, magnetic fields, plasmas

††pagerange: Deciphering solar turbulence from sunspots
records–References††pubyear: 2009

##  1 Introduction

Sunspots observed at the surface of the convection zone of the Sun are usually
understood as the manifestation of solar magnetic activity. Naked-eye and
telescope observations of sunspots are available from AD 1610, providing
reliable records of sunspot numbers (SSN). Several sets of data exist, varying
in how they have been sampled, averaged (daily or monthly), whether they
concern sunspots or sunspots groups, and on the scientific societies who have
compiled the records. Here we consider the American daily SSN (D)
111ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/
AMERICAN_NUMBERS/RADAILY.PLT, the American daily group SSN (G)
222ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/
GROUP_SUNSPOT_NUMBERS/dailyrg.dat, and the International monthly averaged SSN
(M) 333http://solarscience.msfc.nasa.gov/greenwch/spot_num.txt (figure 1). In
addition SSN at earlier times have been reconstructed from proxies, based on
the concentration of cosmogenic isotopes 14C in tree rings
444ftp://ftp.ncdc.noaa.gov/pub/data/paleo/climate_forcing/solar_variability/
solanki2004-ssn.txt or 10Be in ice core bubbles. The production rate of such
isotopes increases with the cosmic-ray flux, which is higher when the solar
magnetic activity is low. Plotting the SSN versus time reveals a cycle of
about 11 years known as the Schwabe cycle. This cyclical solar magnetic
activity is sufficiently robust to be detected in 10Be concentration records,
even during some long periods with almost no visible sunspots, like the
Maunder minimum (1645-1715) (Beer et al., 1998). Analyzing long time series of
14C and 10Be, it has been shown that the solar activity of the last 70 years
has been exceptionally high (Usoskin et al., 2003; Solanki et al., 2004), and
that a decline is expected within the next two or three cycles (Abreu et al.,
2008).

Figure 1: Data sets of 10 year-averaged SSN from 14C (top blue) and 10Be (top
red), International SSN monthly averaged M (second), American daily group SSN
G (third), and American daily SSN D (bottom).

The occurrence of sunspots is of course an important diagnostic that must be
reproduced by any solar dynamo model. The fact that it is irregular (in spite
of the Schwabe cycle) reflects the complexity inherent in the nonlinear
coupling between the turbulent flow and the magnetic field in the solar
convection zone (Browning et al., 2006). With the long time series of direct
SSN observations and 14C and 10Be data available, it is tempting to calculate
the corresponding frequency spectra, to infer some signature of the underlying
turbulence in the convection zone. Similar attempts have been made for the
Earth (Courtillot & Le Mouel, 1988; Constable & Johnson, 2005; Sakuraba &
Hamano, 2007); the frequency spectrum of the geomagnetic dipole moment
obtained from paleomagnetic data is consistent with an inertial range of
$f^{-5/3}$ scaling, and a dissipation range of $f^{-11/3}$ as expected in
magnetohydrodynamic turbulence (Alemany et al., 2000). For the Sun, analysis
of International daily SSN have led to a $f^{-2/3}$ scaling (Morfill et al.,
1991; Lawrence et al., 1995), and this has been attributed to some sequential
sampling of the field upon arrival at the photosphere on top of a Kolmogorov
spatial scaling due to the underlying turbulence. Here we extend this analysis
to the other SSN records mentioned above and test how robust the $f^{-2/3}$
scaling is, in particular during minima and maxima of sunspot activity. In
addition, for time scales smaller than 2 years, the stochastic character of
the SSN records suggests strong intermittency (Lawrence et al., 1995), as
opposed to a low-dimensional chaotic interpretation. We shall characterize
this intermittency by calculating the corresponding scaling exponents.

##  2 Wavelet spectra

Figure 2: Wavelet coefficients (logarithm of absolute value) of 14C (top), M
(second), G (third), and D (bottom).

In order to filter out the noise, we use a wavelet decomposition of the
signal. In figure 2, the wavelet coefficients are plotted versus time (in
years) and frequency (in month-1). The light (resp. dark) colors correspond to
low (resp. high) values of these coefficients. For a given year, the curve
giving the wavelet coefficient versus frequency corresponds to a power
spectral density of the signal. The dark horizontal stripe for $f\sim 0.01$
month-1 corresponds to the Schwabe cycle. It is not visible in the 14C data
set due to the coarse sampling of this data (averaged over 10 years). On the
other hand a clear dark stripe is visible for $f\sim 3.5\times 10^{-5}$
month-1, corresponding to a $\sim$2400 years cycle. In the figure for M data
the dark stripe of the Schwabe cycle almost disappears during the Dalton
minimum (1790-1820) (Frick et al., 1997). In the G and D data figures we
identify another horizontal stripe at $f\sim 1$ month-1, corresponding to the
solar mean rotation rate. This indicates that the latitudinal repartition of
sunspots is not homogeneous. Finally the vertical light stripes correspond to
minima of magnetic activity.

Figure 3: Wavelet spectral density versus frequency, on average (top), and for
magnetic activity minima (bottom, dashed curves) and maxima (bottom, solid
curves). The dashed straight lines correspond to $-2/3$ and $-8/3$ slopes
(top), and $-2/3$ and $-1$ slopes (bottom).

The time average of the wavelet coefficients are shown in figure 3 (top) for
the five data sets. As mentioned earlier the two peaks around $0.01$ month-1
and $1$ month-1 correspond to the Schwabe cycle and solar mean rotation rate.
Between them the three spectra of the directly observed SSN (M, G, D) present
a common scaling in $f^{-2/3}$. The other data sets 14C and 10Be are
compatible with a $f^{-2/3}$ scaling as well, although the 14C PSD is
overestimated by roughly a factor 10, probably due to the proxy used in the
reconstruction of the SSN from the 14C data.

Figure 4: Same data as in figure 1 (the last three panels) but subdivided into
sets of maximum (dark) and minimum (light) magnetic activity.

For the three data sets (M), (G) and (D), instead of averaging on all times,
we now average on periods corresponding to either maximum or minimum magnetic
activity as shown in figure 4. For the maximum (resp. minimum) activity
subset, the excluded data is that centred around the times of the Schwabe
minima (resp. maxima) The corresponding spectral densities are plotted in
figure 3 (bottom). The slopes for minima are systematically steeper than those
for maxima, indicating two different regimes. To estimate these slopes we vary
both the range of frequency $\left[f_{\min},f_{\max}\right]$ on which they are
calculated, and the way the data sets are split into subsets of maximum and
minimum activity. For the former we take $f_{\max}=0.7$ month-1 for the data
sets (G) and (D) in order to escape from the influence of the peak $f=1$
month-1, and $f_{\max}=0.5$ month-1 for the data set (M). When changing
$f_{\min}$ the slopes change. We vary $f_{\min}$ such that the ratio
$f_{\max}/f_{\min}$ is about 10, and the standard deviation of the slope
remains 10 % or less of its average value. This leads to
$f_{\min}\in\left[0.05,0.07\right]$ for (G),
$f_{\min}\in\left[0.06,0.08\right]$ for (D), and
$f_{\min}\in\left[0.03,0.05\right]$ for (M). In addition we consider at least
three different degrees of splitting for each data set, this splitting degree
being related to the time length of the subsets of maximum and minimum
activity. The choices of this splitting degree are made such that both subsets
are long enough to provide good statistics, but remain separated by sufficient
time lags that the SSN values falling into each subset do not overlap too
much. Then for both the maximum and minimum subsets, we calculate the mean
slope and the standard deviation obtained when varying both the frequency
range and the degree of splitting. The corresponding slope estimates are given
in table 1. They are consistent with power spectra in $f^{-2/3}$ and $f^{-1}$.
The corresponding dashed lines are plotted in figure 3 to guide the eye. The
standard deviations are small, showing that these slopes are robust with
respect to the details of our analysis. The formal standard errors from each
of the individual regressions (for specific frequency ranges and degrees of
splitting) are of comparable magnitude.

Activity | (M) | (G) | (D)
---|---|---|---
Max. | $-0.61\pm 0.05$ | $-0.69\pm 0.05$ | $-0.63\pm 0.05$
Min. | $-1.03\pm 0.03$ | $-0.94\pm 0.07$ | $-0.98\pm 0.02$
Table 1: Slope estimates for the PSD curves plotted in figure 3 (bottom). They
correspond to average values plus standard deviation errors when varying both
the frequency range and the degree of splitting of each data set (into the two
subsets of maximum and minimum magnetic activity).

As noted by Lawrence et al. (1995), the question of causality complicates the
interpretation of such temporal data. The difference of spectral slopes
between minima and maxima can be attributed to two different effects: a change
of the underlying turbulence, affecting the spatial structure the magnetic
field; or a change in the frequency of the sequential sampling of the magnetic
field, as suggested by Lawrence et al. (1995). Although the latter effect
cannot be excluded, there is a simple argument in favour of the former. It is
generally accepted that the occurrence of sunspots at the photosphere is due
to the magnetic buoyancy force $\nabla B^{2}$, where $B$ is some magnetic
induction intensity in the convection zone (Tobias et al., 2001). It is then
tempting to interpret the two spectral slopes as the signatures of this
buoyancy, assuming that the frequency of sunspot occurrence at the photosphere
is proportional to this force. Then the $f^{-2/3}$ and $f^{-1}$ SSN spectra
would correspond to buoyancy spectra of $k^{-2/3}$ and $k^{-1}$, where $k$ is
the spatial wave number. During maxima this implies a Kolmogorov magnetic
energy spectrum of $k^{-5/3}$, compatible with inertia-driven turbulence in
the convection zone. During minima it implies a magnetic energy spectrum of
$k^{-2}$, compatible with turbulence dominated by the solar rotation (Zhou,
1995). In the transport scenario proposed by Tobias et al. (2001), the field
which arises at the surface is the strongest part of a poloidal field
generated by cyclonic turbulence in the convection zone. Our interpretation
then suggests two different regimes for such cyclonic turbulence, controlled
by either inertia or rotation.

##  3 Intermittency

The stochastic nature of SSN occurrence for times scales smaller than 2 years
has been shown by Lawrence et al. (1995), suggesting an intermittent
turbulence. Here our goal is to quantify this intermittency for the three sets
(M, G, D), calculating the corresponding scaling exponents. For that we first
calculate the associated generalized structure function (GSF) and look for its
scaling exponents, as usually done in turbulence. We define the SSN increment
by

$\delta y(t,\tau)=S(t+\tau)-S(t)\;,$ (1)

where $S(t)$ denotes the SSN at time $t$. Assuming statistical stationarity in
the frequency range of interest, the $t$ dependence in $\delta y(t,\tau)$ can
be dropped and the GSF is then given by (Nicol et al., 2008)

$S_{m}(\tau)=\left\langle\left|\delta
y\right|^{m}\right\rangle=\int^{\infty}_{-\infty}\left|\delta
y\right|^{m}P(\delta y,\tau)d(\delta y)\;,$ (2)

where $P$ is the probability density function (PDF) of $\delta y$, the angle
brackets $\left\langle\cdot\right\rangle$ denote time averaging, and $m$ is a
positive integer.

Figure 5: Probability density functions of $S(t+\tau)-S(t)$ for M (top), G
(middle) and D (bottom). The labels indicate the value of $\tau$ in years.

In figure 5 PDFs for the three data sets (M, G, D) are plotted for selected
values of $\tau$. The PDFs of the 14C data are poorly defined, and so we drop
this data set for the rest of the study. The PDFs of the other data sets show
peaks at $\tau=11/2$ years, corresponding to the Schwabe cycle. For other
values of $\tau$ they exhibit tails containing a higher number of rare events
than for a gaussian distribution, suggesting intermittency. Similar results
were shown in Lawrence et al. (1995).

To quantify this intermittency we first check whether the GSF obey a scaling
law in the form

$S_{m}(\tau)\sim\tau^{\zeta(m)}.$ (3)

We find (not shown) that this is clearly the case. In homogeneous and
isotropic fully developed turbulence, intermittency corresponds to
$\zeta(m)<m/3$. The ratio $\zeta(m)/\zeta(3)$ is calculated, estimating the
scaling power of $S_{m}(\tau)/S_{m}(3)$. Plotting the ratio
$\zeta(m)/\zeta(3)$ versus $m$ for the three data sets (figure 6) we see a
clear departure from the Kolmogorov straight line $\zeta(m)/\zeta(3)=m/3$, and
clear indications of intermittency. It is remarkable that the (D) set, which
has the best sampling, leads to the largest intermittency. It is also
remarkable that the (M) and (G) sets lead to similar scaling exponents,
supporting the equivalence between averaging over space and time.

Figure 6: Scaling exponents $\zeta(m)/\zeta(3)$ plotted versus $m$ for the
three sets of data M, G, D, and for 2 months $<\tau<$ 14 months. The dashed
line corresponds to a Kolmogorov scaling $\zeta(m)/\zeta(3)=m/3$.

The exponents can be fitted to the standard $p$-model derived for hydrodynamic
(Meneveau & Sreenivasan, 1987) and magnetohydrodynamic turbulence (Carbone,
1993). This model is defined by

$\zeta(m)=1-\log_{2}\left[p^{m/3}+(1-p)^{m/3}\right].$ (4)

We find $p_{\rm(G)}=0.68$, $p_{\rm(M)}=0.68$, and $p_{\rm(D)}=0.83$. The last
value compares surprising well with those for the solar wind measured by the
Ulysses spacecraft (Pagel & Balogh, 2002; Nicol et al., 2008), and for the
magnetospheric cusp measured by the Polar satellite (Yordanova et al., 2004),
even though the frequencies differ by several orders of magnitude.

##  4 Summary

In conclusion, the wavelet spectral analysis of sunspot records has revealed
two different behaviors, depending on whether the Sun is quiescent or active.
This suggests two different kinds of turbulence in the convection zone,
controlled either by inertia or by rotation. The signature of such fully
developed turbulence is confirmed by the calculation of the GSF scaling
exponents, which indicate strong intermittency.

## Acknowledgments

We are grateful for support from the Dynamo Program at KITP (supported in part
by the National Science Foundation under Grant No. PHY05-51164), during which
this work was started. We thank Prof. Steve Tobias for helpful comments, and
Prof. Ilya Usoskin for providing the 10Be data. Finally F.P. and R.S. are
grateful for support from a RFBR/CNRS 07-01-92160 PICS grant.

## References

  * Abreu et al. (2008) Abreu J. A., Beer J., Steinhilber F., Tobias S. M., Weiss N. O., 2008, Geophys. Research Lett., 35, 20109
  * Alemany et al. (2000) Alemany A., Marty P., Plunian F., Soto J., 2000, Journal of Fluid Mechanics, 403, 263
  * Beer et al. (1998) Beer J., Tobias S., Weiss N., 1998, Solar Physics, 181, 237
  * Browning et al. (2006) Browning M. K., Miesch M. S., Brun A. S., Toomre J., 2006, ApJ Lett., 648, L157
  * Carbone (1993) Carbone V., 1993, Physical Review Letters, 71, 1546
  * Constable & Johnson (2005) Constable C., Johnson C., 2005, Physics of the Earth and Planetary Interiors, 153, 61
  * Courtillot & Le Mouel (1988) Courtillot V., Le Mouel J. L., 1988, Annual Review of Earth and Planetary Sciences, 16, 389
  * Frick et al. (1997) Frick P., Galyagin D., Hoyt D. V., Nesme-Ribes E., Schatten K. H., Sokoloff D., Zakharov V., 1997, A&A, 328, 670
  * Lawrence et al. (1995) Lawrence J. K., Cadavid A. C., Ruzmaikin A. A., 1995, ApJ, 455, 366
  * Meneveau & Sreenivasan (1987) Meneveau C., Sreenivasan K. R., 1987, Physical Review Letters, 59, 1424
  * Morfill et al. (1991) Morfill G. E., Scheingraber H., Voges W., Sonett C. P., 1991, in Sonett C. P., Giampapa M. S., Matthews M. S., eds, The Sun in Time Sunspot number variations - Stochastic or chaotic. pp 30–58
  * Nicol et al. (2008) Nicol R. M., Chapman S. C., Dendy R. O., 2008, ApJ, 679, 862
  * Pagel & Balogh (2002) Pagel C., Balogh A., 2002, Journal of Geophysical Research (Space Physics), 107, 1178
  * Sakuraba & Hamano (2007) Sakuraba A., Hamano Y., 2007, Geophys. Research Lett., 34, 15308
  * Solanki et al. (2004) Solanki S. K., Usoskin I. G., Kromer B., Schüssler M., Beer J., 2004, Nat, 431, 1084
  * Tobias et al. (2001) Tobias S. M., Brummell N. H., Clune T. L., Toomre J., 2001, ApJ, 549, 1183
  * Usoskin et al. (2003) Usoskin I. G., Solanki S. K., Schüssler M., Mursula K., Alanko K., 2003, Physical Review Letters, 91, 211101
  * Yordanova et al. (2004) Yordanova E., Grzesiak M., Wernik A., Popielawska B., Stasiewicz K., 2004, Annales Geophysicae, 22, 2431
  * Zhou (1995) Zhou Y., 1995, Physics of Fluids, 7, 2092


# An Entropy Weighted Nonnegative Matrix Factorization Algorithm for Feature
Representation

Jiao Wei  Can Tong  Bingxue Wu  Qiang He  Shouliang Qi  Yudong Yao  Fellow,
IEEE  and Yueyang Teng* This work was supported by the Fundamental Research
Funds for the Central Universities of China (N2019006 and N180719020).
J. Wei, C. Tong, B. Wu, Q He and S. Qi are with the College of Medicine and
Biological Information Engineering, Northeastern University, Shenyang 110169,
China.
Y. Yao is with the Department of Electrical and Computer Engineering, Stevens
Institute of Technology, Hoboken, NJ 07030, USA.
Y. Teng is with the College of Medicine and Biological Information
Engineering, Northeastern University, Shenyang 110169, China, and also with
the Key Laboratory of Intelligent Computing in Medical Image, Ministry of
Education, Shenyang 110169, China (email: tengyy@bmie.neu.edu.cn).

###### Abstract

Nonnegative matrix factorization (NMF) has been widely used to learn low-
dimensional representations of data. However, NMF pays the same attention to
all attributes of a data point, which inevitably leads to inaccurate
representation. For example, in a human-face data set, if an image contains a
hat on the head, the hat should be removed or the importance of its
corresponding attributes should be decreased during matrix factorizing. This
paper proposes a new type of NMF called entropy weighted NMF (EWNMF), which
uses an optimizable weight for each attribute of each data point to emphasize
their importance. This process is achieved by adding an entropy regularizer to
the cost function and then using the Lagrange multiplier method to solve the
problem. Experimental results with several data sets demonstrate the
feasibility and effectiveness of the proposed method. We make our code
available at https://github.com/Poisson-EM/Entropy-weighted-NMF.

clustering, entropy regularizer, low-dimensional representation, nonnegative
matrix factorization (NMF).

##  1 Introduction

With the rapid development of data acquisition technology, large amounts of
data, such as online documents, medical images, various video series, traffic
data, health data and other high-dimensional data, are accumulating.
Therefore, dimensionality reduction [1] has become an essential step in data
mining. Vector quantization (VQ) [2], singular value decomposition (SVD) [3],
principal component analysis (PCA) [4], independent component analysis (ICA)
[5] and concept factorization (CF) [6], [7] are some most commonly used
dimensionality reduction methods. However, due to negative components, these
methods typically cannot be reasonably explained in some practical problems.
Therefore, developing a nonnegative factorization method is valuable for
research. Thus, researchers have investigated the nonnegative matrix
factorization (NMF) proposed by Lee and Seung [8]. This method represents a
nonnegative matrix by a product of two low rank nonnegative matrices and then
learns a parts-based representation. In recent years, NMF has been applied in
many fields, including cluesting, dimensionality reduction, blind source
separation, etc. [9, 10, 11, 12].
In recent years, many variants of NMF have been proposed to extend its
applicable range. For example, Ding _et al._ [13] learned new low-dimensional
features from data with convenient clustering interpretation using Semi-NMF,
which allows the data matrix and the base matrix to have mixed signs [14, 15].
They also developed Convex-NMF by restricting the base vectors to a convex
combination of the data points. Considering that the orthogonality constraint
leads to sparsity, Pompili _et al._ [16] proposed an orthogonal NMF (ONMF)
method, which adds orthogonality constraints to the base and representation
matrices. There are two ways to solve ONMF: the Lagrange method [16] and the
natural gradient on the Stiefel manifold [17].

Blondel _et al._ [18] conducted an extended study for NMF by incorporating
predetermined weights to each attribute of each data point, demonstrating that
weights can produce important flexibility by better emphasizing certain
features in image approximation problems. Other researchers then conducted a
series of studies with the weighted NMF (WNMF). For example, Kim and Choi [19]
proposed a new WNMF method to process an incomplete data matrix with missing
entries, which combined the binary weights into the NMF multiplication update.
Li and Wu [20] proposed a weighted nonnegative matrix tri-factorization method
for co-clustering that weights each row and column of the original matrix in a
specific way, and the normalized cut information is combined into the
optimization model. Then, to assign labels to images, Kalayeh _et al._ [21]
proposed a weighted expansion method of multiview NMF, which imposes a
consistent constraint on the representation matrices between different
features, and a weight matrix mitigates data set imbalance. Dai _et al._ [22]
proposed applying the WNMF method to image recovery, and experimental results
showed that, particularly for data affected by salt and pepper noise, the
method could remove noise effectively and could also provide a more accurate
subspace representation.

We refer these methods as “hard WNMF”; however, one primary drawback is that
the weights, which they rely on, must be predetermined. Variable weights for
NMF are challenging to implement that can be resolved by designing
interpretable and computable weights.

This paper presents an entropy weighted NMF (EWNMF) method that assigns a
weight to indicate the importance of each attribute of each data point in
matrix factorizing. Then, the entropy of these weights is used to regularize
the cost function for obtaining an easy computable solution, which makes the
range of weights fall within [0, 1] with a summation of one; thus, the weights
can be explained as the probability of the contribution of an attribute of one
data point to NMF. Experimental results with several real data sets show that
EWNMF performs better than other NMF variants.

##  2 Methodology

###  2.1 Related research

NMF is a matrix factorization method that focuses on data matrices with
nonnegative elements, which can reveal hidden structures and patterns from
generally redundant data. We will review the standard NMF as follows.
$\mathbf{Notations}~{}$ In this paper, matrices are denoted as capital
letters. For a matrix $A$, $A_{*i}$, $A_{i*}$ and $A_{ij}$ denote the $i-th$
column, the $i-th$ row and $(i,j)-th$ element of $A$, respectively; the
Frobenius norm is represented as $\left\|A\right\|_{F}$; $\odot$ and $./$ mean
the item-by-item multiplication and division of two matrices, respectively;
$A^{T}$ denotes the transpose of $A$; $A\geq 0$ means that all elements of $A$
are equal to or larger than 0.
The expression of NMF is:

$X\approx WH$ (1)

where the matrix $X\subseteq R^{M\times N}$ denotes the given nonnegative
matrix in which each column is a data point. The goal of NMF is to find two
low-dimensional nonnegative matrices: $W\subseteq R^{M\times K}$ is called the
base matrix, and $H\subseteq R^{K\times N}$ is called the representation
matrix, whose product can approximate the original matrix [23, 24], where
$K<<min\\{M,N\\}$.
There are different standards to measure the quality of decomposition. Lee and
Seung proposed using the square of the Euclidean distance and the Kullback-
Leibler divergence. In this paper, the Euclidean distance is used, and the
formula is expressed as:

$\displaystyle\min$ $\displaystyle~{}F_{1}(W,H)=\left\|X-WH\right\|_{F}^{2}$
(2) $\displaystyle s.~{}t.$ $\displaystyle~{}W\geq 0,~{}H\geq 0$

To alternatively minimize $W$ and $H$ in Eq. (2), the construction of the
auxiliary function is important to determine the iterative update rule.
$\mathbf{Definition~{}1}~{}(Auxiliary~{}function)$ If the function
$G(h,h^{\prime})$ satisfies the following conditions:

$\displaystyle G(h,h^{\prime})\geq F(h)$ (3) $\displaystyle
G(h^{\prime},h^{\prime})=F(h^{\prime})$

where $h^{\prime}$ is a given value; then, $G(h,h^{\prime})$ is the auxiliary
function of $F(h)$ on $h^{\prime}$.

Then, we can draw the following conclusion.
$\mathbf{Lemma~{}1.}~{}$ If $G(h,h^{\prime})$ is an auxiliary function of
$F(h)$, then under the update rule:

$h^{*}=arg\mathop{\min}_{h}G(h,h^{\prime})$ (4)

the function $F(h)$ does not increase.
$\mathbf{Proof.~{}}$ The conditions satisfied by the auxiliary function make
this proof marked because:

$F(h^{*})\leq G(h^{*},h^{\prime})\leq G(h^{\prime},h^{\prime})\leq
F(h^{\prime})$ (5)

Thus, if the auxiliary function reaches the minimum, the original function
should decrease.
Then, we construct the update rule for the NMF problem. We consider $W$ first,
where $W^{t}>0$ and $H>0$ are given. Let
$\xi_{ijk}=W^{t}_{ik}H_{kj}/(W^{t}H)_{ij}$, of course, $\xi_{ijk}\geq 0$ and
$\sum_{k=1}^{K}\xi_{ijk}=1$. Therefore, the auxiliary function of standard NMF
is:

$f_{1}(W,W^{t})=\sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{K}\xi_{ijk}(X_{ij}-\frac{W_{ik}H_{kj}}{\xi_{ijk}})^{2}$
(6)

Because the function is separable, it can be easily minimized. We thus take
the partial derivative of Eq. (6) and set it to zero so that we can obtain the
following update rule as:

$W\leftarrow W\odot(XH^{T})~{}./~{}(WHH^{T})$ (7)

The method of constructing the auxiliary function of $H$ is similar to that of
$W$. Then, the following update of $H$ is obtained:

$H\leftarrow H\odot(W^{T}X)~{}./~{}(W^{T}WH)$ (8)

Figure 1: Diagram of the proposed method: a hat in the first image destroys
the accuracy of feature representation, which should be assigned no or little
importance.

###  2.2 Proposed method

Different from previous methods, an optimizable weight matrix is used to
measure the importance of the attributes in matrix factorizing. Fig. 1 shows
the explanation of the proposed idea that, if there is a hat in the only
image, the corresponding attributes of the hat must destroy the result of NMF,
and we should certainly eliminate or weaken their importances. Thus, it may be
necessary to provide different importances to the attributes of each data
point based on the following constraint, which will lead to a new optimization
model:

$\displaystyle\min$
$\displaystyle~{}F_{2}(W,H,T)=\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}[X_{ij}-(WH)_{ij}]^{2}$
(9) $\displaystyle s.~{}t.$ $\displaystyle~{}W\geq 0,~{}H\geq 0,~{}T\geq
0,~{}\sum_{i=1}^{M}T_{ij}=1$

We first consider the new variable $T$, which can be solved in an alternative
optimization manner. For fixed $W$ and $H$, $T_{ij}$ is very easy to solve as
$T_{ij}=1$ if $E_{ij}=min\\{|E_{1j}|,|E_{2j}|,...,|E_{Mj}|\\}$ or 0
otherwise111The update rule can be easily explained by an example as:
$\displaystyle\min~{}\\{3,1,2\\}=min$ $\displaystyle~{}3T_{1}+1T_{2}+2T_{3}$
$\displaystyle s.~{}t.$ $\displaystyle~{}T_{1}\geq 0,T_{2}\geq 0,T_{3}\geq 0$
$\displaystyle~{}T_{1}+T_{2}+T_{3}=1$ The solution is that $T_{1}=0$,
$T_{2}=1$ and $T_{3}=0$, in which $T_{2}$ corresponds to the minimum value of
{3, 1, 2}. This process is similar to the computation of the weights in the
k-means algorithm., where $E=X-WH$. And it demonstrates the simple fact that
only one element of $T_{*j}$ is 1, and the others are 0, which is incompatible
with the real problem.
To address this issue, we apply an entropy regularizer to penalize the cost
function of NMF to obtain a weight in the range of [0, 1] instead of 0 or 1,
which uses information entropy to calculate the uncertainty of weights. The
new optimization problem is rewritten as follows:

$\displaystyle\min$
$\displaystyle~{}F_{3}(W,H,T)=\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}[X_{ij}-(WH)_{ij}]^{2}$
(10)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\gamma\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}ln(T_{ij})$
$\displaystyle s.~{}t.$ $\displaystyle~{}W\geq 0,~{}H\geq 0,~{}T\geq
0,~{}\sum_{i=1}^{M}T_{ij}=1$

where $\gamma\geq 0$ is a given hyperparameter. The first term in Eq. (10) is
the sum of errors, and the second term is the negative entropy of the weights.
The original cost function in Eq. (9) results in only one attribute of each
data point being involved in feature representation, and the entropy
regularizer will stimulate more attributes to help feature representation.
This equation can be solved by a simple algorithm, which is based on the
following proposition:
$\mathbf{Proposition~{}1.~{}}$Given the matrices $W$ and $H$, $T_{ij}$ in Eq.
(10) is minimized when:

$T_{ij}=\frac{e^{-\frac{[X_{ij}-(WH)_{ij}]^{2}}{\gamma}}}{\sum_{l=1}^{M}e^{-\frac{[X_{lj}-(WH)_{lj}]^{2}}{\gamma}}}$
(11)

$\mathbf{Proof.~{}}$ We construct the Lagrange function of Eq. (10) with
respective to $T$ as:

$\displaystyle L(T,\lambda)=$
$\displaystyle\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}[X_{ij}-(WH)_{ij}]^{2}$ (12)
$\displaystyle+\gamma\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}ln(T_{ij})-\sum_{j=1}^{N}\lambda_{j}(\sum_{i=1}^{M}T_{ij}-1)$

where $[\lambda_{1},\lambda_{2},...,\lambda_{N}]$ is a vector containing the
Lagrange multipliers corresponding to the constraints.
By setting the gradient of Eq. (12) with respect to $\lambda_{j}$ and $T_{ij}$
to zero, we obtain the following equation system:

$\displaystyle\frac{\partial L}{\partial\lambda_{j}}$
=$\sum_{i=1}^{M}T_{ij}-1=0$ (13) $\displaystyle\frac{\partial L}{\partial
T_{ij}}$ =$[X_{ij}-(WH)_{ij}]^{2}+\gamma lnT_{ij}+\gamma-\lambda_{j}=0$ (14)

From Eq. (14), we know that:

$\displaystyle
T_{ij}=e^{\frac{\lambda_{j}-\gamma}{\gamma}}e^{-\frac{[X_{ij}-(WH)_{ij}]^{2}}{\gamma}}$
(15)

Substituting Eq. (15) into Eq. (13), we have:

$\displaystyle\sum_{i=1}^{M}T_{ij}=e^{\frac{\lambda_{j}-\gamma}{\gamma}}\sum_{i=1}^{M}e^{{-\frac{[X_{ij}-(WH)_{ij}]^{2}}{\gamma}}}=1$
(16)

It follows that:

$\displaystyle
e^{\frac{\lambda_{j}-\gamma}{\gamma}}=\frac{1}{\sum_{l=1}^{M}e^{{-\frac{[X_{lj}-(WH)_{lj}]^{2}}{\gamma}}}}$
(17)

Substituting this expression to Eq. (15), we find that:

$T_{ij}=\frac{e^{-\frac{[X_{ij}-(WH)_{ij}]^{2}}{\gamma}}}{\sum_{l=1}^{M}e^{-\frac{[X_{lj}-(WH)_{lj}]^{2}}{\gamma}}}$
(18)

Then, we can solve $W$ and $H$ with fixed $T$, which is similar to the
standard NMF method. For example, we can construct the following auxiliary
function about $W$:

$\displaystyle f_{3}(W,W^{t})=$
$\displaystyle\sum_{i=1}^{M}\sum_{j=1}^{N}\sum_{k=1}^{K}T_{ij}\xi_{ijk}(X_{ij}-\frac{W_{ik}H_{kj}}{\xi_{ijk}})^{2}$
(19) $\displaystyle+\gamma\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}lnT_{ij}$

Setting the partial derivative of $f_{3}(W,W^{t})$ to zero yields the
following update rule:

$W\leftarrow W\odot(T\odot X)H^{T}~{}./~{}\\{[T\odot(WH)]H^{T}\\}$ (20)

Similarly, we can also easily obtain the update rule for $H$ as follows:

$H\leftarrow H\odot W^{T}(T\odot X)~{}./~{}\\{W^{T}[T\odot(WH)]\\}$ (21)

The update rules to $W$ and $H$ are similar to the existing WNMF methods in
[21], [22]. Optimizing EWNMF is summarized as follows in
$\mathbf{Algorithm~{}1}$:

Algorithm 1 EWNMF

1:Given the input nonnegative matrix $X\subseteq R^{M\times N}$, the number
_K_ of reduced dimensions and hyperparameter $\gamma$;

2:the weight matrix $T$, the base matrix $W$ and the representation matrix
$H$;

3:Randomly initialize $W\subseteq R^{M\times K}\textgreater 0$ and $H\subseteq
R^{K\times N}\textgreater 0$;

4:while not convergence do

5: Update $T$ by Eq. (18);

6: Update $W$ by Eq. (20);

7: Update $H$ by Eq. (21);

8:end while

9:return $T$, $W$ and $H$.

###  2.3 Extensions

The proposed entropy weighted method can also be applied to standard NMF using
the square of the Euclidean distance and to those that use KL divergence,
$\alpha$-divergence and other existing NMFs, such as ONMF, Semi-NMF, Convex-
NMF, etc., which demonstrates its good compatibility.
For example, incorporating the entropy regularizer into KL-NMF, its cost
function is expressed as:

$\displaystyle F_{4}(T,W,H)=$
$\displaystyle\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}[X_{ij}log\frac{X_{ij}}{(WH)_{ij}}-X_{ij}$
(22)
$\displaystyle+(WH)_{ij}]+\gamma\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}ln(T_{ij})$

where $T$ can be derived from the above method, and then the update is shown
below:

$T_{ij}=\frac{e^{-\frac{X_{ij}log\frac{X_{ij}}{(WH)_{ij}}-X_{ij}+(WH)_{ij}}{\gamma}}}{\sum_{l=1}^{M}e^{-\frac{X_{lj}log\frac{X_{lj}}{(WH)_{lj}}-X_{lj}+(WH)_{lj}}{\gamma}}}$
(23)

Applying the entropy regularizer to the $\alpha$-divergence-based NMF method,
the cost function is:

$\displaystyle F_{5}(T,W,H)=$
$\displaystyle\frac{1}{\alpha(\alpha-1)}\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}[X_{ij}^{\alpha}(WH)_{ij}^{1-\alpha}-\alpha
X_{ij}$ (24)
$\displaystyle+(\alpha-1)(WH)_{ij}]+\gamma\sum_{i=1}^{M}\sum_{j=1}^{N}T_{ij}ln(T_{ij})$

where $\alpha\in R$ is a given value. The update method of $T$ is as follows:

$T_{ij}=\frac{e^{-\frac{[X_{ij}^{\alpha}(WH)_{ij}^{1-\alpha}-\alpha
X_{ij}+(\alpha-1)(WH)_{ij}]}{\alpha(\alpha-1)\gamma}}}{\sum_{l=1}^{M}e^{-\frac{[X_{lj}^{\alpha}(WH)_{lj}^{1-\alpha}-\alpha
X_{lj}+(\alpha-1)(WH)_{lj}]}{\alpha(\alpha-1)\gamma}}}$ (25)

The entropy weighted method can also be popularized to many other NMFs, which
is similar to the above derivations, and thus, we omit them.

Figure 2: Illustrative figures of EWNMF: (a) the source signals, (b) the mixed
signals, in which the first one is corrupted at the beginning part and (c) the
obtained weights corresponding to (b).

##  3 experiments

###  3.1 Experimental description

Experiments were performed on a HP Compaq PC with a 3.40-GHz Core i7-6700 CPU
and 16 GB memory, and all the methods were implemented in MATLAB. We compare
the performance of the proposed methods with NMF [8], ONMF [16], Semi-NMF [13]
and Convex-NMF [13] on five public data sets, including the Yale, UMISTface,
Caltech101, GTFD and TDT2 data sets.
The Yale face data set [25] was created by Yale University and contains data
describing 15 people. Each person has 11 face images with different
expressions, postures and lighting. The data set has a total of 150 images,
and each has a size of $32\times 32$.
The UMISTface data set [26] was established by the University of Manchester,
UK. The data set has a total of 1012 images, including 20 people, each with
different angles and different poses. In this experiment, the pixel number of
each image is $32\times 32$.
The Caltech101 data set [27] contains 9144 images split between 101 different
object categories, as well as an additional background/clutter category. This
data set has approximately 40 to 800 images per category, and most categories
have approximately 50 images. The size of each image is $32\times 32$.
The GTFD [28] contains 750 images taken in two different sessions and includes
50 people. All people in the data set are represented with cluttered
backgrounds. The images show faces with different facial expressions and
lighting conditions, and each has a size of $32\times 32$.
The TDT2 Audio Corpus [29] contains six sources, including two news special
lines (APW and NYT), two radio programs (VOA and PRI) and two TV programs (CNN
and ABC). In this experiment, only the largest 30 categories are used (a total
of 9394 documents).
Important details of these data sets are shown in Table 1.

Table 1: Statistics of the used standard data sets. Data set | Yale | UMISTface | Caltech101 | GTFD | TDT2
---|---|---|---|---|---
Points | 165 | 1012 | 9144 | 750 | 9394
Dimensions | 1024 | 1024 | 1024 | 1024 | 36771
Class | 15 | 20 | 101 | 50 | 30
Figure 3: Clustering performance versus the hyperparameter $\gamma$ on the
Yale data set: (a) ACC and (b) NMI. Figure 4: Clustering performance versus
the hyperparameter $\gamma$ on the UMISTface data set: (a) ACC and (b) NMI.

After obtaining a new feature representation, we use k-means to cluster them
and then compare it with the label to evaluate the clustering results.
Clustering accuracy (ACC) [30], [31] and normalized mutual information (NMI)
[32], [33] are used to evaluate the performance of these clustering results.
Given a set of the ground true class labels $y$ and the obtained cluster
labels $y^{\prime}$, the clustering accuracy is defined as:

$ACC=\frac{\sum_{i=1}^{N}\delta(y_{i},map(y^{\prime}_{i}))}{N}$ (26)

where:

$\delta(a,b)=\begin{cases}1&\text{\emph{a = b}}\\\
0&\text{otherwise}\end{cases}$

and $map(\cdot)$ is a permutation mapping function that maps the obtained
cluster labels to the real labels. The higher the ACC value is, the better the
clustering performance.
NMI is used to calculate the agreement between the two data distributions and
is defined as follows:

$NMI(y,y^{\prime})=\frac{MI(y,y^{\prime})}{max(H(y),H(y^{\prime}))}$ (27)

where $H(y)$ is the entropy of $y$. $MI(y,y^{\prime})$ quantifies the amount
of information between two random variables (i.e., $y$ and $y^{\prime}$) and
is defined as:

$MI(y,y^{\prime})=\sum_{y_{i}\in y,y^{\prime}_{j}\in
y^{\prime}}p(y_{i},y^{\prime}_{j})log(\frac{p(y_{i},y^{\prime}_{j})}{p(y_{i})p(y^{\prime}_{j})})$
(28)

where $p(y_{i})$ and $p(y^{\prime}_{j})$ are the probabilities that a data
point selected from the data set belongs to the clusters $y_{i}$ and
$y^{\prime}_{j}$, respectively; and $p(y_{i},y^{\prime}_{j})$ is the joint
probability that an arbitrarily selected data point belongs to clusters
$y_{i}$ and $y^{\prime}_{j}$ concurrently. The NMI score ranges from 0 to 1,
and the larger NMI is, the better the clustering performance.

Table 2: ACC and NMI on the Yale data set (%). Method | NMF | ONMF | Semi-NMF | Convex-NMF | EWNMF
---|---|---|---|---|---
ACC | 39.03 | 41.21 | 42.48 | 31.88 | 46.06
NMI | 44.74 | 47.36 | 48.48 | 37.91 | 53.27
Table 3: ACC and NMI on the UMISTface data set (%). Method | NMF | ONMF | Semi-NMF | Convex-NMF | EWNMF
---|---|---|---|---|---
ACC | 51.06 | 50.99 | 52.15 | 23.52 | 64.05
NMI | 69.63 | 68.34 | 67.70 | 29.13 | 76.01
Figure 5: Clustering performance versus cluster number on the Yale data set:
(a) AC and (b) NMI. Figure 6: Clustering performance versus cluster number on
the UMISTface data set: (a) AC and (b) NMI. Figure 7: Clustering performance
versus cluster number on the Caltech101 data set: (a) AC and (b) NMI.

Before the experiment, we normalized all the data sets to scale the minimum
and maximum values of each data point to 0 and 1, respectively. All the
methods use the same random distribution for initialization of $W$ and $H$ to
make them uniformly distributed on [0.1 1.1] and perform 300 iterations to
ensure sufficient convergence.

###  3.2 Experimental results

####  3.2.1 Signal unmixing on synthetic data

We use a synthetic data set to simulate a signal-mixture process, which mixes
two source signals by a uniformly distributed matrix on [0 1] and then
destroys the beginning part of the first one. The destroyed part in the mixed
signal will hinder signal recovery. Then, we use EWNMF to unmix the signals to
demonstrate the usefulness of the proposed entropy weighted method. Fig. 2
shows the source signals, mixed signals (including the destroyed signal) and
the obtained weights. The weights of the destroyed elements in the first
signal are small and can even be considered as zero, and the other weights are
similar. Thus, we can conclude that EWNMF can provide correct weights to the
importance of attributes in the data set.

####  3.2.2 EWNMF compared with the standard NMF

We investigate the ability of the entropy weighted strategy to improve the
performance of the standard NMF. We apply NMF and EWNMF to reduce the
dimension number of the Yale and UMISTface data sets and then use k-means to
cluster these new representations. Note that the number of reduced dimensions
is equal to that of the clusters. Figs. 3 and 4 show the clustering results
evaluated by ACC and NMI. EWNMF indeed provides better performance than the
standard NMF and can achieve a consistently superior performance to the
standard NMF on a wide range of the hyperparameter $\gamma$, demonstrating the
robustness of EWNMF.
The optimal clustering results on the entire Yale and UMISTface data sets are
shown in Tabs. 2 and 3. An interesting observation is highlighted and shows a
consistent result with the above figures that EWNMF still provides the best
evaluation standards.

Figure 8: Clustering performance versus cluster number on the GTFD data set:
(a) AC and (b) NMI. Figure 9: Clustering performance versus cluster number on
the TDT2 data set: (a) AC and (b) NMI.

Figure 10: Cost function on the Yale data set. From left to right, they are
obtained by (a) NMF, (b) ONMF, (c) Semi-NMF, (d) Convex-NMF and (e) EWNMF.

####  3.2.3 Clustering results with different cluster numbers

We study the relationship between the evaluation standards and cluster number.
The hyperparameter $\gamma$ is selected in $\\{10^{i},i=-8,-7,...,7,8\\}$ to
obtain the optimal results in a large range. Because the NMF problem does not
have a sole solution, we randomly initialized 10 times to obtain a credible
averaged ACC and NMI. Different cluster numbers ranging from 2 to 10 are
selected. In a certain data set with _k_ clusters, the experimental details
are described as follows:
1) Randomly select _k_ categories as a subset for the following experiment.
2) Randomly initialize $W$ and $H$, obtain new representations, and cluster
them by k-means. Note that EWNMF uses the selected hyperparameter $\gamma$
according to the above instruction.
3) Repeat 1) and 2) 10 times to obtain an average result.

Figure 11: Some base images with respect to the Yale data set, in which the
summation of every column is normalized to one. From left to right, they are
obtained by (a) NMF, (b) ONMF, (c) Semi-NMF, (d) Convex-NMF and (e) EWNMF.

The clustering results, ACC and NMI versus the number of clusters, are
reported in Figs. 5-9. The proposed method can generally yield more accurate
clustering results, and only for two clusters of the GTFD data set, Semi-NMF
occasionally exceeds the accuracy of the proposed method. Also, none of the
other methods performs better in all aspects than EWNMF.

####  3.2.4 Convergence speed and part-based learning

From the theoretical analysis described above, we can conclude that the
proposed method is monotonically decreasing, but due to the nonconvexity of
the cost function, it cannot be guaranteed to be strictly convergent [34].
Thus, we investigate the convergence speed of the NMF methods.

Fig. 10 shows the cost functions of NMF, ONMF, Semi-NMF, Convex-NMF and EWNMF
on the Yale data set. Except for the Convex-NMF method, which converges
slowly, the other methods reached a stable point within 200 iterations,
demonstrating that they have similar convergence speeds. However, ONMF cannot
guarantee the monotonic decrease of the cost function because it uses a
proximal Lagrange multiplier method. Fig. 11 shows the base images on the Yale
data set and demonstrates that Semi-NMF and Convex-NMF identify more global
faces. EWNMF has poor locality compared to ONMF; however, the former has a
better clustering effect than the latter. It is visually difficult to evaluate
the performance of the NMF methods, even though they are all markedly
different.

##  4 Conclusion

This paper proposes a new NMF method, which adds weights to each attribute of
each data to emphasize their importance. We introduce an entropy regularizer
to consider these weights as the probabilities of importance within [0 1],
which mimics the process of human reasoning more accurately. These weights can
be solved by the Lagrange multiplier method, and a simple update is achieved.
The experimental results show that the proposed method produces performance
that is competitive with those of existing methods.
The entropy regularizer requires an additional hyperparameter to control the
certainty of the weights. In the future, we plan to develop an auto-adjustment
strategy for this hyperparameter.

## Acknowledgment

The authors would like to thank the editor and anonymous reviewers for their
constructive comments and suggestions.

## References

  * [1] I. K. Fodor, “A survey of dimensionality reduction techniques,” in _Proc. Center Appl. Sci. Comput., Lawrence Livermore Nat. Lab._ , 2002, pp. 1-24.
  * [2] A. Gersho and R. M. Gray, _Vector Quantization and Signal Compression._ Springer Science and Business Media, 1992.
  * [3] C. Hu, X. Lu, M. Ye, and W. Zeng, “Singular value decomposition and local near neighbors for face recognition under varying illumination,” _Pattern Recognit._ , vol. 64, pp. 60-83, 2017.
  * [4] I. T. Jolliffe, _Principal Component Analysis._ Springer, 2002.
  * [5] A. Hyvarinen, E. Oja, “Independent component analysis: algorithms and applications,” _Neural Netw._ vol. 13, no. 4-5, pp. 411-430, 2000.
  * [6] S. Peng, W. Ser, B. Chen, L. Sun, and Z. Lin, “Correntropy based graph regularized concept factorization for clustering,” _Neurocomput._ , vol. 316, pp. 34-48, 2018.
  * [7] W. Yan, B. Zhang, S. Ma, and Z. Yang, “A novel regularized concept factorization for document clustering,” _Knowl. Based Syst._ , vol. 135, pp. 147-158, 2017.
  * [8] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” _Nature_ , vol. 401, no. 6755, pp. 788-791, 1999.
  * [9] F. Shahnaz, M. W. Berry, V. P. Pauca, and R. J. Plemmons, “Document clustering using nonnegative matrix factorization,” _Inf. Manag._ , vol. 42, no. 2, pp. 373-386, 2006.
  * [10] V. P. Pauca, J. Piper, and R. J. Plemmons, “Nonnegative matrix factorization for spectral data analysis,” _Linear Algebra Appl._ , vol. 416, no. 1, pp. 29-47, 2006.
  * [11] S. Jia and Y. Qian, “Constrained nonnegative matrix factorization for hyperspectral unmixing,” _IEEE Trans. Geosci. Remote Sens._ , vol. 47, no. 1, pp. 161-173, 2009.
  * [12] H. Kim and H. Park, “Sparse non-negative matrix factorizations via alternating non-negativity constrained least squares for microarray data analysis,” _Bioinform._ , vol. 23, no. 12, pp. 1495-1502, 2007.
  * [13] C. H. Q. Ding, T. Li, and M. I. Jordan, “Convex and semi-nonnegative matrix factorizations,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 32, no. 1, pp. 45-55, 2010.
  * [14] A. Suleman, “A convex semi-nonnegative matrix factorization approach to fuzzy k-means clustering,” _Fuzzy Sets Syst._ , vol. 270, pp. 90-110, 2015.
  * [15] R. Hamon, V. Emiya, and C. Fevotte, “Convex nonnegative matrix factorization with missing data,” in _Proc. IEEE 26th Int. Workshop Mach. Learn. Signal Process._ , 2016, pp. 1-6.
  * [16] C. Ding, T. Li, W. Peng, and P. Haesun, “Orthogonal nonnegative matrix t-factorizations for clustering,” in _Proc. 12th ACM SIGKDD Int. Conf. Knowl. Discov. Data Mining_ , 2006, pp. 126-135.
  * [17] S. Choi, “Algorithms for orthogonal nonnegative matrix factorization,” in _Proc. Int. Jt. Conf. Neural Netw._ , 2008, pp. 1828-1832.
  * [18] V. D. Blondel, N.-D. Ho, and P. V. Dooren, “Weighted nonnegative matrix factorization and face feature extraction,” _Image. Vis. Comput._ , 2007.
  * [19] Y. Kim and S. Choi, “Weighted nonnegative matrix factorization,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Proces_ , 2009, pp. 1541-1544.
  * [20] Z. Li and X. Wu, “Weighted nonnegative matrix tri-factorization for co-clustering,” in _Proc. IEEE Int. Conf. Tools Arti. Intell._ , 2011, pp. 811-816.
  * [21] M. M. Kalayeh, H. Idrees, and M. Shah, “NMF-KNN: image annotation using weighted multi-view non-negative matrix factorization,” in _Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2014, pp. 184-191.
  * [22] X. Dai, K. Zhang, J. Xiong, X. Zhang, Z. Tu, and N. Zhang, “Weighted non-negative matrix factorization for image recovery and representation,” in _Int. Conf. Inf. Sci. Tech. (ICIST)_ , 2020, pp. 288-293.
  * [23] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in _Proc. Int. Conf. Neural Inform. Process. Syst._ , 2001, pp. 556-562.
  * [24] P. Paatero and U. Tapper, “Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values,” _Environ._ , vol. 5, no. 2, pp. 111-126, 1994.
  * [25] P. N. Bellhumer, J. P. Hespanha, and D. J. Kriegman, “Eigenfaces vs. fisherfaces: recognition using class specific linear projection,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 19, no. 7, pp. 711-720, 1997.
  * [26] Sheffield/UMIST face data set. [Online]. Available: http://www.sheffield.ac.uk/eee/research/iel/research/face/.
  * [27] F. Li, R. Fergus, and P. Perona, “Learning generative visual models from few training examples: an incremental Bayesian approach tested on 101 object categories,” in _Comput. Vis. Image Underst._ , 2004, pp. 178-178.
  * [28] Georgia tech face data set. [Online]. Available: Ftp://ftp.ee.gatech.edu/pub/users/hayes/facedb/.
  * [29] D. Graff, TDT2 English audio LDC99s84, DVD, Linguistic Data Consortium, Philadelphia, 1999.
  * [30] H. Liu, Z. Wu, X. Li, D. Cai, and T. S. Huang, “Constrained nonnegative matrix factorization for image representation,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 34, no. 7, pp. 1299-1311, 2012.
  * [31] L. Lovasz and M. Plummer, _Matching Theory._ Akademiai Kiado, 1986.
  * [32] D. Cai, X. He, and J. Han, “Document clustering using locality preserving indexing,” _IEEE Trans. Knowl. Data Eng._ , vol. 17, no. 12, pp. 1624-1637, 2005.
  * [33] F. Shahnaz, M. W. Berry, V. Pauca, and R. J. Plemmons, “Document clustering using nonnegative matrix factorization,” _Inf. Process. Manag._ , vol. 42, no. 2, pp. 373-386, 2004.
  * [34] C. J. Lin, “On the convergence of multiplicative update algorithms for nonnegative matrix factorization,” _IEEE Trans. Neural Netw._ , vol. 18, no. 6, pp. 1589-1596, 2007.

[]Jiao Wei is pursuing her M.S. degree in the School of Medical and
Bioinformatics Engineering at Northeastern University and received her B.S.
degree from Qufu Normal University in 2014. Her research interests include
machine learning and nonnegative matrix decomposition algorithms. []Tong Can
is a Ph.D. student and received a bachelor’s degree in mathematics and applied
mathematics from Northeast University in 2014 and a master’s degree in
computational mathematics from Northeastern University in 2018. Her primary
research interests include basic algorithm theory and acceleration methods of
machine learning. []Bingxue Wu is pursuing her M.S. degree in the School of
Medical and Bioinformatics Engineering at Northeastern University and received
her B.S. degree from Harbin Business University in 2015. Her research
interests are machine learning and compressive perception in medical imaging.
[]Qiang He received the Ph.D. degree in computer application technology from
the Northeastern University, Shenyang, China, in 2020. He is currently an
Associate Professor with the College of Medicine and Biological Information
Engineering, Northeastern University, Shenyang, China. His research interests
include social network analytic, machine learning, etc. []Shouliang Qi is an
associate professor at Northeastern University, China. After he received a
Philosophy degree doctor from Shanghai Jiao Tong University in 2007, he joined
the GE Global Research Center and was responsible for designing an innovative
magnetic resonance imaging (MRI) system. In 2014-2015, he worked as a visiting
scholar at Eindhoven University of Technology and Epilepsy Center
Kempenhaeghe, the Netherlands. In recent years, he has been conducting
productive studies in advanced MRI technology, intelligent medical imaging
computing and modeling, and the convergence of Nano-Bio-Info-Cogn (NBIC). He
has published more than 60 papers in peer-reviewed journals and international
conferences. He has won many academic awards, such as Chinese Excellent PH. D.
Dissertation Nomination Award, and the Award for Outstanding Achievement in
Scientific Research from Ministry of Education.

[]Yudong Yao (S’88-M’88-SM’94-F’11) received B.Eng. and M.Eng. degrees from
Nanjing University of Posts and Telecommunications, Nanjing, in 1982 and 1985,
respectively, and a Ph.D. degree from Southeast University, Nanjing, in 1988,
all in electrical engineering. He was a visiting student at Carleton
University, Ottawa, in 1987 and 1988. Dr. Yao has been with Stevens Institute
of Technology, Hoboken, New Jersey, since 2000 and is currently a professor
and department director of electrical and computer engineering. He is also a
director of the Stevens’ Wireless Information Systems Engineering Laboratory
(WISELAB). Previously, from 1989 to 2000, Dr. Yao worked for Carleton
University, Spar Aerospace Ltd., Montreal, and Qualcomm Inc., San Diego. He
has been active in a nonprofit organization, WOCC, Inc., which promotes
wireless and optical communications research and technical exchange. He served
as WOCC president (2008-2010) and chairman of the board of trustees
(2010-2012). His research interests include wireless communications and
networking, cognitive radio, machine learning, and big data analytics. He
holds one Chinese patent and thirteen U.S. patents. Dr. Yao was an Associate
Editor of IEEE Communications Letters (2000-2008) and IEEE Transactions on
Vehicular Technology (2001-2006) and an Editor for IEEE Transactions on
Wireless Communications (2001-2005). He was elected an IEEE Fellow in 2011 for
his contributions to wireless communications systems and was an IEEE ComSoc
Distinguished Lecturer (2015-2018). In 2015, he was elected a Fellow of
National Academy of Inventors. []Yueyang Teng received his Bachelor and Master
degrees from the Department of Applied Mathematics, Liaoning Normal University
and Dalian University of Technology, China, in 2002 and 2005, respectively.
From 2005 to 2013, he was a software engineer at Neusoft Positiron Medical
Systems Co., Ltd. In 2013, he received his doctorate degree in Computer
Software and Theory, Northeastern University. Since 2013, he has been a
lecturer at the Sino-Dutch Biomedical and Information Engineering School,
Northeastern University. His research interests include image processing and
machine learning.