# Totally non-negativity of a family of change-of-basis matrices

David Galvin and Yufei Zhang Department of Mathematics, University of Notre
Dame, Notre Dame IN; dgalvin1<EMAIL_ADDRESS>Galvin in part supported by a
Simons Collaboration Grant for Mathematicians.

###### Abstract

Let ${\bf a}=(a_{1},a_{2},\ldots,a_{n})$ and ${\bf
e}=(e_{1},e_{2},\ldots,e_{n})$ be real sequences. Denote by $M_{{\bf
e}\rightarrow{\bf a}}$ the $(n+1)\times(n+1)$ matrix whose $(m,k)$ entry
($m,k\in\\{0,\ldots,n\\}$) is the coefficient of the polynomial
$(x-a_{1})\cdots(x-a_{k})$ in the expansion of $(x-e_{1})\cdots(x-e_{m})$ as a
linear combination of the polynomials
$1,x-a_{1},\ldots,(x-a_{1})\cdots(x-a_{m})$. By appropriate choice of ${\bf
a}$ and ${\bf e}$ the matrix $M_{{\bf e}\rightarrow{\bf a}}$ can encode many
familiar doubly-indexed combinatorial sequences, such as binomial
coefficients, Stirling numbers of both kinds, Lah numbers and central
factorial numbers.

In all four of these examples, $M_{{\bf e}\rightarrow{\bf a}}$ enjoys the
property of total non-negativity — the determinants of all its square
submatrices are non-negative. This leads to a natural question: when, in
general, is $M_{{\bf e}\rightarrow{\bf a}}$ totally non-negative?

Galvin and Pacurar found a simple condition on ${\bf e}$ that characterizes
total non-negativity of $M_{{\bf e}\rightarrow{\bf a}}$ when ${\bf a}$ is non-
decreasing. Here we fully extend this result. For arbitrary real sequences
${\bf a}$ and ${\bf e}$, we give a condition that can be checked in $O(n^{2})$
time that determines whether $M_{{\bf e}\rightarrow{\bf a}}$ is totally non-
negative. When $M_{{\bf e}\rightarrow{\bf a}}$ is totally non-negative, we
witness this with a planar network whose weights are non-negative and whose
path matrix is $M_{{\bf e}\rightarrow{\bf a}}$. When it is not, we witness
this with an explicit negative minor.

Keywords: Total non-negativity, matrices, combinatorial sequences, planar
networks, Lindström’s lemma

##  1 Introduction

Let ${\bf a}=(a_{1},a_{2},\ldots,a_{n})$ and ${\bf
e}=(e_{1},e_{2},\ldots,e_{n})$ be real sequences. Each gives rise to a basis
for the space of polynomials in one variable ($x$, say) of degree at most $n$,
in the following way: set

${\mathcal{B}}_{\bf
a}=\\{1,x-a_{1},(x-a_{1})(x-a_{2}),\ldots,\prod_{i=1}^{n}(x-a_{i})\\}$

and define ${\mathcal{B}}_{\bf e}$ analogously.

Associated with ${\mathcal{B}}_{\bf a}$ and ${\mathcal{B}}_{\bf e}$ there is
an $(n+1)\times(n+1)$ change-of-basis matrix $M_{{\bf e}\rightarrow{\bf a}}$,
whose $(m,k)$ entry ($0\leq m\leq n$, $0\leq k\leq n$) is the coefficient of
the polynomial $(x-a_{1})\cdots(x-a_{k})$ in the expansion of
$(x-e_{1})\cdots(x-e_{m})$ as a linear combination of the polynomials in
${\mathcal{B}}_{\bf a}$. (We adopt the standard convention here that the empty
product evaluates to $1$).

Evidently $M_{{\bf e}\rightarrow{\bf a}}$ is lower-triangular, with $1$’s down
the main diagonal. There are numerous combinatorial matrices that are
instances of $M_{{\bf e}\rightarrow{\bf a}}$, for suitable choices of ${\bf
a}$ and ${\bf e}$. We present some examples now (see [10, Section 1] for the
mostly standard justifications).

  1. 1.

Taking $e_{i}=-1$ and $a_{i}=0$ for all $i$ yields $M_{{\bf e}\rightarrow{\bf
a}}=\left[\binom{m}{k}\right]$, the matrix of binomial coefficients.

  2. 2.

Taking $e_{i}=0$ and $a_{i}=i-1$ for all $i$ yields $M_{{\bf e}\rightarrow{\bf
a}}=\left[\genfrac{\\{}{\\}}{0.0pt}{}{m}{k}\right]$, the matrix of Stirling
numbers of the second kind.

  3. 3.

Taking $e_{i}=-(i-1)$ and $a_{i}=0$ for all $i$ yields $M_{{\bf
e}\rightarrow{\bf a}}=\left[\genfrac{\\{}{\\}}{0.0pt}{}{m}{k}\right]$, the
matrix of (unsigned) Stirling numbers of the first kind.

  4. 4.

Taking $e_{i}=-(i-1)$ and $a_{i}=i-1$ for all $i$ yields $M_{{\bf
e}\rightarrow{\bf a}}=\left[L(m,k)\right]$, the matrix of Lah numbers.

  5. 5.

Taking $e_{i}=i-1-b_{i}$ and $a_{i}=i-1$, the $(m,k)$ entry of $M_{{\bf
e}\rightarrow{\bf a}}$ is the Rook number $R_{m-k}(B_{m})$, where
$b_{1},b_{2},\ldots$ is any non-decreasing sequence of non-negative integers,
$B_{m}$ is the Ferrers board with $m$ columns that has $b_{i}$ cells in column
$i$, and $R_{k}(B_{m})$ is the number of ways of placing $k$ non-attacking
rooks on $B_{m}$.

A matrix is totally non-negative (TNN) if all its minors (determinants of
square sub-matrices) are non-negative. Totally non-negative matrices occur
frequently in combinatorics and have been the subject of much investigation.
See for example [5, 6, 9, 11, 15] for an overview. In particular, all of the
matrices listed above are known to be totally non-negative. Motivated by the
first four of the examples above, in [10] Galvin and Pacurar asked the
following question.

######  Question 1.1.

For which pairs of real sequences $({\bf a},{\bf e})$ is $M_{{\bf
e}\rightarrow{\bf a}}$ TNN?

A partial answer was given in [10]. Specifically, a characterization was given
of those ${\bf e}$-sequences for which $M_{{\bf e}\rightarrow{\bf a}}$ is TNN,
when ${\bf a}$ is weakly increasing. We now briefly describe that
characterization, which depends on a notion of restricted growth (introduced,
as far as we are aware, by Arndt [4]).

######  Definition 1.2.

For weakly increasing ${\bf a}$, say that ${\bf e}$ is a restricted growth
sequence relative to ${\bf a}$ if for each $i\geq 1$ it holds that $e_{i}\leq
a_{f(i)}$, where $f(1)=1$ and for $i\geq 1$

$f(i+1)=\left\\{\begin{array}[]{cl}f(i)&\mbox{if $e_{i}<a_{f(i)}$}\\\
f(i)+1&\mbox{if $e_{i}=a_{f(i)}$}.\end{array}\right.$

In other words, each $e_{i}$ is at most a certain cap. The cap for $e_{1}$ is
$a_{1}$. If $e_{1}<a_{1}$ then the cap for $e_{2}$ is also $a_{1}$, while if
$e_{1}=a_{1}$ then the cap for $e_{2}$ is $a_{2}$. In general, the cap for
$e_{i}$ is some $a_{i^{\prime}}$, and if $e_{i}<a_{i^{\prime}}$ then the cap
for $e_{i+1}$ is also $a_{i^{\prime}}$, while if $e_{i}=a_{i^{\prime}}$ then
the cap for $e_{i+1}$ is $a_{i^{\prime}+1}$. If ${\bf
a}=(0,1,\ldots,n-1,\ldots)$ then a non-negative integer sequence ${\bf e}$ is
a restricted growth sequence relative to ${\bf a}$ exactly if it is a
restricted growth sequence in the usual sense, that is, one satisfying
$e_{1}=0$ and $e_{i+1}\leq 1+\max_{j=1,\ldots,i}e_{j}$ for $i\geq 1$.

######  Theorem 1.3.

([10, Theorem 1.2]) Let ${\bf a}$ be a weakly increasing sequence, and let
${\bf e}$ be an arbitrary sequence. Then

  1. 1.

the matrix $M_{{\bf e}\rightarrow{\bf a}}$ is TNN if and only if ${\bf e}$ is
a restricted growth sequence relative to ${\bf a}$, and

  2. 2.

if ${\bf e}$ is not a restricted growth sequence relative to ${\bf a}$ then
the failure of $M_{{\bf e}\rightarrow{\bf a}}$ to be TNN is witnessed by a
negative entry in $M_{{\bf e}\rightarrow{\bf a}}$.

Galvin and Pacurar left open the major part of Question 1.1, that of
determining when $M_{{\bf e}\rightarrow{\bf a}}$ is TNN when no restriction is
placed on ${\bf a}$ (or ${\bf e}$). In this note, we completely resolve this
question. Our resolution involves extending the notion of restricted growth to
the case when ${\bf a}$ is arbitrary, via the following algorithm.

######  Algorithm 1.4.

The input is a pair of real sequences ${\bf a}=(a_{1},\ldots,a_{n})$ and ${\bf
e}=(e_{1},\ldots,e_{n})$.

  1. 1.

Initialize $X=(e_{1},\ldots,e_{n})$ (an ordered list), and $i=1$.

  2. 2.

While $i\leq n$, locate the first element, $e^{\prime}$ say, of $X$ that is
greater than or equal to $a_{i}$.

    1. (a)

If there is no such element, then delete the last element of $X$, and
increment $i$ by $1$.

    2. (b)

If $e^{\prime}=a_{i}$, then delete $e^{\prime}$ from $X$, and increment $i$ by
$1$.

    3. (c)

If $e^{\prime}>a_{i}$, then STOP and report the current value of $i$.

  3. 3.

When $i=n+1$, STOP and report the value $n+1$.

######  Definition 1.5.

For arbitrary real sequences ${\bf a}$ and ${\bf e}$, say that ${\bf e}$ is a
restricted growth sequence relative to ${\bf a}$ if Algorithm 1.4 (on input
${\bf a}$ and ${\bf e}$) terminates with the report $n+1$.

######  Example 1.6.

With ${\bf a}=(3,8,7,5,2,7)$ and ${\bf e}=(2,1,3,7,3,4)$ (so $n=6$), Algorithm
1.4 proceeds as follows:

  * •

Initially $X=(2,1,3,7,3,4)$, $i=1$ and $a_{1}=3$. The first element of $X$
which is greater than or equal to $3$ is the third element, which has value
$3$, so $e^{\prime}=3$.

  * •

Since $e^{\prime}=a_{1}=3$, we update $X$ to $(2,1,7,3,4)$ and set $i=2$.

  * •

Since there is no $e^{\prime}$ (in updated $X$) with $e^{\prime}\geq a_{2}=8$,
we update $X$ to $(2,1,7,3)$ and set $i=3$.

  * •

Now $e^{\prime}=7$, so we update $X$ to $(2,1,3)$ and set $i=4$.

  * •

No element of $X$ is greater than of equal to $a_{4}=5$, so $X$ becomes
$(2,1)$ and $i$ becomes $5$.

  * •

Now $e^{\prime}=2$, $X$ updates to $(1)$ and $i$ to $6$.

  * •

In this final iteration there is no valid $e^{\prime}$ so $X$ updates to the
empty list, $i$ to $7~{}(=n+1)$, and the algorithm terminates with report $7$.
It follows that $(2,1,3,7,3,4)$ is a restricted growth sequence relative to
$(3,8,7,5,2,7)$.

######  Example 1.7.

With ${\bf a}=(11,8,3,1)$ and ${\bf e}=(10,9,2,1)$ (so $n=4$), Algorithm 1.4
proceeds as follows:

  * •

Initially $X=(10,9,2,1)$ and $i=1$; there is no $e^{\prime}$ with
$e^{\prime}\geq a_{1}$, so $X$ updates to $(10,9,2)$.

  * •

At $i=2$ ($a_{2}=8$) we have $e^{\prime}=10>8$, so the algorithm terminates
with report $2~{}(<n+1)$. It follows that $(10,9,2,1)$ is not a restricted
growth sequence relative to $(11,8,3,1)$.

It is easy to check that if ${\bf a}$ is weakly increasing then the notions of
restricted growth given in Definitions 1.2 and 1.5 coincide, so that the
following theorem, the main point of this note, is an extension of item 1 of
Theorem 1.3.

######  Theorem 1.8.

Let ${\bf a}$ and ${\bf e}$ be arbitrary real sequences of length $n$. The
matrix $M_{{\bf e}\rightarrow{\bf a}}$ is TNN if and only if ${\bf e}$ is a
restricted growth sequence relative to ${\bf a}$.

######  Example 1.9.

(Continuing Example 1.6) With ${\bf a}=(3,8,7,5,2,7)$ and ${\bf
e}=(2,1,3,7,3,4)$ we have, by Example 1.6 and Theorem 1.8, that

$M_{{\bf e}\rightarrow{\bf a}}=\left[\begin{array}[]{ccccccc}1&0&0&0&0&0&0\\\
1&1&0&0&0&0&0\\\ 2&8&1&0&0&0&0\\\ 0&42&12&1&0&0&0\\\ 0&42&42&10&1&0&0\\\
0&210&210&62&9&1&0\\\ 0&840&840&272&44&12&1\end{array}\right]$

is TNN (and this may be readily verified directly).

######  Example 1.10.

(Continuing Example 1.7) With ${\bf a}=(11,8,3,1)$ and ${\bf e}=(10,9,2,1)$ we
have that

$M_{{\bf e}\rightarrow{\bf a}}=\left[\begin{array}[]{ccccc}1&0&0&0&0\\\
1&1&0&0&0\\\ 2&0&1&0&0\\\ 18&2&1&1&0\\\ 180&32&4&1&1\end{array}\right]$

is not TNN. This may be readily verified directly, for example by noting that
the minor corresponding to rows 1 and 2, columns 0 and 1 is $-2$ (recall that
our indexing of rows and columns starts at $0$). Note that this is an example
of a non-TNN matrix of the form $M_{{\bf e}\rightarrow{\bf a}}$ where the
failure to be TNN is not witnessed by a negative entry in the matrix; in the
setting of [10], where ${\bf a}$ is weakly increasing, such an example cannot
occur.

Some remarks are in order. Firstly, note that an $(n+1)\times(n+1)$ matrix has
$\Theta\left(4^{n}/\sqrt{n}\right)$ square submatrices and that calculating
the determinant of a $k\times k$ matrix requires $\Omega(k^{2})$ operations.
So the obvious naive algorithm for determining whether such a matrix is TNN
requires $\Omega(n^{3/2}4^{n})$ operations. Ando [3, Corollary 2.2]
established that to test the total non-negativity of a lower-triangular matrix
it suffices to test the $\Theta\left(n2^{n}\right)$ square submatrices whose
columns are consecutive and include the first column. Since a proportion
$1-o(1)$ of these submatrices have dimension $\Theta(n)$, using Ando’s
criterion still requires $\Omega(2^{n})$ operations. Standing in contrast to
this, Theorem 1.8 allows one to check whether $M_{{\bf e}\rightarrow{\bf a}}$
is TNN using $O(n^{2})$ arithmetic operations.

Our next remark is that, unsurprisingly, the proof of Theorem 1.8 involves the
theory of planar networks, and crucially uses the pivot operation of [10]. We
postpone a detailed discussion of this to Section 2. For the moment, let us
note that our proof passes through an intermediate algorithm that takes as
input arbitrary real sequences ${\bf a}$ and ${\bf e}$ of length $n$, and
outputs either

  * •

a planar network with all non-negative weights whose path matrix is $M_{{\bf
e}\rightarrow{\bf a}}$, witnessing that $M_{{\bf e}\rightarrow{\bf a}}$ is TNN

or

  * •

a square submatrix of $M_{{\bf e}\rightarrow{\bf a}}$ (with consecutive rows
and columns) whose determinant is negative, witnessing that $M_{{\bf
e}\rightarrow{\bf a}}$ is not TNN.

(See Algorithm 2.4 and Claim 2.5). We will see that when $M_{{\bf
e}\rightarrow{\bf a}}$ is not TNN, the witnessing negative minor can be
produced in $O(n^{2})$ steps.

Our final remark is that there are various alternate ways of describing the
entries of the matrix $M_{{\bf e}\rightarrow{\bf a}}$; these were all observed
in [10]. Denote by $M^{{\bf a},{\bf e}}(m,k)$ the $(m,k)$ entry of the matrix
$M_{{\bf e}\rightarrow{\bf a}}$.

  1. 1.

The numbers $M^{{\bf a},{\bf e}}(m,k)$ satisfy the recurrence

$M^{{\bf a},{\bf e}}(m,k)=M^{{\bf a},{\bf e}}(m-1,k-1)+(a_{k+1}-e_{m})M^{{\bf
a},{\bf e}}(m-1,k)~{}~{}~{}\mbox{for $m,k>0$}$ (1)

with initial conditions $M^{{\bf a},{\bf e}}(0,0)=1$, $M^{{\bf a},{\bf
e}}(0,k)=0$ for $k>0$ and $M^{{\bf a},{\bf
e}}(m,0)=\prod_{i=1}^{m}(a_{1}-e_{i})$ for $m>0$. Various forms of this
recurrence have appeared in the literature. As observed in [12], with suitable
choices of ${\bf a}$ and ${\bf e}$ the recurrence (1) can encode

     * •

some generalizations of the classical rook numbers [7],

     * •

the normal order coefficients of the word $(VU)^{n}$ in the Weyl algebra
generated by symbols $V,U$ satisfying $UV-VU=hV^{s}$ [7],

     * •

Hsu and Shiue’s generalized Stirling numbers [13],

     * •

the Jacobi-Stirling numbers (coefficients of the Jacobi differential operator)
[7, 8],

as well as encoding Binomial coefficients, Stirling numbers of both kinds, Lah
numbers and rook numbers. Also, with $e_{m}=0$ for all $m$ and
$a_{k}=(k-1)^{2}$ for all $k$, (1) encodes the central factorial numbers of
Riordan and Carlitz (see e.g. [2]).

  2. 2.

We have the explicit expression

$M^{{\bf a},{\bf e}}(m,k)=\sum\prod_{i=1}^{m-k}(a_{s_{i}-i+1}-e_{s_{i}})$

where the sum is over all
$S=\\{s_{1},\ldots,s_{m-k}\\}\subseteq\\{1,\ldots,m\\}$ with
$s_{1}<\ldots<s_{m-k}$.

  3. 3.

We also have the explicit expression

$M^{{\bf a},{\bf
e}}(m,k)=\sum_{\ell=0}^{m-k}(-1)^{\ell}h_{m-k-\ell}(a_{1},\ldots,a_{k+1})s_{\ell}(e_{1},\ldots,e_{m})$

where $h_{\ell}(x_{1},\ldots,x_{t})$ is the degree $\ell$ complete symmetric
polynomial in $x_{1},\ldots,x_{t}$ and $s_{\ell}(x_{1},\ldots,x_{t})$ is the
degree $\ell$ elementary symmetric polynomial in $x_{1},\ldots,x_{t}$.

In Section 2.1 we review the necessary background material from the theory of
planar networks, and describe the pivot operation from [10]. In Section 2.2 we
give the intermediate planar network algorithm. The proof that this algorithm
is correct is given in Section 3.1, while the derivation of Theorem 1.8 from
all that has preceded appears in Section 3.2.

##  2 Planar networks and the intermediate algorithm

###  2.1 Setup

We require some well-known results from the theory of totally non-negative
matrices. (Note that much of what follows is taken from [10, Section 3]).

A planar network $P$ is a directed acyclic planar graph with a subset
$\\{s_{i}:i=0,1,2,\ldots\\}$ of vertices designated as sources and a subset
$\\{t_{i}:i=0,1,2,\ldots\\}$ of vertices designated as sinks. We assume
throughout that the sources and sinks are arranged around the outer face of
the graph, and are ordered cyclically (clockwise) as
$\ldots,s_{2},s_{1},s_{0},t_{0},t_{1},t_{2},\ldots$. A weighted planar network
$(P,w)$ is a planar network $P$ together with a function
$w:\vec{E}(P)\rightarrow{\mathbb{R}}$, which we think of as an assignment of
weights to the edges of $P$.

Figure 1 shows a particular weighted planar network. The vertices set of this
network comprises the sources, the sinks, and the intersection points between
the vertical and horizontal lines. All horizontal lines are oriented to the
right and all vertical lines are oriented upward. All horizontal edges have
weight $1$, while the weights of the vertical edges are given by the
$x_{ij}$’s.

$s_{0}$$s_{1}$$s_{2}$$s_{3}$$\vdots$$s_{n-1}$$s_{n}$$t_{0}$$t_{1}$$t_{2}$$t_{3}$$\vdots$$t_{n-1}$$t_{n}$$\\!x_{11}$$\\!x_{21}$$\\!x_{31}$$\\!x_{n1}$$\\!x_{22}$$\\!x_{32}$$\\!x_{n2}$$\\!x_{33}$$\\!x_{n3}$$\\!x_{n4}$$\\!x_{nn}$$\ddots$
Figure 1: A weighted planar network. The labels on the vertical edges indicate
the value of $w$ for those edges; $w$ takes the value $1$ for all horizontal
edges.

The weighted planar network shown in Figure 1 is the only one that will be
used in our proofs (with various choices of $w$). In figures we will use the
compact array representation shown in Figure 2.

$\begin{array}[]{ccccccc}x_{11}&&&&&&\\\ x_{21}&x_{22}&&&&&\\\
x_{31}&x_{32}&x_{33}&&&&\\\ x_{41}&x_{42}&x_{43}&x_{44}&&&\\\
x_{51}&x_{52}&x_{53}&x_{54}&x_{55}&&\\\
\vdots&\vdots&\vdots&\vdots&&\ddots&\\\
x_{n1}&x_{n2}&x_{n3}&x_{n4}&\cdots&x_{n(n-1)}&x_{nn}.\end{array}$ Figure 2: A
compact array representation of the weighted planar network of Figure 1. We
will use the shorthand $A$ for this array.

We will say that the edge of the planar network that is given weight $x_{mk}$
in Figure 2 is the $[m,k]$ edge of the array (using square brackets to
distinguish this from an entry in a matrix). For each $k=1,\ldots,n$, the
collection of $[m,k]$ edges as $m$ ranges from $k$ to $n$ is the $k$-th column
of the planar network (so, for example, the set of edges with weights
$x_{33},\ldots,x_{n3}$ is the third column of the array in Figure 2).

The path matrix $M(P,w)$ of a generic weighted planar network $(P,w)$ is the
$(n+1)\times(n+1)$ matrix whose $(m,k)$ entry (with rows and columns indexed
by $\\{0,1,2,\ldots,n\\}$) is the sum of the weights of all the directed paths
from $s_{m}$ to $t_{k}$, where the weight of any such path is the product,
over all edges traversed, of the weight of the edge. For example, the $(3,1)$
entry of the path matrix of the weighted planar network shown in Figure 1 is
$x_{31}x_{21}+x_{31}x_{22}+x_{32}x_{22}$. Notice that the path matrix of the
weighted planar network shown in Figure 1 is lower-triangular with $1$’s down
the main diagonal. The following [14] (see also, for example, [1, Chapter 29],
[15]) is a standard result in the theory of totally non-negative matrices.

######  Lemma 2.1.

(Lindström’s Lemma) If $w\geq 0$ (i.e., all weights are non-negative) then the
matrix $M(P,w)$ is TNN.

Moreover, the minor corresponding to selecting the rows indexed by $I$ and
columns indexed by $J$ (with indexing of rows and columns starting from $0$)
equals the sum of the weights of all the collections of $|I|$ vertex disjoint
paths from the sources $\\{s_{i}:i\in I\\}$ to the sinks $\\{t_{j}:j\in J\\}$,
where the weight of a collection of paths is the product of the weights of the
individual paths in the collection.

As mentioned earlier, all the weighted planar networks that appear in our
proofs have the form shown in Figure 1. Such networks are completely specified
by the array of weights shown in Figure 2. Denoting that array by $A$, we use
the shorthand $M(A)$ for $M(P,w)$ (where $(P,w)$ is the weighted planar
network shown in Figure 1). Abusing notation slightly, we will refer to $M(A)$
as the path matrix of the array $A$.

All of the arrays that appear in our proofs will be generated from real
sequences ${\bf a}=(a_{1},\ldots,a_{n})$ and ${\bf e}=(e_{1},\ldots,e_{n})$.
Figure 3 shows the generic array $A({\bf a},{\bf e})$ generated from ${\bf a}$
and ${\bf e}$, while Figure 4 shows a specific instance.

$\begin{array}[]{ccccccc}a_{1}-e_{1}&&&&&&\\\ a_{1}-e_{2}&a_{2}-e_{1}&&&&&\\\
a_{1}-e_{3}&a_{2}-e_{2}&a_{3}-e_{1}&&&&\\\
a_{1}-e_{4}&a_{2}-e_{3}&a_{3}-e_{2}&a_{4}-e_{1}&&&\\\
a_{1}-e_{5}&a_{2}-e_{4}&a_{3}-e_{3}&a_{4}-e_{2}&a_{5}-e_{1}&&\\\
\vdots&\vdots&\vdots&\vdots&&\ddots&\\\
a_{1}-e_{n}&a_{2}-e_{n-1}&a_{3}-e_{n-2}&a_{4}-e_{n-3}&\cdots&a_{n-1}-e_{2}&a_{n}-e_{1}\end{array}$
Figure 3: The array $A({\bf a},{\bf e})$.
$\begin{array}[]{clc}\begin{array}[]{ccccccc}3-2&&&&&&\\\ 3-1&8-2&&&&&\\\
3-3&8-1&7-2&&&&\\\ 3-7&8-3&7-1&5-2&&&\\\ 3-3&8-7&7-3&5-1&2-2&&\\\
3-4&8-3&7-7&5-3&2-1&7-2\end{array}&\mbox{or}{}{}{}&\begin{array}[]{rrrrrrr}1&&&&&&\\\
2&6&&&&&\\\ 0&7&5&&&&\\\ -4&5&6&3&&&\\\ 0&1&4&4&0&&\\\
-1&5&0&2&1&5\end{array}\end{array}$ Figure 4: The array $A({\bf a},{\bf e})$
in the specific case ${\bf a}=(3,8,7,5,2,7)$, ${\bf e}=(2,1,3,7,3,4)$ (Example
1.6).

We need the following lemma ([10, Lemma 3.5]):

######  Lemma 2.2.

For arbitrary ${\bf a}$ and ${\bf e}$, we have $M(A({\bf a},{\bf e}))=M_{{\bf
e}\rightarrow{\bf a}}$.

If it happens that all weights in $A({\bf a},{\bf e})$ are non-negative, then
the total non-negativity of $M_{{\bf e}\rightarrow{\bf a}}$ follows
immediately from Lemmas 2.1 and Lemma 2.2. If $A({\bf a},{\bf e})$ has some
negative weights, then we will give a procedure that modifies the array by
iteratively replacing negative weights with $0$ weights, while not changing
the path matrix of the array. We will show that this procedure succeeds in
replacing all negative weights with $0$ weights exactly when $M_{{\bf
e}\rightarrow{\bf a}}$ is totally non-negative. To describe the procedure we
need ([10, Lemma 3.9]), whose statement requires a little more notation. For
an array $A$ of the kind shown in Figure 2, we refer to the triangle
consisting of the $[m+\ell_{1},k+\ell_{2}]$ edges for $0\leq\ell_{1}\leq n-m$
and $0\leq\ell_{2}\leq\ell_{1}$ as the triangle headed at the $[m,k]$ edge
(see Figure 5).

$\begin{array}[]{ccccccccc}&[1,1]&&&&&&&\\\ &[2,1]&[2,2]&&&&&&\\\
&[3,1]&{\bf[3,2]}&[3,3]&&&&&\\\ &[4,1]&{\bf[4,2]}&{\bf[4,3]}&[4,4]&&&&\\\
&[5,1]&{\bf[5,2]}&{\bf[5,3]}&{\bf[5,4]}&[5,5]&&&\\\
&\vdots&\vdots&&&\ddots&\ddots&&\\\
&[n,1]&{\bf[n,2]}&{\bf[n,3]}&{\bf[n,4]}&{\bf\cdots}&{\bf[n,n-1]}&[n,n]&\end{array}$
Figure 5: The triangle headed at the $[3,2]$ position (bolded entries).

Suppose that the triangle headed at the $[m,k]$ edge of $A$ is an instance of
$A({\bf f},{\bf g})$ for some sequences ${\bf f}=(f_{1},\ldots,f_{n-m+1})$,
${\bf g}=(g_{1},\ldots,g_{n-m+1})$. Denote by $A^{\rm p}(m,k)$ the array
obtained from $A$ by the following operations:

  * •

Leave unchanged all weights of $A$ that are on edges not in the the triangle
headed at the $[m,k]$ edge.

  * •

Also leave unchanged the weight ($f_{1}-g_{1}$) on the $[m,k]$ edge. Replace
the weights on the $[m+1,k]$ and $[m+1,k+1]$ edges ($f_{1}-g_{2}$ and
$f_{2}-g_{1}$, respectively), with $f_{1}-g_{1}$ and $f_{2}-g_{2}$,
respectively.

  * •

More generally, replace the weights on the
$[m+\ell,k],[m+\ell,k+1],[m+\ell,k+2],\ldots,[m+\ell,k+\ell]$ edges
($f_{1}-g_{\ell+1},f_{2}-g_{\ell},\ldots,f_{\ell+1}-g_{1}$, respectively),
with $f_{1}-g_{1},f_{2}--g_{\ell+1},\ldots,f_{\ell+1}-g_{2}$, respectively,
for $1\leq\ell\leq n-m$.

We refer to $A^{\rm p}(m,k)$ as the array obtained from $A$ by pivoting on the
$[m,k]$ edge. Note that in each row of the triangle headed at the $[m,k]$
edge, $A^{\rm p}({\bf f},{\bf g})$ is obtained from $A$ by moving the
$-g_{1}$’s from the last position in the row to the first, and then shifting
all other $-g_{j}$’s in the row one place to the right. Note also that if
$f_{1}=g_{1}$, then after pivoting the weights on the $[m^{\prime},k]$ edges,
for all $m^{\prime}\geq m$, are $0$. (See Figure 6 for an illustration.)

$\begin{array}[]{clc}\begin{array}[]{ccccccc}3-2&&&&&&\\\ 3-1&8-2&&&&&\\\
\underline{{\bf 3-3}}&8-1&7-2&&&&\\\
\underline{3-3}&\underline{8-7}&7-1&5-2&&&\\\
\underline{3-3}&\underline{8-3}&\underline{7-7}&5-1&2-2&&\\\
\underline{3-3}&\underline{8-4}&\underline{7-3}&\underline{5-7}&2-1&7-2\end{array}&\mbox{or}{}{}{}&\begin{array}[]{rrrrrrr}1&&&&&&\\\
2&6&&&&&\\\ 0&7&5&&&&\\\ 0&1&6&3&&&\\\ 0&5&0&4&0&&\\\
0&4&4&-2&1&5\end{array}\end{array}$ Figure 6: The array from Figure 4 (Example
1.6), after pivoting on the $[3,1]$ edge (bolded). The edges of the triangle
headed at the $[3,1]$ edge (the only edges whose weights might change by
pivoting) are underlined.
$\begin{array}[]{clc}\begin{array}[]{ccccccc}3-2&&&&&&\\\ 3-1&8-2&&&&&\\\
3-3&8-1&7-2&&&&\\\ 3-3&8-7&7-1&5-2&&&\\\ 3-3&8-3&\underline{{\bf
7-7}}&5-1&2-2&&\\\
3-3&8-4&\underline{7-7}&\underline{5-3}&2-1&7-2\end{array}&\mbox{or}{}{}{}&\begin{array}[]{ccccccc}1&&&&&&\\\
2&6&&&&&\\\ 0&7&5&&&&\\\ 0&1&6&3&&&\\\ 0&5&0&4&0&&\\\
0&4&0&2&1&5\end{array}\end{array}$ Figure 7: The array from Figure 6 (Example
1.6), after pivoting on the $0$ at the $[5,3]$ position (bolded).

######  Lemma 2.3.

([10, Lemma 3.9]) For an array $A$ of the kind shown in Figure 2, suppose that
the triangle headed at the $[m,k]$ edge of $A$ is an instance of $A({\bf
f},{\bf g})$ for some sequence ${\bf f}=(f_{1},\ldots,f_{n-m+1})$, ${\bf
g}=(g_{1},\ldots,g_{n-m+1})$.

If $f_{1}=g_{1}$ then $M(A^{\rm p}(m,k))=M(A)$ (that is, the path matrix is
unchanged after pivoting on the $[m,k]$ edge, as long as the entry on the
$[m,k]$ edge is $0$).

It should be noted that when Lemma 2.3 was stated in [10], some assumptions
were made about the weights in $A$ on edges other than those edges in the
triangle headed at the $[m,k]$ edge. However, those assumptions were not
actually used in the proof, so the proof goes through unchanged in the setting
of Lemma 2.3.

###  2.2 The intermediate algorithm

We are now in a position to describe the intermediate planar network algorithm
alluded to in the introduction.

######  Algorithm 2.4.

The input is a pair of real sequences ${\bf a}=(a_{1},\ldots,a_{n})$ and ${\bf
e}=(e_{1},\ldots,e_{n})$. The algorithm proceeds by modifying the weighted
planar network $A({\bf a},{\bf e})$ (see Figure 3).

  * •

Initially, set $A=A({\bf a},{\bf e})$ and $k=1$.

  * •

As long as $k\leq n$, scan the weights in the $k$th column of $A$, from the
$[k,k]$ edge down.

    * –

If no zero or negative weights are encountered in the column, leave $A$
unchanged and increase $k$ by $1$.

    * –

If the first non-positive weight that is encountered is negative then STOP,
output the current $A$ and report that $M_{{\bf e}\rightarrow{\bf a}}$ is not
TNN.

    * –

If the first non-positive weight that is encountered is zero, say on the
$[m,k]$ edge of the array $A$:

      * *

if the triangle headed at the $[m,k]$ edge of $A$ is an instance of $A({\bf
f},{\bf g})$ for some sequences ${\bf f}=(f_{1},\ldots,f_{n-m+1})$, ${\bf
g}=(g_{1},\ldots,g_{n-m+1})$, then update $A$ by pivoting on the $[m,k]$
position — that is, replace $A$ by $A^{\rm p}(m,k)$ — and increase $k$ by $1$.
(See Figure 3 for an illustration of an instance of $A({\bf f},{\bf g})$.)

      * *

if not, then STOP and report that the algorithm is inconclusive. (We will see
in the analysis that this clause of the algorithm will never be invoked).

  * •

If $k=n+1$, then STOP, output the current $A$ and report that $M_{{\bf
e}\rightarrow{\bf a}}$ is TNN.

For the input ${\bf a}=(3,8,7,5,2,7)$, ${\bf e}=(2,1,3,7,3,4)$ (Example 1.6),
the evolution of the planar network in Algorithm 2.4 is illustrated in Figures
4, 6 and 7. In this case the algorithm reports that $M_{{\bf e}\rightarrow{\bf
a}}$ is TNN, and outputs the final network that is illustrated in Figure 7.
For the input ${\bf a}=(11,8,3,1)$ and ${\bf e}=(10,9,2,1)$ (Example 1.7),
Algorithm 2.4 terminates at $k=2$ (without having done a pivot), and in this
case the algorithm reports that $M_{{\bf e}\rightarrow{\bf a}}$ is not TNN.
The final network that the algorithm outputs is illustrated in Figure 8.

$\begin{array}[]{clc}\begin{array}[]{cccc}11-10&&&\\\ 11-9&8-10&&\\\
11-2&8-9&3-10&\\\
11-1&8-2&3-9&1-10\end{array}&\mbox{or}{}{}{}&\begin{array}[]{rrrr}1&&&\\\
2&-2&&\\\ 9&-1&-7&\\\ 10&7&-6&-9\end{array}\end{array}$ Figure 8: The final
network in Algorithm 2.4 when ${\bf a}=(11,8,3,1)$ and ${\bf e}=(10,9,2,1)$
(Example 1.7).

######  Claim 2.5.

For every input ${\bf a}$, ${\bf e}$:

  1. 1.

Algorithm 2.4 terminates after at most $O(n^{3})$ operations.

  2. 2.

It never terminates with an inconclusive report.

  3. 3.

If it terminates with the report that $M_{{\bf e}\rightarrow{\bf a}}$ is TNN,
then the outputted weighted planar network $A$ has all non-negative entries
and has path matrix $M_{{\bf e}\rightarrow{\bf a}}$ (and so serves as a
witness that $M_{{\bf e}\rightarrow{\bf a}}$ is TNN).

  4. 4.

If it terminates with the report that $M_{{\bf e}\rightarrow{\bf a}}$ is not
TNN, then there is an $O(n^{2})$ algorithm that identifies, in the outputted
weighted planar network $A$, a collection $S$ of consecutive sources and a
collection $T$ of consecutive sinks with $|S|=|T|$, such that all collections
of vertex disjoint paths from $S$ to $N$ have non-positive weight, and at
least one has negative weight (which serves, via Lemma 2.1, to witness that
$M_{{\bf e}\rightarrow{\bf a}}$ is not TNN).

##  3 Proofs of Claim 2.5 and Theorem 1.8

###  3.1 Proof of Claim 2.5

####  3.1.1 Item 1

Constructing the initial weighted planar network $A$ requires $O(n^{2})$
operations. The algorithm requires scanning $O(n^{2})$ entries. At most $O(n)$
of these scans require performing a pivot, and each pivot requires $O(n^{2})$
recalculations of weights. So in total, the algorithm requires $O(n^{3})$
operations.

####  3.1.2 Item 2

This is an immediate corollary of the following straightforward observation:

######  Observation 3.1.

Suppose that at some point in the algorithm the triangle headed at the $[k,k]$
edge is of the form $A({\bf f},{\bf g})$ for some sequences ${\bf
f}=(f_{k},\ldots,f_{n}),{\bf g}=(g_{1},\ldots,g_{n-k+1})$. If we then pivot on
the $[m,k]$ edge for some $m\geq k$ (so necessarily $f_{k}=g_{m-k+1}$), then
after the pivot

  * •

all the entries in positions $[m^{\prime},k]$ for $m^{\prime}\geq m$ are $0$
(more precisely, they are all $f_{k}-g_{m-k+1}$), and

  * •

the triangle headed at the $[k+1,k+1]$ position is of the form $A({\bf
f}^{\prime},{\bf g}^{\prime})$ where ${\bf f}^{\prime}=(f_{k+1},\ldots,f_{n})$
and ${\bf
g}^{\prime}=(g_{1},\ldots,g_{m-k},\widehat{g_{m-k+1}},g_{m-k+2},\ldots,g_{n-k+1})$.
(The $\widehat{g_{m-k+1}}$ here indicates the removal of that entry from the
sequence).

See Figures 6 and 7 for illustrations of Observation 3.1.

####  3.1.3 Item 3

Here Observation 3.1 is also a key ingredient. The proof of the statement in
item 3 is by induction (on the column $k$ being scanned).

By Lemma 2.2, the initial network has path matrix $M_{{\bf e}\rightarrow{\bf
a}}$, and the triangle headed at the $[1,1]$ edge is of the form $A({\bf
f},{\bf g})$ (specifically with ${\bf f}={\bf a}$ and ${\bf g}={\bf e}$).

Suppose that we are at the point in the algorithm where we are scanning the
$k$th column. Let us assume (by induction) that in the current weighted planar
network $A$, the triangle headed at the $[k,k]$ edge is of the form $A({\bf
f},{\bf g})$, and that all the entries in the first $k-1$ columns are non-
negative.

If column $k$ has only positive entries, then we move on to column $k+1$
without changing $A$; so as we scan column $k+1$, in the current weighted
planar network $A$, the triangle headed at the $[k+1,k+1]$ edge is of the form
$A({\bf f}^{\prime},{\bf g}^{\prime})$ (with ${\bf f}^{\prime}$ obtained from
${\bf f}$ by removing the first element of ${\bf f}$, and ${\bf g}^{\prime}$
obtained from ${\bf g}$ by removing the last element of ${\bf g}$). Also, all
the entries in the first $k$ columns are non-negative.

If column $k$ has its first zero at position $[m,k]$, for some $m\geq k$,
then, by the observation above, after pivoting the $k$th column has all its
entries non-negative, and (in the updated network) the triangle headed at the
$[k+1,k+1]$ edge is of the form $A({\bf f}^{\prime},{\bf g}^{\prime})$. Also
(crucially), by Lemma 2.3, the new network still has path matrix $M_{{\bf
e}\rightarrow{\bf a}}$.

So by induction (on $k$), when we get to the point where $k=n+1$ (which we do
reach, by assumption), the final weighted planar network $A$ has all non-
negative entries and has path matrix $M_{{\bf e}\rightarrow{\bf a}}$.

####  3.1.4 Item 4

Suppose that the first negative weight that is encountered in running the
algorithm is at the $[m,k_{0}]$ edge.

Mark the $[m,k_{0}]$ edge. Consider the path $p_{m}$, from source $s_{m}$ to
sink $t_{k}$, where $k=(k_{0})=k_{0}-1$, constructed as follows: start by
going horizontally from $s_{m}$ to the lower vertex of the marked edge, then
go vertically until the path currently being constructed intersects with the
unique path from $s_{k_{0}-1}$ to $t_{k_{0}-1}$, and then go horizontally to
$t_{k_{0}-1}$. This path only uses edges with non-zero weight (the
$[m^{\prime},k_{0}]$ edges, for $m^{\prime}<m$, have positive weight), and in
fact uses only edges with positive weight, except for the $[m,k_{0}]$ edge,
which has negative weight.

Next, let $k_{1}$ be maximal subject to the conditions that firstly
$k_{1}<k_{0}$, and secondly that the $[m-1,k_{1}]$ edge has positive weight
(assuming such a $k_{1}$ exists — if not, the process terminates, as described
in more detail in a moment), and mark the $[m-1,k_{1}]$ edge. Note that there
is a path, $p_{m-1}$, from source $s_{m-1}$ to sink $t_{k-1}$ — horizontal to
the marked edge, then vertical to the unique path from $s_{k-1}$ to $t_{k-1}$
(the vertical stretch of the network goes at least as high as $k-1$, since it
is strictly to the left of a vertical stretch that goes at least as far as the
unique path from $s_{k}$ to $t_{k}$), then horizontal. This path is vertex
disjoint from $p_{k}$. Also, this path only uses edges with positive weight —
that the marked edge has positive weight forces all edges above it to have
positive weight.

Repeat — mark (if it exists) the $[m-2,k_{2}]$ edge, where $k_{2}<k_{1}$ is
maximal subject to the that edge having positive weight. Note that there is a
path, $p_{m-2}$ say, from source $s_{m-2}$ to sink $t_{k-2}$, that is vertex-
disjoint from $p_{m},p_{m-1}$, and again only uses edges with positive weight.

Repeat, until no markable edge is found. Suppose that the last marked edge is
the $[m-\ell,k_{\ell}]$ edge, where $\ell\geq 0$ (so $\ell+1$ edges have been
marked).

There is a collection of vertex disjoint paths in the network $A$ from
$\\{s_{m-\ell},\ldots,s_{m}\\}$ to $\\{t_{k-\ell},\ldots,t_{k}\\}$ — namely
$\\{p_{m-\ell},\ldots,p_{m}\\}$ — and it has negative weight, since all edges
on it have positive weight, except the first marked edge at position
$[m,k_{0}]$, which is negative. Note also that the construction of this
collection can be done by looking at the weight on each edge of the network
only once, and so in time $O(n^{2})$.

There may be more such collections of vertex disjoint paths from
$\\{s_{m-\ell},\ldots,s_{m}\\}$ to $\\{t_{k-\ell},\ldots,t_{k}\\}$, that avoid
edges with weight zero, but we assert that all of these will use the
$[m,k_{0}]$ edge and no other edge with negative weight. If this were true,
then all such collections of paths would have negative weight, and so the sum
of the weights of all the collections of vertex disjoint paths from
$\\{s_{m-\ell},\ldots,s_{m}\\}$ to $\\{t_{k-\ell},\ldots,t_{k}\\}$ would be
negative. From this, it would follow (by Lindström’s lemma, Lemma 2.1) that
the minor of the path matrix of $A$ corresponding to rows $m-\ell$ through
$m$, columns $k-\ell$ through $k$, would be negative. Since that path matrix
is $M_{{\bf e}\rightarrow{\bf a}}$ (this follows from the discussion of item
3), item 4 would follow.

######  Example 3.2.

For the final network that comes from Algorithm 2.4 with the input ${\bf
a}=(11,8,3,1)$ and ${\bf e}=(10,9,2,1)$ (Example 1.7, see Figure 8), we mark
the $[2,2]$ and $[1,1]$ edges, and the identified set of sources is
$\\{2,1\\}$ while the sinks are $\\{0,1\\}$; this recovers the negative minor
pointed out in the introduction.

######  Example 3.3.

For a more substantial example, consider the weighted planar network whose
sign-pattern in shown in Figure 9.

$\begin{array}[]{ccccccc}+&&&&&&\\\ +&+&&&&&\\\ +&+&+&&&&\\\ +&+&+&+&&&\\\
+&0&+&+&+&&\\\ +&0&+&+&0&+&\\\ +&0&+&0&0&-&\star\end{array}$ Figure 9: An
example of the possible sign pattern on a final planar network from Algorithm
2.4, when the algorithm terminates with the report that $M_{{\bf
e}\rightarrow{\bf a}}$ in not TNN.

In this case the marked edges are those at positions $[7,6]$, $[6,4]$,
$[5,3]$, $[4,2]$ and $[3,1]$, the identified set of sources is
$\\{3,4,5,6,7\\}$ and the identified set of sinks is $\\{1,2,3,4,5\\}$. There
is more than one collection of vertex disjoint paths from these sources to
these sinks that avoids edges with weight $0$, but all of them use exactly one
edge with negative weight (the one at position $[7,6]$), so all have negative
weight.

So it remains to prove the assertion that all collections of vertex disjoint
path systems from $\\{s_{m},\ldots,s_{m-\ell}\\}$ to
$\\{t_{k},\ldots,t_{k-\ell}\\}$, that avoid edges with weight zero, use the
$[m,k_{0}]$ edge and no other edge with negative weight.

For the proof, first observe that any collection of vertex disjoint paths from
$\\{s_{m},\ldots,s_{m-\ell}\\}$ to $\\{t_{k},\ldots,t_{k-\ell}\\}$ that avoids
edges with weight zero, and that uses $p_{m}$, must have exactly one edge with
negative weight (there are no edges with negative weight above or to the left
of the the $[m,k]$ edge).

Suppose we have a path $p^{\prime}_{m}$ from $s_{m}$ to $t_{k}$ that avoids
edges with weight zero and that makes its first vertical turn before the
$[m,k_{0}]$ edge, i.e, before column $k_{0}$, say at column
$k_{0}^{\prime}<k_{0}$. Necessarily, $k_{0}^{\prime}\leq k_{1}$. So any path
from $s_{m-1}$ to $t_{k-1}$ that avoids edges with weight zero and that is
disjoint from $p_{m}^{\prime}$ must make its first vertical turn earlier than
$p_{m-1}$ does. Let $p_{m-1}^{\prime}$ be such a path, turning vertical first
at column $k_{1}^{\prime}<k_{1}$. We have $k_{1}^{\prime}\leq k_{2}$, so any
path from $s_{m-2}$ to $t_{k-2}$ that avoids edges with weight zero and that
is disjoint from $p_{m-1}^{\prime}$ must make its first vertical turn earlier
than $p_{m-2}$ does. Iterating, we eventually get that if
$p_{m}^{\prime},\ldots,p_{m-\ell+1}^{\prime}$ is any vertex disjoint path
system from $\\{s_{m},\ldots,s_{m-\ell+1}\\}$ to
$\\{t_{k},\ldots,t_{k-\ell+1}\\}$ that avoids edges with weight zero, then it
must be the case that $p_{m-\ell+1}^{\prime}$ makes its first vertical turn
earlier than $p_{m-\ell+1}$ does. So any path from $s_{k-\ell}$ to
$t_{k-\ell}$ that is disjoint from $p_{m-\ell+1}^{\prime}$ and avoids edges
with weight zero must make its first vertical turn before it reaches the
$[m-\ell,k_{\ell}]$ edge. But no such vertical turn exists — if it did, the
process of marking edges would not have stopped with the $[m-\ell,k_{\ell}]$
edge. This contradiction finishes the proof of the assertion that all the
vertex disjoint path systems from $\\{s_{m},\ldots,s_{m-\ell}\\}$ to
$\\{t_{k},\ldots,t_{k-\ell}\\}$, that avoid edges with weight zero, use the
$[m,k_{0}]$ edge and no other negative edges, and finishes the verification of
item 4.

###  3.2 Deriving Theorem 1.8 from Claim 2.5

First suppose that ${\bf e}$ is a restricted growth sequence relative to ${\bf
a}$. Then it is easy to see that Algorithm 2.4 terminates with the report that
$M_{{\bf e}\rightarrow{\bf a}}$ is TNN, and moreover that we can recover the
final weighted planar network $A$ of Algorithm 2.4, by the following process:
for $i=1,\ldots,n$, in the $i$th iteration of Algorithm 1.4,

  * •

if no $e^{\prime}\in X$ is found with $e^{\prime}=a_{i}$ (so all
$e^{\prime}<a_{i}$), then for each $m=i,\ldots,n$, at the $[m,i]$ edge of the
weighted planar network put the weight $a_{i}-(X)_{m-i+1}$. Here $(X)_{j}$
indicates the $j$th entry of $X$ (in its $i$th incarnation); and

  * •

if there is an $e^{\prime}\in X$ with $e^{\prime}=a_{i}$, say
$e^{\prime}=(X)_{m^{\prime}-i+1}$ for some $m^{\prime}$, then for each
$m=i,\ldots,m^{\prime}$, at the $[m,i]$ edge of the weighted planar network
put the weight $a_{i}-(X)_{m-i+1}$, and for each $m=m^{\prime}+1,\ldots,n$, at
the $[m,i]$ edge put the weight $0$.

Note moreover that using this process, the weighted planar network that
witnesses that $M_{{\bf e}\rightarrow{\bf a}}$ is TNN can be constructed from
${\bf a}$ and ${\bf e}$ in only $O(n^{2})$ operations.

On the other hand, suppose that ${\bf e}$ is not a restricted growth sequence
relative to ${\bf a}$. If we execute the procedure described above (that
associates a weighted planar network with a running of Algorithm 1.4),
stopping it at the point where an $e^{\prime}>a_{i}$ is first found, then the
resulting partial weighted planar network agrees with the final network of
Algorithm 2.4 at all the edges that are involved in the collections of vertex
disjoint paths that witness a negative minor of $M_{{\bf e}\rightarrow{\bf
a}}$, and so the collection of sources and sinks (rows and columns of $M_{{\bf
e}\rightarrow{\bf a}}$) that witnesses a negative minor of $M_{{\bf
e}\rightarrow{\bf a}}$ can be recovered from the output of Algorithm 1.4.

## Acknowledgement

We thank an anonymous referee for a very careful reading and many helpful
suggestions to improve the clarity of the presentation.

## References

  * [1] Martin Aigner and Günter M. Ziegler, Proofs from THE BOOK (4th edition), Springer Verlag, 2010.
  * [2] Feryal Alayont and Nicholas Krzywonos, Rook polynomials in three and higher dimensions, Involve 6 (2013), 35–52.
  * [3] T. Ando, Totally positive matrices, Linear Algebra Appl. 90 (1987), 165–219.
  * [4] Jörg Arndt, Matters Computational: Ideas, Algorithms, Source Code, Springer, 2011.
  * [5] Francesco Brenti, Combinatorics and total positivity, J. Combin. Theory Ser. A 71 (1995), 175–218.
  * [6] F. Brenti, The applications of total positivity to combinatorics, and conversely, in Total Positivity and its Applications, (M. Gasca, C. Micchelli, eds.), Kluwer Academic Pub., Dordrecht, The Netherlands, 1996, 451–473.
  * [7] Richell O. Celeste, Roberto B. Corcino and Ken Joffaniel M. Gonzales, Two approaches to normal order coefficients, J. Integer Seq. 20 (2017), article 17.3.5.
  * [8] W. N. Everitt, K. H. Kwon, L. L. Littlejohn, R. Wellman, G. J. Yoon, Jacobi-Stirling numbers, Jacobi polynomials, and the left-definite analysis of the classical Jacobi differential expression, J. Comput. Appl. Math. 208 (2007), 29–56.
  * [9] Sergey Fomin and Andrei Zelevinsky, Total positivity: tests and parameterizations, Math. Intelligencer 22 (2000), 23–33.
  * [10] David Galvin and Adrian Pacurar, Total non-negativity of some combinatorial matrices, J. Combin. Th. Ser. A 172 (2020), article 105179.
  * [11] Mariano Gasca and Charles A. Micchelli (eds.), Total Positivity and its Applications, Kluwer Academic Pub., Dordrecht, The Netherlands, 1996.
  * [12] Ken Joffaniel M. Gonzales, Total positivity of a class of combinatorial matrices, unpublished manuscript.
  * [13] Leetsch C. Hsu and Peter Jau-Shyong Shiue, A unified approach to generalized Stirling numbers, Adv. Appl. Math. 20 (1998), 366–384.
  * [14] Bernt Lindström, On the vector representations of induced matroids, Bull. Lond. Math. Soc. 5 (1973), 85–90.
  * [15] Mark Skandera, Introductory notes on total positivity, retrieved from http://www.math.lsa.umich.edu/~fomin/eidma.html, June 1 2023.


# A Low-Complexity Machine Learning Design for mmWave Beam Prediction

Muhammad Qurratulain Khan , Abdo Gaber , Mohammad Parvini , , Philipp Schulz ,
and Gerhard Fettweis This work was funded by the European Union’s SEMANTIC ITN
project under the Marie Skodowska-Curie grant agreement No. 861165.Muhammad
Qurratulain Khan and Abdo Gaber are with the National Instruments Corporation,
Dresden, Germany. (email: {muhammad.qurratulain.khan,
abdo.nasser.ali.gaber}@ni.com).Mohammad Parvini, Philipp Schulz, and Gerhard
Fettweis are with the Vodafone Chair for Mobile Communications Systems,
Technische Universität Dresden, Germany. (email: {mohammad.parvini,
philipp.schulz2, gerhard.fettweis}@tu-dresden.de). 0000-0001-5779-0101
0000-0001-5779-0101 0000-0002-1315-7635 0000-0002-0738-556X
0000-0003-4622-1311

###### Abstract

The 3rd Generation Partnership Project (3GPP) is currently studying machine
learning (ML) for the fifth generation (5G)-Advanced New Radio (NR) air
interface, where spatial and temporal-domain beam prediction are important use
cases. With this background, this letter presents a low-complexity ML design
that expedites the spatial-domain beam prediction to reduce the power
consumption and the reference signaling overhead, which are currently
imperative for frequent beam measurements. Complexity analysis and evaluation
results showcase that the proposed model achieves state-of-the-art accuracy
with lower computational complexity, resulting in reduced power consumption
and faster beam prediction. Furthermore, important observations on the
generalization of the proposed model are presented in this letter.

###### Index Terms:

beam prediction, machine learning (ML), millimeter-wave (mmWave), supervised
learning (SL).

mmWave

    millimeter-wave
5G

    fifth generation
THz

    terahertz
BS

    base station
UE

    user equipment
3GPP

    3rd Generation Partnership Project
SSB

    synchronization sequence block
RSRP

    reference signal received power
EBS

    exhaustive beam scan
HBS

    hierarchical beam scan
AI

    artificial intelligence
ML

    machine learning
FC

    fully-connected
NN

    neural network
DNN

    deep neural network
CNN

    convolutional neural network
NR

    New Radio
MIMO

    multi-in multi-out
LOS

    line-of-sight
NLOS

    non-line-of-sight
AoA

    angle-of-arrival
AoD

    angle-of-departure
UPA

    uniform planar array
RF

    radio frequency
AWGN

    additive white Gaussian noise
KPI

    key performance indicator
CDL

    clustered delay line
FLOP

    floating-point operation

##  I Introduction

The availability of abundant bandwidth at millimeter-wave (mmWave) bands makes
it a requisite for higher throughput. However, to achieve an adequate link
margin, beamforming via large antenna arrays is essential [1]. Consequently,
the evaluation of beam qualities through frequent beam measurements and beam
qualities reporting is imperative to help the base station (BS) and the user
equipment (UE) decide the optimal beam pair for link establishment. Within the
3GPP this is referred to as beam management procedure.

In order to enable the UE to measure the beam qualities, beamformed reference
signals (synchronization sequence blocks) are sequentially transmitted from
the BS in the form of an SSB burst. This allows the UE to measure the
qualities of all the BS transmit beams in terms of their reference signal
received powers through one of its receive beams. Further, to measure the
qualities of all possible transmit-receive beam pairs, several SSB bursts are
transmitted. This procedure of beam qualities measurement is known as
exhaustive beam scan (EBS), which suffers from large beam measurement
overhead, increased latency, and higher power consumption [2, 3]. To overcome
this, a two-level hierarchical beam scan (HBS) consisting of parent (wide) and
child (narrow) beams is employed [4]. Nevertheless, it suffers from increased
latency and inaccuracy of beam selection.

Recently, ML methods have been extensively applied to wireless communications
to solve the non-linear problems that were burdensome to be resolved by
conventional signal processing techniques. Consequently, several studies
propose the use of ML for beam prediction and selection [2]. A straightforward
approach to reduce the beam measurement overhead is to utilize the UE location
information [5] to train an ML model for beam prediction. However,
transmission of UE location information, which may not necessarily be
available always to the BS, poses an additional feedback overhead. To avoid
this issue, the study in [6] fuses the concept of HBS with a supervised ML
model and exploits the spatial correlation among the parent and the child beam
qualities to predict the optimal child beam. A similar approach in [7]
utilizes the received signal vector of parent beams as an input to a
convolutional neural network (CNN). Another approach in [8] proposes to reduce
the beam measurement overhead by transmitting a subset of child beams and then
utilizes a CNN that predicts the optimal beam by learning the spatial
correlation among child beams.

Starting from 2022, the study of ML for the 5G-Advanced NR air interface is an
important project at 3GPP. Here, the focus is to explore the benefits of
augmenting the NR air interface with ML models for enhanced performance and/or
reduced overhead and complexity [9]. An important study item in this project
is the evaluation of ML for beam management, where spatial and temporal-domain
beam prediction are the sub use cases [10]. Following 3GPP guidelines,
companies report their proposed evaluation methodology and results on ML-based
beam prediction [11]. A recent proposal for spatial-domain beam prediction is
presented in [12], where based on the received power of a subset of the
transmit beams, a CNN is trained to predict the RSRPs of the non-transmitted
beams resulting in reduced overhead.

Though most of the discussed ML solutions reduce the beam measurement overhead
while achieving a performance closer to EBS, no significant attention has been
paid to the model computational complexity, model training time and its
generalization capabilities. To bridge this research gap, this letter presents
a low-complexity ML beam prediction approach that achieves the performance
closer to the optimal EBS but with lower computational complexity as compared
to other ML approaches, resulting in faster beam prediction. Additionally, to
investigate the generalization capabilities of our model, we evaluate its
performance over 3GPP specified scenarios.

##  II System Model

This section details channel and beam steering models, followed by an overview
of the beam management procedure.

###  II-A Channel Model

We consider a downlink mmWave multi-in multi-out (MIMO) communication system,
where the BS and the UE are equipped with $N_{\mathrm{T}}$ and
$N_{\mathrm{R}}$ antenna elements, respectively. Using the clustered channel
model, the channel is assumed to be the sum of the line-of-sight (LOS) path
and $C$ non-line-of-sight (NLOS) clusters with $L$ paths per cluster. The
channel matrix $\textbf{H}\in\mathbb{C}^{N_{\mathrm{R}}\times N_{\mathrm{T}}}$
can then be written as [13]

$\textbf{H}=\underbrace{\sqrt{\frac{K\Lambda}{K+1}}\alpha_{\text{LOS}}\textbf{a}_{\mathrm{R}}(\phi^{\mathrm{R}}_{\text{LOS}},\theta^{\mathrm{R}}_{\text{LOS}})\textbf{a}_{\mathrm{T}}^{H}(\phi^{\mathrm{T}}_{\text{LOS}},\theta^{\mathrm{T}}_{\text{LOS}})}_{\textbf{H}_{\text{LOS}}}\\\
+\underbrace{\sqrt{\frac{\Lambda}{L(K+1)}}\sum_{c=1}^{C}\sum_{l=1}^{L}\alpha_{c,l}\textbf{a}_{\mathrm{R}}(\phi^{\mathrm{R}}_{c,l},\theta^{\mathrm{R}}_{c,l})\textbf{a}_{\mathrm{T}}^{{H}}(\phi^{\mathrm{T}}_{c,l},\theta^{\mathrm{T}}_{c,l})}_{\textbf{H}_{\text{NLOS}}}.$
(1)

Here, the $l$-th path of the $c$-th cluster has azimuth (elevation) angle-of-
arrival (AoA) $\phi^{\mathrm{R}}_{c,l}(\theta^{\mathrm{R}}_{c,l})$ and azimuth
(elevation) angle-of-departure (AoD)
$\phi^{\mathrm{T}}_{c,l}(\theta^{\mathrm{T}}_{c,l})$, while $\alpha_{c,l}$ is
the complex path gain. The same variables are analogously defined for the LOS
path and are indicated by the LOS index. Furthermore,
$\textbf{a}_{\mathrm{R}}(\cdot)\in\mathbb{C}^{N_{\mathrm{R}}\times 1}$ and
$\textbf{a}_{\mathrm{T}}(\cdot)\in\mathbb{C}^{N_{\mathrm{T}}\times 1}$ denote
the UE and the BS array response, respectively, $(\cdot)^{{H}}$ denotes
conjugate transpose, $K$ is the Ricean factor, and $\Lambda$ indicates the
pathloss.

We assume a uniform planar array (UPA) in the $\mathrm{y}$-$\mathrm{z}$ plane
at the BS and the UE with $N_{\mathrm{y}}$ and $N_{\mathrm{z}}$ antenna
elements ($N_{\mathrm{y}}N_{\mathrm{z}}=N$) on $\mathrm{y}$ and $\mathrm{z}$
axis, respectively. Here, for ease of notation we drop the subscript for the
BS and UE. The array response vector for the UPA can then be written as

$\textbf{a}(\phi,\theta)=\frac{1}{\sqrt{N}}[1,\cdots,e^{j\frac{2\pi}{\lambda}d(y^{\prime}\sin(\phi)\sin(\theta)+z^{\prime}\cos(\theta))},\cdots,\\\
e^{j\frac{2\pi}{\lambda}d((N_{\mathrm{y}}-1)\sin(\phi)\sin(\theta)+(N_{\mathrm{z}}-1)\cos(\theta))}]^{{T}},$
(2)

where $y^{\prime}\in\\{0,1,\cdots,N_{\mathrm{y}}-1\\}$,
$z^{\prime}\in\\{0,1,\cdots,N_{\mathrm{z}}-1\\}$, while $\lambda$ and
$d=\frac{\lambda}{2}$ indicate the wavelength and antenna element spacing,
respectively.

###  II-B Beam Steering Model

We consider phase shifter based analog beamforming with one radio frequency
(RF) chain. At the BS the transmit signal is beamformed by a beamforming
vector
$\textbf{f}=[{f}_{1},{f}_{2},\cdots,{f}_{N_{\mathrm{T}}}]^{{T}}$$\in\mathbb{C}^{N_{\mathrm{T}}\times
1}$ and at the UE the received signals are combined with a receive combining
vector
$\textbf{w}=[{w}_{1},{w}_{2},\cdots,{w}_{N_{\mathrm{R}}}]^{{T}}$$\in\mathbb{C}^{N_{\mathrm{R}}\times
1}$. Here, $f_{i}$ and $w_{j}$ denote the complex weight on the $i$-th
transmit and $j$-th receive antenna element, respectively. The transmit and
receive beams are selected from the predefined codebooks $\mathcal{F}$ and
$\mathcal{W}$, consisting of $F$ and $W$ candidate beams, respectively. The
codebooks are designed on the following beam steering scheme.

$\displaystyle\textbf{f}\in\mathcal{F}$
$\displaystyle=\\{\textbf{a}_{\mathrm{T}}(\bar{\phi}^{\mathrm{T}}_{1},\bar{\theta}^{\mathrm{T}}_{1}),\textbf{a}_{\mathrm{T}}(\bar{\phi}^{\mathrm{T}}_{2},\bar{\theta}^{\mathrm{T}}_{2}),\cdots,\textbf{a}_{\mathrm{T}}(\bar{\phi}^{\mathrm{T}}_{F},\bar{\theta}^{\mathrm{T}}_{F})\\}$
(3) $\displaystyle\textbf{w}\in\mathcal{W}$
$\displaystyle=\\{\textbf{a}_{\mathrm{R}}(\bar{\phi}^{\mathrm{R}}_{1},\bar{\theta}^{\mathrm{R}}_{1}),\textbf{a}_{\mathrm{R}}(\bar{\phi}^{\mathrm{R}}_{2},\bar{\theta}^{\mathrm{R}}_{2}),\cdots,\textbf{a}_{\mathrm{R}}(\bar{\phi}^{\mathrm{R}}_{W},\bar{\theta}^{\mathrm{R}}_{W})\\}$
(4)

Here, $\bar{\phi}_{m}^{\mathrm{T}}(\bar{\theta}_{m}^{\mathrm{T}})$ for the
$m$-th transmitting beam $\textbf{f}_{m}$, $m\in\\{1,2,\cdots,F\\}$ and
$\bar{\phi}_{n}^{\mathrm{R}}(\bar{\theta}_{n}^{\mathrm{R}})$ for the $n$-the
receiving beam $\textbf{w}_{n}$, $n\in\\{1,2,\cdots,W\\}$ are the quantized
azimuth (elevation) AoD and AoA, respectively. Given the channel matrix H, the
transmit signal $x$, the $m$-th transmitting beam $\textbf{f}_{m}$ and the
$n$-th receiving beam $\textbf{w}_{n}$, the received signal $y_{m,n}$ is

$y_{m,n}=\sqrt{P}\textbf{w}_{n}^{H}\textbf{H}\textbf{f}_{m}x+\textbf{w}_{n}^{H}\boldsymbol{\eta},$
(5)

where ${P}$ is the transmit power and
$\boldsymbol{\eta}\in\mathbb{C}^{{N}_{\mathrm{R}}\times 1}$ is the additive
white Gaussian noise (AWGN).

###  II-C Beam Management in 5G NR

The 3GPP beam management procedure is based on the EBS and aims to find the
optimal beam pair $\\{\textbf{f}_{{m}^{*}},\textbf{w}_{{n}^{*}}\\}$ that
maximizes the RSRP given as: RSRP${}_{m,n}=|y_{m,n}|^{2}$. The optimization
problem can be formulated as

$\\{m^{*},n^{*}\\}=\underset{{\begin{subarray}{c}m\in\\{1,2,\cdots,F\\},\\\
n\in\\{1,2,\cdots,W\\}\end{subarray}}}{\mathrm{argmax}}\,\text{RSRP}_{m,n}.$
(6)

EBS solves this optimization problem by exhaustively searching over all
possible beamforming and combining vectors leading to an excessively huge beam
training overhead of $F\cdot W$ beam measurements.

To reduce this beam measurement overhead, HBS utilizes a multi-resolution
codebook and the problem of beam selection is divided into two levels. The
first-level search identifies the best parent beam by solving

$\\{m^{*}_{\mathrm{p}},n^{*}_{\mathrm{p}}\\}=\underset{{\begin{subarray}{c}m_{\mathrm{p}}\in\\{1,2,\cdots,F_{\mathrm{p}}\\},\\\
n_{\mathrm{p}}\in\\{1,2,\cdots,W_{\mathrm{p}}\\}\end{subarray}}}{\mathrm{argmax}}\,\text{RSRP}^{\mathrm{p}}_{m_{\mathrm{p}},n_{\mathrm{p}}}.$
(7)

Here, $F_{\mathrm{p}}=\frac{F}{s_{\mathrm{T}}}$ and
$W_{\mathrm{p}}=\frac{W}{s_{\mathrm{R}}}$ indicate the number of parent beams
at the BS and UE, respectively. Further, $s_{\text{T}}$ and $s_{\text{R}}$
defines the number of child beams within each parent beam at the BS and UE,
respectively. After identifying the best parent beam pair, the second-level
search confirms the optimal child beam pair within the range of the selected
parent beam pair (7), by

$\\{m^{*},n^{*}\\}=\underset{{\begin{subarray}{c}m\in\\{(m^{*}_{\mathrm{p}}-1)s_{\mathrm{T}}+1,\cdots,m^{*}_{\mathrm{p}}s_{\mathrm{T}}\\},\\\
n\in\\{(n^{*}_{\mathrm{p}}-1)s_{\mathrm{R}}+1,\cdots,n^{*}_{\mathrm{p}}s_{\mathrm{R}}\\}\end{subarray}}}{\mathrm{argmax}}\,\text{RSRP}^{\mathrm{c}}_{m,n}.$
(8)

Notably, the first and the second-level search requires $F_{\mathrm{p}}\cdot
W_{\mathrm{p}}$ and $s_{\mathrm{R}}\cdot s_{\mathrm{T}}$ beam measurements,
respectively, resulting in reduced beam measurement overhead. However, the
multi-level search incurs increased latency.

##  III Low-Complexity machine learning Design for mmWave Beam Prediction

In this section, we leverage the angular domain spatial correlation to propose
a low-complexity beam prediction model for fast beam training. Motivated by
the fact that very large antenna arrays can only be employed at the BS due to
size constraints, in the following sections, we limit our discussion to the
identification of the optimal transmit beam, i.e., the assumption of the
knowledge of the optimal receive beam.

###  III-A Algorithm Framework

Motivated by the two-level beam search, we propose to cover the whole angular
region with the first-level parent beams. By doing so, we observe that there
exists a strong angular spatial correlation among parent and child beams in a
certain environment. As an example, Fig. 1 shows the angular spatial
correlation between the RSRPs of the parent and the child beams, where each
parent beam contains four child beams. Here, it can be observed that the
parent beam has a stronger correlation with a limited number of child beams.
Consequently, we assume that the $\text{RSRP}^{\mathrm{c}}$ of the child beams
is a function $f_{1}(\cdot)$ of the parental RSRP values, i.e.,

$\text{RSRP}^{\mathrm{c}}=f_{1}(\text{RSRP}^{\mathrm{p}}).$ (9)

In particular, we aim on probing the parent beams and obtaining their
corresponding RSRPs from the received signal vector
$\textbf{y}^{\mathrm{p}}=[y_{1}^{\mathrm{p}},y_{2}^{\mathrm{p}},\cdots,y_{F_{\mathrm{p}}}^{\mathrm{p}}]^{{T}}$
and by intelligently merging these parent RSRPs with the strong correlation
among parent and child beams, we can predict the optimal child beam index
$m^{*}$. Due to the discrete number of candidate beams, the beam prediction
problem can be formulated as multiclass-classification problem and can be
written as

$m^{*}=f_{2}(\text{RSRP}^{\mathrm{p}}),\quad m^{*}\in\\{1,2,\cdots,F\\}$ (10)

where $f_{2}(\cdot)$ is the function that learns the correlation between
parent and child RSRPs for optimal beam index prediction. Further, due to the
highly non-linear relationship between RSRPs and channel directivity, the
prediction is difficult to be estimated by conventional signal processing
methods. With this background, we propose a low-complexity ML design for beam
prediction in the following section.

Figure 1: Spatial correlation among RSRPs of parent and child beams.

###  III-B Model Design

In this section, we introduce our ML model and its corresponding inputs and
outputs as shown in Fig. 2.

####  III-B1 Input Layer

Based on our previous discussions, the RSRPp of the parent beams obtained via
the first level of traditional HBS is provided as an input to the model. This
indicates that the input layer consists of $F_{\mathrm{p}}$ nodes. As an
example, considering $F=64$ beams and selecting $s_{\mathrm{T}}=4$ results in
$F_{\mathrm{p}}=16$ parent beams which means that a beam measurement overhead
reduction of $1-\frac{16}{64}=75\%$ is achieved as compared to the EBS.

####  III-B2 Output Layer

For the prediction of the optimal child beam from all the candidate child
beams, a fully-connected (FC) layer, consisting of $F$ nodes is introduced,
which learns the spatial correlation between RSRPp and RSRPc and transforms it
to the candidate child beams. Finally, a non-linear softmax activation layer
is introduced that returns the probabilities of all the child beams. The
output of the proposed low-complexity neural network (NN) can be written as

$\hat{\mathbfcal{P}}=\text{softmax}(\textbf{A}^{{T}}\text{\acs{RSRP}}^{\mathrm{p}}+\textbf{b}).$
(11)

Here, $\hat{\mathbfcal{P}}\in\mathbb{C}^{F\times 1}$ is the predicted output
probability vector of all the child beams, while
$\textbf{A}\in\mathbb{C}^{F_{\mathrm{p}}\times F}$ and
$\textbf{b}\in\mathbb{C}^{F\times 1}$ are the weights and the biases,
respectively. Finally, the child beam with maximum probability
$\hat{\mathcal{P}}_{m}$ is selected, i.e.,

$\hat{m}^{*}=\underset{{\begin{subarray}{c}m\in\\{1,2,\cdots,F\\}\end{subarray}}}{\mathrm{argmax}}\,\hat{\mathcal{P}}_{m}.$
(12)

pp^P∈$\Cstr[0]_{0}$$\Cstr[0]_{0}$$\vdots$Softmax Figure 2:  Proposed low-
complexity ML design for beam prediction.

##  IV Performance Evaluation

This section details detaset generation, model training, complexity analysis,
and performance evaluation over specified key performance indicators. For
reproducibility of results, our simulation dataset and source code is publicly
available [14].

###  IV-A Dataset Generation and Model Training

For dataset collection, we utilize the EBS approach in combination with HBS.
Our dataset consists of parent RSRP measurements, i.e., RSRPp obtained via the
traditional HBS and are provided as input features to the ML model. In
addition, the offline training labels, i.e., optimal beam indices are obtained
via the traditional EBS [15]. Table I lists default simulation parameters. The
location of the UE is drawn based on a uniform spatial distribution in the
cell coverage area. The noise power $\sigma^{2}$ is computed as
$(-174+10\text{log}_{10}B+N_{\text{F}})$ dBm and the path loss is given as
$(20\text{log}_{10}d+20\text{log}_{10}f_{\mathrm{c}}-147.56)$ dB, where $d$
indicates distance. Finally, the channel model is considered as a clustered
delay line (CDL) model [13]. Further, to investigate the generalization
capabilities of our ML model, we consider following scenarios with different
combinations of channel profiles [15].

  * •

Scenario ${1}$: The ML model is trained based on a training dataset
constructed by utilizing the CDL-D channel profile and performs inference on
the UE with same channel profile but with unknown location.

  * •

Scenario ${2}$: The ML model is trained based on a training dataset
constructed by utilizing the CDL-D channel profile and performs inference on a
UE with the CDL-E channel profile and with unknown location.

  * •

Scenario ${3}$: The ML model is trained on the mixed dataset from above
scenarios and performs inference on the UE of both channel profiles but with
unknown location.

Our dataset consists of $25,000$ samples, where the training, validation, and
testing data split is $70\%$, $10\%$, and $20\%$, respectively. Further, the
ML model is trained for $n_{\mathrm{e}}=100$ epochs, the model parameters are
optimized by the Adam optimizer [16] with the mean square error as loss
function.

TABLE I: List of simulation parameters. Parameters | Values
---|---
No. of BS antennas | $N_{\mathrm{T}}$ | $64$ |
UE antenna configuration | $N_{\mathrm{R}}$ | $8$ |
BS codebook size | $F$ | $64$ |
BS parent codebook size | $F_{\mathrm{p}}$ | $16$ |
Transmit power | $P$ | $30$ | dBm
BS antenna gain |  | $8$ | dBi [13]
UE codebook size | $W$ | $8$ |
UE noise figure | $N_{\mathrm{F}}$ | $6$ | dB
UE antenna gain |  | $5$ | dBi [17]
Center frequency | $f_{\mathrm{c}}$ | $28$ | GHz
Bandwidth | $B$ | $100$ | MHz
Sub-carrier spacing |  | $120$ | kHz
Cell radius |  | $200$ | m

###  IV-B Key Performance Indicators

For performance evaluation in terms of beam measurement overhead, the KPI is
selected as reference signalling overhead reduction (%) $1-\frac{N}{M}$, where
$N$ is the number of beams (SSBs) required as input by the ML model and $M$ is
the total number of beams to be predicted [15]. For beam prediction accuracy,
the KPI Top-$K$ (%) is defined as the percentage that the truly optimal genie-
aided transmit beam is among the $K$ best beams predicted by the ML model and
the beam prediction error (%) is calculated as $1-\text{beam prediction
accuracy}$. Here, the Top-$1$ genie-aided transmit beam is obtained via EBS
[15]. Further, the beam prediction accuracy is also evaluated in terms of
achieved average RSRP. Finally, for complexity analysis, we compare the model
complexity in terms of number of trainable parameters and number of floating-
point operations.

###  IV-C Complexity Analysis

An important measure of ML model complexity is the number of trainable
parameters ($n_{\mathrm{l}}$), which for an FC-NN layer with $n_{\mathrm{i}}$
inputs and $n_{\mathrm{o}}$ outputs can be computed as
$n_{\mathrm{l}}=(n_{\mathrm{i}}+1)n_{\mathrm{o}}$. Consequently, for proposed
model the number of parameters are $n_{\mathrm{l}}=(F_{p}+1)F=1088$. Further,
the number of trainable parameters for a convolutional layer can be obtained
as
$n_{\mathrm{l}}=n_{\mathrm{f}}(f_{\mathrm{h}}f_{\mathrm{w}}f_{\mathrm{d}}+1)$,
where $n_{\mathrm{f}},f_{\mathrm{h}},f_{\mathrm{w}},$ and $f_{\mathrm{d}}$
indicate the number of filters, filter height, width, and depth, respectively.
We evaluate the complexity in terms of model size with $32$-bit precision.
Table II indicates that due to a smaller number of trainable parameters the
proposed model has the smallest size as compared to other models.

TABLE II: Computational complexity comparison. |  No. of Trainable Parameters |  Model Size (Mbits) |  No. of FLOPs
---|---|---|---
FC-NN in [5] |  $17,728$ |  $0.5$ |  $17,728$
CNN in [6] |  $352,034$0 |  $11.2$0 |  $1,370,000$
CNN in [7] |  $67,008$ |  $2.1$ |  $332,000$
CNN in [12] |  $739,073$0 |  $23.6$0 |  $47,300,000$
Proposed model |  0$1,088$ |  0$0.04$ |  $1,088$

The time complexity of our proposed ML model is compared in terms of number of
required FLOPs using Big-$\mathcal{O}$ notation. During training, the ML model
performs forward and backward pass and it is useful to analyze the training
and inference time complexity. In both forward and backward pass, the
trainable parameters of a layer with $w$ nodes are updated by a matrix-vector
multiplication resulting in a time complexity of $\mathcal{O}(w^{2})$ FLOPs.
Furthermore, considering an NN with $l$ layers, $w$ nodes per layer, and
training the network with $n_{\mathrm{d}}$ data samples, and for
$n_{\mathrm{e}}$ epochs requires
$\mathcal{O}(n_{\mathrm{e}}n_{\mathrm{d}}lw^{2})$ FLOPs during training, while
the inference requires only $\mathcal{O}(n_{\mathrm{d}}lw^{2})$ FLOPs as only
forward pass is performed during inference. Similarly, the time complexity of
a CNN, with $l_{\mathrm{c}}$ convolutional and $l$ FC layers during training
is
$\mathcal{O}(n_{\mathrm{e}}n_{\mathrm{d}}(n_{\mathrm{f}}l_{\mathrm{c}}i_{\mathrm{h}}i_{\mathrm{w}}(f_{\mathrm{h}}f_{\mathrm{w}}f_{\mathrm{d}}))+lw^{2})$.
Here, in addition to the parameters defined above, $i_{\mathrm{h}}$ and
$i_{\mathrm{w}}$ indicate input height and width, respectively. Further, the
inference time complexity is then given as
$\mathcal{O}(n_{\mathrm{d}}(n_{\mathrm{f}}l_{\mathrm{c}}i_{\mathrm{h}}i_{\mathrm{w}}(f_{\mathrm{h}}f_{\mathrm{w}}f_{\mathrm{d}}))+lw^{2})$.

Table II summarizes the complexity comparison with the state of the art. For a
fair comparison, the number of estimated FLOPs are for one epoch and one data
sample, i.e., $n_{\mathrm{e}}=n_{\mathrm{d}}=1$. Here, it can be seen that the
proposed ML model achieves significantly lower computational complexity and
benefits from lower power consumption. Further, the execution of the proposed
ML model on an Intel i7-1185G7 processor indicates that the training time per
epoch and per data sample is $9\text{\,}\mathrm{\SIUnitSymbolMicro s}$, which
allows efficient and less time consuming model retraining. Besides, the
execution time for each prediction is around
$2\text{\,}\mathrm{\SIUnitSymbolMicro s}$ allowing faster beam prediction.

###  IV-D Simulation Results

For performance evaluation, in addition to the two-level HBS, CNN from [6],
and the FC-NN from [5], the EBS based beam selection is selected as a baseline
for comparison [15]. During inference the input to all ML models are the RSRPp
measurements of the parent beams and the outputs are the predicted
probabilities of each child beam being the best.

In terms of beam measurement overhead, the baseline EBS requires $64$ beam
measurements, resulting in $100\%$ beam measurement overhead. HBS requires
$16$ parent and $4$ child beam measurements, resulting in beam measurement
overhead of $32\%$. For the ML models, during inference, the measurement
overhead depends on the value of $K\in\\{4,2,1\\}$, reflecting the necessity
of probing the remaining $K$ beams for final selection, resulting in beam
measurement overhead of around $32\%$, $28\%$, and $25\%$, respectively, as
shown in Fig. 3. In terms of beam prediction error for $K=1$, our proposed
approach reduces the error by around $2$, $1.4$, and $1$ percentage points as
compared to HBS, [5], and [6], respectively. Similar observations can be made
from Fig. 4, where the performance is compared in terms of the average RSRP.
Here it can be noticed that the mean RSRP achieved by all ML approaches is
well within a $0.15$ dB margin of the genie-aided (EBS) transmit beam.
However, it is worth mentioning that the HBS achieves similar performance at
the cost of increased latency.

Figure 3: Beam prediction error and beam measurement overhead for scenario 1.
Figure 4: Comparison in terms of achieved average RSRP [dBm] for scenario 1.

Fig. 5 showcases the generalization capabilities of our proposed model over
three different scenarios as discussed in Section IV-A. We observe that for ML
Top-1 the prediction error of the model increases by around $5$ percentage
points for scenario $2$, due to different channel profiles used in training
and testing. Further, the error can be reduced when the model is trained on a
mixed data set from different channel profiles, i.e., scenario $3$. However,
the error in scenario $3$ is still around $1.2$ percentage points higher as
compared to scenario $1$. An important observation made here is that training
a model for a large number of scenarios results in reduced inference
performance for a specific scenario. Thus, there exists a trade-off between ML
model accuracy performance and its generalization capabilities.

Figure 5: Generalization capabilities of proposed beam prediction model over
3GPP specified scenarios discussed in Section IV-A.

##  V Conclusion

This letter proposes an ML-based beam prediction design that reduces the
reference signaling overhead and predicts the transmit beam with higher
accuracy and much lower computational complexity as compared to the state-of-
the-art. Specifically, we formulated the beam prediction problem as a
multiclass-classification task and proposed a low-complexity ML design to
learn the spatial angular correlation between parent and child beams to
predict the optimal beam. Due to lower computational complexity, the proposed
model reduces the power consumption at the UE and the beam prediction time
making it suitable for faster beam prediction. Further, through simulation
results, we showed that there exists a trade-off between ML model performance
and its generalization capabilities. These 3GPP compliant evaluation results
indicate the feasibility of ML-based mmWave beam prediction for 5G-Advanced NR
and beyond 5G communication networks.

## References

  * [1] T. S. Rappaport _et al._ , “Wideband millimeter-wave propagation measurements and channel models for future wireless communication system design,” _IEEE Trans. Commun._ , vol. 63, no. 9, pp. 3029–3056, 2015.
  * [2] M. Qurratulain Khan _et al._ , “Machine learning for millimeter wave and terahertz beam management: A survey and open challenges,” _IEEE Access_ , vol. 11, pp. 11 880–11 902, 2023.
  * [3] Q. Li _et al._ , “Machine learning based time domain millimeter-wave beam prediction for 5G-advanced and beyond: Design, analysis, and over-the-air experiments,” _IEEE J. Sel. Areas Commun._ , vol. 41, no. 6, pp. 1787–1809, 2023.
  * [4] Z. Xiao _et al._ , “Hierarchical codebook design for beamforming training in millimeter-wave communication,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 5, pp. 3380–3392, 2016.
  * [5] Y. Heng and J. G. Andrews, “Machine learning-assisted beam alignment for mmwave systems,” _IEEE Trans. on Cogn. Commun. Netw._ , vol. 7, no. 4, pp. 1142–1155, 2021.
  * [6] H. Echigo _et al._ , “A deep learning-based low overhead beam selection in mmWave communications,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 1, pp. 682–691, 2021.
  * [7] K. Ma _et al._ , “Deep learning assisted calibrated beam training for millimeter-wave communication systems,” _IEEE Trans. Commun._ , vol. 69, no. 10, pp. 6706–6721, 2021.
  * [8] Z. Wang _et al._ , “Intelligent beam training with deep convolutional neural network in mmwave communications,” in _Proc. IEEE Global Commun. Conf._ , 2022.
  * [9] 3GPP RP-213599, “New study item: Study on artificial intelligence (AI)/machine learning (ML) for NR air interface,” Qualcomm, 3GPP TSG RAN Meeting #94e, Dec. 2021.
  * [10] 3GPP, “RAN1 chair’s notes,” 3GPP TSG RAN WG1 #109-e, May. 2022.
  * [11] ——, “RAN1 chair’s notes,” 3GPP TSG RAN WG1 #114, Aug. 2023.
  * [12] ——, “Evaluation methodology and results on AI/ML for beam managment,” Keysight Technologies, TSG RAN WG1 Meeting, Aug. 2023.
  * [13] 3GPP, “Study on channel model for frequencies from 0.5 to 100 GHz,” TR 38.901 V17.0, 2022.
  * [14] A low-complexity machine learning design for mmWave beam prediction. [online]. Available: https://github.com/MuhammadKhan86/Low-Complexity-Beam-Prediction.
  * [15] 3GPP, “Study on artificial intelligence (AI)/machine learning (ML) for NR air interface,” TR 38.843 V1.0, 2023.
  * [16] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [17] 3GPP, “Study on NR coverage enhancements,” TR 38.830, V.17.0.0, 2020.

  *[3GPP]: 3rd Generation Partnership Project
  *[ML]: machine learning
  *[5G]: fifth generation
  *[NR]: New Radio
  *[mmWave]: millimeter-wave
  *[BS]: base station
  *[UE]: user equipment
  *[EBS]: exhaustive beam scan
  *[HBS]: hierarchical beam scan
  *[CNN]: convolutional neural network
  *[RSRPs]: reference signal received power
  *[MIMO]: multi-in multi-out
  *[LOS]: line-of-sight
  *[NLOS]: non-line-of-sight
  *[AoA]: angle-of-arrival
  *[AoD]: angle-of-departure
  *[UPA]: uniform planar array
  *[RF]: radio frequency
  *[AWGN]: additive white Gaussian noise
  *[RSRP]: reference signal received power
  *[FC]: fully-connected
  *[NN]: neural network
  *[CDL]: clustered delay line
  *[KPI]: key performance indicator
  *[SSBs]: synchronization sequence block
  *[FLOPs]: floating-point operation