cell-space-limits = 2pt

# Generalized unistochastic matrices

Ion Nechita<EMAIL_ADDRESS>Laboratoire de Physique Théorique,
Université de Toulouse, CNRS, UPS, France  ,  Zikun Ouyang
<EMAIL_ADDRESS>Institut de Mathématiques, Université de Toulouse,
CNRS, UPS, France  and  Anna Szczepanek<EMAIL_ADDRESS>Institut de Mathématiques, Université de Toulouse, CNRS, UPS, France

###### Abstract.

We introduce and study a class of bistochastic matrices generalizing
unistochastic matrices. Given a complex bipartite unitary operator, we
construct a bistochastic matrix having as entries the normalized squares of
Frobenius norm of the blocks. We show that the closure of the set of
generalized unistochastic matrices is the whole Birkhoff polytope. We
characterize the points on the edges of the Birkhoff polytope that belong to a
given level of our family of sets, proving that the different (non-convex)
levels have a rich inclusion structure. We also study the corresponding
generalization of orthostochastic matrices. Finally, we introduce and study
the natural probability measures induced on our sets by the Haar measure of
the unitary group. These probability measures interpolate between the natural
measure on the set of unistochastic matrices and the Dirac measure supported
on the van der Waerden matrix.

###### Contents

  1. 1 Introduction
  2. 2 Generalized unistochastic matrices
  3. 3 Generalizing the bracelet framework
    1. 3.1 Generalized bracelet conditions
    2. 3.2 Generalized bracelet matrices
  4. 4 Generalized orthostochastic matrices
  5. 5 Random generalized unistochastic matrices

##  1\. Introduction

Bistochastic and unistochastic matrices is a classical topic that comes up
repeatedly in various domains of mathematics and mathematical physics. To set
the stage and introduce notation, let us recall that a square real matrix of
size $d$, the set of which we shall denote by $\mathcal{M}_{d}(\mathbb{R})$,
is called _bistochastic_ (or _doubly stochastic_) if it has non-negative
entries that add up to one in every row and column. That is,
$B=(B_{ij})_{i,j=1}^{d}$ is bistochastic when its entries satisfy the
following conditions:

$\forall{i,j}\quad B_{ij}\geq
0\qquad\forall{i}\quad\sum_{j=1}^{d}B_{ij}=1\quad\text{ and
}\quad\forall{j}\quad\sum_{i=1}^{d}B_{ij}=1.$

We shall denote by $\mathsf{B}_{d}$ the set of all bistochastic matrices of
order $d$.

A bistochastic matrix is called _unistochastic_ when its entries are the
squared absolute values of some unitary matrix of the same size. More
formally, we consider the map

$\begin{aligned}
\varphi_{d}\colon\mathcal{U}(d)\ni\left(U_{ij}\right)_{i,j=1}^{d}\longmapsto\left(|U_{ij}|^{2}\right)_{i,j=1}^{d}\in\mathcal{M}_{d}(\mathbb{R})\end{aligned}.$

The image of $\mathcal{U}(d)$ under $\varphi_{d}$ constitutes the set
$\mathsf{U}_{d}$ of unistochastic matrices of order $d$. Alternatively, using
the Hadamard (entrywise) product $\circ$ of matrices, we can write

$\mathsf{U}_{d}:=\varphi_{d}(\mathcal{U}(d))=\\{U\circ\bar{U}\>|\>U\in\mathcal{U}(d)\\}.$

By unitarity, we have $\mathsf{U}_{d}\subseteq\mathsf{B}_{d}$.

One of the reasons behind the prominence of bistochastic matrices is the fact
that its entries can be regarded as the probabilities that some (classical)
physical system evolves from one state to another. If the bistochastic matrix
is also unistochastic, then the system under consideration can be quantized.
There are many references related to the applications of unistochastic
matrices in various areas, e.g., in quantum information theory and in particle
physics, see [Ben04] and references therein.

It is well known that the _Birkhoff polytope_ $\mathsf{B}_{d}$ is convex and
compact. The extreme points of $\mathsf{B}_{d}$ are permutation matrices, so
the Birkhoff polytope has $d!$ vertices. A bistochastic matrix lies at the
boundary of the Birkhoff polytope iff it has a zero entry. There are $d^{2}$
faces and they correspond to the inequalities that the matrix entries must
satisfy. For instance, consider $d=3$ and

$B=\begin{bmatrix}a&b&*\\\ c&d&*\\\ *&*&*\end{bmatrix}$

with $a,b,c,d\in\mathbb{R}$. The $9$ inequalities corresponding to the faces
of the Birkhoff polytope are

$\begin{cases}a,b,c,d\geq 0\\\ a+b\leq 1,\quad c+d\leq 1,\quad a+c\leq 1,\quad
b+d\leq 1\\\ (1-a-c)+(1-b-d)\leq 1\end{cases}$

An obvious example of a unistochastic matrix is a permutation matrix. Another
well-known example is the van der Waerden (flat) matrix, i.e., the matrix
whose entries are all equal to $1/d$. It is worth noting that the unitary
matrices that induce the van der Waerden matrix are precisely the renowned
complex Hadamard matrices, for instance the Fourier matrix
$\frac{1}{\sqrt{d}}(\exp(\frac{2\pi\mathrm{i}}{d}jk))_{j,k=0}^{d-1}$.

One immediately sees that for $d=2$ every bistochastic matrix is
unistochastic, i.e., we have $\mathsf{U}_{2}=\mathsf{B}_{2}$. This, however,
is a sole exception as for every dimension $d$ higher than two we have
$\mathsf{U}_{d}\subsetneq\mathsf{B}_{d}$ and $\mathsf{U}_{d}$ is known to be
non-convex. A lot of effort has been put into characterizing unistochastic
matrices and one of the key tools turned out to be the _bracelet condition_.
It allows us to distinguish the set of _bracelet matrices_ , which is a
superset of unistochastic matrices. Namely, let
$\alpha=(\alpha_{1},\ldots,\alpha_{d})$ and
$\beta=(\beta_{1},\ldots,\beta_{d})$ be probability vectors (i.e., in each
vector the entries are non-negative and sum up to one). We say that
$(\alpha,\beta)$ satisfies the bracelet condition if

$2\max_{j=1,\ldots,d}\sqrt{\alpha_{j}\beta_{j}}\leq\sum_{j=1}^{d}\sqrt{\alpha_{j}\beta_{j}}.$
(1)

Now, a bistochastic matrix is said to be a _bracelet matrix_ if every pair of
its rows and every pair of its columns, regarded as pairs of probability
vectors, satisfies the bracelet condition; we shall denote the set of bracelet
matrices of order $d$ by $\mathsf{L}_{d}$. The bracelet condition plays an
instrumental role in the study of unistochastic matrices because it
characterizes the first non-trivial case $d=3$ [AYP79], i.e.,

$\mathsf{U}_{3}=\mathsf{L}_{3}\subsetneq\mathsf{B}_{3}$

and, as we already mentioned, it provides a necessary (but not sufficient)
condition for unistochasticity in higher dimensions (see [RMKPŻ22]):

$\forall d\geq
4,\qquad\mathsf{U}_{d}\subsetneq\mathsf{L}_{d}\subsetneq\mathsf{B}_{d}.$

In the present paper we introduce the notion of _generalized unistochastic
matrices_ , denoted by $\mathsf{U}_{d,s}$; here, $s$ is an integer parameter.
The idea is to replace $\mathcal{U}(d)$ with $\mathcal{U}(ds)$ and regard a
unitary matrix from $\mathcal{U}(ds)$ as a $d\times d$ matrix consisting of
$s\times s$ submatrices (blocks). Then it suffices to replace the absolute
values of entries by the normalized squares of Frobenius (or Schatten-2) norms
of blocks to arrive at a $d\times d$ bistochastic matrix again. Formally, we
consider the map

$\varphi_{d,s}\colon\mathcal{U}(ds)\ni\big{(}U_{ij}(k,l)\big{)}_{\begin{subarray}{c}1\leq
i,j\leq d\\\ 1\leq k,l\leq s\end{subarray}}\
\longmapsto\big{(}\tfrac{1}{s}||U_{ij}||_{F}^{2}\big{)}_{1\leq i,j\leq
d}\in\mathsf{B}_{d},$

where $U_{ij}(k,l)$ is the $(k,l)$-th entry in the $(i,j)$-th block, and we
define $\mathsf{U}_{d,s}$ as the image of $\varphi_{d,s}$. Let us recall that
the Frobenius norm of a square complex matrix
$X\in\mathcal{M}_{n}(\mathbb{C})$ is given by

$\|X\|_{F}:=\Big{(}\sum_{i,j=1}^{n}|X_{ij}|^{2}\Big{)}^{1/2}=\Tr(XX^{*})^{1/2}.$

In Propositions 2.4 & 2.5 we show that generalized unistochastic matrices do
indeed generalize the notion of unistochastic matrices, i.e., for every $s$ we
have

$\mathsf{U}_{d}\subseteq\mathsf{U}_{d,s}\subseteq\mathsf{B}_{d}.$

One of the main results of the present paper is Theorem 2.8, where we show
that every bistochastic matrix can be arbitrarily well approximated by a
generalized unistochastic matrix of some order. The key ingredient in proving
this result is the convexity-type property of generalized unistochastic
matrices:

$\forall{s,t}\quad\tfrac{s}{s+t}\mathsf{U}_{d,s}+\tfrac{t}{s+t}\mathsf{U}_{d,t}\subseteq\mathsf{U}_{d,s+t},$

see Proposition 2.6. Then in Corollary 2.7 we investigate further non-trivial
inclusion relations between the sets of generalized unistochastic matrices of
different orders.

Let us point our that an alternative generalization of unistochastic matrices
was proposed by Gutkin in [Gut13]. Unfortunately, as we show at the end of
Section 2, the proposed generalization yields only stochastic (and generally
not bistochastic) matrices, so it is quite far from the usual unistochastic
matrices.

The main tool we develop to investigate $\mathsf{U}_{d,s}$ is the _generalized
bracelet condition_. A pair of probability vectors $\alpha,\beta$ is said to
satisfy the generalized bracelet condition of order $s$ if they correspond to
the normalized squares of Frobenius norms of the blocks of some unitary matrix
$U\in\mathcal{U}(ds)$, i.e., if there exist matrices
$A_{1},\ldots,A_{d},B_{1},\ldots,B_{d}\in\mathcal{M}_{s}(\mathbb{C})$
satisfying

$\forall i\quad\tfrac{1}{s}\|A_{i}\|^{2}_{F}=\alpha_{i}\quad\text{ and
}\quad\tfrac{1}{s}\|B_{i}\|^{2}_{F}=\beta_{i}$

as well as

$\sum_{i=1}^{d}A_{i}A_{i}^{*}=\sum_{i=1}^{d}B_{i}B_{i}^{*}=I_{s},\quad\sum_{i=1}^{d}A_{i}B_{i}^{*}=0,$

which means that the $2$-row block matrix

$\begin{bmatrix}A_{1}&\cdots&A_{d}\\\ B_{1}&\cdots&B_{d}\end{bmatrix}$

(of size $2s\times ds$) can be expanded to a unitary matrix of size $ds\times
ds$. See Proposition 3.4 for the proof that these conditions do indeed
generalize the standard bracelet condition (1).

Since the generalized unistochastic matrices $\mathsf{U}_{d,s}$ are defined as
the image of $\mathcal{U}(ds)$ under $\varphi_{d,s}$, it is natural to equip
$\mathsf{U}_{d,s}$ with the probability measure obtained by pushing forward
the Haar measure from $\mathcal{U}(ds)$ via $\varphi_{d,s}$. We compute the
first few joint moments of the elements of a random matrix
$B\in\mathsf{U}_{d,s}$. In particular, the expected value of $B_{ij}$ equals
$\ 1/d$, while its variance decreases as $s$ grows (and $d$ is fixed), which
means that the probability distribution on $\mathsf{U}_{d,s}$ tends to
concentrate around the van der Waerden matrix. We also draw some conclusions
regarding the covariance and correlation of the elements of $B$.

The paper is organized as follows. In Section 2 we introduce generalized
unistochastic matrices and present their basic properties. Section 3 contains
a suitable generalization of the bracelet conditions for unistochastic
matrices; these conditions allow us to showcase the complexity of the
different levels of the generalized unistochastic sets. In Section 4 we
discuss the corresponding generalizations of orthostochastic matrices.
Finally, in Section 5 we explore the properties of the probability measures
induced on the set of generalized unistochastic matrices by the Haar
distribution on the unitary group.

##  2\. Generalized unistochastic matrices

The main idea of this work can be summarized in the following table:

Source | Operation | Result
---|---|---
| Unitary group
---
$U\in\mathcal{U}(d)$
$\mathbb{C}\ni u_{ij}\mapsto|u_{ij}|^{2}$ |  | Unistochastic matrix
---
$B\in\mathsf{U}_{d}$
| Larger unitary group
---
$U\in\mathcal{U}(ds)\subseteq\mathcal{M}_{d}(\mathcal{M}_{s}(\mathbb{C}))$
$\mathcal{M}_{s}(\mathbb{C})\ni U_{ij}\mapsto\frac{1}{s}||U_{ij}||_{F}^{2}$ |  | Generalized unistochastic matrix
---
$B\in\mathsf{U}_{d,s}$

Let $d\geq 2$ and $s\geq 1$ be integers. In what follows we regard
$B\in\mathcal{M}_{ds}(\mathbb{C})$ as a $d\times d$ block matrix consisting of
$s\times s$ blocks. We shall write $B_{ij}$ for the the $(i,j)$-th block and
$B_{ij}(k,l)$ for the $(k,l)$-th coefficient inside this block, where
$k,l\in[s]$ and $i,j\in[d]$. For brevity, we put $[n]:=\\{1,2,\ldots,n\\}$ and
$\mathfrak{S}_{n}$ for the group of permutations of $[n]$.

######  Definition 2.1.

Consider the map

$\displaystyle\varphi_{d,s}\colon$ $\displaystyle\mathcal{U}(ds)$
$\displaystyle\xlongrightarrow{}$ $\displaystyle\ \mathcal{M}_{d}(\mathbb{R})$
$\displaystyle\left(U_{ij}(k,l)\right)_{i,j\in[d];k,l\in[s]}$
$\displaystyle\longmapsto$
$\displaystyle\>\big{(}\tfrac{1}{s}||U_{ij}||_{F}^{2}\big{)}_{i,j\in[d]}$

We define $\mathsf{U}_{d,s}:=\varphi_{d,s}(\mathcal{U}(ds))$ to be the set of
_generalized unistochastic matrices_.

######  Example 2.2.

For $d=s=2$,
$P=\bNiceMatrix[c,margin]\CodeBefore\rectanglecolor{blue!10}{1-1}{2-2}\rectanglecolor{red!10}{3-1}{4-2}\rectanglecolor{red!10}{1-3}{2-4}\rectanglecolor{blue!10}{3-3}{4-4}\Body
0&001\\\ 0010\\\ 0100\\\
1000\stackrel{{\scriptstyle\varphi_{2,2}}}{{\longmapsto}}\bNiceMatrix[c,margin]\CodeBefore\chessboardcolors{blue!10}{red!10}\Body
0&1\\\ 10\in\mathsf{B}_{2}.$

Importantly, there exist generalized unistochastic matrices which are not
unistochastic. This makes the definition above interesting and justifies the
study of generalized unistochastic matrices.

######  Example 2.3.

For $d=3$,
$P=\bNiceMatrix[c,margin]\CodeBefore\rectanglecolor{blue!10}{1-1}{2-2}\rectanglecolor{red!10}{1-3}{2-4}\rectanglecolor{blue!10}{1-5}{2-6}\rectanglecolor{red!10}{3-1}{4-2}\rectanglecolor{blue!10}{3-3}{4-4}\rectanglecolor{red!10}{3-5}{4-6}\rectanglecolor{blue!10}{5-1}{6-2}\rectanglecolor{red!10}{5-3}{6-4}\rectanglecolor{blue!10}{5-5}{6-6}\Body
0&00010\\\ 000100\\\ 100000\\\ 000001\\\ 010000\\\
001000\stackrel{{\scriptstyle\varphi_{3,2}}}{{\longmapsto}}B=\bNiceMatrix[c,margin]\CodeBefore\chessboardcolors{blue!10}{red!10}\Body
0&\frac{1}{2}\frac{1}{2}\\\ \frac{1}{2}0\frac{1}{2}\\\
\frac{1}{2}\frac{1}{2}0\in\mathsf{U}_{3,2}\setminus\mathsf{U}_{3}$.

Indeed, $P$ is a permutation matrix corresponding to the permutation
$(1\,3\,6\,4\,2\,5)\in\mathfrak{S}_{6}$, hence $B\in\mathsf{U}_{3,2}$.
However, $B\notin\mathsf{U}_{3}$ is not unistochastic since it does not
satisfy the bracelet conditions, see Section 3.

In the next two propositions, we show that generalized unistochastic matrices
are bistochastic and that they contain, for every value of the parameter $s$,
the set of (usual) unistochastic matrices $\mathsf{U}_{d}$, which coincides
with the generalized family at $s=1$.

######  Proposition 2.4.

For every $d\geq 2$ and $s\geq 1$, we have
$\mathsf{U}_{d,s}\subseteq\mathsf{B}_{d}$.

###### Proof.

Let $d\geq 2$ and $s\geq 1$. By unitarity,

$\sum_{i=1}^{d}\sum_{k=1}^{s}|U_{ij}(k,l)|^{2}=\sum_{j=1}^{d}\sum_{l=1}^{s}|U_{ij}(k,l)|^{2}=1.$

Therefore, for all $i\in[d]$ we have

$\sum_{j=1}^{d}B_{ij}=\tfrac{1}{s}\sum_{j=1}^{d}\|U_{ij}\|_{F}^{2}=\tfrac{1}{s}\sum_{k=1}^{s}\sum_{j=1}^{d}\sum_{l=1}^{s}|U_{ij}(k,l)|^{2}=1$

and, analogously, $\sum_{i=1}^{d}B_{ij}=1$, which concludes the proof. ∎

######  Proposition 2.5.

For every $d\geq 2$ and $s\geq 1$, we have
$\mathsf{U}_{d}=\mathsf{U}_{d,1}\subseteq\mathsf{U}_{d,s}$.

###### Proof.

Let $d\geq 2$ and $s\geq 1$, and let $B\in\mathsf{U}_{d}$, i.e., there exists
$U\in\mathcal{U}(d)$ such that $B=\varphi_{d}(U)$. Consider $V:=U\otimes
I_{s}\in\mathcal{U}(ds)$. Then for all $i,j\in[d]$ we have
$V_{ij}=U_{ij}\otimes I_{s}$, which implies that
$\tfrac{1}{s}\|V_{ij}\|_{F}^{2}=|U_{ij}|^{2}=B_{ij}$; hence,
$B\in\mathsf{U}_{d,s}$, as desired. ∎

Next, we show that the sets $\mathsf{U}_{d,s}$ satisfy a kind of convexity
property. This result will be key in showing one of our main results, Theorem
2.8.

######  Proposition 2.6.

For every $d\geq 2$ and $s,t\geq 1$, we have
$\frac{s}{s+t}\mathsf{U}_{d,s}+\frac{t}{s+t}\mathsf{U}_{d,t}\subseteq\mathsf{U}_{d,s+t}$.

###### Proof.

Fix $d\geq 2$ and $s,t\geq 1$, and let $B\in\mathsf{U}_{d,s}$ and
$C\in\mathsf{U}_{d,t}$. There exist $V\in\mathcal{U}(ds)$ and
$W\in\mathcal{U}(dt)$ such that $B_{ij}=\frac{1}{s}\|V_{ij}\|_{F}^{2}$ and
$C_{ij}=\frac{1}{t}\|W_{ij}\|_{F}^{2}$. Consider
$U\in\mathcal{M}_{d(s+t)}(\mathbb{C})$ defined as

$U_{ij}:=\begin{bmatrix}V_{ij}&0\\\ 0&W_{ij}\end{bmatrix}.$

In particular, up to a permutation of blocks, $U$ coincides with $V\oplus W$.
Thus, $U\in\mathcal{U}(d(s+t))$ and

$\tfrac{1}{s+t}\|U_{ij}\|_{F}^{2}=\tfrac{1}{s+t}(\|V_{ij}\|_{F}^{2}+\|W_{ij}\|_{F}^{2})=\tfrac{1}{s+t}(sB_{ij}+tC_{ij}).$

That is,
$\varphi_{d,t+s}(U)=\frac{s}{s+t}B+\frac{t}{s+t}C\in\mathsf{U}_{d,s+t}$, as
desired.

∎

######  Corollary 2.7.

Let $d\geq 2$. From Proposition 2.6 we easily conclude that

  1. (1)

For all orders $s_{1},\ldots,s_{k}\geq 1$, we have

$\frac{s_{1}}{s_{1}+\ldots+s_{k}}\mathsf{U}_{d,s_{1}}+\ldots+\frac{s_{k}}{s_{1}+\ldots+s_{k}}\mathsf{U}_{d,s_{k}}\subseteq\mathsf{U}_{d,s_{1}+\ldots+s_{k}}.$

  2. (2)

For all $s,n\geq 1$, we have $\mathsf{U}_{d,s}\subseteq\mathsf{U}_{d,ns}$.

  3. (3)

For all $s,t\geq 1$, we have
$\mathsf{U}_{d,s}\cap\mathsf{U}_{d,t}\subseteq\mathsf{U}_{d,s+t}$.

In relation to the second point of the corollary above, note that, in general,
we do not have

$s\leq t\implies\mathsf{U}_{d,s}\subseteq\mathsf{U}_{d,t},$

see Corollary 3.13 for counterexamples in this direction.

We now prove the main theorem of this section: the closed union of all
generalized unistochatic matrices constitutes the whole set of bistochastic
matrices.

######  Theorem 2.8.

For every dimension $d\geq 2$, we have

$\overline{\bigcup_{s\geq 1}\mathsf{U}_{d,s}}=\mathsf{B}_{d}.$

###### Proof.

We only need to prove the “$\supseteq$” inclusion. Fix $d\geq 2$ and
$\varepsilon>0$. Let $B\in\mathsf{B}_{d}$. We shall construct $N\in\mathbb{N}$
and $B_{\varepsilon}\in\mathsf{U}_{d,N}$ such that
$\|B-B_{\varepsilon}\|_{F}\leq\varepsilon$. As a bistochastic matrix, $B$ can
be written as a convex combination of permutation matrices, i.e., there exists
a family $\\{t_{\sigma}\,|\,\sigma\in\mathfrak{S}_{d}\\}$ of non-negative
coefficients such that $\sum_{\sigma\in\mathfrak{S}_{d}}t_{\sigma}=1$ and
$B=\sum_{\sigma\in\mathfrak{S}_{d}}t_{\sigma}P_{\sigma}$, where $P_{\sigma}$
is the permutation matrix corresponding to $\sigma\in\mathfrak{S}_{d}$. Take
$\delta:={\varepsilon}/[{2(d!-1)\sqrt{d}}]$ and let $N$ be large enough so
that

$\max_{\sigma\neq id}(t_{\sigma}-\tfrac{k_{\sigma}}{N})\leq\delta,$

where $k_{\sigma}:=\lfloor Nt_{\sigma}\rfloor$. Define
$k_{id}:=N-\sum_{\sigma\neq id}\ {k_{\sigma}}$. Then

$\tfrac{k_{id}}{N}-t_{id}=\sum_{\sigma\neq
id}(t_{\sigma}-\tfrac{k_{\sigma}}{N})\in[0,\delta(d!-1)].$

Let us now consider
$B_{\varepsilon}:=\sum_{\sigma\in\mathfrak{S}_{d}}\frac{k_{\sigma}}{N}P_{\sigma}$.
Since $P_{\sigma}\in\mathsf{U}_{d}\subset\mathsf{U}_{d,k_{\sigma}}$ if
$k_{\sigma}\neq 0$, from Corollary 2.7 it follows that
$B_{\varepsilon}\in\mathsf{U}_{d,N}$. Therefore,

$\|B-B_{\varepsilon}\|_{F}=\Big{\|}\sum_{\sigma\in\mathfrak{S}_{d}}(t_{\sigma}-\tfrac{k_{\sigma}}{N})P_{\sigma}\Big{\|}_{F}\leq\sum_{\sigma\in\mathfrak{S}_{d}}|t_{\sigma}-\tfrac{k_{\sigma}}{N}|\cdot\sqrt{d}\leq
2(d!-1)\delta\sqrt{d}=\varepsilon,$

as desired. ∎

Next, we present some numerical simulations regarding the set
$\mathsf{U}_{3,1}$ of $3\times 3$ unistochastic matrices and its generalized
version $\mathsf{U}_{3,2}$, see Figure 1. To decide whether a bistochastic
matrix $B$ is an element of $\mathsf{U}_{d,s}$, we use the NMinimize function
of Wolfram Mathematica to try finding a unitary matrix $U\in\mathcal{U}(ds)$
such that $\varphi_{d,s}(U)=B$. This method does not guarantee finding the
global minimum of non-convex functions, so the results in Figure 1 are
empirical; note however that there is a perfect fit with the theory in the
case $(d,s)=(3,1)$.

Figure 1. In the simplex defined by the identity, $(123)$, $(321)$ permutation
matrices, we plot 1000 uniformly sampled random bistochastic matrices. In the
left panel, we plot in blue unistochastic elements, i.e., samples
$B\in\mathsf{U}_{3,1}$, and in red samples outside this set. In the right
panel, we use the same colours to plot the elements inside and outside the set
$\mathsf{U}_{3,2}$. The gray curves correspond to the bracelet conditions (4)
characterizing unistochastic $3\times 3$ matrices, see Section 3.

We end this section by discussing Gutkin’s generalization of unistochastic
matrices [Gut13]. In that paper, the author generalizes unistochastic matrices
starting from an isometry

$V\colon\mathbb{C}^{d}\to\mathbb{C}^{d}\otimes\mathbb{C}^{s}.$

This isometry can be seen as a $\mathbb{C}^{s}$-valued $d\times d$ matrix:

$\mathcal{M}_{ds\times d}(\mathbb{C})\cong\mathcal{M}_{d}(\mathbb{C}^{s}).$

Denoting the vector elements of $V$ as $v_{ij}\in\mathbb{C}^{s}$, where
$i,j\in[d]$, Gutkin defines

$B_{ij}:=\|v_{ij}\|^{2}.$

Unfortunately, in general the resulting matrix $B$ is then only column
stochastic, and not row stochastic. This can be seen, for example, by
considering the case $d=s=2$ and the isometry

$V=\begin{bmatrix}1&0\\\ 0&1\\\ 0&0\\\ 0&0\end{bmatrix}$

which corresponds to the vectors

$v_{11}=\big{[}\begin{smallmatrix}1\\\ 0\end{smallmatrix}\big{]},\quad
v_{12}=\big{[}\begin{smallmatrix}0\\\ 1\end{smallmatrix}\big{]},\quad
v_{21}=v_{22}=\big{[}\begin{smallmatrix}0\\\ 0\end{smallmatrix}\big{]},$

which, in turn, lead to the matrix

$B=\begin{bmatrix}1&1\\\ 0&0\end{bmatrix}.$

The error in [Gut13] seems to stem from Lemma 1 being wrong, and thus
equations (1) and (2) in that paper not being equivalent.

##  3\. Generalizing the bracelet framework

In this section we generalize, in the same spirit as the main Definition 2.1,
the notions of _bracelet conditions_ and _bracelet matrices_ , ideas
originating in [AYP79]. These notions play an important role in the study of
unistochastic matrices, since bracelet conditions fully characterize the first
non-trivial case, that of dimension $d=3$. Indeed, a $3\times 3$ bistochastic
matrix is unistochastic if and only if it satisfies the bracelet conditions
[AYP79], i.e. iff it is a bracelet matrix.

We first review the standard notion of bracelet condition. The intuitive idea
behind it is that the elements of two rows of a unistochastic matrix
corresponding to the same column cannot be too large simultaneously with
respect to the other row elements. This is because the scalar product of the
corresponding rows of the unitary matrix needs to be zero, so each individual
term of the sum cannot be too large in magnitude with respect to the others.
This intuition is encoded in the following definition:

$\mathsf{Brac}_{d}:=\bigg{\\{}(\alpha,\beta)\in\Delta_{d}^{2}\,:\,\forall
i\in[d],\,\sqrt{\alpha_{i}\beta_{i}}\leq\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{d}\sqrt{\alpha_{j}\beta_{j}}\bigg{\\}}.$ (2)

In the formula above, recall that the $(d-1)$-dimensional probability simplex
is the set of all probability vectors in $\mathbb{R}^{d}$:

$\Delta_{d}:=\Big{\\{}\alpha\in\mathbb{R}^{d}\,:\,\forall
i\in[d],\,\alpha_{i}\geq 0\text{ and }\sum_{i=1}^{d}\alpha_{i}=1\Big{\\}}.$

We have the following important definition; note that the term “bracelet
matrix / condition” was introduced in [RMKPŻ22].

######  Definition 3.1.

A bistochastic matrix $B\in\mathsf{B}_{d}$ is said to be a _bracelet matrix_
if all pairs of rows and all pairs of columns of $B$ satisfy the bracelet
condition from (2). We introduce the set of bracelet matrices

$\mathsf{L}_{d}:=\Big{\\{}B\in\mathsf{B}_{d}:\,\forall i_{1}\neq
i_{2}\,(B_{i_{1},\cdot},B_{i_{2},\cdot})\in\mathsf{Brac}_{d}\,\text{ and
}\,\forall j_{1}\neq
j_{2}\,(B_{\cdot,j_{1}},B_{\cdot,j_{2}})\in\mathsf{Brac}_{d}\big{\\}}.$ (3)

It was observed in [RMKPŻ22] that being bracelet is a necessary condition for
unistochasticity. We give here the proof of this claim for the sake of
completeness.

######  Proposition 3.2.

For all dimensions $d\geq 2$, we have $\mathsf{U}_{d}\subseteq\mathsf{L}_{d}$.

###### Proof.

Let $B=(B_{ij})_{i,j}\in\mathsf{U}_{d}$ be such that $B_{ij}=|U_{ij}|^{2}$,
where $U\in\mathcal{U}(d)$ is a corresponding unitary matrix. Fix two row
indices $i_{1}\neq i_{2}$. By unitarity, for all $k\in[d]$ we have

$-U_{i_{1}k}\overline{U_{i_{2}k}}=\sum\limits_{\begin{subarray}{c}j=1\\\ j\neq
k\end{subarray}}^{d}U_{i_{1}j}\overline{U_{i_{2}j}}.$

Taking norm and applying the triangle inequality, for all $k\in[d]$ we obtain

$|U_{i_{1}k}\overline{U_{i_{2}k}}|\leq\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq k\end{subarray}}^{d}|U_{i_{1}j}\overline{U_{i_{2}j}}|,$

which translates into

$\sqrt{B_{i_{1}k}B_{i_{2}k}}\leq\sum\limits_{\begin{subarray}{c}j=1\\\ j\neq
k\end{subarray}}^{d}\sqrt{B_{i_{1}j}B_{i_{2}j}}.$

A similar computation shows that the columns of $B$ also satisfy the bracelet
condition from Eq. (2); hence, $B\in\mathsf{L}_{d}$, as claimed. ∎

The bracelet conditions characterize unistochasticity for $3\times 3$ matrices
(i.e., $\mathsf{U}_{3}=\mathsf{L}_{3}$, see [AYP79]), while being only
necessary for $d\geq 4$, see [RMKPŻ22]. Note that the complete description of
the non-convex set $\mathsf{U}_{3}$ was obtained, thanks to the
characterization in terms of bracelet conditions, in [Nak96]. For example, for
the bistochastic matrices studied in Figure 1, which are of the form

$B=\begin{bmatrix}\lambda_{1}&\lambda_{2}&\lambda_{3}\\\
\lambda_{3}&\lambda_{1}&\lambda_{2}\\\
\lambda_{2}&\lambda_{3}&\lambda_{1}\end{bmatrix},$

the bracelet conditions read

$\sqrt{\lambda_{i}\lambda_{j}}\leq\sqrt{\lambda_{i}\lambda_{k}}+\sqrt{\lambda_{j}\lambda_{k}},$
(4)

for any permutation $(i,j,k)$ of the set $\\{1,2,3\\}$. The matrices
satisfying these conditions are precisely the unistochastic matrices, and they
correspond to the region delimited by the gray curves in Figure 1.

###  3.1. Generalized bracelet conditions

Since the idea behind the bracelet condition from Eq. (2) was to use the
orthogonality of the rows/columns of a unitary operator, we generalize this
insight to our setting in the following definition, encoding in it the block-
orthogonality of block unitary matrices.

######  Definition 3.3.

A pair of probability vectors $\alpha,\beta\in\Delta_{d}$ is said to satisfy
the _generalized bracelet condition_ of order $s$ if they correspond to the
normalized squares of Frobenius norms of the blocks of some unitary matrix
$U\in\mathcal{U}(ds)$:

$\displaystyle\mathsf{Brac}_{d,s}:=\Big{\\{}(\alpha,\beta)\in\Delta_{d}^{2}\,:\,$
$\displaystyle\exists
A_{1},\ldots,A_{d},B_{1},\ldots,B_{d}\in\mathcal{M}_{s}(\mathbb{C})\text{ such
that }$
$\displaystyle\sum_{i=1}^{d}A_{i}A_{i}^{*}=\sum_{i=1}^{d}B_{i}B_{i}^{*}=I_{s},\quad\sum_{i=1}^{d}A_{i}B_{i}^{*}=0_{s},$
(5) $\displaystyle\tfrac{1}{s}\|A_{i}\|_{F}^{2}=\alpha_{i}\>\text{ and }\
\tfrac{1}{s}\|B_{i}\|_{F}^{2}=\beta_{i}\quad\forall i\in[d]\Big{\\}}.$

The conditions above mean precisely that we can expand the 2-row block matrix

$\begin{bmatrix}A_{1}&\cdots&A_{d}\\\ B_{1}&\cdots&B_{d}\end{bmatrix}$

to a full unitary matrix $U\in\mathcal{U}(ds)$. Note that the case $d=1$ of
the bracelet conditions is empty for every order, i.e.,
$\mathsf{Brac}_{1,s}=\emptyset$ for every $s\geq 1$; indeed, if
$A_{1}A_{1}^{*}=B_{1}B_{1}^{*}=I_{s}$, then $A_{1},B_{1}\in\mathcal{U}(s)$,
which implies that $A_{1}B_{1}^{*}\neq 0$.

We first show that the newly introduced conditions from Eq. (3.3) do indeed
generalize the standard bracelet conditions from Eq. (2).

######  Proposition 3.4.

For all dimensions $d\geq 2$, the generalized bracelet conditions of order
$s=1$ are precisely the usual bracelet conditions:
$\mathsf{Brac}_{d,1}=\mathsf{Brac}_{d}$.

###### Proof.

For $s=1$ the generalized bracelet condition takes the form

$\mathsf{Brac}_{d,1}=\left\\{\begin{aligned}
&(\alpha,\beta)\in\Delta_{d}^{2}\,:\,\exists
a_{1},\ldots,a_{d},b_{1},\ldots,b_{d}\in\mathbb{C}\text{ such that}\\\
&\sum\nolimits_{i}|a_{i}|^{2}=\sum\nolimits_{i}|b_{i}|^{2}=1,\,\sum\nolimits_{i}a_{i}\bar{b}_{i}=0,\,\text{
and }\,|a_{i}|^{2}=\alpha_{i},\,|b_{i}|^{2}=\beta_{i}\ \ \forall
i\in[d]\end{aligned}\>\right\\}.$

The inclusion $\mathsf{Brac}_{d,1}\subseteq\mathsf{Brac}_{d}$ follows by
mimicking the proof of Proposition 3.2. The converse inclusion can be thought
of as a generalized version of the triangle inequality: for
$l_{1}\geq\ldots\geq l_{d}\geq 0$ satisfying $l_{1}\leq\sum_{j=2}^{d}l_{j}$,
there exist $\theta_{1},\ldots,\theta_{d}\in[0,2\pi)$ such that
$\sum_{j=1}^{d}l_{j}e^{\mathrm{i}\theta_{j}}=0$. Therefore, for any
$(\alpha,\beta)\in\mathsf{Brac}_{d}$ we can choose the phases
$\theta_{1},\ldots,\theta_{d}$ so that
$a_{j}:=\sqrt{\smash[b]{\alpha_{j}}}e^{{\mathrm{i}\theta_{j}}}$ and
$b_{j}:=\sqrt{\smash[b]{\beta_{j}}}$ satisfy
$\sum\nolimits_{j=1}^{d}a_{j}\bar{b}_{j}=0$. The other conditions follow
trivially, and so $(\alpha,\beta)\in\mathsf{Brac}_{d,1}$, as claimed. ∎

As it is the case for the sets $\mathsf{U}_{d,s}$ (see Proposition 2.6), the
sets $\mathsf{Brac}_{d,s}$ satisfy the following “convexity” relation:

######  Proposition 3.5.

For every $d\geq 2$ and $s,t\geq 1$, we have
$\frac{s}{s+t}\mathsf{Brac}_{d,s}+\frac{t}{s+t}\mathsf{Brac}_{d,t}\subseteq\mathsf{Brac}_{d,s+t}.$

###### Proof.

Consider pairs of probability vectors $(\alpha,\beta)\in\mathsf{Brac}_{d,s}$
and $(\mu,\nu)\in\mathsf{Brac}_{d,t}$ together with the generating matrices
$A_{1},\ldots,A_{d},B_{1},\ldots,B_{d}\in\mathcal{M}_{s}(\mathbb{C})$ and
$C_{1},\ldots,C_{d},D_{1},\ldots,D_{d}\in\mathcal{M}_{t}(\mathbb{C})$, i.e.,

$\begin{bmatrix}A_{1}&\cdots&A_{d}\\\
B_{1}&\cdots&B_{d}\end{bmatrix}\xmapsto{\frac{1}{s}\|\cdot\|^{2}_{F}}\begin{bmatrix}\alpha\\\
\beta\end{bmatrix}=\begin{bmatrix}\alpha_{1}&\cdots&\alpha_{d}\\\
\beta_{1}&\cdots&\beta_{d}\end{bmatrix}$ $\begin{bmatrix}C_{1}&\cdots&C_{d}\\\
D_{1}&\cdots&D_{d}\end{bmatrix}\xmapsto{\frac{1}{t}\|\cdot\|^{2}_{F}}\begin{bmatrix}\mu\\\
\nu\end{bmatrix}=\begin{bmatrix}\mu_{1}&\cdots&\mu_{d}\\\
\nu_{1}&\cdots&\nu_{d}\end{bmatrix}$

Then, using a direct sum construction, we have:

$\begin{bmatrix}A_{1}\oplus C_{1}&\cdots&A_{d}\oplus C_{d}\\\ B_{1}\oplus
D_{1}&\cdots&B_{d}\oplus
D_{d}\end{bmatrix}\xmapsto{\frac{1}{s+t}\|\cdot\|^{2}_{F}}\begin{bmatrix}\frac{s\alpha_{1}+t\mu_{1}}{s+t}&\cdots&\frac{s\alpha_{d}+t\mu_{d}}{s+t}\\\\[3.30002pt]
\frac{s\beta_{1}+t\nu_{1}}{s+t}&\cdots&\frac{s\beta_{d}+t\nu_{d}}{s+t}\end{bmatrix}=\begin{bmatrix}\frac{s\alpha+t\mu}{s+t}\\\\[3.30002pt]
\frac{s\beta+t\nu}{s+t}\end{bmatrix}.$

Hence,
$\frac{s}{s+t}(\alpha,\beta)+\frac{t}{s+t}(\mu,\nu)\in\mathsf{Brac}_{d,s+t}$,
as desired. ∎

The result above easily generalizes to more than two summands.

######  Corollary 3.6.

For all dimensions $d\geq 2$ and all orders $s_{1},\ldots,s_{k}\geq 1$, we
have

$\frac{s_{1}}{s_{1}+\ldots+s_{k}}\mathsf{Brac}_{d,s_{1}}+\ldots+\frac{s_{k}}{s_{1}+\ldots+s_{k}}\mathsf{Brac}_{d,s_{k}}\subseteq\mathsf{Brac}_{d,s_{1}+\ldots+s_{k}}.$

In particular, for all $s\geq 1$ we have
$\mathsf{Brac}_{d,1}=\mathsf{Brac}_{d}\subseteq\mathsf{Brac}_{d,s}$.

In what follows, we analyze the sets $\mathsf{Brac}_{d,s}$ for different
values of the dimension $d$ and of the generalization parameter $s$. We will
show that for $d=2$ nothing new happens when the value of $s$ changes, i.e.,
$\mathsf{Brac}_{2,s}=\mathsf{Brac}_{2}$ for every $s\geq 1$; however, for all
dimensions $d\geq 3$ we will obtain strict inclusion
$\mathsf{Brac}_{d,s}\supsetneq\mathsf{Brac}_{d}$ for every $s>1$. Let us start
with an auxiliary lemma.

######  Lemma 3.7.

Consider two real diagonal matrices
$X=\operatorname{diag}(x_{1},\cdots,x_{s})$ and
$Y=\operatorname{diag}(y_{1},\cdots,y_{s})$. The following conditions are
equivalent:

  1. (i)

There exist $U,V\in\mathcal{U}(s)$ such that $XUY=V$.

  2. (ii)

There exists $\sigma\in\mathfrak{S}_{s}$ such that $|x_{i}y_{\sigma(i)}|=1$
for every $i\in[s]$.

###### Proof.

We shall prove the double implication.
$(i)\Leftarrow(ii)$: Let $\sigma\in\mathfrak{S}_{s}$ be such that
$|x_{i}y_{\sigma(i)}|=1$ for all $i\in[s]$. Taking $U$ to be the phase-
permutation matrix corresponding to $\sigma$ and phases
$\operatorname{arg}(x_{i}y_{\sigma(i)})$, we have $XUY=I_{s}$.
$(i)\Rightarrow(ii)$: If $XUY=V$ for some unitary matrices $U$ and $V$, then
the real diagonal matrices $X,Y$ are both invertible; furthermore, we have

$I_{s}=V^{*}V=Y^{*}U^{*}X^{*}XUY=YU^{*}X^{2}UY,$

which implies that

$X^{2}U=UY^{-2}.$

The latter equation means that the $s$ linearly independent column vectors of
the unitary matrix $U$ form $s$ eigenvectors of $X^{2}$, and so the $s$
eigenvalues of $X^{2}$ are exactly the $s$ diagonal elements of $Y^{-2}$.
Therefore, the set of $s$ diagonal elements of $X^{2}$ (counting
multiplicities) is equal to the set of $s$ diagonal elements of $Y^{-2}$.
Hence, there exists $\sigma\in\mathfrak{S}_{s}$ such that
$x_{i}^{2}=1/y_{\sigma(i)}^{2}$ for all $i\in[s]$, which concludes the proof.
∎

######  Proposition 3.8.

For $d=2$ and all orders $s\geq 1$, we have

$\mathsf{Brac}_{2,s}=\mathsf{Brac}_{2}=\left\\{\big{(}(p,1-p),(1-p,p)\big{)}\,:\,p\in[0,1]\right\\}.$

###### Proof.

Fix $s\geq 2$ and let
$\big{(}(\alpha_{1},1-\alpha_{1}),(\beta_{1},1-\beta_{1})\big{)}\in\mathsf{Brac}_{2,s}$.
It suffices to show that $\beta_{1}=1-\alpha_{1}$. Let
$A_{1},A_{2},B_{1},B_{2}\in\mathcal{M}_{s}(\mathbb{C})$ be such that
$\tfrac{1}{s}\Tr(A_{1}A_{1}^{*})=\alpha_{1}$ and
$\tfrac{1}{s}\Tr(B_{1}B_{1}^{*})=\beta_{1}$ as well as

$\displaystyle A_{1}A_{1}^{*}+A_{2}A_{2}^{*}$ $\displaystyle=I_{s}$
$\displaystyle B_{1}B_{1}^{*}+B_{2}B_{2}^{*}$ $\displaystyle=I_{s}$
$\displaystyle A_{1}B_{1}^{*}+A_{2}B_{2}^{*}$ $\displaystyle=0_{s}.$

Since $A_{1}A_{1}^{*}$ and $A_{2}A_{2}^{*}$ commute, there exists
$U\in\mathcal{U}(s)$ such that

$\displaystyle A_{1}A_{1}^{*}$
$\displaystyle=U\operatorname{diag}(a_{1},\ldots,a_{s})U^{*}$ $\displaystyle
A_{2}A_{2}^{*}$
$\displaystyle=U\operatorname{diag}(1-a_{1},\ldots,1-a_{s})U^{*}$

for some real numbers $a_{1},\ldots,a_{s}\in[0,1]$. Moreover, using the
singular value decomposition, there exist $U_{1},U_{2}\in\mathcal{U}(s)$ such
that

$\displaystyle A_{1}$
$\displaystyle=U\operatorname{diag}(\sqrt{a_{1}},\ldots,\sqrt{a_{s}})U_{1}^{*}$
$\displaystyle A_{2}$
$\displaystyle=U\operatorname{diag}(\sqrt{1-a_{1}},\ldots,\sqrt{1-a_{s}})U_{2}^{*}.$

Similarly, there exist $V,V_{1},V_{2}\in\mathcal{U}(s)$ such that

$\displaystyle B_{1}$
$\displaystyle=V\operatorname{diag}(\sqrt{\smash[b]{b_{1}}},\ldots,\sqrt{\smash[b]{b_{s}}})V_{1}^{*}$
$\displaystyle B_{2}$
$\displaystyle=V\operatorname{diag}(\sqrt{\smash[b]{1-b_{1}}},\ldots,\sqrt{\smash[b]{1-b_{s}}})V_{2}^{*}$

with $b_{1},\ldots,b_{s}\in[0,1]$. Denoting $W_{1}:=U_{1}^{*}V_{1}$ and
$W_{2}:=U_{2}^{*}V_{2}$, we have $W_{1},W_{2}\in\mathcal{U}(s)$ and the
condition $A_{1}B_{1}^{*}+A_{2}B_{2}^{*}=0$ now reads

$\displaystyle\operatorname{diag}(\sqrt{a_{1}},\ldots,\sqrt{a_{s}})W_{1}\operatorname{diag}($
$\displaystyle\sqrt{\smash[b]{b_{1}}},\ldots,\sqrt{\smash[b]{b_{s}}})$
$\displaystyle=$
$\displaystyle-\operatorname{diag}(\sqrt{1-a_{1}},\ldots,\sqrt{1-a_{s}})W_{2}\operatorname{diag}(\sqrt{\smash[b]{1-b_{1}}},\ldots,\sqrt{\smash[b]{1-b_{s}}}).$

For non-degenerate $A_{1},A_{2},B_{1},B_{2}$ (i.e., if none of $\alpha_{i}$’s
or $\beta_{i}$’s equals $0$ or $1$), the condition
$A_{1}B_{1}^{*}+A_{2}B_{2}^{*}=0$ can be written as

$\operatorname{diag}\left(\sqrt{\frac{a_{1}}{1-a_{1}}},\ldots,\sqrt{\frac{a_{s}}{1-a_{s}}}\right)W_{1}\operatorname{diag}\left(\sqrt{\frac{b_{1}}{1-b_{1}}},\ldots,\sqrt{\frac{b_{s}}{1-b_{s}}}\right)=-W_{2}.$

Using Lemma 3.7, there exists a permutation $\sigma\in\mathfrak{S}_{s}$ such
that $\frac{a_{i}}{1-a_{i}}\frac{b_{\sigma(i)}}{1-b_{\sigma(i)}}=1$, that is,
$a_{i}=1-b_{\sigma(i)}$, for every $i\in[s]$. Therefore,

$\alpha_{1}=\tfrac{1}{s}\Tr(A_{1}A_{1}^{*})=\tfrac{1}{s}\sum_{i=1}^{s}a_{i}=\tfrac{1}{s}\sum_{i=1}^{s}(1-b_{\sigma(i)})=1-\beta_{1};$

thus, for $p:=\alpha_{1}\in[0,1]$ we have

$\begin{bmatrix}A_{1}&A_{2}\\\
B_{1}&B_{2}\end{bmatrix}\xmapsto{\frac{1}{s}\|\cdot\|^{2}_{F}}\begin{bmatrix}p&1-p\\\
1-p&p\end{bmatrix}.$

For general $A_{1},A_{2},B_{1},B_{2}$, because the general linear group
$\mathcal{GL}_{s}(\mathbb{C})$ is dense in $\mathcal{M}_{s}(\mathbb{C})$ and
the map $\frac{1}{s}\|\cdot\|^{2}_{F}$ is continuous, we again obtain
$\alpha_{1}=1-\beta_{1}$, which finishes the proof. ∎

We now move on to the case $d=3$ and we show that
$\mathsf{Brac}_{3,s}\supsetneq\mathsf{Brac}_{3,1}=\mathsf{Brac}_{3}$ for every
$s\geq 2$. This fact will actually imply that

$\mathsf{Brac}_{d,s}\supsetneq\mathsf{Brac}_{d,1}=\mathsf{Brac}_{d}\ \text{
for all }\ d\geq 3,\,s\geq 2$

since for all $d,n,s$ we have
$\mathsf{Brac}_{d+n,s}\cap\\{\alpha_{d+1}=\cdots=\alpha_{d+n}=\beta_{d+1}=\cdots=\beta_{d+n}=0\\}\cong\mathsf{Brac}_{d,s}.$

Back to the claim about $d=3$, we shall focus on the slice

$\displaystyle\mathsf{Brac}_{3,s}\cap\\{\alpha_{3}=\beta_{2}=0\\}=\big{\\{}$
$\displaystyle\big{(}(\alpha_{1},1-\alpha_{1},0),(\beta_{1},0,1-\beta_{1})\big{)}\in\Delta_{3}^{2}\,\colon\exists
A_{1},A_{2},B_{1},B_{3}\in\mathcal{M}_{s}(\mathbb{C})\text{ s.t.}$
$\displaystyle
A_{1}A_{1}^{*}+A_{2}A_{2}^{*}=B_{1}B_{1}^{*}+B_{3}B_{3}^{*}=I_{s},\,A_{1}B_{1}^{*}=0_{s},$
(6)
$\displaystyle\tfrac{1}{s}\\!\Tr(A_{1}A_{1}^{*})=\alpha_{1},\,\tfrac{1}{s}\\!\Tr(A_{2}A_{2}^{*})=1-\alpha_{1},$
$\displaystyle\tfrac{1}{s}\\!\Tr(B_{1}B_{1}^{*})=\beta_{1},\,\tfrac{1}{s}\\!\Tr(B_{3}B_{3}^{*})=1-\beta_{1}\,\big{\\}}.$

######  Proposition 3.9.

For $d=3$ and every $s\geq 1$, we have

$\mathsf{Brac}_{3,s}\cap\\{\alpha_{3}=\beta_{2}=0\\}=\left\\{\big{(}(\alpha_{1},1-\alpha_{1},0),(\beta_{1},0,1-\beta_{1})\big{)}\in\Delta_{3}^{2}\,:\,\lceil\alpha_{1}s\rceil+\lceil\beta_{1}s\rceil\leq
s\right\\}.$

###### Proof.

Let $s\geq 1$. We shall prove the double inclusion.

“$\subseteq$”: Consider matrices $A_{1},A_{2},B_{1},B_{3}$ as in Eq. (3.1) and
let $\sqrt{a_{1}}\geq\ldots\geq\sqrt{a_{s}}$ be the singular values of
$A_{1}$. The eigenvalues of $A_{1}A_{1}^{*}$ are therefore
$a_{1}\geq\ldots\geq a_{s}\geq 0$, and $A_{1}A_{1}^{*}+A_{2}A_{2}^{*}=I_{s}$
guarantees that $a_{1}\leq 1$. Since
$\operatorname{rk}(A_{1})=\operatorname{rk}(A_{1}A_{1}^{*})=\>\text{the number
of nonzero }a_{i}$’s, it follows that

$\alpha_{1}s=\Tr(A_{1}A_{1}^{*})=\sum_{i=1}^{s}a_{i}=\sum_{i=1}^{\operatorname{rk}(A_{1})}a_{i}\leq\operatorname{rk}(A_{1})a_{1}\leq\operatorname{rk}(A_{1}).$

In consequence, $\lceil\alpha_{1}s\rceil\leq\operatorname{rk}(A_{1})$ and,
similarly, $\lceil\beta_{1}s\rceil\leq\operatorname{rk}(B_{1})$. Finally, from
$A_{1}B_{1}^{*}=0$ we obtain $\operatorname{rk}B_{1}^{*}\leq\dim\ker
A_{1}=s-\operatorname{rk}A_{1}$, and so
$\operatorname{rk}(A_{1})+\operatorname{rk}(B_{1})\leq s$, which proves the
first inclusion.

“$\supseteq$”: Conversely, given $\alpha_{1},\beta_{1}\in[0,1]$ satisfying
$\lceil\alpha_{1}s\rceil+\lceil\beta_{1}s\rceil\leq s$, let us consider

$\displaystyle A_{1}$
$\displaystyle:=\operatorname{diag}\Big{(}\sqrt{a_{1}},\ldots,\sqrt{\smash[b]{a_{\lceil\alpha_{1}s\rceil}}},0,\ldots,0\Big{)}$
$\displaystyle A_{2}$
$\displaystyle:=\operatorname{diag}\Big{(}\sqrt{1-a_{1}},\ldots,\sqrt{\smash[b]{1-a_{\lceil\alpha_{1}s\rceil}}},1,\ldots,1\Big{)}$
$\displaystyle B_{1}$
$\displaystyle:=\operatorname{diag}\Big{(}0,\ldots,0,\sqrt{\smash[b]{b_{1}}},\ldots,\sqrt{\smash[b]{b_{\lceil\beta_{1}s\rceil}}}\Big{)}$
$\displaystyle B_{3}$
$\displaystyle:=\operatorname{diag}\Big{(}1,\ldots,1,\sqrt{\smash[b]{1-b_{1}}},\ldots,\sqrt{\smash[b]{1-b_{\lceil\beta_{1}s\rceil}}}\Big{)},$

where

$a_{1}=\cdots=a_{\lfloor\alpha_{1}s\rfloor}=1$

and, in the case when $\alpha_{1}s$ is not an integer,

$a_{\lceil\alpha_{1}s\rceil}=a_{\lfloor\alpha_{1}s\rfloor+1}=\alpha_{1}s-\lfloor\alpha_{1}s\rfloor\in(0,1).$

The $b_{i}$’s are defined analogously, using $\beta_{1}s$. The assumption
$\lceil\alpha_{1}s\rceil+\lceil\beta_{1}s\rceil\leq s$ guarantees that the
non-zero elements of $A_{1}$ and $B_{1}$ do not overlap, which implies that
$A_{1}B_{1}^{*}=0_{s}$. One can easily verify that all the other conditions
from Eq. (3.1) hold as well, and so the proof is finished. ∎

Let us consider

$E(s):=\left\\{(\alpha_{1},\beta_{1})\in[0,1]^{2}\colon\lceil\alpha_{1}s\rceil+\lceil\beta_{1}s\rceil\leq
s\right\\}.$

The sets $E(2),E(3),E(4),E(5)$ are displayed in Figure 2. Note that every
$E(s)$ contains the axes:

$E(1)=\big{\\{}(\alpha_{1},0)\,:\,\alpha_{1}\in[0,1]\big{\\}}\cup\big{\\{}(0,\beta_{1})\,:\,\beta_{1}\in[0,1]\big{\\}}\subseteq
E(s)\ \text{ for every}\ s\geq 1.$

Moreover, we have

$\big{\\{}(\alpha_{1},\beta_{1})\,:\,\alpha_{1}+\beta_{1}\leq
1-{2}/{s}\big{\\}}\subseteq
E(s)\subseteq\big{\\{}{(\alpha_{1},\beta_{1})\,:\,\alpha_{1}+\beta_{1}\leq
1\big{\\}}}\ \text{ for every}\ s\geq 1;$

hence,

$\overline{\lim_{s\to\infty}E(s)}=\big{\\{}{(\alpha_{1},\beta_{1})\in[0,1]^{2}\,:\,\alpha_{1}+\beta_{1}\leq
1\big{\\}}}.$

This shows that, after taking the limit $s\to\infty$ and the closure, the
bracelet conditions corresponding to the slice considered in Section 3.1
become trivial; this result is in the spirit of Theorem 2.8.

Figure 2. The (interior of the) sets $E(s)$ for $s=2,3,4,5$.

###  3.2. Generalized bracelet matrices

Similarly to Definition 3.1, we introduce the set of generalized bracelet
matrices.

######  Definition 3.10.

A bistochastic matrix $B\in\mathsf{B}_{d}$ is called a _generalized bracelet
matrix_ of order $s$ if all pairs of rows and all pairs of columns of $B$
satisfy the generalized bracelet condition of order $s$ from Eq. (3.3). That
is, the set of generalized bracelet matrices is defined as

$\mathsf{L}_{d,s}:=\big{\\{}B\in\mathsf{B}_{d}:\,\forall i_{1}\neq
i_{2}\,(B_{i_{1}\cdot},B_{i_{2}\cdot})\in\mathsf{Brac}_{d,s}\text{ and
}\forall j_{1}\neq j_{2}\,(B_{\cdot j_{1}},B_{\cdot
j_{2}})\in\mathsf{Brac}_{d,s}\big{\\}}.$ (7)

The following result is a generalization of Proposition 3.2:

######  Proposition 3.11.

For all dimension $d\geq 2$ and all order $s\geq 1$ we have

$\mathsf{U}_{d,s}\subseteq\mathsf{L}_{d,s}\subseteq\mathsf{B}_{d}.$

In particular, the closure of the set of all generalized bracelet matrices is
the full Birkhoff polytope:

$\overline{\bigcup_{s\geq 1}\mathsf{L}_{d,s}}=\mathsf{B}_{d}.$

###### Proof.

The only fact that needs checking is
$\mathsf{U}_{d,s}\subseteq\mathsf{L}_{d,s}$, as then the claim about
$\bigcup_{s\geq 1}\mathsf{L}_{d,s}$ will follow from Theorem 2.8. Fix
$B\in\mathsf{U}_{d,s}$ and let $U\in\mathcal{U}(ds)$ be the corresponding
unitary matrix, i.e., $\frac{1}{s}\operatorname{Tr}(U_{ij}U_{ij}^{*})=B_{ij}$
for all $i,j\in[d]$. Using the unitarity of $U$, one easily verifies that the
pair of different rows $(B_{i_{1}\cdot},B_{i_{2}\cdot})$ satisfies the
generalized bracelet condition $\mathsf{Brac}_{d,s}$ with matrices
$U_{i_{1}1},\ldots,U_{i_{1}d};U_{i_{2}1},\ldots,U_{i_{2}d}\in\mathcal{M}_{s}(\mathbb{C})$.
Analogously for columns; hence $B\in\mathsf{L}_{d,s}$, as claimed.

∎

We show next that the (non-convex) sets of generalized bracelet matrices
$\mathsf{L}_{d,s}$ have a very intricate inclusion structure. To this end, we
study the intersections of $\mathsf{L}_{d,s}$ and of $\mathsf{U}_{d,s}$ with
segments connecting two extremal points of the Birkhoff polytope.

######  Proposition 3.12.

Let $d\geq 3$ and consider permutations $\pi,\sigma\in\mathfrak{S}_{d}$ such
that $\pi^{-1}\sigma$ has a $p$-cycle for some $p\geq 3$. Then the
intersection of the segment $[\pi,\sigma]$ with $\mathsf{U}_{d,s}$ is equal to
the intersection of $[\pi,\sigma]$ with $\mathsf{L}_{d,s}$, and coincides with
the discrete set of convex mixtures of $\pi$ and $\sigma$ with rational
weights with denominator $s$:

$\\{\lambda\in[0,1]\,:\,(1-\lambda)\pi+\lambda\sigma\in\mathsf{U}_{d,s}\\}=\\{\lambda\in[0,1]\,:\,(1-\lambda)\pi+\lambda\sigma\in\mathsf{L}_{d,s}\\}=\\{{k}/{s}\colon
k=0,\ldots,s\\}.$

In particular, this holds for the $1$-faces (i.e. edges) of the Birkhoff
polytope $\mathsf{B}_{d}$ (for which $\pi^{-1}\sigma$ is a $d$-cycle).

###### Proof.

The inclusion of the first set in the second one is trivial. We shall prove
that the second set is contained in the third, and then that the third set is
contained in the first.

For the first inclusion, without loss of generality we may assume that the
decomposition of $\pi^{-1}\sigma$ contains the cycle $(1,2,\ldots,p)$ for some
$p\geq 3$. Furthermore, we may also assume that $\pi$ is the identity
permutation; hence, $(1-\lambda)\pi+\lambda\sigma$ has the form:

$\begin{bmatrix}1-\lambda&\lambda&0&\cdots&0&0\\\
0&1-\lambda&\lambda&\cdots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\cdots&1-\lambda&\lambda\\\
\lambda&0&0&\cdots&0&1-\lambda\end{bmatrix}_{p\times p}\bigoplus\
\Big{[}\text{another matrix of size }(d-p)\times(d-p)\Big{]}.$

Considering any two rows of the $p\times p$ submatrix, we obtain, after
permuting the columns, the slice analysed in Prop. 3.9: if
$(1-\lambda)\pi+\lambda\sigma\in\mathsf{L}_{d,s}$, there exist
$A_{1},A_{2},B_{1},B_{3}\in\mathcal{M}_{s}(\mathbb{C})$ such that

$\begin{bmatrix}A_{1}&A_{2}&0&0&\cdots\\\
B_{1}&0&B_{3}&0&\cdots\end{bmatrix}_{2s\times
ps}\xmapsto{\frac{1}{s}\|\cdot\|^{2}_{F}}\begin{bmatrix}\lambda&1-\lambda&0&0&\cdots\\\
1-\lambda&0&\lambda&0&\cdots\end{bmatrix}_{2\times p}.$

As in the “$\subseteq$” part of the proof of Prop. 3.9, we obtain
$\operatorname{rk}(A_{1})\geq\lceil\lambda s\rceil$ and
$\operatorname{rk}(B_{1})\geq\lceil(1-\lambda)s\rceil$, and then
$\lceil\lambda s\rceil+\lceil(1-\lambda)s\rceil\leq s$; thus,
$\lambda={k}/{s}$ for some $k\in\\{0,1,\ldots,s\\}$, proving the first claim.

As for the second inclusion, use Corollary 2.7, together with the fact that
permutation matrices are unistochastic, to conclude that

$\frac{s-k}{s}\pi+\frac{k}{s}\sigma\in\mathsf{U}_{d,s}.$

Finally, the claim about the face structure of the Birkhoff polytope can be
found in, e.g., [BG77, Theorem 2.2]. ∎

######  Corollary 3.13.

The family of sets $\mathsf{U}_{d,s}$ is not increasing in $s$.

As a final remark, note that the result above was the motivation behind the
study of the slice of the bracelet condition set considered in Eq. (3.1).
Indeed, if we consider the segment $[(1)(2)(3),(123)]$ connecting the identity
and the full cycle permutations of $\mathfrak{S}_{3}$, we have

$(1-\lambda)\cdot(1)(2)(3)+\lambda\cdot(123)=\begin{bmatrix}1-\lambda&\lambda&0\\\
0&1-\lambda&\lambda\\\ \lambda&0&1-\lambda\end{bmatrix}.$

We see that any two rows (or columns) of the matrix above have the zero
pattern found in Eq. (3.1).

##  4\. Generalized orthostochastic matrices

Much of the theory developed in the previous sections for generalized
unistochastic matrices can be carried out to the case of _generalized
orthostochastic_ matrices, which is what we do in this section. However, since
many things are very similar, we shall only present the main definitions and
some observations. We leave the detailed study of generalized orthostochastic
matrices (and that of their quaternionic counterpart, the _qustochastic_
matrices) to future work.

Recall that the function $\varphi_{d,s}$ from Definition 2.1 maps a $d\times
d$ block-matrix (with blocks of size $s\times s$) to the matrix of normalized
squares of Frobenius norms of the blocks. We denote by $\mathcal{O}(n)$ the
group of $n\times n$ orthogonal matrices.

######  Definition 4.1.

We define

$\mathsf{O}_{d,s}:=\varphi_{d,s}(\mathcal{O}(ds))$

to be the set of _generalized orthostochastic matrices_.

As in the complex case, for $s=1$ we recover the usual orthostochastic
matrices, which have received, along with the unistochastic matrices, a lot of
attention in the literature [AYP79, AYC91, BEK+05, ŻKSS03, CD08]. Clearly,
$\mathsf{O}_{d,s}\subseteq\mathsf{U}_{d,s}$, with the inclusion being strict
for $d\geq 3$.

Importantly, the van der Waerden matrix $J_{d}/d$ is orthostochastic if and
only if there exists a real Hadamard matrix of order $d$ [Had93, HW78, KTR05].
This can only happen if $d=2$ or if $d$ is a multiple of four, and it has long
been conjectured that these conditions are also sufficient. In particular, the
distance between $J_{3}/3$ and the set $\mathsf{O}_{d,1}$ is equal to
$\sqrt{2}/3$ [CD08, Proposition 3.2].

We prove now the main result of this section.

######  Proposition 4.2.

For all dimensions $d\geq 2$ and all orders $s\geq 1$, we have

$\mathsf{U}_{d,s}\subseteq\mathsf{O}_{d,2s}.$

###### Proof.

The result follows from the standard embedding
$\mathcal{U}(n)\subseteq\mathcal{O}(2n)$, obtained by replacing a complex
entry $z_{ij}$ by the $2\times 2$ block
$\big{[}\begin{smallmatrix}\operatorname{Re}z_{ij}&\operatorname{Im}z_{ij}\\\
-\operatorname{Im}z_{ij}&\operatorname{Re}z_{ij}\end{smallmatrix}\big{]}$. ∎

######  Corollary 4.3.

For all $d\notin\\{2\\}\cup 4\mathbb{N}$, we have

$J_{d}/d\in\mathsf{O}_{d,2}\setminus\mathsf{O}_{d,1}.$

###### Proof.

We have $J_{d}/d\in\mathsf{U}_{d,1}\subseteq\mathsf{O}_{d,2}$. On the other
hand, since there cannot exist a real Hadamard matrix of order $d$, we have
$J_{d}/d\notin\mathsf{O}_{d,1}$, as claimed. ∎

##  5\. Random generalized unistochastic matrices

In the previous sections we have introduced and discussed generalized
unistochastic matrices, which form the set

$\mathsf{U}_{d,s}=\varphi_{d,s}(\mathcal{U}(ds)),$

where $\varphi_{d,s}$ is the map from Eq. (2.1). As the unitary group
$\mathcal{U}(ds)$ comes equipped with the (normalized) Haar measure
$\mathfrak{h}_{ds}$, it is natural to introduce and examine its image measure
via $\varphi_{d,s}$.

######  Definition 5.1.

We endow the set of generalized unistochastic matrices $\mathsf{U}_{d,s}$ with
the probability measure

$\mu_{d,s}=(\varphi_{d,s})_{\\#}\mathfrak{h}_{ds},$

that is, the image measure of the normalized Haar distribution
$\mathfrak{h}_{ds}$ on $\mathcal{U}(ds)$ through the map $\varphi_{d,s}$. In
other words, if $U\in\mathcal{U}(ds)$ is Haar-distributed, then
$B:=\varphi_{d,s}(U)\in\mathsf{U}_{d,s}$ is $\mu_{d,s}$-distributed.

We recall the following result about the first few joint moments of the
entries of a Haar-distributed random unitary matrix.

######  Lemma 5.2 ([HP00, Proposition 4.2.3]).

Let $U=(U_{ij})_{i,j\in[n]}\in\mathcal{U}(n)$ be Haar-distributed. We have

$\displaystyle\mathbb{E}\Big{[}|U_{ij}|^{2}\Big{]}$
$\displaystyle=\frac{1}{n}\qquad(1\leq i,j\leq n)$
$\displaystyle\mathbb{E}\Big{[}|U_{ij}|^{4}\Big{]}$
$\displaystyle=\frac{2}{n(n+1)}\qquad(1\leq i,j\leq n)$
$\displaystyle\mathbb{E}\Big{[}|U_{ij}|^{2}|U_{i^{\prime}j}|^{2}\Big{]}=\mathbb{E}\Big{[}|U_{ij}|^{2}|U_{ij^{\prime}}|^{2}\Big{]}$
$\displaystyle=\frac{1}{n(n+1)}\qquad(i\neq i^{\prime},j\neq j^{\prime})$
$\displaystyle\mathbb{E}\Big{[}|U_{ij}|^{2}|U_{i^{\prime}j^{\prime}}|^{2}\Big{]}$
$\displaystyle=\frac{1}{n^{2}-1}\qquad(i\neq i^{\prime},j\neq j^{\prime}).$

We leverage now this result to obtain the first moments of a random
generalized unistochastic matrix.

######  Proposition 5.3.

For a Haar-distributed random unitary matrix $U\in\mathcal{U}(ds)$, consider
the corresponding $\mu_{d,s}$-distributed bistochastic matrix

$B:=\varphi_{d,s}(U)=\left(\frac{1}{s}||U_{ij}||_{F}^{2}\right)_{i,j\in[d]}\in\mathsf{B}_{d}.$

For all $1\leq i\neq i^{\prime},j\neq j^{\prime}\leq d$ and $n:=ds$ we have:

$\displaystyle\mathbb{E}\Big{[}B_{ij}\Big{]}$ $\displaystyle=\frac{1}{d}$
$\displaystyle\mathbb{E}\Big{[}B_{ij}^{2}\Big{]}$
$\displaystyle=\frac{d(s^{2}+1)-2}{d(n^{2}-1)}$
$\displaystyle\mathbb{E}\Big{[}B_{ij}B_{i^{\prime}j}\Big{]}=\mathbb{E}\Big{[}B_{ij}B_{ij^{\prime}}\Big{]}$
$\displaystyle=\frac{ds^{2}-1}{d(n^{2}-1)}$
$\displaystyle\mathbb{E}\Big{[}B_{ij}B_{i^{\prime}j^{\prime}}\Big{]}$
$\displaystyle=\frac{s^{2}}{n^{2}-1}.$

###### Proof.

We show the different claims one by one, using Lemma 5.2.

$\mathbb{E}\Big{[}B_{ij}\Big{]}=\frac{1}{s}\mathbb{E}\Big{[}\Tr(U_{ij}U_{ij}^{*})\Big{]}=\frac{1}{s}\mathbb{E}\Big{[}||U_{ij}||_{F}^{2}\Big{]}=\frac{1}{s}\sum_{k,l=1}^{s}\mathbb{E}\Big{[}|U_{ij}(k,l)|^{2}\Big{]}=\frac{1}{s}\cdot
s^{2}\cdot\frac{1}{n}=\frac{1}{d}.$

$\displaystyle\mathbb{E}\Big{[}B_{ij}^{2}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\mathbb{E}\Big{[}\Tr(U_{ij}U_{ij}^{*})^{2}\Big{]}=\frac{1}{s^{2}}\mathbb{E}\Big{[}\big{(}\sum_{k,l=1}^{s}|U_{ij}(k,l)|^{2}\big{)}\big{(}\sum_{p,q=1}^{s}|U_{ij}(p,q)|^{2}\big{)}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\mathbb{E}\Big{[}\big{(}\sum_{k,l=1}^{s}|U_{ij}(k,l)|^{2}\big{)}\big{(}|U_{ij}(k,l)|^{2}+\sum_{\begin{subarray}{c}q=1\\\
q\neq l\end{subarray}}^{s}|U_{ij}(k,q)|^{2}+\sum_{\begin{subarray}{c}p=1\\\
p\neq k\end{subarray}}^{s}|U_{ij}(p,l)|^{2}+\sum_{\begin{subarray}{c}p,q=1\\\
p\neq k\\\ q\neq l\end{subarray}}^{s}|U_{ij}(p,q)|^{2}\big{)}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\cdot
s^{2}\cdot\Big{[}\frac{2}{n(n+1)}+2(s-1)\frac{1}{n(n+1)}+(s-1)^{2}\frac{1}{n^{2}-1}\Big{]}=\frac{d(s^{2}+1)-2}{d(n^{2}-1)}.$

$\displaystyle\mathbb{E}\Big{[}B_{ij}B_{ij^{\prime}}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\mathbb{E}\Big{[}\big{(}\sum_{k,l=1}^{s}|U_{ij}(k,l)|^{2}\big{)}\big{(}\sum_{p,q=1}^{s}|U_{ij^{\prime}}(p,q)|^{2}\big{)}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\mathbb{E}\Big{[}\big{(}\sum_{k,l=1}^{s}|U_{ij}(k,l)|^{2}\big{)}\big{(}\sum_{q=1}^{s}|U_{ij^{\prime}}(k,q)|^{2}+\sum_{\begin{subarray}{c}p,q=1\\\
p\neq k\end{subarray}}^{s}|U_{ij^{\prime}}(p,q)|^{2}\big{)}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\cdot
s^{2}\cdot\Big{[}s\frac{1}{n(n+1)}+s(s-1)\frac{1}{n^{2}-1}\Big{]}.$

$\displaystyle\mathbb{E}\Big{[}B_{ij}B_{i^{\prime}j^{\prime}}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\mathbb{E}\Big{[}\big{(}\sum_{k,l=1}^{s}|U_{ij}(k,l)|^{2}\big{)}\big{(}\sum_{p,q=1}^{s}|U_{i^{\prime}j^{\prime}}(p,q)|^{2}\big{)}\Big{]}$
$\displaystyle=\frac{1}{s^{2}}\cdot
s^{4}\cdot\frac{1}{n^{2}-1}=\frac{s^{2}}{n^{2}-1}.$

∎

######  Remark 5.4.

We could get the same results using the (graphical) Weingarten calculus [CŚ06,
CN10] used to compute general integrals over the unitary group with respect to
the Haar measures.

Let us now analyze the correlations between different matrix elements of $B$.
Recall that _Pearson’s correlation coefficient_ of a pair of random variables
$(X,Y)$ is defined as

$\rho(X,Y):=\frac{\operatorname{Cov}(X,Y)}{\sqrt{\operatorname{Var}(X)}\cdot\sqrt{\operatorname{Var}(Y)}}.$

######  Corollary 5.5.

For all $d\geq 2$, $s\geq 1$, $n:=ds$, and $1\leq i\neq i^{\prime},j\neq
j^{\prime}\leq d$, we have:

$\displaystyle\operatorname{Var}(B_{ij})$
$\displaystyle=\frac{(d-1)^{2}}{d^{2}(n^{2}-1)}$
$\displaystyle\operatorname{Cov}(B_{ij},B_{i^{\prime}j})=\operatorname{Cov}(B_{ij},B_{ij^{\prime}})$
$\displaystyle=-\frac{d-1}{d^{2}(n^{2}-1)}$
$\displaystyle\rho(B_{ij},B_{i^{\prime}j})=\rho(B_{ij},B_{ij^{\prime}})$
$\displaystyle=-\frac{1}{d-1}$
$\displaystyle\operatorname{Cov}(B_{ij},B_{i^{\prime}j^{\prime}})$
$\displaystyle=\frac{1}{d^{2}(n^{2}-1)}$
$\displaystyle\rho(B_{ij},B_{i^{\prime}j^{\prime}})$
$\displaystyle=\frac{1}{(d-1)^{2}}.$

###### Proof.

The formulas follow by direct computation from Proposition 5.3; we leave the
details to the reader. ∎

######  Remark 5.6.

Note that the covariance (and the correlation coefficient) of elements of $B$
situated on the same row (or column) is negative; this anti-correlation is
explained by the normalization conditions of bistochastic matrices. However,
the correlation of elements not belonging to the same row and column is
positive.

Let us also note that while the covariance of the matrix elements of $B$
decreases (in absolute value) with the parameter $s$, the correlation
coefficient is constant (at fixed matrix dimension $d$).

The fact that the variance of $B_{ij}$ decreases with the parameter $s$ at
fixed matrix dimension $d$ means that the distribution $\mu_{d,s}$
concentrates, as $s\to\infty$, around the van der Waerden matrix. The same
phenomenon can be seen at the level of spectra, see Figure 3: the non-trivial
eigenvalues of the random matrices $B\sim\mu_{d,s}$ tend to concentrate around
the origin as $s$ grows. We point the reader interested in spectral properties
of bistochastic and unistochastic matrices to the papers [ŻKSS03, CSBŻ09].
Finally, note the pair of complex eigenvalues outside the gray-bounded region
in the middle top panel of Figure 3; they are a signature of the fact that
$\mathsf{U}_{3,2}\supsetneq\mathsf{U}_{3,1}$, see Example 2.3.

Figure 3. Spectra of random generalized bistochastic matrices. We plot the
(complex) eigenvalues of 10 000 samples from the measure $\mu_{d,s}$
introduced in Def. 5.1. On the top row, we have $d=3$, and, respectively,
$s=1,2,3$. On the bottom row, $d=4$ and $s=1,2,4$. We also plot in gray the
hypocycloid curves which are conjectured to bound the spectra in the
unistochastic case $s=1$, see [ŻKSS03, Section 4.3].

Acknowledgements. I.N. was supported by the ANR projects ESQuisses, grant
number ANR-20-CE47-0014-01 and STARS, grant number ANR-20-CE40-0008, and by
the PHC program _Star_ (Applications of random matrix theory and abstract
harmonic analysis to quantum information theory). A.S. was supported by the
ANR project Quantum Trajectories, grant number ANR-20-CE40-0024-01. Z.O. would
like to extend my sincere gratitude to Professor Ion Nechita and Anna
Szczepanek for their invaluable guidance and insightful discussions.
Additionally, I also wish to express my appreciation to all professors at IMT
for their exceptional guidance and support, which makes my M1 year in Toulouse
a truly fulfilling and enriching experience.

## References

  * [AYC91] Yik-Hoi Au-Yeung and Che-Man Cheng.  Permutation matrices whose convex combinations are orthostochastic.  Linear Algebra and its Applications, 150:243–253, 1991.
  * [AYP79] Yik-Hoi Au-Yeung and Yiu-Tung Poon.  3$\times$ 3 orthostochastic matrices and the convexity of generalized numerical ranges.  Linear Algebra and its Applications, 27:69–79, 1979.
  * [BEK+05] Ingemar Bengtsson, Åsa Ericsson, Marek Kuś, Wojciech Tadej, and Karol Życzkowski.  Birkhoff’s polytope and unistochastic matrices, n= 3 and n= 4.  Communications in mathematical physics, 259:307–324, 2005.
  * [Ben04] Ingemar Bengtsson.  The importance of being unistochastic.  arXiv:quant-ph/0403088, 2004.
  * [BG77] Richard A Brualdi and Peter M Gibson.  Convex polyhedra of doubly stochastic matrices. i. applications of the permanent function.  Journal of Combinatorial Theory, Series A, 22(2):194–230, 1977\.
  * [CD08] Oleg Chterental and Dragomir Ž Doković.  On orthostochastic, unistochastic and qustochastic matrices.  Linear algebra and its applications, 428(4):1178–1201, 2008.
  * [CN10] Benoît Collins and Ion Nechita.  Random quantum channels i: Graphical calculus and the bell state phenomenon.  Communications in Mathematical Physics, 297(2):345–370, feb 2010\.
  * [CŚ06] Benoît Collins and Piotr Śniady.  Integration with respect to the haar measure on unitary, orthogonal and symplectic group.  Communications in Mathematical Physics, 264(3):773–795, mar 2006\.
  * [CSBŻ09] Valerio Cappellini, Hans-Jürgen Sommers, Wojciech Bruzda, and Karol Życzkowski.  Random bistochastic matrices.  Journal of Physics A: Mathematical and Theoretical, 42(36):365209, 2009.
  * [Gut13] Eugene Gutkin.  On a multi-dimensional generalization of the notions of orthostochastic and unistochastic matrices.  Journal of Geometry and Physics, 74:28–35, 2013.
  * [Had93] Jacques Hadamard.  Résolution d’une question relative aux déterminants.  Bull. sci. math, 17(1):240–246, 1893.
  * [HP00] F. Hiai and D. Petz.  The Semicircle Law, Free Random Variables and Entropy.  Mathematical surveys and monographs. American Mathematical Society, 2000\.
  * [HW78] A Hedayat and Walter Dennis Wallis.  Hadamard matrices and their applications.  The annals of statistics, pages 1184–1238, 1978.
  * [KTR05] Hadi Kharaghani and Behruz Tayfeh-Rezaie.  A hadamard matrix of order 428.  Journal of Combinatorial Designs, 13(6):435–440, 2005.
  * [Nak96] Hiroshi Nakazato.  Set of 3$\times$ 3 orthostochastic matrices.  Nihonkai mathematical journal, 7(2):83–100, 1996.
  * [RMKPŻ22] Grzegorz Rajchel-Mieldzioć, Kamil Korzekwa, Zbigniew Puchała, and Karol Życzkowski.  Algebraic and geometric structures inside the birkhoff polytope.  Journal of Mathematical Physics, 63(1), 2022.
  * [ŻKSS03] Karol Życzkowski, Marek Kuś, Wojciech Słomczyński, and Hans-Jürgen Sommers.  Random unistochastic matrices.  Journal of Physics A: Mathematical and General, 36(12):3425, 2003\.


# Emotion-Conditioned Melody Harmonization with Hierarchical Variational
Autoencoder

Shulei Ji∗ and Xinyu Yang∗ ∗Shulei Ji and Xinyu Yang are with the School of
Computer Science and Technology, Xi’an Jiaotong University, Xi’an, 710049,
China (e-mail<EMAIL_ADDRESS>yxyphd@mail.xjtu.edu.cn)(Corresponding author: Xinyu Yang)

###### Abstract

Existing melody harmonization models have made great progress in improving the
quality of generated harmonies, but most of them ignored the emotions beneath
the music. Meanwhile, the variability of harmonies generated by previous
methods is insufficient. To solve these problems, we propose a novel LSTM-
based Hierarchical Variational Auto-Encoder (LHVAE) to investigate the
influence of emotional conditions on melody harmonization, while improving the
quality of generated harmonies and capturing the abundant variability of chord
progressions. Specifically, LHVAE incorporates latent variables and emotional
conditions at different levels (piece- and bar-level) to model the global and
local music properties. Additionally, we introduce an attention-based melody
context vector at each step to better learn the correspondence between
melodies and harmonies. Objective experimental results show that our proposed
model outperforms other LSTM-based models. Through subjective evaluation, we
conclude that only altering the types of chords hardly changes the overall
emotion of the music. The qualitative analysis demonstrates the ability of our
model to generate variable harmonies.

##  I Introduction

Melody harmonization is a conditional music generation task to generate chord
progression conditioned on the given melody [1, 2]. This task aims at
generating human-like harmonies to enrich monotonous melodies and enhance
musical expression. In addition to improving computer creativity,
understanding human composition mechanisms, and bringing inspiration to
composers, melody harmonization can also be applied in areas such as
entertainment and education.

Many previous neural-network-based melody harmonization methods used recurrent
neural networks [3, 4], e.g., Long Short-Term Memory (LSTM), to generate the
corresponding chords for the melody input, which hardly generated variable
results for the same input and ignore the learning of chord transitions.
Recently, the variational auto-encoder (VAE) [5, 6] has been utilized to
generate coherent and varied harmonies. However, these studies all neglected
the emotion conveyed by the harmonies.

Generating music conditioned on emotions using deep learning techniques, also
known as emotion-conditioned music generation, has gradually become a trend
[7, 8, 9], where researchers tend to generate music from scratch based on the
given emotions. Nevertheless, few studies explore the impact of emotion on
conditional music generations, i.e., generating partial music given other
musical contexts, such as melody harmonization, the focus of this paper.

The role of emotional condition in the melody harmonization task warrants
careful consideration. There are two cases: firstly, the harmony exhibits
distinct emotions, so it maybe feasible to change the emotion of the overall
lead sheet by generating harmony with specified emotion; secondly, in cases
where discerning the emotion of the harmony is challenging, the emotional
conditions can serve to assist in generating better harmonies.

Our previous work focused on the second case, which proposed an emotion-
conditioned hierarchical transformer VAE named EmoMusicTV to explore the
impact of emotional conditions on multiple tasks on lead sheets. It has been
proven that introducing emotions into the model can benefit the generated
results for two conditional generation tasks, i.e., melody harmonization and
melody generation given harmony. Note that the above practice is based on our
empirical findings that melody seemingly dominates in conveying the emotion of
music rather than harmony. However, this assumption has not been
experimentally verified.

In addition, the transformer-based EmoMusicTV was compared with two LSTM-based
methods [11, 5] for the melody harmonization task in our previous work , one
of which was also based on VAE generative model. However, it is difficult for
such a comparison to highlight the advantages of the hierarchical hidden
variable structure, since the transformer appears to be superior in sequence
learning when compared to LSTM. Therefore, it is hard to say whether the
advantages of EmoMusicTV come from the transformer or the hierarchical hidden
variable structure.

To make up for the deficiency in previous work, a novel LSTM-based
Hierarchical VAE (LHVAE) is proposed to generate harmonies with emotional
conditions at different levels and improve the diversity of the generated
harmonies. Following EmoMusicTV , LHVAE also incorporates a hierarchical
latent variable structure, where a piece-level latent variable learns holistic
properties and a chain of bar-level latent variables captures short-term
variations. Emotional conditions at different levels are embedded in their
corresponding latent spaces neatly. The comparison between LHVAE and other
LSTM-based models can better highlight the effectiveness of hierarchical
latent variable structure since all the models adopt the same backbone, i.e.,
LSTM.

(a) Event-based music representation.

(b) The architecture of the LSTM-based hierarchical VAE. The model first
produces a piece-level latent variable, and then iteratively generates bar-
level latent variables at the beginning of each bar.

Figure 1: The adopted music representation and the proposed model.

We conduct a series of experiments to verify the effectiveness of our proposed
model for melody harmonization and to answer the following question of
concern: Whether the harmonies generated under different emotional conditions
for the same melody can change the overall emotion of music. Experimental
results show that our proposed model outperforms other LSTM-based methods via
objective evaluation and brings music variability through qualitative
analysis. The subjective evaluation shows that it is hard to change the
musical emotion perceived by humans via only altering the chord
progressions111Note that in this paper we only change the types of chords,
regardless of the density of chords, i.e., each bar contain two chords..

##  II Related Work

Deep neural networks have been widely applied for automatic melody
harmonization. Lim et al. [3] utilized a two-layer bidirectional LSTM (BiLSTM)
network to learn the correspondence between melody and chord sequences. Yeh et
al. [4] proposed an extended BiLSTM model to predict not only chord labels but
also chord functions. Sun et al. [11] input the melody and the partially
masked chord sequence into the BiLSTM to learn the missing parts in the chord
sequence. Then they adopted blocked Gibbs sampling for inference. Chen et al.
[5] proposed a VAE-based melody harmonization model conditioned on surprise
contours. Rhyu et al. [6] proposed the regularized variational Transformer to
control the attribute of generated chords. Zeng et al. [12] proposed a novel
reinforcement learning method to first learn a refined structured melody
representation, and then to generate chords for each learned segment
representation. Our previous work [13] also proposed a novel melody
harmonization system based on reinforcement learning, which contains three
well-designed reward modules for introducing other guidance, such as music
theory knowledge. However, none of the above methods explore melody
harmonization conditioned on emotions.

To the best of our knowledge, only Takahashi et al. [14] and EmoMusicTV
generated harmonic arrangements for the given melodies with emotional
constraints. Takahashi et al. [14] adopted crowd-sourced mood tags as emotion
labels, which are considered as “weak labels”. They found that the errors
between the perceived and the target emotions remain relatively large and it
is difficult to prove that the model consistently expresses the desired
emotions. EmoMusicTV empirically assumed that melody dominates in conveying
the emotion of music and took emotional conditions as auxiliary information to
improve the quality of the generated harmonies. However, neither of the above
researches addressed our proposed question, i.e., the possibility of changing
the overall emotion of music via altering chord progressions1.

##  III Methodology

###  III-A Music Representation

Same as , we adopt the event-based music representation [1] to represent the
lead sheets, as illustrated in Fig. 1(a). Each melody note is represented
using a pair of consecutive tokens, namely the pitch and duration tokens. The
beginning of a bar is indicated by the bar token. Each chord is represented by
a super token [15] consisting of chord type and root tone. The numbers of
pitch (plus REST) and duration events are 61 and 37, respectively. There are 7
distinct chord type events (including REST) and 41 different root tone events.
Consequently, the melody event is encoded as a 99-dimensional one-hot vector,
while the chord event is represented by a multi-hot vector of 48 dimensions.

###  III-B LSTM-based Hierarchical VAE

We propose a novel LSTM-based Hierarchical VAE (LHVAE) for emotion-conditioned
melody harmonization to improve the variability of the generated chord
progressions and to explore the impact of emotional conditions on the
generated harmonies, as shown in Fig. 1(b). The hierarchical latent variable
structure in LHVAE follows the design of EmoMusicTV .

We represent the melody sequence as $Y=\\{y_{1},y_{2},...,y_{m}\\}$, the chord
sequence as $X=\\{x_{1},x_{2},...,x_{n}\\}$, where $m$ and $n$ are the number
of melody events and chords. The proposed model contains latent variables at
different levels, i.e., piece-level latent variable $z_{P}$ and bar-level
latent variables $Z_{B}=\\{z_{B}^{1},\ldots,z_{B}^{I}\\}$, $I$ is the number
of bars. The piece-level emotion $s_{P}$ and emotion expressed in each bar
$S_{B}=\\{s_{B}^{1},\ldots,s_{B}^{I}\\}$ are embedded in their corresponding
latent spaces to guide the chord generation. We use $S$ to represent all the
emotions in the follow-up.

The latent variables exhibit structured dependencies that each bar-level
latent variable $z_{B}^{i}$ is conditioned not only on $z_{P}$ but also on the
previous latent trajectories $Z_{B}^{\textless i}$. By variational inference,
we derive the objective function which is the Evidence Lower Bound ($ELBO$) of
the conditional log likelihood of $p(X|Y,S)$

$\begin{split}\log
p(X|Y,S)&\geq\textstyle\sum_{i}\mathbb{E}_{q_{\phi}(Z_{B}^{\textless
i}|X,z_{P},s_{P},S_{B},Y)}[{ELBO_{i}}]-\\\ &\quad\
KL[q_{\phi}(z_{P}|s_{P},X,Y)\|p_{\varphi}(z_{P}|s_{P},Y)]\\\
&=-\mathcal{L}(X;Y,S,\theta,\phi,\varphi)\end{split}$ (1)

where KL is the Kullback-Leibler divergence. The loss
$\mathcal{L}(X;Y,S,\theta,\phi,\varphi)$ consists of a piece-level KL term and
the sum of negative bar-level $ELBO_{i}$. The $ELBO_{i}$ can be factorized as:

$\begin{split}ELBO_{i}=&\mathbb{E}_{q_{\phi}(z_{B}^{i}|X,U)}[\log
P_{\theta}(X^{i}|X^{\textless i},z_{B}^{i},U)]-KL_{i}\\\
KL_{i}=&KL[q_{\phi}(z_{B}^{i}|X,U)\|p_{\varphi}(z_{B}^{i}|X^{\textless
i},U)]\end{split}$

(2)

where $U$ refers to the union of $Y,z_{P},s_{P},Z_{B}^{\textless i},$ and
$S_{B}^{\leq i}$, $X^{i}$ and $X^{\textless i}$ represent the chords within
$i$th bar and the chords before $i$th bar, respectively. $\log P_{\theta}$ is
the log-likelihood of the predicted probabilities, and $p_{\varphi}$ and
$q_{\phi}$ are the prior distribution and the approximate posterior
distribution, respectively.

We assume that the latent variables at different levels obey the isotropic
Gaussian distribution, like previous work [16] did. The mean $\mu$ and
standard deviation $\sigma$ are generated by the Multi-Layer Perceptron
(MLP)-based prior model and recognition model [17], as depicted in Fig. 1(b).
The piece-level $z_{P}$ is computed only once at the beginning of generating
the entire harmony, which initializes the first hidden state $h_{0}$ of the
LSTM. While the bar-level latent variables $Z_{B}$ are computed at the
beginning of each bar. Note that each bar contains two chords in this paper.

Referring to Eq. (1), the piece-level prior model $p_{\varphi}$ takes as input
the melody and the piece-level emotion to obtain $z_{P}$, and the piece-level
recognition model $q_{\phi}$ additionally needs the chord. According to the
Eq. (2), the bar-level prior model $p_{\varphi}$ takes historical output
$X^{\textless i}$ and other elements as input to get $z_{B}^{i}$, while the
bar-level recognition model $q_{\phi}$ requires the entire output sequence
$X$. Thus, we add a reverse chord encoder to model output sequences
$X^{\geqslant i}$ that have not been generated yet. For code implementation,
the last hidden state $h_{2i}$ of LSTM, implicitly containing $X^{\textless
i},Y,z_{P},s_{P},Z_{B}^{\textless i},$ and $S_{B}^{\textless i}$, are input
into the bar-level prior and recognition network, as illustrated in Fig. 1(b).

Figure 2: Three ways to inform the model bar-level hidden variable and melody
context vector. $z_{B}^{i}$ and $s_{B}^{i}$ are latent variable and emotional
label of the $i$th bar, respectively. $c_{2i}$ is the melody context vector at
time step $2i$.

Additionally, to better learn the correspondence between melody and chord, we
calculate the melody context vector $c_{t}$ at each step by making the hidden
state $h$ of LSTM account for the melody hidden states $h^{Y}$ within the
current bar.

$c_{t}=\sum_{j=m^{i}_{start}}^{m^{i}_{end}}\alpha_{tj}h_{j}^{Y},\
\alpha_{tj}=f_{attention}(h_{t-1},h_{j}^{Y})$ (3)

where $h_{j}^{Y}$ is the $j$th hidden states of melody encoder,
$m^{i}_{start}$ and $m^{i}_{end}$ indicate start and end of melody sequences
within the $i$th bar. $f_{attention}$ is a function that performs $softmax$
after the inner product, and $\alpha_{tj}$ is the attention weight.

At each step, we feed the bar-level latent variable, the context vector and
the emotional label together into the model to generate chords. As shown in
Fig. 2, we provide three ways to inform the model bar-level latent variable
$z_{B}$ and melody context vector $c$. The experimental results will prove
that the second one is the best choice, so the model diagram in Fig. 1(b) is
drawn according to the second type.

TABLE I: Holistic emotion distribution of NMD dataset.

Negative | Neutral | Positive
---|---|---
624 (5.52%) | 1272 (11.25%) | 9408 (83.23%)

TABLE II: The results of three ways informing model latent variables and
context vectors. The metric values closest to GT are marked in bold. We adopt
two-tailed t-test scores to determine whether there is a significant
difference between the objective results of the model and GT. The objective
results of the models are expected to be not significantly different from
those of GT, i.e., p $\geqslant 0.05$. All p-values equal to or greater than
0.05 are marked in bold.

Models | CHE | p-value | CC | p-value | CTD | p-value | CTnCTR | p-value | PCS | p-value | MCTD | p-value
---|---|---|---|---|---|---|---|---|---|---|---|---
GT | 0.7493$\pm$0.20 | - | 6.3120$\pm$2.04 | - | 0.1034$\pm$0.07 | - | 0.4214$\pm$0.11 | - | 0.1801$\pm$0.40 | - | 0.0999$\pm$0.03 | -
LHVAE1 | 0.7083$\pm$0.20 | 1.2e-05 | 7.9930$\pm$2.78 | 2.9e-31 | 0.0257$\pm$0.04 | 6.7e-56 | 0.4214$\pm$0.11 | 0.1670 | 0.1687$\pm$0.37 | 0.0032 | 0.1007$\pm$0.03 | 1.2e-36
LHVAE2 | 0.7993$\pm$0.22 | 4.5e-21 | 8.8050$\pm$3.09 | 1.1e-47 | 0.0912$\pm$0.05 | 0.3107 | 0.4212$\pm$0.11 | 0.2322 | 0.1872$\pm$0.41 | 0.0530 | 0.0995$\pm$0.03 | 1.2e-14
LHVAE3 | 0.8122$\pm$0.23 | 3.1e-25 | 8.9900$\pm$2.68 | 2.4e-61 | 0.0788$\pm$0.04 | 2.5e-07 | 0.4214$\pm$0.11 | 0.6367 | 0.1905$\pm$0.42 | 0.0415 | 0.0992$\pm$0.03 | 6.2e-19

TABLE III: Model comparison results.

Models | CHE | p-value | CC | p-value | CTD | p-value | CTnCTR | p-value | PCS | p-value | MCTD | p-value
---|---|---|---|---|---|---|---|---|---|---|---|---
GT | 0.7493$\pm$0.20 | - | 6.3120$\pm$2.04 | - | 0.1034$\pm$0.07 | - | 0.4214$\pm$0.11 | - | 0.1801$\pm$0.40 | - | 0.0999$\pm$0.03 | -
BGS[11]* | 0.6994$\pm$0.17 | 3.2e-10 | 5.5820$\pm$1.77 | 0.0746 | 0.1407$\pm$0.06 | 3.8e-28 | 0.4216$\pm$0.11 | 0.0879 | 0.1875$\pm$0.40 | 0.0038 | 0.0994$\pm$0.03 | 2.1e-27
SurpriseNet[5]* | 0.7964$\pm$0.21 | 1.8e-10 | 7.9980$\pm$1.79 | 2.6e-06 | 0.1193$\pm$0.05 | 5.7e-07 | 0.4212$\pm$0.11 | 0.1717 | 0.1834$\pm$0.41 | 0.1238 | 0.0997$\pm$0.03 | 8.3e-05
LHVAE2 | 0.7993$\pm$0.22 | 1.0e-12 | 8.8050$\pm$3.09 | 1.1e-47 | 0.0912$\pm$0.05 | 0.3107 | 0.4212$\pm$0.11 | 0.2322 | 0.1872$\pm$0.41 | 0.0530 | 0.0995$\pm$0.03 | 1.2e-14

  * *

The results of BGS[11] and SurpriseNet[5] are derived from [].

###  III-C Loss function

The output of LSTM at each step will be fed into two feed-forward heads to
predict chord type and root tone, respectively. To alleviate the unbalanced
emotion distribution caused by the uneven number of chord types (see Table I),
we use the chord balancing trick proposed by Sun et al. [11]. Specifically, we
assign a weight $w_{c}$ for each chord type $c$ in the negative log-likelihood
(NLL) loss of chord type.

$\begin{split}L\left(X;Y,S,\theta\right)&=-\sum_{c}{w_{c}x_{t,c}\log
P_{\theta}\left(x_{t,c}|x_{<t},Y,S\right)}\\\
w_{c}&=\left|C\right|\times\frac{{1}/{\left(1e5+n_{c}\right)}}{\sum_{c^{\prime}}{1}/{\left(1e5+n_{c^{\prime}}\right)}}\end{split}$
(4)

where $x_{t,c}$ is the $c$-th element of the one-hot vector $x_{t}$ of chord
type and $\theta$ represents the model parameter. $\left|C\right|$ is the
number of chord types (i.e., 7), $n_{c}$ is the number of chord type $c$ in
the training data and $1e5$ is a temperature-like parameter for count
smoothing. Note that the reconstruction loss for the root tone is the
unweighted NLL loss.

##  IV Experiments

###  IV-A Experimental Setup

####  IV-A1 Dataset

We adopt the preprocessed Nottingham Music Dataset222abc.sourceforge.net/NMD/
(NMD) by with the annotated piece- and bar-level emotional labels. One
difference is that we simplify five emotion categories (very negative,
moderate negative, neutral, moderate positive, and very positive) in into
three categories: negative (0), neutral (1), and positive (2). The holistic
emotion distribution of the dataset is presented in Table I, from which we see
that the number of musical pieces with positive emotion is much larger than
that of neutral and negative emotions, and the number of negative pieces is
the least. Since the emotional labels are computed from chords using the
method of Makris et al.[18], the unbalanced distribution of emotions is caused
by the uneven number of different chord types. Thus, we used chord balancing
[9] to alleviate this problem in this paper (Section III-C). Note that each
bar contains two chords.

####  IV-A2 Implementation Details

The encoders and decoder in the LHVAE model are all 2-layer LSTMs. The size of
the hidden state, latent space and batch are 256, 128 and 24, respectively.
The learning rate is 3e-4 and the model is trained with the Adam optimizer. We
use the KL annealing trick [19, 20] to alleviate the posterior collapse
problem [21] in VAE training. Specifically, the piece- and bar-level KL
weights start from zero and linearly increase 1e-5 and 1e-4 per epoch,
respectively. The dataset is randomly divided into 90% training set and 10%
validation set for objective evaluation. The emotional labels in the
evaluation stage are derived from the validation set.

####  IV-A3 Compared Models

We compare our model with Ground-Truth (GT) and other two LSTM-based melody
harmonization models.

BGS[11]: A BiLSTM model with orderless NADE training process and chord
balancing, then using blocked Gibbs sampling (BGS) for inference.

SurpriseNet[5]: A BiLSTM-based CVAE taking surprise contours as conditions to
generate chord sequences.

####  IV-A4 Evaluation Metrics

To evaluate the generated chords quantitatively, we use six music metrics
defined by Yeh et al. [4], i.e., Chord histogram entropy (CHE), Chord coverage
(CC), Chord tonal distance (CTD), Chord tone to non-chord tone ratio (CTnCTR),
Pitch consonance score (PCS) and Melody-chord tonal distance (MCTD). Among
them, CC, CHE and CTD evaluate chords in terms of the number of chord types,
chord distribution and chord transition, while the rest metrics evaluate the
consistency and harmonicity between melody and harmony. Note that the
generated results are generally better when the metrics values are closer to
the GT. Besides, we use four criteria [11] to evaluate the generated chords
subjectively, i.e., Chord progression, Coherence, Interestingness, and
Overall.

TABLE IV: Emotional control accuracy.

Accuracy | LHVAE1 | LHVAE2 | LHVAE3
---|---|---|---
Piece-level | 0.98 | 0.95 | 0.89
Bar-level | 0.96 | 0.93 | 0.88

TABLE V: Ablation study results of emotional conditions and latent variables.
$S$ denotes the emotional conditions. $z_{P}$ is the piece-level latent
variable, $Z_{B}$ is the chain of bar-level latent variables, and $Z_{A}$
represents all the latent variables.

Models | CHE | p-value | CC | p-value | CTD | p-value | CTnCTR | p-value | PCS | p-value | MCTD | p-value
---|---|---|---|---|---|---|---|---|---|---|---|---
GT | 0.7493$\pm$0.20 | - | 6.3120$\pm$2.04 | - | 0.1034$\pm$0.07 | - | 0.4214$\pm$0.11 | - | 0.1801$\pm$0.40 | - | 0.0999$\pm$0.03 | -
LHVAE2 | 0.7993$\pm$0.22 | 1.0e-12 | 8.8050$\pm$3.09 | 1.1e-47 | 0.0912$\pm$0.05 | 0.3107 | 0.4212$\pm$0.11 | 0.2322 | 0.1872$\pm$0.41 | 0.0530 | 0.0995$\pm$0.03 | 1.2e-14
w/o $S$ | 0.7395$\pm$0.21 | 0.0637 | 7.8660$\pm$2.99 | 1.6e-18 | 0.0522$\pm$0.06 | 5.9e-21 | 0.4211$\pm$0.11 | 0.0392 | 0.1704$\pm$0.38 | 0.0089 | 0.1006$\pm$0.03 | 2.8e-24
w/o $z_{P}$ | 0.7746$\pm$0.21 | 1.9e-12 | 8.6660$\pm$3.13 | 1.5e-44 | 0.0836$\pm$0.05 | 0.0343 | 0.4213$\pm$0.11 | 0.4765 | 0.1880$\pm$0.40 | 0.0047 | 0.0999$\pm$0.03 | 0.0657
w/o $Z_{B}$ | 0.7562$\pm$0.20 | 1.8e-10 | 7.0400$\pm$2.13 | 2.6e-27 | 0.0782$\pm$0.05 | 0.0001 | 0.4215$\pm$0.11 | 0.0326 | 0.1774$\pm$0.40 | 0.6148 | 0.1000$\pm$0.03 | 0.9385
w/o $Z_{A}$ | 0.7445$\pm$0.20 | 0.0008 | 6.6810$\pm$2.08 | 1.9e-15 | 0.0586$\pm$0.05 | 1.4e-24 | 0.4214$\pm$0.11 | 0.0201 | 0.1756$\pm$0.39 | 0.9358 | 0.1001$\pm$0.03 | 3.8e-05

###  IV-B Experimental Results and Analyses

Figure 3: Histogram of human rating results. M1, M2, M3 and M4 represent chord
progression, coherence, interestingness and overall, respectively.
($s_{GT}$)$\rightarrow$($s_{Given}$) in the legend denotes music with ground-
truth annotated emotion $s_{GT}$ is re-harmonized with given emotional
conditions $s_{Given}$. The shadow coverage indicates significant differences
from human-composed chords. Figure 4: Emotion recognition results given by
humans. The x-axis represents the given emotional labels, and the y-axis
denotes the frequency of emotions perceived by the subjects.

####  IV-B1 Objective Evaluation

##### Comparative study

We first investigate three ways (shown in Fig. 2) of feeding bar-level latent
variable and melody context vector into the model, and the results are
demonstrated in Table II. We conclude that LHVAE2 gets the best results,
especially on the metric evaluating chord transition (CTD). LHVAE1 perform
better on CHE and CC, indicating better learning of chord distribution and the
number of chord type. It seems that the concatenation of input and latent
variables leads to better results of the melody-chord harmonicity (PCS, MCTD),
while concatenating input with context vector leads to better CHE and CC.

We also explored the emotional control capabilities of these three ways.
Specifically, we generate harmonies for the melodies in the validation set
given the annotated emotional labels. Then we calculate the piece- and bar-
level emotions of the generated harmonies using the method proposed in [18].
The emotional control accuracy is the number of generated harmonies conveying
given emotions divided by the total number of pieces or bars, as shown in
Table IV. LHVAE1 achieves the highest accuracy in both piece- and bar-level
emotion control followed by LHVAE2.

After making a trade-off between the objective results and the emotional
control accuracy, we select LHVAE2 as the final version of the model. Then we
compare it with other two LSTM-based models, as shown in Table III. We can see
that LHVAE2 excels at learning the chord transitions, showing the best CTD.
And LHVAE2 is on par with SurpriseNet on learning the harmonicity between
melody and chord (CTnCTR, PCS, MCTD). However, LHVAE2 and SurpriseNet perform
poorly on CHE and CC, presenting higher values than BGS, which has no latent
variable. We conjecture that more latent variables enable the generated chords
more diverse than GT.

It is worth noting that although LHVAE2 outperforms the other two LSTM models,
it is inferior to the EmoMusicTV (compared with the 5th row in Table IV(a) in
). This result shows that the advantages of EmoMusicTV come from the joint
effect of the transformer and hierarchical hidden variable structure, and
combining a stronger backbone with this structure can achieve better results.

##### Ablation study

We first conduct an ablation study on LHVAE2 by eliminating the emotions, the
results are shown in Table V. We observe that removing emotions (w/o $S$)
makes the results on CHE and CC improved, but the results of other metrics
become worse. The conclusion comes that emotional conditions assist in
improving the quality of generated chords in terms of chord transition and
melody-chord harmonicity. This conclusion is slightly different from the
results of EmoMusicTV w/o $S$ on the TheoryTab
Dataset333https://www.hooktheory.com/theorytab (TTD), which perform worse on
almost all metrics (refer to 4th row in Table V(a) in ). However, both
experimental results demonstrate that emotional conditions can improve the
generated harmonies more or less.

Additionally, we also show the effects of hierarchical latent variable
structure in Table V. Referring to the p-value, the lack of piece-level latent
variables (w/o $z_{P}$) damages CTD and PCS, and removing bar-level latent
variables (w/o $Z_{B}$ or $Z_{A}$) enables the results of CTD and CTnCTR far
away from the GT. However, eliminating latent variables makes CHE, CC and MCTD
closer to the GT than the complete LHVAE2. By and large, the combination of
latent variables at different levels achieves better performance in learning
chord transitions (CTD) and melody-chord harmonicity (CTnCTR and PCS),
presenting insignificant differences with GT.

####  IV-B2 Subjective Evaluation

To answer the question introduced in Section I (i.e., whether the chord
progressions generated under different emotional conditions for the same
melody can change the overall emotion of music) and evaluate the generated
chords subjectively, we invite human subjects to participate in two listening
tests. In the first test, subjects are asked to identify the emotions for 9
pieces of music, which are generated by first selecting three lead sheets
whose ground-truth holistic emotions are negative, neutral and positive
respectively, then generating three harmonies for each melody given three
different emotional conditions. Specifically, the piece- and bar-level
emotions of the three emotional conditions are set to the same value at a time
(i.e., all 0, all 1 or all 2)444The accuracy of bar-level positive, neutral
and negative emotions for the generated chords are 97.98%, 100% and 88.89%,
respectively, which is computed as the number of generated emotions that are
consistent with given emotions divided by the total number.. For the second
test, subjects are asked to rate 12 chord progressions according to the
subjective criteria, including the above-mentioned nine and their ground-
truth. A total of 20 subjects containing 10 males and 10 females participated
in these two surveys. Their ages range from 18 to 30 years old, and six of
them have learning experiences in musical instruments (e.g., piano and Chinese
zither) and are familiar with harmonic theory more or less.

The human emotion recognition results are shown in Fig. 4. The three sub-
figures correspond to three melodies with negative, neutral and positive
ground-truth emotions. We observe that although the melody is harmonized under
different emotions, the music is still identified by the majority as its
ground-truth emotion. It seems hard to change the overall emotion of music
only by altering the simple pillar chords. On another front, the ground-truth
emotional labels annotated using the method of [18] from the harmony align
with human annotation to some extent.

The human rating results for the generated chord progressions are shown in
Fig. 3. The three sub-figures still correspond to the three melodies with
negative, neutral and positive ground-truth emotions. On the whole, human-
composed harmonies are slightly better than the generated harmonies, and the
quality of chords generated under different emotions has no obvious
regularity. Through the two-tailed t-test, we find that in most cases, there
are no significant differences between the generated and human-composed
harmonies, and cases with significant differences are shaded in Fig. 3.

Figure 5: The visualizations of the generated chord progressions. For each
piece of music, we show the original melody (top), the original harmony
(bottom), and the harmonies generated by LHVAE2 (middle three).

####  IV-B3 Qualitative Analysis

Given the inherent nature of VAE to sample from the latent spaces and
introduce variability, it becomes possible to perform latent variable
samplings multiple times and re-harmonize a given melody. Especially, the
chain of bar-level latent variables in our proposed model makes the generated
results more varied. Fig. 5 shows two pieces with chord progressions created
by humans and the LHVAE2 model. It can be seen that our model generates chord
progressions corresponding well to the ground-truth emotional labels and
generates varied harmonies for the same melody while maintaining similar
harmonic patterns.

##  V Conclusion and Future Work

In this paper, we make up for the deficiency of prior work on melody
harmonization by proposing a novel LSTM-based VAE (LHVAE) model with a
hierarchical latent variable structure and time-varying emotional conditions.
We combine latent variables and emotions at different levels to produce
harmonious, coherent, and diverse chord progressions. Experimental results
show that LHVAE outperforms other LSTM-based melody harmonization models. We
conclude that the guidance of emotional conditions can improve the quality of
generated harmonies in terms of chord transition and harmonicity, while only
altering chords hardly changes the overall emotion of music. For future work,
we will explore the power of transformer-based hierarchical VAE and
investigate music elements that have an impact on emotion (such as tempo, note
density) to realize emotion-controllable music generation.

## References

  * [1] S. Ji, X. Yang, and J. Luo, “A survey on deep learning for symbolic music generation: Representations, algorithms, evaluations, and challenges,” ACM Comput. Surv., 2023.  Just Accepted.
  * [2] S. Ji, J. Luo, and X. Yang, “A comprehensive survey on deep music generation: Multi-level representations, algorithms, evaluations, and future directions,” arXiv preprint arXiv:2011.06801, 2020.
  * [3] H. Lim, S. Rhyu, and K. Lee, “Chord generation from symbolic melody using BLSTM networks,” in Proceedings of the 18th International Society for Music Information Retrieval Conference, ISMIR 2017, Suzhou, China, October 23-27, 2017, pp. 621–627, 2017.
  * [4] Y.-C. Yeh, W.-Y. Hsiao, S. Fukayama, T. Kitahara, B. Genchel, H.-M. Liu, H.-W. Dong, Y. Chen, T. Leong, and Y.-H. Yang, “Automatic melody harmonization with triad chords: A comparative study,” Journal of New Music Research, vol. 50, no. 1, pp. 37–51, 2021.
  * [5] Y. Chen, H. Lee, Y. Chen, and H. Wang, “Surprisenet: Melody harmonization conditioning on user-controlled surprise contours,” in Proceedings of the 22nd International Society for Music Information Retrieval Conference, ISMIR 2021, pp. 105–112, 2021.
  * [6] S. Rhyu, H. Choi, S. Kim, and K. Lee, “Translating melody to chord: Structured and flexible harmonization of melody with transformer,” IEEE Access, vol. 10, pp. 28261–28273, 2022.
  * [7] L. Ferreira and J. Whitehead, “Learning to generate music with sentiment,” in Proceedings of the 20th International Society for Music Information Retrieval Conference, ISMIR, pp. 384–390, 2019.
  * [8] L. Ferreira, L. Lelis, and J. Whitehead, “Computer-generated music for tabletop role-playing games,” in Proceedings of the AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment, vol. 16, pp. 59–65, AAAI Press, 2020.
  * [9] L. N. Ferreira, L. Mou, J. Whitehead, and L. H. Lelis, “Controlling perceived emotion in symbolic music generation with monte carlo tree search,” in Proceedings of the AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment, vol. 18, pp. 163–170, 2022.
  * [10] S. Ji and X. Yang, “Emomusictv: Emotion-conditioned symbolic music generation with hierarchical transformer vae,” IEEE Transactions on Multimedia, pp. 1–13, 2023.
  * [11] C.-E. Sun, Y.-W. Chen, H.-S. Lee, Y.-H. Chen, and H.-M. Wang, “Melody harmonization using orderless nade, chord balancing, and blocked gibbs sampling,” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 4145–4149, 2021.
  * [12] T. Zeng and F. C. Lau, “Automatic melody harmonization via reinforcement learning by exploring structured representations for melody sequences,” Electronics, vol. 10, no. 20, p. 2469, 2021.
  * [13] S. Ji, X. Yang, J. Luo, and J. Li, “Rl-chord: Clstm-based melody harmonization using deep reinforcement learning,” IEEE Transactions on Neural Networks and Learning Systems, pp. 1–14, 2023.
  * [14] T. Takahashi and M. Barthet, “Emotion-driven harmonisation and tempo arrangement of melodies using transfer learning,” in ISMIR 2022 Hybrid Conference, 2022.
  * [15] W.-Y. Hsiao, J.-Y. Liu, Y.-C. Yeh, and Y.-H. Yang, “Compound word transformer: Learning to compose full-song music over dynamic directed hypergraphs,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, pp. 178–186, 2021.
  * [16] X. Yang, X. Lin, S. Suo, and M. Li, “Generating thematic chinese poetry using conditional variational autoencoders with hybrid decoders,” in IJCAI, pp. 4539–4545, 2018.
  * [17] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” arXiv preprint arXiv:1312.6114, 2013.
  * [18] D. Makris, K. R. Agres, and D. Herremans, “Generating lead sheets with affect: A novel conditional seq2seq framework,” in 2021 International Joint Conference on Neural Networks (IJCNN), pp. 1–8, IEEE, 2021.
  * [19] S. R. Bowman, L. Vilnis, O. Vinyals, A. M. Dai, R. Jozefowicz, and S. Bengio, “Generating sentences from a continuous space,” in 20th SIGNLL Conference on Computational Natural Language Learning, CoNLL 2016, pp. 10–21, ACL, 2016.
  * [20] E. Schonfeld, S. Ebrahimi, S. Sinha, T. Darrell, and Z. Akata, “Generalized zero-and few-shot learning via aligned variational autoencoders,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8247–8255, 2019.
  * [21] C. K. Sønderby, T. Raiko, L. Maaløe, S. K. Sønderby, and O. Winther, “Ladder variational autoencoders,” Advances in neural information processing systems, vol. 29, 2016.