# Universality classes of the Anderson transition in three-dimensional
symmetry classes AIII, BDI, C, D and CI

Tong Wang  International Center for Quantum Materials, School of Physics,
Peking University, Beijing 100871, China  Collaborative Innovation Center of
Quantum Matter, Beijing 100871, China  Tomi Ohtsuki  Physics Division, Sophia
University, Chiyoda-ku, Tokyo 102-8554, Japan  Ryuichi Shindou
<EMAIL_ADDRESS>International Center for Quantum Materials, School of
Physics, Peking University, Beijing 100871, China  Collaborative Innovation
Center of Quantum Matter, Beijing 100871, China

###### Abstract

We clarify universal critical properties of delocalization-localization
transitions in three-dimensional (3D) unitary and orthogonal classes with
particle-hole and/or chiral symmetries (classes AIII, BDI, D, C and CI). We
first introduce tight-binding models on cubic lattice that belong to these
five nonstandard symmetry classes respectively. Unlike the Bogoliubov-de
Gennes Hamiltonian for superconductors, all the five models have finite areas
of Fermi surfaces in the momentum space in the clean limit. Thereby, the
scaling theory of the Anderson transition guarantees the presence of the
delocalization-localization transitions at finite disorder strength in these
models. Based on this expectation, we carry out extensive transfer matrix
calculations of the Lyapunov exponents for zero-energy eigenstates of the
disordered tight-binding models with quasi-one-dimensional geometry. Near the
Anderson transition point, the correlation length diverges with a universal
critical exponent $\nu$. Using finite-size scaling analysis of the
localization length, we determine the critical exponent of the transitions and
scaling function of the (normalized) localization length in the three non-
standard unitary symmetry classes: $\nu_{\text{AIII}}=1.06\pm 0.02$,
$\nu_{\text{D}}=0.87\pm 0.03$, and $\nu_{\text{C}}=0.996\pm 0.012$. Our result
of the class C is consistent with a previous study of classical network model
[M. Ortuño et al, Phys. Rev. Lett. 102, 070603 (2009)]. The critical exponents
of the two non-standard orthogonal classes are estimated as
$\nu_{\text{CI}}=1.17\pm 0.02$ and $\nu_{\text{BDI}}=1.12\pm 0.06$. Our result
of the class CI is consistent with another previous study of a lattice model,
while the exponent of the class BDI is at variance with the previous
evaluation of nodal Dirac ring model [X. L. Luo et al, Phys. Rev. B. 101,
020202(R) (2020)].

##  I Introduction

The quantum interference of two counter-propagating waves form a standing wave
that does not move in space, suppressing particle diffusions completely in
random media. After Anderson’s seminal proposal of the localization of
electron wavefunctions in disordered solids [1, 2, 3], theoretical
understandings of the localization phenomena have been elaborated by scaling
theories [4, 5], field theories [6, 7, 8, 9, 10], numerical simulations [11,
12, 13, 14], and symmetry classifications of random matrices [15, 16, 17, 18].
It is widely acknowledged that delocalization-localization (Anderson)
transition occurs in a variety of physical systems, including many-body
electronic systems [19, 20, 21], Bose–Einstein condensates [22, 23, 24], and
classical optical [25, 26, 27, 28, 29, 30] and acoustic systems [31, 32, 33].
The Anderson transition is a continuous quantum phase transition. Like other
second-order phase transitions in statistical physics, the quantum phase
transitions are categorized by universality classes. Each universality class
is characterized by universal critical exponent and scaling functions, the
information of which is encoded into scaling properties of an effective field
theory around its fixed points. It is widely believed that the universality
class of the Anderson transition is determined only by symmetries of
disordered Hamiltonian (random matrix) and spatial dimension of the system [4,
5].

| class | TRS | PHS | CS | $\nu$ | Ref.
---|---|---|---|---|---|---
Unitary | A | No | No | No | $1.443\pm 0.003$ | [34]
Orthogonal | AI | $1$ | No | No | $1.572\pm 0.003$ | [35]
Symplectic | AII | $-1$ | No | No | $1.37\pm 0.01$ | [36]
Chiral Unitary | AIII | No | No | Yes | $1.06\pm 0.02$ | *
Chiral Orthogonal | BDI | $1$ | $1$ | Yes | $1.12\pm 0.06$ | *
Chiral Symplectic | CII | $-1$ | $-1$ | Yes | – | –
| D | No | 1 | No | $0.87\pm 0.03$ | *
BdG | C | No | $-1$ | No | $0.996\pm 0.012$ | *
| DIII | $-1$ | $1$ | Yes | $1.1\pm 0.05$ | [37]
| CI | 1 | $-1$ | Yes | $1.17\pm 0.02$ | *
Table 1: Classification according to time reversal symmetry (TRS), particle-
hole symmetry (PHS), chiral symmetry (CS), and their universal critical
exponents $\nu$ in 3D Anderson transition. TRS is defined as
$\mathbb{T}\mathbb{H}^{T}\mathbb{T}^{-1}=\mathbb{H}$ and PHS as
$\mathbb{C}\mathbb{H}^{T}\mathbb{C}^{-1}=-\mathbb{H}$. For TRS (PHS) 1 mean
$\mathbb{T}^{T}=\mathbb{T}$ ($\mathbb{C}^{T}=\mathbb{C})$ whereas $-1$ means
$\mathbb{T}^{T}=-\mathbb{T}$ $(\mathbb{C}^{T}=-\mathbb{C})$. * indicates the
values estimated in this paper, whereas – indicates yet to be determined.

The classification of Hermitian random matrices leads to ten symmetry classes
known as the Altland-Zirnbauer (AZ) classes. They consist of standard Wigner-
Dyson symmetry classes (A, AI, AII), three chiral symmetry classes (AIII, BDI,
CII) [9, 10], and four Bogoliubov-de Gennes (BdG) symmetry classes (D, C,
DIII, CI) [18, 38]. The random matrices and disordered Hamiltonian in the non-
standard symmetry classes not only describe universal aspects of the
localization of Bogoliubov quasiparticle wavefunctions in disordered
superconductors [18, 38] but also they are closely related to the localization
of classical waves in random dissipative media [39], non-Hermitian disordered
systems, which are attracting a lot of research interests during the last
couple of years [40, 41, 42, 43, 44, 45, 46, 47]. As some will be demonstrated
for the first time in this paper, three-dimensional (3D) models in each
symmetry class show the delocalization-localization transition and its quantum
criticality is characterized by its own universal scaling functions and the
critical exponent. So far, the critical exponent of the 3D Anderson transition
in seven out of the ten AZ classes have been studied numerically. They are the
three Wigner-Dyson classes A, AI and AII [34, 35, 36]; the chiral orthogonal
class BDI [48]; and three BdG classes C [49], CI [48], and DIII [37], while
the other three (AIII, CII, D) remain to be clarified. Numerical values of the
critical exponent of these classes, including those evaluated in this paper,
are listed in TABLE 1.

In this article, we clarify the universal critical exponents and scaling
functions associated with the localization length in the 3D Anderson
transition of the symmetry classes AIII, BDI, C, D and CI, using transfer
matrix method. The chiral class AIII has topological phases classified by
integer winding number in $d=1,3$ dimension(s). The study of the class AIII
was rather established both numerically and analytically in 1D [50, 51, 52,
53, 54], but most of them focused on topological quantum critical point and
few were directed to localization phenomena in 3D class AIII systems. Class C
and D in 3D are also less studied than in 2D, where there exists 2D spin and
thermal quantum Hall insulator phases in these two symmetry classes,
respectively. Low-energy Bogoliubov excitations in spin-triplet
superconductors with broken time-reversal symmetry Sr2RuO4 [55, 56] and UPt3
[57] potentially belong to these two non-standard symmetry classes in 3D.

In this paper, we propose two-band tight-binding models on the cubic lattice
that belong to symmetry classes AIII, BDI C, D and CI respectively, which
enable precise simulation studies of the 3D Anderson transition. Critical
exponents of the Anderson transition in the 3D classes AIII and D are
estimated for the first time, whose values are distinct from those of the
other known symmetry classes. Our estimates of exponents for the classes C and
CI agree with previous works of a classical network model [49] and a lattice
model [48].

##  II Unitary models with particle-hole or chiral symmetries

Figure 1:  Fermi surfaces of the two-bands models in the five non-standard
symmetry classes at $E=0$ in the clean limit. The yellow surfaces are of the
upper band, and the blue surfaces are of the lower band.

Let us first introduce the following two-orbital tight-binding model on the
cubic lattice that is shared by all the three unitary models studied in this
paper:

$\displaystyle H_{0}$
$\displaystyle\equiv\sum_{{\bm{i}},{\bm{j}}}\sum_{d,d^{\prime}=a,b}|{\bm{i}},d\rangle\big{[}\mathbb{H}_{0}\big{]}_{({\bm{i}},d|{\bm{j}},d^{\prime})}\langle{\bm{j}},d|$
$\displaystyle=\sum_{\bm{i}}\left(\epsilon_{\bm{i}}+\Delta\right)\left(|{\bm{i}},a\rangle\langle{\bm{i}},a|-|{\bm{i}},b\rangle\langle{\bm{i}},b|\right)\\!\
+$
$\displaystyle\sum_{\bm{i}}\bigg{\\{}t_{\perp}\Big{[}\mathrm{i}\left(|{\bm{i}}+{\bm{e}}_{x},a\rangle\langle{\bm{i}},b|+|{\bm{i}}+{\bm{e}}_{x},b\rangle\langle{\bm{i}},a|\right)$
$\displaystyle\ \ \
+\left(|{\bm{i}}+{\bm{e}}_{y},a\rangle\langle{\bm{i}},b|-|{\bm{i}}+{\bm{e}}_{y},b\rangle\langle{\bm{i}},a|\right)\Big{]}$
$\displaystyle\
+t_{\parallel}\left(|{\bm{i}}+{\bm{e}}_{z},a\rangle\langle{\bm{i}},a|+|{\bm{i}}+{\bm{e}}_{z},b\rangle\langle{\bm{i}},b|\right)+\mathrm{H.c.}\bigg{\\}}.$
(1)

Here ${\bm{i}}\equiv(i_{x},i_{y},i_{z})$ specifies the cubic-lattice site;
${\bm{e}}_{x}\equiv(1,0,0)$, ${\bm{e}}_{y}\equiv(0,1,0)$, and
${\bm{e}}_{z}\equiv(0,0,1)$. The lattice constant is taken to be the unit
length. $a$ and $b$ are the two orbitals: $d,d^{\prime}=a,b$. The nearest
neighbor (NN) hopping within the $x$-$y$ plane is an inter-orbital hopping
$t_{\perp}$, and the NN hopping along $z$ direction is an intra-orbital
hopping $t_{\parallel}$. $\pm\Delta$ is an onsite energy for the $a,b$
orbitals respectively, and $\epsilon_{\bm{i}}$ is a random potential
independently distributed on each lattice site. In numerical simulations, we
choose the distribution of the random potential to be uniform:
$\epsilon_{\bm{i}}\in[-W/2,W/2],~{}\overline{\epsilon_{\bm{i}}\epsilon_{\bm{j}}}=\delta_{\bm{ij}}W^{2}/12$.
The tight binding Hamiltonian has the following chiral symmetries (Eqs. (2,3))
and particle-hole symmetries (Eqs. (4),(5)):

$\displaystyle\mathbb{P}_{1}\mathbb{H}_{0}\mathbb{P}^{-1}_{1}$
$\displaystyle=-\mathbb{H}_{0},\ \
\mathbb{P}_{1}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{x}+i_{z}}[\sigma_{2}]_{d,d^{\prime}},$
(2) $\displaystyle\mathbb{P}_{2}\mathbb{H}_{0}\mathbb{P}^{-1}_{2}$
$\displaystyle=-\mathbb{H}_{0},\ \
\mathbb{P}_{2}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{y}+i_{z}}[\sigma_{1}]_{d,d^{\prime}},$
(3) $\displaystyle\mathbb{C}_{1}\mathbb{H}^{T}_{0}\mathbb{C}^{-1}_{1}$
$\displaystyle=-\mathbb{H}_{0},\ \
\mathbb{C}_{1}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{z}}[\sigma_{1}]_{d,d^{\prime}},$
(4) $\displaystyle\mathbb{C}_{2}\mathbb{H}^{T}_{0}\mathbb{C}^{-1}_{2}$
$\displaystyle=-\mathbb{H}_{0},\ \
\mathbb{C}_{2}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{x}+i_{y}+i_{z}}[\sigma_{2}]_{d,d^{\prime}},$
(5)

with $\mathbb{P}_{2}\equiv\mathbb{P}_{1}\mathbb{C}_{1}\mathbb{C}_{2}$.
$\sigma_{1,2,3}$ are the Pauli matrices in the $a,b$ orbital space and
$\sigma_{0}$ is the unit matrix. The combination of particle-hole symmetries
and chiral symmetries gives the time-reversal symmetries, e.g.
$(\mathbb{P}_{1}\mathbb{C}_{1})^{T}\mathbb{H}^{T}_{0}\mathbb{P}_{1}\mathbb{C}_{1}=\mathbb{H}_{0}$.
$\mathbb{P}_{1}\mathbb{P}_{2}\equiv\mathbb{C}_{1}\mathbb{C}_{2}=(-1)^{i_{x}+i_{y}}\mathrm{i}\sigma_{z}$
commutes with $\mathbb{H}_{0}$, where the tight-binding Hamiltonian can be
block-diagonalized in a basis of the real-valued eigenvectors of
$\mathbb{P}_{1}\mathbb{P}_{2}$. Since $\mathbb{P}_{1}\mathbb{C}_{1}$ and
$\mathbb{P}_{1}\mathbb{C}_{2}$ commute with $\mathbb{P}_{1}\mathbb{P}_{2}$,
$\mathbb{P}_{1}\mathbb{C}_{1}$ and $\mathbb{P}_{1}\mathbb{C}_{2}$ are block-
diagonalized in the same basis, too. Accordingly, the two blocks of the tight-
binding Hamiltonian are time-reversal symmetric. Since
$(\mathbb{P}_{1}\mathbb{C}_{1})^{*}(\mathbb{P}_{1}\mathbb{C}_{1})=+1$, $H_{0}$
belongs to the orthogonal class. On the basis of such $H_{0}$, we will next
introduce other parts of the hopping terms, which lead to the three unitary
models belonging to the three non-standard symmetry classes one by one.

###  II.1 Class AIII model

To break the two particle-hole symmetries, we add the following nearest
neighbor hoppings within the $x$-$y$ plane;

$\displaystyle H_{\rm AIII}$
$\displaystyle\equiv\sum_{{\bm{i}},{\bm{j}}}\sum_{d,d^{\prime}=a,b}|{\bm{i}},d\rangle\big{[}\mathbb{H}_{\rm
AIII}\big{]}_{({\bm{i}},d|{\bm{j}},d^{\prime})}\langle{\bm{j}},d^{\prime}|$
$\displaystyle=H_{0}+\sum_{\bm{i}}\Big{[}t_{1}\left(|{\bm{i}}+{\bm{e}}_{x},a\rangle\langle{\bm{i}},a|-|{\bm{i}}+{\bm{e}}_{x},b\rangle\langle{\bm{i}},b|\right).$
$\displaystyle+t_{2}\left(|{\bm{i}}+{\bm{e}}_{y},a\rangle\langle{\bm{i}},a|+|{\bm{i}}+{\bm{e}}_{y},b\rangle\langle{\bm{i}},b|\right)+\mathrm{H.c.}\Big{]}.$
(6)

The model respects the chiral symmetry in Eq. (3), while the two particle-hole
symmetries are broken by $t_{2}$ and $t_{1}$, respectively. We set
$\Delta=0.0$, $t_{\perp}=0.6$, $t_{\parallel}=0.4$, and $t_{1}=t_{2}=0.5$. The
Hamiltonian in the clean limit is Fourier-transformed to

$\displaystyle H_{\mathrm{AIII}}(\bm{\mathrm{k}})=$ $\displaystyle
2(t_{\parallel}\cos k_{z}+t_{2}\cos k_{y})\sigma_{0}+(\Delta+2t_{1}\cos
k_{x})\sigma_{3}$ $\displaystyle-2t_{\perp}(\sin k_{x}\sigma_{1}+\sin
k_{y}\sigma_{2}),$ (7)

where the two bands are particle-hole symmetric under the translation of the
momentum by $(0,\pi,\pi)$. The upper (lower) energy band forms an electron
(hole) pocket at $E=0$ around the zone boundary (center) axis of
$(k_{y},k_{z})=(\pi,\pi)$ ($(0,0)$) (Fig. 1(a)).

###  II.2 Class D model

To break the chiral symmetries and one of the two particle-hole symmetries
(Eq. (5)), we add to $H_{0}$ the following NN hoppings within the $x$-$y$
plane:

$\displaystyle H_{\rm D}$
$\displaystyle\equiv\sum_{{\bm{i}},{\bm{j}}}\sum_{d,d^{\prime}=a,b}|{\bm{i}},d\rangle\big{[}\mathbb{H}_{\rm
D}\big{]}_{({\bm{i}},d|{\bm{j}},d^{\prime})}\langle{\bm{j}},d^{\prime}|$
$\displaystyle=H_{0}+\sum_{\bm{i}}\Big{[}t_{1}\left(|{\bm{i}}+{\bm{e}}_{x},a\rangle\langle{\bm{i}},a|-|{\bm{i}}+{\bm{e}}_{x},b\rangle\langle{\bm{i}},b|\right).$
$\displaystyle+t_{2}\left(|{\bm{i}}+{\bm{e}}_{y},a\rangle\langle{\bm{i}},a|-|{\bm{i}}+{\bm{e}}_{y},b\rangle\langle{\bm{i}},b|\right)+\mathrm{H.c.}\Big{]}.$
(8)

The model respects the particle-hole symmetry in Eq. (4), while the two chiral
symmetries are broken by $t_{1}$ and $t_{2}$, respectively. We set
$\Delta=0.1$, $t_{\perp}=0.3$, $t_{\parallel}=0.2$, and $t_{1}=t_{2}=0.5$.
Since $\mathbb{C}^{T}_{1}=\mathbb{C}_{1}$, the model belongs to the symmetry
class D. The tight-binding model can be regarded as layered Chern insulator
models [58, 59]. In momentum space, the clean-limit Hamiltonian takes the
following form

$\displaystyle H_{\mathrm{D}}(\bm{\mathrm{k}})$
$\displaystyle=2t_{\parallel}\cos k_{z}\sigma_{0}-2t_{\perp}(\sin
k_{x}\sigma_{1}+\sin k_{y}\sigma_{2})$ $\displaystyle+(\Delta+2t_{1}\cos
k_{x}+2t_{2}\cos k_{y})\sigma_{3}.$ (9)

A finite $t_{\parallel}$ closes the band gap at $E=0$, where electron and hole
pockets appears at $k_{z}=0$ and $\pi$, respectively, as shown in Fig. 1(b).

###  II.3 Class C model

To construct a class C model from Eq. (1), we add to $H_{0}$ the following
nearest neighbor hoppings within the $x$-$y$ plane:

$\displaystyle H_{\rm C}$
$\displaystyle\equiv\sum_{{\bm{i}},{\bm{j}}}\sum_{d,d^{\prime}=a,b}|{\bm{i}},d\rangle\big{[}\mathbb{H}_{\rm
C}\big{]}_{({\bm{i}},d|{\bm{j}},d^{\prime})}\langle{\bm{j}},d^{\prime}|$
$\displaystyle=H_{0}+\sum_{\bm{i}}\Big{[}t_{1}\left(|{\bm{i}}+{\bm{e}}_{x},a\rangle\langle{\bm{i}},a|+|{\bm{i}}+{\bm{e}}_{x},b\rangle\langle{\bm{i}},b|\right).$
$\displaystyle+t_{2}\left(|{\bm{i}}+{\bm{e}}_{y},a\rangle\langle{\bm{i}},a|+|{\bm{i}}+{\bm{e}}_{y},b\rangle\langle{\bm{i}},b|\right)+\mathrm{H.c.}\Big{]}.$
(10)

The model respects the particle-hole symmetry Eq. (5) with
$\mathbb{C}^{T}_{2}=-\mathbb{C}_{2}$, while finite $t_{1}$ and $t_{2}$ break
Eq. (3) and (2), respectively. Therefore the model belongs to symmetry class
C. In this paper we set $\Delta=0.1$, $t_{\perp}=0.6$, $t_{\parallel}=0.4$,
and $t_{1}=t_{2}=0.5$. In the clean limit, this Hamiltonian is Fourier-
transformed to

$\displaystyle H_{\mathrm{C}}(\bm{\mathrm{k}})$
$\displaystyle=2(t_{\parallel}\cos k_{z}+t_{1}\cos k_{x}+t_{2}\cos
k_{y})\sigma_{0}+\Delta\sigma_{3}$ $\displaystyle-2t_{\perp}(\sin
k_{x}\sigma_{1}+\sin k_{y}\sigma_{2}).$ (11)

The Fermi surface at $E=0$ is shown in Fig. 1(c).

The above models respect only one of the four (particle-hole or chiral)
symmetries, Eqs. (2), (3), (4), or (5). If there were time reversal symmetry
$\mathbb{T}$, $\mathbb{T}\mathbb{H}^{T}\mathbb{T}^{-1}=\mathbb{H}$, combining
$\mathbb{T}$ and the particle-hole or chiral symmetry should give rise to
another chiral or particle-hole symmetry. Having only one particle-hole or
chiral symmetry, therefore, means TRS is broken. In Appendix A, we explicitly
demonstrate the broken-time-reversal symmetries in the these three models, by
showing the non-zero Hall conductivity calculated from the Berry phase of
Bloch bands.

##  III Orthogonal models with particle-hole or chiral symmetries

In this chapter, we introduce two-orbital tight-binding models that belong to
two orthogonal classes CI and BDI with chiral and particle-hole symmetries.
The two models are real-valued and symmetric, hence chiral and particle-hole
symmetries are equivalent. In the following two models, the random potential
$\epsilon_{\boldsymbol{i}}$ takes the same uniform distribution as in Eq. (1).
The class CI model is given as follows,

$\displaystyle H_{\mathrm{CI}}$
$\displaystyle\equiv\sum_{\boldsymbol{i},{\boldsymbol{j}}}\sum_{d,d^{\prime}=a,b}|{\boldsymbol{i}},d\rangle\big{[}\mathbb{H}_{\mathrm{CI}}\big{]}_{({\boldsymbol{i}},d|{\boldsymbol{j}},d^{\prime})}\langle{\boldsymbol{j}},d^{\prime}|$
$\displaystyle=\sum_{\boldsymbol{i}}\left(\epsilon_{\boldsymbol{i}}+\Delta\right)\left(|{\boldsymbol{i}},a\rangle\langle{\boldsymbol{i}},b|+|{\boldsymbol{i}},b\rangle\langle{\boldsymbol{i}},a|\right)\\!\
$
$\displaystyle+\sum_{\boldsymbol{i}}\bigg{\\{}\sum_{\mu=x,y}t_{\perp}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{\mu},a\rangle\langle{\boldsymbol{i}},a|+|{\boldsymbol{i}}+{\boldsymbol{e}}_{\mu},b\rangle\langle{\boldsymbol{i}},b|\big{)}$
$\displaystyle+t_{\parallel}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},a\rangle\langle{\boldsymbol{i}},a|-|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},b\rangle\langle{\boldsymbol{i}},b|\big{)}$
$\displaystyle+t_{\parallel}^{\prime}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},a\rangle\langle{\boldsymbol{i}},b|+|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},b\rangle\langle{\boldsymbol{i}},a|\big{)}+\mathrm{H.c.}\bigg{\\}}.$
(12)

It has only the following particle-hole (chiral) symmetry,

$\displaystyle\mathbb{C}_{3}\mathbb{H}^{T}_{0}\mathbb{C}^{-1}_{3}$
$\displaystyle=-\mathbb{H}_{0},\ \
\mathbb{C}_{3}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{x}+i_{y}}[\sigma_{2}]_{d,d^{\prime}},$
(13)

with $\mathbb{C}_{3}^{T}=-\mathbb{C}_{3}$. Note that there is no other
particle-hole symmetries. Namely, $\mathbb{C}_{3}$ must be diagonal with
respect to the cubic-lattice site index, because of the on-site random terms.
From the on-site term, it has to take a form of either $\sigma_{2}$ or
$\sigma_{3}$ in the orbital space. From the nearest-neighbor hoppings within
the $x$-$y$ plane, it comes with a U(1) phase of
$(-1)^{i_{x}+i_{y}}e^{i\theta(i_{z})}$. From the hopping along $z$, such
$\mathbb{C}_{3}$ is uniquely determined by $(-1)^{i_{x}+i_{y}}\sigma_{2}$,
hence the model belongs to the symmetry class CI. In this paper, we take
$\Delta=t_{\perp}=t_{\parallel}=1,\ t_{\parallel}^{\prime}=2$. The Fermi
surface at $E=0$ is shown in Fig. 1(d).

The class BDI model is given by

$\displaystyle H_{\mathrm{BDI}}$
$\displaystyle\equiv\sum_{\boldsymbol{i},{\boldsymbol{j}}}\sum_{d,d^{\prime}=a,b}|{\boldsymbol{i}},d\rangle\big{[}\mathbb{H}_{\mathrm{BDI}}\big{]}_{({\boldsymbol{i}},d|{\boldsymbol{j}},d^{\prime})}\langle{\boldsymbol{j}},d^{\prime}|$
$\displaystyle=\sum_{\boldsymbol{i}}\left(\epsilon_{\boldsymbol{i}}+\Delta\right)\left(|{\boldsymbol{i}},a\rangle\langle{\boldsymbol{i}},a|-|{\boldsymbol{i}},b\rangle\langle{\boldsymbol{i}},b|\right)\\!\
$
$\displaystyle+\sum_{\boldsymbol{i}}\bigg{\\{}\sum_{\mu=x,y}t_{\perp}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{\mu},a\rangle\langle{\boldsymbol{i}},a|+|{\boldsymbol{i}}+{\boldsymbol{e}}_{\mu},b\rangle\langle{\boldsymbol{i}},b|\big{)}$
$\displaystyle+t_{\parallel}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},a\rangle\langle{\boldsymbol{i}},a|-|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},b\rangle\langle{\boldsymbol{i}},b|\big{)}$
$\displaystyle+t_{\parallel}^{\prime}\big{(}|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},a\rangle\langle{\boldsymbol{i}},b|-|{\boldsymbol{i}}+{\boldsymbol{e}}_{z},b\rangle\langle{\boldsymbol{i}},a|\big{)}+\mathrm{H.c.}\bigg{\\}}.$
(14)

It has only the following chiral (particle-hole) symmetry;

$\displaystyle\mathbb{P}_{3}\mathbb{H}_{\rm BDI}\mathbb{P}^{-1}_{3}$
$\displaystyle=-\mathbb{H}_{\rm BDI},\ \
\mathbb{P}_{3}=\delta_{{\bm{i}},{\bm{j}}}(-1)^{i_{x}+i_{y}}[\sigma_{1}]_{d,d^{\prime}}.$
(15)

with $\mathbb{P}_{3}^{T}=\mathbb{P}_{3}$. Namely, from the first two terms of
Eq. (14), $\mathbb{P}_{3}$ must take a form of either
$(-1)^{i_{x}+i_{y}}e^{i\theta(i_{z})}\sigma_{1}$ or
$(-1)^{i_{x}+i_{y}}e^{i\theta(i_{z})}\sigma_{2}$. The last two terms uniquely
determine $\mathbb{P}_{3}$ to be $(-1)^{i_{x}+i_{y}}\sigma_{1}$; the model
belongs to the symmetry class BDI. Here we take
$\Delta=t_{\perp}=t_{\parallel}^{\prime}=1,\ t_{\parallel}=0.5$. The Fermi
surface at $E=0$ is shown in Fig. 1(e)

We emphasize that the BdG symmetry classes describe the localization of
Bogoliubov quasiparticle wavefunctions in superconductors [18], where
disordered single-particle Hamiltonian is given by a mixed basis in particle
and hole space. In the class C, D and CI models introduced above, the particle
and hole degrees of freedom are undertaken by the two orbitals. Thereby, the
concept of a finite area of Fermi surface applies to single-particle
wavefunctions of these two-orbitals models. We expect the critical behaviors
of the Anderson transition of the single-particle wavefunctions are the same
as those of the superconducting quasiparticle wavefunctions.

##  IV Method and Result

Figure 2:  Plots of normalized localization length (or its inverse or its
logarithm) of the five models near the critical point of the Anderson
transition for different system size in the transverse direction (see text).
$\Lambda$ is the dimensionless localization length and $\Gamma=\Lambda^{-1}$.
Black dots are numerical data points with error bar and colored curves are the
fitting results with the largest GOF in TABLE 2. The quasi-one-dimensional
samples are typically of length $L_{z}=10^{6}-10^{7}$ to ensure a 0.2 %
precision of each data point. Figure 3:  The single-parameter scaling function
of the normalized localized length for the five symmetry classes.
$\Lambda_{\mathrm{corrected}}$ is the numerical data $\Lambda$ subtracted by
all terms of irrelevant scaling variable in Eq. (18). The upper branch is for
the delocalized phase and the lower branch is for the insulator phase. We note
that $\Lambda_{c}$ depends on the anisotropy of the model.

The scaling theory is at the heart of studying quantum criticality of the
Anderson localization. On approaching the critical point, various physical
quantities diverge according to a power law with the universal critical
exponent $\nu$. One of them is the correlation length
$\xi\sim|x-x_{c}|^{-\nu}$, which in a localized phase, characterises a decay
length of a density-density correlation function in disordered medium. Here
$x$ represents the coordinate in a certain parameter space of the system and
$x_{c}$ stands for the critical point. To determine the critical exponent, we
calculate a quasi-one-dimensional (Q1D) localization length $\lambda$ of a
single-particle wavefunction amplitude, using transfer matrix method [11, 12].
To be specific, we consider the Q1D geometry of the 3D cubic lattice
($L_{x}\times L_{y}\times L_{z}$) with a cross-section $L_{x}=L_{y}=L\ll
L_{z}$. In such geometry, all the eigenstates of the disordered Hamiltonian
are localized in $z$ direction, $|\psi|\sim\exp(-|z-z_{0}|/\lambda)$, with the
localization length $\lambda$. The scaling argument around the critical point
suggests that the localization length normalized by the linear dimension of
the cross-section, $\Lambda\equiv\lambda/L$ or $\Gamma=\Lambda^{-1}$, must be
scale-invariant (independent of $L$) at the critical point (see Fig. 2).
Namely, numerical data of $\Lambda$ for different cross-sectional size $L$
must fit in so-called the one-parameter scaling function near the transition
point [11, 12, 13],

$\displaystyle\Lambda(W,L)=f(L/\xi(W))$ (16)

where $f(x)$ is a universal scaling function for the normalized localization
length. $\xi(W)$ is the correlation length that depends on disorder strength:
$\xi(W)\sim|W-W_{c}|^{-\nu}$. In practical calculations, one often encounters
deviations of the data sets from the single-parameter scaling form. The
deviations are attributed to finite-size effects associated with irrelevant
scaling variables around the critical point in the framework of a
renormalization group (RG) theory of general critical phenomena. In the RG
theory, generic continuous phase transition is characterized by a saddle-point
fixed point of RG equations for a certain effective theory with multiple
system parameters. When high-energy degrees of freedom are integrated out,
low-energy system parameters run away from the fixed point only along one
direction in the parameter space, while they flow into the fixed point along
all the other directions. The unstable direction around the fixed point is
characterized by relevant scaling variable, while its complementary directions
are characterized by a number of irrelevant scaling variables. Scaling
dimension of the relevant scaling variable is positive and it defines the
universal critical exponent, $1/\nu$. Scaling dimensions of the irrelevant
scaling variables are all negative, where the irrelevant variable with the
largest negative scaling dimension is the least irrelevant. Empirically, we
attribute the deviation of the numerical data from the single-parameter
scaling to the least irrelevant scaling variable [14, 60]. The scaling
argument that includes the effect of the least irrelevant scaling variable
modifies the single-parameter scaling function form into

$\displaystyle\Lambda(W,L)=F\left(u_{1}(w)L^{1/\nu},u_{2}(w)L^{-y}\right),$
(17)

where $w=(W-W_{c})/W_{c}$, and $1/\nu>0$ and $-y<0$ are scaling dimensions of
the relevant scaling variable $u_{1}$ and the least irrelevant scaling
variable $u_{2}$, respectively. The omission of the irrelevant scaling
variables with smaller negative scaling dimensions are justified a posterior
by large $y$ obtained from the fitting (see TABLE 2).

To obtain the exponents $\nu$, $y$ as well as the scaling function $f(x)\equiv
F(x^{1/\nu},0)$, we expand the scaling function $F$ as Taylor series of its
arguments

$\displaystyle\Lambda=\sum_{i=0}^{n_{1}}\sum_{j=0}^{n_{2}}a_{ij}\left(u_{1}L^{1/\nu}\right)^{i}\left(u_{2}L^{-y}\right)^{j}.$
(18)

When the system is close to the critical point, the scaling variables can be
also expanded in small $w$

$\displaystyle u_{1}=\sum_{k=1}^{m_{1}}b_{k}w^{k},\quad
u_{2}=\sum_{k=0}^{m_{2}}c_{k}w^{k}.$ (19)

To fix the ambiguity of fitting, we at $a_{1,0}=a_{0,1}=1$ [14, 60]. Having
the polynomial form of the scaling function $F$ in $w$, we use the least-
squares fitting procedure and minimize $\chi^{2}$ defined as follows,

$\displaystyle\chi^{2}=\sum_{i=1}^{N_{D}}\frac{(\Lambda_{i}-F_{i})^{2}}{\sigma_{i}^{2}}.$
(20)

Here $\Lambda_{i},\,\sigma_{i},\,F_{i}$ are the numerical mean value of
$\Lambda$, numerical standard error of $\Lambda$, and a fitting value of
$\Lambda$ from the polynomial fitting function. They are given at every data
point specified by $L$ and $W$ with the data point index $i=1,\cdots,N_{D}$.
The fitting results give estimate of $\nu,\,y,\,W_{c},\,\Lambda_{c}$ and
expansion coefficients $\\{a_{ij}\\},\,\\{b_{k}\\},\,\\{c_{k}\\}$. We choose
$n_{1}=3$, $n_{2}=1$ and perform the polynomial fitting for different values
of $m_{1}$ and $m_{2}$. Results are shown in TABLE 2. They are stable against
changes of the expansion order $m_{1}$, $m_{2}$, as well as the range of
system sizes. The goodness of fit (GOF) is a probability that the $N_{D}$ data
points sampled from the fitting function would give a $\chi^{2}$ larger than
the minimized value $\chi^{2}_{\text{min}}$ of the fitting [61]. We use
$\Gamma$ as data inputs for the evaluations in the class AIII, C and CI
models, and $\ln\Lambda$ for the class D and BDI models, in order to get
better fitting qualities. The resulting critical exponents of the five
symmetry classes are

$\displaystyle\nu_{\text{AIII}}\ =$ $\displaystyle\ 1.06\pm 0.02,$
$\displaystyle\nu_{\text{D}}\ =$ $\displaystyle\ 0.87\pm 0.03,$
$\displaystyle\nu_{\text{C}}\ =$ $\displaystyle\ 0.996\pm 0.012,$
$\displaystyle\nu_{\text{CI}}\ =$ $\displaystyle\ 1.17\pm 0.02,$
$\displaystyle\nu_{\text{BDI}}\ =$ $\displaystyle\ 1.12\pm 0.06.$ (21)

The numbers after $\pm$ sign are the standard deviation determined from 1000
fittings of synthetic data sets.

The values of the evaluated critical exponents in the three non-standard
symmetry classes with broken time-reversal symmetry are distinct from each
other and also from previously evaluated critical exponents in other symmetry
classes. The symmetry class C in 3D has been studied by network model of both
quantum systems [62] and the classical counterpart [49]. The evaluated
critical exponent in Ref. [49] is $\nu=0.9985\pm 0.0015$, which is in good
agreement with our result. The symmetry class CI in 3D has also been studied
with a lattice model [48] as $\nu=1.16\pm 0.02$, also in good agreement with
our result. These consistencies demonstrate that the universality class of the
Anderson transitions are free from details of the models, and it is determined
only by the symmetry of the random matrices and the spatial dimension.

In Fig. 3, we show the scaling functions of these five non-standard symmetry
classes. We define $\Lambda_{\text{corrected}}$ from raw numerical data
subtracted by the contributions of the irrelevant scaling variable in the
fitting function $F$. $\Lambda_{\text{corrected}}$ should obey the single-
parameter scaling form in Eq. (16) around the critical point, with the scaling
argument $L/\xi\sim L|u_{1}|^{\nu}$. When plotted as a function of
$L|u_{1}|^{\nu}$, $\Lambda_{\text{corrected}}$ all collapse onto two branches
of the scaling functions, where the upper (lower) branch corresponds to the
metallic (localized) phase regime.

Table 2:  Finite-size scaling analyses of the normalized localization length
around the Anderson transition in the five non-standard symmetry classes. The
scaling function of the normalized localization length is expanded in power of
$w\equiv(W-W_{c})/W_{c}$ (see text). The expansion orders $n_{1}=3,\ n_{2}=1$
are fixed, while $m_{1},m_{2}$ are variables. The values inside the square
brackets are 95% confidence intervals from 1000 Monte Carlo simulations. Only
fitting samples with good of fit (GOF) greater than 0.1 are shown. We note
that our model is anisotropic, and $\Lambda_{c}$ depends on the anisotropy.

##  V Summary and discussion

In conclusion, we have clarified the universal critical exponents and scaling
functions associated with the Q1D localization length of the 3D conventional
localization-delocalization transition in three non-standard unitary symmetry
classes (the chiral unitary classes AIII, the BdG class D and C), and in two
non-standard orthogonal classes (the chiral orthogonal class BDI, the BdG
class CI).

We introduced two-orbital tight-binding models on the cubic lattice that
belong to these symmetry classes. In these models, the particle-hole degrees
of freedom in usual BdG Hamiltonians for superconductors [18] and the
sublattice degrees of freedom in usual chiral symmetry [9, 63, 64, 10, 65] are
overtaken by an artificial orbital degrees of freedom. One of the advantages
of these tight-binding models is that all of them have finite area of Fermi
surfaces in the clean limit. The presence of the finite Fermi surfaces in the
clean limit guarantees the presence of 3D Anderson transitions at finite
disorder strength. Using comprehensive numerical analyses, we (re)determine
the critical exponents of the Anderson transition in these five non-standard
symmetry classes. The exponents of the 3D class AIII and D are evaluated in
this paper for the first time. The exponents of class C and CI are consistent
with value in previous literature.

On the contrary, our exponent of the class BDI model is significantly
different from a previous evaluation in the nodal Dirac semimetal (NDS) model
[48]. The disordered NDS model has two topologically different types of
delocalization-localization transitions in its phase diagram. One is a
‘topological’ phase transition line between diffusive metal (DM) phase and
topological insulator (TI) phase with 1D topological winding number in the BDI
class. The other is a ‘conventional’ phase transition line between the DM
phase and ordinary (topologically trivial) band insulator or Anderson
localized phase. The previous evaluation of the critical exponents at these
two types of the phase transition lines are consistent with each other. On the
one hand, the zero-energy density of state (DOS) suggests that the dynamical
exponent at these two transitions could be different; the zero-energy DOS
takes a finite constant value in the former ‘topological’ phase transition
line, while the DOS shows a weak singularity at the zero energy around the
‘conventional’ phase transition line [48]. The singularity of the zero-energy
DOS indicates a possible deviation of the dynamical exponent from the spatial
dimension. It is an interesting open issue to study how the zero-energy DOS
behaves at the Anderson transition point in the present BDI model. Such
information could provide comprehensive understanding on the 3D Anderson
transition in the BDI class.

Recently, the critical exponents and scaling functions in the non-standard
symmetry classes acquire a lot of research interests from a view point of
conventional and unconventional delocalization-localization transitions in
non-Hermitian disordered systems. A previous dimensional regularization study
of the non-linear sigma models in $d=2+\epsilon$ shows that the $\beta$
function in the three chiral symmetry classes is identical to zero [66, 9,
10]. The vanishing $\beta$ function implies an unusual nature of
delocalization-localization transition in these chiral symmetry classes [67].
Our evaluation of the universal critical properties provide solid information
to these research community and will shed a new light in future experiments on
quasiparticle heat transports in superconductors as well as the localization
phenomena in random dissipative systems.

###### Acknowledgements.

T. W. and R. S. thanks the fruitful discussion and correspondence with Xunlong
Luo. T. O. was supported by JSPS KAKENHI Grants 19H00658. T. W. and R. S. was
supported by the National Basic Research Programs of China (No.
2019YFA0308401) and by National Natural Science Foundation of China
(No.11674011 and No. 12074008).

##  Appendix A Berry Curvature and Time-reversal Symmetry

$H$ is time-reversal invariant, when it commutes with the time reversal
operator $\Theta$,

$\displaystyle[\Theta,H]=0,\quad\Theta=\mathbb{T}\cdot K\,,$ (22)

where $K$ is the complex conjugation and $\mathbb{T}$ a unitary operator.
Since the tight-binding models studied in this paper contain on-site random
potentials, the unitary transformation $\mathbb{T}$ must be diagonal with
respect to the lattice sites. To prove the absence of such $\mathbb{T}$ in the
class AIII and C models, we calculate the Hall conductivity at $E=0$ in the
clean limit.

Two-bands Hamiltonians in the momentum space are written in the following
form,

$\displaystyle
H(\bm{\mathrm{k}})=E_{0}(\bm{\mathrm{k}})\sigma_{0}+E_{1}(\bm{\mathrm{k}})\bm{n}({\bm{k}})\cdot\bm{\sigma}.$
(23)

with a three-component unit vector ${\bm{n}}(\bm{k})$ and the two by two Pauli
matrices ${\bm{\sigma}}$. For the 2 by 2 Hamiltonian, the two energy bands and
the Berry curvature for each band are given by

$\displaystyle E_{\pm}(\bm{\mathrm{k}})$
$\displaystyle=E_{0}(\bm{\mathrm{k}})\pm E_{1}(\bm{\mathrm{k}}),$
$\displaystyle\Omega_{xy}^{\pm}(\bm{\mathrm{k}})$
$\displaystyle=\pm\left(\frac{\partial\bm{n}}{\partial
k_{x}}\times\frac{\partial\bm{n}}{\partial k_{y}}\right)\cdot\bm{n}.$ (24)

In terms of this, the Berry curvature of the class AIII model is calculated
from Eq. (7),

$\displaystyle\Omega_{xy}^{\pm}(\bm{\mathrm{k}})=\pm\frac{4t_{\perp}^{2}}{E_{1}^{3}}(\Delta\cos
k_{x}+2t)\cos k_{y},$ (25)

with $E_{1}(\bm{\mathrm{k}})\equiv[(\Delta+2t\cos
k_{x})^{2}+4t_{\perp}^{2}(\sin^{2}k_{x}+\sin^{2}k_{y})]^{1/2}$. The Berry
curvature of the class D model is calculated from Eq. (9),

$\displaystyle\Omega_{xy}^{\pm}(\bm{\mathrm{k}})=\pm\frac{4t_{\perp}^{2}}{E_{1}^{3}}\left[\Delta\cos
k_{x}\cos k_{y}+2t(\cos k_{x}+\cos k_{y})\right],$ (26)

with $E_{1}(\mathbf{k})\equiv[(\Delta+2t\cos k_{x}+2t\cos
k_{y})^{2}+4t_{\perp}^{2}(\sin^{2}k_{x}+\sin^{2}k_{y})]^{1/2}$. Likewise, the
Berry curvature of the class C model is calculated from Eq. (11),

$\displaystyle\Omega_{xy}^{\pm}(\bm{\mathrm{k}})=\pm\frac{4t_{\perp}^{2}}{E_{1}^{3}}\Delta\cos
k_{x}\cos k_{y},$ (27)

with
$E_{1}(\bm{\mathrm{k}})\equiv[\Delta^{2}+4t_{\perp}^{2}(\sin^{2}k_{x}+\sin^{2}k_{y})]^{1/2}$.

An integral of the Berry curvature over the occupied bands ($E_{\pm}<0$) is
nothing but the Hall conductivity at $E=0$ at certain $k_{z}$ [68],

$\displaystyle\sigma_{xy}(k_{z})\equiv$
$\displaystyle\int_{E_{+}<0}\mathrm{d}^{2}k\\!\
\Omega_{xy}^{+}(\bm{\mathrm{k}})+\int_{E_{-}<0}\mathrm{d}^{2}k\\!\
\Omega_{xy}^{-}(\bm{\mathrm{k}}).$ (28)

Throughout this paper we take lattice constant to be one, and the first
Brillouin zone is $k_{x,y,z}\in[-\pi,\pi]$. By integrating over $k_{z}$, one
gets the total Hall conductivity
$\sigma_{xy}=\int_{-\pi}^{\pi}\mathrm{d}k_{z}\sigma_{xy}(k_{z})$. A non-zero
Hall conductivity indicates the broken time-reversal symmetry. With the chosen
parameters for numerical simulations of the three model, we get
$\sigma_{xy}/8\pi^{2}\approx 0.04$ for the class AIII model,
$\sigma_{xy}/8\pi^{2}\approx 0.73$ for the class D model, and
$\sigma_{xy}/8\pi^{2}\approx 0.58$ for the class C model.

## References

  * Anderson [1958] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * Abrahams [2010] E. Abrahams, ed.,  _50 Years of Anderson Localization_ (World Scientific Publishing Company, 2010).
  * Wegner [1976] F. J. Wegner, Zeitschrift fr Physik B Condensed Matter and Quanta 25, 327 (1976).
  * Abrahams  _et al._ [1979] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
  * Effetov  _et al._ [1980] K. B. Effetov, A. I. Larkin, and D. E. Khmel’nitsukii, Soviet Phys. JETP 52, 568 (1980).
  * Hikami [1981a] S. Hikami, Phys. Rev. B 24, 2671 (1981a).
  * Hikami [1981b] S. Hikami, Physics Letters B 98, 208 (1981b).
  * Gade and Wegner [1991] R. Gade and F. Wegner, Nucl. Phys. B 360, 213 (1991).
  * Gade [1993] R. Gade, Nucl. Phys. B 398, 499 (1993).
  * MacKinnon and Kramer [1981] A. MacKinnon and B. Kramer, Phys. Rev. Lett. 47, 1546 (1981).
  * Mackinnon and Kramer [1983] A. Mackinnon and B. Kramer, Zeitschrift Fur Physik B-Condensed Matter 53, 1 (1983).
  * Pichard and Sarma [1981] J.-L. Pichard and G. Sarma, J. Phys. C: Solid State Phys. C14, L127 (1981).
  * Slevin and Ohtsuki [1999] K. Slevin and T. Ohtsuki, Phys. Rev. Lett. 82, 382 (1999).
  * Wigner [1951] E. P. Wigner, Mathematical Proceedings of the Cambridge Philosophical Society 47, 790–798 (1951).
  * Dyson [1962a] F. J. Dyson, Journal of Mathematical Physics 3, 140 (1962a).
  * Dyson [1962b] F. J. Dyson, Journal of Mathematical Physics 3, 1199 (1962b).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Physical Review B 55, 1142 (1997).
  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annals of Physics 321, 1126 (2006).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Billy  _et al._ [2008] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Nature 453, 891 (2008).
  * Modugno [2010] G. Modugno, Reports on Progress in Physics 73, 102401 (2010).
  * Sanchez-Palencia  _et al._ [2008] L. Sanchez-Palencia, D. Clément, P. Lugan, P. Bouyer, and A. Aspect, New Journal of Physics 10 (2008), 10.1088/1367-2630/10/4/045019.
  * John [1987] S. John, Phys. Rev. Lett. 58, 2486 (1987).
  * Wiersma  _et al._ [1997] D. S. Wiersma, P. Bartolini, A. Lagendijk, and R. Righini, Nature 390, 671 (1997).
  * Segev  _et al._ [2013] M. Segev, Y. Silberberg, and D. N. Christodoulides, Nature Photonics 7, 197 (2013).
  * Mafi [2015] A. Mafi, Advances in Optics and Photonics 7, 459 (2015).
  * Skipetrov and Sokolov [2015] S. E. Skipetrov and I. M. Sokolov, Phys. Rev. Lett. 114, 053902 (2015).
  * Skipetrov [2018] S. E. Skipetrov, Phys. Rev. Lett. 121, 093601 (2018).
  * Kirkpatrick [1985] T. R. Kirkpatrick, Phys. Rev. B 31, 5746 (1985).
  * Weaver [1990] R. Weaver, Wave Motion 12, 129 (1990).
  * Hu  _et al._ [2008] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nature Physics 4, 945 (2008).
  * Slevin and Ohtsuki [2016] K. Slevin and T. Ohtsuki, Journal of the Physical Society of Japan 85, 104712 (2016).
  * Slevin and Ohtsuki [2018] K. Slevin and T. Ohtsuki, Journal of the Physical Society of Japan 87, 094703 (2018), https://doi.org/10.7566/JPSJ.87.094703 .
  * Asada  _et al._ [2005] Y. Asada, K. Slevin, and T. Ohtsuki, J. Phys. Soc. Jpn. 74, 238 (2005).
  * Roy  _et al._ [2017] B. Roy, Y. Alavirad, and J. D. Sau, Phys. Rev. Lett. 118, 227002 (2017).
  * Heinzner  _et al._ [2005] P. Heinzner, A. Huckleberry, and M. R. Zirnbauer, Communications in Mathematical Physics 257, 725 (2005).
  * Kawabata and Ryu [2021] K. Kawabata and S. Ryu, Phys. Rev. Lett. 126, 166801 (2021).
  * Xu  _et al._ [2016] B. Xu, T. Ohtsuki, and R. Shindou, Phys. Rev. B 94, 220403 (2016).
  * Tzortzakakis  _et al._ [2020] A. F. Tzortzakakis, K. G. Makris, and E. N. Economou, Phys. Rev. B 101, 014202 (2020).
  * Wang and Wang [2020] C. Wang and X. R. Wang, Phys. Rev. B 101, 165114 (2020).
  * Huang and Shklovskii [2020a] Y. Huang and B. I. Shklovskii, Phys. Rev. B 101, 014204 (2020a).
  * Huang and Shklovskii [2020b] Y. Huang and B. I. Shklovskii, Phys. Rev. B 102, 064212 (2020b).
  * Luo  _et al._ [2021a] X. Luo, T. Ohtsuki, and R. Shindou, Phys. Rev. Lett. 126, 090402 (2021a).
  * Luo  _et al._ [2021b] X. Luo, T. Ohtsuki, and R. Shindou, “Transfer matrix study of the anderson transition in non-hermitian systems,” (2021b), arXiv:2103.05239 [quant-ph] .
  * Cao  _et al._ [1999] H. Cao, Y. G. Zhao, S. T. Ho, E. W. Seelig, Q. H. Wang, and R. P. H. Chang, Phys. Rev. Lett. 82, 2278 (1999).
  * Luo  _et al._ [2020] X. Luo, B. Xu, T. Ohtsuki, and R. Shindou, Phys. Rev. B 101, 020202 (2020).
  * Ortuño  _et al._ [2009] M. Ortuño, A. M. Somoza, and J. T. Chalker, Phys. Rev. Lett. 102, 070603 (2009).
  * Mondragon-Shem  _et al._ [2014] I. Mondragon-Shem, T. L. Hughes, J. Song, and E. Prodan, Phys. Rev. Lett. 113, 046802 (2014).
  * Song and Prodan [2014] J. Song and E. Prodan, Phys. Rev. B 89, 224203 (2014).
  * Li  _et al._ [2017] L. Li, C. Yin, S. Chen, and M. A. N. Araújo, Phys. Rev. B 95, 121107 (2017).
  * Sbierski  _et al._ [2020] B. Sbierski, J. F. Karcher, and M. S. Foster, Phys. Rev. X 10, 021025 (2020).
  * Claes and Hughes [2020] J. Claes and T. L. Hughes, Phys. Rev. B 101, 224201 (2020).
  * Luke  _et al._ [1998] G. M. Luke, Y. Fudamoto, K. M. Kojima, M. I. Larkin, J. Merrin, B. Nachumi, Y. J. Uemura, Y. Maeno, Z. Q. Mao, Y. Mori, H. Nakamura, and M. Sigrist, Nature 394, 558 (1998).
  * Xia  _et al._ [2006] J. Xia, Y. Maeno, P. T. Beyersdorf, M. M. Fejer, and A. Kapitulnik, Phys. Rev. Lett. 97, 167002 (2006).
  * Schemm  _et al._ [2014] E. R. Schemm, W. J. Gannon, C. M. Wishne, W. P. Halperin, and A. Kapitulnik, Science 345, 190 (2014).
  * Liu  _et al._ [2016] S. Liu, T. Ohtsuki, and R. Shindou, Phys. Rev. Lett. 116, 066401 (2016).
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Phys. Rev. Lett. 107, 127205 (2011).
  * Slevin and Ohtsuki [2014] K. Slevin and T. Ohtsuki, New J. Phys. 16, 015012 (2014).
  * Bevington and Robinson [2003] P. Bevington and D. Robinson,  _Data Reduction and Error Analysis for the Physical Sciences_ (McGraw-Hill Education, 2003).
  * Kagalovsky  _et al._ [2004] V. Kagalovsky, B. Horovitz, and Y. Avishai, Phys. Rev. Lett. 93, 246802 (2004).
  * Verbaarschot and Zahed [1993] J. J. M. Verbaarschot and I. Zahed, Phys. Rev. Lett. 70, 3852 (1993).
  * Verbaarschot and Zahed [1994] J. J. M. Verbaarschot and I. Zahed, Phys. Rev. Lett. 73, 2288 (1994).
  * Slevin and Nagao [1993] K. Slevin and T. Nagao, Phys. Rev. Lett. 70, 635 (1993).
  * Wegner [1989] F. Wegner, Nuclear Physics B 316, 663 (1989).
  * König  _et al._ [2012] E. J. König, P. M. Ostrovsky, I. V. Protopopov, and A. D. Mirlin, Phys. Rev. B 85, 195130 (2012).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).


program

# Repetitive Transients Extraction Algorithm
for Detecting Bearing Faults111Preprint

Wangpeng He  School of Aerospace Science and Technology, Xidian University,
Xi’an, China  Yin Ding email<EMAIL_ADDRESS>Tandon School of Engineering, New
York University, 6 Metrotech Center, Brooklyn, NY, USA  Yanyang Zi  Ivan W.
Selesnick  Tandon School of Engineering, New York University, 6 Metrotech
Center, Brooklyn, NY, USA

###### Abstract

This paper addresses the problem of noise reduction with simultaneous
components extraction in vibration signals for faults diagnosis of bearing.
The observed vibration signal is modeled as a summation of two components
contaminated by noise, and each component composes of repetitive transients.
To extract the two components simultaneously, an approach by solving an
optimization problem is proposed in this paper. The problem adopts convex
sparsity-based regularization scheme for decomposition, and non-convex
regularization is used to further promote the sparsity but preserving the
global convexity. A synthetic example is presented to illustrate the
performance of the proposed approach for repetitive feature extraction. The
performance and effectiveness of the proposed method are further demonstrated
by applying to compound faults and single fault diagnosis of a locomotive
bearing. The results show the proposed approach can effectively extract the
features of outer and inner race defects.

##  1 Introduction

Rolling bearings are one of the most prevalent components in rotating machines
and reciprocating machines [35]. Vibration-based fault detection has become
the preferred technique for bearing fault diagnosis [44, 18]. Bearing
vibrations are usually modeled as cyclostationary signals (or group-sparse
periodic) [9, 10, 4]. When a localized defect occurs on the bearing,
repetitive transients will be generated due to the passing of rollers over the
defect [26, 38]. These transients have repetitive structure and are usually
submerged in background noise. Many denoising methods have been introduced to
extract fault features for the purpose of detecting faults in machines, such
as wavelet transform [43, 12], singular value decomposition (SVD) [27], time-
frequency analysis methods [21, 23], empirical mode decomposition (EMD) [30],
methods for blind extraction of a cyclostationary signal and spectral kurtosis
(SK) [9, 10, 5]. If compound faults exist, then the observed vibration signals
are rather complex and it is difficult to identify each fault using
traditional signal processing methods.

A number of approaches have been developed for the multiple fault or compound
fault diagnosis. Principal component analysis (PCA) based method has been used
to classify faults [29, 31, 33]. Support vector machine (SVM) based methods
have been used for multi-fault diagnosis [46, 1] as well. Some other
techniques such as neural network and independent component analysis (ICA),
have been introduced to assist fault detection and classification [42, 8, 7].
Some of these methods require large data collection as a training set and
further require an off-line training phase; and some methods rely on careful
selection of features (e.g., wavelet packet sub-bands).

Adopting sparsity in the field of fault detection was initially illustrated in
Ref. [45], where basis pursuit denoising (BPD) [14] was used to exploit sparse
features in various domains to detect faults. Some recent works consider other
sparse representations for fault diagnosis [44, 16, 25, 17]. The fault
frequencies of potential bearing fault features can be simply obtained using
the geometry of the components in many cases or directly obtained from the
user operation manual. Many works have considered the fault frequencies as
priori information [44, 32, 40, 34, 24].

In this work, a method using the temporal periodicity (namely fault
characteristic or fundamental frequency) directly in the time domain is
proposed to extract fault features while simultaneously denoising. The
proposed method is based on convex optimization using non-convex
regularization. Specifically, this paper aims to address the problem of
extracting compound features caused by faults in vibration signals, where the
features exhibit repetitive group sparsity.

In particular, the observed signal is modeled as

$\displaystyle y=x_{1}+x_{2}+w,$ (1)

where $w$ denotes additive white Gaussian noise (AWGN), and $x_{1}$ and
$x_{2}$ are both repetitively group-sparse signals with periods $T_{1}$ and
$T_{2}$ respectively. Note that here the repetitively group-sparse signal
means the useful features (group-sparse) appear repetitively. In other words,
the features have group-sparse property and appear repetitively. To avoid
confusion, here the “period” is referred to as the cycle for cyclostationary
signal in Ref. [2]. In the case of compound faults detection of bearings,
useful features $x_{1}$ and $x_{2}$ also satisfy the following two conditions.

  1. 1.

The periods $T_{1}$ and $T_{2}$ are different.

  2. 2.

Each period is not close to an integer multiple of the other.

Note that, in real applications, the fault features are not strictly group-
sparse periodic. In other words, the real fault frequencies have a variation
from the calculated frequencies (up to 1-2%), i.e., there are not strict
periods $T_{1}$ and $T_{2}$ [35]. In such cases, the proposed approach can
still work owing to the overlapping of group structure, which will be shown in
the following sections. The proposed approach is not an improvement of the
cutting edge techniques developed in the cyclostationary framework, instead,
it is an alternative to these techniques. A work closely related to the
proposed approach is the periodic overlapping group sparsity (POGS) problem
[24], which assumes only one repetitive group-sparse component is present in
the vibration signal. Notice that although the observed signal is modeled with
two components, the proposed method also works when there exist only one fault
component. Therefore, the proposed method generalizes POGS, and is useful when
there exists multiple components.

Another related work is group-sparse signal denoising (GSSD), which is also
known as overlapping group sparsity (OGS) with non-convex regularization [13],
where mathematical derivations and proofs have been given in detail to show
that non-convex regularization can be used to promote group-sparsity, while
maintaining convexity of the problem as a whole.

The method proposed in this paper also uses the concept of morphological
component analysis (MCA) [39], which is a method to decompose signals based on
sparse representations. In contrast to MCA, the proposed method does not
utilize any transform (e.g., Fourier or wavelet transform), i.e., the sparse
features are in the signal domain. Moreover, non-convex regularization is used
to strongly induce sparsity while maintaining convexity of the proposed
problem formulation. The proposed approach reduces to the OGS method without
prior knowledge, i.e., we can utilize the sparsity-based OGS approach if we do
not have prior knowledge of the characteristic frequencies.

##  2 Preliminaries

###  2.1 Notation

In this paper, the elements of a vector $x$ are denoted as $x_{n}$ or
$[x]_{n}$. The norms of $x$ are defined as

$\displaystyle\lVert x\rVert_{1}:=\sum_{n}\left|x_{n}\right|,\quad\lVert
x\rVert_{2}:=\bigg{(}\sum_{n}\left|x_{n}\right|^{2}\bigg{)}^{1/2}.$ (2)

A function of $x$ determined by parameter $a$ is denoted as $f(x;a)$, and to
distinguish it from a function with two ordered arguments, e.g., $f(x,y)$.

###  2.2 Review of majorization-minimization

In this paper, the majorization-minimization (MM) approach is used to derive a
fast-converging algorithm. This subsection briefly describes the MM approach
for minimizing a convex cost function. The MM is an approach to simplify a
complicated optimization problem into a sequence of simpler ones [22]. More
specifically, consider an optimization problem

$\displaystyle u^{\mathsf{opt}}=\arg\min_{u}F(u).$ (3)

Using MM, the problem can be solved iteratively by

$\displaystyle u^{(i+1)}=\arg\min_{u}F^{\mathsf{M}}(u,u^{(i)}),$ (4)

where $F^{\mathsf{M}}:{}^{N}\times{}^{N}\to\real$ is an upper bound
(majorizer) of the objective function $F$, satisfying

$\displaystyle F^{\mathsf{M}}(u,v)\geq F(u),\quad F^{\mathsf{M}}(u,u)=F(u).$
(5)

Note that the majorizer $F^{\mathsf{M}}(u,v)$ touches $F(u)$ for $u=v$, as
shown in (5). Figure 1 illustrates the majorizer (red line) of a penalty
function (blue line), which will be described in the following subsection. The
proof of convergence for MM has been given in Ref. [28, Chapter 10]. More
details about the MM procedure can be found in [28, 22] and references
therein.

###  2.3 Non-convex penalty functions

Non-convex penalty functions can promote sparsity more strongly than convex
penalty functions [36, 37]. This subsection briefly describes the non-convex
penalty functions which will be used in the proposed approach. The smoothed
non-convex penalty function $\phi_{\epsilon}:\real\to{}_{+}$ is used in this
work. Table 1 gives several examples of the functions defined by

$\displaystyle\phi_{\epsilon}(u;a):=\phi(\sqrt{u^{2}+\epsilon};a),\quad\epsilon>0,$
(6)

where $\phi$ is a non-smooth penalty function satisfying the following
properties:

  1. 1.

$\phi(u;a)$ is continuous on .

  2. 2.

$\phi(u;a)$ is twice continuously differentiable on
$\mathbb{R}\backslash\\{0\\}$.

  3. 3.

$\phi(u;a)$ is even symmetric: $\phi(u)=\phi(-u)$.

  4. 4.

$\phi(u;a)$ is increasing and concave on +.

  5. 5.

$\phi(u;a)=\left|x\right|$ when $a=0$.

Note that, for both $\phi$ and $\phi_{\epsilon}$, the parameter $a\geq 0$
controls the concavity of the function.

The parameter $\epsilon$ controls the smoothness of the functions. As a
special case, when $\epsilon=0$, $\phi_{\epsilon}(u;0)=\phi(u)$, then the
penalty function is non-differentiable at 0. In practice, $\epsilon$ is
specified very small, e.g. $10^{-10}$, so that the function is differentiable.
Figure 1 gives two specific examples of $\phi_{\epsilon}$.

Figure 1:  (a) Smoothed $\ell_{1}$-norm (abs) penalty function (blue) and
majorizer (red). (b) Smoothed arctangent (atan) penalty function and majorizer
(gray). For all the figures, the parameters are set to
$v=0.5,\epsilon=0.02,a=0.5$.

A majorizer of penalty function $\phi_{\epsilon}$ is given by

$\displaystyle\phi_{\epsilon}^{\mathsf{M}}(u,v;a):=\frac{u^{2}}{2\psi(v;a)}-\underbrace{\bigg{(}\frac{v^{2}}{2\psi(v;a)}-\phi_{\epsilon}(v;a)\bigg{)}}_{\text{only
depends on $v$}},$ (7)

where $\psi(v;a)$ is listed in the third column of Table 1. Note that
$\psi(v;a)>0$ for all $v\in\real$. Also note that
$\phi_{\epsilon}^{\mathsf{M}}(u,v)$ is quadratic in $u$. In Ref. [19], a
detailed proof has been given to show that when $\phi$ satisfies the above
properties, the majorizer (7) satisfies the condition (5) for
$\phi_{\epsilon}$.

###  2.4 Review of POGS

The periodic overlapping group sparsity (POGS) problem [24] considers the
signal model

$\displaystyle y=x+w,$ (8)

where $x$ is a repetitive group-sparse signal, and $w$ is additive white
Gaussian noise. The POGS method extracts $x$ by solving a convex problem

$\displaystyle
x^{\mathsf{opt}}=\arg\min_{x}\Big{\\{}P_{0}(x)=\frac{1}{2}\lVert
y-x\rVert_{2}^{2}+\lambda\Phi(x,b;a)\Big{\\}},$ (9)

where $\Phi:{}^{N}\to\real$ is defined as

$\displaystyle\Phi(x,b;a):=\sum_{n}\phi_{\epsilon}\Big{(}\Big{[}\sum_{k}[b]_{k}[x]_{n+k}^{2}\Big{]}^{1/2};a\Big{)}.$
(10)

The function $\Phi$ is a regularization term that promotes repetitive group-
sparsity. Moreover, $b\in\\{0,1\\}^{K}$ is a binary weight array designated
according to the period. To simplify following derivation, denote that

$\displaystyle\Phi(x,b;0)=\Phi(x;b),$ (11)

where the penalty function is strictly convex as shown in Figure 1(a).

Table 1: Sparsity-promoting penalty functions.

Penalty | $\phi(u;a)$ | $\phi_{\epsilon}(u;a)$ | $\psi(u;a)=u/\phi_{\epsilon}^{\prime}(u;a)\quad$
---|---|---|---
abs ($a=0$) | $\left|u\right|$ | $\sqrt{u^{2}+\epsilon}$ | $\sqrt{u^{2}+\epsilon}$
log | $\displaystyle\frac{1}{a}\log(1+a|u|)$ | $\displaystyle\frac{1}{a}\log(1+a\sqrt{u^{2}+\epsilon})$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}\right)$
rat | $\displaystyle\frac{\left|u\right|}{1+a\left|u\right|/2}$ | $\displaystyle\frac{\sqrt{u^{2}+\epsilon}}{1+a\sqrt{u^{2}+\epsilon}/2}$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}/2\right)^{2}$
atan | $\displaystyle\frac{2}{a\sqrt{3}}\left(\tan^{-1}\left(\frac{1+2a|u|}{\sqrt{3}}\right)-\frac{\pi}{6}\right)$ | $\displaystyle\frac{2}{a\sqrt{3}}\left(\tan^{-1}\left(\frac{1+2a\sqrt{u^{2}+\epsilon}}{\sqrt{3}}\right)-\frac{\pi}{6}\right)$ | $\displaystyle\sqrt{u^{2}+\epsilon}\left(1+a\sqrt{u^{2}+\epsilon}+a^{2}(u^{2}+\epsilon)\right)$

##  3 Repetitive transients extraction algorithm

In this section, the proposed algorithm termed repetitive transients
extraction algorithm (RTEA) is presented.

###  3.1 Problem formulation

To extract two repetitive group-sparse components, an optimization problem is
formulated as

$\displaystyle\\{x_{1}^{\mathsf{opt}},x_{2}^{\mathsf{opt}}\\}$
$\displaystyle=\arg\min_{x_{1},x_{2}}\Big{\\{}P(x_{1},x_{2})=\frac{1}{2}\lVert
y-(x_{1}+x_{2})\rVert_{2}^{2}+\lambda_{0}R(x_{1},x_{2};a_{0})+\sum_{i\in\\{1,2\\}}\lambda_{i}\Phi(x_{i};b_{i})\Big{\\}}$
(12)

where $R:{}^{N}\times{}^{N}\to\real$ is defined as

$\displaystyle R(x_{1},x_{2};a_{0})$
$\displaystyle:=\sum_{n}\phi_{\epsilon}\Big{(}\Big{[}\sum_{k=0}^{K_{0}-1}[x_{1}+x_{2}]_{n+k}^{2}\Big{]}^{1/2};a_{0}\Big{)}.$
(13)

The function $R$ is an overlapping group sparsity (OGS) regularization
function with group size $K_{0}$222Using the notation of [24], it is denoted
$\sum_{i}\phi_{\epsilon}(\lVert x\rVert_{i,K_{0}};a_{0})$. . There are two
more regularizers in (12) promoting the repetitive group-sparsity of $x_{1}$
and $x_{2}$ respectively, and formulated in (11).

Furthermore, in problem (12), $b_{1}\in\\{0,1\\}^{K_{1}}$ and
$b_{2}\in\\{0,1\\}^{K_{2}}$ are two binary-weighting arrays as

$\displaystyle b_{i}=[\ \underbrace{\underbrace{1~{}1\cdots 1}_{N_{i1}}\
\underbrace{0~{}0\cdots 0}_{N_{i0}}\ \dots\ \underbrace{1~{}1\cdots
1}_{N_{i1}}\ \underbrace{0~{}0\cdots 0}_{N_{i0}}\ \underbrace{1~{}1\cdots
1}_{N_{i1}}}_{\text{spanning }M_{i}\text{ periods}}\ ],$ (14)

for $i=1,2$ and $M_{1}$ and $M_{2}$ defines the number of periods included in
$b_{1}$ and $b_{2}$ respectively.

Moreover, in contrast to MCA, which has two regularizers, problem (12) has
three. The regularization term $R$ is introduced because according to the
signal model, the summation of the two components is also sparse. Note that,
$(x_{1}+x_{2})$ might be sparse when $x_{1}$ and $x_{2}$ are not, but in (12),
the regularizers with $\Phi$ does force $x_{1}$ and $x_{2}$ to be sparse.

###  3.2 Convexity of the objective function

The idea of using non-covex regularization with “maximizing concavity” in a
convex problem has been illustrated in [36], wherein the quadratic formulation
of data fidelity term can be used to compensate the non-convexity in the
regularization so that the objective function is still convex. In this work,
we adopt this idea for an extraction algorithm with simultaneous denoising.
Entirely there are three regularizers in (12), and we allow one of them to be
non-convex to promote the global sparsity more strongly. Moreover, the
specific condition to assure the convexity of problem (12) is derived as the
following proposition.

######  Proposition 1.

Suppose the parameterized penalty function $\phi_{\epsilon}$ is defined by
formula (6) and $\lambda_{0}>0$. If

$0\leq a_{0}<\frac{1}{K_{0}\lambda_{0}},$ (15)

then the objective function $P:{}^{N}\times{}^{N}\to\real$ in (12) is strictly
convex.

A proof of the above proposition is given in Appendix A.

##  4 Algorithm derivation

In this section, an algorithm is derived to solve (12) based on MM. The
majorizer of $R$ in (13) is
$R^{\mathsf{M}}:{}^{N}\times{}^{N}\times{}^{N}\times{}^{N}\to\real$, and
written explicitly as

$\displaystyle R$ ${}^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2};a_{0})$
$\displaystyle=\sum_{n}\Big{\\{}r_{0}(n,z_{1}+z_{2})([x_{1}]_{n}^{2}+[x_{2}]_{n}^{2}-[z_{1}-z_{2}]_{n}[x_{1}]_{n}-[z_{2}-z_{1}]_{n}[x_{2}]_{n})\Big{\\}}+C(z_{1},z_{2}),$
(16)

where $C(z_{1},z_{2})$ is a constant only dependent on $z_{1}$ and $z_{2}$. In
(4), $r_{0}$ is a function $r_{0}:\mathbb{Z}\times{}^{N}\to\real$ given by

$\displaystyle
r_{0}(n,z_{1}+z_{2})=\sum_{j=0}^{K_{0}-1}{\psi^{-1}\Big{(}\Big{[}\sum_{k=0}^{K_{0}-1}[z_{1}+z_{2}]_{n-j+k}^{2}\Big{]}^{1/2};a_{0}\Big{)}}.$
(17)

The majorizer of function $\Phi$ has been derived in Ref. [24, Section 3.3].
Here, it can be rewritten using a slightly different notation, that
$\Phi^{\mathsf{M}}:{}^{N}\times{}^{N}\to\real$ is

$\displaystyle\Phi^{\mathsf{M}}(x,z;b,a)=\frac{1}{2}\sum_{n}r(n,z)x_{n}^{2}+C(z),$
(18)

where $r:\mathbb{Z}\times{}^{N}\to\real$ is defined as

$\displaystyle
r(n,z):=\sum_{j=0}^{K-1}\frac{[b]_{j}}{\psi\Big{(}\Big{[}\sum_{k}[b]_{k}[z]_{n-j+k}^{2}\Big{]}^{1/2};a\Big{)}}$
(19)

Using the above results, the majorizer of the objective function $P$ in (12)
can be written as

$\displaystyle P^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2})=$
$\displaystyle\frac{1}{2}\lVert
y-(x_{1}+x_{2})\rVert_{2}^{2}+\frac{1}{2}\lVert(x_{1}-z_{1})-(x_{2}-z_{2})\rVert_{2}^{2}$
$\displaystyle\quad+\lambda_{0}R^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2};a_{0})+\\!\sum_{i\in\\{1,2\\}}\\!\lambda_{i}\Phi^{\mathsf{M}}(x_{i},z_{i};b_{i},a_{i}),$
(20)

where the function
$R^{\mathsf{M}}:{}^{N}\times{}^{N}\times{}^{N}\times{}^{N}\to\real$ is a
majorizer of $R$ in (13), wherein the derivation in detail is given in
Appendix B.

Table 2: Explicit steps of proposed method (RTEA).
$\displaystyle\text{Input:}~{}y\in{}^{N},~{}K_{0},~{}\lambda_{0},~{}\lambda_{i},~{}b_{i}\in\\{0,1\\}^{K_{i}},\text{
for }i\in\\{1,2\\}.$ $\displaystyle\text{Initialization: }a_{0}\geq
0,~{}a_{i}\geq 0,~{}x_{i}=y,\text{ for }i=1,2.$ Repeat:
$\displaystyle\quad[r_{0}]_{n}=\sum_{j=0}^{K_{0}-1}\frac{1}{\displaystyle\psi_{0}\Big{(}\Big{[}\sum_{k=0}^{K_{0}-1}[x_{1}+x_{2}]_{n-j+k}^{2}\Big{]}^{1/2};a_{0}\Big{)}}$
$\displaystyle\quad[r_{i}]_{n}=\sum_{j=0}^{K_{i}-1}\frac{[b_{i}]_{j}}{\displaystyle\psi_{i}\Big{(}\Big{[}\sum_{k=0}^{K_{i}-1}[b_{i}]_{k}[x_{i}]_{n-j+k}^{2}\Big{]}^{1/2};a_{i}\Big{)}},\text{
for }i\in\\{1,2\\}$
$\displaystyle\quad[p_{i}]_{n}=2+2\lambda_{0}[r_{0}]_{n}+\lambda_{i}[r_{i}]_{n},\text{
for }i\in\\{1,2\\}$
$\displaystyle\quad[q_{1}]_{n}=y_{n}+(1+[r_{0}]_{n})[x_{2}-x_{1}]_{n}$
$\displaystyle\quad[q_{2}]_{n}=y_{n}+(1+[r_{0}]_{n})[x_{1}-x_{2}]_{n}$
$\displaystyle\quad[x_{i}]_{n}=[q_{i}]_{n}/[p_{i}]_{n},\text{ for
}i\in\\{1,2\\}$ Until convergence $\displaystyle\text{Return: }x_{1},x_{2}$

Note that in (4), an extra term
$\frac{1}{2}\lVert(x_{1}-z_{1})-(x_{2}-z_{2})\rVert_{2}^{2}$ is introduced,
which does not break the property of majorizer, but helps to cancel the term
$x_{1}^{\mathsf{T}}x_{2}$ in the data fidelity term. Then
$P^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2})$ can be written as

$\displaystyle P^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2})=$
$\displaystyle\sum_{n}\frac{1}{2}[p_{1}(z_{1},z_{2})\odot
x_{1}]_{n}^{2}+\frac{1}{2}[p_{2}(z_{1},z_{2})\odot x_{2}]_{n}^{2}$
$\displaystyle-[q_{1}(z_{1},z_{2})\odot x_{1}]_{n}-[q_{2}(z_{1},z_{2})\odot
x_{2}]_{n}+C(z_{1},z_{2}),$ (21)

where $\odot$ denotes element-wise multiplication, and the coefficients
$p_{1},p_{2},q_{1},q_{2}\in{}^{N}$ can be written explicitly as

$\displaystyle[p_{1}(z_{1},z_{2})]_{n}$
$\displaystyle=2+2\lambda_{0}r_{0}(n,z_{1}+z_{2})+\lambda_{1}r_{1}(n,z_{1})$
(22a) $\displaystyle[p_{2}(z_{1},z_{2})]_{n}$
$\displaystyle=2+2\lambda_{0}r_{0}(n,z_{1}+z_{2})+\lambda_{2}r_{2}(n,z_{2})$
(22b) $\displaystyle[q_{1}(z_{1},z_{2})]_{n}$
$\displaystyle=y_{n}+[1+r_{0}(n,z_{1}+z_{2})][z_{1}-z_{2}]_{n}$ (22c)
$\displaystyle[q_{2}(z_{1},z_{2})]_{n}$
$\displaystyle=y_{n}+[1+r_{0}(n,z_{1}+z_{2})][z_{2}-z_{1}]_{n}.$ (22d)

Note that $r_{1}$ and $r_{2}$ are functions defined by (19) dependent on the
binary weighting arrays $b_{1}$ and $b_{2}$ in (14) respectively.

Consequently, using MM, the problem (12) is solved iteratively by

$\displaystyle\\{x_{1}^{(i+1)},x_{2}^{(i+1)}\\}=\arg\min_{x_{1},x_{2}}P^{\mathsf{M}}(x_{1},x_{2},x_{1}^{(i)},x_{2}^{(i)}),$
(23)

and the results of $x_{1}$ and $x_{2}$ in each iteration can be written
separately as

$\displaystyle[x_{1}^{(i+1)}]_{n}$
$\displaystyle=\frac{[q_{1}(x_{1}^{(i)},x_{2}^{(i)})]_{n}}{[p_{1}(x_{1}^{(i)},x_{2}^{(i)})]_{n}},$
(24a) $\displaystyle[x_{2}^{(i+1)}]_{n}$
$\displaystyle=\frac{[q_{2}(x_{1}^{(i)},x_{2}^{(i)})]_{n}}{[p_{2}(x_{1}^{(i)},x_{2}^{(i)})]_{n}}.$
(24b)

Table 2 gives the specific steps of the proposed algorithm RTEA. Note that, in
this algorithm, the computation of $r_{0}$, $r_{1}$ and $r_{2}$ can be
implemented directly with convolution, and the rest steps are element-wised
independent, where practically multi-threading computation can be directly
adopted. The systematic structure of the proposed RTEA for fault detection of
rolling element bearings is presented in Figure 2. Since this algorithm is
derived using MM procedure, the convergence can be guaranteed and converges to
the optimal minimizer in this case. The detailed proof of MM for convex
problems can be found in Ref. [28, Chapter 10].

Figure 2: Procedure of the proposed RTEA for fault detection of bearings.
Figure 3: Example 1: Test signal. (a) Component 1 with period $T_{1}=32$
samples. (b) Component 2 with period $T_{2}=53$ samples. (c) Summation of the
two compoennts. (d) Noisy observation.  Figure 4: Example 1: Results of
proposed method. Figure 5: Example 1: Cost function history.

##  5 Synthetic data example

Figure 3 shows the synthetic test signal. It consists of two sparse transient
sequences exhibiting distinct periods, where $T_{1}=32$ and $T_{2}=53$
samples. In this example, each transient is generated by adding a random
number of sinusoids with random frequencies and initial phases. More
specifically, each transient can be written as

$\displaystyle g(n)=\sum_{j=1}^{J}A_{j}\sin(\omega_{j}n+\theta_{j}),\quad
n\in\\{0,1,2\dots 9\\},$ (25)

where $1\leq J\leq 10$ is a random integer, and for each $j$, $A_{j}$ is a
random amplitude, $\omega_{j}$ is a random frequency, and $\theta_{j}$ is a
random phase. The sequences of transients are shown in Figure 3(a) and Figure
3(b), and their summation is shown in Figure 3(c).

In this example, the parameters are set to $N_{11}=N_{21}=K_{0}=3$ and $M=4$,
to determine the binary arrays $b_{1}$ and $b_{2}$ by the priori known
periods. The root-mean-square error (RMSE) is used as an evaluation metric.
Moreover, the penalty functions are selected by Proposition 1 assuring the
problem is strictly convex. Figure 4 shows the results from the proposed
method, and Figure 5 shows the convergence behavior of the proposed algorithm.

Figure 6: Example 1: Results by using (26) with convex regularizations. Figure
7: Example 1: Results by using some popular and/or state-of-the-art methods.
(a) Undecimated wavelet based denoising. (b) Fused lasso with non-convex
regularization. (c) group-sparse signal denoising (GSSD) with group size 3.

As a comparison, considering another problem, based strictly on MCA, whose
objective function is

$\displaystyle\arg\min_{x_{1},x_{2}}\frac{1}{2}\lVert
y-(x_{1}+x_{2})\rVert_{2}^{2}+\sum_{i\in\\{1,2\\}}\lambda_{i}\Phi(x_{i};b_{i}).$
(26)

Note that (26) has one regularization term less than (12), where the global
sparsity is not promoted by function $R$.

Although non-convex penalties can help to promote sparsity, this will break
the convexity of problem (26), so that a global optimal solution is not
assured. Moreover, experiment results find that convex formulation (12)
obtains as good a performance as the non-convex formulation (26). As a
consequence, the formulation (12) is preferred, because there is no need to
sacrifice convexity.

In addition, the proposed method with non-convex formulation is also
evaluated, where $a_{1}$ and/or $a_{2}$ are greater than zero. In this case,
losing the convexity, the resulting sparsity can be even further promoted,
then the formulation (12) is still preferred to the MCA formulation (26).

Further comparisons to some denoising methods are also presented. Firstly,
wavelet-based denoising method is adopted to the test signal, more
specifically, a 6-scale undecimated wavelet transform [15] using Haar wavelet
filter. For denoising, hard-thresholding is applied and the threshold value is
chosen by $3\sigma$-rule for each subband. As shown in Figure 7(a), the
denoising results adheres the shape of impulse response of the wavelet, where
some of the signal does not have a zero baseline.

Secondly, the non-convex regularized fused lasso (FL) approach proposed in [6]
is adopted, which is an improved version of conventional FL [41], allowing one
regularizer to be non-convex and preserving the global convexity. The result
is shown in Figure 7(b). Thirdly, group-sparse signal denoising (GSSD) also
known as OGS with non-convex regularization [13] is adopted. More
specifically, ‘atan’ penalty with group size to be 3 samples is used, and the
regularization parameter is chosen to optimize the RMSE. The result is shown
in Figure 7(c). Both of the above methods also have worse recoveries in terms
of RMSE. In addition, these methods are all denoising methods only, they
cannot decompose the signal into two distinct sequences of transients, each
exhibiting its own period.

###  5.1 Parameter selection

Table 3: Selection of $\beta_{i}$ for $i=0,1,2$. $M_{i}$ $N_{i1}$ | 1 | 2 | 3 | 4
---|---|---|---|---
1 $(K_{0})$ | 3.700 | 1.700 | 1.150 | 0.925
2 | 1.700 | 0.850 | 0.625 | 0.475
3 | 1.150 | 0.625 | 0.450 | 0.375
4 | 0.925 | 0.475 | 0.375 | 0.325

In Section 4.1 of [24], the schemes to set the binary weight array $b$ (14)
and the regularization parameter $\lambda$ for POGS problem have been
discussed in detail. Moreover, a look-up table has been given as a guide to
choose the regularization parameter, where using the given multiplier in [24,
Table 3], the regularization parameter can be chosen by
$\lambda=\beta\sigma_{w}$, where $\sigma_{w}$ is the deviation of the additive
noise. Here, the table is quoted with a slight change of notation.

Figure 8: Example 1: RMSE values as a function of $\eta$ using (27).

As an extension of POGS with simultaneously decomposition, this parameter
selection scheme can be used with a modification wherein the weights are
shared among the three regularizers. In this case, the regularization
parameters $\\{\lambda_{0},\lambda_{1},\lambda_{2}\\}$ are suggested to be set
by

$\displaystyle\lambda_{0}$ $\displaystyle=\eta\beta_{0}$ (27a)
$\displaystyle\lambda_{i}$ $\displaystyle=0.5(1-\eta)\beta_{i}~{}\text{ for
}i=1,2$ (27b)

where $0<\eta<1$ is a parameter to balance sparsity of the sum $(x_{1}+x_{2})$
and sparsity of the individual signal component $x_{i}$.

To run the method, firstly the binary weight arrays $b_{1},b_{2}$ need to be
set by (14) using the priori periods. Then $\beta_{i}$ can be determined by
Table 3 using $N_{i1}$ and $M_{i}$ for $i=1,2$. Note that it is necessary to
chose $K_{0}=\min\\{N_{11},N_{21}\\}$ to induce group-sparsity. Hence in
practice $K_{0}$ is not necessary to be chosen, and so as $\beta_{0}$.

In practice, the parameters are set to be $K_{0}=N_{11}=N_{21}$, and
$M_{1}=M_{2}=4$. The regularization parameters can be straight-forwardly
determined by Table 3 and (27).

Setting parameter $\eta$. The parameter $0<\eta<1$ in (27) is to balance the
sparsity, as mentioned above. As special cases:

  1. 1.

If $\eta\to 1$, then $\lambda_{1}=\lambda_{2}\to 0$, and the approach promotes
sparsity of $x_{1}+x_{2}$, but leads to $x_{1}=x_{2}$.

  2. 2.

If $\eta\to 0$, then the problem (12) reduces to (26), which is a conventional
MCA problem. It is able to seperate the components, but promotes sparsity
weakly.

To obtain both of the benefits, through experiments, parameter $\eta$ is
suggested to be set around 0.5.

Figure 8 illustrates the RMSE values of $x_{1},x_{2}$ and $(x_{1}+x_{2})$ as
functions of $\eta$. As shown, when $\eta$ is greater than 0.5, although the
RMSE of the global signal $(x_{1}+x_{2})$ is low, the performance of the
decomposition is worse. Moreover, $\eta$ should not be too small, for then the
results will be similar to conventional MCA shown in Figure 6.

##  6 Engineering Examples

Figure 9: Outer race defect (left) and inner race (right) in the bearing. Table 4: Parameters of 552732QT bearing Inner Race (mm) | Outer Race(mm) | Roller (mm) | Number of rollers | Contact angle (degree)
---|---|---|---|---
160 | 290 | 34 | 17 | $0$

In this section, the proposed RTEA is applied to analyze vibration signals
collected from a rolling element bearing with compound defects on a
locomotive, using a SONY EX data acquisition system operating at a sampling
rate of $f_{s}=12.8$ kHz. The locomotive bearing with faults in inner and out
races is shown in Figure 9. The bearing (552732QT) parameters are given in
Table 4.

###  6.1 Compound faults detection

In this example, the vibration signals were collected at a constant shaft
speed of 360 r/min. Thus, based on the geometric parameters and rotational
frequency, the ball-pass frequency of the outer race is about $43.3$ Hz, and
that of the inner race is about $58.7$ Hz.

In practice, the regularization parameters in (12) can be estimated from the
noise level, or the deviation of the vibration signal without any fault under
an approximately identical experiment environment. Moreover, when the healthy
data is not available, the ‘noise’ level can still be determined by the
formula

$\displaystyle\hat{\sigma}=\mathsf{MAD}(y)/0.6745$ (28)

which is a conventional estimator of noise level used for wavelet-based
denoising [20], where $\mathsf{MAD}$ is the median absolute deviation defined
as

$\displaystyle\mathsf{MAD}(y):=\mathsf{median}(\left|y_{n}-\mathsf{median}(y)\right|).$
(29)

In this example, the formula (28) is used to estimate $\sigma$ directly from
the observation data, and then that value is used to choose the regularization
parameters from Table 3.

Figure 10 shows the results using RTEA, where $x_{1}$ is the transient
sequence generated by the fault in the outer race, and $x_{2}$ is the
transient sequence generated by the fault in the inner race. During the
recording, the outer race is stable and the inner race is rotating. Hence, the
amplitude of transients in $x_{2}$ (inner race fault) exhibit modulating
effect.

To further reveal the characteristic frequencies, the Hilbert envelope
spectrums of the extracted components are shown in Figure 11. The smoothed (by
lowpass filtering) profiles of the Hilbert envelope spectrum is also presented
to indicate the characteristic frequencies more clearly. Figure 11 shows that
the characteristic frequencies of outer and inner race are at about 45 and 60
Hz, and their harmonics are revealed by the peaks of the smoothed profiles.
Thus, the fault features of the two defects are clearly extracted by the
proposed approach.

As a comparison, we include the Welch’s estimate of the cyclic spectral
coherence for the analysis of cyclostationary signals [2, 3, 4]. The result is
illustrated in Figure 12. Note that, the signal length is 6400, the window
length is 32, the signal is divided into 14 overlapping blocks and the block
overlap is 21 with a Hanning window, which is suggested in Ref. [2]. When
observing the coherence, the cyclic frequency of the outer race (approximate
43.3 Hz) and its harmonic components can be identified. However, the
characteristic frequencies of inner race 58.7 Hz and its harmonic components
cannot be observed from Figure 12. This may be caused by the fact that the
component of inner race in this case is weak.

Figure 10: Example 2: Output of RTEA (proposed method) from data with compound
faults. Figure 11: Example 2: Hilbert envelope spectrums of the extracted
repetitive transient sequences. Figure 12:  Example 2: Cyclic spectral
coherence of the test signal given in Figure 10.

###  6.2 Example 3: Single fault detection

Figure 13: Example 3: Output of RTEA (proposed method) from data with outer
race fault only. Note that the Y-axis has a different scale to the previous
example.

In this example, the proposed method also works when the bearing has only one
fault. Figure 13 shows the measured data and the results from the bearing in
Figure 9, but there was only one defect on the outer race. The acceleration
signals were collected at a constant shaft speed of 481 r/min. Thus, the
characteristic fault frequencies of the outer and inner races are about $57.8$
Hz and $78.4$ Hz respectively.

Using the given information of the periods, the repetitive transient sequences
$x_{1}$, $x_{2}$ and the residual noisy signals are extracted simultaneously
from the measured vibration data, where $x_{1}$ is corresponds to outer race
defect. Repetitive transients can be observed from the extracted $x_{1}$,
where the fault frequency can be directly observed around 58 Hz, because there
are about 29 sparse transients evenly distributed with the 0.5 second. The
extracted characteristic frequency 58 Hz is approximately in accordance with
the outer race fault frequency of 57.8 Hz.

Moreover, component $x_{2}$ in Figure 13 is almost purely zero and it exhibits
almost no repetitive transients. This implies that the status of the inner
race is healthy.

##  7 Conclusion

This paper proposes a novel approach for the extraction of repetitive
transients with group-sparse structure in vibration signals, and for detecting
faults in rolling element bearings. To simultaneously extract both sparse
components and perform denoising, a repetitive group-sparsity based
optimization problem is formulated. To solve the problem, a computationally
efficient iterative algorithm, termed repetitive transients extraction
algorithm(RTEA) is derived. The non-convex penalty function is used to
strongly promote sparsity, and a condition is given so that the objective
function is strictly convex. Moreover, for practical problems, an approach to
select regularization and non-convexity parameters is provided.

## Appendices

##  Appendix A Proof of Proposition 1

###### Proof.

The first two terms of the objective function in (12) can be rewritten as

$\displaystyle F(u)=\frac{1}{2}\lVert
y-u\rVert_{2}^{2}+\lambda_{0}\sum_{n}\phi_{\epsilon}\Big{(}\Big{[}\sum_{k=0}^{K_{0}-1}[u]_{n+k}^{2}\Big{]}^{1/2};a_{0}\Big{)}$
(A.1)

where $u\in{}^{N}$ and $u-(x_{1}+x_{2})=\mathbf{0}$. Note that $F(u)$ is
exactly identical to GSSD problem {see Equation (20) in [13, Theorem 1]}.
Adopting Theorem 1 and Corollary 2 in [13], it can be shown that: when

$0\leq a_{0}<\frac{1}{K_{0}\lambda_{0}},$ (A.2)

$F(u)$ is strictly convex. Moreover, it is immediate that when $a=0$, function
(11) is convex. As a consequence, an equivalent problem of (12) is considered,
namely

$\displaystyle\\{u^{\mathsf{opt}},x_{1}^{\mathsf{opt}},x_{2}^{\mathsf{opt}}\\}$
$\displaystyle=\arg\min_{u,x_{1},x_{2}}F(u)+\sum_{i\in\\{1,2\\}}\lambda_{i}\Phi(x_{i};b_{i})$
$\displaystyle\qquad\text{ such that }u-x_{1}-x_{2}=\mathbf{0},$ (A.3)

which satisfies the convexity condition of equality constrained problem in
[11], and which implies that: when (12) satisfies (15), it is a convex
problem. ∎

##  Appendix B Derivation of majorizing function $R$ (13)

Since $\phi_{\epsilon}^{\mathsf{M}}$ in (7) majorizes $\phi_{\epsilon}$, an
upper bound of function $R$ can be found as

$\displaystyle\bar{R}(x_{1},x_{2},v;a_{0})$
$\displaystyle=\sum_{n}\phi_{\epsilon}^{\mathsf{M}}\Big{(}\Big{[}\sum_{k}[x_{1}+x_{2}]_{n+k}^{2}\Big{]}^{1/2},v_{n};a_{0}\Big{)}$
$\displaystyle=\sum_{n}\Bigg{\\{}\frac{1}{2\psi(v_{n};a_{0})}\sum_{k}[x_{1}+x_{2}]_{n+k}^{2}\Bigg{\\}}+C$
(B.1)

where $C$ is a constant does not depend on $x_{1}$ or $x_{2}$, and for any
$v\in{}^{N}$,

$\displaystyle R(x_{1},x_{2};a_{0})$
$\displaystyle\leq\bar{R}(x_{1},x_{2},v;a_{0}).$ (B.2)

After algebraic manipulations, $\bar{R}$ can be further expressed as

$\displaystyle\bar{R}(x_{1},x_{2},v(z);a_{0})=\frac{1}{2}\sum_{n}r_{0}(n,z)[x_{1}+x_{2}]_{n}^{2}+C,$
(B.3)

where vector $v$ is dependent on another vector $z\in{}^{N}$,as

$\displaystyle[v(z)]_{n}=\Big{[}\sum_{k=0}^{K_{0}-1}[z]_{n+k}^{2}\Big{]}^{1/2},$
(B.4)

and $r_{0}:\mathbb{Z}\times{}^{N}\to\real,$ is defined by

$\displaystyle
r_{0}(n,z):=\sum_{j=0}^{K_{0}-1}\frac{1}{\psi\Big{(}[v(z)]_{n-j};a_{0}\Big{)}}=\sum_{j=0}^{K_{0}-1}\frac{1}{\psi\bigg{(}\Big{[}\sum_{k}[z]_{n-j+k}^{2}\Big{]}^{1/2};a_{0}\bigg{)}}$
(B.5)

which is similar to Equation (37) of [13].

Note that, in this case the inequality of (B.2) is still valid and the
equality holds when

$\displaystyle z=x_{1}+x_{2}.$ (B.6)

Moreover, considering a simple inequality that: for any
$\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}\in\real$,

$\displaystyle(\alpha_{1}+\alpha_{2})^{2}\leq(\alpha_{1}+\alpha_{2})^{2}+\big{[}(\alpha_{1}-\beta_{1})-(\alpha_{2}-\beta_{2})\big{]}^{2},$
(B.7)

whose right side can be re-written as

$\displaystyle
2\alpha_{1}^{2}+2\alpha_{2}^{2}-2(\beta_{1}-\beta_{2})\alpha_{1}-2(\beta_{2}-\beta_{1})\alpha_{2}+(\beta_{1}-\beta_{2})^{2},$
(B.8)

where $\alpha_{1}$ and $\alpha_{2}$ are de-coupled and the equality holds when
$\alpha_{1}=\beta_{1}$ and $\alpha_{2}=\beta_{2}$. Then, using the inequality
(B.7) element-wise, an upper bound of $\bar{R}$ can be found as

$\displaystyle R^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2};a_{0})$ $\displaystyle=$
$\displaystyle\sum_{n}\Big{\\{}r_{0}(n,z_{1}+z_{2})([x_{1}]_{n}^{2}+[x_{2}]_{n}^{2}-[z_{1}-z_{2}]_{n}[x_{1}]_{n}-[z_{2}-z_{1}]_{n}[x_{2}]_{n})\Big{\\}}+C(z_{1},z_{2}),$
(B.9)

where $C(z_{1},z_{2})$ is a constant only dependent on $z_{1}$ and $z_{2}$.
Note that since (B) is an upper bound of $\bar{R}$, consequently

$\displaystyle R^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2};a_{0})\geq
R(x_{1},x_{2};a_{0}).$ (B.10)

Furthermore, when $z_{1}=x_{1}$ and $z_{2}=x_{2}$, using (B.7) it can be seen
that
$R^{\mathsf{M}}(x_{1},x_{2},z_{1},z_{2};a_{0})=\bar{R}(x_{1},x_{2},v(z_{1}+z_{2});a_{0})$,
and this implies that

$\displaystyle
R^{\mathsf{M}}(x_{1},x_{2},x_{1},x_{2};a_{0})=R(x_{1},x_{2};a_{0}).$ (B.11)

Therefore, equation (B) is a majorizer of function $R$ in (13).

## References

  * [1] S. Abbasion, A. Rafsanjani, A. Farshidianfar, and N. Irani.  Rolling element bearings multi-fault classification based on the wavelet denoising and support vector machine.  Mech. Syst. Signal Process., 21(7):2933–2945, 2007.
  * [2] J. Antoni.  Cyclic spectral analysis in practice.  Mech. Syst. Signal Process., 21(2):597–630, 2007.
  * [3] J. Antoni.  Cyclic spectral analysis of rolling-element bearing signals: Facts and fictions.  304(3–5):497–529, 2007.
  * [4] J. Antoni.  Cyclostationarity by examples.  Mech. Syst. Signal Process., 23(4):987–1036, 2009.
  * [5] J. Antoni and R. B. Randall.  The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines.  Mech. Syst. Signal Process., 20(2):308–331, 2006.
  * [6] I. Bayram, P.-Y. Chen, and I. W. Selesnick.  Fused lasso with a non-convex sparsity inducing penalty.  In Proc. ICASSP 2014, pages 4156–4160, May 2014.
  * [7] G. Bin, J. Gao, X. Li, and B. S. Dhillon.  Early fault diagnosis of rotating machinery based on wavelet packets -empirical mode decomposition feature extraction and neural network.  Mech. Syst. Signal Process., 27:696–711, 2012.
  * [8] T. Boukra, A. Lebaroud, and G. Clerc.  Statistical and neural-network approaches for the classification of induction machine faults using the ambiguity plane representation.  IEEE Trans. Ind. Electron., 60(9):4034–4042, 2013.
  * [9] R. Boustany and J. Antoni.  A subspace method for the blind extraction of a cyclostationary source: Application to rolling element bearing diagnostics.  Mech. Syst. Signal Process., 19(6):1245–1259, 2005.
  * [10] R. Boustany and J. Antoni.  Blind extraction of a cyclostationary signal using reduced-rank cyclic regression—a unifying approach.  Mech. Syst. Signal Process., 22(3):520–541, 2008.
  * [11] S. Boyd and L. Vandenberghe.  Convex Optimization.  Cambridge University Press, 2004.
  * [12] B. Chen, Z. Zhang, C. Sun, B. Li, Y. Zi, and Z. He.  Fault feature extraction of gearbox by using overcomplete rational dilation discrete wavelet transform on signals measured from vibration sensors.  Mech. Syst. Signal Process., 33:275–298, 2012.
  * [13] P.-Y. Chen and I. W. Selesnick.  Group-sparse signal denoising: Non-convex regularization, convex optimization.  IEEE Trans. Signal Process., 62(13):3464–3478, July 2014.
  * [14] S. Chen and D. L. Donoho.  Basis pursuit.  In 1994 Conference Record of the Twenty-Eighth Asilomar, Conference on Signals, Systems and Computers, 1994., volume 1, pages 41–44, October 1994.
  * [15] R. R. Coifman and D. L. Donoho.  Translation-invariant de-noising.  In Wavelet and statistics, pages 125–150. Springer-Verlag, 1995\.
  * [16] L. Cui, J. Wang, and S. Lee.  Matching pursuit of an adaptive impulse dictionary for bearing fault diagnosis.  J. Sound Vib., 333(10):2840–2862, 2014.
  * [17] L. Cui, N. Wu, C. Ma, and H. Wang.  Quantitative fault analysis of roller bearings based on a novel matching pursuit method with a new step-impulse dictionary.  Mech. Syst. Signal Process., 68:34–43, 2016.
  * [18] L. Cui, Y. Zhang, F. Zhang, J. Zhang, and S. Lee.  Vibration response mechanism of faulty outer race rolling element bearings for quantitative analysis.  J. Sound Vib., 364:67–76, 2016.
  * [19] Y. Ding and I. W. Selesnick.  Sparsity-based correction of exponential artifacts.  Signal Process., 120:236–248, March 2016.
  * [20] D. Donoho, I. Johnstone, and I. M. Johnstone.  Ideal spatial adaptation by wavelet shrinkage.  Biometrika, 81:425–455, 1993.
  * [21] Z. Feng, M. Liang, and F. Chu.  Recent advances in time–frequency analysis methods for machinery fault diagnosis: A review with application examples.  Mech. Syst. Signal Process., 38(1):165–205, 2013.
  * [22] M. Figueiredo, J. Bioucas-Dias, and R. Nowak.  Majorization-minimization algorithms for wavelet-based image restoration.  IEEE Trans. Image Process., 16(12):2980–2991, December 2007.
  * [23] Q. He and X. Wang.  Time–frequency manifold correlation matching for periodic fault identification in rotating machines.  J. Sound Vib., 332(10):2611–2626, 2013.
  * [24] W. He, Y. Ding, Y. Zi, and I. W. Selesnick.  Sparsity-based algorithm for detecting faults in rotating machines.  Mech. Syst. Signal Process., 72–-73:46–64, May 2016.
  * [25] W. He, Y. Zi, B. Chen, S. Wang, and Z. He.  Tunable Q-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis.  Sci. China Technol. Sci., 56(8):1956–1965, 2013.
  * [26] W. He, Y. Zi, B. Chen, F. Wu, and Z. He.  Automatic fault feature extraction of mechanical anomaly on induction motor bearing using ensemble super-wavelet transform.  Mech. Syst. Signal Process., 54-55:457–480, 2015.
  * [27] H. Jiang, J. Chen, G. Dong, T. Liu, and G. Chen.  Study on Hankel matrix-based SVD and its application in rolling element bearing fault diagnosis.  Mech. Syst. Signal Process., 52:338–359, 2015.
  * [28] K. Lange, D. Hunter, and I. Yang.  Optimization transfer using surrogate objective functions.  J. of Comp. Graph. Statist., 9:1–20, 2000.
  * [29] Y. Lei, Z. He, Y. Zi, and X. Chen.  New clustering algorithm-based fault diagnosis using compensation distance evaluation technique.  Mech. Syst. Signal Process., 22(2):419–435, 2008.
  * [30] Y. Lei, J. Lin, Z. He, and M. J Zuo.  A review on empirical mode decomposition in fault diagnosis of rotating machinery.  Mech. Syst. Signal Process., 35(1):108–126, 2013.
  * [31] Z. Li, X. Yan, C. Yuan, Z. Peng, and L. Li.  Virtual prototype and experimental research on gear multi-fault diagnosis using wavelet-autoregressive model and principal component analysis method.  Mech. Syst. Signal Process., 25(7):2589–2607, 2011.
  * [32] M. Liang and I. Soltani Bozchalooi.  An energy operator approach to joint application of amplitude and frequency-demodulations for bearing fault detection.  Mech. Syst. Signal Process., 24(5):1473 – 1494, 2010.
  * [33] A. Malhi and R. X. Gao.  PCA-based feature selection scheme for machine defect classification.  IEEE Trans. Instrum. Meas., 53(6):1517–1525, 2004.
  * [34] Y. Qin, Y. Mao, and B. Tang.  Vibration signal component separation by iteratively using basis pursuit and its application in mechanical fault detection.  J. Sound Vib., 332(20):5217–5235, 2013.
  * [35] R. B. Randall and J. Antoni.  Rolling element bearing diagnostics – a tutorial.  Mech. Syst. Signal Process., 25(2):485–520, 2011.
  * [36] I. W. Selesnick and I. Bayram.  Sparse signal estimation by maximally sparse convex optimization.  IEEE Trans. Signal Process., 62(5):1078–1092, 2014.
  * [37] I. W. Selesnick, A. Parekh, and I. Bayram.  Convex 1-d total variation denoising with non-convex regularization.  IEEE Signal Processing Letters, 22(2):141–144, 2015.
  * [38] W. A. Smith and R. B. Randall.  Rolling element bearing diagnostics using the case western reserve university data: A benchmark study.  Mech. Syst. Signal Process., 64:100–131, 2015.
  * [39] J.-L. Starck, M. Elad, and D. Donoho.  Redundant multiscale transforms and their application for morphological component analysis.  Advances in Imaging and Electron Physics, 132:287–348, 2004.
  * [40] W. Su, F. Wang, H. Zhu, Z. Zhang, and Z. Guo.  Rolling element bearing faults diagnosis based on optimal morlet wavelet filter and autocorrelation enhancement.  Mech. Syst. Signal Process., 24(5):1458–1472, 2010.
  * [41] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight.  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society Series B, pages 91–108, 2005.
  * [42] Z. Wang, J. Chen, G. Dong, and Y. Zhou.  Constrained independent component analysis and its application to machine fault diagnosis.  Mech. Syst. Signal Process., 25(7):2501–2512, 2011.
  * [43] R. Yan and R. X. Gao.  Harmonic wavelet-based data filtering for enhanced machine defect identification.  J. Sound Vib., 329(15):3203–3217, 2010.
  * [44] R. Yan, R. X. Gao, and X. Chen.  Wavelets for fault diagnosis of rotary machines: A review with applications.  Signal Process., 96:1–15, 2014.
  * [45] H. Yang, J. Mathew, and L. Ma.  Fault diagnosis of rolling element bearings using basis pursuit.  Mech. Syst. Signal Process., 19(2):341–356, 2005.
  * [46] X. Zhang and J. Zhou.  Multi-fault diagnosis for rolling element bearings based on ensemble empirical mode decomposition and optimized support vector machines.  Mech. Syst. Signal Process., 41(1):127–140, 2013.