# Minimal sized generating sets of permutation groups

Derek F. Holt and Gareth Tracey

(5th April 2024)

###### Abstract

We present a randomised variant of an algorithm of Lucchini and Thakkar [10]
for finding a smallest sized generating set in a finite groupd, which has
polynomial time expected running time in finite permutation groups.

##  1 Introduction

Lucchini and Thakkar prove in [10] that a generating set of smallest possible
size in a finite group $G$ can be found in time polynomial in the group order
$|G|$ and, in the case when $G$ is given as a subgroup of the symmetric group
$\mathrm{Sym}(n)$, this can be done using $n^{2}\lambda(G)^{13/5}$ _generating
tests_ , where $\lambda(G)$ is the maximum of the orders of the nonabelian
composition factors of $G$, with $\lambda(G)=1$ for solvable $G$. In this note
we prove the following result.

######  Theorem 1.

There is a randomised variant of the algorithm described in [10] for subgroups
of $\mathrm{Sym}(n)$ with expected running time polynomial in $n$. The
expected number of generating tests is $O(n^{2}\log n)$, and
$O(n^{2}/\sqrt{\log n})$ if $G$ is transitive.

We have implemented this variant in Magma [2] as a function
SmallestGeneratingSet, and it runs fast in practice.

A generating test consists of testing whether a given subset $S$ of $G$ (which
in this paper will always have size bounded by the permutation degree $n$),
generates the whole of $G$. The _Schreier-Sims_ algorithm, which is easily
seen to run in polynomial time, can be used for this purpose. There is an
extensive discussion and detailed analysis of the complexity of variants of
this fundamental algorithm in [13, Chapters 4 and 8].

For the proof of the theoretical result on expected running time, we assume
that a deterministic version of the Schreier-Sims algorithm, which we denote
by SS, is used. There is however a randomised Monte-Carlo version of the
Schreier-Sims algorithm described in [13, page 64] that often runs much faster
in practice, but has a small probability of returning an incorrect negative
answer; that is, it falsely reports that $\langle S\rangle\neq G$. We shall
call this version RSS. In some (but not all) of the situations in which we
apply generating tests in the algorithm presented in this paper, we repeatedly
choose random elements in a search for certain generators of $G$ that we know
exist, and in these situations we can safely use RSS, because an incorrect
negative answer would merely prolong the search. Experiments with
implementations indicate that the increased speed that we gain from using RSS
rather than SSmore than compensates for the time wasted by the occasional
wrong answer. In our descriptions below we shall indicate in brackets when we
can use RSS in an implementation, and when we must use SS.

##  2 Definitions and preliminary results

We shall assume that the reader is familiar with the results and proofs in
[10] and with the proof of Proposition 6 of that paper in particular.

The notion of $G$-equivalent chief factors of $G$ is described in [10, Section
4] and, for a nonabelian chief factor $A$, the total number of chief factors
of $G$ that are $G$-equivalent to $A$ is denoted by $\delta_{G}(A)$. There is
a normal subgroup $R_{G}(A)$ of $G$ such that the socle of $G/R_{G}(A)$ is the
direct product of these $\delta_{G}(A)$ chief factors.

For a finite group $X$, let $P(X)$ be the smallest degree of a faithful
permutation representation of $X$; that is, the smallest $n$ for which $G$
embeds in $\mathrm{Sym}(n)$. It is proved in [8] that, if $X$ is a finite
group and $X/N$ has no abelian normal subgroups, then $P(X/N)\leq P(X)$. In
[5, Theorem 3.1] it is shown that, for a direct product
$N=S_{1}\times\cdots\times S_{r}$ of finite nonabelian simple groups, we have
$P(N)=\sum_{i=1}^{r}P(S_{i})$. So we have the following result.

######  Lemma 2.

If $N\cong S^{b}$ is a nonabelian chief factor of the group $G$ with $S$
simple, and $\delta_{G}(N)$ is as defined above, then $P(G)\geq
b\delta_{G}(N)P(S)$.

We will also need a bound on the number of nonabelian chief factors in a
transitive permutation group.

######  Lemma 3.

Let $G\leq S_{n}$ be transitive. Then $G$ has at most $\log{n}$ nonabelian
chief factors.

###### Proof.

When $G$ is primitive, the claim is [12, Theorem 2.10]. So assume that $G$ is
imprimitive, and let $\Sigma$ be a block system with $|\Sigma|$ minimal, so
$G^{\Sigma}$ is primitive. Let $\Delta\in\Sigma$, and write $|\Delta|=r$ and
$|\Sigma|=s$, so that $n=rs$. Then by [14, Lemma 5.8], $G$ embeds as a
subgroup of $\mathrm{Stab}_{G}(\Delta)^{\Delta}\wr G^{\Sigma}$ in such a way
so that if $E:=X/Y$ is a nonabelian chief factor of
$R:=\mathrm{Stab}_{G}(\Delta)^{\Delta}$, then $(G\cap X^{s})/(G\cap Y^{s})$ is
either trivial or a nonabelian chief factor of $R$. It follows that, if $b(G)$
denotes the number of nonabelian chief factors of $G$, then $b(G)\leq
b(R)+b(G^{\Sigma})$. The result then follows by induction. ∎

##  3 The algorithm

The input to the algorithm is a permutation group $G\leq S_{n}$ defined by a
generating set. We start by computing a chief series
$1=N_{0}<N_{1}<...<N_{u}=G$ of $G$. This can be done in polynomial time; in
fact in _nearly linear_ time [13, Section 6.2.7]. We have $u\leq n-1$ by [11,
Theorem 10.0.5].

Now we compute a generating set of the top factor $G/N_{u-1}$, which is
simple. If it is abelian then it is cyclic and any element of $G\setminus
N_{u-1}$ generates $G$ modulo $N_{u-1}$. Otherwise it can be generated by two
elements, and by [11, Theorem 3.1.7] a random pair of elements of $G$ has a
probability of at least $53/90$ of generating $G$ modulo $N_{u-1}$. So we find
such a pair by repeatedly choosing random pairs and doing generating tests,
for which we can use RSS.

In the main “lifting” step of the algorithm, we have found elements
$g_{1},\ldots,g_{d}$ that generate $G$ modulo $N_{k}$ for some $k$ with $1\leq
k\leq u-1$, and we want to lift them to a generating set of $G$ modulo
$N_{k-1}$. We know from Proposition 9 and the proof of Corollary 13 of [10]
that either (Lifting Case 1) there exist $n_{1},\ldots,n_{d}$ in $N_{k}$ such
that $g_{1}n_{1},\ldots,g_{d}n_{d}$ generate $G$ modulo $N_{k-1}$, or (Lifting
Case 2) there exist $n_{1},\ldots,n_{d+1}$ in $N_{k}$ such that
$g_{1}n_{1},\ldots,g_{d}n_{d},n_{d+1}$ generate $G$ modulo $N_{k-1}$. Note in
particular that we always have $d\leq u+1-k$ and so $d\leq n-1$.

Let $N:=N_{k}/N_{k-1}$. If $N$ is abelian and hence isomorphic to $C_{p}^{l}$
for some prime $p$ and $l>0$, then we proceed as in [10, Remark 11]. That is,
we find a generating set $e_{1},\ldots,e_{l}$ of $N_{k}$ module $N_{k-1}$ and,
for each $i=1,\ldots,d$ and $j=1,\ldots,l$, we test whether
$g_{1},\ldots,g_{i-1},g_{i}e_{j},g_{i+1},\ldots,g_{d}$ generate $G$ modulo
$N_{k-1}$. If so then we are in Lifting Case 1 and we are done. If not, then
by [10, Proposition 9] we are in Lifting Case 2, and furthermore $G$ is
generated modulo $N_{k-1}$ by $g_{1},\ldots,g_{d},x$ for any $x\in
N_{k}\setminus N_{k-1}$.

Otherwise $N$ is nonabelian, so $N\cong S^{b}$ for some nonabelian simple
group $S$ and $b\geq 1$. Let $\delta:=\delta_{G/N_{k-1}}(N)$. As is observed
in [10] there is no known efficient method of calculating $\delta$, so in our
algorithm we use the possibly larger number
$\delta^{\prime}:=\delta^{\prime}_{G/N_{k-1}}(N)$, which we define to be the
number of chief factors of $G/N_{k-1}$ having order $|N|$. This could
potentially increase the running time, but will not affect the correctness of
the algorithm.

If $d=1$ then, since $N$ is not cyclic, we know that we must be in Lifting
Case 2, so we keep testing random pairs $n_{1},n_{2}$ (for which we may RSS in
the implementation) until $g_{1}n_{1}$ and $n_{2}$ generate $G$ modulo
$N_{k-1}$ and, as we shall demonstrate in the next section, the expected
number of pairs that we need test is at most $18n/53$.

So we assume henceforth that $d\geq 2$. There are two cases. Define
$t:=\lceil{8/5+\log_{|N|}(\delta)}\rceil$ as in [10, Lemma 12] and
$t^{\prime}:=\lceil{8/5+\log_{|N|}(\delta^{\prime})}\rceil$. Then, if
$t^{\prime}\leq d$ (and hence also $t\leq d$), we are in Lifting Case 1 by
[10, Lemma 12], and we repeatedly test sequences
$g_{1}e_{1},\ldots,g_{d}e_{d}$ for random $e_{1},\dots,e_{d}\in N$ (using RSS
in the implementation) until we find the required generators. As we shall show
in the next section, the expected number of sequences that we need to choose
is at most $18n/53$.

Otherwise (if $t^{\prime}>d$), we do not know whether we are in Lifting Case 1
or 2. In that case, we choose (at most) $90d\delta^{\prime}\log|N|/53$ random
sequences $n_{1},\ldots,n_{d}$, and test them all (using RSS in the
implementation). If we find generators $g_{1}n_{1},\ldots,g_{d}n_{d}$ of $G$
modulo $N_{k-1}$, then we are done.

If not, then we do an exhaustive test of all $|N|^{d}$ such sequences using
SS. If we still fail to find them then we know that we are in Lifting Case 2,
so we keep testing random sequences $n_{1},\ldots,n_{d+1}$ (using RSS in the
implementation) until we find one such that
$g_{1}n_{1},\ldots,g_{d}n_{d},n_{d+1}$ generate $G$ modulo $N_{k-1}$. The
expected number of trials is again at most $18n/53$.

##  4 Proof of polynomial running time

We shall prove that the expected number of generating tests applied in the
algorithm that we have described is $O(n^{2}\log{n})$, and $O(n^{2}/\sqrt{\log
n})$ when $G$ is transitive, thereby proving Theorem 1. Many of the processes
involved are easily seen to be polynomial-time, usually linear in $n$, and
will not require comment. When we refer to the _cost_ of part of the
algorithm, we shall mean the (expected) number of generating tests required.

When $N$ is abelian, the elements $e_{1},\ldots,e_{l}$ from the previous
section are easily found using random choice and up to $l$ generating tests.
We then apply up to $dl$ generating tests (for which we must use SS) to see if
we are in Lifting Case 1. It is proved in [6, Theorem 1.2] that a composition
series of a subgroup of $\mathrm{Sym}(n)$ with $s$ orbits has length at most
$4(n-s)/3$, so the sum of the values of $l$ arising from all elementary
abelian chief factors is at most $4(n-1)/3$ and, since $d\leq n-1$, the total
cost of the algorithm for all elementary abelian composition factors is at
most $4(n-1)^{2}/3$ generating tests. Moreover, when $G$ is transitive we have
$d=O(n/\sqrt{\log n})$ by [9], and so the number of generating tests required
is $O(n^{2}/\sqrt{\log n})$.

So assume now that $N$ is nonabelian. We first justify our claims above that
the expected number of choices of random sequences of elements of $N$ in order
to find the required generating set of $G$ modulo $N_{k-1}$ when it exists is
at most $18n/53$. For this we need to look more carefully at the proof of [10,
Proposition 6]. In the notation of that proof, there is a group $\Gamma$
acting semiregularly on a set $\Omega$ that consists of certain sequences
$e_{1},e_{2},\ldots,e_{d}$ of elements of $N$ where, by [10, Lemma 5], we have
$|\Omega|/|N|^{d}\geq 53/90$. Furthermore, the property of whether or not
$g_{1}e_{1},\ldots,g_{d}e_{d}$ generate $N_{k}$ modulo $N_{k-1}$ is the same
for all elements in an orbit of this action, and there are at most $\delta-1$
orbits for which this sequence fails to generate $G$ modulo $N_{k-1}$.

Now if we are in Lifting Case 1, then there must be at least one orbit for
which the sequence does generate $G$ modulo $N_{k-1}$. In that case, it
follows that the probability that $G$ is generated modulo $N_{k-1}$ by
$g_{1}e_{1},\ldots,g_{d}e_{d}$ for a random sequence
$e_{1},e_{2},\ldots,e_{d}$ of elements of $N$ is at least $53/(90\delta)$.
Since $\delta\leq n/5$ by Lemma 2, this probability is at least $53/(18n)$. In
Lifting Case 2, by the same reasoning, the probability that
$g_{1}e_{1},\ldots,g_{d}e_{d},e_{d+1}$ generate $G$ modulo $N_{k-1}$ for a
random sequence $e_{1},e_{2},\ldots,e_{d+1}$ of elements of $N$ is also at
least $53/(18n)$. So in both of these situations, if we are looking for
generators of $G$ modulo $N_{k-1}$ by repeatedly choosing random sequences of
elements of $N$, then the expected number of choices necessary is at most
$18n/53$.

In the case $t^{\prime}\leq d$ the expected cost is at most $18n/53$. But when
$t^{\prime}>d$, we carry out up to $90d\delta^{\prime}\log|N|/53$ generating
tests. We have $d<13/5+\log_{|N|}(\delta^{\prime})$ in this case, and so

$90d\delta^{\prime}\log|N|/53<234\delta^{\prime}\log|N|/53+90\delta^{\prime}\log{\delta^{\prime}}/53.$

Since the sum of $\log|N|$ over all chief factors of $G$ is at most
$\log(n!)=O(n\log(n))$, the sum of the first of these terms over all chief
factors is $O(n^{2}\log n)$ and, since there are at most $n$ such chief
factors and $\delta^{\prime}<n$, the same applies to the sum of the second
terms. So the total expected cost of tests of this type over all chief factors
is $O(n^{2}\log n)$. If $G$ is transitive, then $G$ has at most $\log{n}$
nonabelian chief factors by Lemma 3, so the total expected cost we get above
is $O(n(\log{n})^{2})$.

It remains to estimate the expected cost of the exhaustive tests over all
sequences $e_{1},\ldots,e_{d}\in N^{d}$ that we need to carry out. We do the
exhaustive test in two situations:

1\. We are in Lifting Case 1, but we fail to find the required elements
$n_{1},\ldots,n_{d}$. Since we perform $90d\delta^{\prime}\log|N|/53$ tests of
random sequences, each of which has a probability of at least $53/(18n)$ of
succeeding but fails, and the cost of the exhaustive test itself is
$|N|^{d}=e^{d\log|N|}$, the expected cost for each chief factor is at most

$(1-53/(90\delta))^{90d\delta^{\prime}\log|N|/53}e^{d\log|N|}\leq
e^{-d\log|N|},$

giving a total expected cost over all chief factors of at most $ne^{-2\log
60}$.

2\. We are in Lifting Case 2. From the proofs of Lemma 12 and Proposition 6 of
[10], we see that their conclusions must hold provided that $|N|^{t}\geq
90b|N||\mathrm{Out}(S)|\delta/53$, and so we can replace the value of $t$ in
these results by $\lceil 1+\log_{|N|}(f(N)\delta)\rceil$, where
$f(N):=90b|\mathrm{Out}(S)|/53$. By Lemma 2, we have $n\geq b\delta P(S)$ and
we have $|\mathrm{Out}(S)|\leq 2P(S)/3$ by Lemma 4 below, so $f(N)\delta\leq
60n/53$.

So, if we in Lifting Case 2, then $d<\lceil 1+\log_{|N|}(f(N)\delta)\rceil$
and, since $d$ is an integer, we have

$d<1+\log_{|N|}(f(N)\delta)\leq 1+\log_{|N|}(60n/53)$

so $|N|^{d-1}<60n/53$, and the cost $|N|^{d}$ of the test is less than
$(60n/53)^{d/(d-1)}<(60/53)^{2}n^{2}$ (maximised when $d=2$).

Let $k$ be the number of nonabelian chief factors for which Lifting Case 2
occurs. Then on the final occurrence we must have $d\geq 2+k-1$, and so
$2+k-1<1+\log_{|N|}(60n/53)$ and hence $k+1<\log n$. Then the total cost over
all such factors is at most

$\sum_{d=2}^{k+1}n^{d/(d-1)}<\sum_{d=2}^{\log n}n^{d/(d-1)}<n^{2}+n^{3/2}\log
n=O(n^{2}).$

Suppose now that $G$ is transitive. From the estimate above, we see that to
verify the claimed bound of $O(n^{2}/\sqrt{\log n})$ generating tests, it
suffices to do this for a nonabelian factor with $d=2$. In that case we have
$|N|=|N|^{d-1}\leq 60n/53$. From Lemma 2.3.8 (b) of [6] we have
$|\mathrm{Out}(S)|\leq\log_{2}|S|$ for all nonabelian finite simple groups
$S$, and so

$\displaystyle f(N)=\ $ $\displaystyle
90b|\mathrm{Out}(S)|/53\leq(90/53)b\log_{2}|S|=(90/53)\log_{2}|N|\leq$
$\displaystyle(90/53)\log_{2}(60n/53).$

But $2=d<1+\log_{|N|}(f(N)\delta)$ implies $|N|<f(N)\delta$, and
$\delta\leq\log n$ by Lemma 3, so we get $|N|=O((\log n)^{2})$ and cost is
$O(|N|^{2})=O((\log n)^{4})$ generating tests, which completes the proof.

##  5 Implementation and performance

Our Magma implementation works for finite groups of matrices (type GrpMat and
solvable groups defined by a power-conjugate presentation (type GrpPC) as well
as permutation groups (type GrpPerm) but since this article relates to
permutation groups, we shall use these only in our sample run-times below. We
note however that for soluble like $3^{k}:2$, it appears to be faster to work
with a group of type GrpPC.

In the implementation we deal separately first with cyclic groups, then with
$p$-groups, and then nilpotent groups. In these cases we can easily calculate
$d(G)$ and find generators as inverse images of a smallest sized generating
set of the abelian group $G/\Phi(G)$.

Otherwise, putting $k:=\min(2,d(G/[G,G])$, which we can easily calculate, we
know that $d(G)\geq k$, and we start by choosing a small number of random
subsets of $G$ of size $k$, and test (using RSS) whether any of them generate
$G$. This enables us to find a smallest generating set very quickly in a large
number of groups.

If this fails, then we apply the algorithm that we have described in Section
3. For the exhaustive test in the case when the chief factor $N$ is non-
abelian, we provide the option of computing the conjugation action of $N$ on
the set of $|N|^{d}$ candidates for generating sets, and only carrying out
generating tests on orbits representatives of this action. This results in
significant improvements in run-times.

Immediately after processing each chief factor we check whether the current
generators generate the whole group, and we can stop early if they do. This
results in significant variation from run to run of the times for the examples
in the table below, depending on whether or not we are lucky, and find
generators early.

We found that, for groups with $d(G)>2$ with both abelian and nonabelian chief
factors, if there is more than one chief series, then a series in which the
abelian factors are nearer the top of the group is likely to result in faster
run-times. This is because processing the abelian factors is generally faster,
and if the current generating set becomes larger early in the complete
calculation, then the exhaustive tests with the non-abelian factors are less
likely to be necessary.

Unfortunately, the Magma ChiefSeries intrinsic does the opposite, and
constructs series with the chief factors in the soluble radical at the bottom.
It is difficult to write an alternative version using existing Magma
functionality, but we did manage a compromise in which the factors in the non-
abelian socle of the group are at the bottom. This results in much faster
times in the final three examples in the table below, where the two reported
times are using the default ChiefSeries function, and our modified version.

It is not easy to find examples that take more than a few seconds to run. For
example the transitive groups of degree 48 almost all take less than $0.1$
seconds.

The examples in the table below have been chosen because they are more
difficult, and some of them involve applying the exhaustive test for a non-
abelian chief factor. The split extension in the table are crowns: that is,
the complement of the specified normal subgroup acts in the same way on each
of the direct factors of the normal subgroup. The two irreducible
$4$-dimensional ${\mathbb{F}}_{2}A_{5}$-modules are denoted by $2^{4}_{a}$ and
$2^{4}_{b}$, where the former is from the ntural module for $\mathrm{SL}(2,4)$
(which is isomorpic to $A_{5}$) and the latter is the deleted permutation
module.

Group | Degree | $d(G)$ | Time
---|---|---|---
$A_{5}^{19}$ | 95 | 2 | 0.6
$A_{5}^{20}$ | 100 | 3 | 3
$A_{5}^{100}$ | 500 | 3 | 252
$L_{3}(2)^{57}$ | 399 | 2 | 75
$L_{3}(2)^{58}$ | 506 | 3 | 236
$L_{3}(2)^{100}$ | 700 | 3 | 350
$L_{3}(2)^{100}:2$ | 800 | 3 | 513
$A_{6}^{53}$ | 318 | 2 | 68
$A_{6}^{54}$ | 324 | 3 | 804
$A_{6}^{100}$ | 600 | 3 | 1295
$3^{30}:2$ | 90 | 31 | 0.51
$3^{50}:2$ | 150 | 51 | 8
$3^{100}:2$ | 300 | 101 | 297
$(2^{4}_{a})^{20}:A_{5}$ | 320 | 12 | 2.8
$(2^{4}_{a})^{50}:A_{5}$ | 800 | 27 | 172
$(2^{4}_{b})^{50}:A_{5}$ | 500 | 14 | 48
$((2^{4}_{b})^{5}:A_{5})\times A_{5}^{50}$ | 300 | 3 | 27, 16
$((2^{4}_{b})^{5}:A_{5})\times(L_{3}(2)^{58})$ | 456 | 3 | 259, 77
$(3^{10}:2)\times(L_{3}(2)^{58})$ | 436 | 11 | 224, 34

##  6 Appendix

Lemma 2.3.8 (a) of [6] asserts that $P(S)>2|\mathrm{Out}(S)|$ for all finite
simple groups $S$ except for $S=A_{6}$, when $P(S)=6$ and
$|\mathrm{Out}(S)|=4$. For the proof, the authors cite [1, Lemma 2.7] which
claims that the result is true “by inspection”. But in fact
$S=\mathrm{PSL}(3,4)$ with $P(S)=21$ and $|\mathrm{Out}(S)|=12$ is another
exception to the general result, so we decided to re-prove it here. Since
$P(S)$ and $|\mathrm{Out}(S)|$ are known for all such $S$, the proof is indeed
a straightforward verification but, in view of the error, it seems preferable
to write it down.

In fact the only place in the paper [6] where Lemma 2.3.8 (a) of [6] is used
is in Case 2 (simple diagonal) of their proof of Theorem 1.6 (page 568) where
what they actually need is $|\mathrm{Out}(S)|<|S|/16$ for all nonabelian
simple groups $S$. This follows from Lemma 2.3.8 (b) of [6], which (correctly)
states that $|\mathrm{Out}(S)|\leq\log_{2}|S|$, and $\log_{2}|S|<|S|/16$ for
$|S|>108$. So to complete the proof of [6, Theorem 6.1] we need only observe
that $|\mathrm{Out}(A_{5})|=2<|A_{5}|/16=3.75$.

######  Lemma 4.

Let $S$ be a finite simple group. Then $P(S)>3|\mathrm{Out}(S)|$ except when
$(S,P(S),|\mathrm{Out}(S)|)=(A_{6},6,4)$, $(\mathrm{PSL}(3,4),21,12)$,
$(A_{5},5,2)$, or $(\mathrm{PSL}(2,8),9,3)$.

###### Proof.

It is well-known that $P(A_{n})=n$ and $|\mathrm{Out}(A_{n})|=2$ when $n\neq
6$, and $|\mathrm{Out}(A_{6})|=4$. For $S$ sporadic, we have
$|\mathrm{Out}(S)|\leq 2$ and $P(S)\geq 11$. So the result holds in these
cases.

For $S$ of Lie type, we refer to [4, Page xvi] for $|\mathrm{Out}(S)|$, and
[7, Table 4] (which is derived from earlier results of Mazurov, Vasilyev and
others) for $P(S)$. Let $q=p^{e}$ with $p$ prime in the following list of
cases.

  1. 1.

$S=\mathrm{PSL}(2,q)$ with $q\geq 7$ and $q\neq 9$:
$|\mathrm{Out}(S)|=\gcd(2,q-1)e$ and $P(S)=q+1$ (or $q$ when $q=7$ or $11$).
We check the result directly for $q\leq 16$ and note that $P(S)=p^{e}+1>6e\geq
3|\mathrm{Out}(S)|$ for $q>16$.

  2. 2.

$S=\mathrm{PSL}(d,q)$ with $d\geq 3$ and $(d,q)\neq(4,2)$
($\mathrm{PSL}(4,2)\cong A_{8}$): $|\mathrm{Out}(S)|=2e\gcd(d,q-1)$ and
$P(S)=(q^{d}-1)/(q-1)$. For $d\geq 4$ and $(d,q)\neq(4,2)$, we have
$P(S)>q^{d-1}>6eq>3|\mathrm{Out}(S)|$ and the same holds for $d=3$ (i.e.
$q>6e$) when $q\not\in\\{2,3,4,5,8,9,16\\}$. For those cases we can check that
$P(S)=q^{2}+q+1>3|\mathrm{Out}(S)|=6e\gcd(d,q-1)$ except for
$S=\mathrm{PSL}(3,4)$.

  3. 3.

$S=\mathrm{PSU}(d,q)$ with $d\geq 3$ and $(d,q)\neq(3,2)$:
$|\mathrm{Out}(S)|=2e\gcd(d,q+1)$. When $d=3$ (and $q>2$) we have
$P(S)=q^{3}+1>18e\geq 3|\mathrm{Out}(S)|$ when $q\neq 5$, and
$P(S)=50>18=3|\mathrm{Out}(S)|$ when $q=5$. There are various cases for $P(S)$
when $d\geq 4$, but it is easily checked that $P(S)>q^{4}>6(q+1)e\geq
3|\mathrm{Out}(S)|$ when $q>2$, and $P(S)\geq 27>3|\mathrm{Out}(S)|$ when
$q=2$.

  4. 4.

$S=\mathrm{PSp}(d,q)$ with $d\geq 4$ and $(d,q)\neq(4,2)$:
$|\mathrm{Out}(S)|=\gcd(d,q-1)e$. In all cases we have $P(S)\geq
q^{3}>3qe>3|\mathrm{Out}(S)|$.

  5. 5.

$S=\mathrm{P}\Omega^{\epsilon}(d,q)$ with $d\geq 7$ and $d(q-1)$ even (note
that $\mathrm{P}\Omega^{\circ}(d,q)\cong\mathrm{PSp}(d-1,q)$ when $d$ is odd
and $q$ is even): in all cases $|\mathrm{Out}(S)|\leq 24e$ with
$|\mathrm{Out}(S)|\leq 6$ when $q=2$, and $P(S)>q^{5}>3|\mathrm{Out}(S)|$.

  6. 6.

$S$ is an exceptional group of Lie type. If $G={{}^{2}}B_{2}(q)$ with $p=2$
and $e\geq 3$ odd, then $P(S)=q^{2}+1>3e=3|\mathrm{Out}(S)|$. If
$G={{}^{2}}G_{2}(q)$ with $p=3$ and $e\geq 3$ odd, then
$P(S)=q^{3}+1>3e=3|\mathrm{Out}(S)|$. If $G=G_{2}(4)$ then
$P(S)=416>6=3|\mathrm{Out}(S)|$. In all other cases we have
$|\mathrm{Out}(S)|\leq 6e$, and $P(S)>q^{5}>18e\geq 3|\mathrm{Out}(S)|$.

∎

## References

  * [1] M. Aschbacher and R. M. Guralnick. On abelian quotients of primitive groups. Proc. Amer. Math. Soc. 107 (1989), 89–95.
  * [2] Wieb Bosma, John Cannon, and Catherine Playoust.  The Magma algebra system I: The user language.  J. Symbolic Comput. 24 (1997), 235–265.
  * [3] P. J. Cameron. _Permutation Groups_. London Mathematical Society Student Texts, 45. CUP, Cambridge, 1999.
  * [4] J. H. Conway, R. T. Curtis, S. P. Norton, R. A. Parker, and R. A. Wilson.  An $\mathbb{ATLAS}$ of Finite Groups. Clarendon Press, Oxford, 1985; reprinted with corrections 2003.
  * [5] D. Easdown and C. E. Praeger, On minimal faithful permutation representations of finite groups. Bull. Austral. Math. Soc. 38 (1988), 207–220.
  * [6] S. P. Glasby, C. E. Praeger, K. Rosa, and G. Verret. Bounding the composition length of primitive permutation groups and completely reducible linear groups. J. Lond. Math. Soc. (2) 98 no. 3 (2018), 557–572.
  * [7] S. Guest, J. Morris, C. E. Praeger, and P. Spiga P. On the maximum orders of elements of finite almost simple groups and primitive permutation groups. Trans. Amer. Math. Soc. 367 (2015), 7665–7694.
  * [8] L. G. Kovács and C. E. Praeger. On minimal faithful permutation representations of finite groups. Bull. Aust. Math. Soc. 62 no. 2 (2000), 311–317.
  * [9] A. Lucchini, F. Menegazzo, and M. Morigi. Asymptotic results for transitive permutation groups. Bull. London. Math. Soc. 32 (2000) 191–195.
  * [10] Andrea Lucchini and Dhara Thakkar. The Minimum Generating Set Problem, arXiv:2306.07633, 2023.
  * [11] Nina Menezes. Random generation and chief length of finite groups. PhD. Thesis, `http://research-repository.st-andrews.ac.uk/handle/10023/3578`
  * [12] László Pyber. Asymptotic results for permutation groups. Groups and Computation DIMACS Ser. Discrete Math. Theoret. Computer Sci. 11 (ed. Finkelstein, L. and Kantor, W.M., Amer. Math. Soc., Providence, 1993) 197-219.
  * [13] Ákos Seress. Permutation Group Algorithms. Cambridge Tracts in Mathematics 152. Cambridge University Press, Cambridge, 2003\.
  * [14] Gareth Tracey. Minimal generation of transitive permutation groups. _J. Algebra_ 509 (2018) 40–100.


# Baseline Drift Tolerant Signal Encoding for ECG Classification with Deep
Learning

Robert O’Shea  King’s College London
Email<EMAIL_ADDRESS>Prabodh Katti  King’s College London
Email<EMAIL_ADDRESS>Bipin Rajendran  King’s College London
Email<EMAIL_ADDRESS>
###### Abstract

Common artefacts such as baseline drift, rescaling, and noise critically limit
the performance of machine learning-based automated ECG analysis and
interpretation. This study proposes Derived Peak (DP) encoding, a non-
parametric method that generates signed spikes corresponding to zero crossings
of the signal’s first and second-order time derivatives. Notably, DP encoding
is invariant to shift and scaling artefacts, and its implementation is further
simplified by the absence of user-defined parameters. DP encoding was used to
encode the 12-lead ECG data from the PTB-XL dataset ($n$=18,869 participants)
and was fed to 1D-ResNet-18 models trained to identify myocardial infarction,
conductive deficits and ST-segment abnormalities. Robustness to artefacts was
assessed by corrupting ECG data with sinusoidal baseline drift, shift,
rescaling and noise, before encoding. The addition of these artefacts resulted
in a significant drop in accuracy for seven other methods from prior art,
while DP encoding maintained a baseline AUC of 0.88 under drift, shift and
rescaling. DP achieved superior performance to unencoded inputs in the
presence of shift (AUC under $1\,$mV shift: 0.91 vs 0.62), and rescaling
artefacts (AUC 0.91 vs 0.79). Thus, DP encoding is a simple method by which
robustness to common ECG artefacts may be improved for automated ECG analysis
and interpretation.

##  I Introduction

Deep neural networks have demonstrated promising results for automated
electrocardiograph (ECG) classification and monitoring [1]. However, intensive
energy requirements have been identified as a potential barrier to clinical
implementation within the resource constraints of implanted and wearable
devices [1]. Event-driven ECG processing limits communication to salient
events such as waveform complexes, thereby reducing energy requirements [2, 3,
4, 5, 6]. “Spike”-based encoding optimises efficiency and latency further
through ECG representation as sparse sequences of “all-or-none” spike signals.
Spiking neural network models (SNNs) employ spike encoding throughout,
facilitating accurate ECG classification on microjoule-range energy
budgets[6]. Thus, SNNs present an opportunity to both improve battery-life and
response latency in implanted cardiac defibrillators and pacemakers and
wearable monitors [2, 3, 4, 5]. However, ECG signals are subject to
significant distributional variability under normal operating conditions [7,
8], presenting a major reliability obstacle for spike-based cardiovascular
monitoring and diagnosis. Physiological factors such as respiration and limb
movement cause baseline drift, a low-frequency oscillatory artefact[9]. In
contrast, diagnostic ECG features are subtle, as a $0.2\,$mV ST segment
elevation may signify myocardial infarction [10]. This problem is considered
one of the most significant challenges for robust ECG classification [8, 7, 1,
11]. Event-driven ECG encoding techniques typically convert analogue input to
spikes according to fine-tuned threshold criteria [12, 6, 4]. Consequently,
these methods are particularly vulnerable to commonly encountered artefacts
such as baseline drift [11, 13]. To improve the robustness of spike-based ECG
encoding to these factors, we propose a parameter-free method based on zero-
crossings of first and second-order time derivatives.

###  I-A Related Work

“Rate” or “count” encoding [14, 3], encodes inputs over multiple channels,
determining the proportion of active channels by input discretisation. “Field
response” encoding activates a single channel corresponding to the discretised
input, improving efficiency by limiting instantaneous spike output to one-per-
timestep in the channel dimension. “Level crossing” encoding [12, 6] generates
spikes when discretised input values change, improving efficiency by omitting
sequential spikes at similar levels [12, 6]. “BitString” encoding or “virtual
neuron” [15] encoding represents the discretised input as a signed $m$-bit
integer over the channel dimension. BitString encoding encodes inputs as
binary integer strings. As discretisation-based methods require pre-specified
thresholds, parameterisation must balance adequate signal representation with
energy efficiency and channel depth. All discretisation methods are
theoretically vulnerable to baseline drift, which effectively shifts spiking
thresholds relative to the ECG waveform. Scaling transformations affect the
ratio of signal variance to discretisation interval widths, potentially
leading to incorrect discretisation. Lastly, discretisation methods require
multi-channel spiking, increasing model complexity. To address this problem,
delta encoding[2, 5], also known as “temporal contrast” [14, 16], yields a
signed spike when the rate of change in the signal exceeds a predefined
threshold. However, parameterisation may be challenging, and fixed thresholds
may misrepresent data under gain changes.

The “Crossing Time-Encoding Machine” (CTEM) emits spikes when the input
waveform intersects a regular reference waveform [17]. The “Crossing Time
Encoding Machine with Feedback” (FTEM) introduces resetting of the reference
waveform by the input to CTEM [17]. The integrate-and-fire time encoding
machine (IFTEM) [18, 19] applies an integrate-and-fire neuron to the input
waveform, emitting spikes when the membrane potential reaches a certain
threshold. Due to the bipolar nature of ECG signal, for our baseline
comparison studies, we extend the IFTEM to have dual threshold such that both
positive and negative spikes can be emitted. Some spiking neural network
classifiers have employed recurrent architectures receiving “direct input” of
the unencoded analogue current to the first layer of integrate-and-fire
neurons [20, 21, 22]. This approach is functionally equivalent to real-time
IFTEM, with empirical fitting of integrate and fire parameters.

##  II Derived Peak Encoding

To address the vulnerabilities of threshold and discretisation-based encoding
schemes, we propose “Derived Peak” (DP) encoding for event-driven signalling.
DP is physiologically inspired by the peaks and nadirs which the P-wave, QRS
complex and T-wave represent. Accordingly, DP generates signed spikes at zero-
crossings of the signal’s time derivatives such that:

$S_{DP}(t)_{k}=\begin{cases}1&\frac{d^{k}x}{dt}(t)>0,\frac{d^{k}x}{dt}(t-1)\leq
0\\\ -1&\frac{d^{k}x}{dt}(t)<0,\frac{d^{k}x}{dt}(t-1)\geq 0\\\
0&\text{otherwise}\end{cases}$ (1)

Consequently, DP may be considered a variant of the CTEM class of encoding
methods, with the distinctions that 1) DP operates on the derivatives of the
input to achieve location invariance, and 2) the zero constant is used as the
crossing reference function, removing requirements for user-defined
parameterisation. The zero crossings of the first derivative correspond
temporally to peaks and nadirs of the input signal, while those of the second
derivative correspond temporally to peaks and nadirs in the rate of change of
the input signal (Fig. 1). Thus, it is observed that the information encoded
by DP is directly related to the inflections after which the P-wave, QRS
complex and T-wave are named[23]. Furthermore, as
$\operatorname{sgn}(\frac{d^{k}x}{dt})=\operatorname{sgn}(\frac{d^{k}(ax+b)}{dt})$
for all $k\geq 1$ and $a>0$, DP encoding is invariant to gain changes and
shifts in the ECG waveform.

We propose two DP encoding schemes:

  * $\bullet$

_DP(1)_ : DP encoding utilising the first order derivative of the signal
($k=1$ in equation (1)). The signal is presented to the network with 1 encoded
channel per input channel.

  * $\bullet$

_$DP(1,2)$_ : DP encoding utilising both the first ($k=1)$ and the second
order derivatives ($k=2$) of the signals. The signal is presented to the
network with 2 encoded channels per input channel.

Figure 1: DP encoding scheme, for first and second-order derivatives, for a
sample ECG.

##  III Experiments and Data

Electrocardiography data was extracted from the PTB-XL dataset via the
PhysioNet repository [24, 25]. 21,837 12-lead ECG records sampled at 100 Hz
were used in the study. 11,354 male and 10445 female patients were included,
with median age 62 years (IQR: 50-72 years), median height 166 cm (IQR:
160-174 cm) and median weight 70 kg (IQR: 60-80 kg). 285 participants had
pacemakers. Three (non-mutually exclusive) classes of cardiovascular
conditions relating to ECG waveform morphology were annotated: ST-
segment/T-wave changes (STTC), myocardial infarction (MI) and conduction
deficits (CD). 5,235 samples (24%) demonstrated STTC, 5,469 (25%) demonstrated
MI and 4,898 (23%) demonstrated CD. Following the dataset authors’
instructions, predefined folds 1-9 were used for model development, and fold
10 was reserved for testing.

Figure 2: Experimental flow to test ResNet classification performance of
spike-encoded ECG signals. To model the continuous low-latency ECG monitoring
setting, subseries of 2,560ms duration were extracted from random timepoints.
Signal pre-processing was limited to 3-point moving average (MA) for reduction
of high-frequency noise.

The objective was to predict STTC, MI and CD from a 256-timestep, 12-lead ECG
segment. A 256-timestep window was employed to assess the classification
performance in the low latency setting, including approximately three beats in
the temporal field. To recreate the most challenging and clinically useful
scenario, models were provided with a random ECG segment, without waveform
annotation or centering (Fig. 2). Preprocessing was limited to a 3-point
moving average, which served to reduce high-frequency noise. Event-driven
encoding methods were applied to the ECG signal before it was passed to a 1D
ResNet18 classifier for training/inference.

Delta threshold parameterisation followed Corradi et al. [2]. Discretisation-
based methods (Count, LC, FR, BitString) employed 31 bins of width $0.05\,$mV
over the range $[-0.55,1]$. CTEM was implemented with sine reference wave
centered at $0\,$mV with $2\,$mV amplitude and 10-timestep wavelength
($10\,$Hz), for a target spike probability of $0.2$. FTEM was parameterised
for a target spike probability of $0.2$ and IFTEM was implemented with
threshold, leakage and null bias targeting spiking probability of $0.05-0.1$.

Figure 3: ECG sample corrupted with artefacts used in the study

The training was performed with Tensorflow version 2.10 [26]. Encoded inputs
were passed to a 1D ResNet-18 architecture with a base feature depth of 64
convolutional filters. Global max-pooling was applied in the penultimate layer
and passed to a 3-channel dense layer with sigmoid activation [27]. 3-channel
binary cross-entropy was minimised over batches of size 32, using the Adam
optimizer with an initial learning rate of $10^{-3}$ and a decay factor of
$10^{-7}$ per batch. Discrimination and calibration were evaluated on the test
set using Area under Curve (AUC) and accuracy metrics. To assess robustness,
further model testing was performed with four common artefacts (see Fig. 3):

  1. 1.

Shift: $x(t)+a,\hskip 5.69054pta=\pm 1\,$mV

  2. 2.

Rescale: $ax(t),\hskip 5.69054pta\in\\{\frac{1}{3},3\\}$

  3. 3.

Drift: $x(t)+\sin\left(\frac{(256+t+a)\pi}{512}\right),\hskip
5.69054pta\in[0,256]$

  4. 4.

Noise: $x(t)+\epsilon,\hskip 5.69054pt\epsilon\sim\mathcal{N}(0,0.02)$

Shift represents a change in the location of the ECG waveform baseline
corresponding to an approximate QRS complex amplitude. Rescale represents a
three-fold gain change. Drift is low frequency ($0.512\,$Hz) baseline drift
oscillation with $1\,$mV amplitude. Noise represents spontaneous Gaussian
interference with $10\%$ of signal standard deviation.

##  IV Results

AUC and accuracy for each encoding method and artefact condition (averaged
over diagnostic class) are provided in Table I. We highlight the following
salient observations from our experiments:

  * $\bullet$

In the absence of any artefact, the $DP(1,2)$ encoding schema yielded
performance similar to all other methods, including unencoded inputs.

  * $\bullet$

In the presence of shift artefacts, models with DP(1) and $DP(1,2)$ encoded
inputs maintain stable calibrative (accuracy) and discriminative (AUC)
performance. Encoding methods relying on predefined discretisation thresholds
– Rate, BitString, FR, LC – as well as TEM-based methods failed under this
perturbation. Delta encoding discriminative performance was robust to shift
artefact.

  * $\bullet$

DP-encoded inputs maintained predictive performance in the presence of rescale
artefact. All other methods, except for the one with unencoded inputs,
suffered significant AUC and accuracy loss upon input rescaling.

  * $\bullet$

DP-encoded inputs maintained predictive performance in the presence of drift
artefact. The model receiving unencoded inputs was moderately robust to drift,
as was delta encoding, but model performance with discretisation-based
encoding and TEM-based methods reduced.

  * $\bullet$

Noise artefact induced mild deterioration in performance of models with DP(1)
encoded inputs. However in the $DP(1,2)$-encoded case, model calibration
decreased significantly.

_Metric_ |  _Encoding_ |  _Original_ |  _Shift_ |  _Rescale_ |  _Drift_ |  _Noise_
---|---|---|---|---|---|---
AUC |  None |  0.92 |  0.62 |  0.79 |  0.89 |  0.92
|  Rate |  0.91 |  0.50 |  0.77 |  0.79 |  0.91
|  BitString |  0.91 |  0.50 |  0.77 |  0.79 |  0.91
|  Delta |  0.89 |  0.89 |  0.75 |  0.89 |  0.89
|  FR |  0.90 |  0.53 |  0.74 |  0.69 |  0.90
|  LC |  0.91 |  0.52 |  0.81 |  0.81 |  0.91
|  IFTEM |  0.87 |  0.62 |  0.77 |  0.83 |  0.86
|  CTEM |  0.87 |  0.50 |  0.8 |  0.71 |  0.87
|  FTEM |  0.87 |  0.49 |  0.81 |  0.80 |  0.87
|  DP(1) |  0.90 |  0.90 |  0.90 |  0.90 |  0.89
|  DP(1,2) |  0.91 |  0.91 |  0.91 |  0.90 |  0.88
Acc. |  None |  0.87 |  0.76 |  0.80 |  0.85 |  0.87
|  Rate |  0.87 |  0.76 |  0.78 |  0.79 |  0.87
|  BitString |  0.87 |  0.55 |  0.78 |  0.75 |  0.87
|  Delta |  0.86 |  0.86 |  0.76 |  0.85 |  0.86
|  FR |  0.86 |  0.68 |  0.73 |  0.76 |  0.86
|  LC |  0.87 |  0.52 |  0.78 |  0.75 |  0.87
|  IFTEM |  0.84 |  0.60 |  0.77 |  0.81 |  0.84
|  CTEM |  0.84 |  0.67 |  0.80 |  0.77 |  0.84
|  FTEM |  0.85 |  0.75 |  0.80 |  0.80 |  0.85
|  DP(1) |  0.86 |  0.86 |  0.86 |  0.85 |  0.84
|  DP(1,2) |  0.86 |  0.86 |  0.86 |  0.86 |  0.81
TABLE I: Test AUC and accuracy of encoding schemes under each artefact
condition

Spike output rates produced by each method in the test dataset are provided in
Table II. Rate and BitString encoding required multiple spikes per timestep.
DP, TEM methods and LC generated relatively sparse outputs. DP based methods
achieved low spike counts throughout. Of the methods that outperform DP in
efficiency, FR and LC have high output channel ratio, thereby increasing the
model complexity.

_Encoding_ |  _Output Channel Ratio_ |  _Total Channel Output_ |  _Mean Channel Output_
---|---|---|---
Delta |  1 |  0.11 |  0.11
Rate |  31 |  11.56 |  0.37
FR |  31 |  1 |  0.03
LC |  31 |  0.51 |  0.02
BitString |  7 |  3.88 |  0.55
CTEM |  1 |  0.21 |  0.21
FTEM |  1 |  0.19 |  0.19
IFTEM |  1 |  0.06 |  0.06
DP(1) |  1 |  0.2 |  0.2
DP(1,2) |  2 |  0.57 |  0.29
TABLE II: Spike output per 256-timestep ECG segment for various encoding
methods, averaged over all test segments. The output channel ratio represents
the total number of output channels per input channel. Total channel output is
the total number of spikes across all channels, and mean channel output is the
mean spikes per channel.

##  V Discussion

Without artefacts, the classification accuracy obtained with a wide variety of
spike-based encoding schemes is similar to that with analogue inputs. Test
discrimination performance was comparable to that observed in a benchmark
survey of deep learning model performance on the dataset $(AUC:0.87-0.93)$,
where the entire 1000-timestep series was provided to the model [28]. This
suggests that information loss that results due to event-driven encoding had a
negligible effect on classification performance.

However, previously reported spike-encoding methods performed poorly in the
presence of corrupting artefacts. Furthermore, models operating on the
unencoded signal were also vulnerable to such corruption. DP-encoded methods
proved resilient to shift, rescale and drift artefacts, maintaining accuracy
and AUC within a narrow 1-2% range. An added advantage of DP encoding is its
non-parametric nature, requiring minimal domain knowledge or extensive data
distribution analysis to anticipate encoding behaviour and tuning.

However, high-frequency noisy interference presented a challenge for DP(1 ,2)
based encoding, with DP(1) showing slightly better performance. High-frequency
noise is known to complicate the computation of derivatives, and
regularisation techniques have been proposed to address this problem [29].
Oscillations due to such noise in low amplitude regions in the signal may have
caused artefactual zero crossings in the second-order time derivative.
Furthermore, diagnostic tasks which involve amplitude-based diagnostic
criteria, such as hypertrophy [30], may be better served by amplitude-based
encoding schemes. Notably, DP methods achieved high performance with
relatively low spike rates, making them suitable candidates for deployment on
energy-constrained diagnostic systems at the edge.

##  VI Conclusion

DP is a nonparametric, event-driven, and resource-friendly encoding technique
which guarantees invariance to shifts and rescaling of the input signal with
minimal pre-processing. DP encoded inputs yield similar classification
performance to the unencoded input with improved robustness to three common
artefacts. DP encoding has strengths and vulnerabilities complementary to
those of currently available techniques, presenting an opportunity to improve
cross-artefact robustness through hybrid implementation.

## Acknowledgment

This research was supported in part by the EPSRC Open Fellowship EP/X011356/1.

## References

  * [1] Z. Ebrahimi, M. Loni, M. Daneshtalab, and A. Gharehbaghi, “A review on deep learning methods for ECG arrhythmia classification,” _Expert Systems with Applications: X_ , vol. 7, p. 100033, Sep. 2020.
  * [2] F. Corradi, S. Pande, J. Stuijt, N. Qiao, S. Schaafsma, G. Indiveri, and F. Catthoor, “ECG-based Heartbeat Classification in Neuromorphic Hardware,” in _2019 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2019, pp. 14–19.
  * [3] Z. Yan, J. Zhou, and W.-F. Wong, “Energy efficient ECG classification with spiking neural network,” _Biomed. Signal Process. Control_ , vol. 63, p. 102170, Jan. 2021.
  * [4] D. Banerjee, S. Dey, and A. Pal, “An snn based ecg classifier for wearable edge devices,” in _NeurIPS 2022 Workshop on Learning from Time Series for Health_ , 2022.
  * [5] F. C. Bauer, D. R. Muir, and G. Indiveri, “Real-Time Ultra-Low Power ECG Anomaly Detection Using an Event-Driven Neuromorphic Processor,” _IEEE Trans. Biomed. Circuits Syst._ , vol. 13, no. 6, pp. 1575–1582, Nov. 2019.
  * [6] H. Chu, Y. Yan, L. Gan, H. Jia, L. Qian, Y. Huan, L. Zheng, and Z. Zou, “A Neuromorphic Processing System With Spike-Driven SNN Processor for Wearable ECG Classification,” _IEEE Trans. Biomed. Circuits Syst._ , vol. 16, no. 4, pp. 511–523, Jul. 2022.
  * [7] G. Lenis, N. Pilia, A. Loewe, W. H. W. Schulze, and O. Dössel, “Comparison of Baseline Wander Removal Techniques considering the Preservation of ST Changes in the Ischemic ECG: A Simulation Study,” _Comput. Math. Methods Med._ , vol. 2017, Mar. 2017.
  * [8] V. S. Chouhan and S. S. Mehta, “Total Removal of Baseline Drift from ECG Signal,” in _2007 International Conference on Computing: Theory and Applications (ICCTA’07)_. IEEE, 2007, pp. 05–07.
  * [9] A. R. Pérez-Riera, R. Barbosa-Barros, R. Daminello-Raimundo, and L. C. de Abreu, “Main artifacts in electrocardiography,” _Annals of Noninvasive Electrocardiology : The Official Journal of the International Society for Holter and Noninvasive Electrocardiology, Inc_ , vol. 23, no. 2, Mar. 2018.
  * [10] B. Ibanez, S. James, S. Agewall, M. J. Antunes, C. Bucciarelli-Ducci, H. Bueno, A. L. P. Caforio, F. Crea, J. A. Goudevenos, S. Halvorsen, G. Hindricks, A. Kastrati, M. J. Lenzen, E. Prescott, M. Roffi, M. Valgimigli, C. Varenhorst, P. Vranckx, P. Widimský, and ESC Scientific Document Group, “2017 ESC Guidelines for the management of acute myocardial infarction in patients presenting with ST-segment elevation: The Task Force for the management of acute myocardial infarction in patients presenting with ST-segment elevation of the European Society of Cardiology (ESC),” _Eur. Heart J._ , vol. 39, no. 2, pp. 119–177, Jan. 2018.
  * [11] A. Das, P. Pradhapan, W. Groenendaal, P. Adiraju, R. T. Rajan, F. Catthoor, S. Schaafsma, J. L. Krichmar, N. Dutt, and C. Van Hoof, “Unsupervised heart-rate estimation in wearables with Liquid states and a probabilistic readout,” _Neural Networks_ , vol. 99, pp. 134–147, Mar. 2018.
  * [12] M. Saeed, Q. Wang, O. Märtens, B. Larras, A. Frappé, B. Cardiff, and D. John, “Evaluation of Level-Crossing ADCs for Event-Driven ECG Classification,” _IEEE Trans. Biomed. Circuits Syst._ , vol. 15, no. 6, pp. 1129–1139, Dec. 2021.
  * [13] S. Somani, A. J. Russak, F. Richter, S. Zhao, A. Vaid, F. Chaudhry, J. K. De Freitas, N. Naik, R. Miotto, G. N. Nadkarni, J. Narula, E. Argulian, and B. S. Glicksberg, “Deep learning and the electrocardiogram: review of the current state-of-the-art,” _Europace_ , vol. 23, no. 8, p. 1179, Aug. 2021.
  * [14] D. Auge, J. Hille, E. Mueller, and A. Knoll, “A survey of encoding techniques for signal processing in spiking neural networks,” _Neural Processing Letters_ , vol. 53, no. 6, p. 4693–4710, Jul 2021. [Online]. Available: https://link.springer.com/article/10.1007/s11063-021-10562-2
  * [15] P. Date, S. Kulkarni, A. Young, C. Schuman, T. Potok, and J. Vetter, “Encoding integers and rationals on neuromorphic computers using virtual neuron,” _Sci. Rep._ , vol. 13, no. 10975, pp. 1–18, Jul. 2023.
  * [16] B. Petro, N. Kasabov, and R. M. Kiss, “Selection and Optimization of Temporal Spike Encoding Methods for Spiking Neural Networks,” _IEEE Trans. Neural Networks Learn. Syst._ , vol. 31, no. 2, pp. 358–370, Apr. 2019.
  * [17] D. Gontier and M. Vetterli, “Sampling based on timing: Time encoding machines on shift-invariant subspaces,” _Appl. Comput. Harmon. Anal._ , vol. 36, no. 1, pp. 63–78, Jan. 2014.
  * [18] A. A. Lazar, “Time encoding with an integrate-and-fire neuron with a refractory period,” _Neurocomputing_ , vol. 58-60, pp. 53–58, Jun. 2004.
  * [19] H. Naaman, S. Mulleti, and Y. C. Eldar, “FRI-TEM: Time Encoding Sampling of Finite-Rate-of-Innovation Signals,” _IEEE Trans. Signal Process._ , vol. 70, pp. 2267–2279, Apr. 2022.
  * [20] P. Kovács and K. Samiee, “Arrhythmia Detection Using Spiking Variable Projection Neural Networks,” in _2022 Computing in Cardiology (CinC)_. IEEE, 2022, pp. 04–07.
  * [21] H. De Melo Ribeiro, A. Arnold, J. P. Howard, M. J. Shun-Shin, Y. Zhang, D. P. Francis, P. B. Lim, Z. Whinnett, and M. Zolgharni, “ECG-based real-time arrhythmia monitoring using quantized deep neural networks: A feasibility study,” _Comput. Biol. Med._ , vol. 143:105249., Apr. 2022.
  * [22] X. Liang, H. Li, A. Vuckovic, J. Mercer, and H. Heidari, “A Neuromorphic Model With Delay-Based Reservoir for Continuous Ventricular Heartbeat Detection,” _IEEE Trans. Biomed. Eng._ , vol. 69, no. 6, pp. 1837–1849, Jun. 2022.
  * [23] Y. Sattar and L. Chhabra, “Electrocardiogram,” in _StatPearls [Internet]_. StatPearls Publishing, Jun. 2023. [Online]. Available: https://www.ncbi.nlm.nih.gov/books/NBK549803
  * [24] P. Wagner, N. Strodthoff, R.-D. Bousseljot, D. Kreiseler, F. I. Lunze, W. Samek, and T. Schaeffter, “PTB-XL, a large publicly available electrocardiography dataset,” _Sci. Data_ , vol. 7, no. 154, pp. 1–15, May 2020.
  * [25] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C. K. Peng, and H. E. Stanley, “PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals,” _Circulation_ , vol. 101, no. 23, pp. 215–220, Jun. 2000.
  * [26] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/
  * [27] K. He, X. Zhang, S. Ren, and J. Sun, “Deep Residual Learning for Image Recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2016, pp. 27–30.
  * [28] N. Strodthoff, P. Wagner, T. Schaeffter, and W. Samek, “Deep learning for ecg analysis: Benchmarks and insights from ptb-xl,” _IEEE Journal of Biomedical and Health Informatics_ , vol. 25, no. 5, pp. 1519–1528, 2020.
  * [29] F. van Breugel, J. N. Kutz, and B. W. Brunton, “Numerical differentiation of noisy data: A unifying multi-objective optimization framework,” _IEEE access : practical innovations, open solutions_ , vol. 8, p. 196865, 2020.
  * [30] J. G. Peguero, S. Lo Presti, J. Perez, O. Issa, J. C. Brenes, and A. Tolentino, “Electrocardiographic Criteria for the Diagnosis of Left Ventricular Hypertrophy,” _J. Am. Coll. Cardiol._ , vol. 69, no. 13, pp. 1694–1703, Apr. 2017.


# Convex ordering of Pólya random variables and monotonicity of the error
estimate of Bernstein-Stancu operators

Diana Meleşteu  Faculty of Mathematics and Computer Science, Transilvania
University of Braşov, Str. Iuliu Maniu 50, Braşov – 500091, Romania.
<EMAIL_ADDRESS>,  Mihai N. Pascu  Faculty of Mathematics and
Computer Science, Transilvania University of Braşov, Str. Iuliu Maniu 50,
Braşov – 500091, Romania<EMAIL_ADDRESS>and  Nicolae R. Pascu
Department of Mathematics, Southern Polytechnic State University, 1100 South
Marietta Parkway, Marietta, GA 30060-2896, U.S.A<EMAIL_ADDRESS>
###### Abstract.

In the present paper we show that in Pólya’s urn model, for an arbitrarily
fixed initial distribution of the urn, the corresponding random variables
satisfy a convex ordering with respect to the replacement parameter. As an
application, we show that in the class of convex functions, the absolute value
of the error of Bernstein-Stancu operators is a non-decreasing (strictly
increasing under an additional hypothesis) function of the corresponding
parameter.

The proof relies on two results of independent interest: an interlacing lemma
of three sets and the monotonicity of the (partial) first moment of Pólya
random variables with respect to the replacement parameter.

###### Key words and phrases:

Convex ordering, Pólya urn model, error estimate of Bernstein-Stancu
operators, Bernstein operator, approximation theory

###### 2020 Mathematics Subject Classification:

Primary 60E15, 41A36, 41A10, 47A63, 47A58

##  1\. Introduction

More than $100$ years ago, in a beautiful and short paper Serge Bernstein
([1]) gave a simple, constructive proof of Weierstrass’s theorem on uniform
approximation of continuous functions by polynomials, known nowadays as
Bernstein polynomials.

About $50$ years later, D. D. Stancu noticed that the binomial distribution
used by Bernstein is a particular case of the Pólya’s urn distribution (the
case of the replacement parameter being equal to zero), and he introduced
([6], [7]) a more general class of polynomials/operators, known in the
literature as the Pólya-Stancu operators (the operator $P_{n}^{c}$ defined by
(2.7)).

Aside from a passing remark that for a particular choice (negative real
number) of the replacement parameter $c$ one obtains the Lagrange
interpolation polynomial (which cannot be used for uniform approximation), in
his work Stancu considered only non-negative values $c\geq 0$ of the
replacement parameter in Pólya’s urn model, this choice being adopted by
subsequent researchers in the field.

Recently (in [3], see also [4] and [9]), the last two authors introduced the
operator $R_{n}$ (given by (2.8)), corresponding to a negative choice
(pointwise minimal) of the replacement parameter of Pólya’s urn distribution,
and showed that this leads to better approximation results. To be precise, we
showed that the upper bounds of the error estimates (in terms of the
first/second order modulus of continuity of the function, etc) are smaller
than the corresponding estimates for the Bernstein operator, and we also
provided numerical evidence (for various choices of the function, smooth,
continuous, even discontinuous - see [3]) which indicated that among all
Bernstein-Stancu type operators, the newly introduced operator $R_{n}$ gives
the best approximation.

A criticism received while publishing these results was that even though the
upper bounds for the error of the new operator are smaller than the
corresponding ones for the Bernstein operator, this is not a proof that the
new operator is a better approximation operator.

In the present paper, we fill this gap (at least partially), by showing that
for a sufficiently large class of functions (convex functions) the absolute
value of the error of approximation of Bernstein-Stancu type operators
$P_{n}^{c}$ is (pointwise) a monotonically increasing function of the
replacement parameter $c$. This shows that in the case of operators considered
by Stancu (non-negative replacement parameter), the Bernstein operator $B_{n}$
gives the best approximation of convex functions, and that the choice of the
operator $R_{n}$ considered by the last two authors (minimal admissible choice
of the replacement parameter) further improves this approximation, giving the
best approximation of convex functions.

The proof of the main result (Theorem 4.1) relies on a result of independent
interest (Theorem 3.3), which shows that the Pólya random variables satisfy a
convex ordering with respect to the replacement parameter. In turn, the proof
of this result relies on two other results with interest on their own: a
result concerning the ordering (interlacing) of three sets (Lemma 3.1) and the
monotonicity of the partial centered moment of the Pólya distribution (Lemma
3.2).

The structure of the paper is the following. In Section 2 we introduce the
notation and the main results needed in the sequel.

In Section 3 we prove the convex ordering of Pólya random variables, and, as
an application, in Section 4 we prove that in the case of convex functions,
the absolute value of the error of approximation of the Bernstein-Stancu
operators is a non-decreasing function of the replacement parameter (strictly
increasing under an additional hypothesis).

##  2\. Preliminary results

Recall that for given integer parameters $a,b\geq 0$, $c$, and $n\geq 1$
satisfying the compatibility condition

(2.1) $a+\left(n-1\right)c\geq 0\text{\qquad and\qquad}b+\left(n-1\right)c\geq
0\text{,}$

Pólya’s urn model gives the number of white balls (”successes”) extracted in
$n$ trials from an urn containing initially $a$ white balls and $b$ black
balls, where after each extraction, the extracted ball is replaced in the urn
together with $c$ balls of the same color. Denoting by $X_{n}^{a,b,c}$ the
random variable representing the number of successes in this experiment (which
will be referred to as a Pólya random variable with parameters $a,b,c$ and
$n$), we have

(2.2)
$p_{n,k}^{a,b,c}=P\left(X_{n}^{a,b,c}=k\right)=C_{n}^{k}\frac{a^{\left(k,c\right)}b^{\left(n-k,c\right)}}{(a+b)^{\left(n,c\right)}},\qquad
k\in\\{0,1,\ldots,n\\},$

where for $x,h\in\mathbb{R}$ and $n\in\mathbb{N}$ we denoted by

(2.3)
$x^{\left(n,h\right)}=x\left(x+h\right)\left(x+2h\right)\cdot\ldots\cdot\left(x+\left(n-1\right)h\right)$

the generalized (rising) factorial with increment $h$. We are using the
convention that an empty product is equal to $1$, that is
$x^{\left(0,h\right)}=1$ for any $x,h\in\mathbb{R}$. The binomial theorem for
the rising factorial shows that (2.2) still defines a distribution
$X_{n}^{a,b,c}$ for real values of the parameters $a,b\geq 0$ and $c$,
satisfying the compatibility condition (2.1).

It is known (see [2]) that the partial first absolute centered moment of the
Pólya distribution $X_{n}^{x,1-x,c}$ is given by

(2.4) $E\left(\left(nx-X_{n}^{x,1-x,c}\right)1_{X_{n}^{x,1-x,c}\leq
s-1}\right)=\sum_{k=0}^{s-1}\left(nx-k\right)p_{n,k}^{x,1-x,c}=sp_{n,s}^{x,1-x,c}\left(1-x+\left(n-s\right)c\right)$

for all $s\in\left\\{1,\ldots,n\right\\}$.

Denoting by $\mathcal{F}\left(\left[0,1\right]\right)$ the set of real-valued
functions defined on $\left[0,1\right]$, for an integer $n\geq 1$ denote by
$P_{n}^{a,b,c}:\mathcal{F}\left(\left[0,1\right]\right)\rightarrow\mathcal{F}\left(\left[0,1\right]\right)$
the operator defined by

(2.5)
$P_{n}^{a,b,c}\left(f;x\right)=E\left(\frac{1}{n}X_{n}^{a,b,c}\right)=\sum_{k=0}^{n}f\left(\frac{k}{n}\right)p_{n}^{a,b,c}$

where the parameters $a,b\geq 0$ and
$c\geq-\frac{\min\left\\{x,1-x\right\\}}{n-1}$ may depend on $n$ and
$x\in\left[0,1\right]$ and satisfy the compatibility condition (2.1) (see
[3]). In particular, consider the Bernstein operator

(2.6) $B_{n}\left(f;x\right)=P_{n}^{x,1-x,0}\left(f;x\right),$

the Bernstein-Stancu operator

(2.7) $P_{n}^{c}\left(f;x\right)=P_{n}^{x,1-x,c}\left(f;x\right),$

and the operator $R_{n}$ introduced by the last two authors in [3]

(2.8)
$R_{n}\left(f;x\right)=P_{n}^{x,1-x,-\min\left\\{x,1-x\right\\}/\left(n-1\right)},$

which corresponds to the minimal value of the replacement parameter $c$
(satisfying the compatibility condition (2.1)) in the case $a=x$ and $b=1-x$.

Finally, recall that a random variable $X$ is said to be smaller in the convex
order than a random variable $Y$ (denoted by $X\leq_{\text{cx}}Y$) iff

(2.9) $E\left(\phi(X)\right)\leq E\left(\phi(Y)\right)$

for any convex function $\phi:\mathbb{R}\rightarrow\mathbb{R}$ for which the
above expectations exist.

If $X$ and $Y$ are random variables for which the means $EX$, $EY$ exist and
are equal, it is known (see e.g. [5], Theorem 3.A.1) that $X\leq_{\text{cx}}Y$
iff

(2.10) $E\left((t-X)_{+}\right)\leq E\left((t-Y)_{+}\right),\qquad\text{for
all }t\in\mathbb{R},$

where $x_{+}=\max\\{x,0\\}$ denotes the positive part of $x\in\mathbb{R}$.

##  3\. Main results

In order to prove the main result of this section, we begin with the following
auxiliary result, of independent interest.

######  Lemma 3.1.

For any real number $x\in\left(0,1\right)$, integers $n\geq 3$ and
$k\in\left\\{1,\ldots,n-2\right\\}$, there exists a disjoint partition
$\left\\{n_{1},\ldots,n_{k}\right\\}\bigsqcup\left\\{m_{1},\ldots,m_{n-k-1}\right\\}$
of the set $\left\\{1,2,\ldots,n-1\right\\}$ such that

(3.1) $n_{i}\leq\frac{i}{x},\quad
i\in\left\\{1,\ldots,k\right\\}\qquad\text{and}\qquad
m_{i}\leq\frac{i}{1-x},\quad i\in\left\\{1,\ldots,n-k-1\right\\}.$

###### Proof.

First note that for $i\geq 1$ the interval $\left(i,i+1\right)$ cannot contain
two distinct elements of the set
$\left\\{\frac{1}{x},\ldots,\frac{k}{x}\right\\}$. This is so for otherwise,
there would exist indices $1\leq j_{1}<j_{2}\leq k$ such that
$i<\frac{j_{1}}{x}<\frac{j_{2}}{x}<i+1$, and therefore

$\frac{1}{x}\leq\frac{j_{2}-j_{1}}{x}<i+1-i=1,$

which implies $x>1$, a contradiction. A similar proof shows that the interval
$\left(i,i+1\right)$ cannot contain two distinct elements of the set
$\left\\{\frac{1}{1-x},\ldots,\frac{n-k-1}{1-x}\right\\}$.

Secondly, note that the interval $\left(i,i+1\right)$ cannot contain an
element of the set $\left\\{\frac{1}{x},\ldots,\frac{k}{x}\right\\}$ and an
element of the set $\left\\{\frac{1}{1-x},\ldots,\frac{n-k-1}{1-x}\right\\}$.
This is so for otherwise there would exist indices $1\leq j_{1}\leq k$ and
$1\leq j_{2}\leq n-k-1$ such that $i<\frac{j_{1}}{x}<i+1$ and
$i<\frac{j_{2}}{1-x}<i+1$, and therefore

$i=ix+i\left(1-x\right)<j_{1}+j_{2}<\left(i+1\right)x+\left(i+1\right)\left(1-x\right)=i+1,$

which is impossible since $j_{1}+j_{2}$ is an integer.

A similar proof shows that if the interval $\left[i,i+1\right]$ contains an
element of the set $\left\\{\frac{1}{x},\ldots,\frac{k}{x}\right\\}$ and an
element of the set $\left\\{\frac{1}{1-x},\ldots,\frac{n-k-1}{1-x}\right\\}$,
then they both must be equal either to $i$ or to $i+1$. Note that in the
latter case the interval $(i+1,i+2]$ cannot contain any element of the set
$\left\\{\frac{1}{x},\ldots,\frac{k}{x}\right\\}\cup\left\\{\frac{1}{1-x},\ldots,\frac{n-k-1}{1-x}\right\\}$
(since $x<1$ and $1-x<1$).

Define the sequences $\left(n_{i}^{\prime}\right)_{i=\overline{1,k}}$ and
$\left(m_{i}^{\prime}\right)_{i=\overline{1,n-k-1}}$ by $n_{i}^{\prime}=j$ if
$\frac{i}{x}\in(j,j+1]$ for some $j\in\mathbb{N}$,
$i\in\left\\{1,\ldots,k\right\\}$, and

$m_{i}^{\prime}=\left\\{\begin{array}[]{ll}j,&\text{if
}\frac{i}{1-x}\in(j,j+1)\text{ for some }j\in\mathbb{N}\\\ j+1,&\text{if
}\frac{i}{1-x}=j+1\text{ for some }j\in\mathbb{N}\end{array}\right.,\qquad
i\in\left\\{1,\ldots,n-k-1\right\\},$

thus $n_{i}^{\prime}\leq\frac{i}{x}$ for $i\in\left\\{1,\ldots,k\right\\}$ and
$m_{i}^{\prime}\leq\frac{i}{1-x}$ for $i\in\left\\{1,\ldots,n-k-1\right\\}$.

The first part of the proof shows that $n_{1}^{\prime},\ldots
n_{k}^{\prime},m_{1}^{\prime},\ldots,m_{n-k-1}^{\prime}$ are all distinct, and
in particular this shows that the set $N=\left\\{n_{1}^{\prime},\ldots
n_{k}^{\prime},m_{1}^{\prime},\ldots,m_{n-k-1}^{\prime}\right\\}$ has $n-1$
elements.

Denoting by $r$ the number of distinct elements of the set
$N\cap\left(n-1,\infty\right)$, it follows that the set
$\left\\{1,\ldots,n-1\right\\}\backslash\left(N\cap[0,n-1]\right)$ also has
$r$ elements, and therefore there exists a bijection
$f:N\cap\left(n-1,\infty\right)\rightarrow\left\\{1,\ldots,n-1\right\\}\backslash\left(N\cap[0,n-1]\right)$.

Define

$n_{i}=\left\\{\begin{array}[]{ll}n_{i},&\text{if }n_{i}^{\prime}\leq n-1\\\
f\left(n_{i}^{\prime}\right),&\text{if
}n_{i}^{{}^{\prime}}>n-1\end{array}\right.,\qquad i\in\\{1,\ldots,k\\},$

and

$m_{i}=\left\\{\begin{array}[]{ll}m_{i}^{\prime},&\text{if }m_{i}^{\prime}\leq
n-1\\\ f\left(m_{i}^{\prime}\right),&\text{if
}m_{i}^{\prime}>n-1\end{array}\right.,\qquad
i\in\left\\{1,\ldots,n-k-1\right\\}.$

It is not difficult to see that $n_{1},\ldots,n_{k},m_{1},\ldots,m_{n-k-1}$
are distinct,
$\left\\{1,\ldots,n-1\right\\}=\left\\{n_{1},\ldots,n_{k},m_{1},\ldots,m_{n-k-1}\right\\}$
(recall the definition of the function $f,$ in particular $f\leq n-1$), and
they satisfy

$n_{i}\leq\frac{i}{x},\text{\quad}i\in\left\\{1,\ldots,k\right\\}\qquad\text{and}\qquad
m_{i}\leq\frac{i}{1-x},\text{\quad}i\in\left\\{1,\ldots,n-k-1\right\\},$

concluding the proof. ∎

A second auxiliary result of independent interest is the following
monotonicity of the (partial) first centered moment of Pólya random variables.

######  Lemma 3.2.

For any $x\in\left[0,1\right]$ and any integers $n\geq 1$ and
$k\in\left\\{0,1,\ldots,n\right\\}$, the sum

(3.2) $\sum_{i=0}^{k}\left(x-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}$

is a non-decreasing function of $c\geq-\frac{1}{n-1}\min\\{x,1-x\\}$.

Moreover, for $x\in\left(0,1\right)$, $n\geq 2,$ and
$k\in\left\\{0,\ldots,n-1\right\\}$, the above sum is increasing with respect
to $c\geq-\frac{1}{n-1}\min\\{x,1-x\\}$.

###### Proof.

For $k=n$ we have

(3.3)
$\sum_{i=0}^{n}\left(x-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}=x-\frac{1}{n}EX_{n}^{x,1-x,c}=x-x=0,$

and for $k=0$ we have

(3.4)
$\sum_{i=0}^{0}\left(x-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}=xp_{n,0}^{x,1-x,c}=x\left(1-x\right)\left(1-x+c\right)\cdot\ldots\cdot\left(1-x+(n-1)c\right),$

thus the claim of the lemma is true in these cases.

Also, since for $x=0$ and $x=1$ the probabilities $p_{n,i}^{x,1-x,c}$ are
independent on the value of $c$ ($p_{n,0}^{0,1,c}=p_{n,n}^{1,0,c}=1$, and
$p_{n,i}^{x,1-x,c}=0$ for $i\in\left\\{1,\ldots,n-1\right\\}$), the claim of
the lemma is also true in these cases.

Without loss of generality we may therefore assume that $x\in(0,1)$,
$k\in\\{1,\ldots,n-1\\}$, and $n\geq 2$. Since the expression in (3.2) is a
continuous, differentiable function of $c\geq-\frac{1}{n-1}\min\\{x,1-x\\}$,
in order to prove the claim of the lemma, using (2.4), it suffices to show
that for any $c>-\frac{1}{n-1}\min\\{x,1-x\\}$ we have

(3.5) $\frac{\partial}{\partial
c}\sum_{i=0}^{k}\left(nx-i\right)p_{n,i}^{x,1-x,c}=\frac{\partial}{\partial
c}\left(\left(k+1\right)p_{n,k+1}^{x,1-x,c}\left(1-x+\left(n-\left(k+1\right)\right)c\right)\right)>0.$

Taking logarithms and using (2.2), we have left to prove that

$\frac{\partial}{\partial
c}\left(\ln\left(\left(k+1\right)C_{n}^{k+1}\frac{x^{\left(k+1,c\right)}\left(1-x\right)^{\left(n-k-1,c\right)}\left(1-x+\left(n-k-1\right)c\right)}{1^{\left(n,c\right)}}\right)\right)>0,$

or equivalent

(3.6)
$\sum_{i=0}^{k}\frac{i}{x+ic}+\sum_{i=0}^{n-k-1}\frac{i}{1-x+ic}>\sum_{i=0}^{n-1}\frac{i}{1+ic},$

for any $c>-\frac{1}{n-1}\min\\{x,1-x\\}$.

For $k=n-1$ the above inequality is readily satisfied (recall that
$x\in\left(0,1\right)$), so we have left to consider the case $n\geq 3$ and
$k\in\left\\{1,\ldots,n-2\right\\}$.

For arbitrarily fixed $n\geq 3$, $k\in\left\\{1,\ldots,n-2\right\\}$,
$x\in\left(0,1\right)$, and $c>-\frac{1}{n-1}\min\\{x,1-x\\}$, the function
$\varphi:\left[1,\frac{n-1}{\min\left\\{x,1-x\right\\}}\right]\rightarrow\mathbb{R}$
defined by $\varphi\left(u\right)=\frac{u}{1+uc}$ is increasing.

Lemma 3.1 shows that we can find a partition
$\left\\{n_{1},\ldots,n_{k}\right\\}\bigsqcup\left\\{m_{1},\ldots,m_{n-k-1}\right\\}$
of the set $\left\\{1,2,\ldots,n-1\right\\}$ such that $n_{i}\leq\frac{i}{x}$
for $i\in\left\\{1,\ldots,k\right\\}$ and $m_{i}\leq\frac{i}{1-x}$ for
$i\in\left\\{1,\ldots,n-k-1\right\\}$, and therefore we obtain

$\sum_{i=1}^{k}\varphi\left(\frac{i}{x}\right)+\sum_{i=1}^{n-k-1}\varphi\left(\frac{i}{1-x}\right)>\sum_{i=1}^{k}\varphi\left(n_{i}\right)+\sum_{i=1}^{n-k-1}\varphi\left(m_{i}\right)=\sum_{i=1}^{n-1}\varphi\left(i\right),$

which is equivalent to (3.6), thus concluding the proof. ∎

With the above preparation we can now proceed to prove the convex ordering of
Pólya random variables $X_{n}^{a,b,c}$ with respect to the replacement
parameter $c$. The precise statement is the following.

######  Theorem 3.3.

The Pólya random variables $X_{n}^{x,1-x,c}$ satisfy the following convex
ordering

(3.7) $X_{n}^{x,1-x,c}\leq_{\text{cx}}X_{n}^{x,1-x,c^{\prime}},$

for any integer $n\geq 1$, $x\in\left[0,1\right]$, and any $c^{\prime}\geq
c\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$.

###### Proof.

Since $EX_{n}^{x,1-x,c}=nx$ for any value of $c$ satisfying the compatibility
condition $c\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$, the claim of the
theorem is equivalent by (2.10) to

(3.8) $\frac{\partial}{\partial
c}E\left(\left(t-X_{n}^{x,1-x,c}\right)_{+}\right)\geq 0,\qquad
t\in\mathbb{R}.$

Since $X_{n}^{x,1-x,c}$ takes values in $\left\\{0,1,\ldots,n\right\\}$ and
$EX_{n}^{x,1-x,c}=nx$ is independent of $c$, it is easy to see that the above
inequality is satisfied for $t<0$ and $t>n$, thus it remains to prove it for
$t\in[0,n]$.

Replacing $t$ by $nt$ with $t\in[0,1]$, the above inequality is thus
equivalent to

(3.9) $\frac{\partial}{\partial
c}E\left(t-\frac{1}{n}X_{n}^{x,1-x,c}\right)_{+}=\sum_{i=0}^{\left[nt\right]}\left(\left(t-\frac{i}{n}\right)\frac{\partial}{\partial
c}p_{n,i}^{x,1-x,c}\right)\geq 0,\qquad t\in\left[0,1\right].$

The above inequality is further equivalent to the apparent weaker inequality

(3.10)
$\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial}{\partial
c}p_{n,i}^{x,1-x,c}\right)\geq 0,\qquad k\in\left\\{0,1,\ldots,n\right\\},$

the reason being the following.

For $t=1=\frac{n}{n}$ the inequality (3.9) follows immediately from (3.10).
For $t\in[0,1)$, we can write
$t=\alpha\frac{k}{n}+\left(1-\alpha\right)\frac{k+1}{n}$ as a convex
combination of $\frac{k}{n}$ and $\frac{k+1}{n}$, where $k=\left[nt\right]$
and $\alpha=k+1-nt\in[0,1)$. If the inequality (3.10) holds true, then

$\displaystyle\frac{\partial}{\partial
c}\sum_{i=0}^{[nt]}\left(\left(t-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
c}\left(\alpha\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
c}\left(\alpha\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k+1}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)\right)$
$\displaystyle=$
$\displaystyle\alpha\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial
c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k+1}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle\geq$ $\displaystyle 0,$

by Lemma 3.2, thus proving (3.9).

In order to prove (3.10), first note that the claim is trivial if $x=1$, for
in this case $p_{n,n}^{1,0,c}=1$ and $p_{n,i}^{1,0,c}=0$ for
$i\in\left\\{0,\ldots,n-1\right\\}$ (thus $p_{n,i}^{1,0,c}$ is independent of
the value of the replacement parameter $c$). Without loss of generality we can
therefore assume that $x\neq 1$.

Suppose that (3.10) does not hold for a certain $x\in[0,1)$ and
$k\in\\{0,1,\ldots,n\\}$, that is

(3.11)
$\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial}{\partial
c}p_{n,i}^{x,1-x,c}\right)<0.$

We distinguish the following cases.

  * i)

$nx>k$

Note that in this case we cannot have $k=0$, for in this case the sum in
(3.11) is equal to $0$. Since $k\geq 1$ and $x\in[0,1)$ we can write
$x=\alpha\frac{k-1}{n}+\left(1-\alpha\right)\frac{k}{n}$, where
$\alpha=k-nx<0$.

Using again Lemma 3.2 and the same argument as above we obtain

$\displaystyle\alpha\sum_{i=0}^{k-1}\left(\left(\frac{k-1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial
c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle=$
$\displaystyle\alpha\sum_{i=0}^{k-1}\left(\left(\frac{k-1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial
c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k-1}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
c}\sum_{i=0}^{k-1}\left(\left(\alpha\left(\frac{k-1}{n}-\frac{i}{n}\right)+\left(1-\alpha\right)\left(\frac{k}{n}-\frac{i}{n}\right)\right)p_{n,i}^{x,1-x,c}\right)$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{k-1}\left(\left(x-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle\geq$ $\displaystyle 0,$

and therefore (the second sum on the first line above being assumed to be
strictly negative, and since $\alpha<0$) we deduce that

(3.12)
$\sum_{i=0}^{k-1}\left(\left(\frac{k-1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)<0.$

We showed that if the sum in (3.11) is strictly negative, then so is the sum
in (3.12). Proceeding inductively on $k$, this implies that

$0\equiv\sum_{i=0}^{0}\left(\left(\frac{0}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)<0,$

a contradiction.

  * ii)

$nx\leq k<k+1$

Note that in this case we cannot have $k=n$, for

(3.13)
$\sum_{i=0}^{n}\left(\left(\frac{n}{n}-\frac{i}{n}\right)\frac{\partial}{\partial
c}p_{n,i}^{x,1-x,c}\right)=\frac{\partial}{\partial
c}E\left(1-\frac{1}{n}X_{n}^{x,1-x,c}\right)=\frac{\partial}{\partial
c}\left(1-x\right)\equiv 0.$

Since $k\leq n-1$ and $x\in[0,1)$ we can write
$x=\alpha\frac{k}{n}+\left(1-\alpha\right)\frac{k+1}{n}$, where $k+1\leq n$
and $\alpha=k+1-nx\geq 1$.

We have

$\displaystyle\alpha\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial
c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k+1}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle=$
$\displaystyle\alpha\sum_{i=0}^{k}\left(\left(\frac{k}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial
c}\right)+\left(1-\alpha\right)\sum_{i=0}^{k}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
c}\sum_{i=0}^{k}\left(\left(\alpha\frac{k}{n}+\left(1-\alpha\right)\frac{k+1}{n}-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
c}\sum_{i=0}^{k}\left(\left(x-\frac{i}{n}\right)p_{n,i}^{x,1-x,c}\right)$
$\displaystyle\geq$ $\displaystyle 0$

by Lemma 3.2.

Using (3.11) and the fact that $\alpha\geq 1$, from the above inequality we
conclude that

(3.14)
$\sum_{i=0}^{k+1}\left(\left(\frac{k+1}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)<0.$

We showed that if the sum in (3.11) is strictly negative, then so is the sum
in (3.14). Proceeding inductively on $k$, we obtain

$0\equiv\frac{\partial}{\partial c}\left(1-x\right)=\frac{\partial}{\partial
c}E\left(1-\frac{1}{n}X_{n}^{x,1-x,c}\right)=\sum_{i=0}^{n}\left(\left(\frac{n}{n}-\frac{i}{n}\right)\frac{\partial
p_{n,i}^{x,1-x,c}}{\partial c}\right)<0,$

a contradiction.

In both cases above we showed that (3.11) leads to a contradiction. This shows
that (3.10) holds true, thus concluding the proof of the theorem. ∎

##  4\. An application to the error estimate of Pólya-Stancu operators

As an application of Theorem 3.3, we have the following.

######  Theorem 4.1.

For any convex function $f:\left[0,1\right]\rightarrow\mathbb{R}$ the absolute
value of the error of approximation of the Pólya-Stancu operator $P_{n}^{c}$
is a non-decreasing function of $c$, that is

(4.1)
$\left|P_{n}^{c_{2}}f\left(x\right)-f\left(x\right)\right|\geq\left|P_{n}^{c_{1}}f\left(x\right)-f\left(x\right)\right|,$

for any integer $n\geq 2$, $x\in\left[0,1\right]$, and any
$c_{2}>c_{1}\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$.

If moreover

(4.2) $B_{n}f\left(x\right)\neq B_{1}f\left(x\right)$

for certain values of $n\geq 2$ and $x\in\left[0,1\right]$, then the above
monotonicity is strict, that is

(4.3)
$\left|P_{n}^{c_{2}}f\left(x\right)-f\left(x\right)\right|>\left|P_{n}^{c_{1}}f\left(x\right)-f\left(x\right)\right|,$

for any $c_{2}>c_{1}\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$.

###### Proof.

Jensen’s inequality shows that

$P_{n}^{c}f\left(x\right)=Ef\left(\frac{1}{n}X_{n}^{x,1-x,c}\right)\geq
f\left(\frac{1}{n}EX_{n}^{x,1-x,c}\right)=f\left(x\right),\qquad
c\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\},$

thus the claim (4.1) is equivalent to $P_{n}^{c_{2}}f\left(x\right)\geq
P_{n}^{c_{1}}f\left(x\right)$, and it follows immediately from Theorem 3.3 and
the definition (2.9) of convex ordering.

To prove the second part of the theorem, note that if for certain values
$n\geq 2$ and $x\in\left[0,1\right]$ we have
$P_{n}^{c_{2}}f\left(x\right)=P_{n}^{c_{1}}f\left(x\right)$ for some
$c_{2}>c_{1}\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$, then by the first
part of the proof we have that

(4.4) $P_{n}^{c}f\left(x\right)=C,\qquad c\in\left[c_{1},c_{2}\right],$

is a constant function of $c$ (the constant $C$ may depend on $n$ and $x$).

The above still holds if we replace $f$ by
$f+1+\max_{x\in\left[0,1\right]}f\left(x\right)$ (being convex, $f$ is also
continuous), thus without loss of generality we may assume that
$f\left(x\right)>0$ for $x\in\left[0,1\right]$.

Next note that for $x\in\left(0,1\right)$, from the definition (2.3) of the
rising factorial, it follows that
$x^{\left(k,c\right)}\left(1-x\right)^{\left(n-k,c\right)}$ is a polynomial of
degree $n-2$ in the variable $c$ if $k\in\left\\{1,\ldots,n-1\right\\}$ (with
leading coefficient $x\left(1-x\right)\left(k-1\right)!\left(n-k-1\right)!$),
and a polynomial of degree $n-1$ in the variable $c$ if
$k\in\left\\{0,n\right\\}$ (with leading coefficient
$\left(n-1\right)!\left(1-x\right)$ for $k=0$ and $\left(n-1\right)!x$ for
$k=n$). Using this and the definition (2.7) of the operator $P_{n}^{c}$, it
can be seen that for fixed values $n\geq 2$, $x\in\left[0,1\right]$ and $f$,

$\displaystyle P_{n}^{c}f\left(x\right)$ $\displaystyle=$
$\displaystyle\frac{\sum_{k=0}^{n}\left(f\left(\frac{k}{n}\right)C_{n}^{k}x^{\left(k,c\right)}\left(1-x\right)^{\left(n-k,c\right)}\right)}{1^{\left(n,c\right)}}$
$\displaystyle=$
$\displaystyle\frac{\left(\left(n-1\right)!f\left(0\right)\left(1-x\right)+\left(n-1\right)!f\left(1\right)x\right)c^{n-1}+\ldots+\sum_{k=0}^{n}f\left(\frac{k}{n}\right)x^{k}\left(1-x\right)^{n-k}}{\left(n-1\right)!c^{n-1}+\ldots+1}$

is the ratio of two polynomials of degree $n-1$ in the variable
$c\geq-\frac{1}{n-1}\min\left\\{x,1-x\right\\}$ (recall that $f>0$ and the
special cases $k=0$ and $k=n$ above).

From (4.4) we conclude that these two polynomials (in the variable $c$) are a
constant multiple of each other; in particular this shows that their leading
and free term coefficients are proportional, thus

$\frac{\left(\left(n-1\right)!f\left(0\right)\left(1-x\right)+\left(n-1\right)!f\left(1\right)x\right)}{\left(n-1\right)!}=\frac{\sum_{k=0}^{n}f\left(\frac{k}{n}\right)x^{k}\left(1-x\right)^{n-k}}{1},$

or equivalent

$f\left(0\right)\left(1-x\right)+f\left(1\right)x=\sum_{k=0}^{n}f\left(\frac{k}{n}\right)x^{k}\left(1-x\right)^{n-k}.$

Finally, note that
$\sum_{k=0}^{n}f\left(\frac{k}{n}\right)x^{k}\left(1-x\right)^{n-k}=B_{n}f\left(x\right)$
is just the Bernstein polynomial of degree $n$ corresponding to $f$ evaluated
at $x$, and
$f\left(0\right)\left(1-x\right)+f\left(1\right)x=B_{1}f\left(x\right)$, thus
we have equivalent

(4.5) $B_{n}f\left(x\right)=B_{1}f\left(x\right).$

We have shown that (4.4) implies the condition (4.5). Therefore if (4.5) does
not hold, then (4.4) cannot hold, thus concluding the proof of the theorem. ∎

## References

  * [1] S. N. Bernstein, Démonstration du Théoreme de Weierstrass fondée sur le calcul des Probabilités, Comm. Soc. Math. Kharkov 2 (1912), Series XIII, No.1, pp. 1 – 2.
  * [2] I. Koźniewska, The first absolute central moment for Pólya’s distribution. (Polish) Zastos. Mat. 1 (1954), pp. 206 – 211.
  * [3] M. N. Pascu, N. R. Pascu, F. Tripşa, A new Bernstein-Stancu type operator with negative parameter, Proceedings of the Romanian Academy, Series A, 20 (2019), No. 1, pp. 19 – 28.
  * [4] M. N. Pascu, N. R. Pascu, F. Tripşa, An error estimate for a Bernstein-Stancu operator with negative parameter, Results in Mathematics, 74 (2019), No. 1, Art. 39, 11 pp.
  * [5] M. Shaked, J. G. Shanthikumar, Stochastic orders, Springer Series in Statistics, Springer, New York, 2007.
  * [6] D. D. Stancu, On a new positive linear polynomial operator, Proc. Japan Acad. 44 (1968), pp. 221 – 224.
  * [7] D. D. Stancu, On a generalization of the Bernstein polynomials, Studia Univ. Babeş-Bolyai Ser. Math.-Phys. 14 (1969), No. 2, pp. 31 – 45.
  * [8] D. D. Stancu, Approximation of functions by means of some new classes of positive linear operators, Numerische Methoden der Approximationstheorie, Proc. Conf. Oberwolfach 1971, ISNM 16, pp. 187 – 203, Birkhäuser, Verlag, Basel, 1972.
  * [9] F. Tripşa, N. R. Pascu, Stochastic ordering of Pólya random variables and monotonicity of the Bernstein-Stancu operator for a negative parameter, J. Ineq. Appl. 2019, Paper No. 47, 10 pp.