# Bohemian Matrix Geometry

Robert M. Corless<EMAIL_ADDRESS>0000-0003-0515-1572 Cheriton School
of Computer Science
University of WaterlooCanada ,  George Labahn<EMAIL_ADDRESS>Cheriton
School of Computer Science
University of WaterlooCanada ,  Dan Piponi<EMAIL_ADDRESS>0000-0002-6599-2835 Epic GamesUSA and  Leili Rafiee Sevyeri
0000-0002-7819-9764<EMAIL_ADDRESS>Cheriton School of
Computer Science
University of WaterlooCanada

(2022)

###### Abstract.

A Bohemian matrix family is a set of matrices all of whose entries are drawn
from a fixed, usually discrete and hence bounded, subset of a field of
characteristic zero. Originally these were integers—hence the name, from the
acronym BOunded HEight Matrix of Integers (BOHEMI)—but other kinds of entries
are also interesting. Some kinds of questions about Bohemian matrices can be
answered by numerical computation, but sometimes exact computation is better.
In this paper we explore some Bohemian families (symmetric, upper Hessenberg,
or Toeplitz) computationally, and answer some open questions posed about the
distributions of eigenvalue densities.

Bohemian matrix, height, upper Hessenberg, Toeplitz, complex symmetric

††journalyear: 2022††copyright: acmlicensed††conference: Proceedings of the
2022 International Symposium on Symbolic and Algebraic Computation; July 4–7,
2022; Lille, France††booktitle: Proceedings of the 2022 International
Symposium on Symbolic and Algebraic Computation (ISSAC ’22), July 4–7, 2022,
Lille, France††price: 15.00††doi: 10.1145/xxxx.yyyy††isbn: not-an-isbn-
number††ccs: Computing methodologies Hybrid symbolic-numeric methods††ccs:
Computing methodologies Linear algebra algorithms††ccs: Computing
methodologies Symbolic calculus algorithms

##  1\. Introduction

Figure 1. Eigenvalues of $5\cdot 10^{6}$ dimension $m=8$ matrices with entries
chosen uniformly at random among the population $P=(-1,0,1)$. We see that the
uniform distribution as $m\to\infty$ result of (tao2006random, ) is evident
already. As $m$ increases, the “holes” close up, and the relative percentage
of the real eigenvalues (which extend past the $\sqrt{m}$ disk radius) becomes
negligible.

A fuzzy-edged light-grey disk, bisected with a dark line with dark bumps at
$0$ and $\pm 1$. Reminiscent of the moon occulted by a dark shadow. Smaller
holes are placed symmetrically on on near-circular arcs.

A Bohemian matrix family is a set of matrices all of whose entries are drawn
from a fixed finite population, usually integers, algebraic integers, or
Gaussian integers. The name “Bohemian” was invented in 2015 at the Fields
Institute Thematic Year in Symbolic Computation; the mnemonic is useful
because it highlights searching for commonality among features of matrices
with discrete populations. Our original interest was for software testing, and
as a testing ground for optimization over (in search of improved computational
bounds for certain quantities, such as the growth factor in Gaussian
elimination with complete pivoting, or the departure from normality). Bohemian
matrices are a specialization in the sense of Pólya, and have led now to
several workshops, at Manchester in 2018, at ICIAM in 2019, and at SIAM in
2021. There have been several publications since, including (Chan2020, ),
(fasi2020determinants, ), (DudaFonsecaXuYe+2021+505+514, ), (Sendra2022, ),
and the very interesting (bogoya2022upper, ) which explores a connection to
the asymptotic spectral theory of Toeplitz matrices (Schmidt1960, ), which is
very much alive today: see e.g. (bottcher2012introduction, ) and
(barrera2018eigenvalues, ).

The study of matrices with rational integer entries is very old, and the
literature too vast to survey coherently here. We instead point to the early
survey by Olga Taussky–Todd (Taussky1960, ) as an entry point. We are also
going to be working with Gaussian integer and algebraic integer entries; see
for instance (butson1962generalized, ) for important work on generalized
Hadamard matrices where the entries are roots of unity.

The study of random matrices where the entries are drawn from discrete
distributions is also very advanced; see (tao2006random, ; tao2017random, )
for instance. Those papers established that dense square matrices of dimension
$m$ whose entries are drawn from a discrete population, say $-1$, $0$, and
$1$, have eigenvalues that are asymptotically uniformly distributed on a disk
of radius $\sqrt{m}$. See Figure 1.

However, if the matrices are structured, other pictures arise, and little is
known about the asymptotic distributions of their eigenvalues. For instance,
in (Corless:2021:WhatCan, ) we find skew-symmetric tridiagonal matrices with
unexpectedly square distributions, or even diamond-shaped, as in Figure 2. We
will explain some of these shapes in this paper, and prove that they will not
fill out to disk shapes as $m\to\infty$. We will also explain some other
interesting features that arise in certain Bohemian families, including upper
Hessenberg Toeplitz Bohemians.

Figure 2. Eigenvalues of all dimension $m=31$ skew-symmetric tridiagonal
matrices with entries drawn from $P=(-1,i,1,-i)$, the fourth roots of unity.
Picture courtesy Aaron Asner. We explain the unexpected eigenvalue geometry in
this paper.

A pinkish fractal diamond shape with a hole in the center outlined in yellow
and red with a green and blue cross surrounding the hole. The fractal edges on
the outside are mostly green. The fractal edges are reminiscent of a saw
blade. The sharpest peaks are at the four main corners $(\pm\sqrt{2},0)$,
$(0,\pm\sqrt{2})$.

There is a significant connection to number theoretic works, as well. Kurt
Mahler (mahler1963two, ) was interested in the distribution of zeros of
polynomials with given length (the one-norm of the vector of coefficients) and
height (the infinity-norm of the vector of coefficients). This is connected
with the Littlewood conjecture for polynomials (littlewood1968, ) (How large
on the unit circle must a polynomial with $-1$,$1$ coefficients be?). The
numerical visualization of zeros of polynomials with coefficients $0$ or $1$
was apparently first done in (Odlyzko1992, ), who proved that the limiting set
was connected; later work explained the “holes" (borwein1997polynomials, ) and
visualizations by Peter Borwein and Loki Jörgenson made several other
questions clearer (borwein2001visible, ). See also the web pages of John
Carlos Baez and of Dan Christensen. Their article, The Beauty of Roots,
published on Baez’ website at the previous link, explains quite a few of the
visible structures. Then we will see a connection to Kate Stange’s work on
Schmidt Tessellations, https://math.katestange.net/illustration/schmidt-
arrangements/. See also (harriss2020algebraic, ) and (dorfsman2021searching, )
who connect Galois theory and visualization of roots of polynomials (and
therefore, although they do not point this out, of eigenvalues of Bohemian
matrices).

It is only a small jump from polynomials of bounded height to matrices of
bounded height; but the questions become (to our minds) even more interesting.

###  1.1. Organization of the paper

In section 2 we mention a few research questions about Bohemian matrices that
may be interesting to the computer algebra community; in 3 we discuss specific
matrix structures and give our main theorems, which describe and explain the
constraints that these matrix structures place on the spectra. In particular,
we give a new theorem about eigenvalues of upper Hessenberg Toeplitz matrices
which gives an essentially complete explanation of the “fractal” edges seen in
some of the figures. This is based on the well-known asymptotic spectral
theory of Toeplitz matrices, but extended to the case where we have an
uncountable number of such matrices in the limit as the dimension goes to
infinity. This new theorem extends a result of Schmidt and Spitzer, which is
concerned with Toeplitz matrices whose “symbol” is a Laurent polynomial and
with certain semi-algebraic curves that arise from that Laurent polynomial, to
matrices whose symbol is a Laurent series.

Together these theorems explain the appearance of some of these figures. We
also explain some of the “algebraic number starscape" appearance
(harriss2020algebraic, ) and connect to Schmidt tesselations (Stange2017, ) by
making an approximate computation of eigenvalues, and displaying the results
in Figure 5.

##  2\. Some questions of interest

Every polynomial written in the monomial basis can be embedded as a Frobenius
companion matrix into a matrix of the same height (the height of a matrix
$\mathbf{A}$, as opposed to a polynomial, is the infinity norm of the matrix
reshaped into a vector). Therefore every question about roots of polynomials
of bounded height translates directly into a question about eigenvalues of
Bohemian matrices. It will become clear as we go that this is one-way, that
is, there are questions of Bohemian matrices that do not translate into
questions about bounded height polynomials.

One question is “which matrices in the family have the largest characteristic
height?” The characteristic height is the height of the characteristic
polynomial; as previously noted, the characteristic height might be
exponentially larger than the matrix height. This is so for certain upper
Hessenberg Toeplitz matrices (Chan2020, ), where a lower bound containing a
Fibonacci number is given for the maximum characteristic height in the family
studied in that paper.

The characteristic height is connected to the numerical conditioning of the
characteristic polynomial; indeed one may take the Lebesgue constant for the
polynomial (CorlessFillion2013, , ch. 8) on the interval $-1\leq x\leq 1$ to
be the characteristic height.

Typically, one asks these questions in an asymptotic form; one wants answers
valid in the limit as the dimension $m$ goes to infinity. Examples of this
include work by Tao and Vu, who show that for general unstructured matrices
the distribution of eigenvalues divided by $\sqrt{m}$ is asymptotically
uniform on the unit disk (tao2006random, ; tao2017random, ). This is not true
for structured matrices; see e.g. (Corless:2021:WhatCan, ; Corless2021, )
where skew-symmetric tridiagonal matrices (independently of dimension, so no
scaling is needed) are seen to be confined to a diamond shape. We explain that
mystery in this paper, using a century-old theorem, which deserves to be
better-known.

##  3\. The Effect of Matrix Structure

The first Bohemian results were on real symmetric matrices or Hermitian
matrices (Tracy1993, ). We look at other structures here, in order to get
another view. We begin with complex symmetric matrices.

###  3.1. Complex Symmetric Matrices

A complex symmetric matrix $\bm{A}$ satisfies $\bm{A}^{T}=\bm{A}$, where T is
the real transpose operation. These occur, for instance in Bézout matrices for
polynomials with complex coefficients. Unlike Hermitian matrices, the
eigenvalues of complex symmetric matrices need not be real. Indeed, any matrix
may be brought by similarity transformation to a complex symmetric matrix
(horn2012matrix, , Thm 4.4.9). In many cases this can be done by unitary
similarity; see for instance, the characterizations of when this can be done,
in (Liu2013, ).

Here let us examine a specific complex symmetric family, with population
$-1\pm i$ where $i=(0,1)$ is the square root of $-1$. At dimension $m$, such a
matrix has $m(m+1)/2$ free entries, each of which can be either of $-1\pm i$.
This gives $2^{m(m+1)/2}$ such matrices; this growth is (much) faster than
exponential. Still, examining eigenvalues of small dimension examples can tell
us much. If we take $m=6$, then the number of such matrices is only $2^{21}$,
slightly more than $2$ million. The eigenvalues of all these matrices can be
computed in a reasonable time, and plotted. As depicted in Figure 3, they seem
confined to a strip in the left-half plane.

We now prove that this will always be true at any dimension.

###### Theorem 3.1.

If the symmetric matrix $\bm{A}$ has dimension $m$, entries drawn from $-1\pm
i$, and eigenvalue $\lambda$, then $-m\leq\Re(\lambda)\leq 0$ and
$-m\leq\Im(\lambda)\leq m$.

###### Proof.

Write $\bm{A}=-\bm{E}+i\bm{M}$ where $\bm{E}=\mathbf{e}\mathbf{e}^{T}$ is the
rank-one matrix that has all $1$s, and $\bm{M}$ is symmetric and has entries
only $\pm 1$. We use the Bendixon–Bromwich–Hirsch theorem (hogben2013handbook,
, Fact 5, p. 16-2) (original references (Bendixson1902, ; 1902, ;
Bromwich1906, )) as follows. Note that $(\bm{A}+\bm{A}^{*})/2=-\bm{E}$ is the
Hermitian part of $\bm{A}$, while $(\bm{A}-\bm{A}^{*})/(2i)=\bm{M}$ is the
skew-Hermitian part. The Bendixon–Bromwich–Hirsch theorem says that if the
eigenvalues of $-\bm{E}$ are written $\mu_{1}\geq\mu_{2}\geq\ldots\geq\mu_{m}$
(they are all real because $\bm{E}$ is Hermitian) and the eigenvalues of
$\bm{M}$ are written $\nu_{1}\geq\nu_{2}\geq\ldots\geq\nu_{m}$, then
$\mu_{m}\leq\Re(\lambda)\leq\mu_{1}$ and $\nu_{m}\leq\Im(\lambda)\leq\nu_{1}$.
But a short computation shows that the eigenvalues of $-\bm{E}$ are $-m$ and
$0$ with multiplicity $m-1$, and the Gerschgorin disk theorem shows that the
eigenvalues of $\bm{M}$ lie in the union of circles centred at $1$ of radius
$m-1$ and centred at $-1$ of radius $m-1$. This establishes the theorem. ∎

###### Remark 1.

We could have stated and proved that theorem with greater generality. That is,
the population could equally well have been $a\pm bi$ for real numbers $a$ and
$b$ and the conclusions would have been the same, apart from scaling. However,
we only need the idea, for the next theorem.

Figure 3. All eigenvalues with $\Im(\lambda)\geq 0$ of symmetric dimension
$m=6$ matrices with entries $-1\pm i$. The set is symmetric about the real
axis. We see the eigenvalues apparently confined to the strip
$-6\leq\Re(\lambda)\leq 0$, and bounded below and above by
$-6\leq\Im(\lambda)\leq 6$. There are several other unexplained features of
this eigenvalue distribution.

We now use this method to explain why skew-symmetric tridiagonal matrices with
population $1$ and $i$ (the population considered in (Corless:2021:WhatCan, )
and (Corless2021, )) are confined to a diamond shape.

####  3.1.1. Squares, Diamonds, and Lozenges

###### Theorem 3.2.

Let $\bm{A}$ be a square skew-symmetric matrix of dimension $m$ and population
$-1\pm i$. Then its eigenvalues $\lambda$ satisfy $-2\leq\Re(\lambda)\leq 2$
and $-2\leq\Im(\lambda)\leq 2$, and are thus confined to a square.

###### Proof.

A matrix $\bm{A}$ from this family can be written as $\bm{A}=\bm{S}+i\bm{T}$
where the superdiagonal of $\bm{S}$ is $-1$ and the superdiagonal of $\bm{T}$
is $\pm 1$. Both are skew-symmetric: $\bm{S}^{T}=-\bm{S}$ and
$\bm{T}^{T}=-\bm{T}$. The Hermitian part of $\bm{A}$ is
$(\bm{A}+\bm{A}^{H})/2=i\bm{T}$ and the skew-Hermitian part is
$(\bm{A}-\bm{A}^{H})/(2i)=-i\bm{S}$. [Both are Hermitian, and have real
eigenvalues.] Application of the Gerschgorin circle theorem to each of these
shows that the eigenvalues of either matrix are confined to the interval
$-2\leq\mu\leq 2$. By the theorem of Bendixon–Bromwich–Hirsch cited earlier,
the eigenvalues of $\bm{A}$ are confined to the square $-2\leq\Re(\lambda)\leq
2$, $-2\leq\Im(\lambda)\leq 2$. ∎

###### Remark 2.

Skew-symmetric tridiagonal matrices with population $1$ and $i$ are confined
to a diamond $|\Re(\lambda)|+|\Im(\lambda)|\leq\sqrt{2}$. To see this,
multiply the matrix by $-1+i$, which rotates its eigenvalues by $\pi/4$ and
stretches them by $\sqrt{2}$; but now the matrix population is $-1\pm i$ and
it’s still skew-symmetric, and hence confined to the square as described
above. Rotate the square back by $\pi/4$ and shrink by $\sqrt{2}$, and the
result follows.

We have therefore explained the non-round shape of the eigenvalue distribution
of these matrices.

###  3.2. Upper Hessenberg Matrices

An upper Hessenberg matrix is a matrix of the form

(3.1) $\begin{bmatrix}h_{1,1}&h_{1,2}&h_{1,3}&\ldots&h_{1,m}\\\
h_{2,1}&h_{2,2}&\ddots&&\vdots\\\ 0&h_{3,2}&h_{3,3}&\ddots&\vdots\\\
\vdots&&\ddots&\ddots&\vdots\\\ 0&\ldots&0&u_{m,m-1}&u_{m,m}\end{bmatrix}\>.$

That is, it is zero below the first subdiagonal. If any entry of the first
subdiagonal is zero, then the matrix is said to be reducible, because the
matrix then separates into blocks containing distinct eigenvalues. Here we
restrict our attention to irreducible matrices and indeed we specify that all
the subdiagonal entries $h_{j+1,j}=-1$. If the subdiagonal entries all have
$|h_{j+1,j}|=1$ we say that the matrix is unit upper Hessenberg. If all the
diagonal entries are zero, we say that it is Zero Diagonal.

###### Remark 3.

Suppose that all $h_{i,j}$ are roots of unity. Consider the $5$ by $5$ case,
for definiteness. Take four roots of unity, as yet unspecified: call them
$s_{2}$, $s_{3}$, $s_{4}$, $s_{5}$. Form the diagonal matrix
$\bm{D}=\mathrm{diag}{(}1,s_{2},s_{3},s_{4},s_{5})$ and perform the similarity
transform $\bm{D}\bm{H}\bm{D}^{-1}$. The result is

(3.2)
$\left[\begin{array}[]{ccccc}h_{1,1}&\frac{h_{1,2}}{s_{2}}&\frac{h_{1,3}}{s_{3}}&\frac{h_{1,4}}{s_{4}}&\frac{h_{1,5}}{s_{5}}\\\
s_{2}h_{2,1}&h_{2,2}&\frac{s_{2}h_{2,3}}{s_{3}}&\frac{s_{2}h_{2,4}}{s_{4}}&\frac{s_{2}h_{2,5}}{s_{5}}\\\
0&\frac{s_{3}h_{3,2}}{s_{2}}&h_{3,3}&\frac{s_{3}h_{3,4}}{s_{4}}&\frac{s_{3}h_{3,5}}{s_{5}}\\\
0&0&\frac{s_{4}h_{4,3}}{s_{3}}&h_{4,4}&\frac{s_{4}h_{4,5}}{s_{5}}\\\
0&0&0&\frac{s_{5}h_{5,4}}{s_{4}}&h_{5,5}\end{array}\right]\>.$

Now choose $s_{2}=\overline{h_{2,1}}$, $s_{3}=s_{2}\overline{h_{3,2}}$,
$s_{4}=s_{3}\overline{h_{4,3}}$, and $s_{5}=s_{4}\overline{h_{5,4}}$. This
forces the subdiagonal entries to be $1$, leaves the diagonal as it was
before, and shuffles the upper triangle to be possibly different roots of
unity to what they were before. Since the elements of the upper triangle were
independent, the resulting formulae will range over the entire set of
possibilities as we make the $h_{i,j}$ range over the entire set of
possibilities. It is for this specific population that we started studying
unit upper Hessenberg matrices.

Many populations have zero as their mean value (e.g. $\\{-1,1\\}$ or indeed
roots of unity; or $\\{-1,0,1\\}$). In that case, the eigenvalues are
typically symmetric about zero and a simplified picture is obtained simply by
setting all the diagonal entries to zero, in which case the Gerschgorin
circles are all centred at $0$. Sometimes we will have to transform back to
the nonzero diagonal case, but a surprising amount of information is retained
even with the simplification of insisting on a zero diagonal.

Figure 4. Eigenvalues of a sample of $5$ million upper Hessenberg matrices of
dimension $m=5$ matrices with population cube roots of unity. The image is
visually indistinguishable from the density plot of all $14,348,907$ unit
upper Hessenberg zero diagonal matrices of dimension $m=6$.

####  3.2.1. Rayleigh Quotients

Recall the Rayleigh Quotient:

(3.3) $r=\frac{\mathbf{y}^{T}\bm{A}\mathbf{x}}{\mathbf{y}^{T}\mathbf{x}}\>.$

If $\mathbf{x}$ is an approximate eigenvector, and $\mathbf{y}^{T}$ an
approximate left eigenvector, then this quotient is a least-squares
approximation to an eigenvalue of $\bm{A}$. If we replace $\bm{A}$ by
$\bm{A}^{-1}$ above, then this is an approximation to an eigenvalue of
$\bm{A}^{-1}$, and typically the largest one; this of course is the reciprocal
of the smallest eigenvalue of $\bm{A}$.

We will consider this not as an eigenvalue approximation, but as a process in
its own right, and plot the results of a single iteration of this on a
Bohemian family, with both $\mathbf{y}$ and $\mathbf{x}$ taken to be the first
elementary vectors. Thus the result is the top left corner of the inverse of
our Bohemian matrix. See Figure 5.

Figure 5. Density plot of the upper left corner element of $\bm{A}^{-1}$ where
$\bm{A}$ is sampled randomly from the $8$ by $8$ upper Hessenberg Bohemian
(not unit or zero diagonal) family with population cube roots of unity. There
are $3^{43}>3.2\cdot 10^{20}$ such matrices; we sampled only $2\cdot 10^{7}$
of these and gave a density plot on a $2048$ by $2048$ grid, enhanced by anti-
aliased point rendering.

This can also be computed by the recurrence relation (3.7), as the ratio of
two determinants, from Cramer’s rule:

(3.4)
$r=\frac{Q_{m}(0;t_{1},t_{2},\ldots,t_{m-1})}{Q_{m-1}(0;t_{1},t_{2},\ldots,t_{m-2})}$

Perhaps surprisingly to a numerical analyst, this offers an effective way to
perform this computation when the population consists of small Gaussian
integers, which can be represented as complex “flints” and are not subject to
rounding error when ring arithmetic is carried out in floats. In this case the
only division occurs at the end, and so rounding error is trivial. Even with
roots of unity, the rounding errors are generally not of serious consequence.

What we are computing here is representable in (if cube roots of unity are
used) $\mathbb{Q}(\sqrt{-3})$, and moreover the rational numbers involved will
not have overly large values. This suggests that the appearance in Figure 5
can be explained as Schmidt arrangements, as done by Katherine Stange. See her
blog at https://math.katestange.net/, but note in particular (Stange2017, )
and her papers previously cited.

###  3.3. Unit upper Hessenberg zero diagonal matrices

We will use the following Theorem repeatedly, to ensure that we see all the
eigenvalues of the family in the window of the plot.

###### Theorem 3.3.

Suppose that every entry of a unit upper Hessenberg zero diagonal matrix
$\mathbf{H}$ has magnitude at most $B$: that is, $|h_{i,j}|\leq B$. Then every
eigenvalue $\lambda$ of $\mathbf{H}$ is bounded by

(3.5) $|\lambda|\leq 1+2\sqrt{B}\>.$

This bound is independent of the dimension. The proof uses an idea already
present in (Schmidt1960, ), and doubtless in other places. Namely one chooses
a diagonal matrix $\bm{D}=\mathrm{diag}{(}1,r,r^{2},\ldots,r^{m-1})$ with a
free parameter $r>1$ and considers the similar matrix
$\bm{D}\bm{H}\bm{D}^{-1}$ which has the same eigenvalues as $\bm{H}$. Suppose,
without loss of generality, the subdiagonal entries of the zero diagonal unit
upper Hesssenberg $\bm{H}$ are all $-1$. Then

(3.6)
$\bm{D}\bm{H}\bm{D}^{-1}=\begin{bmatrix}0&h_{1,2}/r&h_{1,3}/r^{2}&\ldots&h_{1,m}/r^{m-1}\\\
-r&0&h_{2,3}/r&&h_{2,m}/r^{m-2}\\\ 0&-r&0&h_{3,4}/r&\vdots\\\
\vdots&&\ddots&\ddots&\vdots\\\ 0&\ldots&0&-r&0\end{bmatrix}\>.$

The Gerschgorin disk for the first row has radius at most $B/(r-1)$ by
comparison with a geometric series; the second and all subsequent rows have
the bound $r+B/(r-1)$ for the radius which is larger because $r>1$. To
minimize this bound, we write it as $1+r-1+B/(r-1)$ and use the AGM inequality
to say that this is minimized when $r-1=B/(r-1)$ or $r=1+\sqrt{B}$; this gives
the value of the Gershgorin radius as $1+2\sqrt{B}$, as desired.

###### Remark 4.

This is the only Gerschgorin-like theorem that we are aware of that gives a
bound for eigenvalues which is independent of the dimension and depends on the
square root of the bound for the entries in the matrix instead of the more
usual linear power of the bound. Of course, if one multiplies a matrix
$\bm{A}$ by a constant factor, then the eigenvalues must also be multiplied by
that factor; but we cannot perform such a multiplication here and remain in
the class of unit upper Hessenberg matrices.

###  3.4. Unit upper Hessenberg zero diagonal Toeplitz Matrices

A Toeplitz matrix $\bm{T}$ has constant elements on every diagonal, that is,
$t_{i,j}=t_{0,j-i}$. The authors of (Chan2020, ) found that the Toeplitz
subset of unit upper Hessenberg matrices maximized the characteristic height,
and so decided to study that subset directly; it contains only exponentially
many elements ($\\#P^{m-1}$ entries, rather than $\\#P^{O(m^{2})}$) and has
several other interesting features. For large dimension, the spectral theory
connects to the well-known asymptotic spectral theory for Toeplitz matrices.
Although the population does not decay, we can use the method of 3.3 to
transform it to one for which we may compute the so-called “symbol".

Figure 6. Density plot of eigenvalues of all $1,048,576$ upper Hessenberg
Toeplitz matrices of dimension $11$ with zero diagonal, $-1$ subdiagonal, and
population $\pm 1$, $\pm i$ (fourth roots of unity) otherwise. Brighter
colours correspond to higher density. Figure 7. Density plot of eigenvalues of
all $262,144$ upper Hessenberg Toeplitz matrices of dimension $10$ with zero
diagonal, $-1$ subdiagonal, and population $\pm 1\pm i$ (four corners of a
square) otherwise. Brighter colours correspond to higher density.

####  3.4.1. The Fractal Edges

At the edges of Figures 6 and 7 we see clear indication of fractal gasket-like
structures. When the population is third roots of unity, we see Sierpinski
gaskets; when the population only has two elements, we see pairs, and pairs of
pairs, recursively; with five-element populations we see recursive pentagonal
structures. The phenomenon seems universal for upper Hessenberg Toeplitz
matrices. We now give an explanation for this behaviour in this section.

The keys to understanding this are the spectral theory of Toeplitz matrices,
the improved Gerschgorin bound of Theorem 3.3, and to consider the
characteristic polynomials. For this structure, there is a recurrence relation
for the characteristic polynomial, say $Q_{n}(z)$, of a unit upper Hessenberg
zero diagonal matrix, namely

$\displaystyle Q_{n+1}(z;t_{1},t_{2},\ldots,t_{n})=$ $\displaystyle
zQ_{n}(z;t_{1},t_{2},\ldots,t_{n-1})$ (3.7)
$\displaystyle{}-\sum_{k=1}^{n}(-1)^{k}t_{k}Q_{n-k}(z;t_{1},t_{2},\ldots,t_{n-k-1})\>.$

There is a somewhat more complicated recurrence relation for a general upper
Hessenberg matrix; see (Chan2020, ). Proofs can be found in many places, for
instance (Cahill2002, ). The final term is $\pm t_{n}Q_{0}(z)$ and
$Q_{0}(z)=1$.

Using Theorem 3.3 with $r=2$ (this is also used in (bogoya2022upper, )) we
find the Toeplitz symbol to be one of the (uncountably many)

(3.8) $-\frac{2}{e^{i\theta}}+\sum_{k\geq
1}\frac{t_{k}}{2^{k}}e^{ik\theta}\>.$

All of these share the same bound, however; and the eigenvalues will lie in
the union of the convex hulls of the zeros of these symbols.

More than that, though, the eigenvalues of any Toeplitz matrix with a finite
symbol are known to converge, as the matrix dimension goes to infinity, to the
union of several algebraic curves described by the Schmidt–Spitzer theorem
(hogben2013handbook, , Ch. 22). An improved algorithm for computing these
curves is given in (Bottcher2021, ), which we have implemented. For
convenience, we state the Schmidt–Spitzer theorem here.

###### Theorem 3.4.

Given the Laurent polynomial

(3.9) $f(z)=\sum_{\nu=-q}^{h}a_{\nu}z^{\nu}$

where $q>0$, $h>0$ and $a_{-q}a_{h}\neq 0$, and given $\lambda\in\mathbb{C}$,
define the polynomial

(3.10) $Q(\lambda;z)=z^{q}(f(z)-\lambda)\>.$

Denote the moduli of its $q+h$ zeros counted according to multiplicity by
$\alpha_{j}$ for $1\leq j\leq q+h$. Assume that they have been ordered so that
$\alpha_{1}\leq\alpha_{2}\leq\cdots\leq\alpha_{q+h}$. Then
$\lambda\in\Lambda$, that is the set of Schmidt–Spitzer points, if and only if
the $q$th largest and $(q+1)$st largest roots have the same magnitude:
$\alpha_{q}=\alpha_{q+1}$. Recall that $q$ is fixed by the Laurent polynomial.

For us, $q=1$; we will comment further on this below.

What we observed experimentally, as demonstrated in Figure 8, is that the
Schmidt–Spitzer curves themselves seem to converge rapidly—for upper
Hessenberg Toeplitz matrices, which have very special Toeplitz symbols—as we
add new sequence elements $t_{k}$.

This convergence is seen to be general by the following theorem.

###### Theorem 3.5.

Suppose a fixed infinite sequence $t_{k}\in P$is given, where each
$|t_{k}|\leq B$. The Schmidt–Spitzer curves for the truncated Toeplitz symbols
$a_{m,r}(z)=-r/z+\sum_{k=1}^{m-1}t_{k}z^{k}/r^{k}$ converge in the Haußdorff
metric to the set $\Lambda$ of piecewise analytic arcs inside
$|z|<1+2\sqrt{B}$ in the limit as $m\to\infty$.

###### Proof.

The curves are all contained in the disk $|z|\leq 1+2\sqrt{B}$ by Theorem 3.3,
and in the union over $r>1$ of the images $a_{r}(e^{i\theta})$ of the unit
circle. The Laurent polynomial symbols $a_{m,r}(z)$ obtained by truncation to
degree $m-1$ converge uniformly to an analytic function $a_{r}(z)$ as
$m\to\infty$ in $|z|<1+2\sqrt{B}$, by Theorem 2.7a in (Henrici, , Vol I). By
Hurwitz’ theorem (titchmarsh1939theory, ) the zeros of the sequence
$a_{m,r}(z)-\lambda$ also converge to a (possibly multiple) zero of
$a_{r}(z)-\lambda$, for any fixed $\lambda$, and every zero of
$a_{r}(z)-\lambda$ has the appropriate number of zeros of the sequence
converge to it. If $\xi$ is a zero of $a_{r}(z)-\lambda$, then $\xi/r$ is a
zero of $a_{r}(z)-\lambda$ with $r=1$; call this $a(z)-\lambda$. If the two
smallest roots ($q=1$) of $a(z)-\lambda$ have equal magnitude, then
$\lambda\in\Lambda$, the limiting set of eigenvalues of the Toeplitz matrices;
moreover, the set $\Lambda$ is not empty because the case all $t_{j}=0$ for
$j>h$ is included in the convergence above, and one may perturb from that
using analyticity. Hurwitz’ theorem then gives the desired result, when
applied to any $\lambda\in\Lambda$. ∎

###### Remark 5.

The Schmidt–Spitzer theorem entails that, for a Laurent polynomial with a pole
of order $q$, it is the zeros of $q$th smallest magnitude that matter. For our
application, $q=1$ because the matrix is not only Toeplitz but also unit upper
Hessenberg with zero diagonal. So it is the smallest magnitude zeros that
matter, and these are the first ones to converge. Indeed it can happen that
the largest roots “converge" to a complicated object at the boundary, which is
the image of the circle of radius $r$. Consider the example where all
$t_{k}=1$. The symbol is then, when $r=1$,

(3.11) $a(z)=-\frac{1}{z}+\sum_{k\geq 1}z^{k}$

or $-1/z+z/(1-z)$. Only the smallest roots of $a(z)-\lambda$ converge. The
compactness, and the scaling by $r>1$, really matter.

Applying the improved algorithm of (Bottcher2021, ), we solve
$a(z)-a(e^{i\theta}z)=0$ and find

(3.12)
$z=\frac{{\mathrm{e}}^{\mathrm{i}\theta}+1\pm\sqrt{\left({\mathrm{e}}^{\mathrm{i}\theta}\right)^{2}-6\,{\mathrm{e}}^{\mathrm{i}\theta}+1}}{4\,{\mathrm{e}}^{\mathrm{i}\theta}}$

Drawing $\lambda=a(z)$ in the complex plane for $-\pi\leq\theta\leq\pi$ gives
the correct half-circle determining the curve $\Lambda$ describing the
asymptotic location of the eigenvalues of the Toeplitz matrices with this
symbol.

###### Remark 6.

We have not investigated the case $q>1$, that is, Toeplitz matrices with more
subdiagonals. We suspect that the above theorem will go through, mutatis
mutandis, perhaps with slower convergence the more nonzero subdiagonals the
matrices have.

(a) $t_{1}=1$

(b) $t_{1}=1$, $t_{2}=-1$

(c) $t_{3}=1$

(d) $t_{4}=-1$

(e) $t_{5}=-1$

(f) $t_{6}=-1$

(g) $t_{7}=-1$

(h) $t_{8}=1$

(i) $t_{9}=1$

Figure 8. The Schmidt–Spitzer curves for Toeplitz matrices with symbol
$a(z)=-1/z+\sum_{k=1}^{m-1}t_{k}z^{k}$ for a fixed sequence of $t_{k}$ and
various $m$. As we increase $m$, and thus see a new member of the sequence,
the curves are seen to converge quite rapidly. After $m=10$ they are visually
indistinguishable from the final case shown here. The red dots are the
eigenvalues of the dimension $m$ Toeplitz matrix with those entries. In the
limit as $m\to\infty$ these are guaranteed to converge to the Schmidt–Spitzer
curves but we see that they are essentially there from the beginning.

We have

(3.13) $Q_{n+1}(z)=F_{n}(z)+t_{n}$

where $F_{n}(z)=zQ_{n}(z)-\sum_{k=1}^{n-1}(-1)^{k}t_{k}Q_{n-k}$ is a fixed
polynomial depending only on previous $t_{k}$. This final term $t_{n}$
perturbs that fixed polynomial in one of (in this case) three ways. Now use a
homotopy argument, replacing $t_{n}$ by $st_{n}$ where $0\leq s\leq 1$, we see
the roots of $Q_{n+1}$ arising by paths emanating in three directions from the
roots of $F(z)$. That is, for each root of $F(z)$, three nearby roots of
different $Q_{n+1}(z)$ arise. The roots of $F(z)$ are also near the zeros of
the symbol. This explains the fractal structure.

This recursive construction is not a linear one: the structures resembling
Sierpinski gaskets seen in the close-up in Figure 9 are clearly not rigid
triangles, but rather have been distorted into curved shapes. Nonetheless we
believe the above explanation is one way of understanding why this structure
arises.

Figure 9. A close-up (window $-2.5\leq\Re(\lambda)\leq-1.5$,
$0\leq\Im(\lambda)\leq 1$) of an $800$ by $800$ density plot of the
eigenvalues of all $531,441$ upper Hessenberg zero diagonal matrices with
population cube roots of unity. The resemblance to a Sierpinski gasket is
striking.

##  4\. On Visualization

Most of the figures in this present paper use only the simplest techniques of
visualization: Figures 1 and 2 show colourized density plots of eigenvalues in
the complex plane, about which more in a moment. Figure 3 is a simple plot.
Figure 9 is a greyscale density plot on an $800$ by $800$ grid. Figures 6 and
7 are colourized density plots, where the colouring scheme was chosen by using
a cumulative frequency count in order to attempt to equalize the apparent
density of eigenvalues using colour; this verges on true computer imaging
techniques but is actually very crude. The technique has some value because it
is relatively faithful to the underlying mathematics: brighter colours
correspond to higher eigenvalue density, and when the cividis or viridis
colour palette is used, the colours are relatively even perceptually. The
copper palette of Figure 4 retains the correlation of brightness to density,
but has a smaller colour range.

The appearance of Figure 5, however, depends on some rather more professional
techniques, as described in (Piponi2009, ). The basic idea is to estimate
spatial derivative information (using TensorFlow Gradient) and use that to
enhance the figure, making the density visible even with relatively sparse
data (for this figure, only $500,000$ matrices were used, and the
computational cost was substantially lower than for the other figures). The
technique is called “anti-aliased point rendering.” The picture remains
faithful to the underlying mathematics, however.

Figure 10 also uses this enhancement, this time because without it some
features of the eigenvalues (relative increase in density near the edge
$\Re(\lambda)=0$, for instance) are not so easy to see.

##  5\. Concluding remarks

The notion of a Bohemian matrix seems to be a remarkably productive one, with
substantial connections to very active areas of research, including
visualizations in number theory, combinatorial design, random matrices in
physics, numerical analysis, and computer algebra.

Many combinatorial questions about Bohemians, such as “how many different
characteristic polynomials are there” for a given Bohemian family (say unit
upper Hessenberg with population $(-1,0,1)$ for concreteness), can be
addressed computationally, but floating-point error is an issue. One is forced
to look at polynomials (and thus computer algebra, however implemented)
because of the multiple-eigenvalue problem: in those circumstances, numerical
computation of eigenvalues is ill-conditioned and one cannot really count
things by “clustering” nearby eigenvalues. One is often tempted, when thinking
of random matrices, to say that “multiple eigenvalues never happen” but of
course this is not true. Supplying constraints (either on the population or
the matrix structure) significantly enhances the probability that multiplicity
will be encountered.

One topic some of us have looked at briefly is that of stable matrices. Which
Bohemian matrices have all their eigenvalues strictly in the left half-plane?
For those matrices $\bm{A}$, and those matrices only, the solutions to the
linear differential equation $\dot{y}=\bm{A}y$ will ultimately decay to zero.
If the dimensions are large, then one may have to consider pseudospectra (and
thus matrix non-normality) as well.

It can be unsatisfactory to compute the eigenvalues of a matrix $\bm{A}$ and
check to see if they are all in the left half plane; rounding errors may drift
some of them into the right half plane. Computation of the characteristic
polynomial, and subsequent use of the Routh–Hurwitz criterion, seems in order.
One would like to take advantage of the compression seen for several families:
rather than computing eigenvalues of several million matrices, instead compute
the roots of the (equivalent) several thousand characteristic polynomials.
Better yet, apply the Routh–Hurwitz criterion, which is a rational criterion,
to make the decision in an arena uncontaminated by rounding errors.

As an example, consider the symmetric matrices with population $-1\pm i$ of
dimension $m=6$. We already know from Theorem 3.1 that all eigenvalues lie in
$\Re(\lambda)\leq 0$. But are there any of the $4,970$ characteristic
polynomials of these $2^{21}$ matrices which have all of their roots strictly
in the left half plane? Yes. By applying Maple’s Hurwitz tool to the
characteristic polynomials (which were actually computed using Python and
exported in a JSON container to Maple) we identified 1328 of these
polynomials, all of whose roots were strictly in $\Re(\lambda)<0$. Indeed, the
maximum real part was approximately $-1.03\cdot 10^{-5}$. Corresponding to
these $1328$ polynomials were $966,240$ matrices, or about $46\%$ of the
total.

For upper Hessenberg matrices with population $(-1,-1\pm i)$ and dimension
$m=4$, out of $1,594,323$ matrices we find $365,307$ distinct characteristic
polynomials. Of these, only $14,604$ (associated with $66,782$ matrices, about
$4.2\%$ of the total) have all their roots strictly in the left half plane.
The maximum real part is about $-7.1\cdot 10^{-5}$. See Figure 10. We enhanced
the figure to show more clearly the increase in density near $\Re(\lambda)=0$.
For comparison, we plot the eigenvalues of the whole family in Figure 11. One
sees immediately that the majority of matrices in this collection have
eigenvalues in the right half-plane.

Figure 10. A $1024$ by $1024$ grid density plot of the roots of all $14,604$
stable characteristic polynomials, enhanced by anti-aliased point rendering.
The Bohemian family is upper Hessenberg, population $-1-i$, $-1$, $-1+i$, and
dimension $m=4$. The plot is on $-L-1\leq\Re(\lambda)<L-1$,
$-L\leq\Im(\lambda)\leq L$, where $L=1+2\cdot 2^{1/4}$.  Figure 11. A $1200$
by $1200$ grid unenhanced density plot of the eigenvalues of the whole family,
to compare with the plot of Figure 10. The Bohemian family is upper
Hessenberg, population $-1-i$, $-1$, $-1+i$, and dimension $m=4$. The plot is
on $-L-1\leq\Re(\lambda)<L-1$, $-L\leq\Im(\lambda)\leq L$, where $L=1+2\cdot
2^{1/4}$. The spirals are completely unexplained.

###### Acknowledgements.

The Python code we used for many of the experiments was written largely by
Eunice Y. S. Chan (with some help by Rob Corless) for another project; we
thank her for her work on that code. We also thank her for many discussions on
Bohemian matrices and colouring algorithms. We thank Mark Giesbrecht for hours
of discussion on this topic. We thank Neil J. Calkin for valuable discussions
about finding good questions to direct research. We thank Ilias Kotsireas for
pointing out several references relevant to the Hadamard matrix literature,
and for helpful discussions. We thank Nick Higham for an independent proof of
Theorem 3.1, and for searching out the original references to Hirsch,
Bromwich, and Bendixon. We also thank Owen Maresh (@graveolens) for pointing
out the connection to Kate Stange’s work. Partially supported by NSERC grant
RGPIN-2020-06438, and partially supported by the grant PID2020-113192GB-I00
(Mathematical Visualization: Foundations, Algorithms and Applications) from
the Spanish MICINN.

## References

  * [1] Mauricio Barrera, Albrecht Böttcher, Sergei M Grudsky, and Egor A Maximenko.  Eigenvalues of even very nice Toeplitz matrices can be unexpectedly erratic.  In The diversity and beauty of applied operator theory, pages 51–77. Springer, 2018.
  * [2] Ivar Bendixson.  Sur les racines d'une équation fondamentale.  Acta Mathematica, 25(0):359–365, 1902.
  * [3] Manuel Bogoya, Stefano Serra-Capizzano, and Ken Trotti.  Upper Hessenberg and Toeplitz Bohemian matrix sequences: a note on their asymptotical eigenvalues and singular values.  Elec. Trans. Numer. Anal., 55:76–91, 2022.
  * [4] Peter Borwein and Loki Jörgenson.  Visible structures in number theory.  The American Mathematical Monthly, 108(10):897–910, 2001.
  * [5] Peter Borwein and Christopher Pinner.  Polynomials with $\\{$0,+ 1,-1$\\}$ coefficients and a root close to a given point.  Canad. J. Math., 49(5):887–915, 1997.
  * [6] Albrecht Böttcher, Juanita Gasca, Sergei M. Grudsky, and Anatoli V. Kozak.  Eigenvalue clusters of large tetradiagonal Toeplitz matrices.  Integral Equations and Operator Theory, 93(1), February 2021.
  * [7] Albrecht Böttcher and Bernd Silbermann.  Introduction to large truncated Toeplitz matrices.  Springer Science & Business Media, 2012.
  * [8] T. J. I’a Bromwich.  On the roots of the characteristic equation of a linear substitution.  Acta Mathematica, 30(0):297–304, 1906.
  * [9] A. T Butson.  Generalized Hadamard matrices.  Proceedings of the American Mathematical Society, 13(6):894–898, 1962.
  * [10] Nathan D. Cahill, John R. D'Errico, Darren A. Narayan, and Jack Y. Narayan.  Fibonacci determinants.  College Math. Jour., 33(3):221–225, May 2002.
  * [11] Eunice Y. S. Chan, Robert M. Corless, Laureano González-Vega, J. Rafael Sendra, and Juana Sendra.  Inner Bohemian inverses.  Applied Mathematics and Computation, accepted January 11, 2022, 2022\.
  * [12] Eunice Y. S. Chan, Robert M. Corless, Laureano Gonzalez-Vega, J. Rafael Sendra, Juana Sendra, and Steven E. Thornton.  Upper Hessenberg and Toeplitz Bohemians.  Linear Algebra and its Applications, 601:72–100, September 2020\.
  * [13] Robert M. Corless.  Skew-symmetric tridiagonal Bohemian matrices.  Maple Transactions, 1(2), October 2021.
  * [14] Robert M. Corless.  What can we learn from Bohemian matrices?  Maple Transactions, 1:31, 7 2021.
  * [15] Robert M. Corless and Nicolas Fillion.  A Graduate Introduction to Numerical Methods.  Springer, 2013.
  * [16] Gabriel Dorfsman-Hopkins and Candy Xu.  Searching for rigidity in algebraic starscapes.  arXiv preprint arXiv:2107.06328, 2021.
  * [17] Zhibin Du, Carlos M. da Fonseca, Yingqiu Xu, and Jiahao Ye.  Disproving a conjecture of Thornton on Bohemian matrices.  Open Math., 19(1):505–514, 2021.
  * [18] Massimiliano Fasi and Gian Maria Negri Porzio.  Determinants of normalized Bohemian upper Hessenberg matrices.  The Electronic Journal of Linear Algebra, 36:352–366, 2020.
  * [19] Edmund Harriss, Katherine E. Stange, and Steve Trettel.  Algebraic number starscapes.  arXiv preprint arXiv:2008.07655, 2020.
  * [20] Peter Henrici.  Applied and computational complex analysis, volume I.  Wiley, 1974.
  * [21] M. A. Hirsch.  Sur les racines d'une équation fondamentale.  Acta Mathematica, 25(0):367–370, 1902.
  * [22] Leslie Hogben.  Handbook of linear algebra.  CRC Press, 2nd edition, 2013.
  * [23] Roger A. Horn and Charles R. Johnson.  Matrix analysis.  Cambridge university press, 2012.
  * [24] J. E. Littlewood.  Some problems in real and complex analysis.  1968\.
  * [25] Xuhua Liu, Brice M. Nguelifack, and Tin-Yau Tam.  Unitary similarity to a complex symmetric matrix and its extension to orthogonal symmetric Lie algebras.  Linear Algebra and its Applications, 438(10):3789–3796, May 2013\.
  * [26] Kurt Mahler.  On two extremum properties of polynomials.  Illinois Journal of Mathematics, 7(4):681–701, 1963.
  * [27] Andrew Odlyzko.  Zeros of polynomials with $0$–$1$ coefficients.  In Bruno Salvy, editor, Algorithms Seminar 1992–1993. INRIA, 1992\.  Summary by Xavier Gourdon.
  * [28] Dan Piponi.  Two tricks for the price of one: Linear filters and their transposes.  Journal of Graphics, GPU, and Game Tools, 14(1):63–72, January 2009.
  * [29] Palle Schmidt and Frank Spitzer.  The Toeplitz matrices of an arbitrary Laurent polynomial.  Mathematica Scandinavica, 8(1):15–38, 1960.
  * [30] Katherine E. Stange.  Visualizing the arithmetic of imaginary quadratic fields.  International Mathematics Research Notices, 2018(12):3908–3938, February 2017.
  * [31] Terence Tao and Van Vu.  On random $\pm$1 matrices: singularity and determinant.  Random Structures & Algorithms, 28(1):1–23, 2006.
  * [32] Terence Tao and Van Vu.  Random matrices have simple spectrum.  Combinatorica, 37(3):539–553, 2017.
  * [33] Olga Taussky.  Matrices of rational integers.  Bulletin of the American Mathematical Society, 66(5):327–345, 1960\.
  * [34] Edward Charles Titchmarsh.  The theory of functions.  Oxford university press, 1939.
  * [35] Craig A. Tracy and Harold Widom.  Level-spacing distributions and the Airy kernel.  Physics Letters B, 305(1-2):115–118, May 1993.

##  Appendix A Maple code for the Schmidt–Spitzer theorem

The Maple code to implement our version of the algorithm of [6] can be found
in the Maple workbook at https://github.com/rcorless/Bohemian-Matrix-Geometry.
The algorithm is laid out in Algorithm 1.

Figure 12. An example Schmidt–Spitzer curve, together with the eigenvalues of
the smallest matrix in the Toeplitz family which the Laurent polynomial symbol
(A.14) applies. The closed curve is the graph of $a(e^{i\theta}/\rho)$.
Algorithm 1 Schmidt–Spitzer, specialized to unit upper Hessenberg zero-
diagonal case

0: Scale factor $\rho>1$. Vector $\mathbf{t}$ of Toeplitz matrix entries, of
length $m-1$. $t_{k}\in\mathbb{C}$.

1: Construct the scaled Laurent polynomial
$a(z)=-\rho/z+\sum_{k=1}^{m-1}t_{k}(z/\rho)^{k}$

2: Choose a vector of $\phi$ values in $-\pi\leq\phi_{j}\leq\pi$. More vector
entries mean a finer resolution of the Schmidt–Spitzer curves. We ignore the
case $\phi=0$ which requires special handling but only adds isolated points.

3: for $\phi_{\ell}\in\phi$ do

4: Solve the polynomial equation $za(z)-za(e^{i\phi_{\ell}}z)=0$. There are
$m$ roots $u_{j}$.

5: for $j$ to $m$ do

6: Compute $\lambda=a(u_{j})$.

7: solve $za(z)-z\lambda$. There are $m$ roots $v_{k}$ again, two of which are
$u_{j}$ and $e^{i\phi_{\ell}}u_{j}$.

8: If $|u_{j}|=|e^{i\phi_{\ell}}u_{j}|$ are the smallest roots in magnitude,
then $\lambda$ is on a Schmidt–Spitzer curve. Record it and continue.

9: end for

10: end for

As an example, we chose the vector $\mathbf{t}=[1,-1,1,0,1]$, so $m=6$. This
means that the symbol is the Laurent polynomial

(A.14) $a(z)=-\frac{1}{z}+\sum_{k=1}^{m-1}t_{k}z^{k}\>.$

For Figure 12 we chose $101$ values of $\phi$ equally-spaced on $[-\pi,\pi]$.
We chose the scale factor $\rho=1.75$ because it gives a reasonably tight
bound on the Schmidt–Spitzer curves when we draw the image $a(e^{i\psi}/\rho)$
of the unit circle under the map defined by the symbol. We also plotted the
eigenvalues of the dimension $m$ matrix with that population—this is the
smallest matrix in the family that this symbol and Schmidt–Spitzer curve apply
to. The visible agreement of eigenvalues and Schmidt–Spitzer curve is
satisfactory.


# On the Evaluation of Generative Models in
Distributed Learning Tasks

Zixiao Wang∗,1, Farzan Farnia${}^{~{},1}$, Zhenghao Lin2, Yunheng Shen3, Bei
Yu1
1 Chinese University of Hong Kong
2 South China Agricultural University 3 Tsinghua University
{zxwang22, farnia<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>Equal Contribution.

###### Abstract

The evaluation of deep generative models including generative adversarial
networks (GANs) and diffusion models has been extensively studied in the
literature. While the existing evaluation methods mainly target a centralized
learning problem with training data stored by a single client, many
applications of generative models concern distributed learning settings, e.g.
the federated learning scenario, where training data are collected by and
distributed among several clients. In this paper, we study the evaluation of
generative models in distributed learning tasks with heterogeneous data
distributions. First, we focus on the Fréchet inception distance (FID) and
consider the following FID-based aggregate scores over the clients: 1) FID-avg
as the mean of clients’ individual FID scores, 2) FID-all as the FID distance
of the trained model to the collective dataset containing all clients’ data.
We prove that the model rankings according to the FID-all and FID-avg scores
could be inconsistent, which can lead to different optimal generative models
according to the two aggregate scores. Next, we consider the kernel inception
distance (KID) and similarly define the KID-avg and KID-all aggregations.
Unlike the FID case, we prove that KID-all and KID-avg result in the same
rankings of generative models. We perform several numerical experiments on
standard image datasets and training schemes to support our theoretical
findings on the evaluation of generative models in distributed learning
problems.

##  1 Introduction

Deep generative models including diffusion models (Sohl-Dickstein et al.,
2015) and generative adversarial networks (GANs) (Goodfellow et al., 2014)
have attained impressive results over a wide array of machine learning tasks
(Karras et al., 2019, Ho et al., 2020, Ramesh et al., 2022). This success can
be attributed to the enormous capacity of multi-layer neural networks in
modeling complex distributions of image and text data as well as the intricate
design of the training mechanisms in GANs and diffusion models. The promising
results of these frameworks have inspired the development of several
methodologies for the training and evaluation of generative models in the
literature.

While the existing literature on deep generative models has mostly focused on
centralized settings with training data stored by a single learner, many
modern applications of deep learning algorithms are aimed at _distributed
scenarios_ where training data are collected by multiple agents in a network.
A well-known instance of such a distributed setting is the _federated learning
task_ (McMahan et al., 2017), where several clients are connected to a server
and aim to train a decentralized model through their communications with the
server node while preserving the privacy of their collected data. A
significant challenge in such distributed learning settings is the
heterogeneous data distributions across clients, since the background features
of every client could lead to a different data distribution. Especially, in
training a deep generative model over a distributed network, the heterogeneity
of the clients’ distributions could highly impact the performance and
evaluation of the trained model.

In this work, we focus on the evaluation of deep generative models in
heterogeneous distributed learning settings. Our primary goal is to highlight
the challenges of extending standard evaluation metrics for generative models
from the centralized setting to the heterogeneous distributed case. To do
this, we consider and analyze the following two sensible extensions of an
evaluation score: 1) the average score over clients (score-avg), i.e. the mean
of the evaluation scores for individual clients, 2) the score with respect to
the aggregate distribution of all clients (score-all) considering all the
samples in the network. We consider standard assessment scores and analyze the
consistency of model rankings suggested by the described extensions of the
score to the distributed learning problem.

In our theoretical analysis, we specifically focus on the Fréchet inception
distance (FID) (Heusel et al., 2017) and kernel inception distance (KID)
(Bińkowski et al., 2018) scores, which have been widely used in the evaluation
of generative models. For the FID score, not only do we prove that the FID-all
and FID-avg take different values, but further we show that the model rankings
according to the two scores can be different. Indeed, we prove that the
generative models with the optimal FID-avg and FID-all scores are different
under heterogeneous client distributions, revealing the discrepancy between
the two aggregate assessments.

On the other hand, in the case of KID score we prove that the aggregate scores
KID-all and KID-avg result in the same rankings of generative models although
they could take different values. This result shows the consistency of the KID
evaluation between the mean aggregation in KID-avg and the collective data-
based aggregation in KID-all. Therefore, our theoretical results suggest the
different consistency behaviors of FID and KID scores when aggregated over a
heterogeneous distributed learning setting: while averaging individual KID
scores results in the same ranking as the KID score with respect to the
distribution of collective client’s data, the same conclusion does not hold
for FID and more generally Wasserstein distance-based evaluation metrics.

Finally, we discuss several numerical results on standard image datasets and
generative model architectures which support our theoretical comparison of the
aggregate evaluation scores in federated learning problems. Our empirical
results suggest that the FID-all and FID-avg aggregations could lead to
inconsistent rankings of the trained generative models in standard GANs and
diffusion models. In our experiments, we observed that FID-avg often takes
higher values when generated samples look sharper. On the other hand, the FID-
all score seemed to be higher scores for models with higher diversity in their
generated samples. In contrast, KID-all and KID-avg resulted in consistent
rankings of the trained models where the evaluated score seems to aggregate
both quality and diversity performance. In our numerical experiments, we also
evaluated the precision, recall (Sajjadi et al., 2018, Kynkäänniemi et al.,
2019), density, and coverage (Naeem et al., 2020) scores under the two
introduced aggregations which, except in the case of recall score, could lead
to different rankings of the trained model under the discussed score
aggregations. In the following, we summarize the main contributions of our
study:

  * •

Highlighting the challenges of evaluating generative models in heterogeneous
distributed learning settings;

  * •

Analyzing two types of aggregate FID scores in distributed learning scenarios
and proving the inconsistencies between the optimal models under the two
scores;

  * •

Demonstrating the consistent rankings suggested by the client-wise averaged
KID and the aggregate-data-based KID evaluations;

  * •

Presenting numerical results on the evaluation of generative models in
distributed settings and empirically supporting the theoretical claims.

##  2 Related Work

A large body of related works (Theis et al., 2016) has focused on the
evaluation of generative models in standard centralized learning settings. The
existing evaluation scores can be categorized into two general groups: 1)
Distance-based metrics defining a distance between the distribution of
training data and learnt generative model. The distance between the real and
fake distributions is usually computed after passing the samples through a
pre-trained neural net offering a proper embedding of image data. The well-
known evaluation scores in this category are the FID (Heusel et al., 2017) and
KID (Bińkowski et al., 2018) scores. 2) Quality and diversity-based scores
which output a score based on the sharpness and variety of the generated
samples. The widely-used evaluation metrics in this category are the Inception
score (Salimans et al., 2016), precision and recall metrics (Sajjadi et al.,
2018, Kynkäänniemi et al., 2019), and density and coverage scores (Naeem et
al., 2020). We note that in this work our goal is not to introduce a new
evaluation metric, and our aim is to analyze the extensions of these scores to
distributed and federated learning problems under heterogeneous data
distributions across clients.

In another set of related works, extensions of generative model training
methods including GANs and diffusion models to distributed federated learning
have been studied. Rasouli et al. (2020) propose Fed-GAN to train GANs in a
federated learning setting. In (Hardy et al., 2019), the gradient from sample
generation for the generator is exchanged on the server, while each client
possesses a personalized discriminator. According to (Yonetani et al., 2019),
different weights are assigned to local discriminators in the non-i.i.d.
setting. Conversely, Wu et al. (2022) ensure client privacy by sharing the
discriminator across clients, while keeping the generator private.
Additionally, Su et al. (2023) explore a dual diffusion paradigm to extend
diffusion-based models into the federated learning setting, addressing
concerns related to data leakage.

##  3 Preliminaries

###  3.1 Deep Generative Models

In a deep generative model framework, a neural network generator $G$ is used
to map a hidden random vector $\mathbf{Z}\in\mathbb{R}^{d}$ drawn according to
a fixed distribution, e.g. isotropic multivariate Gaussian
$\mathcal{N}(\mathbf{0},\sigma^{2}I)$ to a real-like sample $G(\mathbf{Z})$.
Several deep learning approaches have been proposed to train such a generator
network, including maximum-likelihood-based methods such as variational
autoencoder (Kingma & Welling, 2013) and flow-based models (Dinh et al.,
2016), generative adversarial networks (GANs) (Goodfellow et al., 2014), and
denoising diffusion models (Ho et al., 2020). In our numerical evaluation, we
mainly concentrated on the latter two methods, GANs and diffusion models, due
to their state-of-the-art performance in computer vision applications.

In GANs, the training of generative models is framed as a min-max game between
a generator network $G$ mapping latent vector $\mathbf{Z}$ to a real-like
output and a discriminator network $D$ attempting to differentiate $G$’s
generated samples from real training data. The GAN game is typically
formulated as the following min-max optimization problem where
$\bm{\theta},\,\bm{\omega}$ represent the parameters of generator and
discriminator neural nets, and
$f\bigl{(}G_{\bm{\theta}},D_{\bm{\omega}}\bigr{)}$ is the min-max objective
representing $D$’s dissimilarity score for the generated and real samples:

$\min_{\bm{\theta}}\;\max_{\bm{\omega}}\;f\bigl{(}G_{\bm{\theta}},D_{\bm{\omega}}\bigr{)}.$
(1)

The training of GANs in a distributed learning problem aims at solving the
above problem via a distributed optimization method. For example, in a
federated learning setting where the local clients are connected to a single
server node, the training of GAN players can be achieved by a federated min-
max optimization algorithm as discussed in the related work section.

In the case of diffusion models, the generative model performs by multi-step
denoising of a Gaussian input. The training of this approach is typically done
by reversing the denoising process where the training data are turned to a
Gaussian input via an iterative addition of independent Gaussian noise
vectors. To extend diffusion models in federated learning, we follow the
simple FedAvg (McMahan et al., 2017) method and average the locally updated
diffusion networks at the server followed by synchronizing the clients with
the averaged model.

###  3.2 Distance-based Evaluation of Generative Models

In order to assess the performance of a generative model, a standard approach
is to measure the distance between the distribution of real and generated
data. Due to the high-dimensionality of standard image data, the evaluation of
image-based generative models is typically performed after passing the data
point through a pre-trained Inception model on the ImageNet dataset.

Specifically, a standard distance-based metric is the Fréchet inception
distance (FID) defined as the $2$-Wasserstein distance between two Gaussian
distributions with the mean and covariance parameters of the data distribution
$P_{X}$, denoted by $\bm{\mu}_{X},C_{X}$, and with the mean and covariance of
the generative model $P_{G}$, denoted by $\bm{\mu}_{G},C_{G}$:

$\mathrm{FID}(P_{X},P_{G}):=\bigl{\|}\bm{\mu}_{X}-\bm{\mu}_{G}\bigr{\|}^{2}_{2}+\mathrm{Tr}\Bigl{(}C_{X}+C_{G}-2\bigl{(}C_{X}C_{G}\bigr{)}^{1/2}\Bigr{)}.$

Note that we can interpret the FID score as an approximation of the
2-Wasserstein distance given the first and second-order moments of the
distributions.

Another widely-used distance-based score for the evaluation of generative
models is the kernel inception distance (KID), which measures the maximum mean
discrepancy (MMD) between the two distributions which is calculated using a
kernel similarity function
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. Here, the
definition of the MMD distance between $P_{X}$ and $P_{G}$ based on kernel $k$
follows from

$\mathrm{MMD}\bigl{(}P_{X},P_{G}\bigr{)}:=\mathbb{E}_{X,X^{\prime}\sim
P_{X}}\bigl{[}k(X,X^{\prime})\bigr{]}+\mathbb{E}_{Y,Y^{\prime}\sim
P_{G}}\bigl{[}k(Y,Y^{\prime})\bigr{]}-2\,\mathbb{E}_{X\sim P_{X},\,Y\sim
P_{G}}\bigl{[}k(X,Y)\bigr{]},$

where we suppose samples $X,X^{\prime}\sim P_{X}$ and $Y,Y^{\prime}\sim P_{G}$
are independently drawn.

##  4 Evaluation of Generative Models in Distributed Learning Settings

In this section, we discuss two extensions of distance-based evaluation scores
from a centralized case to heterogeneous distributed learning settings. In our
analysis, we use $\mathcal{D}(P_{X},P_{G})$ to denote a general distance
between data distribution $P_{X}$ and the generative model $P_{G}$. For
example, $\mathcal{D}$ can be chosen to be the FID score or KID score, which
we will analyze later in the section.

In a standard centralized setting, we have only a single distribution $P_{X}$
for real data. However, the main characteristic of a heterogeneous distributed
learning problem is the multiplicity of the involved clients’ distribution.
Here, we suppose a distributed setting $k$ clients and use
$P_{X_{1}},\ldots,P_{X_{k}}$ to denote their underlying distributions, i.e.
$P_{X_{i}}$ stands for the data distribution at client $i$. In addition, we
assume that every client $i$ contributes a fraction $0<\lambda_{i}<1$ of the
data in the network, that is $\lambda_{i}=\frac{n_{i}}{n}$ with
$n=\sum_{j=1}^{k}n_{j}$ is the total number of samples in the network and
$n_{i}$ is the number of samples at client $i$.

As a result of multiple input distributions, we need to define an aggregate
evaluation score that is based on distance measure $\mathcal{D}$. The
aggregate distance is supposed to summarize the performance of the generative
model $P_{G}$ in only one score. To do this, we consider and analyze two
reasonable ways of defining the aggregate score:

  1. 1.

Average Score $\mathcal{D}_{\mathrm{avg}}$: The score-avg is the mean of the
client’s individual distance measures, i.e.

$\mathcal{D}_{\mathrm{avg}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\,:=\,\sum_{i=1}^{k}\lambda_{i}\mathcal{D}\bigl{(}P_{X_{i}},P_{G}\bigr{)}.$
(2)

  2. 2.

Collective-data-based Score $\mathcal{D}_{\mathrm{all}}$: The score-all with
respect to the collective data of the clients is the distance between $P_{G}$
to the averaged distribution
$\widehat{P}_{X}:=\sum_{i=1}^{k}\lambda_{i}P_{X_{i}}$:

$\mathcal{D}_{\mathrm{all}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\,:=\,\mathcal{D}\Bigl{(}\widehat{P}_{X}\,,\,P_{G}\Bigr{)}.$
(3)

In the above, note that $\sum_{i=1}^{k}\lambda_{i}P_{X_{i}}$ is indeed a
mixture distribution with $k$ components $P_{X_{1}},\ldots,P_{X_{k}}$ with
frequency weights $\lambda_{1},\ldots,\lambda_{k}$. To relate the above
aggregate scores, we first observe that when $\mathcal{D}$ is a convex
function of the input distributions, which applies to both FID and KID scores,
the score-avg $\mathcal{D}_{\mathrm{avg}}$ will upper-bound the score-all
$\mathcal{D}_{\mathrm{all}}$:

######  Observation 1.

If $\mathcal{D}(P_{X},P_{G})$ is a convex function of $P_{X}$, then

$\mathcal{D}_{\mathrm{all}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\leq\mathcal{D}_{\mathrm{avg}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}$

######  Remark 1.

The convexity assumption on distance $D$ in the above observation applies to
standard divergence scores, including Wasserstein distances, $f$-divergence
measures, total variation distance, and the maximum mean discrepancy.
Consequently, the result applies to the FID and KID scores.

While the mentioned observation shows how the two aggregate scores are
compared with one another, it does not imply a monotonic relationship between
$\mathcal{D}_{\mathrm{all}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}$
and
$\mathcal{D}_{\mathrm{avg}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}$.
Therefore, this observation does not provide a comparison of the ranking of
generative models according to the two aggregate scores. In the following
subsections, we study this question for the particular FID and KID scores.

###  4.1 Aggregate FID Scores in Distributed Learning

In the case of FID score, we utilize the formulation of FID as the
2-Wassesrtein-distance between the Gaussian-fitted model that leads to a
Riemannian geometry. This observation results in the following theorem on FID-
all and FID-avg aggregations. We defer the proofs to the Appendix.

######  Theorem 1.

Suppose that $P_{X_{1}},\ldots,P_{X_{k}}$ are the clients’ distrbutions with
the mean parameters $\bm{\mu}_{1},\ldots,\bm{\mu}_{k}$ and covariance matrices
$C_{1},\ldots,C_{k}$, respectively, in the semantic space of the Inception
net. Then, the followings hold for a generative model $P_{G}$ with mean
$\bm{\mu}_{G}$ and covariance $C_{G}$.

  1. 1.

For FID-all, if we define a random variable $\widehat{X}$ with the average
mean $\widehat{\bm{\mu}}=\sum_{i=1}^{k}\lambda_{i}\bm{\mu}_{i}$ and covariance
matrix
$\widehat{C}=\sum_{i=1}^{k}\lambda_{i}\bigl{(}C_{i}+\bm{\mu}_{i}\bm{\mu}_{i}^{\top}-\widehat{\bm{\mu}}\widehat{\bm{\mu}}^{\top}\bigr{)}$
in the Inception-based semantic space, we have

$\mathrm{FID}_{\mathrm{all}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\,=\,\mathrm{FID}\bigl{(}P_{\widehat{X}},P_{G}\bigr{)}.$

  2. 2.

For FID-avg, if we define a random vector $\widetilde{X}$ with the average
mean $\widehat{\bm{\mu}}=\sum_{i=1}^{k}\lambda_{i}\bm{\mu}_{i}$ and covariance
matrix $\widetilde{C}$ as the unique solution to
$\widetilde{C}=\sum_{i=1}^{k}\lambda_{i}\bigl{(}\widetilde{C}^{1/2}C_{i}\widetilde{C}^{1/2}\bigr{)}^{1/2}$
in the Inception-based semantic space we will have

$\mathrm{FID}_{\mathrm{avg}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\,=\,\mathrm{FID}\bigl{(}P_{\widetilde{X}},P_{G}\bigr{)}+\sum_{i=1}^{k}\lambda_{i}\mathrm{FID}\bigl{(}P_{\widetilde{X}},P_{X_{i}}\bigr{)}.$

As a result, the FID-all score as a function of $G$ is changing monotonically
with $\mathrm{FID}\bigl{(}P_{\widetilde{X}},P_{G}\bigr{)}$ for defined
variable $\widetilde{X}$.

######  Remark 2.

In Theorem 1, $\widetilde{C}$ follows from the Wasserstein barycenter of the
approximate Gaussian distributions and under the assumption that
$C_{1},\ldots,C_{k}$ commute, i.e. $C_{i}C_{j}=C_{j}C_{i}$ for all $i,j$’s,
will simplify to

$\widetilde{C}=\Bigl{(}\sum_{i=1}^{k}\lambda_{i}C_{i}^{{1}/{2}}\Bigr{)}^{2}.$

######  Remark 3.

In Theorem 1, the optimal mean vectors for FID-all and FID-avg aggregations
are the same. In contrast, the optimal covariance matrix of FID-avg denoted by
$\widetilde{C}$ has no dependence on the choice of $\bm{\mu}_{i}$’s, while the
optimal covariance matrix of FID-all $\widetilde{C}$ will be affected by the
difference between $\bm{\mu}_{i}$’s due to the term
$\sum_{i=1}^{k}\lambda_{i}(\bm{\mu}_{i}\bm{\mu}_{i}^{\top}-\widehat{\bm{\mu}}\widehat{\bm{\mu}}^{\top})$.

As explained in the above remarks, Theorem 1 shows that the optimal covariance
matrices under FID-all and FID-avg could be significantly different in
heterogeneous settings with different $\bm{\mu}_{i}$’s. Therefore, since the
FID-all and FID-avg scores can be interpreted as the distance to covariance
matrices $\widehat{C}$ and $\widetilde{C}$, respectively, the rankings
suggested by the aggregations will be different if two generators’ covariances
$C_{G}$ and $C_{G^{\prime}}$ have different ordering of distances to
$\widehat{C}$ and $\widetilde{C}$.

###  4.2 KID-based Evaluation in Distributed Learning

After showing the possibility of inconsistent rankings by FID-all and FID-avg,
we consider the KID score and analyze the consistency of KID-all and KID-avg
aggregations. The following theorem proves that unlike the FID-case, the KID-
all and KID-avg will result in a consistent ordering of the models and there
is a monotonic relationship between the two aggregate scores.

######  Theorem 2.

Consider a kernel function
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ and the resulting
KID score. Then for the clients’ distributions $P_{X_{1}},\ldots,P_{X_{k}}$
with frequency parameters $\lambda_{1},\ldots,\lambda_{k}$, we will have the
following for the average distribution
$\widehat{P}_{X}=\sum_{j=1}^{k}\lambda_{j}P_{X_{j}}$:

$\mathrm{KID}_{\mathrm{avg}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}\,=\,\mathrm{KID}_{\mathrm{all}}\Bigl{(}P_{X_{1}},\ldots,P_{X_{k}}\,;\,P_{G}\Bigr{)}+\sum_{i=1}^{k}\lambda_{i}\,\mathrm{KID}\bigl{(}\widehat{P}_{X},P_{X_{i}}\bigr{)},$

which implies a monotonic relationship between KID-all and KID-avg as a
function of $P_{G}$.

Therefore, the above theorem shows the consistent rankings implied by the
aggregate KID-avg and KID-all scores, as the difference between the two scores
remains constant while changing generative model $P_{G}$.

##  5 Numerical Results

###  5.1 Evaluation on Synthetic Gaussian Mixture Datasets

As discussed in Remark 3, the optimal selection of the covariance matrix
differs for the FID-all and FID-avg aggregate scores. To illustrate this
distinction, we performed a toy experiment, revealing that FID-avg attains its
minimum value when the generator’s variance closely approximates that of an
individual client, whereas FID-all reaches its minimum value when the
generator’s variance aligns with that of the aggregated distribution.

Setup. Our experimental setup involves two clients, denoted by $C_{1}$ and
$C_{2}$. $C_{1}$ possesses a dataset consisting of 50,000 samples drawn from
the Gaussian distribution $\mathcal{N}([1,0]^{\top},\bm{\Sigma})$, while
$C_{2}$ holds a dataset with 50,000 samples drawn from
$\mathcal{N}([-1,0]^{\top},\bm{\Sigma})$, where
$\bm{\Sigma}=\text{diag}([1,1]^{T})$. We introduce a generator, denoted as
$G_{\mathrm{var}_{x}}$, which is parameterized by $\mathrm{var}_{x}$. The
value of $\mathrm{var}_{x}$ regulates the variance of the generator along the
X-axis. Specifically, $G_{\mathrm{var}_{x}}$ generates 50,000 data points
following a Gaussian distribution $\mathcal{N}([0,0]^{T},\bm{\Sigma_{G}})$,
where $\bm{\Sigma_{G}}=\text{diag}([\mathrm{var}_{x},1]^{T})$. The
relationship between the two clients and the generator is visually depicted in
Figure 1. Additionally, we introduce an "ideal estimator" denoted as
$\hat{E}=C_{1}\cup C_{2}$. This ideal estimator possesses the unique ability
to replicate the distribution of the training dataset perfectly. We employ the
ideal estimator as a reference for our analysis.

Evaluation Metrics. We measure the similarity between samples generated by
clients and generators using the Fréchet distance (FD), which follows from the
Wasserstein-based definition of FID-all and FID-avg without the application of
the pre-trained Inception network. We consider the aggregate scores FD-avg and
FD-all as defined in Equation 2 and Equation 3. Note that the FD-all for the
ideal estimator is zero and we use $\text{\rm FD-
ref}=\frac{1}{2}\sum_{i=1}^{2}\text{\rm FD}(\hat{E},C_{i})$ as a reference for
FD-avg. We also measure the Kernel distance (KD), which follows the definition
of KID-all and KID-avg without Inception network.

Results. By increasing $\mathrm{var}_{x}$ from 0 to 4, we get a sequence of
FD-avg / FD-all pairs and we plot them with the $\mathrm{var}_{x}$ in Figure
1. Our experimental results highlight the following conclusions. First, we
observed that the minimum of FD-all occurs at $\mathrm{var}_{x}=2$, while that
of FD-avg occurs at $\mathrm{var}_{x}=1$, which indicates that the optimal
solutions of $\mathrm{var}_{x}$ to minimize FD-all and FD-avg are
inconsistent. In this case, FD-all and FD-avg lead to different rankings of
the models with $\mathrm{var}_{x}=1$ and $\mathrm{var}_{x}=2$. Additionally,
we observed that, counterintuitively, the ‘ideal estimator’ did not reach the
minimum average of the Fréchet distances. The distance between KD-avg and KD-
all remains the same with the change of $\mathrm{var}_{x}$ and both of which
reach minimum at $\mathrm{var}_{x}=2$.

(a)  (b)  (c)

Figure 1: Experimental results of Gaussian mixture dataset. (a): The optimal
$\mathrm{var}_{x}$ parameters are different under FD-avg and FD-all
evaluations. (b): Distance between KD-avg and KD-all remains the same. (c):
the clients’ and generator’s samples.

###  5.2 Evaluation on Real Image Datasets

We evaluated our theoretical results on standard image datasets including
CIFAR-10, CIFAR-100 (Krizhevsky et al., 2009), and ImageNet-32 (Deng et al.,
2009, Chrabaszcz et al., 2017). In our experiments, we simulated heterogeneous
federated learning experiments consisting of non-i.i.d. data at different
clients: for CIFAR-10, we considered 10 clients, each owning 5000 samples
exclusively from a single class of the image dataset. Therefore, every
client’s dataset contains 5000 images having the same label. We put the
experiment on CIFAR-100 and ImageNet in Appendix.

Figure 2: Left: Images generated by a generative model obtain a lower FID-all.
Right: Images from real datasets with class ’plane’ obtain a lower FID-avg.
FID-avg and FID-all lead to inconsistent rankings, while KID-avg and KID-all
result in the same ranking.

Neural Net-based Generators. We trained WGAN-GP and DDPM models using the
standard training protocols outlined in Gulrajani et al. (2017) and Ho et al.
(2020) under the federated CIFAR-10 setting described above. To train the
generative models in the federated learning setting, we employed the standard
FedAvg approach McMahan et al. (2017). Detailed information regarding the
experiment setting can be found in the Appendix. We trained the WGAN-GP
generative models multiple times using different random states, and we set
different training lengths for every training procedure. We saved the models
at different checkpoints every $10$ epochs, which is common in training
generative models to select the best-performing saved model according to an
evaluation metric such as FID and KID.

Perfect Data-simulating Generators. In our CIFAR-10 experiments, we also
simulated and evaluated an "ideal generator" capable of perfectly replicating
all samples belonging to the ’plane’ class in CIFAR10. In this scenario, the
samples "generated" by the ideal generator exhibit impeccable fidelity but
lack diversity since no samples from other categories can be produced.

FID-based and KID-based Evaluation of Generative Models. We evaluated the
described generative models according to FID-all, FID-avg, KID-all, and KID-
avg as defined in Section 4. In several cases, we observed that FID-all / FID-
avg could assign inconsistent rankings to the generators. Specficially, we
computed FID-all and FID-avg for the ideal ’plane’-class-based generator and
neural net-based DDPM generators under the distributed CIFAR10 setting. We
present some examples generated from the two generators in Figure 2 and report
their scores according to the four metrics. The results suggest that FID-avg
assigns a considerably higher score to the ideal ’plane’-based generator,
whose images preserve perfect details but lack diversity in image categories.
Conversely, FID-all assigns a relatively higher value to the DDPM model
because its images possess greater diversity. On the other hand, we also
observed that KID-avg and KID-all give consistent rankings. Both of them led
to the evaluation that the ideal plane generator is slightly better than the
DDPM generator.

To further experiment the ranking of generative models according to the
discussed aggregate scores, we extracted samples from each class of CIFAR-100
and treated them as the output of one hundred distinct generators, each
corresponding to a single class. By assessing these generators on the
federated CIFAR-10 dataset, we obtained one hundred pairs of FID-avg / FID-all
values, and a subset of these pairs with inconsistent rankings according to
FID-all/FID-avg is visualized in the left of Figure 3. The complete set of
evaluation results is available in the Appendix. These results further
highlight that the rankings provided by FID-all and FID-avg can exhibit
inconsistencies in the context of distributed learning. Such inconsistencies
could pose a challenge when selecting from a series of checkpoints or model
architectures during the training of generative models in distributed learning
scenarios, where a distributed computation of FID-all is more challenging than
obtaining FID-avg due to privacy considerations.

Figure 3: FID-based evaluations on federated CIFAR-10 and variance-limited
federated CIFAR-10, FID-avg and FID-all can lead to inconsistent rankings.
Figure 4:  Left: KID evaluations of WGAN-GP checkpoints on federated CIFAR10.
Right: KID Evaluations of variance-controlled generators on variance-limited
federated CIFAR10.

Regarding the KID-based evaluation, our numerical results suggest that the gap
between KID-all and KID-avg remains constant and hence they lead to the same
rankings of the geneartive models. Here, we conducted our evaluations on all
the geneartive models instances as previously described, and the results are
visualized in the left subfigure of Figure 4. These findings reveal that all
distinct generators consistently exhibit a uniform gap between KID-avg and
KID-all. Consequently, our results indicate that the rankings established by
KID-avg consistently align with those of KID-all in distributed learning
settings.

###  5.3 Evaluation on Variance-Limited Federated Datasets

In the federated learning literature, it is relatively common that each client
possesses only a small portion of the collective dataset, and the data
diversity within each client’s holdings is significantly constrained. To
experiment the effect of such distribution heterogeneity, we simulated and
evaluated generative models under _variance-limited federated datasets_. To
obtain a variance-limited federated dataset, for each class in the image
dataset, we kept only a single image and its K-nearest neighbors. To find the
$K$ nearest neighbors, we used the $L_{2}$-distance in the Inception-V3
2048-dimensional semantic space. This approach effectively mimics scenarios
where each client’s data has limited variance. We simulated the variance-
limited federated learning setting for CIFAR-10, CIFAR-100 and a 32$\times$32
version of ImageNet (IN-32) Chrabaszcz et al. (2017). For CIFAR-10 and
CIFAR-100, we utilized all the classes in the dataset and for ImageNet-32 we
utilized the first 100 classes. We chose $K=20$ in the experiments.
Intuitively, a larger $K$ leads to a more significant intra-client variance.

Variance-controlled Generators. To simulate a generator, we initiate the
process by randomly selecting a sample from the dataset. We then gather its
M-nearest neighbors from the original dataset (w/o federated learning
setting). We consider this subset of samples as a set of generated samples
generated by a generator denoted by $G_{M}$. By increasing the value of $M$,
we generated a sequence of generators with progressively higher variance
values. We tried the $M$ from $\\{100,200,500,1000,2000,5000,\\\
10000,15000,20000,25000,30000,35000,40000,45000,50000\\}$.

Numerical Results. We evaluated all the generative models, denoted as $G_{M}$
with the chosen $M$ values, using the Variance-Limited Federated datasets. The
evaluation results on CIFAR10 are in the right subfigure of Figure 3 and
Figure 4. Results on ImageNet-32 are illustrated in Figure 5. Our findings
reveal a distinct pattern in the behavior of FID-avg and FID-all as generator
variance varies while the distance between KID-avg and KID-all remains the
same. The result on CIFAR-100 can be found in the Appendix. Our numerical
results highlight the impact of the choice of FID-all and FID-avg on model
rankings in federated learning settings with limited intra-client variance.

###  5.4 Precision/Recall and Density/Coverege Evaluations

Figure 5:  FID and KID-based Evaluations of variance-controlled generators on
variance-limited federated ImageNet-32.  Figure 6: Precision / Recall and
Density / Coverage evaluation on variance-limited federated CIFAR10.

In addition to FID and KID, we followed the definition in Equation 2 and
Equation 3 and performed similar experiments to evaluate the consistency
between the two aggregate scores for Precision/Recall (Kynkäänniemi et al.,
2019) and Density/Coverage (Naeem et al., 2020).

In the case of Precision/Recall, we utilized the official implementation with
the number of clusters set to 5. We assessed precision and recall under the
variance-limited CIFAR10 setting, with the results presented in Figure 6. The
numerical scores indicate that in a heterogeneous data setting, the two
aggregate precision scores may not consistently rank the generative models. On
the other hand, based on the recall’s definition, it can be seen that Recall-
all and Recall-avg will always take the same value, since the recall score
reduces to an average over the generated data. Our numerical results are also
consistent with this observation.

For the Density/Coverage-based evaluation, the numerical results in Figure 6
suggest that while the density-based aggregate scores lead to more consistent
rankings under data heterogeneity, both density-all and density-avg may still
provide inconsistent rankings in a few cases. On the other hand, the coverage-
based evaluations mimic the recall-based evaluation and consistently rank
models in our experimental settings.

##  6 Conclusion

In this paper, we studied the evaluation of generative models in heterogeneous
distributed learning problems where the clients have different distributions.
We discussed the challenges of evaluating the overall performance of a trained
generative model with only one score and showed the inconsistent rankings of
sensible aggregations of standard FID scores in the network. On the other
hand, we demonstrated that the same extensions of KID offer the same ranking
of generative models. A possible future direction for our work is to extend
the theoretical study to other evaluation criteria such as precision/recall
and density/coverage scores. Also, understanding the behavior of the aggregate
score using non-arithmetic averaging could be useful for evaluating deep
generative models in federated learning contexts.

## Limitations and Broader Impact

We note that our numerical study focuses on the applications of generative
models to image datasets and the empirical conclusions may not apply to other
standard types of data including text and audio data. Regarding the work’s
broader impact, we note that our analysis could be connected to the fairness
evaluation of generative models in distributed learning contexts, as it
suggests evaluation metrics for the assessment of diversity in the generated
data. The study of fairness and diversity for generative models is required
for a principled deployment of generative models in sensitive machine learning
applications.

## References

  * Bińkowski et al. (2018) Mikołaj Bińkowski, Danica J Sutherland, Michael Arbel, and Arthur Gretton.  Demystifying mmd gans.  _arXiv preprint arXiv:1801.01401_ , 2018.
  * Chrabaszcz et al. (2017) Patryk Chrabaszcz, Ilya Loshchilov, and Frank Hutter.  A downsampled variant of imagenet as an alternative to the cifar datasets.  _arXiv preprint arXiv:1707.08819_ , 2017.
  * Deng et al. (2009) Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ , pp. 248–255. Ieee, 2009.
  * Dinh et al. (2016) Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using real nvp.  _arXiv preprint arXiv:1605.08803_ , 2016.
  * Goodfellow et al. (2014) Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  _arXiv preprint arXiv:1406.2661_ , 2014.
  * Gulrajani et al. (2017) Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  _NeurIPS_ , 30, 2017.
  * Hardy et al. (2019) Corentin Hardy, Erwan Le Merrer, and Bruno Sericola.  Md-gan: Multi-discriminator generative adversarial networks for distributed datasets.  In _IPDPS_ , pp. 866–877. IEEE, 2019.
  * Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  _NeurIPS_ , 30, 2017.
  * Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _NeurIPS_ , 33:6840–6851, 2020.
  * Karras et al. (2019) Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In _CVPR_ , pp. 4401–4410, 2019.
  * Kingma & Welling (2013) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Krizhevsky et al. (2009) Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Kynkäänniemi et al. (2019) Tuomas Kynkäänniemi, Tero Karras, Samuli Laine, Jaakko Lehtinen, and Timo Aila.  Improved precision and recall metric for assessing generative models.  _NeurIPS_ , 32, 2019.
  * McMahan et al. (2017) Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas.  Communication-efficient learning of deep networks from decentralized data.  In _Artificial intelligence and statistics_ , pp. 1273–1282. PMLR, 2017.
  * Naeem et al. (2020) Muhammad Ferjad Naeem, Seong Joon Oh, Youngjung Uh, Yunjey Choi, and Jaejun Yoo.  Reliable fidelity and diversity metrics for generative models.  In _ICML_ , pp. 7176–7185. PMLR, 2020.
  * Ramesh et al. (2022) Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.  Hierarchical text-conditional image generation with clip latents.  _arXiv preprint arXiv:2204.06125_ , 2022.
  * Rasouli et al. (2020) Mohammad Rasouli, Tao Sun, and Ram Rajagopal.  Fedgan: Federated generative adversarial networks for distributed data.  _arXiv preprint arXiv:2006.07228_ , 2020.
  * Sajjadi et al. (2018) Mehdi SM Sajjadi, Olivier Bachem, Mario Lucic, Olivier Bousquet, and Sylvain Gelly.  Assessing generative models via precision and recall.  _NeurIPS_ , 31, 2018.
  * Salimans et al. (2016) Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen.  Improved techniques for training gans.  _NeurIPS_ , 29, 2016.
  * Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _ICML_ , pp. 2256–2265. PMLR, 2015.
  * Su et al. (2023) Xuan Su, Jiaming Song, Chenlin Meng, and Stefano Ermon.  Dual diffusion implicit bridges for image-to-image translation.  In _ICLR_ , 2023.
  * Theis et al. (2016) L Theis, A van den Oord, and M Bethge.  A note on the evaluation of generative models.  In _ICLR_ , pp. 1–10, 2016.
  * Wu et al. (2022) Yuezhou Wu, Yan Kang, Jiahuan Luo, Yuanqin He, Lixin Fan, Rong Pan, and Qiang Yang.  Fedcg: Leverage conditional gan for protecting privacy and maintaining competitive performance in federated learning.  In _IJCAI_ , 2022.
  * Yonetani et al. (2019) Ryo Yonetani, Tomohiro Takahashi, Atsushi Hashimoto, and Yoshitaka Ushiku.  Decentralized learning of generative adversarial networks from non-iid data.  _arXiv preprint arXiv:1905.09684_ , 2019.