# Notes on the complexity of coverings for Kronecker powers of symmetric
matrices

Igor S. Sergeev111e-mail<EMAIL_ADDRESS>
###### Abstract

In the present note, we study a new method of constructing efficient coverings
for Kronecker powers of matrices, recently proposed by J. Alman, Y. Guan, A.
Padaki [1]. We provide an alternative proof for the case of symmetric matrices
in a stronger form. As a consequence, the previously known upper bound on the
depth-2 additive complexity of the boolean $N\times N$ Kneser-Sierpinski
matrices is improved to $O(N^{1.251})$. This work can be viewed as a
supplement to [3].

##  1 Introduction

Let us recall necessary concepts. See [3] for a more detailed introduction to
the subject.

A rectangle of size $a\times b$ is an all-1s matrix with $a$ rows and $b$
columns. Further, depending on the context, sometimes under rectangle we will
understand a rank-1 matrix, i.e. consisting of an all-1s submatrix and all 0s
in other entries.

We define the complexity222In [3], we used a term weight instead. Here we
substitute it with complexity to avoid confusing with the spectral weight of a
rectangle. of an $a\times b$ rectangle $R$ as the sum of lengths of its two
sides, $w(R)=a+b$. We introduce the characteristic of the narrowness of a
rectangle as the ratio of the lengths of its larger and smaller sides,
$\rho(R)=\frac{\max(a,\,b)}{\min(a,\,b)}$. The spectral weight of a rectangle
is defined as $\sigma(R)=\sqrt{ab}$.

A set $F=\\{R_{1},\ldots,R_{k}\\}$ of rectangles is a covering of a boolean
matrix $A$, if

$A=R_{1}+\ldots+R_{k}.$ (1)

(Here under $R_{i}$ we mean rank-1 matrices.) If the operation “+” in (1) is
an integer addition, then $F$ is called $\mathsf{SUM}$-covering. If “+” is a
disjunction, then $F$ is called $\mathsf{OR}$-covering. In the case when “+”
is a $\bmod\,2$ addition, then we have an $\mathsf{XOR}$-covering.

The complexity of a covering $F$ is defined as
$w(F)=w(R_{1})+\ldots+w(R_{k})$, and the spectral weight as
$\sigma(F)=\sigma(R_{1})+\ldots+\sigma(R_{k})$. The $\mathsf{L}$-complexity of
a matrix $A$ is defined for
$\mathsf{L}\in\\{\mathsf{SUM},\mathsf{OR},\mathsf{XOR}\\}$ as the minimal
complexity of its $\mathsf{L}$-covering, denoted by $\mathsf{L}_{2}(A)$ (it
means the complexity of computation of $A$ by depth-2 linear circuits of the
corresponding type).

$\blacktriangleright$ These notions may be extended to the case of matrices
over an arbitrary semiring $S$. A rectangle over $S$ is a matrix
$(c_{1},\ldots,c_{a})^{T}\cdot(d_{1},\ldots,d_{b})$, where $c_{i},d_{j}\in
S\setminus\\{0\\}$. A covering of a matrix is conditioned by

$A=e_{1}R_{1}+\ldots+e_{k}R_{k},\qquad e_{i}\in S.$

The results presented below may be applied also to the analogously defined
measure of complexity of computation of matrices by algebraic linear circuits
of depth 2. $\blacktriangleleft$

Let $\sigma(A)$ denote the minimal spectral weight of a matrix $A$. Since
$w(R)\geq 2\sigma(R)$ for any rectangle $R$, spectral weight serves as a
simple lower bound for complexity333Here and below, symbols $\asymp$, $\prec$,
$\preceq$ denote the equality, strict and non-strict inequalities on the order
of growth.: $\mathsf{L}_{2}(A)\succeq\sigma(A)$.

A convenient property of the spectral weight is its multiplicativity with
respect to the Kronecker product. Recall that the Kronecker product of boolean
matrices $A$, $B$ is a matrix $A\otimes B$ obtained by replacing 1-entries of
$A$ by copies of $B$, and 0-entries by all-0s matrices of the same size.

Note that if $F$ and $G$ are coverings of matrices $A$ and $B$, then $F\otimes
G$ is a covering444The Kronecker product of sets of matrices is $F\otimes
G=\\{R\otimes R^{\prime}\mid R\in F,\,R^{\prime}\in G\\}$. of $A\otimes B$,
and $\sigma(F\otimes G)=\sigma(F)\sigma(G)$. In particular, if we construct a
covering of a matrix $A^{\otimes n}$ (a Kronecker power of $A$) by the
product-of-coverings method above using appropriate coverings of $A$, then the
complexity of a resulting covering $H$ satisfies
$w(H)\succeq\sigma(H)\succeq\sigma^{n}(A)$.

In [1], the authors actually pose a question: can we obtain upper bounds like
$\mathsf{L}_{2}(A^{\otimes n})\preceq\sigma^{n+o(n)}(A)$ or at least
$\mathsf{L}_{2}(A^{\otimes n})\preceq\sigma^{n+o(n)}(F)$ for some appropriate
coverings $F$ of a matrix $A$. In general, it is not possible, just consider
an example $A=\begin{bmatrix}1&1\\\ 0&0\end{bmatrix}$.

A general obstacle for the desired bounds is the growing narrowness of the
covering rectangles. Note that the complexity and the spectral weight of a
rectangle $R$ are related as555To be precise,
$w(R)=\left(\sqrt{\rho(R)}+\sqrt{1/\rho(R)}\right)\sigma(R)$.
$w(R)\asymp\sqrt{\rho(R)}\sigma(R)$. Nevertheless, under some conditions, the
authors of [1] were able to overcome the indicated obstacle and to derive the
desired bounds. The first situation is when the covering $F$ satisfies some
asymmetry criteria, the second one is when the matrix $A$ is symmetric, and
the covering $F$ is one-sided (it means that all rectangles are stretched in
the same direction).

Perhaps, the most challenging object (among boolean matrices) to apply the
theory are the Kneser–Sierpinski (or disjointness) matrices. Recall that the
boolean $N\times N$ Kneser–Sierpinski matrix $D_{N}$ is defined for $N=2^{n}$
as follows. Rows and columns of $D_{N}$ are labeled by distinct subsets
$u\in[n]$, and $D[u,v]=1$ iff $u\cap v=\emptyset$. A matrix $D_{N}$ also may
be viewed as a Kronecker product

$D_{N}=D_{2}^{\otimes n}=\underbrace{D_{2}\otimes\ldots\otimes
D_{2}}_{n},\qquad D_{2}=\begin{bmatrix}1&1\\\ 1&0\end{bmatrix}.$ (2)

The problem of complexity of $\mathsf{OR}$-coverings for $D_{N}$ was almost
closed in [2]. There was established that

$N^{1.16}\prec\mathsf{OR}_{2}(D_{N})\prec N^{1.17}$

(the lower bound is from [3]). Moreover, the authors of [2] constructed a
covering of almost minimal complexity, up to a factor of order $(\log
N)^{O(1)}$.

The question about additive ($\mathsf{SUM}$) complexity of matrices is less
clear. In [3], we propose a simple way to show that

$\mathsf{SUM}_{2}(D_{N})\preceq\sigma^{n}(D_{2})=\left(\sqrt{2}+1\right)^{n}\prec
N^{1.272}$

relying on a trivial (and optimal) decomposition of the matrix $D_{2}$ into
rectangles of size $1\times 1$ and $1\times 2$ (or $2\times 1$). This approach
was nontrivially generalized in [1]. Due to limitations inherent in the
analysis of the proposed method, the efficient implementation of the
Kneser–Sierpinski matrices is justified only with the basic coverings of
matrices $D_{4}$ and $D_{8}$. In the latter case, the obtained bound [1] is

$\mathsf{SUM}_{2}(D_{N})\preceq\sigma^{n/3+o(n)}(D_{8})=\left(\sqrt{8}+\sqrt{7}+3\sqrt{3}+3\right)^{n/3+o(n)}\prec
N^{1.258}.$

In the present note, we describe a version of this method for the (most
interesting) case of symmetric matrices. The limits of applicability of the
method are (comparatively) extended, and a more elementary proof is given. As
a consequence, an upper bound for the complexity of the Kneser–Sierpinski
matrices is reduced to $\mathsf{SUM}_{2}(D_{N})\prec N^{1.251}$.

Note that the question about existence of substantially more efficient
$\mathsf{XOR}$-coverings for matrices $D_{N}$ is still open.

##  2 The synthesis method

In this section, we provide a general method for constructing a covering of a
symmetric matrix $A^{\otimes n}$. The method is essentially equivalent to the
method [1]. By the way, we will follow an illustrating example where a non-
trivial covering of $D_{N}$ is obtained from coverings of the matrix $D_{4}$.
This example doesn’t require the method in its full generality.

Assume there exist two coverings of a symmetric matrix $A$, or to be precise,
two pairs of coverings, if we count transposed ones. The first covering $F$
($F^{T}$) is supposed to be efficient in the terms of spectral weight. The
second covering $G$ is one-sided: the longest sides of all its non-square
rectangles are parallel. Coverings $G/G^{T}$ serve to compensate an imbalance
caused by the use of $F$-type coverings.

$\blacktriangleright$ Fig. 1 shows appropriate coverings of $D_{4}$. The
covering $F_{2}$ (on the left) has optimal spectral weight
$\sigma(F_{2})=4+\sqrt{3}$. The covering $G_{2}$ (on the right) has good
correcting qualities. Its spectral weight is slightly higher,
$\sigma(G_{2})=3+2\sqrt{2}$. $\blacktriangleleft$

$1$$1$$1$$1$$1$$1$$1$$1$$1$$0$$0$$0$$0$$0$$0$$0$$F_{2}$$1$$1$$1$$1$$1$$1$$1$$1$$1$$0$$0$$0$$0$$0$$0$$0$$G_{2}$
Figure 1: Weight-minimal and compensating coverings of $D_{4}$

Let $H$ be a covering of a matrix $B$. Then we can build a covering of a
matrix $A\otimes B$ in the form $\\{F_{i}\otimes R_{i}\mid R_{i}\in H\\}$,
where $F_{i}$ are some coverings of $A$. This way we sequentially obtain
coverings for matrices $A,A^{\otimes 2},A^{\otimes 3},\ldots$ In the main
process, we choose $F_{i}\in\\{F,\,F^{T}\\}$. Precisely, we transform an
$a\times b$ rectangle $R$ into $F\otimes R$, if $a\leq b$, and into
$F^{T}\otimes R$, otherwise.

Consider a covering $F=\\{R_{1},\ldots,R_{s}\\}$ consisting of $a_{i}\times
b_{i}$ rectangles $R_{i}$. The characteristic function of $F$ is defined as

$\chi_{F}(x)=\sigma(R_{1})\left(\frac{a_{1}}{b_{1}}\right)^{x}+\ldots+\sigma(R_{s})\left(\frac{a_{s}}{b_{s}}\right)^{x}-\sigma(F).$
(3)

As follows from the definition, $\chi_{F}(0)=0$. We call a covering $F$
compact, if $\chi_{F}(x)$ takes negative values on some negative
arguments666It is a weak analogue of imbalanced covering from [1].. For a
compact covering $F$, let $\lambda_{F}$ denote the minimal real root of
$\chi_{F}$ in whose right semineighbourhood the function is negative777The
compactness of a covering implies that $b_{i}<a_{i}$ for some $i$. Then,
$\lambda_{F}$ is correctly defined since $\chi_{F}(x)\to+\infty$ as
$x\to-\infty$, and any exponential sum of the form (3) has a finite number of
real zeros, see e.g. [4]..

Note that the condition
$\chi^{\prime}_{F}(0)=\sum_{i}\sigma(R_{i})\ln(a_{i}/b_{i})>0$ is sufficient
for $F$ to be compact.

So, we require the compactness of a covering $F$ for our algorithm of
computation of a Kronecker power of a matrix $A$.

$\blacktriangleright$ The covering $F_{2}$ of $D_{4}$ is compact. Its
characteristic function is $\chi_{F_{2}}(x)=2\cdot 4^{x}+\sqrt{3}\cdot
3^{-x}-2-\sqrt{3}$, and the minimal root is $\lambda_{F_{2}}\approx-0.305$.
$\blacktriangleleft$

For compensation, we will use a compact one-sided covering $G$ consisting of
$a^{\prime}_{i}\times b^{\prime}_{i}$ rectangles $R^{\prime}_{i}$ satisfying
$a^{\prime}_{i}\geq b^{\prime}_{i}$. The quality of such covering is
characterized by the coefficient

$\mu_{G}=\frac{1}{\sigma(G)}\sum_{R\in
G}\frac{\sigma(R)}{\sqrt{\rho(R)}}=\frac{1}{\sigma(G)}\sum_{i}b^{\prime}_{i}.$

For $x>1$, define the function

$\pi_{G}(x)=\frac{1}{\sigma(G)}\sum_{R\in G}\sigma(R)\cdot
x^{-\lfloor\log_{x}\rho(R)\rfloor/2}.$ (4)

It easily follows from the definition that $\pi_{G}(x)\geq\mu_{G}$, and
$\pi_{G}(x)\to\mu_{G}$ as $x\to 1$.

$\blacktriangleright$ For the covering $G_{2}$, we have
$\mu_{G_{2}}=\frac{4}{3+2\sqrt{2}}$. $\blacktriangleleft$

######  Theorem 1.

Let $F$ be a compact $\mathsf{L}$-covering, and $G$ be a compact one-sided
$\mathsf{L}$-covering of a symmetric $r\times r$ matrix $A$, and
$\sigma(G)\geq\sigma(F)$. If the condition

$\frac{\sigma(G)}{\sigma(F)}<\mu_{G}^{2\lambda_{F}}$ (5)

is satisfied, then for $N=r^{n}$,

$\mathsf{L}_{2}(A^{\otimes n})\preceq N^{\log_{r}\sigma(F)}.$

$\blacktriangleright$ The matrix $D_{4}$ and its coverings $F_{2},G_{2}$
satisfy the conditions of the theorem, since $\lambda_{F_{2}}<-0.3$ and
$\mu_{G_{2}}=\frac{4}{3+2\sqrt{2}}$. Hence, $\mathsf{SUM}_{2}(D_{N})\preceq
N^{\log_{4}(4+\sqrt{3})}\prec N^{1.26}.$ $\blacktriangleleft$

###### Proof.

The proof strategy is the following. First, we analyze the evolution of
rectangle sizes after multiple application of type-$F$ coverings. Then, do the
same for type-$G$ coverings. Finally, we propose an appropriate combination of
these two types of coverings.

I. By compactness of the covering $F$, for some (small enough)
$\delta,\epsilon>0$, and $\lambda=\lambda_{F}+\delta$ such that
$\chi_{F}(\lambda)<-\epsilon$, the inequality (5) holds true after replacing
$\lambda_{F}$ by $\lambda$. Let us check that for all small enough $\tau>1$,

$\sigma(R_{1})\tau^{\lambda\lfloor\log_{\tau}(a_{1}/b_{1})\rfloor}+\ldots+\sigma(R_{s})\tau^{\lambda\lfloor\log_{\tau}(a_{s}/b_{s})\rfloor}\leq\sigma(F).$
(6)

Indeed, the left side of (6) doesn’t exceed

$\tau^{-\lambda}(\chi_{F}(\lambda)+\sigma(F))<\tau^{-\lambda}(\sigma(F)-\epsilon),$

thus to satisfy (6), it is sufficient to require $\tau^{-\lambda}\leq
1+\frac{\epsilon}{\sigma(F)}$. Therefore, any choice from the interval
$1<\tau\leq\left(1+\frac{\epsilon}{\sigma(F)}\right)^{-1/\lambda}$ is
suitable.

We assign to the parameter $\tau$ the meaning of a discretization step of
changing the ratios between the rectangle’s sides in the classification of
rectangles. The final choice of $\tau$ will be decided later.

Let us introduce a classification on the set $\cal R$ of all rectangles
depending on the ratio between the longer and the shorter sides. Set ${\cal
R}=\bigcup_{k\geq 0}I_{k}$, where $I_{0}$ contains rectangles $R$ satisfying
$\rho(R)\leq r$, and for $k\geq 1$, the set $I_{k}$ contains rectangles with
ratios $r\cdot\tau^{k-1}<\rho(R)\leq r\cdot\tau^{k}$. Recall that $r$ is the
size of the matrix $A$.

$\blacktriangleright$ In the example with the covering $F_{2}$, we set
$\tau=4$. Then, $I_{0}=\\{R\mid\rho(R)\leq 4\\}$, and $I_{k}=\\{R\mid
4^{k}<\rho(R)\leq 4^{k+1}\\}$ for $k>0$. The function $\chi_{F_{2}}(x)$ is
negative in the interval $(\lambda_{F_{2}},\,0)$, thus we are quite free in
choosing $\lambda$. The specific value will be determined later.
$\blacktriangleleft$

When performing compositions with coverings of $A$, we will track the
distribution of the spectral weight of rectangles among the sets $I_{k}$. In
doing so, we will be guided by the principle of error to the right. It means
that: (a) we allow a rectangle to be placed into a set $I_{k}$ with a higher
index $k$, but not otherwise, and (b) we estimate the distribution of the
weight of a covering of $A\otimes R$ only on the basis of information about
assigning $R$ to a certain set $I_{k}$. As a consequence, the estimated
distribution of a covering obtained as a result of a series of iterations
(with possible errors to the right) may differ from the true distribution only
in that some rectangles appear in sets $I_{k}$ with higher indices.

The redistribution of the spectral weight of a set of rectangles under the
composition with type-$F$ coverings (with possible errors to the right) may be
estimated from the coefficients of the Laurent polynomial

$P_{F}(x)=\sum_{i\in\mathbb{Z}}\beta_{i}x^{i}=\frac{\sigma(R_{1})}{\sigma(F)}x^{\lfloor\log_{\tau}(a_{1}/b_{1})\rfloor}+\ldots+\frac{\sigma(R_{s})}{\sigma(F)}x^{\lfloor\log_{\tau}(a_{s}/b_{s})\rfloor}$
(7)

obtained from (6). It means that for $R\in I_{m}$, the weight of rectangles
from $(F/F^{T})\otimes R$ is distributed so that the portion $\beta_{k}$ of
the total weight belongs to $I_{m-k}$ in the case $m-k>0$, and to $I_{0}$,
otherwise.

$\blacktriangleright$ For the covering $F_{2}$, and $\tau=4$, we obtain
$P_{F_{2}}(x)=\frac{2x+2+\sqrt{3}x^{-1}}{4+\sqrt{3}}$. The corresponding
redistribution diagram is shown on Fig. 2. $\blacktriangleleft$

$I_{0}$$I_{1}$$I_{2}$$I_{k}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{\sqrt{3}/\omega}$$\scriptstyle{\sqrt{3}/\omega}$$\scriptstyle{\sqrt{3}/\omega}$$\scriptstyle{\sqrt{3}/\omega}$$\scriptstyle{\sqrt{3}/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{2/\omega}$$\scriptstyle{4/\omega}$……
Figure 2: Diagram of spectral weight redistribution under the action of the
composition with the covering $F_{2}/F_{2}^{T}$ (here
$\omega=\sigma(F_{2})=4+\sqrt{3}$)

Let $p_{k}(t)$ stand for the fraction of the spectral weight of the
constructed covering of $A^{\otimes t}$ associated with the set $I_{k}$. In
the beginning, one has $p_{0}(0)=1$, and $p_{k}(0)=0$ for all $k>0$. Set
$\nu=\tau^{\lambda}$. By (6), $P_{F}(\nu)\leq 1$. Denote $d=\deg
P_{F}=\max_{\beta_{i}>0}|i|$.

We are going to show that the distribution $\\{p_{k}^{*}(t)\\}$ with
$p^{*}_{k}(t)=\nu^{k}$ for $k\leq dt$, and $p^{*}_{k}(t)=0$ for all $k>dt$, is
a majorant for $\\{p_{k}(t)\\}$, meaning that the values $p^{*}_{k}(t)$ upper
bound the components of some distribution $\\{p^{\prime}_{k}(t)\\}$ obtained
from $\\{p_{k}(t)\\}$ by a partial shift of the distribution to the right:
from components with smaller indices to components with greater
indices888Though the distribution $\\{p_{k}(t)\\}$ is probabilistic, we don’t
require the same from the majorant $\\{p_{k}^{*}(t)\\}$ allowing the sum of
its components be greater than 1. Our goal is just deriving upper bounds on
$p_{k}(t)$..

Obviously, in the moment $t=0$, the majorization condition is fulfilled. Let
us prove the induction step: apply the composition with $F/F^{T}$ to a set of
rectangles with the distribution $\\{p_{k}^{*}(t)\\}$. Rectangles from
$I_{k}$, where $k>d(t+1)$, do not appear here. For $0<k\leq d(t+1)$, the
weight of rectangles from $I_{k}$ may be upper bounded as

$\sum_{i\in\mathbb{Z}}\beta_{i}\nu^{k+i}=\nu^{k}P_{F}(\nu)\leq\nu^{k}.$ (8)

In the case $k=0$, this bound is, generally speaking, wrong, since an
essential portion of the total weight remains in $I_{0}$.

However, if $I_{0}$ receives exceptional weight, then other sets
$I_{1},I_{2},\ldots$ (in general) suffer from the weight deficit, as follows
from (8). Therefore, we can redistribute the exceptional weight from $I_{0}$
to other sets, in accordance to the error-to-the-right principle. It is
possible, since the composition preserves the conditional total weight, i.e.
the sum of distribution components.

As a consequence, after $t$ steps of composition with coverings $F/F^{T}$, any
set $I_{k}$ contains a portion at most $\nu^{k}$ of the total weight
$\sigma^{t}(F)$ of the resulting covering of $A^{\otimes t}$, up to some
errors to the right.

$\blacktriangleright$ In our example with $D_{4}$, we may choose
$\nu=\sqrt{3}/2$ (it’s a root of the polynomial $P_{F_{2}}(x)$). The choice
allows us to avoid a redistribution of the weight from $I_{0}$. The
distribution $\\{1,\nu,\nu^{2},\ldots,\nu^{k},\ldots\\}$ is stationary for the
diagram on Fig. 2. Implicitly, we also chose
$\lambda=\log_{4}\nu>\lambda_{F_{2}}$. $\blacktriangleleft$

II. The redistribution of the spectral weight of a set of rectangles under the
composition with type-$G$ coverings (again, with possible errors to the right)
may be described by the polynomial

$P_{G}(x)=\alpha_{l}x^{l}+\ldots+\alpha_{1}x+\alpha_{0}=\frac{1}{\sigma(G)}\sum_{R\in
G}\sigma(R)x^{\lfloor\log_{\tau}\rho(R)\rfloor}.$

For $R\in I_{m}$, the spectral weight of rectangles from $(G/G^{T})\otimes R$
is distributed so that the portion $\alpha_{k}$ belongs to $I_{m-k}$ in the
case $k<m$, and to $I_{0}$, otherwise.

$\blacktriangleright$ Actually, for the covering $G_{2}$, and $\tau=4$, one
can assign (assuming some errors to the right) $P_{G_{2}}(x)=(x+2)/3$. The
weight redistribution under the action of coverings $G_{2}/G_{2}^{T}$ is shown
by the diagram on Fig. 3. $\blacktriangleleft$

$I_{0}$$I_{1}$$I_{2}$$I_{k}$$\scriptstyle{2/3}$$\scriptstyle{2/3}$$\scriptstyle{2/3}$$\scriptstyle{2/3}$$\scriptstyle{2/3}$$\scriptstyle{1/3}$$\scriptstyle{1/3}$$\scriptstyle{1/3}$$\scriptstyle{1}$……
Figure 3: Diagram of spectral weight redistribution under the action of the
composition with the covering $G_{2}/G_{2}^{T}$

Assuming that the initial weight distribution has a form $q_{m}(0)=1$,
$q_{i}(0)=0$ for all $i\neq m$, that is, all rectangles are located in
$I_{m}$, then after $t$ steps of compositions with $G$-type coverings, we
obtain a distribution999Here $C_{t}^{k_{1},\ldots,k_{l}}$ stands for the
multinomial coefficient representing the number of ways to select from $t$
elements $l$ groups, with $k_{i}$ elements in $i$-th group.

$q_{m-k}(t)=\sum_{\begin{array}[]{c}\scriptstyle k_{1}+\ldots+k_{l}\leq t\\\
\scriptstyle
k_{1}+2k_{2}+\ldots+lk_{l}=k\end{array}}C_{t}^{k_{1},\ldots,k_{l}}\alpha_{1}^{k_{1}}\cdot\ldots\cdot\alpha_{l}^{k_{l}}\alpha_{0}^{t-(k_{1}+\ldots+k_{l})}$
(9)

for $0\leq k<m$. For a component associated with $I_{0}$, we use a trivial
estimate $q_{0}(t)\leq 1$. By consideration, if $\alpha_{i}=0$, then
$k_{i}=0$, and the corresponding factor $\alpha_{i}^{k_{i}}$ in (9) should be
replaced by 1. This remark will be implied in further calculations. The total
weight of the (considered part of the) covering under construction will
increase $\sigma^{t}(G)$ times in $t$ steps.

Note that $P_{G}(1/\sqrt{\tau})=\pi_{G}(\tau)$, see (4). Recall that
$\pi_{G}(x)\to\mu_{G}$ as ${x\to 1}$. Our final choice of $\tau$ is such that
the inequality (5) remains valid after replacing $\lambda_{F}$ by $\lambda$,
and $\mu_{G}$ by $\pi_{G}(\tau)$ ($\tau$ should be small enough).

Observe that our choice implies $\alpha_{0}\neq 1$, since otherwise
$\pi_{G}(\tau)\equiv 1$, and the inequality (5) cannot be satisfied. On the
other hand, $\mu_{G}<1$ due to the compactness of $G$: there should exist
rectangles $R\in G$ with $\rho(R)>1$.

III. Now we are ready to state the synthesis algorithm, namely the rule of
combination of the coverings $F$ and $G$. Choose $\gamma$ satisfying the
condition

$-\frac{1}{\lambda}\log_{\tau}\frac{\sigma(G)}{\sigma(F)}<\gamma\leq-2\log_{\tau}\pi_{G}(\tau).$
(10)

Such $\gamma$ does exist, since the inequality between the left and the right
sides of (10) is equivalent to (5), where $\lambda_{F}$ is replaced by
$\lambda$, and $\mu_{G}$ is replaced by $\pi_{G}(\tau)$ (just apply a
base-$\tau$ logarithm to (5) and divide by $\lambda$).

To construct the required covering of $A^{\otimes n}$, we assign two sets of
rectangles, $\cal F$ and $\cal G$.

$(i)$ Before the start of the algorithm, the set $\cal G$ is empty, and the
set $\cal F$ contains a $1\times 1$ rectangle (a trivial covering of the
matrix $A^{\otimes 0}$). $(ii)$ Then, perform $n$ similar steps. A step $t$
does the following: — with the use of a suitable covering $F$ or $F^{T}$,
transform any rectangle $R\in\cal F$ into $(F/F^{T})\otimes R$; — with the use
of a suitable covering $G$ or $G^{T}$, transform any rectangle $R\in\cal G$
into $(G/G^{T})\otimes R$; — relocate from $\cal F$ to $\cal G$ all rectangles
$R\in\cal F$ belonging to the sets $I_{m}$ with $m\geq\gamma(n-t)$. $(iii)$ By
construction, after any $t$ steps, the set ${\cal F}\cup{\cal G}$ is a
covering of the matrix $A^{\otimes t}$. In the end of the algorithm, the set
$\cal F$ is empty, and $\cal G$ is a covering of $A^{\otimes n}$.

Essentially, this is the algorithm [1]. The parameter $\gamma$ controls the
switching between the two stages of the algorithm. Next, we are going to prove
that $\cal G$ has the desired complexity.

For any $m$, we have to relocate rectangles belonging to $I_{m}\cap\cal F$ on
at most $\frac{d}{\gamma}+2=O(1)$ (consecutive) steps, namely while
$m-d<\gamma(n-t)\leq m$, and one more time, when $m-d\geq\gamma(n-t)$. On the
subsequent steps, rectangles in $I_{m}$ don’t appear.

Let us turn to complexity bounds. The complexity of a rectangle $R\in I_{m}$
in relation to its spectral weight is estimated to be higher, when $m$ is
greater. Recall that $w(R)\asymp\sqrt{\rho(R)}\sigma(R)$. Therefore, an
erroneous assignment of a rectangle to a set $I_{m}$ with a higher index $m$
(an error to the right) leads to overestimation of the complexity.

The complexity of the part of the covering $\cal G$ derived from rectangles
that belonged to $I_{m}\cap\cal F$ at the moment of relocation is bounded from
above as

$L_{m}\preceq\sigma^{n}(F)\nu^{m}\left(\frac{\sigma(G)}{\sigma(F)}\right)^{\frac{m}{\gamma}}\cdot\\\
\cdot\left(1+\sum_{k_{1}+\ldots+k_{l}\leq\frac{m}{\gamma}}C_{\frac{m}{\gamma}}^{k_{1},\ldots,k_{l}}\alpha_{1}^{k_{1}}\cdot\ldots\cdot\alpha_{l}^{k_{l}}\alpha_{0}^{\frac{m}{\gamma}-(k_{1}+\ldots+k_{l})}\tau^{\frac{m-(k_{1}+2k_{2}+\ldots+lk_{l})}{2}}\right).$
(11)

Here the first factor $\sigma^{n}(F)$ is the expected weight of the entire
covering under the (optimistic) assumption that we apply only $F$-type
coverings. The second factor $\nu^{m}$ is the upper bound for the weight
portion of the covering $\cal F$ associated with $I_{m}$ at the moment of
relocation. The third factor reflects the weight increase caused by the
application of $G$-type coverings instead of $F/F^{T}$ on the last
$m/\gamma-O(1)$ steps of the algorithm. In brackets, an additional factor
taking into account the final distribution of the spectral weight is written.
Namely, the term 1 represents the complexity of rectangles from $I_{0}$, and
under the sum are written the products of the partial weight portions of
rectangles from $I_{m-k}$ with $k=k_{1}+2k_{2}+\ldots+lk_{l}$ (provided by
(9)), and the estimate $\tau^{\frac{m-k}{2}}$ of the ratio between the
complexity and the spectral weight of rectangles from $I_{m-k}$. The summands
with $m\leq k$ are excess.

Multinomial coefficients satisfy the standard inequality101010Easily follows
by induction on $l$ from the well-known relation $C_{n}^{k}\leq 2^{H(k/n)}$.

$C_{n}^{k_{1},\ldots,k_{l}}\leq
2^{nH(k_{1}/n,\ldots,k_{l}/n)},\qquad\text{where}$
$H(x_{1},\ldots,x_{l})=-\sum_{i=0}^{l}x_{i}\log_{2}x_{i},\qquad
x_{0}=1-\sum_{i=1}^{l}x_{i},$

is the binary entropy function defined on the $l$-dimensional simplex111111On
the boundary, the function is defined by continuity.
${\mathbb{R}}_{+}^{l}\cap\\{x_{1}+\ldots+x_{l}\leq 1\\}$.

It is easy to check that for any $c_{i}>0$,

$H(x_{1},\ldots,x_{l})+\sum_{i=1}^{l}x_{i}\log_{2}c_{i}\leq\log_{2}\left(1+\sum_{i=1}^{l}c_{i}\right).$
(12)

Indeed, it immediately follows from the variant of the Hölder’s inequality

$\prod_{i=0}^{l}a_{i}^{b_{i}}\leq\sum_{i=0}^{l}a_{i}b_{i}$

that holds for $a_{i},b_{i}>0$, and $\sum b_{i}=1$, see e.g. [5, Ch. V] (just
assign $b_{i}=x_{i}$, $a_{i}=c_{i}/x_{i}$, where $c_{0}=1$, and take a
logarithm; in the case $x_{i}=0$ for some $i$, evaluate the limit).

By setting $k_{i}=\frac{x_{i}m}{\gamma}$ for $i=1,\ldots,l$, and applying
(12), we obtain

$\log_{2}\left(C_{\frac{m}{\gamma}}^{k_{1},\ldots,k_{l}}\alpha_{1}^{k_{1}}\cdot\ldots\cdot\alpha_{l}^{k_{l}}\alpha_{0}^{\frac{m}{\gamma}-(k_{1}+\ldots+k_{l})}\tau^{\frac{m-(k_{1}+2k_{2}+\ldots+lk_{l})}{2}}\right)\leq\\\
\frac{m}{\gamma}\left(H(x_{1},\ldots,x_{l})+\sum_{i=1}^{l}x_{i}\log_{2}\frac{\alpha_{i}}{\alpha_{0}\tau^{i/2}}+\log_{2}\left(\alpha_{0}\tau^{\gamma/2}\right)\right)\leq\\\
\frac{m}{\gamma}\cdot\left(\log_{2}\left(\sum_{i=0}^{l}\frac{\alpha_{i}}{\tau^{i/2}}\right)+\log_{2}\left(\tau^{\gamma/2}\right)\right).$

Thus, we continue (11) as

$L_{m}\preceq\sigma^{n}(F)\left(C_{0}^{\frac{m}{\gamma}}+m^{l}C_{1}^{\frac{m}{\gamma}}\right),\quad\text{where}$
$C_{0}=\frac{\sigma(G)\nu^{\gamma}}{\sigma(F)},\qquad C_{1}=C_{0}\cdot
P_{G}\left(1/\sqrt{\tau}\right)\cdot\tau^{\gamma/2}.$

Here the power of $C_{0}$ corresponds to the contribution of rectangles from
$I_{0}$, and the power of $C_{1}$ to the contribution of the remaining
rectangles.

Let us check that $C_{1}\leq C_{0}<1$. Indeed, in the case
$\sigma(G)=\sigma(F)$, the inequality $C_{0}=\nu^{\gamma}<1$ holds trivially.
In the other case $\sigma(G)>\sigma(F)$, from the left part of (10), it
follows that

$C_{0}=\frac{\sigma(G)}{\sigma(F)}\,\tau^{\lambda\gamma}<1.$

Further, the right part of (10) implies $\tau^{\gamma/2}\leq 1/\pi_{G}(\tau)$,
hence $C_{1}\leq C_{0}$.

$\blacktriangleright$ For our example $A=D_{4}$, choose $\gamma=1/5$. Then
$C_{0}<0.99$, and $C_{1}<0.95$. $\blacktriangleleft$

Finally, we conclude

$\mathsf{L}_{2}(A^{\otimes n})\preceq w({\cal G})=\sum_{m\geq
0}L_{m}\preceq\sigma^{n}(F)\sum_{m\geq
0}\left(C_{0}^{\frac{m}{\gamma}}+m^{l}C_{1}^{\frac{m}{\gamma}}\right)\asymp\sigma^{n}(F).$

$\blacksquare$

As a reserve for improving the method, one can suggest a more subtle
combination of the two types of coverings involving reverse relocations of
rectangles from $\cal G$ to $\cal F$.

##  3 Coverings of the Kneser–Sierpinski matrices

For a matrix $D_{r}$, $r=2^{t}$, we propose a covering $F_{t}$ consisting
solely of rectangles of width 1. It generalizes the examples of Fig. 1 (for
$t=2$), and from [1] (for $t=3$).

We exploit a simple gradient-fashion approach. First, put into $F_{t}$ a
column labeled by $\emptyset$, then add a row labeled by $\emptyset$ from the
remaining part of the matrix. Next, we sequentially extract ones from all
columns, and then from all rows labeled by the size-1 subsets of $[t]$. Then,
we do the same with columns and rows labeled by the size-2 subsets, and go on
until we reach $t/2$-size labels. At this point, all ones in $D_{r}$ are
covered.

Let $s(m,k)$ denote the binomial sum

$s(m,k)=C_{m}^{k}+C_{m}^{k+1}+\ldots+C_{m}^{m}.$

By construction, any rectangle corresponding to a column labeled by a size-$k$
subset has height $s(t-k,k)$, and a rectangle corresponding to a row labeled
by a size-$k$ subset has length $s(t-k,k+1)$. Hence,

$\sigma(F_{t})=\sum_{k=0}^{t/2}C_{t}^{k}\left(\sqrt{s(t-k,k)}+\sqrt{s(t-k,k+1)}\right).$

Direct calculation shows that the quantity $\log_{r}\sigma(F_{t})$ attains its
minimum $1.2502...$ when $t=18$, and, as easy to verify, it tends to the limit
$1.259...$ as $t\to\infty$.

As a correcting covering $G_{t}$, we take a covering of all columns of $D_{r}$
by individual rectangles, as shown on Fig. 1. Easy to see that
$\sigma(G_{t})=(\sqrt{2}+1)^{t}$, and
$\mu_{G_{t}}=\left(\frac{2}{\sqrt{2}+1}\right)^{t}$.

With the use of coverings $F_{t},G_{t}$, it is possible to satisfy the
conditions of Theorem 1 only for $t\leq 15$. It can be directly verified that
$\sigma(F_{15})<442412$, and $\lambda_{F_{15}}<-0.04$.

######  Corollary 1.

$\mathsf{SUM}_{2}(D_{N})\preceq N^{\log_{2^{15}}\sigma(F_{15})}\prec
N^{1.251}$.

It’s a kind of surprise, that coverings of matrices $D_{r}$ by width-1
rectangles appear so efficient for the iterative procedure. However, they are
not optimal in terms of spectral weight. At least starting from $t=7$, one can
construct better coverings via uniting common parts of columns or rows.

The author thanks Stasys Jukna for helpful comments.

## References

  * [1] Alman J., Guan Y., Padaki A. Smaller low-depth circuits for Kronecker powers. 2022\. arXiv:2211.05217v1.
  * [2] Chistikov D., Iván S., Lubiw A., Shallit J. Fractional coverings, greedy coverings, and rectifier networks. Proc. STACS (Hannover, 2017). LIPIcs. 2017. 66, Art. 23.
  * [3] Jukna S., Sergeev I. Complexity of linear boolean operators. Foundations and Trends in Theoretical Computer Science. 2013. 9(1), 1–123.
  * [4] Langer R. E. On the zeros of exponential sums and integrals. Bulletin AMS. 1931. 37, 213–239.
  * [5] Mitrinović D. S., Pečarić J. E., Fink A. M. Classical and new inequalities in analysis. Dordrecht–Boston–London: Kluwer, 1993\.


tablenum SIXtablenum

††thanks<EMAIL_ADDRESS>
# Semantic Segmentation of Solar Radio Spikes at Low Frequencies

Pearse C. Murphy⋆ LESIA, Observatoire de Paris, Université PSL, Sorbonne
Université, Université Paris Cité, CNRS, 92190 Meudon, France  Stéphane
Aicardi  DIO, Observatoire de Paris, Université PSL, CNRS, 75014 Paris, France
Baptiste Cecconi  LESIA, Observatoire de Paris, Université PSL, Sorbonne
Université, Université Paris Cité, CNRS, 92190 Meudon, France  Carine Briand
LESIA, Observatoire de Paris, Université PSL, Sorbonne Université, Université
Paris Cité, CNRS, 92190 Meudon, France  Thibault Peccoux  LESIA, Observatoire
de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS,
92190 Meudon, France  IPSA, 94200 Ivry-sur-Seine, France

###### Abstract

Solar radio spikes are short lived, narrow bandwidth features in low frequency
solar radio observations. The timing of their occurrence and the number of
spikes in a given observation is often unpredictable. The high temporal and
frequency of resolution of modern radio telescopes such as NenuFAR mean that
manually identifying radio spikes is an arduous task. Machine learning
approaches to data exploration in solar radio data is on the rise. Here we
describe a convolutional neural network to identify the per pixel location of
radio spikes as well as determine some simple characteristics of duration,
spectral width and drift rate. The model, which we call SpikeNet, was trained
using an Nvidia Tesla T4 16 GPU with $\sim 100\,000$ sample spikes in a total
time of 2.2 hours. The segmentation performs well with an intersection over
union in the test set of $\sim 0.85$. The root mean squared error for
predicted spike properties is of the order of $23\%$. Applying the algorithm
to unlabelled data successfully generates segmentation masks although the
accuracy of the predicted properties is less reliable, particularly when more
than one spike is present in the same $64\times 64$ pixel time-frequency
range. We have successfully demonstrated that our convolutional neural network
can locate and characterise solar radio spikes in a number of seconds compared
to the weeks it would take for manual identification.

##  1 Introduction

Solar radio bursts occur over a wide frequency range with varying spectral and
temporal characteristics. At decametric wavelengths there are five canonical
types of radio burst which can have durations ranging from seconds to hours
and spectral widths up to tens of . Alongside these, there exist a wide
variety of short duration, narrow spectral width bursts. These bursts
generally fall into three categories; solar radio spikes (Tarnstrom & Philip,
1972), S-bursts (Ellis, 1969) and striae (de La Noe & Boischot, 1972), with
individual striae and spikes sharing a number of characteristics (Clarkson et
al., 2021). These short lived bursts can give an insight into energy transfer
that takes place over a number of milliseconds (e.g. Clarkson et al., 2021) as
well as the turbulent nature of plasma in the corona (e.g. Kontar et al.,
2017; Reid & Kontar, 2021). They may also offer a way of remotely measuring
the magnetic field in the corona (e.g. Melnik et al., 2010; Clarke et al.,
2019). This underscores the importance of understanding the origins of short
duration, narrow bandwidth solar radio bursts.

Modern state of the art radio telescopes such as the LOw Frequency ARray
(LOFAR; van Haarlem et al., 2013) and its extension, the New Extension in
Nançay Upgrading LOFAR (NenuFAR; Zarka et al., 2012) have the temporal and
spectral resolution necessary to analyse short duration, narrow bandwidth
solar radio bursts in great detail (Morosan et al., 2015; Clarke et al., 2019;
Kontar et al., 2017; Sharykin et al., 2018; Clarkson et al., 2021, 2023,
Briand et al. in prep). However, the sporadic occurrence of spike bursts
coupled with the vast data rates of the high resolution observations means
that finding them manually can be and arduous and time consuming task. Modern
machine learning methods are the ideal solution for this as, once the model is
trained, they eliminate the need for human identification of spikes and
determining their characteristics. The application of machine learning
algorithms to solar radio data has seen increased popularity in recent
literature. Already, we have seen the application of the YOLO (you only look
once; Redmon et al., 2016) convolutional neural network (CNN) to detecting
solar radio spikes in the [range-phrase=–]1.11.34 range (Lv et al., 2023),
while Scully et al. (2021, 2023) used YOLO to detect Type III radio bursts in
the [range-phrase=–]1090 range.

YOLO performs what is known as object detection, where an object is located in
an image and a bounding box which determines its height and width is drawn
around it. The task of locating an object in an image on a per pixel basis is
known as semantic segmentation. The varied and irregular morphology of solar
radio spikes make it a problem particularly suited for the use of semantic
segmentation. Murphy et al. (2023) used a CNN called UNET (Ronneberger et al.,
2015) to perform semantic segmentation on NenuFAR dynamic spectra to
automatically detect the presence of solar radio emission. Here we adapt the
work of Murphy et al. (2023) to perform segmentation for solar radio spikes
and determine their location, duration, spectral width and drift rate in
NenuFAR dynamic spectra.

The design and training procedure for our model is described in Section 2. We
describe the results of the training in Section 3 and give further discussion
on the application of our CNN in Section 4. We conclude with a summary in
Section 5.

##  2 Method

Our goal is to develop a CNN to automatically generate a mask of where a solar
radio spike occurred in the time-frequency domain as well as predict the
spike’s central location, duration, spectral width and drift rate. In order to
achieve this, we must first perform pre-processing of the NenuFAR dynamic
spectra in order to create a dataset suitable for training the CNN.

###  2.1 Preparing the datasets

The time and frequency location, spectral width, duration and drift rate of
$1000$ spike bursts were identified manually by summer student Thibault
Peccoux, to whom we owe a debt of gratitude for this tedious work. The spikes
were located between 10 and 90 in NenuFAR observations on 2022-02-02,
2022-05-29, 2023-05-02, 2023-06-01 and 2023-07-10. These are the same data
used by Briand et al. (in prep). Figure 1 shows a dynamic spectrum of a number
of radio spikes occurring on 2022-02-02 over ∼9 in the frequency range [range-
phrase=–]2560. The figure also shows radio frequency interference (RFI) which
is a common feature of low frequency observations and thus our dataset must
incorporate this.

Figure 1: A sample of solar radio spikes on 2022-02-02. The colourbar is
intensity relative to the background. The large number of spikes in such a
short time range indicate the difficulty in manually annotating their every
occurrence. Radio frequency interference and the end of a type III solar radio
burst in the frequency range [range-phrase=–]2535MHz also pose a challenge.

These spikes formed the basis for the training, validation and test datasets.
We performed a stratified split based on the date of observation so that each
split had the same proportion of spikes that were identified on each date. The
percentage of spikes in each split were 64% training, 16% validation and 20%
test. The inputs to the CNN are $64\times 64$ pixel Stokes I and Stokes V/I
dynamic spectra which we hereafter refer to interchangeably as tiles. One time
pixel corresponds to ∼21 and one frequency pixel is ∼6.1. This results in a
time duration of 1.34 and a frequency range of 390.625 for each tile. CNNs
often see improved performance when the input values lie within a small range,
thus we normalised these to the background and rescale in the range $01$. The
frequency range of the dynamic spectrum in was also included as an input. In
order to generate the ground truth output for the model, a segmentation mask
was produced for each spike burst using the CUSUM-slope method described by
Murphy et al. (2023).

To increase the number of samples available for training, we applied random
shifts in time and frequency to the spikes in the training and validation
sets. The same shifts were also applied to the segmentation masks. In an
effort to improve the robustness of the model when dealing with unseen data,
we also located times where no spike occurred to produce samples of the
background and added these to the training set. These also included samples in
common RFI bands $\lesssim$ 30 and at ∼70. In total the training set comprised
of $90\,435$ samples and the validation set $16\,000$. The test set did not
undergo any augmentation and contained 200 samples.

For each spike in the dataset we have a list of its characteristics, i.e., its
location in time and frequency, its duration, its spectral width and its drift
rate. We determine the time and frequency location in terms of pixels while
duration, spectral width and drift rate are in physical units of , and
respectively.

Figure 2: Architecture of our model. The top half of the model is the UNET
adapted from Ronneberger et al. (2015). We pass the outputs from UNET’s
bottleneck layer to a fully connected network. This in turn has two outputs;
the spike characteristics and the objectness.

Finally, we determined an objectness score for each sample. This is used to
determine how certain the model is that a spike exists in the sample at all.
An objectness of 0 means the model is certain there is no spike while an
objectness of 1 means the model is certain there is a spike. We determine the
objectness of each sample as the ratio of total number of pixels in the
shifted segmentation mask to that of the segmentation mask before it was
shifted. Thus, any spike that is not fully in a tile will have an objectness
score less than one.

###  2.2 Model architecture and training the model

The primary output of our model, affectionately referred to as SpikeNet, is a
mask of every pixel where a solar radio spike occurred in time and frequency.
As mentioned in Section 1, this type of task is known as semantic segmentation
and a common CNN used in this regard is called UNET. We use the implementation
of UNET in the segmentation_models Python library (Iakubovskii, 2019) as a
starting point.

As well as trying to produce a segmentation mask from these inputs, we require
our model to determine the central location in time and frequency, duration,
spectral width and drift rate of the spike. Thus, in addition to the decoder
part of the UNET, we included a fully connected deep neural network with 5
hidden layers of 30 neurons each to determine these characteristics, as well
as the objectness score. The model architecture is depicted in Figure 2. We
take the values in the bottleneck layer of the original UNET to use as the
input to the fully connected network. The final hidden layer is connect to two
separate output layers. This then gives two separate outputs, the 5 spike
characteristics and the objectness score.

We apply a different loss for each of the outputs; binary cross entropy for
the segmentation performed by UNET, mean squared error for the regression of
spike characteristics and binary cross entropy for the objectness. The total
loss measured during training was the weighted sum of each of these losses. We
give a weight of 10 to the regression loss and a weight of 1 to segmentation
and objectness loss. The model was trained using a linearly increasing
training rate for the first 5 epochs before implementing a cosine decay. The
total validation loss was monitored per epoch and the training was stopped
once it failed to decrease for three epochs in a row. The model was trained
for a total of 26 epochs which took 2.2 hours on an Nvidia Tesla T4 16 GPU.

Figure 3: Evolution of loss with respect to training epoch. The training
(blue) and validation (orange) loss is shown in each panel for the total loss,
segmentation loss, regression loss and objectness loss.

##  3 Results

The values of the individual and total loss functions per epoch are shown in
Figure 3. The loss for the training set (blue) decreased rapidly in all cases
while the validation loss (orange) decreased slowly before forming a plateau.
The gap between the training and validation loss could indicate under fitting
or that the data samples in the training set were not representative of those
in the validation set.

The key metric to determine the accuracy of our segmentation masks to the
ground truth is the intersection over union (IoU) or Jaccard index. It is
defined, as the name suggests, as

$\text{IoU}(A,B)=\frac{|A\cap B|}{|A\cup B|},$ (1)

where $A$ is the ground truth and $B$ is the prediction. The top panel of
Figure 4 shows the IoU for the validation set reaching $\sim 85\%$. The bottom
panel of Figure 4 shows that the root mean squared (RMS) error of the
regression to predict spike characteristics decreases with each epoch.
Comparing the mean RMS for all characteristics of 0.0736 to the mean of the
true characteristics gives an approximate error of $23\%$ for predicted
characteristics.

Figure 4: The evolution of model metrics with increasing training epoch. Panel
a shows the intersection over union for the training (blue) and validation
(orange) sets. Similarly, the root mean squared error is shown in the bottom
panel.

It is more interesting, and perhaps more informative, to take a detailed look
at the model outputs in Figure 5. Here we show a batch of 32 samples from the
validation set and overlay the centre of the spike (red circle/orange
triangle), a bounding box described by the total duration and spectral width,
and the drift rate. We also overlay the predicted segmentation masks in white.
As we expected from the histograms in Figure 6, the predicted characteristics
mostly match well with the ground truth but not for every sample. In
particular, any time there is more than one spike in the tile, the predicted
location is often somewhere between the two spikes. This has likely occurred
because the ground truth segmentation masks and spike characteristics assume
only one spike is present in each tile thus the model hasn’t learned how to
handle multiple spikes. We discuss possible methods to correct for this in the
next section.

Figure 5: A comparison of the model performance compared to the ground truth.
Each panel represents a sample from the validation set and depicts the dynamic
spectrum of the spike burst with an overlay of the predicted segmentation
mask. The ground truth spike position is given by a red circle while an orange
triangle denotes the predicted position. The bounding box, determined by the
spike duration and spectral width, as well as its drift rate is also plotted
in red for ground truth and orange for model prediction. The axes of the
panels denote time (x axis) and frequency (y axis) in pixels.

The predicted segmentation masks are output as polygons in the TFCat (Cecconi
et al., 2023) format. A selection of these are currently available at
https://doi.org/10.25935/m5cq-f460. This can be queried using the Astronomical
Data Query Language (ADQL, Yasuda et al., 2004). The scripts to prepare the
data and train the model are available at
https://gitlab.obspm.fr/pmurphy/spikenet.

##  4 Discussion

Machine learning methods are being applied to solar radio data in a growing
number of use cases, as outlined in Section 1. We have found success with our
CNN in segmenting and characterising solar radio spikes. The high IoU scores
obtained during training indicate an accurate recreation of the ground truth
segmentation masks. Although Figure 5 shows a good overall agreement between
the ground truth and predicted and characteristics the model still gives a
mean RMS error of $\sim 23\%$. In Figure 6 we investigate the effect of this
error on the model’s predictions by comparing the histogram of each spike
characteristics for the ground truth (blue) and predicted value (orange) of
the validation set. The distribution of predicted central times and
frequencies are similar to the ground truths and the predicted durations
follow the same overall shape as the ground truths. The distributions for
spectral width and drift rate however, do not fully agree with those of the
ground truth. From this, we infer that the dominant source of error is related
to the predicted spectral width.

The histograms of the ground truth characteristics can also guide us to where
the dataset selection for our model may need improvement. For example, we see
peaks in the ground truth central time and frequencies at 0, 0.5 and 1. The
peaks at 0 and 1 are likely a result of samples where the central point lies
outside of the tile while the central peak is due to the validation set
including a sample of every spike in the centre of the tile as well as any
samples that have undergone a very small shift during augmentation described
in Section 2. The shape of the central time distribution can be attributed to
the fact that the random shift applied during augmentation was based on the
duration of the spike burst. The frequency shifts, on the other hand, were a
result of a uniform random distribution. Perhaps the accuracy of the
regression would improve if the time shifts were drawn from a uniform
distribution as well as it will be more used to seeing spikes near the edges
of a tile.

Figure 6: A comparison of the histograms for the ground truth and predicted
spike characteristics. The time and frequency position are determined as a
fraction of the width (height) of the dynamic spectrum in terms of pixels. We
see an abundance of spikes in the validation set have centres outside the
width of the tile, resulting in the peaks at 0 and 1. There is also always 1
sample with the spike directly centred in the tile which may explain the peak
at 0.5 in time and frequency location.

As mentioned above, the characterisation of spikes is noticeably less accurate
when multiple spikes are present in a tile. Performing semantic segmentation
for distinct objects is known as instance segmentation. One network that
performs well in these kinds of tasks is called Mask R-CNN (Region based CNN;
He et al., 2018). Mask R-CNN utilises a region proposal network to locate
areas in each image which likely includes the desired object. We borrow the
idea of objectness from Mask R-CNN and implement it in our network as a
measure of how sure the model is that it “sees” a spike. Hou et al. (2020)
have adapted the region proposal network of Faster R-CNN (Ren et al., 2016) to
locate solar radio spikes in the [range-phrase=–]150500 range. It would be a
significant undertaking for the authors of this work to implement a similar
network in the [range-phrase=–]1090 range, particularly because there is no
reference in Hou et al. (2020) as to where the code to train their network
exists.

###  4.1 Application to unseen, unlabelled data

The ultimate test for our machine learning model is to apply it to a dynamic
spectrum that does not exist in the training, validation or test dataset and
thus does not have any ground truth masks or spike characteristics. This is
instructive to qualitatively determine the robustness of the model to unseen
data and inform how the model can be further developed and put into practical
use. Figure 7 shows a Stokes I dynamic spectrum between 11:03:59 and 11:04:05
on 2022-02-02 in the frequency range [range-phrase=–]47.550 over 4 seconds.
Overlayed in transparent white are the segmentation masks as well as the
central locations, marked by orange triangles, duration and spectral width,
denoted by orange bounding boxes, and drift rate, orange lines.

Figure 7: Performance of the model on unseen data. The panel shows a Stokes I
dynamic spectrum at the same frequency and temporal resolution as Figure 1.
Overlayed on the dynamic spectrum is the segmentation mask detected by the
CNN. Orange triangles denote the predicted centre position of the spikes while
their duration, spectral width, and drift rate are shown as the orange
bounding boxes and lines respectively. The segmentation masks capture the core
of each spike well but fail when a spike extends across a tile.

We see that the core of the spikes are well captured by the predicted
segmentation masks but the model fails to predict most of the spikes duration.
Similarly, the central point and bounding boxes do not cover the full spike.
We attribute this the the spikes not falling fully into the tiles that the
original dynamic spectrum was divided into. Similar to the problem of multiple
spikes appearing in one tile, this behaviour could be correct by using a
region proposal network. Nevertheless, the performance on unseen data is
promising and, if used with appropriate caution, could allow for the analysis
of solar radio spikes without having to manually find them.

##  5 Conclusion

We have trained a machine learning model on $\sim 100\,000$ samples of solar
radio spikes using an Nvidia Tesla T4 16 GPU. The training took 2.2 hours in
total and resulted in a high $\sim 80\%$ IoU and relatively low $<0.1$ RMS
error. The model successfully produced segmentation masks for radio spikes in
the validation set and predicted the spike characteristics, i.e. location in
time and frequency, duration, spectral width and drift rate. The model does
have some shortcomings when applied to unseen data. For example, a certain
number of spike masks are abruptly cut off and the model was less successful
in determining the position of spikes when more than one spike appeared in the
$64\times 64$ pixel dynamic spectrum sample. This is due to the original
dynamic spectrum being divided into $64\times 64$ tiles thus, any burst not
fully in one tile is unlikely to be fully recovered in the mask. A region
proposal network to first determine where spikes are likely to be is the best
solution for this.

With the current data rate generation of LOFAR and NenuFAR, not to mention the
ongoing development of next generation low frequency arrays such as the Square
Kilometre Array (SKA) which are set to produce petabytes of data, the ability
to automatically detect fine scale features in solar radio observation has
become increasingly important. This task is ideally suited to machine learning
methods, the latest of which we have presented here. Our CNN drastically
reduces the time to locate and identify solar radio spikes from minutes to
less than a second per spike. While the accuracy of the predicted spike
characteristics could be improved, the use of our CNN trivialises finding
spikes hidden in solar radio data and allows for more time to be dedicated to
their analysis using whichever methods a particular researcher desires. Future
development of SpikeNet could see it used in near real-time on NenuFAR solar
observations to dynamically determine an appropriate time resolution. For
example, should SpikeNet detect a certain number of spikes over a given time
period, keep this data at its native resolution. Otherwise, the time
resolution should be rebinned. This will alleviate pressure on the fixed
storage space allocated to solar observations with NenuFAR. The expected data
rates of the SKA (Scaife, 2020) mean that dynamically adapting the recording
resolution could drastically reduce the storage required for solar
observations in the future.

## Acknowledgements

The authors would like to thank Thibault Peccoux for his work in identifying
and characterising solar radio spikes. This work has been supported by the
Europlanet-2024 Research Infrastructure project, which received funding from
the European Union’s Horizon 2020 research and innovation programme under
grant agreement No 871149 and a grant from Paris Observatory science council
under the MINERVA project. We acknowledge the use of the Nançay Data Center
(CDN – Centre de Données de Nançay) facility. The CDN is hosted by the
Observatoire Radioastronomique de Nançay (ORN) in partnership with the
Observatoire de Paris, the Université d’Orléans, the Observatoire des Sciences
de l’Univers d’Orléans (OSUC) and the French Centre National de la Recherche
Scientifique (CNRS). The CDN is supported by the Région Centre-Val de Loire
(département du Cher). The ORN is operated by the Observatoire de Paris,
associated with the CNRS. The research leading to these results has received
funding from the European Union’s Horizon Europe Programme under the EXTRACT
Project, grant agreement n° 101093110.

## References

  * Cecconi et al. (2023) Cecconi, B., Louis, C. K., Bonnin, X., Loh, A., & Taylor, M. B. 2023, Frontiers in Astronomy and Space Sciences, 9
  * Clarke et al. (2019) Clarke, B. P., Morosan, D. E., Gallagher, P. T., et al. 2019, Astronomy & Astrophysics, 622, A204, arXiv: 1901.07424
  * Clarkson et al. (2021) Clarkson, D. L., Kontar, E. P., Gordovskyy, M., Chrysaphi, N., & Vilmer, N. 2021, The Astrophysical Journal, 917, L32, aDS Bibcode: 2021ApJ…917L..32C
  * Clarkson et al. (2023) Clarkson, D. L., Kontar, E. P., Vilmer, N., et al. 2023, The Astrophysical Journal, 946, 33, aDS Bibcode: 2023ApJ…946…33C
  * de La Noe & Boischot (1972) de La Noe, J. & Boischot, A. 1972, Astronomy and Astrophysics, 20, 55
  * Ellis (1969) Ellis, G. 1969, Australian Journal of Physics, 22, 177, publisher: CSIRO PUBLISHING
  * He et al. (2018) He, K., Gkioxari, G., Dollár, P., & Girshick, R. 2018, Mask R-CNN, arXiv:1703.06870 [cs]
  * Hou et al. (2020) Hou, Y. C., Zhang, Q. M., Feng, S. W., et al. 2020, Solar Physics, 295, 1, publisher: Springer Nature
  * Iakubovskii (2019) Iakubovskii, P. 2019, Segmentation models
  * Kontar et al. (2017) Kontar, E. P., Yu, S., Kuznetsov, A. A., et al. 2017, Nature Communications, 8, 1515, publisher: Nature Publishing Group
  * Lv et al. (2023) Lv, P. R., Hou, Y. C., Feng, S. W., Du, Q. F., & Tan, C. M. 2023, Astrophysics and Space Science, 368, 14
  * Melnik et al. (2010) Melnik, V. N., Konovalenko, A. A., Rucker, H. O., et al. 2010, Solar Physics, 264, 103, publisher: Springer Netherlands
  * Morosan et al. (2015) Morosan, D. E., Gallagher, P. T., Zucca, P., et al. 2015, Astronomy & Astrophysics, 580, A65, publisher: EDP Sciences
  * Murphy et al. (2023) Murphy, P. C., Cecconi, B., Briand, C., & Aicardi, S. 2023, accepted: 2023-10-23T10:58:13Z
  * Redmon et al. (2016) Redmon, J., Divvala, S., Girshick, R., & Farhadi, A. 2016, You Only Look Once: Unified, Real-Time Object Detection, arXiv:1506.02640 [cs]
  * Reid & Kontar (2021) Reid, H. A. S. & Kontar, E. P. 2021, Nature Astronomy, 5, 796, arXiv: 2103.08424 Publisher: Nature Research
  * Ren et al. (2016) Ren, S., He, K., Girshick, R., & Sun, J. 2016, Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks, arXiv:1506.01497 [cs]
  * Ronneberger et al. (2015) Ronneberger, O., Fischer, P., & Brox, T. 2015, IEEE Access, 9, 16591, arXiv: 1505.04597 Publisher: Institute of Electrical and Electronics Engineers Inc.
  * Scaife (2020) Scaife, A. M. M. 2020, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 378, publisher: Royal Society Publishing
  * Scully et al. (2021) Scully, J., Flynn, R., Carley, E., Gallagher, P., & Daly, M. 2021, in 2021 32nd Irish Signals and Systems Conference (ISSC) (IEEE), 1–6, arXiv: 2105.13387
  * Scully et al. (2023) Scully, J., Flynn, R., Gallagher, P. T., Carley, E. P., & Daly, M. 2023, Astronomy and Astrophysics, 674, A218, aDS Bibcode: 2023A&A…674A.218S
  * Sharykin et al. (2018) Sharykin, I. N., Kontar, E. P., & Kuznetsov, A. A. 2018, Solar Physics, 293, 115, publisher: Springer Netherlands
  * Tarnstrom & Philip (1972) Tarnstrom, G. L. & Philip, K. W. 1972, Astronomy and Astrophysics, 16, 21, aDS Bibcode: 1972A&A….16…21T
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, Astronomy & Astrophysics, 556, A2, arXiv: 1305.3550
  * Yasuda et al. (2004) Yasuda, N., Mizumoto, Y., Ohishi, M., et al. 2004, 314, 293, conference Name: Astronomical Data Analysis Software and Systems (ADASS) XIII ADS Bibcode: 2004ASPC..314..293Y
  * Zarka et al. (2012) Zarka, P., Girard, J. N., Tagger, M., & Denis, L. 2012, LSS/NenuFAR: The LOFAR Super Station project in Nançay, conference Name: SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics Pages: 687-694 ADS Bibcode: 2012sf2a.conf..687Z


# Law of Large Numbers for Permanents of Random Constrained Matrices

Ghurumuruhan Ganesan

Institute of Mathematical Sciences E-Mail<EMAIL_ADDRESS>HBNI  Chennai

###### Abstract

Permanents of random matrices with independent and identically distributed
(i.i.d.) entries have extensively studied in literature and convergence and
concentration properties are known under varying assumptions on the
distributions. In this paper we consider constrained $n\times n$ random
matrices where each row has a deterministic number $r=r(n)$ zero entries and
the rest of the entries are independent random variables that are positive
almost surely. The positions of the zeros within each row are also random and
we establish sufficient conditions for the existence of a weak law of large
numbers for the permanent, appropriately scaled and centred. As a special
case, we see that if $r$ grows faster than $\sqrt{n},$ then the permanent of a
randomly chosen constrained $0-1$ matrix is concentrated around the van der
Waerden bound with high probability, i.e. with probability converging to one
as $n\rightarrow\infty.$

Key words: Weak Law of Large Numbers, Row Constrained Random Matrices, van der
Waerden bound, Second moment method.

AMS 2000 Subject Classification: Primary: 60C05, 05A16.

##  1 Introduction

Estimating the permanents of random matrices is of great interest from both
theoretical and application perspectives and for matrices with independent and
identically distributed (i.i.d.) entries, many convergence properties are
known under varying assumptions on the distributions. For example, (Rempala
and Gupta, 1999) have studied random matrices with independent entries and
have obtained central limit theorems for the permanent as an extension of
Girko’s theorem for determinants. Later (Tao and Vu, 2009) used concentration
techniques and showed that the permanent of a random $n\times n$ matrix with
entries in $\\{-1,1\\}$ is non-zero with high probability, i.e. with
probability converging to one as $n\rightarrow\infty$ and also estimated the
typical value of the permanent. Recently, (Kwan and Sauermann, 2020) have
obtained an analogue of the Tao-Vu result for the permanents of random
symmetric matrices and for a recent survey of combinatorial properties of
random matrices, we refer to (Vu, 2020).

In this paper we study permanents of $n\times n$ random matrices where each
row has a deterministic number zero entries and the rest of the entries are
i.i.d. random variables that are positive almost surely. The positions of the
zeros within each row are also random and we use the second moment method to
establish sufficient conditions on the mean and second moment that guarantee
the existence of a weak law of large numbers for the permanent, appropriately
scaled and centred. In particular, for $0-1$ matrices, this also provides
conditions under which the permanent of a randomly chosen constrained $0-1$
matrix is concentrated around the van der Waerden bound with high probability.
(For completeness, we recall that the van der Waerden bound (Egorychev, 1981),
(Falikman, 1981) asserts that the permanent of a deterministic $n\times n$
matrix containing exactly $r$ ones in each row and column must be at least
$r^{n}\frac{n!}{n^{n}}$).

The paper is organized as follows. In the rest of the section, we state our
main result Theorem 1 regarding the $L^{2}-$convergence of the permanent of
constrained random matrices. Next in Section 2, we collect preliminary
estimates regarding the mean and the second moment of the permanent of $X$ and
finally, in Section 3, we prove Theorem 1.

### Constrained Random Matrices

For integer $n\geq 1$ let $Z:=[Z_{i,j}]_{1\leq i,j\leq n}$ be a random matrix
whose entries $Z_{i,j},1\leq i,j\leq n$ are independent and identically
distributed (i.i.d.) random variables that are positive a.s. with mean
$\mathbb{E}Z_{1,1}$ and second moment $\delta=\delta(n)=EZ^{2}_{1,1}.$ The
permanent of $Z$ is defined as

$per(Z):=\sum_{\sigma}\prod_{i=1}^{n}Z_{i,\sigma(i)},$ (1.1)

where the summation is over all permutations $\sigma$ of $\\{1,2,\ldots,n\\}.$

We introduce constraints in the number of non-zero entries of $Z$ as follows.
Let $1\leq r_{1},\ldots,r_{n}\leq n$ be deterministic integers and let ${\cal
M}$ be the set of all $n\times n$ matrices that contain exactly $r_{i}$ ones
in the $i^{th}$ row and zeros in the rest of the positions, for each $1\leq
i\leq n.$ Let $X=[X_{i,j}]_{1\leq i,j\leq n}$ be uniformly randomly chosen
from ${\cal M},$ independent of $Z$ and form the termwise product matrix
$Y=[X_{i,j}\cdot Z_{i,j}]_{1\leq i,j\leq n}.$ For $1\leq i\leq n,$ the matrix
$Y$ has exactly $r_{i}$ non-zero entries in row $i$ and we define the random
matrix $Y$ on the probability space $(\Omega,{\cal F},\mathbb{P}).$

The following result determines sufficient conditions for the
$L^{2}-$convergence of the permanent $T_{n}=per(Y)$ of the matrix $Y,$
appropriately scaled and centred. For two sequences $\\{a_{n}\\}_{n\geq 1}$
and $\\{b_{n}\\}_{n\geq 1},$ we say that $b_{n}=o(a_{n})$ if
$\frac{b_{n}}{a_{n}}\longrightarrow 0$ as $n\rightarrow\infty$ and
$b_{n}=O(a_{n})$ if there exists a constant $C>0$ such that $b_{n}\leq
Ca_{n}.$ Throughout constants do not depend on $n.$

###### Theorem 1

Setting $\mu_{n}:=\prod_{i=1}^{n}r_{i}\cdot\frac{\nu^{n}n!}{n^{n}},$ we have
that $\mathbb{E}T_{n}=\mu_{n}.$ Moreover, if $r_{low}=\min_{1\leq i\leq
n}r_{i}$ and $r_{up}=\max_{1\leq i\leq n}r_{i}$ are such that

$\frac{\delta}{\nu^{2}r_{low}}=o\left(\frac{1}{\sqrt{n}}\right)\text{ and
}\frac{\delta}{\nu^{2}r_{low}}-\frac{1}{r_{up}}=o\left(\frac{1}{n}\right)$
(1.2)

then

$\mathbb{E}\left(\frac{T_{n}}{\mu_{n}}-1\right)^{2}\longrightarrow 0\text{ as
}n\rightarrow\infty.$ (1.3)

For illustration, we consider the example of $0-1$ matrices and set
$Z_{i,j}=1$ for all $i,j.$ In this case $\delta=\nu^{2}=1$ and the conditions
in (1.2) then hold if the matrices in ${\cal M}$ are reasonably homogenous in
terms of the number of ones in each row and is satisfied if for example

$(\sqrt{n}-1)\log{n}\leq r_{i}\leq(\sqrt{n}+1)\log{n},$

for each $1\leq i\leq n.$ In the case of complete homogeneity, i.e. if
$r_{i}=r=r(n)$ for each $1\leq i\leq n,$ then
$\mu_{n}=r^{n}\cdot\frac{n^{n}}{n!}$ is the van der Waerden bound and
condition (1.2) is equivalent to $\frac{r}{\sqrt{n}}\longrightarrow\infty$ as
$n\rightarrow\infty.$

To prove Theorem 1, we use the second moment method: We first obtain first and
second moment bounds for the permanent of $Y$ and then show that if (1.2)
holds, then the permanent variance, appropriately scaled, converges to zero as
$n\rightarrow\infty.$ This obtains the desired convergence (1.3).

##  2 Preliminaries

In the following main result of this section we obtain bounds for the first
and second moment of the permanent $T_{n}$ of the matrix $Y$ defined prior to
Theorem 1. All constants throughout do not depend on $n.$

###### Proposition 1

Let $r_{up},r_{low}$ and $\mu_{n}$ be as in Theorem 1 and recall that
$\delta=\mathbb{E}Z_{1,1}^{2}$ and $\nu=\mathbb{E}Z_{1,1}.$ Defining

$\alpha_{up}:=\left(\frac{n(r_{up}-1)}{r_{up}(n-1)}\right)^{n}\text{ and
}\beta_{up}:=\frac{\delta r_{up}(n-1)}{\nu^{2}r_{low}(r_{up}-1)}$ (2.1)

and

$\alpha_{low}:=\left(\frac{n(r_{low}-1)}{r_{low}(n-1)}\right)^{n}\text{ and
}\beta_{low}:=\frac{\delta r_{low}(n-1)}{\nu^{2}r_{up}(r_{low}-1)},$ (2.2)

we have that $\mathbb{E}T_{n}=\mu_{n}.$ Moreover, if
$r_{low}\geq\frac{6\delta}{\nu^{2}},$ then we also have that

$\alpha_{low}\cdot
e^{\beta_{low}-1}\cdot\left(1-\frac{2e}{n^{2}}\right)\leq\frac{\mathbb{E}T^{2}_{n}}{\mu_{n}^{2}}\leq\alpha_{up}\cdot
e^{\beta_{up}-1}\cdot\left(1+\frac{2e}{n^{2}}\right)$ (2.3)

for all $n$ large.

To illustrate, we again consider the example of $0-1$ matrices and set
$Z_{i,j}=1$ for all $i,j.$ If $r_{i}=r=r(n)$ for all $i$ so that the matrix is
homogenous, then

$\alpha_{up}=\alpha_{low}=\left(\frac{n(r-1)}{r(n-1)}\right)^{n}=e\left(1-\frac{1}{r}\right)^{n}(1+o(1))$

and $\beta_{up}=\beta_{low}=\frac{n-1}{r-1}.$ Plugging these into (2.3) and
setting $\theta=\left(1-\frac{1}{r}\right)\cdot e^{\frac{1}{r-1}},$ we then
get that

$\frac{\mathbb{E}T^{2}_{n}}{\mu_{n}^{2}}=e^{-\frac{1}{r-1}}\cdot\theta^{n}\cdot(1+o(1)).$

We remark that for constant $r,$ the term $\theta>1$ strictly and therefore
the second moment increases exponentially in $\theta.$ However, we show in the
next Section that if $r=r(n)$ grows faster than $\sqrt{n},$ then the scaled
variance of $T_{n}$ in fact converges to zero as $n\rightarrow\infty.$

### Proof of Proposition 1

The rows of $X$ are mutually independent and so for any $1\leq i,j\leq n$ we
have that

$\mathbb{P}(X_{i,j}=1)=\frac{{n-1\choose r_{i}-1}}{{n\choose
r_{i}}}=\frac{r_{i}}{n}=:p_{i}$ (2.4)

and for any $1\leq j_{1}\neq j_{2}\leq n$ we also have that

$\mathbb{P}(X_{i,j_{1}}=1,X_{i,j_{2}}=1)=\frac{{n-2\choose r_{i}-2}}{{n\choose
r_{i}}}=\frac{r_{i}(r_{i}-1)}{n(n-1)}.$ (2.5)

For a permutation $\sigma$ of $\\{1,2,\ldots,n\\},$ we define
$R_{\sigma}:=\prod_{i=1}^{n}X_{i,\sigma(i)}Y_{i,\sigma(i)}$ and get from (1.1)
that $T_{n}=per(Y)=\sum_{\sigma}R_{\sigma}$ where the summation is over all
permutations of $\\{1,2,\ldots,n\\}.$ For any $\sigma$ we get from (2.4) that
$\mathbb{E}R_{\sigma}=\prod_{i=1}^{n}\nu
p_{i}=\frac{\nu^{n}\prod_{i=1}^{n}r_{i}}{n^{n}}$ and so
$\mathbb{E}T_{n}=\nu^{n}n!\frac{\prod_{i=1}^{n}r_{i}}{n^{n}}=\mu_{n}.$

_Proof of the upper bound in ( 2.3)_: Expanding
$T_{n}^{2}=\left(\sum_{\sigma}R_{\sigma}\right)^{2}$ we get that

$\mathbb{E}T_{n}^{2}=\sum_{\sigma_{1},\sigma_{2}}\mathbb{E}R_{\sigma_{1}}R_{\sigma_{2}}=\sum_{\sigma_{1}}\sum_{k=0}^{n}\sum_{{\cal
S}:\\#{\cal S}=k}\sum_{\sigma_{2}:\sigma_{1}\cap\sigma_{2}={\cal
S}}\mathbb{E}R_{\sigma_{1}}R_{\sigma_{2}}$ (2.6)

where $\sigma_{1}\cap\sigma_{2}:=\\{1\leq i\leq
n:\sigma_{1}(i)=\sigma_{2}(i)\\}$ is the set of indices where $\sigma_{1}$ and
$\sigma_{2}$ agree. If $\sigma_{1}\cap\sigma_{2}={\cal S}$ and $\\#{\cal
S}=k,$ then from (2.5) and the definition of $\mu_{n},$ we see that

$\displaystyle\mathbb{E}R_{\sigma_{1}}R_{\sigma_{2}}$ $\displaystyle=$
$\displaystyle\prod_{j\in{\cal S}}\delta p_{j}\cdot\prod_{i\in{\cal
S}^{c}}(\nu
p_{i})^{2}\cdot\left(\left(1-\frac{1}{r_{i}}\right)\cdot\frac{n}{n-1}\right)$
$\displaystyle=$
$\displaystyle\left(\frac{\mu_{n}}{n!}\right)^{2}\cdot\left(\frac{\delta}{\nu^{2}}\right)^{k}\cdot\frac{n^{n}}{(n-1)^{n-k}}\cdot\frac{1}{\prod_{j\in{\cal
S}}r_{j}}\cdot\prod_{i\in{\cal S}^{c}}\left(1-\frac{1}{r_{i}}\right).$

Using $r_{low}\leq r_{i}\leq r_{up}$ in the last two terms of (LABEL:r_first)
we then get that

$\displaystyle\mathbb{E}R_{\sigma_{1}}R_{\sigma_{2}}$ $\displaystyle\leq$
$\displaystyle\frac{\mu_{n}^{2}\alpha}{(n!)^{2}}\cdot\beta^{k}$ (2.8)

where $\alpha=\alpha_{up}$ and $\beta=\beta_{up}$ are as defined in (2.1).

Substituting (2.8) into (2.6) we get

$\mathbb{E}T_{n}^{2}\leq\frac{\mu_{n}^{2}\alpha}{(n!)^{2}}\sum_{\sigma_{1}}\sum_{k=0}^{n}\sum_{{\cal
S}:\\#{\cal S}=k}\sum_{\sigma_{2}:\sigma_{1}\cap\sigma_{2}={\cal
S}}\beta^{k}.$ (2.9)

Given $\sigma_{1}$ and ${\cal S}$ with $\\#{\cal S}=k,$ the number of
permutations $\sigma_{2}$ satisfying
$\sigma_{2}\cap\sigma_{1}={\cal S}$ is equal to the number of derangements of
$\\{1,2,\ldots,n-k\\},$ which in turn equals $(n-k)!\cdot b_{n-k},$ where
$b_{j}:=\sum_{l=0}^{j}\frac{(-1)^{l}}{l!}.$ From (2.9) we therefore get that

$\displaystyle\mathbb{E}T_{n}^{2}$ $\displaystyle\leq$
$\displaystyle\frac{\mu_{n}^{2}\alpha}{(n!)^{2}}\sum_{\sigma_{1}}\sum_{k=0}^{n}\sum_{{\cal
S}:\\#{\cal S}=k}(n-k)!\cdot b_{n-k}\cdot\beta^{k}$ (2.10) $\displaystyle=$
$\displaystyle\frac{\mu_{n}^{2}\alpha}{(n!)^{2}}\sum_{\sigma_{1}}\sum_{k=0}^{n}{n\choose
k}(n-k)!\cdot b_{n-k}\cdot\beta^{k}$ $\displaystyle=$
$\displaystyle\frac{\mu_{n}^{2}\alpha}{n!}\sum_{k=0}^{n}{n\choose
k}(n-k)!\cdot b_{n-k}\cdot\beta^{k}$ $\displaystyle=$
$\displaystyle\mu_{n}^{2}\alpha\sum_{k=0}^{n}\frac{\beta^{k}}{k!}\cdot
b_{n-k}.$

To estimate the final term in (2.10) we split the corresponding summation into
two terms depending on whether $1\leq k\leq n-\log{n}$ or $n-\log{n}\leq k\leq
n.$ Using $|b_{n-k}-e^{-1}|\leq\sum_{j\geq n-k+1}\frac{1}{j!}$ we have for
$0\leq k\leq n-\log{n}$ that

$|b_{n-k}-e^{-1}|\leq\sum_{j\geq\log{n}}\frac{1}{j!}\leq\frac{1}{n^{2}}$
(2.11)

for all $n$ large and so

$\sum_{k=0}^{n-\log{n}}\frac{\beta^{k}}{k!}\cdot|b_{n-k}-e^{-1}|\leq\frac{1}{n^{2}}\sum_{k=0}^{n-\log{n}}\frac{\beta^{k}}{k!}.$
(2.12)

For the other range $n-\log{n}\leq k\leq n,$ we use $|b_{n-k}|\leq e$ to get
that

$\sum_{k=n-\log{n}}^{n}\frac{\beta^{k}}{k!}\cdot|b_{n-k}-e^{-1}|\leq(e+e^{-1})\sum_{k=n-\log{n}}^{n}\frac{\beta^{k}}{k!}.$
(2.13)

Using $k!\geq k^{k}e^{-k}$ and $\beta\leq\frac{2n\delta}{r_{low}\nu^{2}}$ we
get that the final term in (2.13) is at most

$(e+e^{-1})\sum_{k=n-\log{n}}^{n}\left(\frac{2ne}{kr_{low}}\right)^{k}\leq(e+e^{-1})\sum_{k=n-\log{n}}^{n}\left(\frac{2ne\delta}{(n-\log{n})r_{low}\nu^{2}}\right)^{k}.$

Since $r_{low}\geq\frac{6\delta}{\nu^{2}}$ (see statement of Proposition 1) we
have $\frac{2ne\delta}{(n-\log{n})r_{low}\nu^{2}}\leq e^{-2C}$ for all $n$
large and some positive constant $C.$ Thus

$\sum_{k=n-\log{n}}^{n}\left(\frac{2ne\delta}{(n-\log{n})r_{low}\nu^{2}}\right)^{k}\leq\sum_{k\geq
n-\log{n}}e^{-2Ck}\leq e^{-Cn}$ (2.14)

for all $n$ large.

From (2.12), (2.13), (2.14) and the above discussion we then get that

$\displaystyle\sum_{k=0}^{n}\frac{\beta^{k}}{k!}\cdot b_{n-k}$
$\displaystyle\leq$
$\displaystyle\left(e^{-1}+\frac{1}{n^{2}}\right)\sum_{k=0}^{n}\frac{\beta^{k}}{k!}+(e+e^{-1})e^{-Cn}$
(2.15) $\displaystyle\leq$
$\displaystyle\left(e^{-1}+\frac{1}{n^{2}}\right)e^{\beta}+(e+e^{-1})e^{-Cn}$
$\displaystyle\leq$
$\displaystyle\left(e^{-1}+\frac{2}{n^{2}}\right)e^{\beta}$

for all $n$ large. Substituting (2.15) into (2.10) and using (2.6), we obtain
the upper bound in (2.3).

_Proof of the lower bound in ( 2.3)_: Using (LABEL:r_first) and arguing as in
the paragraph preceding (2.10) we have that

$\mathbb{E}T_{n}^{2}\geq\mu_{n}^{2}\cdot\alpha_{low}\sum_{k=0}^{n}\frac{\beta^{k}_{low}}{k!}\cdot
b_{n-k}$ (2.16)

and using (2.11) we get that the term
$\sum_{k=0}^{n}\frac{\beta^{k}_{low}}{k!}\cdot b_{n-k}$ is bounded below by

$\left(e^{-1}-\frac{1}{n^{2}}\right)\sum_{k=0}^{n-\log{n}}\frac{\beta^{k}_{low}}{k!}$
(2.17)

for all $n$ large. Using $r_{up}\geq r_{low}\geq\frac{6\delta}{\nu^{2}}$ and
arguing as in (2.13) and (2.14), we get that

$\sum_{k\geq n-\log{n}}\frac{\beta^{k}_{low}}{k!}\leq e^{-Dn}$ (2.18)

for some constant $D>0$ and all $n$ large. From (2.16), (2.17) and (2.18), we
then get that

$\displaystyle\mathbb{E}T_{n}^{2}$ $\displaystyle\geq$
$\displaystyle\mu_{n}^{2}\alpha_{low}\left(\left(e^{-1}-\frac{1}{n^{2}}\right)e^{\beta_{low}}-\left(e^{-1}-\frac{1}{n^{2}}\right)e^{-Dn}\right)$
$\displaystyle\geq$ $\displaystyle\mu_{n}^{2}\cdot\alpha_{low}\cdot
e^{\beta_{low}-1}\left(1-\frac{2e}{n^{2}}\right)$

for all $n$ large. Substituting this into (2.6), we get the lower bound in
(2.3).

##  3 Proof of Theorem 1

As before $T_{n}$ is the permanent of the constrained random matrix $Y.$ We
first find an upper bound for the variance of $T_{n}$ using Proposition 1. Let
$r_{up}=r$ and let $\alpha_{up}$ and $\beta_{up}$ be as in (2.1). Using
$\left(1-\frac{1}{n}\right)^{-n}=e(1+o(1))$ where $o(1)\longrightarrow 0$ as
$n\rightarrow\infty,$ we then have that

$e^{-1}\cdot\alpha_{up}=e^{-1}\cdot\left(1-\frac{1}{n}\right)^{-n}\cdot\left(1-\frac{1}{r}\right)^{n}=\left(1-\frac{1}{r}\right)^{n}(1+o(1))$

and $\beta_{up}=\frac{\delta r(n-1)}{\nu^{2}r_{low}(r-1)}\leq\frac{\delta
rn}{\nu^{2}r_{low}(r-1)}.$ Thus

$\alpha_{up}\cdot e^{\beta_{up}-1}=e^{-1}\alpha_{up}\cdot
e^{\beta_{up}}\leq\left(e^{\frac{\delta
r}{\nu^{2}r_{low}(r-1)}}\cdot\left(1-\frac{1}{r}\right)\right)^{n}\cdot(1+o(1)).$
(3.1)

For $r\geq 2$ we have that $\frac{r}{r-1}\leq 2$ and moreover, from (1.2) we
know that $\frac{\delta}{\nu^{2}r_{low}}\leq\frac{1}{4}$ for all $n$ large.
This implies that $\frac{\delta
r}{\nu^{2}r_{low}(r-1)}\leq\frac{2\delta}{\nu^{2}r_{low}}\leq\frac{1}{2}$ for
all $n$ large and so using standard $O$ notations we have that

$\displaystyle e^{\frac{\delta r}{\nu^{2}r_{low}(r-1)}}$ $\displaystyle=$
$\displaystyle 1+\frac{\delta
r}{\nu^{2}r_{low}(r-1)}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}$ (3.2)
$\displaystyle=$ $\displaystyle
1+\frac{\delta}{\nu^{2}r_{low}}+\frac{\delta}{\nu^{2}r_{low}(r-1)}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}$
$\displaystyle=$ $\displaystyle
1+\frac{\delta}{\nu^{2}r_{low}}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}$

since $r=r_{up}\geq r_{low}$ and
$\delta=\mathbb{E}Z^{2}\geq(\mathbb{E}Z)^{2}=\nu^{2}.$ Thus

$\displaystyle e^{\frac{\delta
r}{\nu^{2}r_{low}(r-1)}}\cdot\left(1-\frac{1}{r}\right)$ $\displaystyle=$
$\displaystyle\left(1+\frac{\delta}{\nu^{2}r_{low}}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}\right)\cdot\left(1-\frac{1}{r_{up}}\right)$
(3.3) $\displaystyle=$ $\displaystyle
1+\frac{\delta}{\nu^{2}r_{low}}-\frac{1}{r_{up}}-\frac{\delta}{\nu^{2}r_{low}r_{up}}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}$
$\displaystyle=$ $\displaystyle
1+\frac{\delta}{\nu^{2}r_{low}}-\frac{1}{r_{up}}+O\left(\frac{\delta}{\nu^{2}r_{low}}\right)^{2}$
$\displaystyle=$ $\displaystyle 1+o\left(\frac{1}{n}\right)$

where as before, the third expression in (3.3) follows from the fact that
$r_{up}\geq r_{low}$ and $\delta\geq\nu^{2}$ and the final expression in (3.3)
follows from the conditions in (1.2).

Plugging (3.3) into (3.1) we get that

$\alpha_{up}\cdot
e^{\beta_{up}-1}\leq\left(1+o\left(\frac{1}{n}\right)\right)^{n}\cdot(1+o(1))=1+o(1).$

From (2.3) we therefore get that $\frac{\mathbb{E}T_{n}^{2}}{\mu_{n}^{2}}\leq
1+o(1)$ or equivalently that
$var(T_{n})\leq o(\mu_{n}^{2}).$ This obtains (1.3) and therefore completes
the proof of Theorem 1.

## Acknowledgement

I thank Professors Rahul Roy, C. R. Subramanian and the referee for crucial
comments that led to an improvement of the paper. I also thank IMSc for my
fellowships.

## References

  * [1] Egorychev, G. P. (1981).  The Solution of van der Waerden’s Problem for Permanents.  Advances in Mathematics, 42, pp. 299–305.
  * [2] Falikman, D. I. (1981).  Proof of the van der Waerden’s Conjecture on the Permanent of a Doubly Stochastic Matrix.  Mathematical notes of the Academy of Sciences of the USSR, 29, pp. 475–479.
  * [3] Kwan, M. and Sauermann, L. (2020).  On the Permanent of a Random Symmetric Matrix.  arXiv preprint, Link: _https://arxiv.org/abs/2010.08922_
  * [4] Rempala, G. and Gupta, A. (1999).  Some Extensions of Girko’s Limit Theorems for Permanents of Random Matrices.  Random Operators and Stochastic Equations, 8, pp. 305–318.
  * [5] Tao, T. and Vu, V. (2009).  On the Permanent of Random Bernoulli Matrices.  Advances in Mathematics, 220, pp. 657–669.
  * [6] Vu, V. (2020).  Recent Progress in Combinatorial Random Matrix Theory.  arXiv preprint, Link: _https://arxiv.org/pdf/2005.02797.pdf_.