# Level compressibility of certain random unitary matrices

Eugene Bogomolny  Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France

###### Abstract

The value of spectral form factor at the origin, called level compressibility,
is an important characteristic of random spectra. The paper is devoted to
analytical calculations of this quantity for different random unitary matrices
describing models with intermediate spectral statistics. The computations are
based on the approach developed by G. Tanner in [J. Phys. A: Math. Gen. 34,
8485 (2001)] for chaotic systems. The main ingredient of the method is the
determination of eigenvalues of a transition matrix whose matrix elements
equal squared moduli of matrix elements of the initial unitary matrix. The
principal result of the paper is the proof that the level compressibility of
random unitary matrices derived from the exact quantisation of barrier
billiards and consequently of barrier billiards themselves is equal to $1/2$
irrespectively of the height and the position of the barrier.

##  I Introduction

The leading idea behind statistical descriptions of complex deterministic
quantum problems is that quantum characteristics (e.g., eigenenergies) of a
large variety of such problems are so erratic and irregular that their precise
values are irrelevant (like the position of a molecule in the air) and only
their statistical properties are of importance. As matrices are inherent in
quantum mechanics, random matrices occupy a predominant place in the
application of statistics to quantum problems wigner . In a typical setting
one tries to find a random matrix ensemble whose eigenvalues have the same
statistical distributions as (high-excited) eigenenergies of a given
deterministic quantum problem. Till now this query has been figured out only
for two limiting classes of quantum problems: (i) models whose classical limit
is integrable berry_tabor  and (ii) models whose classical limit is chaotic
bgs . For generic integrable models quantum eigenenergies are distributed as
eigenvalues of diagonal matrices with independent identically distributed
(i.i.d.) elements which means that their correlation functions after unfolding
coincide with the ones of the Poisson distribution berry_tabor . For generic
chaotic systems it was conjectured in bgs  that their eigenenergies are
distributed as eigenvalues of standard random matrix ensembles (GOE, GUE, GSE)
depending only on system symmetries whose correlation functions are known
explicitly mehta . The difference between these two cases is clearly seen from
the limiting behaviour of their nearest-neighbour distribution $P_{0}(s)$
which is the probability density that two eigenvalues are separated by a
distance $s$ and there are no other eigenvalues in-between. For the Poisson
statistics there is no level repulsion which means that $\lim_{s\to
0}P_{0}(s)\neq 0$ and for large argument $P_{0}(s)$ decreases exponentially
with $s$. For standard random matrix ensembles levels repel each other,
$\lim_{s\to 0}P_{0}(s)=0$, and $P_{0}(s)\sim\exp(-as^{2})$ when $s\to\infty$.

These two big conjectures form a cornerstone of quantum chaos and have been
successfully applied to various problems from nuclear physics to number
theory. Nevertheless, they do not cover all possible types of models.
Especially intriguing is the class of pseudo-integrable billiards (see, e.g.,
richens_berry ) which are 2-dimensional polygonal billiards whose angles
$\theta_{j}$ are rational multiplies of $\pi$

$\theta_{j}=\frac{m_{j}}{n_{j}}\pi$

with co-prime integers $m_{j}$ and $n_{j}$. A peculiarity of such billiards is
seen in the fact that their classical trajectories belong to a 2-dimensional
surface of genus $g$ related with angles as follows katok

$g=1+\frac{N_{n}}{2}\sum_{j}\frac{m_{j}-1}{n_{j}}$

where $N_{n}$ is the least common multiple of all denominators $n_{j}$.
Consequently, any such model with at least one numerator $m_{j}>1$ is neither
integrable (which would imply that trajectories belong to a 2-dimensional
torus with $g=1$) nor fully chaotic (in which case trajectories should cover
3-dimensional surface of constant energy) and the aforementioned conjectures
cannot be applied to such systems. Numerical calculations show that for many
pseudo-integrable billiards spectral statistical properties of corresponding
quantum problems differ from both the Poisson statistics and the random matrix
statistics mentioned above (see, e.g., billiard_1 ; billiard_2  and references
therein). In particular, for these models (i) $\lim_{s\to 0}P_{0}(s)=0$ as for
standard random matrix ensembles but (ii) $P_{0}(s)\sim\exp(-bs)$ for large
$s$ as for the Poisson statistics. Such hybrid statistics, labeled
intermediate statistics, had been first observed in the Anderson model at the
metal-insulate transition altshuler ; shklovskii  and they constitute a
special, interesting but poorly investigated class of spectral statistics.

Probably the simplest example of pseudointegrable systems is the so-called
barrier billiard which is a rectangular billiard with a barrier inside
sketched in figure 1(a). The quantum problem for this model consists in
solving the Helmholtz equation

$(\Delta+E_{\alpha})\Psi_{\alpha}(x,y)=0$

imposing that eigenfunction $\Psi_{\alpha}(x,y)$ obeys (e.g.) the Dirichlet
conditions on the boundary of the rectangle as well as on the barrier

$\Psi_{\alpha}(x,y)|_{\mathrm{boundary}}=0,\qquad\Psi_{\alpha}(x,y)|_{\mathrm{barrier}}=0.$

(a)

(b)

Figure 1: (a) Barrier billiard. (b) An infinite slab with a half-plane inside.
Numbers indicate 3 possible channels.

Calculating the exact $S$-matrix for the scattering inside the infinite slab
with a barrier depicted in figure 1(b), it has been demonstrated in billiard_1
; billiard_2  that spectral statistics of this model is the same as the
statistics of eigenvalues of the following $N\times N$ random unitary matrix

$B_{m,n}=e^{i\Phi_{m}}\frac{L_{m}L_{n}}{x_{m}+x_{n}},\qquad m,n=1,\ldots,N$
(1)

where $\Phi_{m}$ are i.i.d. random variables uniformly distributed on interval
$[0,2\pi)$ and

$L_{m}=\sqrt{2x_{m}\prod_{k\neq m}\frac{x_{m}+x_{k}}{x_{m}-x_{k}}}$ (2)

where coordinates $x_{m}$ depend on the position of the barrier.

Define the following quantities (momenta) of propagating modes in each of 3
channels indicated in figure 1(b)

$p_{m}^{(1)}=\sqrt{k^{2}-\frac{\pi^{2}m^{2}}{b^{2}}},\qquad
p_{m}^{(2)}=\sqrt{k^{2}-\frac{\pi^{2}m^{2}}{(b-h)^{2}}},\qquad
p_{m}^{(3)}=\sqrt{k^{2}-\frac{\pi^{2}m^{2}}{h^{2}}}$

If the ratio $h/b$ is an irrational number, coordinates $x_{m}$ have the
following form

$\vec{x}=b\Big{(}\underbrace{p_{1}^{(1)},\ldots,p_{N_{1}}^{(1)}}_{\text{$N_{1}$}},\underbrace{-p_{1}^{(2)},\ldots,-p_{N_{2}}^{(2)}}_{\text{$N_{2}$}},\underbrace{-p_{1}^{(3)},\ldots,-p_{N_{3}}^{(3)}}_{\text{$N_{3}$}}\Big{)}.$
(3)

Here $N_{j}$ with $j=1,2,3$ are the numbers of propagating modes in each
channel

$N_{1}=\left[\frac{kb}{\pi}\right],\qquad
N_{2}=\left[\frac{k(b-h)}{\pi}\right],\qquad
N_{3}=\left[\frac{kh}{\pi}\right]$ (4)

where $[x]$ is the largest integer $\leq x$ and the total dimension of the
$B$-matrix is $N=N_{1}+N_{2}+N_{3}$.

When the ratio $h/b$ is a rational number, $h/b=p/q$ with co-prime integers
$p$ and $q$ ($p<q$), there are exist exact plane wave solutions of barrier
billiard equal zero at the whole line passing through the barrier. It is
natural to disregard them and take into account only non-trivial eigenvalues.
In such case coordinates $x_{m}$ have to be chosen as indicated below

$\vec{x}=b\Big{(}\underbrace{p_{1}^{(1)},\ldots,p_{k}^{(1)},\ldots,p_{N_{1}}^{(1)}}_{\text{$k\neq
0\mod
q$}},\underbrace{-p_{1}^{(2)},\ldots,-p_{N_{2}}^{(2)}}_{\text{$N_{2}$}},\underbrace{-p_{1}^{(3)},\ldots,-p_{k}^{(3)},\ldots,-p_{N_{3}}^{(3)}}_{\text{$k\neq
0\mod p$}}\Big{)}\,.$ (5)

The dimension of this vector is

$N^{(\mathrm{r})}=N_{1}+N_{2}+N_{3}-2N_{0},\qquad N_{0}=\left[\frac{kb}{\pi
q}\right].$ (6)

The matrix $B$ can be generalised for arbitrary vector $x_{m}$ provided the
following interlacing conditions are fulfilled

$|x_{1}|>|x_{2}|>\ldots>|x_{N}|,\qquad x_{m}=(-1)^{m+1}|x_{m}|.$

Exact correlation functions for the $B$-matrix are unknown at present. In
billiard_1 ; billiard_2  it was argued that an approximate Wigner-type surmise
for this matrix corresponds to the so-called semi-Poisson distribution gerland
. In particular, it implies that the probability density $P_{n}(s)$ that two
levels are separated by a distance $s$ and there are exactly $n$ levels in-
between (after the standard unfolding) is given by the following expression

$P_{n}(s)=\frac{2^{2n+2}}{(2n+1)!}s^{2n+1}e^{-2s},\qquad n\geq 0\,.$

Numerical calculations presented in billiard_1 ; billiard_2  agree with these
simple formulas.

This paper is devoted to the calculation of another important characteristic
of spectral statistics, namely the level (or spectral) compressibility. This
quantity is determined by the limiting behaviour of the variance of the number
of levels inside a given interval. More precisely, let $N(L)$ be the number of
eigenvalues in an interval $L$ unfolded to the unit mean density which means
that the mean number of levels in interval $L$ equals $L$, $\langle
N(L)\rangle=L$. By definition the number variance is
$\Sigma^{(2)}(L)\equiv\langle(N(L)-L)^{2}\rangle$. If for large $L$

$\Sigma^{(2)}(L)\underset{L\to\infty}{\longrightarrow}\chi L$ (7)

constant $\chi$ is called the level compressibility. The importance of this
quantity follows from the fact that for integrable systems with the Poisson
statistics $\chi=1$ but for standard random matrix statistics typical for
chaotic models $\chi=0$. For all examples of intermediate statistics it was
observed that $0<\chi<1$.

The conventional way of determination the level compressibility for dynamical
systems is the summation over all periodic orbits in the diagonal
approximation initiated in rigidity . For the symmetric barrier billiard with
$h/b=1/2$ and $d/a=1/2$ it has been demonstrated in wiersig  that $\chi=1/2$.
The same value had been obtained in olivier_thesis  for the case $h/b=1/2$ and
arbitrary barrier height $d/a$. Finally for $h/b=p/q$ with co-prime integers
$p$, $q$ and irrational values of $d/a$ it has been proved in giraud  that

$\chi=\frac{1}{2}+\frac{1}{q}$

but in the calculations exact eigenvalues whose eigenfunctions are zero on the
whole line $y=h$ (see figure 1)(a) have not been excluded. When these trivial
eigenvalues are removed the answer is $\chi=1/2$ private_communication .
Therefore direct (and quite tedious) calculations suggest that for (almost)
all positions and heights of the barrier the level compressibility is the same
as for the semi-Poisson statistics gerland

$\chi=\frac{1}{2}$ (8)

but the reason of this universality remains obscure.

The purpose of the paper is to find analytically the spectral compressibility
for barrier billiards and for a few other models directly from the
corresponding random unitary matrices. To achieve the goal it is convenient to
slightly generalise the method developed in tanner  for random unitary
matrices appeared in the quantisation of quantum graphs gnutzmann . The method
is briefly explained in Section II. In Section III this method is applied to
random unitary matrices derived in marklof  by quantisation of a simple
interval-exchange map. In this case the transition matrix is a circulant
matrix whose eigenvalues are known explicitly. The results coincide with exact
level compressibility for these models obtained in map -integrable_ensembles .
This example gives credit to the method and permits to explain its main
features without unnecessary complications. In particular it clarifies the
situation (not covered by Refs. map -integrable_ensembles ) when a parameter
entered the matrix takes an irrational value. Numerically it has been observed
general_case  that in such case the spectral statistics are well described by
the ones of chaotic systems (GOE or GUE) though the Lyapunov exponent of the
underlying classical map is always zero.

The main part of the paper is devoted to the derivation of the level
compressibility for random matrices associated with barrier billiards (1). The
calculation are more complicated as no eigenvalues (except one) are known
analytically. The simplest case of the symmetric barrier billiard with
$h/b=1/2$ is investigated in Section IV. To get tractable expressions a kind
of paraxial approximation is developed which permits to control the largest
terms. By using such approximation the transition matrix is transformed into a
Toeplitz matrix with a quickly decreasing symbol which allows to find its
eigenvalues for large matrix dimension. The result of this Section is that the
level compressibility equals $1/2$ in accordance with periodic orbit
calculations in wiersig  and olivier_thesis .

In Section V random unitary matrices corresponding to barrier billiards with
irrational ratio $h/b$ (with coordinates given by (3)) are considered. In this
case the transition matrix in paraxial approximation contains quickly
oscillating terms and, consequently, has forbidden zones in the spectrum. In
spite of that one can argue that largest moduli eigenvalues are insensitive to
fast oscillations and are determined solely by a matrix averaged over such
oscillations. The final matrix is also a Toeplitz type which allows of
analytical calculations proving that the level compressibility is again equal
$1/2$.

Section VI is addressed to the calculation of level compressibility in the
most complicated case of barrier billiards with rational ratio $h/b=p/q\neq
1/2$, the only one for which the direct calculation of the level
compressibility by the summation over periodic orbits was not yet done. The
computation is cumbersome but in the end one comes to the conclusion that the
level compressibility remains equal $1/2$. In other words the level
compressibility of barrier billiards is universal (i.e., independent of the
barrier position and its height) and coincides with the semi-Poisson
prediction gerland .

Section VII gives a brief summary of the results. A few technical points are
discussed in Appendices A-C.

##  II Generalities

It is well known that the level compressibility (7) is related with the two-
point correlation form factor $K(\tau)$ as follows

$\chi=\lim_{\tau\to 0}K(\tau)\,.$ (9)

For $N\times N$ random unitary matrices $U$ the form factor can, conveniently,
be written in the following concise form (see, e.g.,tanner )

$K(\tau)=\Big{\langle}\frac{1}{N}\left|\mathrm{Tr}\,U^{n}\right|^{2}\Big{\rangle},\qquad\tau=\frac{n}{N}$
(10)

where the average is taken either on different realisations of random
parameters, or over of a small window of $\tau$, or the both.

Unitary matrices considered in the paper all have the product form

$U_{j,k}=e^{i\Phi_{j}}w_{j,k}$ (11)

where $\Phi_{j}$ are i.i.d. random phases uniformly distributed between $0$
and $2\pi$ and matrix $w_{j,k}$ is a fixed unitary matrix.

In tanner  (see gnutzmann  for more detailed discussion) it was shown that for
such unitary matrices the averaging over random phases leads in the diagonal
approximation to the following formula

$K^{(\mathrm{diag})}\Big{(}\frac{n}{N}\Big{)}=g\frac{n}{N}\mathrm{Tr}\,\left(T^{n}\right)$
(12)

where matrix elements of matrix $T$, called below the transition matrix, are
squared moduli of matrix elements of matrix $U$

$T_{j,k}=|U_{j,k}|^{2}=|w_{j,k}|^{2}\,.$ (13)

For systems without time-reversal invariance $g=1$ and for models with time-
reversal invariance $g=2$ (here only cases of $g=1$ are considered). Due to
the unitarily of matrix $U$ the $T$-matrix is double stochastic matrix,
$\sum_{j}T_{j,k}=1$ and $\sum_{k}T_{j,k}=1$, thus having the meaning of
classical transition matrix.

Let $\Lambda_{\beta}$ be eigenvalues of the $T$-matrix

$T_{j,k}u_{k}(\beta)=\Lambda_{\beta}u_{j}(\beta),\qquad\beta=0,\ldots,N-1\,.$

From (12) it follows that

$K^{(\mathrm{diag})}\Big{(}\frac{n}{N}\Big{)}=\frac{n}{N}\sum_{\beta=0}^{N-1}\Lambda_{\beta}^{n}\,.$
(14)

The unitarity imposes that one eigenvalues $\Lambda_{0}=1$. Perron-Frobenius
theorem states that all other eigenvalues $|\Lambda_{\beta}|\leq 1$. Let the
set of eigenvalues be ordered
$1=\Lambda_{0}\geq\Lambda_{1}\geq\ldots\geq\Lambda_{N-1}$.

One has

$K^{(\mathrm{diag})}(\tau)=\tau\left(\mathrm{Tr}\,\left(T^{n}\right)\right)=\tau+\tau\sum_{\beta=1}^{N-1}\Lambda_{\beta}^{N\tau},\qquad\sum_{\beta=1}^{N-1}\Lambda_{\beta}^{N\tau}\leq(N-1)|\Lambda_{1}|^{N\tau}\,.$
(15)

It has been noted in tanner  that if

$\lim_{N\to\infty}|\Lambda_{1}|^{N}=0$ (16)

then the second term in (15) goes to zero for all finite $\tau$ and
$K^{(\mathrm{diag})}(\tau)=\tau$ for small $\tau$. Consequently it is
reasonable to conjecture (as it has been done in tanner ) that the whole
spectral statistics of such matrices will be well described by standard random
matrix formulas.

For matrices discussed below criterium (16) is not fulfilled. Instead, in all
considered cases largest moduli eigenvalues of transition matrices have the
form $1-\mathcal{O}(1/N)$. To calculate the level compressibility from (14)
the summation over all such eigenvalues is performed analytically and then the
limit $\tau\to 0$ is taken.

It is well known that the form factor is not a self-average quantity. It has
strong fluctuations and necessarily requires a smoothing. There exist two
different sources of fluctuations. The first is related with random phases in
matrices (11). Eq. (12) corresponds to the averaging over these random phases
in the diagonal approximation. The second has its roots in non-smoothness of
$K^{(\mathrm{diag})}(n/N)$ for different $n$ and could be removed by a
smoothing over a small interval of $\tau$. (A trivial example is $(-1)^{n}$.)

##  III Interval-exchange matrices

This Section is devoted to the calculation of the level compressibility for
special unitary matrices derived in marklof  by quantisation of a simple
2-dimensional parabolic map. Slightly generalising their result map  one can
write these matrices in the following form

$M_{n,m}=e^{i\Phi_{n}}\frac{1-e^{2\pi i\alpha N}}{N\big{(}1-e^{2\pi
i(n-m+\alpha N)/N}\big{)}},\qquad n,m=1,\ldots,N.$ (17)

Here $\alpha$ is a real parameter and $\Phi_{n}$ are i.i.d. random variables
uniformly distributed between $0$ and $2\pi$ (the case with ’time-reversal
symmetry’ when $\Phi_{N-n+1}=\Phi_{n}$ requires only multiplication the
formulas below by factor $g=2$ as indicated in (12)).

When $\alpha$ is a rational number $\alpha=p/q$ with co-prime integers $p$ and
$q$ the original classical map is an interval-exchange map and, as it was
shown in detail in map -integrable_ensembles , spectral statistics of matrices
$\eqref{map}$ in the limit of $N\to\infty$ are unusual and peculiar. It
appears that the limiting results depend on the residue of $pN\mod q$ (when
$pN\equiv 0\mod q$ matrix (17) have explicit eigenvalues not interesting for
our purposes). For example, if $\alpha=1/5$ there are 2 possibilities:
$N\equiv\pm 1$ mod $5$ and $N\equiv\pm 2\mod 5$. In the first case the
nearest-neighbour distribution is

$P_{0}(s)=\frac{5^{5}}{4!}s^{4}e^{-5s}$

but for the second one the exact result is different

$P_{0}(s)=(a_{2}s^{2}+a_{3}s^{3}+a_{4}s^{4}+a_{5}s^{5}+a_{6}s^{6})e^{-5s}$

where $a_{2}=625/2-275\sqrt{5}/2$, $a_{3}=3125/2-1375\sqrt{5}/2$,
$a_{4}=71875/48+33125\sqrt{5}/48$, $a_{5}=-15625/3+9375\sqrt{5}/4$,
$a_{6}=1015625/288-453125\sqrt{5}/288$.

Though the spectral correlation functions for $\alpha=p/q$ are different for
different residues $pN\mod q$, the calculations show that the spectral
compressibility for all residues remains the same map ; general_case

$\chi=\frac{1}{q}\,.$ (18)

Matrices (17) with irrational $\alpha$ were investigated numerically in
general_case  and it was observed that their spectral statistics are well
described by standard random matrix ensembles (GOE or GUE). In particular, it
implies that in such case

$\chi=0\,.$ (19)

Below it is shown that values (18) and (19) can easily be recovered by the
discussed method.

The transition matrix (13) for the discussed case has the form

$T_{n,m}\equiv|M_{n,m}|^{2}=\frac{\sin^{2}(\pi\alpha
N)}{N^{2}\sin^{2}\big{(}\pi(n-m+\alpha N)/N\big{)}}$ (20)

This is a circulant matrix and its eigenvalues are simply the Fourier
transforms of its matrix elements

$\Lambda_{\beta}=\frac{\sin^{2}(\pi\alpha
N)}{N^{2}}\sum_{s=0}^{N-1}\frac{e^{-2\pi i\beta
s/N}}{\sin^{2}\big{(}\pi(s+\alpha N)/N\big{)}}$ (21)

Differentiating the both sides of the identity with integer $s$ on $z$

$\sum_{m=0}^{N-1}e^{2\pi im(s+z)/N}=\frac{1-e^{2\pi iz}}{1-e^{2\pi
i(s+z)/N}},$

after straightforward transformations one proves that

$\frac{\sin^{2}(\pi z)}{\sin^{2}\big{(}\pi(s+z)/N\big{)}}=-2i\sin(\pi
z)e^{-\pi iz}\sum_{m=0}^{N-1}\left(m-N\frac{e^{\pi iz}}{2i\sin(\pi
z)}\right)e^{2\pi im(s+z)/N}$

from which is follows that eigenvalues (21) are

$\Lambda_{\beta}=\left(1-\frac{\beta}{N}\Big{(}1-e^{-2\pi i\alpha
N}\Big{)}\right)e^{2\pi i\beta\alpha},\qquad\beta=0,1\,\ldots,N-1\,.$ (22)

Notice that $\Lambda_{0}=1$, $\Lambda_{N-\beta}=\Lambda_{\beta}^{*}$ for
$\beta=1,\ldots,N-1$.

Consider first the case of rational $\alpha=p/q$. Assume that $N\neq 0$ mod
$q$ and calculate the form factor from (14) separately for $n\equiv r$ mod $q$
with $r=0,1,\ldots,q-1$ (i.e., $n=qt+r$ with integer $t$)

$K\left(\frac{qt+r}{N}\right)=\frac{qt+r}{N}\left[1+2\mathrm{Re}\,\sum_{\beta=1}^{N/2}\left(1-\frac{\beta}{N}\Big{(}1-e^{-2\pi
i\alpha N}\Big{)}\right)^{qt+r}e^{2\pi i\beta pr/q}\right].$

Here we do not take into account that when $N=$ even the term with $\beta=N/2$
is real and has additional factor $1/2$.

Notice that the phase depends only on $r$. For large $N$ one can put the first
factor in the exponent and sum the geometric progression from 1 to infinity.
The answer is

$K_{r}(\tau)=\tau+2\tau\mathrm{Re}\left(\frac{e^{-\xi\tau+2\pi
ipr/q}}{1-e^{-\xi\tau+2\pi ipr/q}}\right),\qquad\tau=\frac{qt+r}{N}$ (23)

where $\xi=1-e^{-2\pi ipN/q}=1-e^{-2\pi ik/q}$ with $k=pN\mod q$.

When $\tau\to 0$ all terms except the one with $r=0$ tends to zero as they do
not have a pole at $\tau=0$. The remaining term equals

$\lim_{\tau\to
0}K_{0}(\tau)=2\,\mathrm{Re}\left(\frac{1}{\xi}\right)=\frac{2\,\mathrm{Re}\,(1-e^{2\pi
iz})}{2-2\cos(2\pi z)}=1\,.$

It means that after the averaging over random phases the limiting value of the
form factor $K_{r}(\tau)$ strongly depends on small changes of
$\tau=(qt+r)/N$. For $r=0$ $K_{0}(\tau)$ tends to $1$ for small $\tau\ll 1$
but all other terms with $r=1,\ldots,q-1$ $K_{r}(\tau)$ tend to $0$. The
difference between these different values of $\tau$ is very small, of the
order of $1/N$. Therefore after the averaging over any small (but finite)
interval of $\tau$ one gets

$\lim_{\tau\to 0}\langle K(\tau)\rangle=\frac{1}{q}$

which agrees with (18) obtained in general_case  by a different method.

(a)

(b)

Figure 2: Form factor for the matrix (17) with $\alpha=1/5$ and (a)
$N=399\equiv-1\,\mathrm{mod}\,5$ and (b) $N=398\equiv-2\,\mathrm{mod}\,5$
averaged over $1000$ realisations. Points are values of $K\big{(}n/N\big{)}$
for integers $n$ with indicated residues modulo $5$. Thin solid lines are
guide for the eye. Thick solid lines indicate the average over all 5 residues:
$\frac{1}{5}\sum_{r=0}^{4}K\big{(}(qt+r)/n\big{)}$.

For illustration, the results of the direct calculation of the form factor for
$\alpha=1/5$ and $N=399$ and $N=398$ are presented in figure 2. First,
eigenvalues of the matrix (17) were calculated numerically and then using (10)
the form factor for different $n$ was computed. The result is averaged over
$1000$ realisations of random phases. It is clear that, indeed, for different
residues of $n$ modulo $5$ the results are different and when $n\equiv 0$ mod
$5$ the form factor at small argument is close to 1 but for all other residues
it starts at $0$. The average over all 5 residues begins at $1/5$ as expected.

Such clear picture appears when the form factor is calculated at special
values of $\tau$, $\tau=n/N$ with integer $n$. Computing it at arbitrary
arguments leads to an irregular plot but, of course, the average curve remains
unchanged.

Exactly the same formulas can be applied for an irrational value of parameter
$\alpha$. In this case one has

$K(\tau)=\tau+2\,\tau\,\mathrm{Re}\left(\frac{e^{-\zeta\tau}}{1-e^{-\zeta\tau}}\right),\qquad\zeta=1-e^{-2\pi
i\alpha N}-2\pi i\alpha N$ (24)

The exponent $\zeta=2\sin^{2}(\pi\alpha N)+i\big{[}\sin(2\pi\alpha
N)-2\pi\alpha N\big{]}$ has a large imaginary part when $N\to\infty$. It means
that the above expression is a strongly oscillated function of $\tau$. When
averaged over a small interval of $\tau$ one obtains $K(\tau)=\tau$ as it
should be for the ensemble of usual random matrices (GUE). This result follows
without calculations from the fact that the average of all eigenvalues
$\Lambda_{\beta}^{\tau N}$ except $\beta=0$ equals zero as a consequence of
rapidly changing phases. (For even $N$ the term with $\beta=N/2$ is real but
as it tends to zero at large $N$ its contribution is negligible.)

Notice that criterion (16) for matrix (17) with irrational $\alpha$ is not
fulfilled. Nevertheless the spectral statistics of such matrix is close to GUE
statistics. This example illustrates a new mechanism for the appearance of
random matrix statistics. The contribution of higher eigenvalues of the
transition matrix (13) decreases not because a gap between the first and the
second eigenvalues as has been proposed in (16) but due to rapid oscillations
for large matrix dimensions.

##  IV Symmetric barrier billiard with $\mathbf{h/b=1/2}$

The central problem of the paper is the determination of level compressibility
for the $B$-matrices given by (1) and (2) by employing the method proposed in
tanner  and used in the previous Section for matrices derived from the
quantisation of an interval-exchange map. The simplicity of treatment of
interval-exchange matrices comes from the fact that their transition matrices
are circulant matrices whose eigenvalues are known exactly. For the
$B$-matrices calculations are more complicated as no explicit formulas for
eigenvalues of the corresponding transition matrix.

$T_{m,n}=|B_{m,n}|^{2}=\frac{L_{m}^{2}\,L_{n}^{2}}{(x_{m}+x_{n})^{2}},\qquad
m,n=1,\ldots,N$ (25)

are available.

This section is devoted to the investigation of the $B$-matrix corresponding
to the symmetric barrier billiard with ratio $h/b=1/2$. In this case $q=2$,
$b-h=h$, and the second part of the vector $\vec{x}$ in (5) coincides with the
third one. Now trivial eigenfunctions can be removed by considering a
desymmetrised rectangular billiard with height $h=b/2$ and imposing the
Neumann boundary conditions for negative $x$ and $y=h$. It is equivalent of
dropping the second part of vector (5) and taking coordinates $x_{m}$ as
follows billiard_1

$x_{m}=(-1)^{m+1}b\sqrt{k^{2}-\frac{\pi^{2}m^{2}}{b^{2}}},\qquad
m=1,\ldots,N,\qquad N=\left[\frac{kb}{\pi}\right].$ (26)

Odd (resp., even) indices describe the first (resp., the third) part of vector
(5).

The numerically calculated spectrum of the transition matrix in this case is
presented in figure 3(a).

(a)

(b)

Figure 3: (a) Transition matrix spectrum for symmetric barrier billiard with
$h/b=1/2$, $b=\pi$, $k=1000.5$, and $N=1000$. Insert shows the difference
between true eigenvalues $\Lambda_{\beta}$ and the straight line $2\beta/N-1$.
(b) Form factor for symmetric barrier billiard with $h/b=1/2$, $b=\pi$,
$k=400.5$, and $N=400$. Data are averaged over 1000 realisations of random
phases. The upper dots are $K\big{(}2n/N\big{)}$, the lower dots are
$K\big{(}(2n-1)/N\big{)}$, and the middle dots correspond to
$\Big{(}K\big{(}2n/N\big{)}+K\big{(}(2n-1)/N\big{)}\Big{)}/2$. Solid line is
the semi-Poisson prediction $K(\tau)=(2+\pi^{2}\tau^{2})/(4+\pi^{2}\tau^{2})$
gerland .

To calculate this spectrum (or, at least, the behaviour of largest-moduli
eigenvalues) analytically a kind of paraxial approximation has been developed.
It is based on the fact that the main ingredient of matrices with intermediate
statistics is a linear fall-off of matrix elements from the diagonal levitov ;
altshuler_levitov . In the simplest setting it means that

$M_{m,n}\sim\frac{R_{m,n}}{m-n+\mathrm{const.}},\qquad m,n\gg 1,\qquad
m-n=\mathcal{O}(1).$

Therefore it is natural to assume that the most important contributions come
from the pole terms with $R_{m,n}\approx R_{m,m}$. This type of approximation
can be done directly from the definition (2) as it is demonstrated in Appendix
A. According to these results the $T$-matrix in the paraxial approximation is
a block matrix

$T=\left(\begin{array}[]{cc}0_{\mathrm{o,o}}&t_{\mathrm{o,e}}\\\
t^{T}_{\mathrm{o,e}}&0_{\mathrm{e,e}}\end{array}\right),\qquad
t_{\mathrm{o,e}}\equiv t_{2m-1,2n}=\frac{1}{\pi^{2}(n-m+1/2)^{2}}\,.$ (27)

Here subscripts ’o’ and ’e’ indicate odd and even indices respectively.

It is instructive to get this answer without the knowledge of the exact
$B$-matrix. One can achieved it by using the instantaneous approximation used
in quantum mechanics when the interaction changes suddenly. In optics such
approximation is analogous to the Fraunhofer diffraction. In the barrier
billiard it corresponds to the situation when a wave with large momentum
quickly moving in a channel enters into another channel (cf., figure 1(b)). In
the instantaneous approximation eigenfunctions in the new channel are just re-
expansion of initial eigenfunctions into a complete set of eigenfunctions with
correct boundary conditions inside the final channel.

Consider a normalised wave with the Neumann boundary conditions at $y=h=b/2$
and the Dirichlet ones at $y=0$

$\psi^{(1)}_{2m-1}(x,y)=\frac{2}{\sqrt{b}}\sin\Big{(}\frac{(2m-1)\pi}{b}y\Big{)}\exp\big{(}ip_{2m-1}^{(1)}x\big{)},\qquad
x<0,\qquad 0\leq y\leq b/2$

propagating in the desymmetrised barrier billiard at negative $x$. When it
penetrates into the region of positive $x$ it has to be expanded into correct
waves propagating inside that region

$\psi^{(1)}_{2m-1}(x,y)=\sum_{n=1}^{\infty}S_{2m-1,2n}\psi^{(3)}_{2n}(x,y)$

where $\psi^{(3)}_{2n}(x,y)$ are waves obeying the Dirichlet boundary
conditions at $y=0$ and $y=h=b/2$

$\psi^{(3)}_{2n}(x,y)=\frac{2}{\sqrt{b}}\sin\Big{(}\frac{2\pi
n}{b}y\Big{)}\exp\big{(}ip_{n}^{(3)}x\big{)},\qquad x>0,\qquad 0\leq y\leq
b/2\,.$

Coefficients $S_{2m-1,2n}$ are the $S$-matrix for this process. In the
paraxial approximation they are calculated as follows (notice that in the
paraxial approximation $p_{2m-1}^{(1)}\approx p_{n}^{(3)}$)

$S_{2m-1,2n}^{(\mathrm{p.a.})}=\frac{4}{b}\int_{0}^{b/2}\sin\Big{(}\frac{\pi(2m-1)}{b}\Big{)}\sin\Big{(}\frac{2\pi
n}{b}y\Big{)}dy=\frac{(-1)^{m+n}}{\pi^{2}(n-m+1/2)}+\frac{(-1)^{m+n}}{\pi^{2}(n+m-1/2)}\,.$

Taking into account only the pole term (and symmetry of the $S$-matrix) one
obtains for the $T$-matrix exactly the same expression as (52).

Thus the transition matrix (27) is a block Toeplitz matrix. It is plain that
its eigenvalues $\Lambda=\pm\sqrt{\lambda_{\beta}}$ where $\lambda_{\beta}$
are eigenvalues of a $N_{1}\times N_{1}$ matrix (with $N_{1}=[N/2]$)

$P_{m,n}=\sum_{k=1}^{N_{1}}t_{m,k}t_{n,k},\qquad
t_{m,k}=\frac{1}{\pi^{2}(k-m+1/2)^{2}}\,.$

Dominant contributions to the sum come from regions $k\sim m$ and $k\sim n$.
Due to a quick decrease of the summands the finite summation over $k$ can
safely be substitute in the limit $N_{1}\to\infty$ by the sum over all integer
$k$

$P_{m,n}\approx\sum_{k=-\infty}^{\infty}t_{m,k}t_{n,k}\,.$

Using (58) the necessary sum is easily calculated and the result is

$P_{m,n}=t_{m-n},\qquad t_{0}=\frac{1}{3},\qquad t_{r\neq
0}=\frac{2}{\pi^{2}r^{2}}\,.$ (28)

This formula is valid when $m,n\gg 1$ and $m-n=\mathcal{O}(1)$.

Matrix (28) is a Toeplitz matrix with quickly decreasing matrix elements. It
is well known that eigenvalues of $N\times N$ Toeplitz matrix can be
asymptotically calculated as follows (see, e.g.,szego -rambour  and references
therein)

$\lambda_{\beta}=f\left(\frac{\beta}{2N}+\mathcal{O}\left(\frac{1}{N}\right)\right),\qquad\beta=1,\ldots,N$
(29)

where function $f(x)$ called the symbol is the Fourier series of $t_{r}$

$f(x)=\sum_{r=-\infty}^{\infty}t_{r}e^{2\pi irx}.$

(More precise formulas can be found in the above references.)

Using (56) and (57) one finds that the symbol of matrix (28) is

$f(x)=\frac{1}{3}+\frac{2}{\pi^{2}}\sideset{}{{}^{\prime}}{\sum}_{r=-\infty}^{\infty}\frac{e^{2\pi
irx}}{r^{2}}=\frac{1}{3}+4B_{2}(x)=(1-2x)^{2}.$

Therefore eigenvalues of the $P$-matrix for large $N_{1}$ are

$\lambda_{\beta}\approx\left(1-\frac{\beta}{N_{1}}\right)^{2},\qquad\beta=1,\ldots,N_{1}\,.$

Eigenvalues of block matrix (27) $\Lambda=\pm\sqrt{\lambda_{\beta}}$. Taking
into account that the dimension of matrix (27) is $N\approx 2N_{1}$ one
concludes that approximately its eigenvalues are

$\Lambda=\pm\left(1-\frac{2\beta}{N}\right),\qquad\beta=1,\ldots,\tfrac{1}{2}N$
(30)

With the corresponding redefinition of index $\beta$ these eigenvalues can be
rewritten in the form

$\Lambda_{\beta}\approx\frac{2\beta}{N}-1,\qquad\beta=1,\ldots,N$

which agrees well with numerical calculations (see figure 3(a)).

The form factor in the diagonal approximation is related with transition
matrix eigenvalues by (14)

$K^{\mathrm{(diag)}}\left(\frac{n}{N}\right)=\frac{n}{N}\sum_{\beta}\Lambda_{\beta}^{n}\,.$

As $\Lambda_{N-\beta}=-\Lambda_{\beta}$ for $\beta=1,\ldots,N-1$ (which is a
consequence of the block structure of the transition matrix (27)) the form
factor $K^{\mathrm{diag}}(n/N)$ with odd $n$ in the diagonal approximation
tends to zero when $\tau=n/N\to 0$

$K^{(\mathrm{diag})}(\tau)\underset{\tau\to 0}{\longrightarrow}0\,.$ (31)

But for even $n$ one gets a different answer. Eq. (30) may not be accurate for
extreme eigenvalues with small $\beta$. For $\tau=2n/N$ one can separate
contribution of small $\beta<\beta_{0}$ and the rest for which (30) is a good
approximation

$K^{(\mathrm{diag})}(\tau)=\tau\left(\mathrm{const}+2\sum_{\beta=\beta_{0}}^{N/2}\left(1-\frac{2\beta}{N}\right)^{\tau
N}\right)\underset{N\to\infty}{\longrightarrow}\tau\left(\mathrm{const}+\frac{2e^{-2j_{0}\tau}}{1-e^{-2\tau}}\right)\underset{\tau\to
0}{\longrightarrow}1\,.$ (32)

As has been discussed in the previous Section it means that the spectral
compressibility of the $B$-matrix for symmetric barrier billiard coincides
with the semi-Poisson value

$\chi\equiv K(0)=\frac{1}{2}\,.$ (33)

For illustration, the form factor for the symmetric billiard calculated
numerically by direct diagonalisation of $400\times 400$ matrices (1) with
coordinates given by (26) and averaged over $1000$ realisations is shown in
figure 3(b). Two branches corresponding to odd and even $n$ are clearly seen.
The average over odd and even values agrees well with the semi-Poisson
expression for the form factor gerland  and, in fortiori, the level
compressibility is $1/2$ as in (33).

##  V Barrier billiard with irrational ratio $\mathbf{h/b}$

The transition matrices for general barrier billiard with off-centre barrier
remain the same as in (25) but coordinates $x_{m}$ should have the form (3)
for irrational ratio $h/b$ and (5) for rational $h/b=m/q$. The direct
calculations of eigenvalues of these matrices reveal that they are more
complicated that the ones for symmetric billiard with $h/b=1/2$ discussed in
the previous Section. As an example, in figure 4 the spectra of the transition
matrices with $h/b=1/\sqrt{5}$ and $h/b=2/5$ are presented. It is clearly seen
that, though eigenvalues with small moduli are quite irregular and have gaps,
largest moduli eigenvalues are well described by a straight line
$\Lambda_{\beta}=2\beta/N-1$.

(a)

(b)

Figure 4: Spectra of the transition matrix for barrier billiard with (a)
$h/b=1/\sqrt{5}$, $b=\pi$, $k=500.5$, and $N=999$ (cf. (4)) (b) $h/b=2/5$,
$b=\pi$, $k=650.5$,and $N=1040$ (cf. (6)). Straight dashed line in the both
figures is $\Lambda_{\beta}=2\beta/N-1$.

This Section is concentrated on the analytical treatment of billiards with
irrational ratio $h/b$. As in the previous Section the first step consists in
the calculation of paraxial $S$-matrix for the scattering inside the slab with
a barrier as in figure 1(b). It can easily be done in the instantaneous
approximation exactly as above. In such approximation only transitions from
channel $1$ to $2$ and to $3$ and their inverse are non-zero. One has

$\displaystyle S_{m,n}^{1\to 3}$ $\displaystyle=$
$\displaystyle\frac{2}{\sqrt{bh}}\int_{0}^{h}\sin\Big{(}\frac{\pi
m}{b}y\Big{)}\sin\Big{(}\frac{\pi n}{h}y\Big{)}dy$ $\displaystyle=$
$\displaystyle\frac{(-1)^{n}\sin(\pi
hm/b)}{\pi\sqrt{bh}}\left((m/b-n/h)^{-1}-(m/b+n/h)^{-1}\right).$

Similarly

$\displaystyle S_{m,n}^{1\to 2}$ $\displaystyle=$
$\displaystyle\frac{2}{\sqrt{b(b-h)}}\int_{h}^{b}\sin\Big{(}\frac{\pi
m}{b}y\Big{)}\sin\Big{(}\frac{\pi n}{b-h}(b-y)\Big{)}dy$ $\displaystyle=$
$\displaystyle\frac{(-1)^{n}\sin(\pi
mh/b)}{\pi\sqrt{b(b-h)}}\left((m/b-n/(b-h))^{-1}-(m/b+n/(b-h))^{-1}\right)\,.$

The transition matrix $T=|S|^{2}$ also has the same block structure. Retaining
only the pole (the first) terms (and slightly changing the notations) one
obtains

$T=\left(\begin{array}[]{ccc}0_{1,1}&t^{1\to 2}_{n_{1},n_{2}}&t^{1\to
3}_{n_{1},n_{3}}\\\ t^{1\to 2}_{n_{1},n_{2}}&0_{2,2}&0_{2,3}\\\ t^{1\to
3}_{n_{1},n_{3}}&0_{3,2}&0_{3,3}\end{array}\right)$ (34)

where $n_{i}=1,\ldots,N_{i}$ with $N_{i}$ given by (4) and

$t^{1\to 2}_{n_{1},n_{2}}=\frac{z_{2}\sin^{2}(\pi
n_{1}h/b)}{\pi^{2}\big{(}z_{2}n_{1}-n_{2}\big{)}^{2}},\qquad t^{1\to
3}_{n_{1},n_{3}}=\frac{z\sin^{2}(\pi
n_{1}h/b)}{\pi^{2}\big{(}zn_{1}-n_{3}\big{)}^{2}},\qquad z_{2}=1-z,\qquad
z=\frac{h}{b}.$ (35)

Due to the block structure of the transition matrix (34) it follows that its
eigenvalues $\Lambda$ are determined by the relation
$\Lambda^{2}=\lambda_{\beta}$ where $\lambda_{\beta}$ are eigenvalues of
$N_{1}\times N_{1}$ matrix

$P_{m,n}=\sum_{k=1}^{N_{2}}t^{1\to 2}_{m,k}t^{1\to
2}_{n,k}+\sum_{k=1}^{N_{3}}t^{1\to 3}_{m,k}t^{1\to 3}_{n,k}$

For large matrix dimension the summation can be extended over all integer $k$
and the sums can be calculated explicitly by using (59) from Appendix C. The
results are

$P_{m,m}=\frac{z^{2}+(1-z)^{2}}{3}(3-2\sin^{2}(\pi mz)),$ (36)

and for $m\neq n$

$P_{m,n}=\frac{2\big{(}\sin^{2}(\pi zm)+\sin^{2}(\pi
nz)\big{)}}{\pi^{2}(m-n)^{2}}-\frac{2(1-2z)\sin(\pi mz)\sin(\pi
nz)\sin(\pi(m-n)z)}{\pi^{3}z(1-z)(m-n)^{3}}\,.$ (37)

This matrix is a combination of Toeplitz terms depended on the difference
$m-n$ and oscillating terms (which explains the existence of forbidden zones
in its spectrum, see figure 4(a).

Due to the unitarity of the $B$-matrix the exact transition matrix $T=|B|^{2}$
has the largest eigenvalue equals $1$ whose corresponding eigenvector is
$(1,1,\ldots,1)$. It is natural (and is confirmed by calculations) that
eigenvectors of the $P$-matrix corresponding to large moduli eigenvalues are
slowly varying functions. Consequently, all oscillating terms in (36) and (37)
for large $m$ and $n$ could be ignored. These arguments lead to the following
recipe of the next step of approximation. Put $m=n+r$ and average all matrix
elements of the $P$-matrix over quickly changing phase $\pi nz$. The
calculations are straightforward and

$\langle{P}_{m,n}\rangle\equiv\lim_{N\to\infty}\frac{1}{N}\sum_{n=1}^{N}P_{n+r,n}=f_{m-n}$
(38)

where

$f_{0}=\frac{2(z^{2}+(1-z)^{2})}{3},\qquad f_{r\neq
0}=\frac{2}{\pi^{2}r^{2}}-\frac{(1-2z)\sin(2\pi rz)}{2\pi^{3}z(1-z)r^{3}}\,.$

Eigenvalues of such matrix for large $N$ are calculated by the Fourier
transform of this symbol

$\lambda_{\beta}=\sum_{r=-\infty}^{\infty}f_{r}e^{2\pi irx},\qquad
x=\frac{\beta}{2N_{1}}\,.$

The necessary sums are expressed through the Bernoulli polynomials (56), (57)
and the result is

$\lambda_{\beta}=\frac{2}{3}(z^{2}+(1-z)^{2})+4B_{2}(x)-\frac{1-2z}{3z(1-z)}\Big{(}B_{3}(\\{x+z\\})-B_{3}(\\{x-z\\})\Big{)}.$

From the beginning one can assume that $h<h-b$, i.e., $z=h/b<1/2$ (the case
$h/b=1/2$ was discussed in Section IV ). Then

$\lambda_{\beta}=\left\\{\begin{array}[]{cc}1-4x+\dfrac{x^{2}}{z(1-z)},&0\leq
x\leq z\\\ \dfrac{1-2z+2z^{2}}{1-z}-\dfrac{2x(1-x)}{1-z},&z\leq
x\leq\frac{1}{2}\end{array}\right.,\qquad x=\frac{\beta}{2N_{1}}\,.$ (39)

As eigenvalues of the block matrix (34) $\Lambda=\pm\sqrt{\lambda_{\beta}}$ it
follows that close to maximum value (i.e., with small $\beta$)

$\Lambda\approx\pm\left(1-\frac{\beta}{N_{1}}+\mathcal{O}\left(\frac{1}{N_{1}}\right)\right)\approx\pm\left(1-\frac{2\beta}{N}\right)$

As in the calculation of the form factor (14) small moduli eigenvalues are
irrelevant one can ignore higher order terms in the above expression which
gives the same expression as in (30). It means that the level compressibility
of barrier billiards with irrational ratio $h/b$ has the same value as in the
preceding Sections $\chi=\tfrac{1}{2}$.

In figure 5(a) the above formulas are compared with the results of direct
calculations for the $P$-matrix with $h/b=1/\sqrt{5}$. As has been
demonstrated, approximate expression (39) is tangent to the exact spectrum
close to $1$. The form factor computed numerically for the same ratio $h/b$ is
presented in figure 28(b). The agreement with the above result is clearly
seen.

(a)

(b)

Figure 5: (a) Eigenvalues of the $P$-matrix (36), (37) for for
$h/b=1/\sqrt{5}$ for the same parameters as in figure 3(a). Solid line
indicates approximate expressions (39). (b) Form factor for $h/b=1/\sqrt{5}$,
$b=\pi$, $k=200.5$, $N=399$ averaged over 1000 realisations. Other notations
are as in figure 3(b).

##  VI Barrier billiard with rational ratio $\mathbf{h/b=p/q}$

The calculation of transition matrix eigenvalues when the ratio $h/b$ is a
rational number can be done by a similar method. An additional difficulty in
such case is that one has to select special combinations of states in the
second and the third channels to remove trivial eigenvalues equal zero on the
whole line passing through the barrier. It has been discussed in detail in
billiard_2  and briefly reviewed in Appendix B. Combining all terms together
one concludes that the transition matrix when $h/b=p/q$ with $p$ and $q$ being
co-prime integers has the block form similar to (34) but with one more block

$T=\left(\begin{array}[]{cccc}0_{1,1}&t^{1\to 2}_{n_{1},n_{2}}&t^{1\to
3}_{n_{1},n_{3}}&t^{1\to 4}_{n_{1},n4}\\\ t^{1\to
2}_{n_{1},n_{2}}&0_{2,2}&0_{2,3}&0_{2,4}\\\ t^{1\to
3}_{n_{1},n_{3}}&0_{3,2}&0_{3,3}&0_{3,4}\\\ t^{1\to
4}_{n_{1},n_{4}}&0_{4,2}&0_{4,3}&0_{4,4}\end{array}\right).$ (40)

Here indices $n_{j}=1,\ldots N_{j}^{\mathrm{(r)}}$ have the following
restrictions

$n_{1}\not\equiv 0\mod q,\qquad n_{2}\not\equiv 0\mod q-p,\qquad
n_{3}\not\equiv 0\mod p$

and

$N_{1}^{\mathrm{(r)}}=N_{1}-N_{0},\qquad
N_{2}^{\mathrm{(r)}}=N_{2}-N_{0},\qquad
N_{3}^{\mathrm{(r)}}=N_{3}-N_{0},\qquad N_{4}^{\mathrm{(r)}}=N_{0}$ (41)

with $N_{1},N_{2},N_{3}$ given by (4) and $N_{0}$ is determined by (54) or
(6). The total matrix dimension is
$N^{(\mathrm{r})}=\sum_{j=1}^{4}N_{j}^{\mathrm{(r)}}=N_{1}+N_{2}+N_{3}-2N_{0}$
as in (6).

Matrices $t^{1\to 2}$ and $t^{1\to 3}$ are the same as in (35) and $t^{1\to
4}$ given by (55) from Appendix B

$t_{n_{1},n_{4}}^{1\to 4}=\frac{\sin^{2}(\pi
mh/b)}{\pi^{2}p(q-p)(n_{1}/q-n_{4})^{2}}\,.$

The eigenvalues of block matrix (40) $\Lambda$ are
$\Lambda=\lambda_{\beta}^{2}$ where $\lambda_{\beta}$ with
$\beta=1,\ldots,N_{1}$ are eigenvalues of matrix (superscript (res) indicates
that the matrix describes the resonance case $h/b=p/q$)

$P_{m,n}^{\mathrm{(res)}}=\sum_{k\neq 0\mod(q-p)}t_{m,k}^{1\to 2}t_{n,k}^{1\to
2}+\sum_{k\neq 0\mod p}t_{m,k}^{1\to 3}t_{n,k}^{1\to 3}+\sum_{k}t_{m,k}^{1\to
4}t_{n,k}^{1\to 4}\,.$

Using an evident relation

$\sum_{k\neq 0\mod r}f(k)=\sum_{k}f(k)-\sum_{k}f(rk)$

and (59) the above sums can be explicitly calculated.

The results are

$P_{m,m}^{\mathrm{(res)}}=\frac{1}{3q^{2}}\big{(}p^{2}+(q-p)^{2}\big{)}\big{(}3-2\sin^{2}(\pi
pm/q)\big{)}+\frac{2\sin^{4}(\pi pm/q)\big{(}3-2\sin^{2}(\pi
m/q)\big{)}}{3p(q-p)q^{2}\sin^{4}(\pi m/q)}$ (42)

and when $m\neq n$

$P_{m,n}^{\mathrm{(res)}}=\frac{b_{m,n}}{(m-n)^{2}}-\frac{c_{m,n}\sin(\pi
p(m-n)/q)}{(m-n)^{3}}-\frac{d_{m,n}\sin(\pi(m-n)/q)}{(m-n)^{3}}$ (43)

where

$\displaystyle b_{m,n}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi^{2}}\big{(}\sin^{2}(\pi pm/q)+\sin^{2}(\pi
pn/q)\big{)}$ (44) $\displaystyle+$ $\displaystyle\frac{2\sin^{2}(\pi
pm/q)\sin^{2}(\pi pn/q)}{p(q-p)\pi^{2}}\Big{(}\frac{1}{\sin^{2}(\pi
m/q)}+\frac{1}{\sin^{2}(\pi n/q)}\Big{)}$ $\displaystyle c_{m,n}$
$\displaystyle=$ $\displaystyle\frac{2q(q-2p)}{\pi^{3}p(q-p)}\sin(\pi
pm/q)\sin(\pi pn/q),\qquad d_{m,n}=\frac{4q\sin^{2}(\pi pm/q)\sin^{2}(\pi
pn/q)}{\pi^{3}p(q-p)\sin(\pi m/q)\sin(\pi n/q)}\,.$

Though these expressions are indexed by integers $m$ and $n$ this notation is
symbolic. The point is that by construction these integers cannot be arbitrary
but have to be not divisible by $q$. Let us ordered such numbers and let
$\nu(k)$ with $k=1,2,\ldots,$ be the $k^{\mathrm{th}}$ integer $\not\equiv
0\mod q$. Then indices of matrix $P_{m,n}^{\mathrm{(res)}}$ have to be
considered as follows: $m=\nu(j)$, $n=\nu(k)$ with
$j,k=1,2,\ldots,N_{1}^{\mathrm{(r)}}$ with $N_{1}^{\mathrm{(r)}}$ defined in
(41). In such notation matrix $P^{\mathrm{(res)}}$ is
$N_{1}^{\mathrm{(r)}}\times N_{1}^{\mathrm{(r)}}$ matrix

$P_{m,n}^{\mathrm{(res)}}\equiv P_{\nu(j),\nu(k)}^{\mathrm{(res)}},\qquad
j,k=1,\ldots,N_{1}^{\mathrm{(r)}}.$

The next step, as in the previous Section (cf., (38)), consists in the
calculation instead of the above exact expressions their mean values with
fixed difference between the indices

$\langle
P_{m,n}^{\mathrm{(res)}}\rangle(r)=\lim_{N\to\infty}\frac{1}{N}\sum_{n=1}^{N}P_{n+r,n}^{\mathrm{(res)}}$

where the both integers $n$ and $n+r$ have to be not divisible by $q$.

According to (42) and (43) the $P^{\mathrm{(res)}}$ matrix is a mixture of
functions depending explicitly on the differences of indices and certain
coefficients depending on indices modulus $q$. Only the latter requires the
explicit averaging. Using (60)-(63) from Appendix C one obtains that

$\displaystyle\sum_{m=1}^{q-1}P_{m,m}=\frac{(p^{2}+(q-p)^{2})(2q-3)}{3q^{2}}+\frac{2p^{2}(2q-3p)}{3(q-p)q^{2}},$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{m=1}^{q-1}b_{m,m+r}=\frac{2q(q-p+1-2\sin^{2}(\pi
pr/q))}{\pi^{2}(q-p)},$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{m=1}^{q-1}c_{m,m+r}=\frac{q^{2}(q-2p)\cos(\pi
pr/q)}{\pi^{3}p(q-p)},\qquad\sideset{}{{}^{\prime}}{\sum}_{m=1}^{q-1}d_{m,m+r}=\frac{q^{2}\sin(2\pi
pr/q)}{\pi^{3}p(q-p)\sin(\pi r/q)}\,.$

Here it is taken into account that $p/q<1/2$. The superscript $\prime$ in
these sums indicates that the term with $m+r\equiv 0\mod q$ is omitted. The
latter condition implies that the number of independent terms equal $q-1$ if
$r\equiv 0\mod q$ or $q-2$ otherwise. Finally one obtains

$\langle P_{m,n}^{\mathrm{(res)}}\rangle=\alpha_{1}\delta_{m,n}+f_{m-n}$ (45)

with

$f_{m-n}=\left\\{\begin{array}[]{cc}\dfrac{\alpha_{2}}{(m-n)^{2}},&(m-n)\equiv
0\mod q\\\ &\\\ \dfrac{\alpha_{3}}{(m-n)^{2}}+\dfrac{\alpha_{4}\cos(2\pi
pr/q)}{(m-n)^{2}}+\dfrac{\alpha_{5}\sin(2\pi pr/q)}{(m-n)^{3}},&(m-n)\equiv
r\neq 0\mod q\end{array}\right.$ (46)

where constants $\alpha_{j}$ are

$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle\frac{(p^{2}+(q-p)^{2})(2q-3)}{3q^{2}(q-1)}+\frac{2p^{2}(2q-3p)}{3(q-p)q^{2}(q-1)}\qquad\alpha_{2}=\frac{2q(q-p+1)}{\pi^{2}(q-p)(q-1)},$
$\displaystyle\alpha_{3}$ $\displaystyle=$
$\displaystyle\frac{2q}{\pi^{2}(q-2)},\qquad\alpha_{4}=\frac{2q}{\pi^{2}(q-p)(q-2)},\qquad\alpha_{5}=\frac{q^{2}(2p-q-2)}{2\pi^{3}p(q-p)(q-2))}\,.$
(47)

Though this matrix depends only on the difference of indices $m-n$ it is not a
Toeplitz matrix as $m$ and $n$ are not arbitrary numbers but only integers not
divisible by $q$. Nevertheless one can argue that largest eigenvalues for
large matrix dimension can be calculated by a formula similar to Toeplitz
matrices (which is a kind of variational method)

$\lambda_{\beta}=\frac{1}{N_{1}^{\mathrm{(r)}}}\sum_{j,k=1}^{N_{1}^{\mathrm{(r)}}}\langle
P_{\nu(j),\nu(k)}^{\mathrm{(res)}}\rangle e^{2\pi ix(j-k)},\qquad
x=\frac{\beta}{2N_{1}^{\mathrm{(r)}}}\,.$ (48)

Here, as above, $\nu(k)$ is the $k^{\mathrm{th}}$ integer $\not\equiv 0\mod
q$.

In Appendix C (see (65)) it is shown that such sum can be written as follows

$\displaystyle\lambda_{\beta}$ $\displaystyle=$
$\displaystyle\alpha_{1}+\alpha_{2}\sum_{u=-\infty}^{\infty}\frac{e^{2\pi
ix(q-1)u}}{(qu)^{2}}+2\mathrm{Re}\,\sum_{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\times$
$\displaystyle\times$ $\displaystyle\sum_{u=-\infty}^{\infty}e^{2\pi
ix((q-1)u+t)}\left(\frac{\alpha_{3}}{(uq+t)^{2}}+\frac{\alpha_{4}\cos(2\pi
pt/q)}{(uq+t)^{2}}+\frac{\alpha_{5}\sin(2\pi pr/q)}{(uq+t)^{3}}\right)\,.$

The first sum is calculated through the Bernoulli polynomial $B_{2}(x)$ (see
(56)). The last sums are expressed through two functions

$G(x,r)=\sum_{k=-\infty}^{\infty}\frac{e^{2\pi ix(kq+r)}}{(kq+r)^{2}},\qquad
F(x,r)=\sum_{k=-\infty}^{\infty}\frac{e^{2\pi ix(kq+r)}}{(kq+r)^{3}}\,.$

The explicit expressions of these function can be obtained as follows.

Define one more function

$g(x,r)=\sum_{k=-\infty}^{\infty}\frac{e^{2\pi ix(kq+r)}}{kq+r}\,.$

By the differentiation over $x$ one has $G^{\prime}(x,r)=2\pi ig(x,r)$ and
$F^{\prime}(x,r)=2\pi iG(x,r)$. As the differentiation of $g(x,r)$ over $x$
gives the sum of $\delta$-function it is plain that $g(x,r)$ is the piece-wise
constant function in interval $[j/q,(j+1)/q]$. Using (58) one gets

$g(j/q+x,r)=\frac{\pi}{q\sin(\pi r/q)}e^{i\pi r(2j+1)/q}\,.$

Correspondingly, function $G(x,r)$ is a piece-wise linear function in the same
intervals

$G(j/q+x,r)=\frac{\pi^{2}}{q^{2}\sin^{2}(\pi r/q)}e^{2\pi
irj/q}+\frac{2\pi^{2}ix}{q\sin(\pi r/q)}e^{i\pi r(2j+1)/q}\,.$ (49)

In the same way one proves that function $F(x,r)$ is a piece-wise quadratic
function

$F(j/q+x,r)=\frac{\pi^{3}\cos(\pi r/q)}{q^{3}\sin^{3}(\pi r/q)}e^{2\pi
irj/q}+\frac{2\pi^{3}ix}{q^{2}\sin^{2}(\pi r/q)}e^{2\pi
irj/q}-\frac{2\pi^{3}x^{2}}{q\sin(\pi r/q)}e^{i\pi r(2j+1)/q}\,.$ (50)

In all these formulas $j=0,\ldots,q-1$ and $0\leq x\leq 1/q$.

Combining all terms together one finds

$\displaystyle\lambda_{\beta}=\alpha_{1}+\frac{2\pi^{2}\alpha_{2}}{q^{2}}B_{2}\big{(}\big{\\{}x(q-1)\big{\\}}\big{)}+2\mathrm{Re}\,\sum_{t=1}^{q-1}\Big{(}1-\frac{t}{q-1}\Big{)}\times$
(51) $\displaystyle\times$ $\displaystyle
e^{ixt/q}\left[\big{(}\alpha_{3}+\alpha_{4}\cos(2\pi
pt/q)\big{)}G\Big{(}\frac{x(q-1)}{q},t\Big{)}+\alpha_{5}\sin(2\pi
pt/q)F\Big{(}\frac{x(q-1)}{q},t\Big{)}\right]\,.$

The main interest for the calculation of the form factor is the behaviour of
the largest eigenvalues for $x$ close to zero. Using (49) and (50) one
concludes that

$\lambda_{\beta}=C_{0}+C_{1}\frac{\beta}{N_{1}^{\mathrm{(r)}}},\qquad
0\leq\beta\ll N_{1}^{\mathrm{(r)}}\,.$

Here

$C_{0}=\alpha_{1}+\frac{\pi^{2}\alpha_{2}}{3q^{2}}+\frac{2\pi^{2}}{q^{2}}\sum_{t=1}^{q-1}\Big{(}1-\frac{t}{q-1}\Big{)}\Big{(}\frac{\alpha_{3}+\alpha_{4}\cos(2\pi
pt/q)}{\sin^{2}(\pi t/q)}+\frac{\pi\alpha_{5}\sin(2\pi pr/q)\cos(\pi
t/q)}{q\sin^{3}(\pi t/q)}\Big{)}$

and

$C_{1}=-\frac{2\pi^{2}\alpha_{2}(q-1)}{2q^{2}}-\frac{2\pi^{2}(q-1)}{q^{2}}\sum_{t=1}^{q-1}\Big{(}1-\frac{t}{q-1}\Big{)}\big{(}\alpha_{3}+\alpha_{4}\cos(2\pi
pt/q)\big{)}\,.$

The sum over residues is of the form

$\sum_{t=1}^{q-1}\Big{(}1-\frac{t}{q-1}\Big{)}h(t)$

and (as it is easy to check) in the considered case $h(q-t)=h(t)$. Therefore

$\sum_{t=1}^{q-1}t\,h(t)=\sum_{t=1}^{q-1}(q-t)\,h(q-t)=\sum_{t=1}^{q-1}(q-t)\,h(t)\longrightarrow\sum_{t=1}^{q-1}t\,h(t)=\frac{q}{2}\sum_{t=1}^{q-1}h(t)\,.$

Consequently

$C_{0}=\alpha_{1}+\frac{\pi^{2}\alpha_{2}}{3q^{2}}+\frac{2\pi^{2}(q-2)}{q^{2}(q-1)}\sum_{t=1}^{q-1}\Big{(}\frac{\alpha_{3}+\alpha_{4}\cos(2\pi
pt/q)}{\sin^{2}(\pi t/q)}+\frac{\pi\alpha_{5}\sin(2\pi pr/q)\cos(\pi
t/q)}{q\sin^{3}(\pi t/q)}\Big{)}$

and

$C_{1}=-\frac{2\pi^{2}\alpha_{2}(q-1)}{2q^{2}}-\frac{2\pi^{2}(q-2)}{q^{2}}\sum_{t=1}^{q-1}\big{(}\alpha_{3}+\alpha_{4}\cos(2\pi
pt/q)\big{)}\,.$

Using sums indicated in Appendix C and collecting all terms in the end one
finds that

$C_{0}=1,\qquad C_{1}=-2.$

This result signifies that largest moduli eigenvalues of the transition matrix
for the barrier billiard with rational ratio $h/b=p/q$ are (i) independent on
values of integers $p$ and $q$ and (ii) have the same asymptotic expression as
in (30) (taking into account that $N^{(\mathrm{r})}\approx
2N_{1}^{\mathrm{(r)}}$)

$\Lambda=\pm\left(1-\frac{\beta}{N_{1}^{\mathrm{(r)}}}+\mathcal{O}\left(\frac{1}{N_{1}^{\mathrm{(r)}}}\right)\right)\approx\pm\left(1-\frac{2\beta}{N^{\mathrm{(r)}}}\right)\,.$

As it has been explained above it implies that (i) the form factor for barrier
billiards is different for odd and even $n$ and (ii) the spectral
compressibility is exactly equal $1/2$ for all positions of the barrier

$\chi\equiv K(0)=\frac{1}{2}\,.$

The numerical calculations exemplified in figure 6 confirm well these results.

(a)

(b)

Figure 6: (a) Eigenvalues of the $P$-matrix (42), (43) for $h/b=2/5$ with
$N_{1}=520$. Solid line is the spectrum (51) of the asymptotic matrix. (b)
Form factor for $h/b=2/5$, $b=\pi$, $k=250.5$,and $N=400$ averaged over 1000
realisations. Other notations are as in figure 3(b).

##  VII Summary

It is demonstrated that the method of calculation of the level compressibility
proposed by G. Tanner in tanner  for chaotic systems can successfully be
applied for intermediate statistics models. The criterium discussed in tanner
states that if the difference between the dominant eigenvalue of the
transition matrix $\Lambda=1$ and the second in magnitude eigenvalue is big
enough then only the dominant eigenvalues contributes to the form factor and
one gets the usual value of the form factor corresponding to standard random
matrix ensembles. Notably, the level compressibility is zero.

For models considered in the paper no individual transition matrix eigenvalues
dominate and one has to sum over many of them with moduli close to $1$. Two
types of random unitary matrices were investigated. The first corresponding to
a quantisation of an interval-exchange map marklof  has been discussed in
detail in map -integrable_ensembles . In particular, the values of the level
compressibility were derived. The application of the transition matrix
approach for this case serves first of all to check the validity of Tanner’s
method for intermediate statistics models. It appears that interval-exchange
matrices lead to circulant transition matrices whose eigenvalues are
explicitly known and all necessary sums are easily estimated. In the end one
gets the same values of the level compressibility as obtained in map
-integrable_ensembles  but with much simpler and transparent calculations. An
example is of a special interest. It corresponds to interval-exchange matrices
with an irrational value of a parameter (which strictly speaking describes not
an interval-exchange map but only a parabolic one). Numerically it has
observed in general_case  that in such case spectral statistics is of usual
random matrix type (GOE or GUE depending on a symmetry) as for chaotic systems
which looks strange as the Lyapunov exponent of any parabolic map is zero. The
transition matrix approach clearly indicates that, though there is no dominant
eigenvalue as been discussed in tanner , all eigenvalues except $\Lambda=1$
for large matrix dimensions are so quickly oscillating that averaging over a
small interval of the argument effectively removes their contributions
producing the standard random matrix result.

The main part of the paper contains the calculation of the level
compressibility for random unitary matrices derived from the exact
quantisation of barrier billiards in billiard_1 ; billiard_2 . The importance
of such matrices comes from the fact that they have the same spectral
statistics as high-excited states of barrier billiards which are the simplest
examples of pseudo-integrable models for which very little is known
analytically.

The barrier billiard transition matrices are more complicated that the ones
for interval-exchange matrices. Their spectra contain forbidden zones and
their exact eigenvalues, seems, not to be accessible in closed form.
Nevertheless, as the level compressibility requires the control only of
largest moduli eigenvalues of the transition matrix it is possible to find
such eigenvalues for large matrix dimensions precisely. The main
simplification comes from the fact that eigenvectors corresponding to largest
moduli eigenvalues are slow oscillating functions. Therefore quickly
oscillating terms in matrix elements will give negligible contributions on
these eigenvectors and one can substitute instead of exact matrix elements
their average over fast oscillations. The resulting matrices are simpler and
permit to find their large moduli eigenvalues analytically. In the end one
proves that the level compressibility of barrier billiards for all positions
and heights of the barrier is the same and equals $1/2$. This result strongly
indicates that spectral statistics of the $B$-matrices associated with barrier
billiards is universal (i.e., independent on the barrier position) and well
described by the semi-Poisson distribution.

##  Appendix A Approximate expression for the $B$-matrix for the symmetric
billiard

The purpose of this Appendix is the determination of the transition matrix for
the symmetric case (i.e., $x_{m}=(-1)^{m+1}y_{m}$,
$y_{m}=b\sqrt{k^{2}-\pi^{2}m^{2}/b^{2}}$) in the paraxial approximation by
taking into account only the pole terms in the definition (1). From (2) with
odd $m=2j-1$ it follows (for simplicity it is assumed $N$ is even $N=2n$ and
the products is taken from $q=1$ till $q=n$)

$\displaystyle 2y_{2j-1}L_{2j-1}^{2}=(2y_{2j-1})^{2}\prod_{q\neq
j}\frac{y_{2j-1}+y_{2q-1}}{y_{2j-1}-y_{2q-1}}\prod_{q}\frac{y_{2j-1}-y_{2q}}{y_{2j-1}+y_{2q}}$
$\displaystyle=$
$\displaystyle\prod_{q}\Big{(}\frac{y_{2j-1}+y_{2q-1}}{y_{2j-1}+y_{2q}}\Big{)}^{2}\frac{\pi^{2}\prod_{q}[q^{2}-(j-1/2)^{2}]}{\prod_{q\neq
j}[(q-1/2)^{2}-(j-1/2)^{2}]}=g_{j}\,F^{2}(y_{2j-1})$

where

$g_{j}=\lim_{v\to
j-1/2}\frac{\pi^{2}((j-1/2)^{2}-v^{2})\prod_{q}[1-(j-1/2)^{2}/q^{2}]}{\prod_{q}[1-v^{2}/(q-1/2)^{2}]},\qquad
F(z)=\prod_{q}\frac{q(z+y_{2q-1})}{(q-1/2)(z+y_{2q})}\,.$

Exactly in the same way one gets

$\displaystyle 2y_{2j}L_{2j}^{2}=-(2y_{2j})^{2}\prod_{q\neq
j}\frac{y_{2j}+y_{2q}}{y_{2j}-y_{2q}}\prod_{q}\frac{y_{2j}-y_{2q-1}}{y_{2j}+y_{2q-1}}$
$\displaystyle=$
$\displaystyle-\prod_{q}\Big{(}\frac{y_{2j}+y_{2q}}{y_{2j}+y_{2q-1}}\Big{)}^{2}\frac{\,\pi^{2}\prod_{q}[(q-1/2)^{2}-j^{2}]}{\prod_{q\neq
j}[q^{2}-j^{2}]}=\frac{f_{j}}{F^{2}(y_{2j})}$

with

$f_{j}=-\lim_{v\to
j}\frac{\pi^{2}(j^{2}-v^{2})\prod_{q}[1-j^{2}/(q-1/2)^{2}]}{\prod_{q}[1-v^{2}/q^{2}]}\,.$

All products in the above expressions should be taken from $1$ to $n$. If $j$
is not too close to $kb/\pi$ (i.e., the momentum is not close to the threshold
of new propagating modes) the products in $g_{j}$ and $f_{j}$ can be extended
to infinity and these functions can easily be calculated from standard
expressions

$\sin\pi x=\pi
x\prod_{q=1}^{\infty}\Big{(}1-\frac{x^{2}}{q^{2}}\Big{)},\qquad\cos\pi
x=\prod_{q=1}^{\infty}\Big{(}1-\frac{x^{2}}{(q-1/2)^{2}}\Big{)}\,.$

In this way one finds that $g_{j}=2$ and $f_{j}=2\pi^{2}j^{2}$. Finally

$L_{2j-1}=\frac{F(y_{2j-1})}{\sqrt{y_{2j-1}}},\qquad L_{2j}=\frac{\pi
j}{F(y_{2j})\sqrt{y_{2j}}}\,.$

The $T$-matrix elements are

$T_{m,n}=\frac{L_{m}^{2}L_{n}^{2}}{(x_{m}+x_{n})^{2}}\,.$

In the paraxial approximation one should take into account only the terms with
$x_{m}$ of different signs. For symmetric billiard it means that
$T_{2n,2m}\approx 0$, $T_{2n-1,2m-1}=0$ and

$T_{2m-1,2n}=\frac{L_{2m-1}^{2}L_{2n}^{2}}{(y_{2m-1}-y_{2n})^{2}}\approx\frac{F^{2}(y_{2m-1})}{F^{2}(y_{2n})}\frac{\pi^{2}n^{2}}{y_{2m-1}y_{2n}(y_{2m-1}-y_{2n})^{2}}\,.$

As only the pole terms are important one can put $y_{2n}\approx y_{2m-1}$ and
$F^{2}(y_{2m-1})\approx F^{2}(y_{2n})$

$\displaystyle T_{2m-1,2n}$ $\displaystyle\approx$
$\displaystyle\frac{\pi^{2}n^{2}}{y_{2m-1}y_{2n}(y_{2m-1}-y_{2n})^{2}}=\frac{n^{2}}{\pi^{2}(n^{2}-(m-1/2)^{2})^{2}}\frac{(y_{2m-1}+y_{2n})^{2}}{y_{2m-1}y_{2n}}$
(52) $\displaystyle\approx$
$\displaystyle\frac{4n^{2}}{\pi^{2}(n^{2}-(m-1/2)^{2})^{2}}=\frac{1}{\pi^{2}(n-m+1/2)^{2}}+\frac{1}{\pi^{2}(n+m-1/2)^{2}}\,.$

##  Appendix B Instantaneous approximation for the resonance case

When the ratio $h/b$ is a rational number $h/b=p/q$ with co-prime integer $p$
and $q$ it is plain that for the barrier billiard as in figure 1(a) the
following 3 transverse momenta with integer $t$ (and the corresponding
longitudinal momenta $p_{t}=\sqrt{k^{2}-\pi^{2}w_{t}^{2}}$) are equal

$w_{t}=\frac{\pi qt}{b}=\frac{\pi((q-p)t}{b-h}=\frac{\pi pt}{h}\,.$ (53)

Introduce the elementary solutions with these momenta in each of 3 regions
indicated in figure 1(b)

$\displaystyle\phi_{qt}^{(1)}(x,y)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2}{b}}\sin\Big{(}\frac{\pi
qt}{b}y\Big{)}\exp\big{(}ip_{t}x\big{)}\,,\qquad 0\leq y\leq b,\qquad x<0\,,$
$\displaystyle\phi_{(q-p)t}^{(2)}(x,y)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2}{b-h}}\sin\Big{(}\frac{\pi(q-p)t}{b-h}(b-y)\Big{)}\exp\big{(}ip_{t}x\big{)}\,,\qquad
b-h\leq y\leq b,\qquad x>0\,,$ $\displaystyle\phi_{pt}^{(3)}(x,y)$
$\displaystyle=$ $\displaystyle\sqrt{\frac{2}{h}}\sin\Big{(}\frac{\pi
pt}{h}y\Big{)}\exp\big{(}ip_{t}x\big{)}\,,\qquad 0\leq y\leq h,\qquad x>0\,.$

Due to the resonant conditions (53) all these solutions represent exact
solutions for the scattering inside the slab in figure 1(b). The number of
such solutions is

$N_{0}=\left[\frac{kb}{\pi q}\right].$ (54)

When spectral statistics of non-trivial eigenvalues is considered these
solutions should be removed. It has been done in detail in billiard_2 . Below
the derivation of the paraxial approximation for the $T$-matrix in such case
is briefly discussed.

The paraxial approximation for $T$-matrix for the scattering inside the slab
in figure 1(b) for non-resonant waves when in the first region $n_{1}\neq 0\
\mathrm{mod}\ q$, in the second region $n_{2}\neq 0\mod q-p$, and in the third
one $n_{3}\neq 0\mod p$ are given by the same expression as in (35). The first
step consists in removing all wave from region $1$ proportional to $q$. But it
is not enough as waves from the second and the third regions can diffract into
waves in the first region with $n_{1}\equiv 0\mod q$. To remove them notice
that

$\phi_{qt}^{(1)}(x,y)=\left\\{\begin{array}[]{cc}-(-1)^{qt}\sqrt{\dfrac{b}{b-h}}\phi_{(q-p)t}^{(2)}(x,y),&h\leq
y\leq b\\\ \sqrt{\dfrac{b}{h}}\phi_{pt}^{(3)}(x,y),&0\leq y\leq
b\end{array}\right..$

Therefore the following linear combination billiard_2

$\phi_{t}^{(4)}(x,y)=\sqrt{\frac{h}{b}}\phi_{(q-p)t}^{(2)}(x,y)+(-1)^{qt}\sqrt{\frac{b-h}{b}}\phi_{(q-p)t}^{(3)}(x,y)$

is orthogonal to $\phi_{qt}^{(1)}(x,y)$ and cancels undesirable waves with
$n_{1}\equiv 0\mod q$.

The calculation of the scattering into such state can be performed as above

$\displaystyle S_{m,t}^{1\to
4}=2\sqrt{\frac{h}{b^{2}(b-h)}}\int_{h}^{b}\sin\Big{(}\frac{\pi
my}{b}\Big{)}\sin\Big{(}\frac{\pi(q-p)t(b-y)}{b-h}\Big{)}dy$ $\displaystyle+$
$\displaystyle
2(-1)^{qt}\sqrt{\frac{b-h}{b^{2}h}}\int_{0}^{h}\sin\Big{(}\frac{\pi
my}{b}\Big{)}\sin\Big{(}\frac{\pi pty}{h}\Big{)}dy$ $\displaystyle=$
$\displaystyle\frac{(-1)^{(p+q)t}\sin\big{(}\pi
mh/b\big{)}b}{\sqrt{h(b-h)}}\Big{(}\frac{1}{m-qt}-\frac{1}{m+qt}\Big{)}\,.$

The paraxial approximation of the corresponding transition matrix elements is
given by the pole term

$t_{m,t}^{1\to 4}=\frac{\sin^{2}(\pi mh/b)}{\pi^{2}p(q-p)(m/q-t)^{2}},\qquad
m\not\equiv 0\ \mod q\,.$ (55)

##  Appendix C Divers relations

In this Appendix a few formulas used in the text are briefly reviewed.

The known identity (see e.g., bateman )

$\sideset{}{{}^{\prime}}{\sum}_{r=-\infty}^{\infty}\frac{1}{r^{n}}e^{2\pi
irx}=-\frac{(2\pi i)^{n}}{n!}B_{n}(\\{x\\})$ (56)

where $\\{x\\}$ is the fractional part of $x$ and $B_{n}(x)$ are Bernoulli
polynomials. For example,

$B_{2}(x)=x^{2}-x+\frac{1}{6}\,,\qquad
B_{3}(x)=x^{3}-\frac{3}{2}x^{2}+\frac{1}{2}x\,.$ (57)

The following identities are standard and presented for completeness.

$\sum_{n=-\infty}^{\infty}\frac{1}{n+x}=\pi\cot(\pi
x),\qquad\sum_{n=-\infty}^{\infty}\frac{1}{(n+x)(n+y)}=\frac{\pi\big{(}\cot(\pi
x)-\cot(\pi y)\big{)}}{y-x}\,.$ (58)

The first formula is simply the expansion of the right-hand side over poles.
The second is a consequence of the first. Differentiating the above
expressions by $x$ and $y$ shows that

$\displaystyle\sum_{n=-\infty}^{\infty}\frac{1}{(n+x)^{2}}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{\sin^{2}(\pi
x)},\qquad\sum_{n=-\infty}^{\infty}\frac{1}{(n+x)^{4}}=\frac{\pi^{4}(3-2\sin^{2}(\pi
x))}{3\sin^{4}(\pi x)},$ (59)
$\displaystyle\sum_{n=-\infty}^{\infty}\frac{1}{(n+x)^{2}(n+y)^{2}}$
$\displaystyle=$
$\displaystyle\frac{\pi^{2}}{(x-y)^{2}}\Big{(}\frac{1}{\sin^{2}(\pi
x)}+\frac{1}{\sin^{2}(\pi y)}\Big{)}+\frac{2\pi(\cot(\pi x)-\cot(\pi
y))}{(x-y)^{3}}\,.$

In different places of the paper one needs to calculate finite sums
$\sum_{t=1}^{q-1}f(t)$ over residues $\mod q$. A characteristic feature of
such sums is that their summands $f(t)$ can be rewritten as ratio of two
polynomials in variable $z=e^{2\pi it/q}$. The summation over $t$ from $1$ to
$q-1$ corresponds to the calculation of the integral

$\sum_{t=1}^{q-1}f(t)=\frac{1}{2\pi i}\oint_{C}\frac{qz^{q-1}}{z^{q}-1}f(z)dz$

where contour $C$ encircled all roots of $z^{q}-1=0$ except the one with
$z=1$. By deforming the contour and calculating the necessary residues one can
obtain the necessary sums in closed form. Below a few formulas obtained by
this manner are listed

$\displaystyle\sum_{t=1}^{q-1}\frac{1}{\sin^{2}(\pi
t/q)}=\frac{q^{2}-1}{3}\,,\qquad\sum_{t=1}^{q-1}\frac{\cos(2\pi
pt/q)}{\sin^{2}(\pi t/q)}=\frac{q^{2}+6p^{2}-1}{3}-2qp\,,$ (60)
$\displaystyle\sum_{t=1}^{q-1}\frac{\sin(2\pi pt/q)\cos(\pi t/q)}{\sin^{3}(\pi
t/q)}=\frac{2p(q^{2}+2p^{2})}{3}-2qp^{2}\,,$ (61)
$\displaystyle\sum_{m=1}^{q-1}\frac{\sin^{2}(\pi pm/q)\sin^{2}(\pi
p(m+r)/q)}{\sin^{2}(\pi m/q)}=\frac{pq}{2}-p^{2}\sin^{2}(\pi pr/q)\,,$ (62)
$\displaystyle\sum_{m=1}^{q-1}\frac{\sin^{4}(\pi pm/q)}{\sin^{4}(\pi
m/q)}=\frac{pq(2p^{2}+1)}{3}-p^{4}\,.$ (63)

Here it is assumed that $p/q<1/2$.

In Section VI one has to calculate the following sum where $\nu(k)$ is the
$k^{\mathrm{th}}$ integer not divisible by $q$

$S=\frac{1}{N}\sum_{j,k=1}^{N}f(\nu(j)-\nu(k))e^{ix(j-k)}$

in the limit $N\to\infty$ with a certain quickly decreasing function $f(x)$,
$f(x)\underset{|x|\to\infty}{\longrightarrow}x^{-2}$.

To get an explicit expression of that sum notice that the number of integers
from $1$ to $n$ divisible by $q$ is $[n/q]$ where $[x]$ is the largest integer
less of equal $n$. Therefore if $n=ql+r$ with $r=1,\ldots,q-1$ then
$k=(q-1)l+r$. It means that

$\nu((q-1)l+r)=ql+r\,,\qquad r=1,\ldots,q-1\,.$ (64)

As integer $\nu(k)\neq 0\mod q$ the residue $r\neq 0$.

Writing $j=(q-1)m+s$ and $k=(q-1)n+r$ with integer $m$, $n$ and
$r,s=1,\ldots,q-1$ one gets that

$j-k=(q-1)u+t,\qquad\nu(j)-\nu(k)=qu+t,\qquad u=m-n,t=s-r\,.$

The summation over integers $j$ and $k$ is equivalent to the summation over
integers $m,n$ and $s,r$. Fixing the differences $k=m-n$ and $t=s-r$, using
the fact that integers with fixed residue $\mod(q-1)$ are uniformly
distributed

$\lim_{N\to\infty}\frac{1}{N}\sum_{j=1}^{N}F(j\mod(q-1))=\frac{1}{q-1}\sum_{s=1}^{q-1}F(s)$

and that

$\sum_{s=1}^{q-1}\sum_{r=1}^{q-1}f(s-r)=\sum_{t=-(q-2)}^{q-2}(q-1-|t|)f(t)$

one finds that

$S=\sum_{u=-\infty}^{\infty}\sum_{t=-(q-2)}^{q-2}\left(1-\frac{|t|}{q-1}\right)f(qu+t)e^{ix((q-1)u+t)}\,.$

Due to a quick decrease of function $f(x)$ the summation over $u$ can be
extended to the sum over all integers. It is convenient to separate term with
$t=0$, add together terms with $t$ and $-t$, and in the last term change
$u\to-u$

$\displaystyle S$ $\displaystyle=$
$\displaystyle\sum_{u=-\infty}^{\infty}f(qu)e^{ix(q-1)u}$ $\displaystyle+$
$\displaystyle\sum_{u=-\infty}^{\infty}\sum_{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)\left(f(qu+t)e^{ix((q-1)u+t)}+f(-qu-t)e^{-ix((q-1)u+t)}\right)\,.$

The used function $f$ is even $f(-x)=f(x)$ and this expression can be written
as follows

$S=\sum_{u=-\infty}^{\infty}f(qu)e^{ix(q-1)u}+2\mathrm{Re}\,\sum_{u=-\infty}^{\infty}\sum_{t=1}^{q-1}\left(1-\frac{t}{q-1}\right)f(qu+t)e^{ix((q-1)u+t)}\,.$
(65)

## References

  * (1) E. P. Wigner, Random matrices in physics, SIAM Review 9, 1 (1967).
  * (2) M. V. Berry and M. Tabor, Level clustering in the regular spectrum, Proc. R. Soc. Lond. A 356, 375, (1977).
  * (3) O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52, 1 (1984).
  * (4) M. L. Mehta, Random matrices, ?Third edition, Academic Press (2014).
  * (5) P.J. Richens and M.V. Berry, Pseudointegrable systems in classical and quantum mechanics, Physica D: Nonlinear Phenomena 2, 495 (1981).
  * (6) A.N. Zemlyakov and A.B. Katok, Topological transitivity in billiards in polygons, Math. Notes 18, 760 (1975).
  * (7) E. Bogomolny, Barrier billiard and random matrices, J. Phys. A: Math. Theor. 55, 024001 (2022).
  * (8) E. Bogomolny, Random matrices associated with general barrier billiards, arXiv:2111.00198 (2021).
  * (9) B. L. Altshuler, I. Kh. Zharekeshev, S. A. Kotochigava, and B.I. Shklovskii, Repulsion between levels and the metal-insulator transition, Sov. Phys. JETP 67, 625 (1988).
  * (10) B.I. Shklovskii, B. Shapiro, B.R. Sears, P. Lambrianides, and H.B. Shore, Statistics of spectra of disordered systems near the metal-insulator transition, Phys. Rev. B 47, 11487 (1993).
  * (11) E. Bogomolny, U. Gerland, and C. Schmit, Short-range plasma model for intermediate spectral statistics, Eur. Phys. J. 19, 121 (2001).
  * (12) M.V. Berry, Semiclassical theory of spectral rigidity, Proc. R. Soc. Lond., A 400, 229 (1985).
  * (13) J. Wiersig, Spectral properties of quantized barrier billiards, Phys. Rev. E 65, 046217 (2002).
  * (14) O. Giraud, Spectral statistics of diractive systems, PhD thesis (2002).
  * (15) O.Giraud, Periodic orbits and semiclassical form factor in barrier billiards, Commun. Math. Phys. 260, 183 (2005).
  * (16) O. Giraud, private communication (2021).
  * (17) G. Tanner, Unitary-stochastic matrix ensembles and spectral statistics, J. Phys. A: Math. Gen. 34, 8485 (2001).
  * (18) S. Gnutzmann and U. Smilansky, Quantum graphs: applications to quantum chaos and universal spectral statistics, Advances in Physics 55, 527 (2006).
  * (19) O. Giraud, J. Marklof, and S. O’Keefe Intermediate statistics in quantum maps, J. Phys. A: Math. Gen. 37, L303 (2004).
  * (20) E. Bogomolny and C. Schmit, Spectral statistics of a quantum interval-exchange map, Phys. Rev. Lett. 93, 254102, (2004).
  * (21) E. Bogomolny, R. Dubertrand, and C. Schmit, Spectral statistics of a pseudo-integrable map: the general case, Nonlinearity 22, 2101, (2009).
  * (22) E. Bogomolny, O. Giraud, and C. Schmit, Random matrix ensembles associated with Lax matrices, Phys. Rev. Lett. 103, 054103, (2009).
  * (23) E. Bogomolny, O. Giraud, and C. Schmit, Integrable random matrix ensembles, Nonlinearity 24, 3179, (2011).
  * (24) L. S. Levitov, Localization-delocalization transition for one-dimensional alloy potentials, EPL 7, 343 (1988).
  * (25) B. L. Altshuler and S. Levitov, Weak chaos in a quantum Kepler problem, Phys. Rep. 288, 487 (1997).
  * (26) U. Grenander and G. Szego, Toeplitz forms and their applications, Univ. of California Press, Berkeley, Los Angeles (1958).
  * (27) A. Böttcher, S. M. Grudsky, and E. A. Maksimenko, Inside the eigenvalues of certain Hermitian Toeplitz band matrices, J. Comp. Appl. Math. 233, 2245 (2010).
  * (28) P. Deift, A. Its, and I. Krasovsky, Eigenvalues of Toeplitz matrices in the bulk of the spectrum, Bull. Inst. Math. Academia Sinica 7, 437 (2012).
  * (29) J. M. Bogoya, A. Böttcher, S. M. Grudsky, and E. A. Maximenko, Eigenvectors of Hermitian Toeplitz matrices with smooth simple-loop symbols, Lin. Algeb. Appl. 493, 606 (2016).
  * (30) P. Rambour, Asymptotic of the eigenvalues of Toeplitz matrices with even symbol, arXiv:2101.11250 (2021).
  * (31) H. Bateman and A. Erdelyi, Higher transcendental functions, v. I. Mc. Graw-Hill Book Company, Inc. (1953).


# Robust matrix completion via Novel M-estimator Functions

Zhi-Yong Wang and Hing Cheung So, E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>This work is supported by a grant from the Research
Grants Council of the Hong Kong Special Administrative Region, China [Project
No. CityU 11207922].

###### Abstract

M-estmators including the Welsch and Cauchy have been widely adopted for
robustness against outliers, but they also down-weigh the uncontaminated data.
To address this issue, we devise a framework to generate a class of nonconvex
functions which only down-weigh outlier-corrupted observations. Our framework
is then applied to the Welsch, Cauchy and $\ell_{p}$-norm functions to produce
the corresponding robust loss functions. Targeting on the application of
robust matrix completion, efficient algorithms based on these functions are
developed and their convergence is analyzed. Finally, extensive numerical
results demonstrate that the proposed methods are superior to the competitors
in terms of recovery accuracy and runtime.

###### Index Terms:

Low-rank matrix completion, matrix factorization, outlier-robustness, implicit
regularizer.

##  I Introduction

Matrix completion (MC) [1, 2] refers to recovering the missing entries of a
partially-observed matrix. It has numerous applications in signal processing
and machine learning, such as hyperspectral imaging [3] and image inpainting
[4]. MC can be formulated as a constrained rank minimization problem [5], but
it is NP-hard since the rank is discrete. To address this, nuclear norm
minimization is exploited [6] to recast MC as a semi-definite program [7].
Although it can be solved by the interior-point method [8], its computational
complexity is high. On the other hand, computationally efficient algorithms
such as singular value thresholding [9], and accelerated proximal gradient
with linesearch [10], have been proposed. Nevertheless, they still involve
full singular value decomposition (SVD) per iteration, implying an expensive
cost especially for large-size data.

To avoid performing SVD, factorization based MC approach [11, 12] has been
developed, whose idea is to utilize the product of two much smaller matrices
to approximately match the observed matrix and thus the low-rank property is
automatically fulfilled. MC algorithms such as low-rank matrix fitting [13]
and alternating minimization [14], are then proposed. In addition, Zhu $et$
$al$. [11, 12] have shown that this kind of MC problem has no spurious local
minima and obeys the strict saddle property which requires the cost function
to have a directional negative curvature at all critical points but local
minima. That is, although the factorization based MC problem is nonconvex,
global optimality can be achieved under some conditions. Nevertheless, the
above mentioned methods are vulnerable to outliers.

To resist outliers, $\ell_{1}$-norm is suggested, resulting in numerous robust
MC algorithms, including robust matrix factorization by majorization
minimization (RMF-MM) [15] and practical low-rank matrix approximation under
robust $\ell_{1}$-norm (RegL1) [16]. Nevertheless, the $\ell_{1}$-norm is
still vulnerable to outliers with large magnitudes because it is not upper
bounded [17]. To enhance outlier-robustness, the $\ell_{p}$-norm with $0<p<1$
[17, 18] and the nonconvex penalty functions such as the Welsch and Tukey, are
adopted [19, 20]. Although these functions penalize outliers via assigning
small weights to outlier-corrupted data, they also down-weigh the normal data.
Here, normal data refer to observations without noise or with only Gaussian
noise. Compared with these nonconvex functions, the Huber function only
penalizes outlier-contaminated entries but is still sensitive to large
outliers because it employs the $\ell_{1}$-norm for robustness [21].

In this paper, we devise a framework to produce a class of M-estimator
functions, which only down-weigh outlier-contaminated data and can resist even
large outliers. The framework is then applied to commonly-used Welsch, Cauchy
and $\ell_{p}$-norm M-estimators, resulting in the corresponding robust
functions. Since these functions are nonconvex, the Legendre-Fenchel (LF)
transform is exploited to convert the nonconvex problem into a sum of convex
problems with closed-form solutions. Furthermore, we apply the developed
functions to factorization based MC and propose robust MC algorithms with
convergence guarantees.

##  II Problem formulation

Let ${\Omega}\subset\\{1,\cdots,m\\}\times\\{1,\cdots,n\\}$ represent the
index set of the known entries of an incomplete matrix
$\boldsymbol{X}_{\Omega}$, and $(\cdot)_{\Omega}$ is a projection operator,
defined as:

$\left[\boldsymbol{X}_{\Omega}\right]_{ij}=\left\\{\begin{aligned}
&X_{ij},\quad{\rm if}~{}(i,j)\in\Omega\\\ &0,\quad\quad\rm{otherwise}.\\\
\end{aligned}\right.$

Given $\boldsymbol{X}_{\Omega}$, MC is to seek a low-rank matrix
$\boldsymbol{M}$ to match $\boldsymbol{X}_{\Omega}$ and estimate its missing
entries, which can be modeled as a rank minimization problem [5]:

$\begin{split}&\mathop{\min}\limits_{\boldsymbol{M}}~{}\text{rank}(\boldsymbol{M}),~{}\text{s.t.}~{}\boldsymbol{M}_{\Omega}=\boldsymbol{X}_{\Omega}\end{split}$
(1)

However, (1) is an NP-hard problem. Instead, the convex nuclear norm is
suggested [6]:

$\begin{split}&\mathop{\min}\limits_{\boldsymbol{M}}~{}\|{\boldsymbol{M}}\|_{*},~{}\text{s.t.}~{}\boldsymbol{M}_{\Omega}=\boldsymbol{X}_{\Omega}\end{split}$
(2)

where the nuclear norm $\|\boldsymbol{M}\|_{*}$ is the sum of singular values
of $\boldsymbol{M}$. Although it is convex, full SVD calculation is required
per iteration. To address this problem, factorization based MC has been
exploited [13]:

$\begin{split}&\mathop{\min}\limits_{\boldsymbol{U,V}}~{}\left\|{\boldsymbol{X}_{\Omega}}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}\right\|_{F}^{2}\end{split}$
(3)

where $\boldsymbol{U}\in\mathbb{R}^{m\times r}$ and
$\boldsymbol{V}\in\mathbb{R}^{r\times n}$ are the two small size matrices with
rank $r\ll{\rm min}(m,n)$. Although the restored matrix
$\boldsymbol{M}=\boldsymbol{U}\boldsymbol{V}$ is low-rank, the recovery
performance is sensitive to outliers because the $\ell_{2}$-norm is utilized.

To resist outliers, loss functions such as the Huber and Welsch functions are
suggested [19, 20], resulting in:

$\begin{split}&\mathop{\min}\limits_{\boldsymbol{U,V}}~{}l({\boldsymbol{X}_{\Omega}}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega})\end{split}$
(4)

where $l(\cdot)$ refers to a robust loss function. Table I tabulates the
commonly-used loss functions and the corresponding weight functions. It is
known that a good robust loss function should only down-weigh the outlier-
corrupted elements. We see that the quadratic function, namely, the
$\ell_{2}$-norm, has the same weights for all entries, including the outlier-
contaminated data, thus it is suitable to handle normal data but is not robust
against outliers. Although the Huber function only assigns small weights for
noisy observations, it is still vulnerable to large outliers since it employs
the $\ell_{1}$-norm to combat gross errors. To enhance outlier-robustness, the
nonconvex functions such as the Welsch and Cauchy, are suggested [22].
However, as shown in Table I, these functions also reduce the weights of
normal data. Recently, we design a novel robust function called hybrid
ordinary-Welsch (HOW) [21] where ‘ordinary’ refers to the quadratic function,
whose expression is shown in Table I. We see that only the Huber and HOW
functions down-weigh outlier-corrupted entries because they assign the same
weights for $|x|\leq c$, where $c$ is a parameter to differentiate whether an
entry is contaminated by outliers or not. That is, when $|x|>c$, the
corresponding element is considered corrupted by an outlier, and assumed
normal entry otherwise.

TABLE I: Commonly-used loss functions and their weight functions | Quadratic | Huber | Cauchy | Welsch | HOW
---|---|---|---|---|---
$l(x)$ | $\frac{x^{2}}{2}$ | $\begin{cases}x^{2}/2,&|x|\leq c\\\ c|x|-\frac{c^{2}}{2},&|x|\textgreater c\end{cases}$ | $\frac{\gamma^{2}}{2}\ln\left(1+\left(\frac{x}{\gamma}\right)^{2}\right)$ | $\frac{\sigma^{2}}{2}\left(1-e^{\frac{-x^{2}}{\sigma^{2}}}\right)$ | $\begin{cases}x^{2}/2,&|x|\leq c\\\ \frac{\sigma^{2}}{2}\left(1-e^{\frac{c^{2}-x^{2}}{\sigma^{2}}}\right)+\frac{c^{2}}{2},&|x|\textgreater c\end{cases}$
$w(x)=\frac{l^{\prime}(x)}{x}$ | $1$ | $\begin{cases}1,&|x|\leq c\\\ c\cdot{\rm sign}(x)/x,&|x|\textgreater c\end{cases}$ | $\frac{1}{1+x^{2}/\gamma^{2}}$ | $e^{\frac{-x^{2}}{\sigma^{2}}}$ | $\begin{cases}1,&|x|\leq c\\\ e^{\frac{c^{2}-x^{2}}{\sigma^{2}}},&|x|\textgreater c\end{cases}$

##  III Novel M-estimator function and its application to robust matrix
completion

In this section, motivated by the thoughts of the construction for Huber and
HOW, we devise a framework to generate a class of M-estimator functions, which
only down-weights the outlier-corrupted data.

###  III-A Framework to Generate M-estimator Functions

We generalize the expressions of Huber and HOW functions for $|x|>c$ and
develop a new generic function:

$l_{g,c}(x)=\begin{cases}x^{2}/2,&|x|\leq c\\\ a\cdot
g(|x|)+b,&|x|\textgreater c\end{cases}$ (5)

where $g(x)$ is a continuous function and $g^{\prime}(x)\geq 0$ for $x>0$,
while $a$ and $b$ are constants to ensure that $l_{g,c}(x)$ is continuous and
smooth at $x=c$. Thus, $a=c/g^{\prime}(c)>0$ ($g^{\prime}(c)\neq 0$), and
$b=c^{2}/2-ag(c)$. We easily see that $l_{g,c}$ only down-weighs the outlier-
corrupted data because when $|x|\leq c$, $l_{g,c}$ is the quadratic function
and does not reduce the weight of normal data, while when $|x|\textgreater c$,
$l_{g,c}$ employs the robust function $g$ to handle outliers. Compared with
the Huber function where $g(x)=|x|_{1}$, we will focus on nonconvex $g(x)$
because the latter can resist large outliers. Hence, the function $l_{g,c}$ is
assumed nonconvex in our study.

Analogous to [21], the LF transform [23] is applied to $l_{g,c}$, resulting
in:

$\begin{split}l_{g,c}(x)=\mathop{\inf}\limits_{y}~{}\frac{(y-x)^{2}}{2}+\varphi_{g,c}(y)\end{split}$
(6)

where $\varphi_{g,c}$ is the dual function of $l_{g,c}$ and is also called the
implicit regularizer (IR). The value of $y$ that solves (6) is:

$P_{\varphi_{g,c}}(x):=\nabla f(x)={\rm max}\left\\{0,|x|-a\cdot
g^{\prime}(|x|)\right\\}\cdot{\rm sign}(x)$ (7)

The process of obtaining (7) from (5) can be found in [21].

Next, we will specify $g$ as several commonly-used functions, including the
$\ell_{p}$-norm and Cauchy function. Note that we have already specified $g$
as the Welsch M-estimator and have proposed the HOW function in [21].

When $g(x)=|x|^{p}$, we have the smooth hybrid ordinary-$\ell_{p}$ (HOP)
function:

$l_{p,c}(x)=\begin{cases}x^{2}/2,&|x|\leq c\\\
\frac{1}{p}c^{2-p}|x|^{p}+\frac{c^{2}}{2}-\frac{1}{p}c^{2},&|x|\textgreater
c\end{cases}$ (8)

By the LF transform, we obtain

$\begin{split}l_{p,c}(x)=\mathop{\min}\limits_{y}~{}\frac{(x-y)^{2}}{2}+\varphi_{p,c}(y)\end{split}$
(9)

where $\varphi_{p,c}(y)$ is the IR related to $l_{p,c}(x)$, and the solution
to $y$ is:

$P_{{\varphi}_{p,c}}(x)={\rm
max}\left\\{0,|x|-{c}^{2-p}|x|^{p-1}\right\\}\cdot{\rm sign}(x)$ (10)

Note that when $p=1$, (8) becomes the Huber function.

We then replace $g(x)$ by the Cauchy M-estimator shown in Table I, and then
develop the hybrid ordinary-Cauchy (HOC) function:

$l_{\gamma,c}(x)=\begin{cases}x^{2}/2,&|x|\leq c\\\
\frac{\gamma^{2}+c^{2}}{2}\ln\left(1+\left(\frac{x}{\gamma}\right)^{2}\right)+b,&|x|\textgreater
c\end{cases}$

where
$b=\frac{c^{2}}{2}-\frac{\gamma^{2}+c^{2}}{2}\ln\left(1+\left({c}/{\gamma}\right)^{2}\right)$
and $\gamma$ is the scale parameter. Employing the LF transform results in:

$\begin{split}l_{\gamma,c}(x)=\mathop{\min}\limits_{y}~{}\frac{(x-y)^{2}}{2}+\varphi_{\gamma,c}(y)\end{split}$

with the solution to $y$ being:

$p_{\varphi_{\gamma,c}}(x):={\rm
max}\left\\{0,|x|-\frac{\left(\gamma^{2}+c^{2}\right)|x|}{\gamma^{2}+x^{2}}\right\\}\cdot{\rm
sign}(x)$ (11)

###  III-B Algorithms for Robust Matrix Completion

We replace the Frobenius norm in (3) with our M-estimator functions, leading
to:

$\mathop{\min}\limits_{\boldsymbol{U},\boldsymbol{V}}~{}l_{g,c}\left(\boldsymbol{X}_{\Omega}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}\right)$
(12)

where
$l_{g,c}\left(\boldsymbol{X}_{\Omega}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}\right)$
is separable, i.e.,
$l_{g,c}\left(\boldsymbol{X}_{\Omega}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}\right)\\\
=\sum_{i,j\in\Omega}l_{g,c}\left(\boldsymbol{X}_{i,j}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{i,j}\right)$.
According to (6), we have:

$\mathop{\min}\limits_{\boldsymbol{U},\boldsymbol{V},\boldsymbol{S}}~{}\mathcal{L}_{g,c}\left(\boldsymbol{U},\boldsymbol{V},\boldsymbol{S}\right):=\frac{1}{2}\left\|\boldsymbol{X}_{\Omega}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}-\boldsymbol{S}_{\Omega}\right\|_{F}^{2}+\varphi_{g,c}(\boldsymbol{S}_{\Omega})$
(13)

where $\varphi(\boldsymbol{S})=\sum_{i,j}\varphi(\boldsymbol{S}_{i,j})$, and
$\boldsymbol{S}_{\Omega^{c}}=\boldsymbol{0}$. In the $(k+1)$th iteration,
given $\boldsymbol{U}^{k}$ and $\boldsymbol{V}^{k}$, (13) is equal to:

$\mathop{\min}\limits_{\boldsymbol{S}}~{}\frac{1}{2}\left\|\boldsymbol{D}^{k}_{\Omega}-\boldsymbol{S}_{\Omega}\right\|_{F}^{2}+\varphi_{g,c}(\boldsymbol{S}_{\Omega})$
(14)

where
$\boldsymbol{D}^{k}=\boldsymbol{X}-\boldsymbol{U}^{k}\boldsymbol{V}^{k}$. The
solution to (14) via (7) is:

$\boldsymbol{S}^{k+1}_{\Omega}={P}_{\varphi_{g,c}}(\boldsymbol{D}^{k}_{\Omega})$
(15)

Given $\boldsymbol{S}^{k+1}$, (13) amounts to:

$\mathop{\min}\limits_{\boldsymbol{U},\boldsymbol{V}}~{}h(\boldsymbol{U},\boldsymbol{V}):=\frac{1}{2}\left\|\boldsymbol{H}^{k+1}_{\Omega}-\left(\boldsymbol{U}\boldsymbol{V}\right)_{\Omega}\right\|_{F}^{2}$
(16)

where
$\boldsymbol{H}^{k+1}_{\Omega}=\boldsymbol{X}_{\Omega}-\boldsymbol{S}^{k+1}_{\Omega}$,
and it can be efficiently solved by the scaled alternating steepest descent
(SASD) [24]. Then, the scaled gradient descent directions for our case are:

$\displaystyle\widetilde{\nabla}h_{\boldsymbol{V}}(\boldsymbol{U})=\left(\boldsymbol{H}_{\Omega}-(\boldsymbol{U}\boldsymbol{V})_{\Omega}\right)\boldsymbol{V}^{T}(\boldsymbol{V}\boldsymbol{V}^{T})^{-1}$
(17a)
$\displaystyle\widetilde{\nabla}h_{\boldsymbol{U}}(\boldsymbol{V})=(\boldsymbol{U}^{T}\boldsymbol{U})^{-1}\boldsymbol{U}^{T}\left(\boldsymbol{H}_{\Omega}-(\boldsymbol{U}\boldsymbol{V})_{\Omega}\right)$
(17b)

with the corresponding step sizes being:

$\displaystyle\widetilde{\mu}_{\boldsymbol{U}}=\left<\nabla
h_{\boldsymbol{V}}(\boldsymbol{U}),\widetilde{\nabla}h_{\boldsymbol{V}}(\boldsymbol{U})\right>\bigg{/}\left\|\left(\widetilde{\nabla}h_{\boldsymbol{V}}(\boldsymbol{U})\boldsymbol{V}\right)_{\Omega}\right\|_{F}^{2}$
(18a) $\displaystyle\widetilde{\mu}_{\boldsymbol{V}}=\left<\nabla
h_{\boldsymbol{U}}(\boldsymbol{V}),\widetilde{\nabla}h_{\boldsymbol{U}}(\boldsymbol{V})\right>\bigg{/}\left\|\left(\boldsymbol{U}\widetilde{\nabla}h_{\boldsymbol{U}}(\boldsymbol{V})\right)_{\Omega}\right\|_{F}^{2}$
(18b)

Thus, the SASD updates for $\boldsymbol{U}$ and $\boldsymbol{V}$ are:

$\displaystyle\boldsymbol{U}^{k+1}=\boldsymbol{U}^{k}-\widetilde{\mu}_{\boldsymbol{U}}^{k}\widetilde{\nabla}h_{\boldsymbol{V}^{k}}(\boldsymbol{U}^{k})$
(19a)
$\displaystyle\boldsymbol{V}^{k+1}=\boldsymbol{V}^{k}-\widetilde{\mu}_{\boldsymbol{V}}^{k}\widetilde{\nabla}h_{\boldsymbol{U}^{k+1}}(\boldsymbol{V}^{k})$
(19b)

Recall that the value of $c$ in (5) is the boundary to determine whether an
entry is corrupted or not. Similar to [19, 21], its value is set as:

$\begin{split}c^{k}&=\min\left\\{\xi d^{k},c^{k-1}\right\\}\end{split}$ (20)

where $\xi\textgreater 0$ is a user-defined constant, $d^{k}$ is the robust
normalized interquartile range of the vectorized $\boldsymbol{D}_{\Omega}$,
defined as:

$\begin{split}d^{k}={\rm IQR}\left({\rm
vec}(\boldsymbol{D}_{\Omega}^{k})\right)/1.349\end{split}$ (21)

with ${\rm IQR}(\cdot)$ being the sample interquartile range operator [25].

When HOW, HOP and HOC are adopted in (12), the resultant algorithms are
referred to as robust MC via HOW (RMC-HOW), HOP (RMC-HOP) and HOC (RMC-HOC),
respectively.

In addition, SASD dominates the complexity of the proposed approach with
complexity of $\mathcal{O}\left(8|\Omega|r+4(m+n)r^{2}\right)$ per iteration.
Defining
$E_{k}=l_{g,c}\left(\boldsymbol{X}_{\Omega}-\left(\boldsymbol{U}^{k}\boldsymbol{V}^{k}\right)_{\Omega}\right)$,
we terminate our algorithms until the relative error
$rel_{E}^{k}=\left(E_{k}-E_{k-1}\right)/E_{k-1}<\zeta$. On the other hand, the
convergence analysis results are shown in the following theorems, whose proofs
are analogous to those in our previous work [21], and we omit them due to page
limit.

######  Theorem 1.

The generated sequence
$\left\\{\mathcal{L}_{g^{k},c^{k}}\left(\boldsymbol{U}^{k},\boldsymbol{V}^{k},\boldsymbol{S}^{k}\right)\right\\}$
converges.

######  Theorem 2.

Let
$\left\\{{\left(\boldsymbol{U}^{k},\boldsymbol{V}^{k},\boldsymbol{S}^{k}\right)}\right\\}$
be the generated sequence, and suppose that
$\left(\boldsymbol{U}^{k},\boldsymbol{V}^{k}\right)$ are of full rank. Then,
$\left\\{{\left(\boldsymbol{U}^{k},\boldsymbol{V}^{k},\boldsymbol{S}^{k}\right)}\right\\}$
is bounded. Besides, let
$\left\\{{\left(\boldsymbol{U}^{k_{j}},~{}\boldsymbol{V}^{k_{j}},~{}\boldsymbol{S}^{k_{j}}\right)}\right\\}$
be a generated subsequence such that
$\lim_{k_{j}\rightarrow\infty}\left(\boldsymbol{U}^{k_{j}},\boldsymbol{V}^{k_{j}},\boldsymbol{S}^{k_{j}}\right)\\\
=\left(\boldsymbol{U}^{\star},\boldsymbol{V}^{\star},\boldsymbol{S}^{\star}\right)$.
Then,
${\left(\boldsymbol{U}^{\star},\boldsymbol{V}^{\star},\boldsymbol{S}^{\star}\right)}$
is a critical point.

##  IV Experimental Results

We compare our algorithms with the competing methods, including HQ-ASD [19],
$(\boldsymbol{S}+\boldsymbol{L})_{{1}/{2}}$ [18],
$(\boldsymbol{S}+\boldsymbol{L})_{{2}/{3}}$ [18], RMF-MM [15] and $\rm
RegL_{1}$ [16]. All numerical simulations are conducted using a computer with
3.0 GHz CPU and 16 GB memory. We first generate a low-rank matrix
$\boldsymbol{X}=\boldsymbol{U}\boldsymbol{V}$, where the entries of
$\boldsymbol{U}\in\mathbb{R}^{m\times r}$ and
$\boldsymbol{V}\in\mathbb{R}^{r\times n}$ satisfy the standard Gaussian
distribution. Impulsive noise is modeled by the Gaussian mixture model (GMM)
while the signal-to-noise ratio (SNR) is defined as:

${\rm
SNR}=\frac{\|\boldsymbol{X}_{\Omega}\|_{F}^{2}}{|\Omega|\left((1-\tau)\sigma_{1}^{2}+\tau\sigma_{2}^{2}\right)}$
(22)

where $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$ are variances with
$\sigma_{1}^{2}\ll\sigma_{2}^{2}$, and $\tau$ controls the proportion of
outliers. To model outliers, we set $\sigma_{2}^{2}=100\sigma_{1}^{2}$ and
$\tau=0.1$. Besides, the root mean square error (RMSE) defined as ${\rm
RMSE}={\|\boldsymbol{X}-\boldsymbol{M}\|_{F}}/{\sqrt{mn}}$ is utilized to
measure the performance of all algorithms. Similar to the parameter selection
in [21], we set $\xi=2$ and $\zeta=10^{-4}$ in our methods.

Figure 1: RMSE versus SNR. Figure 2: RMSE versus percentage of observations. TABLE II: Average runtime (in seconds) for different matrix sizes.  Method | Case 1 | Case 2 | Case 3 | Case 4
---|---|---|---|---
RMF-MM | $22.98$ | $175.1$ | $460.7$ | $855.2$
RegL1 | $0.161$ | $0.963$ | $2.233$ | $5.346$
$(S+L)_{1/2}$ | $0.663$ | $2.548$ | $6.349$ | $13.30$
$(S+L)_{2/3}$ | $0.573$ | $2.340$ | $5.852$ | $11.59$
HQ-ASD | $0.080$ | $0.372$ | $0.846$ | $1.765$
RMC-HOW | $0.027$ | $0.118$ | $0.268$ | $0.532$
RMC-HOC | $0.029$ | $0.136$ | $0.300$ | $0.608$
RMC-HOP ($p=0.6$) | $0.048$ | $0.199$ | $0.425$ | $0.812$
RMC-HOP ($p=0.3$) | $0.050$ | $0.205$ | $0.454$ | $0.864$

We conduct experiments on data matrices with $m=300$, $n=200$ and $r=5$. Fig.
1 plots the RMSE versus SNR with $30\%$ observations. It is seen that the RMSE
for all methods decreases with SNR, while the proposed algorithms have smaller
recovery error than the competing techniques. In addition, RMC-HOW yields the
best recovery because HOW is bounded. Moreover, the impact of percentage of
observations is investigated, and the results are shown in Fig. 2. Again, our
methods outperform the competitors, with HOW attaining the best recovery
performance.

Finally, we investigate the runtime of all methods, and perform numerical
experiments on four cases with different matrix dimensions. Here, Case 1:
$m=300,n=200,r=5$, Case 2: $m=600,n=400,r=10$, Case 3: $m=900,n=600,r=15$, and
Case 4: $m=1200,n=800,r=20$. The results under SNR = $10$ dB and $50\%$
observations are tabulated in Table II. We see that the runtime of our
algorithms is less than that of the competitors.

##  V Conclusion

In this paper, we provide a framework to generate a new class of robust loss
functions via combining the quadratic and other robust loss functions. The
proposed functions can be used to combat gross errors and only down-weigh the
outlier-contaminated observations. Applying our framework to the Welsch,
Cauchy and $\ell_{p}$-norm functions yields the HOW, HOC and HOP,
respectively, which are then adopted for robust MC. Furthermore, although the
resultant optimization problem is nonconvex, the LF transform is adopted to
transform it into a sum of convex subproblems. Then, efficient robust MC
algorithms based on SASD are developed. Finally, experimental results show the
superiority of the proposed algorithms over the competing methods in terms of
recovery error and runtime.

## References

  * [1] M. A. Davenport and J. Romberg, “An overview of low-rank matrix recovery from incomplete observations,” _IEEE J. Sel. Topics Signal Process._ , vol. 10, no. 4, pp. 608–622, Jun. 2016.
  * [2] Y. Chen and Y. Chi, “Harnessing structures in big data via guaranteed low-rank matrix estimation: Recent theory and fast algorithms via convex and nonconvex optimization,” _IEEE Signal Process. Mag._ , vol. 35, no. 4, pp. 14–31, Jul. 2018.
  * [3] Z.-Y. Wang, X. P. Li and H. C. So, “Robust matrix completion based on factorization and truncated-quadratic loss function,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 33, no. 9, pp. 1521-1534, Apr. 2023.
  * [4] Z.-Y. Wang, X. P. Li, H. C. So and A. M. Zoubir, “Adaptive rank-one matrix completion using sum of outer products.” _IEEE Trans. Circuits Syst. Video Technol._ , Mar. 2023. Early Access.
  * [5] E. J. Candès, and B. Recht, “Exact matrix completion via convex optimization,” _Found. Comut. Math._ , vol. 9, no. 6, pp. 717, Dec. 2009.
  * [6] M. Fazel, “Matrix rank minimization with applications,” Ph.D. dissertation, Dept. Elect. Eng. Stanford Univ. California, USA, 2002.
  * [7] E. J. Candès and Y. Plan, “Matrix completion with noise,” _Proc. IEEE_ , vol. 98, no. 6, pp. 925–936, Jun. 2010.
  * [8] Z. Liu, and L. Vandenberghe, “Interior-point method for nuclear norm approximation with application to system identification,” _SIAM J. Matrix Anal. Appl._ , vol. 31, no. 3, pp. 1235–1256, Aug. 2009.
  * [9] J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” _SIAM J. Opt._ , vol. 20, no. 4, pp. 1956–1982, Mar. 2010.
  * [10] K. C. Toh and S. W. Yun, “An accelerated proximal gradient algorithm for nuclear norm regularized least squares problems,” _Pacific. J. Optim._ , vol. 6, no. 15, pp. 615–640, Sep. 2010.
  * [11] Z. Zhu, Q. Li, G. Tang, and M. B. Wakin, “Global optimality in low-rank matrix optimization,” _IEEE Trans. Signal Process._ , vol. 66, no. 13, pp. 3614–3628, Jul. 2018.
  * [12] Z. Zhu, Q. Li, G. Tang, and M. B. Wakin, “The global optimization geometry of low-rank matrix optimization,” _IEEE Trans. Inf. Theory_ , vol. 67, no. 2, pp. 1308–1331, Feb. 2021.
  * [13] Z. Wen, W. Yin, and Y. Zhang, “Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm,” _Math. Programm. Comput._ , vol. 4, no. 4, pp. 333–361, Dec. 2012.
  * [14] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in _Proc. 45th Annu. ACM Symp. Symp. Theory Comput. (STOC)_ , 2013, pp. 665–674.
  * [15] Z. Lin, C. Xu, and H. Zha, “Robust matrix factorization by majorization minimization,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 40, no. 1, pp. 208–220, Jan. 2018.
  * [16] Y. Zheng, G. Liu, S. Sugimoto, S. Yan, and M. Okutomi, “Practical low-rank matrix approximation under robust $\ell_{1}$-norm,” in _Proc. IEEE Conf. Comput. Vis. Pattern Recognit._ , Providence, USA, Jun. 2012, pp. 1410–1417.
  * [17] W.-J. Zeng and H. C. So, “Outlier-robust matrix completion via $\ell_{p}$-minimization,” _IEEE Trans. Signal Process._ , vol. 66, no. 5, pp. 1125–1140, Mar. 2018.
  * [18] F. Shang, J. Cheng, Y. Liu, Z. Luo, and Z. Lin, “Bilinear factor matrix norm minimization for robust PCA: Algorithms and applications,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 40, no. 9, pp. 2066–2080, Sep. 2018\.
  * [19] Y. He, F. Wang, Y. Li, J. Qin, and B. Chen, “Robust matrix completion via maximum correntropy criterion and half-quadratic optimization,” _IEEE Trans. Signal Process._ , vol. 68, pp. 181–195, Nov. 2019.
  * [20] M. Muma, W.-J. Zeng, and A. M. Zoubir, “Robust M-estimation based matrix completion,” in _Proc. IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP)_ , Brighton, UK, May 2019, pp. 5476–5480.
  * [21] Z.-Y. Wang, H. C. So and A. M. Zoubir, “Robust low-rank matrix recovery via hybrid ordinary-Welsch function.” _IEEE Trans. Signal Process._ , vol. 71, pp. 2548-2563, Jul. 2023.
  * [22] F. D. Mandanas and C. L. Kotropoulos, “Robust multidimensional scaling using a maximum correntropy criterion,” _IEEE Trans. Signal Process._ , vol. 65, no. 4, pp. 919–932, Feb. 2017.
  * [23] R. T. Rockafellar and R. J. B. Wets, _Variational Analysis_ , 2nd ed. Berlin, Germany: Springer, 2004.
  * [24] J. Tanner and K.Wei, “Low rank matrix completion by alternating steepest descent methods,” _Appl. Comput. Harmon. Anal._ , vol. 40, no.2, pp. 417–429, Mar. 2015.
  * [25] R. A. Maronna, D. R. Martin and V. J. Yohai, _Robust Statistics: Theory and Methods_ , England: Wiley, 2006.