# High-performance sampling of generic Determinantal Point Processes

Jack Poulson111<EMAIL_ADDRESS>Hodge Star Scientific Computing

###### Abstract

Determinantal Point Processes (DPPs) were introduced by Macchi [1] as a model
for repulsive (fermionic) particle distributions. But their recent
popularization is largely due to their usefulness for encouraging diversity in
the final stage of a recommender system [2].

The standard sampling scheme for finite DPPs is a spectral decomposition
followed by an equivalent of a randomly diagonally-pivoted Cholesky
factorization of an orthogonal projection, which is only applicable to
Hermitian kernels and has an expensive setup cost. Researchers [3, 4] have
begun to connect DPP sampling to $LDL^{H}$ factorizations as a means of
avoiding the initial spectral decomposition, but existing approaches have only
outperformed the spectral decomposition approach in special circumstances,
where the number of kept modes is a small percentage of the ground set size.

This article proves that trivial modifications of $LU$ and $LDL^{H}$
factorizations yield efficient direct sampling schemes for non-Hermitian and
Hermitian DPP kernels, respectively. Further, it is experimentally shown that
even dynamically-scheduled, shared-memory parallelizations of high-performance
dense and sparse-direct factorizations can be trivially modified to yield DPP
sampling schemes with essentially identical performance.

The software developed as part of this research, Catamari
[hodgestar.com/catamari] is released under the Mozilla Public License v2.0. It
contains header-only, C++14 plus OpenMP 4.0 implementations of dense and
sparse-direct, Hermitian and non-Hermitian DPP samplers.

##  1 Introduction

Determinantal Point Processes (DPPs) were first studied as a distinct class by
Macchi in the mid 1970’s [1, 5] as a probability distribution for the
locations of repulsive (fermionic) particles, in direct contrast with
Permanental – or, bosonic – Point Processes. Particular instances of DPPs,
representing the eigenvalue distributions of classical random matrix
ensembles, appeared in a series of papers in the beginning of the 1960’s [6,
7, 8, 9, 10] (see [11] for a review, and [12] for a computational perspective
on sampling classical eigenvalue distributions). Some of the early
investigations of (finite) DPPs involved their usage for uniformly sampling
spanning trees of graphs [13, 14] (see Figs. 1 and 2), non-intersecting random
walks [15], and domino tilings of the Aztec diamond [16, 17, 18, 19] (see
Figs. 3 and 4).

Let us recall that the class of immanants provides a generalization of the
determinant and permanent of a matrix by using an arbitrary character $\chi$
of the symmetric group $S_{n}$ to determine the signs of terms in the
summation:

$|A|^{\chi}\equiv\sum_{\sigma\in
S_{n}}\chi(\sigma)\prod_{i=1}^{n}A_{i,\sigma(i)}.$

Choosing the trivial character $\chi\equiv 1$ yields the permanent, while
choosing $\chi(\sigma)=\text{sign}(\sigma)$ provides the determinant. The
resulting generalization from Determinantal Point Processes and Permanental
Point Processes to Immanantal Point Processes was studied in [20].

Figure 1: Two samples from the uniform distribution over spanning trees of a
40 x 40 box of $\mathbb{Z}^{2}$. Each has likelihood $\text{exp}(-1794.24)$.

Figure 2: Two samples from the uniform distribution over spanning trees of a
10 x 10 section of hexagonal tiles. Each has likelihood
$\text{exp}(-299.101)$.

Figure 3: Two samples from the Aztec diamond DPP of order 10, defined via a
complex, non-Hermitian kernel based upon Kenyon’s formula over the Kasteleyn
matrix. The tiles all begin on dark squares and are oriented either: left
(blue), right (red), up (yellow), or down (green). Their likelihoods are both
$\text{exp}(-38.1231)$.

Figure 4: Two samples from the Aztec diamond DPP of order 80, defined via a
complex, non-Hermitian kernel based upon Kenyon’s formula over the Kasteleyn
matrix. Their likelihoods are both $\text{exp}(-2245.8)$. The diamond is large
enough to clearly display the arctic circle phenomenon discussed in [21] and
[19]: “when $n$ is sufficiently large, the shape of the central sub-region
becomes arbitrarily close to a perfect circle of radius $n/\sqrt{2}$ for all
but a negligible proportion of the tilings.”

###### Definition 1

A finite Determinantal Point Process is a random variable
$\mathbf{Y}\sim\text{DPP}(K)$ over the power set of a ground set
$\\{0,1,...,n-1\\}=[n]$ such that

$\mathbb{P}[Y\subseteq\mathbf{Y}]=\det(K_{Y}),$

where $K\in\mathbb{C}^{n\times n}$ is called the marginal kernel matrix and
$K_{Y}$ denotes the restriction of $K$ to the row and column indices of $Y$.

We can immediately observe that the $j$-th diagonal entry of a marginal kernel
$K$ is the probability of index $j$ being in the sample, $\mathbf{Y}$, so the
diagonal of every marginal kernel must lie in $[0,1]\subset\mathbb{R}$. A
characteristic requirement for a complex matrix to be admissible as a marginal
kernel is given in the following proposition, due to Brunel [22], which we
will provide an alternative proof of after introducing our factorization-based
sampling algorithm.

######  Proposition 1 (Brunel [22])

A matrix $K\in\mathbb{C}^{n\times n}$ is admissible as a DPP marginal kernel
iff

$(-1)^{|J|}\det(K-\mathds{1}_{J})\geq 0,\;\;\forall J\subseteq[n].$

We can also observe that, because determinants are preserved under similarity
transformations, there is a nontrivial equivalence class for marginal kernels
– which is to say, there are many marginal kernels defining the same DPP:

###### Proposition 2

The equivalence class of a DPP kernel $K\in\mathbb{C}^{n\times n}$ contains
its orbit under the group of diagonal similarity transformations, i.e.,

$\\{D^{-1}KD\,:\,D=\text{diag}(d),d\in(\mathbb{C}\setminus\\{0\\})^{n}\\}.$

If we restrict to the classes of complex Hermitian or real symmetric kernels,
the same statement holds with the entries of $d$ restricted to the circle
group $U(1)=\\{z\in(\mathbb{C},\times):|z|=1\\}$ and the scalar orthogonal
group $O(1)=\\{z\in(\mathbb{R},\times):|z|=1\\}$, respectively.

###### Proof.

Determinants are preserved under similarity transformations, so the
probability of inclusion of each subset is unchanged by global diagonal
similarity. That a unitary diagonal similarity preserves Hermiticity follows
from recognition that it becomes a Hermitian congruence. Likewise, a signature
matrix similarity transformation becomes a real symmetric congruence. ∎

In most works, with the notable exception of Kasteleyn matrix approaches to
studying domino tilings of the Aztec diamond [16, 17, 18, 19], the marginal
kernel is assumed Hermitian. And Macchi [1] showed (Cf. [23, 24]) that a
Hermitian matrix is admissibile as a marginal kernel if and only if its
eigenvalues all lie in $[0,1]$. The viewpoint of these eigenvalues as
probabilities turns out to be productive, as the most common sampling
algorithm for Hermitian DPPs, due to [24] and popularized in the machine
learning community by [2], produces an equivalent random orthogonal projection
matrix by preserving eigenvectors with probability equal to their eigenvalue:

######  Theorem 1 (Theorem 7 of [24])

Given any Hermitian marginal kernel matrix $K\in\mathbb{C}^{n\times n}$ with
spectral decomposition $Q\Lambda Q^{H}$, sampling from
$\mathbf{Y}\sim\text{DPP}(K)$ is equivalent to sampling from a realization of
the random DPP with kernel $Q_{:,\mathbf{Z}}Q_{:,\mathbf{Z}}^{H}$, where
$Q_{:,\mathbf{Z}}$ consists of the columns of $Q$ with
$\mathbb{P}[j\in\mathbf{Z}]=\Lambda_{j}$, independently for each $j$.

Such an orthogonal projection marginal kernel is said to define a
Determinantal Projection Process [24], or elementary DPP [2], and it is known
that the resulting samples almost surely have cardinality equal to the rank of
the projection (e.g., Lemma 17 of [24]). The algebraic specification in Alg. 1
of [2] of the Alg. 18 of [24] involved an $O(nk^{3})$ approach, where $k$ is
the rank of the projection kernel. In many important cases, such as uniformly
sampling spanning trees or domino tilings, $k$ is a large fraction of $n$, and
the algorithm has quartic complexity (assuming standard matrix
multiplication). In the words of [2]:

> Alg. 1 runs in time $O(nk^{3})$, where $k$ is the number of eigenvectors
> selected […] the initial eigendecomposition […] is often the computational
> bottleneck, requiring $O(n^{3})$ time. Modern multi-core machines can
> compute eigendecompositions up to $n\approx 1,000$ at interactive speeds of
> a few seconds, or larger problems up to $n\approx 10,000$ in around ten
> minutes.

A faster, $O(nk^{2})$ algorithm for sampling elementary DPPs was given as Alg.
2 of [25]. We will later show that this approach is equivalent to a small
modification of a diagonally-pivoted, rank-revealing, left-looking, Cholesky
factorization [26], where the pivot index is chosen at each iteration by
sampling from the probability distribution implied by the diagonal of the
remaining submatrix (which is maintained out-of-place).

As a brief aside, we recall that the distinction between up-looking, left-
looking, and right-looking factorizations is based upon their choice of loop
invariant. Given the matrix partitioning

$\begin{pmatrix}A_{1,1}&A_{1,2}\\\ A_{2,1}&A_{2,2}\end{pmatrix},$

where the diagonal of $A_{1,1}$ comprises the list of eliminated pivots: an
up-looking algorithm will have only updated $A_{1,1}$ (by overwriting it with
its factorization), a left-looking algorithm will also have overwritten
$A_{2,1}$ with the corresponding block column of the lower triangular factor
(and $A_{1,2}$ with its upper-triangular factor in the nonsymmetric case), and
a right-looking algorithm will additionally have overwritten $A_{2,2}$ with
the Schur complement resulting from the block elimination of $A_{1,1}$.

Researchers have begun proposing algorithms which directly sample from
Hermitian marginal kernels by sequentially deciding whether each index should
be in the result by sampling the Bernoulli process defined by the probability
of the item’s inclusion, conditioned on the combination of all of the explicit
inclusion and exclusion decisions made so far. For example, when deciding
whether to include index $j$ in the sample, we will have already partitioned
$[0,...,j-1]$ into a set of included indices, $F$, and excluded indices, $G$.
So we must sample index $j$ with probability
$\mathbb{P}[j\in\mathbf{Y}\,|\,F\subseteq\mathbf{Y},\,G\cap\mathbf{Y}=\emptyset]$.
Directly performing such sampling using formulae for marginal kernels of
conditional DPPs is referred to in [3] as sequential sampling (Cf. [4] for
greedy, maximum-likelihood inference).

The primary contribution of this manuscript is to show that sequential
sampling can be performed via a small modification of an unpivoted $LDL^{H}$
factorization process (where $L$ is unit lower-triangular and $D$ is real
diagonal) and to extend them to non-Hermitian marginal kernels using an
unpivoted $LU$ factorization; [27] motivates an algorithm for learning non-
Hermitian DPP kernels by their ability to incorporate both attraction and
repulsion between items.

It is then demonstrated that high-performance factorizations techniques [28,
29, 30] lead to orders of magnitude accelerations, and that sparse-direct
techniques [31, 32, 33] can yield further orders of magnitude speedups.

##  2 Prototype factorization-based DPP sampling

To describe factorization processes, we will extend our earlier notation that,
for an $n\times n$ matrix $K$ and an index subset $Y\subseteq[n]$, $K_{Y}$
refers to the restriction of $K$ to the row and column indices of $Y$. Given a
second index subset $Z\subseteq[n]$, $K_{Y,Z}$ will denote the restriction of
$K$ to the $|Y|\times|Z|$ submatrix consisting of the rows indices in $Y$ and
the column indices in $Z$. And we use the notation $[j:k]$ to represent the
integer range $\\{j,j+1,...,k-1\\}$.

Our derivation will make use of a few elementary propositions on the forms of
marginal kernels for conditional DPPs. These propositions will allow us to
define modifications to the pivots of an LU factorization so that the
resulting Schur complements correspond to the marginal kernel of the DPP over
the remaining indices, conditioned on the inclusion decisions of the indices
corresponding to the eliminated pivots.

###### Proposition 3

Given disjoint subsets $A,B\subseteq[n]$ of the ground set of a DPP with
marginal kernel $K$, almost surely

$\mathbb{P}[B\subseteq\mathbf{Y}|A\subseteq\mathbf{Y}]=\det(K_{B}-K_{B,A}K_{A}^{-1}K_{A,B}).$

###### Proof.

If $A\subseteq\mathbf{Y}$, then
$\det(K_{A})=\mathbb{P}[A\subseteq\mathbf{Y}]>0$ almost surely, so we may
perform a two-by-two block LU decomposition

$\begin{pmatrix}K_{A}&K_{A,B}\\\
K_{B,A}&K_{B}\end{pmatrix}=\begin{pmatrix}I&0\\\
K_{B,A}K_{A}^{-1}&K_{B}-K_{B,A}K_{A}^{-1}K_{A,B}\end{pmatrix}\begin{pmatrix}K_{A}&K_{A,B}\\\
0&I\end{pmatrix}.$

That $\det:GL(n,\mathbb{C})\mapsto(\mathbb{C},\times)$ is a homomorphism
yields

$\det(K_{A\cup B})=\det(K_{A})\det(K_{B}-K_{B,A}K_{A}^{-1}K_{A,B}).$

The result then follows from the definition of conditional probabilities for a
DPP:

$\mathbb{P}[B\subseteq\mathbf{Y}|A\subseteq\mathbf{Y}]=\frac{\mathbb{P}[A,B\subseteq\mathbf{Y}]}{\mathbb{P}[A\subseteq\mathbf{Y}]}=\frac{\det(K_{A\cup
B})}{\det(K_{A})}.$

∎

###### Proposition 4

Given disjoint $a\in[n]$ and $B\subset[n]$ for a ground set $[n]$ of a DPP
with marginal kernel $K$, almost surely

$\mathbb{P}[B\subset\mathbf{Y}|a\notin\mathbf{Y}]=\det(K_{B}-K_{B,a}(K_{a}-1)^{-1}K_{a,B}).$

###### Proof.

$\displaystyle\mathbb{P}[B\subset\mathbf{Y}|a\notin\mathbf{Y}]$
$\displaystyle=\frac{\mathbb{P}[a\notin\mathbf{Y}|B\subset\mathbf{Y}]\mathbb{P}[B\subset\mathbf{Y}]}{\mathbb{P}[a\notin\mathbf{Y}]}$
$\displaystyle=\frac{(1-\mathbb{P}[a\in\mathbf{Y}|B\subset\mathbf{Y}])\mathbb{P}[B\subset\mathbf{Y}]}{1-\mathbb{P}[a\in\mathbf{Y}]}$
$\displaystyle=\frac{\det(K_{B})-\det(K_{a\cup B})}{1-K_{a}}$
$\displaystyle=\det(K_{B})(1-K_{a,B}\frac{{K_{B}}^{-1}}{K_{a}-1}K_{B,a})$
$\displaystyle=\det(K_{B}-K_{B,a}(K_{a}-1)^{-1}K_{a,B}),$

where the last equality makes use of the Matrix Determinant Lemma. The
formulae are well-defined almost surely. ∎

Propositions 3 and 4 are enough to derive our direct, non-Hermitian DPP
sampling algorithm. But, for the sake of symmetry with Proposition 3, we first
generalize to set exclusion:

###### Proposition 5

Given disjoint subsets $A,B\subset\mathcal{Y}$, almost surely

$\mathbb{P}[B\subseteq\mathbf{Y}|A\subseteq\mathbf{Y}^{c}]=\det(K_{B}-K_{B,A}(K_{A}-I)^{-1}K_{A,B}).$

###### Proof.

The claim follows from recursive formulation of conditional marginal kernels
using the previous proposition. The resulting kernel is equivalent to the
Schur complement produced from the block LU factorization

$\begin{pmatrix}K_{A}-I&K_{A,B}\\\
K_{B,A}&K_{B}\end{pmatrix}=\begin{pmatrix}I&0\\\
K_{B,A}(K_{A}-I)^{-1}&K_{B}-K_{B,A}(K_{A}-I)^{-1}K_{A,B}\end{pmatrix}\begin{pmatrix}K_{A}-I&K_{A,B}\\\
0&I\end{pmatrix},$

as the subtraction of 1 from each eliminated pivot commutes with the outer
product updates. ∎

######  Theorem 2 (Factorization-based DPP sampling)

Given a (possibly non-Hermitian) marginal kernel matrix $K$ of order $n$, an
LU factorization of $K$ can be modified as in Algorithm
LABEL:lst:unblocked_dpp to almost surely provide a sample from
$\text{DPP}(K)$. This algorithm involves roughly $\frac{2}{3}n^{3}$ floating-
point operations, and the likelihood of any returned sample will be given by
the product of the absolute value of the diagonal of the result. If the
marginal kernel is Hermitian, the work can be roughly halved by exploiting
symmetry in the Schur complement updates.

###### Proof.

We first demonstrate, by induction, that our factorization algorithm samples
the DPP generated by the marginal kernel $K$. The loop invariant is that, at
the start of the iteration for pivot index $j$, $A_{[j:n]}$ represents the
equivalence class of kernels for the DPP over indices $[j:n]$ conditioned on
the inclusion decisions for indices $0,...,j-1$.

Since the diagonal entries of a kernel matrix represent the likelihood of the
corresponding index being in the sample, the loop invariant implies that index
$j$ is kept with the correct conditional probability. Prop. 3 shows that the
loop invariant is almost surely maintained when the Bernoulli draw is
successful, and Prop. 4 handles the alternative. Thus, the loop invariant
holds almost surely, and, upon completion, the proposed algorithm samples each
subset with the correct probability by sequentially iterating over each index,
making an inclusion decision with the appropriate conditional probability.

The likelihood of a sample produced by the algorithm is thus the product of
the likelihoods of the results of the Bernoulli draws: when a draw for a
diagonal entry $p_{j}$ is successful, its probability was $p_{j}$, and, when
unsuccessful, $1-p_{j}$. In both cases, the multiplicative contribution is the
absolute value of the final state of the $j$’th diagonal entry. ∎

Algorithm 1: Unblocked, right-looking, non-Hermitian DPP sampling. For a
sample $Y$, the returned matrix $A$ will contain the in-place $LU$
factorization of $K-\mathds{1}_{Y^{C}}$, where $\mathds{1}_{Y^{C}}$ is the
diagonal indicator matrix for the entries not in the sample. Symmetry can be
exploited in the outer products when the kernel is Hermitian. ⬇ 1sample := [];
A := K 2for j in range(n):  3 sample.append(j) if Bernoulli($A_{j}$) else
$A_{j}\,-\\!\\!=\,1$ 4 $A_{[j+1:n],\,j}\,\,/\\!\\!=\,\,A_{j}$ 5
$A_{[j+1:n]}\,\,-\\!\\!=\,\,A_{[j+1:n],\,j}\,A_{j,\,[j+1:n]}$ 6return sample,
A

Theorem 2 provides us with another interpretation of the generic DPP kernel
admissibility condition of [22]:

###### Proposition 1 above.

When Algorithm LABEL:lst:unblocked_dpp produces a sample $Y$, the resulting
upper and strictly-lower triangular portions of the resulting matrix
respectively contain the $U$ and strictly-lower portion of the unit-diagonal
$L$ from the LU factorization of $K-\mathds{1}_{Y^{C}}$. And we have:

$\det(K-\mathds{1}_{Y^{C}})=\det(U)=\prod_{j=0}^{n-1}U_{j}=\prod_{j\in
Y}p_{j}\prod_{j\in Y^{C}}(p_{j}-1),$

where $p_{j}$ is the inclusion probability for index $j$, conditioned on the
inclusion decisions of indices $0,...,j-1$.

One can inductively show, working backwards from the last pivot, that
$(-1)^{|Y^{C}|}\det(K-\mathds{1}_{Y^{C}})$ always being non-negative is
equivalent to all potential pivots produced by our sampling algorithm, the set
of conditional inclusion probabilities, living in $[0,1]$. Otherwise, there
would exist an index inclusion decision change which would not change the sign
of $\det(K-\mathds{1}_{Y^{C}})$. ∎

By replacing the Bernoulli samples of algorithm 2 with the maximum-likelihood
result for each index inclusion, we arrive at an analogous (approximate)
maximum-likelihood inference algorithm.

In both cases, the same specializations that exist for modifying an unpivoted
LU factorization into a Cholesky or (unpivoted) $LDL^{H}$ or $LDL^{T}$
factorization apply to our DPP sampling algorithms. And as we will see in the
next two sections, so do high-performance dense and sparse-direct
factorization techniques.

As a brief aside, for real matrices, assuming standard matrix multiplication
algorithms, the highest-order terms for the operation counts of dense Gaussian
Elimination and an MRRR-based [34] Hermitian eigensolver are
$\frac{2}{3}n^{3}$ and $\frac{10}{3}n^{3}$ [35]. But the coefficient for the
Hermitian eigensolver is misleadingly small, as the initial phase of a
Hermitian eigensolver traditionally involves a unitary reduction to Hermitian
tridiagonal form that, due to only modest potential for data reuse, executes
significantly less efficiently than traditional dense factorizations. So-
called Successive Band Reduction techniques [36, 37, 38] were therefore
introduced as a gambit for trading higher operation counts for decreased data
movement and – as a consequence – increased performance.

######  Theorem 3 (Factorization-based elementary DPP sampling)

Given a marginal kernel matrix $K$ of order $n$ which is an orthogonal
projection of rank $k$, performing $k$ steps of diagonally-pivoted Cholesky
factorization, where each pivot index is chosen by sampling from the diagonal
as in Algorithm LABEL:lst:unblocked_elementary_dpp, is equivalent to sampling
from $\text{DPP}(K)$. This approach has complexity $O(nk^{2})$ and almost
surely completes and returns a sample of cardinality $k$; the likelihood of
the resulting sample is the square of the product of the first $k$ diagonal
entries of the partially factored matrix.

###### Proof.

Alg. LABEL:lst:unblocked_elementary_dpp is, up to permutation, algebraically
equivalent to the sampling phase of Alg. 2 of [25]: for the sake of
simplicity, it assumes the Gramian is preformed rather than forming it on the
fly from the factor. That the probability of a cardinality $k$ sample $Y$ is,
almost surely, equal to the product of the squares of the first $k$ diagonal
entries of the result follows from recognizing that the lower triangle of the
top-left $k\times k$ submatrix will be the Cholesky factor of $K_{Y}$, and

$\mathbb{P}[Y=\mathbf{Y}]=\mathbb{P}[Y\subseteq\mathbf{Y}]=\det(K_{Y})=\det(L_{Y}L_{Y}^{H})=\det(L_{Y})\det(L_{Y}^{H})=\prod_{j=0}^{k-1}A_{j}^{2}.$

∎

Algorithm 2: Unblocked, left-looking, diagonally-pivoted, Cholesky-based
sampling of a Hermitian Determinantal Projection Process. For a sample $Y$,
the returned matrix will contain the in-place Cholesky factorization of
$K_{Y}$. The computational cost is $O(nk^{2})$. ⬇ 1$A\,:\\!=\,K$;
$d\,:\\!=\,\text{diag}(K)$; orig_indices := [0:n]  2for j in range(k):  3 #
Sample pivot index and permute 4 Draw index t from $[j:n]$ with probability
$d_{t}/(k-j)$ 5 Perform Hermitian swap of indices j and t of $A$ 6 Swap
positions j and t of orig_indices and $d$ 7 $A_{j}:=\sqrt{d_{j}}$ 8 if j == k
- 1:  9 break 10 # Form new column and update diagonal 11
$A_{[j+1:n],\,j}\,\,-\\!\\!=\,A_{[j+1:n],\,[0:j]}A_{j,\,[0:j]}^{H}$ 12 for t
in range(j+1, n):  13 $A_{t,j}\,\,/\\!\\!=\,\,A_{j}$ 14
$d_{t}\,\,-\\!\\!=\,\,|A_{t,j}|^{2}$ 15return orig_indices[0:k], $A_{[0:k]}$

The permutations in Alg. LABEL:lst:unblocked_elementary_dpp were introduced to
solidify the connection to a traditional diagonally-pivoted Cholesky
factorization. There is the additional benefit of simplifying the usage of
Basic Linear Algebra Subprograms (BLAS) [39, 40, 41] calls for the
matrix/vector products. But the bulk of the work of this approach will be in
rank-one updates, which have essentially no data reuse, and are therefore not
performant on modern machines, where the peak floating-point performance is
substantially faster than what can be read directly from main memory. The
Linear Algebra PACKage (LAPACK) [28] was introduced in 1990 as proof that
dense factorizations can be recast in terms of matrix/matrix multiplications;
we present analogues in the following section, as well as multi-core, tiled
extensions similar to [29, 30].

##  3 High-performance, dense, factorization-based DPP sampling

The main idea of LAPACK [28] is to reorganize the computations within dense
linear algebra algorithms so that as much of the work as possible is recast
into composing matrices with nontrivial minimal dimensions. Basic operations
rich in such matrix compositions – typically referred to as Level 3 BLAS [41]
– are the building blocks of LAPACK. In most cases, such block sizes range
from roughly 32 to 256, with 64 to 128 being most common.

In the case of triangular factorizations, such as Cholesky, $LDL^{H}$, and LU,
a blocked algorithm can be produced from the unblocked algorithm by carefully
matricizing each of the original operations. In the case of a right-looking LU
factorization without pivoting, one arrives at Alg. LABEL:lst:blocked_lu,
where $A_{J_{1}}$ takes the place of the scalar pivot and must be factored –
typically, using the unblocked algorithm being generalized – before its
components are used to solve against $A_{J_{2},\,J_{1}}$ and
$A_{J_{1},\,J_{2}}$.

Algorithm 3: Blocked LU factorization without pivoting. The functions triu and
unit_tril mirror MATLAB notation and respectively return the upper-triangular
and unit-diagonal lower-triangular restrictions of their input matrices. ⬇ 1j
:= 0  2while j < n:  3 $\text{bsize}\,:=\,\min(\text{blocksize},n-j)$ 4
$J_{1}=[j:j+\text{bsize}]$; $J_{2}=[j+\text{bsize}:n]$ 5
$A_{J_{1}}=\text{unblocked\\_lu}(A_{J_{1}})$ 6
$A_{J_{2},\,J_{1}}\,\,:=\,\,A_{J_{2},\,J_{1}}\,\,\text{triu}(A_{J_{1}})^{-1}$
7
$A_{J_{1},\,J_{2}}\,\,:=\,\,\text{unit\\_tril}(A_{J_{1}})^{-1}\,\,A_{J_{1},\,J_{2}}$
8 $A_{J_{2}}\,\,-\\!\\!=\,\,A_{J_{2},\,J_{1}}\,\,\,A_{J_{1},\,J_{2}}$ 9
$j\,\,+\\!\\!=\,\,\text{bsize}$ 10return sample, A

The correct form of these solves can be derived via the relationship:

$\begin{pmatrix}A_{1}&A_{1,\,2}\\\
A_{2,\,1}&A_{2}\end{pmatrix}=\begin{pmatrix}L_{1}&0\\\
L_{2,\,1}&L_{2}\end{pmatrix}\begin{pmatrix}U_{1}&U_{1,\,2}\\\
0&U_{2}\end{pmatrix}=\begin{pmatrix}L_{1}U_{1}&L_{1}U_{1,\,2}\\\
L_{2,\,1}U_{1}&L_{2,\,1}U_{1,\,2}+L_{2}U_{2}\end{pmatrix},$

where we used shorthand of the form $A_{1,2}$ to represent
$A_{J_{1},\,J_{2}}$. An unblocked factorization can be used to compute the
triangular factors $L_{1}$ and $U_{1}$ of the diagonal block $A_{1}$, and then
triangular solves of the form $L_{2,1}:=A_{2,1}U_{1}^{-1}$ and
$U_{1,2}:=L_{1}^{-1}A_{1,2}$ yield the two panels. After forming the Schur
complement $S_{2}:=A_{2}-L_{2,1}U_{1,2}$, the problem has been reduced to an
LU factorization with fewer variables – that of $L_{2}U_{2}=S_{2}$.
Asymptotically, all of the work is performed in the outer-product updates
which form the Schur complements.

The conversion of Alg. LABEL:lst:unblocked_dpp, an unblocked DPP sampler, into
Alg. LABEL:lst:blocked_dpp, a blocked DPP sampler, is essentially identical to
the formulation of Alg. LABEL:lst:blocked_lu from an unblocked LU
factorization. We emphasize that, while it is well-known that $LU$
factorizations without pivoting fail on large classes of nonsingular matrices
– for example, any matrix with a zero in the top-left position – the analogue
for DPP sampling succeeds almost surely for any marginal kernel.

Algorithm 4: Blocked, factorization-based, non-Hermitian DPP sampling. For a
sample $Y$, the returned matrix $A$ will contain the in-place $LU$
factorization of $K-\mathds{1}_{Y^{C}}$, where $\mathds{1}_{Y^{C}}$ is the
diagonal indicator matrix for the entries not in the sample, $\mathbf{Y}$. The
usual specializations from $LU$ to $LDL^{H}$ factorization applies if the
kernel is Hermitian. ⬇ 1sample := []; A := K; j := 0  2Replace line 5 of main
loop of Alg. 3 with:  3
$\text{subsample},A_{J_{1}}=\text{unblocked\\_dpp}(A_{J_{1}})$ 4
sample.append(subsample + j)  5return sample, A

Beyond the order-of-magnitude improvement provided by such algorithms, even on
a single core of a modern computer, they also simplify the incorporation of
parallelism. Lifting algorithms into blocked form allows for each core to be
dynamically assigned tasks corresponding to individual updates of a tile [29,
30], a roughly $\text{tile\\_size}\times\text{tile\\_size}$submatrix which, in
our experiments, was typically most performant for $\text{tile\\_size}=256$.

Within the context of Alg. LABEL:lst:blocked_dpp, our tiled algorithm assigns
each `unblocked_dpp` call an individual OpenMP 4.0 [42] task which depends
upon the last Schur complement update of its tile. The triangular solves
against $U_{J_{1}}$ are split into submatrices mostly of the form
$\text{tile\\_size}\times\text{block\\_size}$, those of the solves against
$L_{J_{1}}$ are roughly of the transpose dimensions, while the Schur
complement tasks are roughly of size
$\text{tile\\_size}\times\text{tile\\_size}$. In each case, tile reads and
writes are scheduled using dependencies on previous updates of the same tile.

The performance of such a dynamically-scheduled parallelization of the
Hermitian specialization of Alg. LABEL:lst:blocked_dpp is demonstrated for
arbitrary dense, Hermitian marginal kernels on a 16-core Intel i9-7960x in
Fig. 5. The performance of a similarly parallelized unpivoted $LDL^{H}$
factorization is similarly plotted, and one readily observes that their
runtimes are essentially identical. The runtime of Alg. LABEL:lst:blocked_dpp
for arbitrary, complex, non-Hermitian marginal kernels is similarly shown in
Fig. 6. For purpose of comparison, we note that the double-precision High-
Performance LINPACK benchmark [43] achieves roughly 1 TFlop/second on this
machine.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0$$0.5$$1$$1.5$Matrix size (in
thousands)TFlop/s32-bit DPP32-bit $LDL^{H}$64-bit DPP64-bit $LDL^{H}$

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10^{-2}$$10^{-1}$Matrix size (in
thousands)Seconds32-bit DPP32-bit $LDL^{H}$64-bit DPP64-bit $LDL^{H}$

Figure 5: Dense, real $LDL^{H}$-based DPP sampling performance; tile sizes of
$128$ and $256$ were used for matrix sizes less than or equal to, and greater
than, $2000$, respectively.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0.5$$1$$1.5$$2$Matrix size (in
thousands)TFlop/s32-bit DPP64-bit DPP

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10^{-2}$$10^{-1}$$10^{0}$Matrix size (in
thousands)Seconds32-bit DPP64-bit DPP

Figure 6: Dense, complex $LU$-based DPP sampling performance. For single-
precision, a tile size of 128 was used up to matrix sizes of 3000, and tile
sizes of 256 were used thereafter. For double-precision, the switch occurred
above matrices of size 4000.

It is worth emphasizing the compounding performance gains from both the
formulation of DPP sampling as a small modification of a level-3 BLAS focused
dense matrix factorization and from the 16-way parallelism. When combined, a
factor of 2500x speedup is observed relative to the timings on the same
machine of DPPy v0.1.0 [44] when sampling a spectrally-preprocessed arbitrary,
dense Hermitian $5000\times 5000$ matrix. Similar speedups exist relative to
the timings of both the “sequentially thinned” and spectrally-preprocessed
algorithms of [3]. In the case of DPPy, it would be fair to attribute at least
an order of magnitude of the performance difference to its implementation
being standard Python code rather than optimized C++.

In the case of DPPs where the expected number of samples,
$k\equiv\text{trace}(K)$, is much less than the ground set size, $n$,
Derezinski et al. [45] have proposed a rejection-sampling approach with a
setup cost of $n\cdot\text{polylog}(n)\cdot\text{poly}(k)$ and a subsequent
sampling cost of $\text{poly}(k)$. Ignoring polylogarithmic terms, their setup
cost is $O(nk^{6}+k^{9})$ and the subsequent sampling cost is $O(k^{6})$.
Large speedups can therefore be expected when $k\ll\sqrt{n}$. But the examples
in Figs. 1–4 of this article all involve $k$ being a modest fraction of $n$,
e.g., $k\approx n/4$ in the case of the Aztec diamond, so the rejection
sampler would involve a prohibitive $O(n^{9})$ setup and $O(n^{6})$ subsequent
sampling cost.

Timings for the $\mathbb{Z}^{2}$ and hexagonal-tiling Uniform Spanning Tree
DPPs can be extracted from Fig. 5 via the formulae $n\approx 2d^{2}$ and
$n\approx 6d^{2}$, respectively, where $d$ is the order of the grid. And
timings for the uniform domino tilings can be recovered from Fig. 6 via the
formula $n\approx 4d^{2}$. When computing results for the latter, it was
noticed that sampling domino tilings from the Kenyon formula DPP in single-
precision, the results are essentially always inconsistent once the diamond
size exceeds 60, but no inconsistencies have yet been observed with double-
precision sampling. See Fig. 7 for an example.

Figure 7: Corrupted single-precision sample from an Aztec diamond of size 80.
Missing tiles are clearly visible in the image.

Such a numerical instability in an unpivoted dense matrix factorization would
lead any numerical analyst to wonder if dynamic pivoting can mitigate error
accumulation. Indeed, our factorization-based approach frees us to perform
arbitrary diagonal pivoting, as long as the pivot is decided before its
corresponding Bernoulli draw. While it will be the subject of future work, the
author conjectures that maximum-entropy pivoting, that is, pivoting towards a
diagonal entry with conditional probability as close to $\frac{1}{2}$ as
possible, would be the most beneficial, as it maximizes the magnitude of the
smallest possible pivot.

Before moving on to sparse-direct DPP sampling, we demonstrate the multiple
orders of magnitude speedup that are possible for Determinantal Projection
Processes of sufficiently low rank. Unlike our other experiments, Fig. 8 only
executes on a single core, as the unblocked, left-looking rank-revealing
elementary DPP sampling approach of Alg. LABEL:lst:unblocked_elementary_dpp
spends the majority of its time in matrix/vector multiplication. Studying the
benefits of parallelizations of elementary DPP samplers is left for future
work.

$10$$50$$100$$250$$500$$0$$20$$40$RankGFlop/s64-bit generic DPP sampler64-bit
elementary DPP sampler

$10$$50$$100$$250$$500$$10^{-2}$$10^{0}$RankSeconds64-bit generic DPP
sampler64-bit elementary DPP sampler

Figure 8: Sequential, dense, real $LDL^{H}$-based elementary DPP sampling
performance for a ground set of size $5000$ and varying rank. Multithreaded
BLAS was explicitly disabled to ensure execution on a single core of the Intel
i9-7960x.

##  4 Sparse-direct DPP sampling

As was demonstrated in Fig. 5, the performance of a generic Hermitian DPP
sampler can be made to match that of a high-performance Hermitian matrix
factorization. This correspondence should not be a surprise, as we have shown
that matrix factorizations can be trivially modified to yield Hermitian and
non-Hermitian DPP samplers alike. The same insight applies to the conversion
of sparse-direct, unpivoted $LDL^{H}$ factorization into sparse-direct
Hermitian DPP samplers. Such an implementation, mirroring many of the
techniques from CholMod [33], was implemented within Catamari.

The high-level approach is to use a dynamically scheduled – via OpenMP 4.0
task scheduling – multifrontal method [46, 32] when the computational
intensity of the factorization is deemed high enough after determining the
number of nonzeros and of the triangular factor and the number of operations
required to compute it. If the arithmetic count and intensity is not
sufficiently high, an up-looking, scalar algorithm [33] is used instead.

The multifrontal method performs the bulk of its work processing an
(elimination) tree of dense frontal matrices in a post-ordering. Each such
frontal matrix for a supernode $s$ initially takes the form

$\begin{pmatrix}A_{s}&*\\\ A_{\text{struct}(s),s}&0\end{pmatrix},$

where $A_{s}$ is the diagonal block of the input matrix corresponding to
supernode $s$, $A_{\text{struct}(s),s}$ is the submatrix of $A$ below the
supernodal diagonal block that contains at least one nonzero in each row, and
the $*$ denotes the upper-right quadrant not needing to be accessed due to
Hermiticity. Processing such a front involves adding any child Schur
complements onto the front, factoring $A_{s}=L_{s}D_{s}L_{s}^{H}$, solving
against $D_{s}L_{s}^{H}$ to replace $A_{\text{struct}(s),s}$ with
$L_{\text{struct}(s),s}$, and replacing the bottom-right quadrant with the
Schur complement $-L_{\text{struct}(s),s}D_{s}L_{\text{struct}(s),s}^{H}$.

This multifrontal process is then parallelized using a tile-based, dynamically
scheduled DPP sampler or dense factorization on the diagonal blocks, combined
with nested OpenMP tasks launched for each (relaxed) supernode’s [32] subtree
[47].

While the author is not aware of any well-studied sparse marginal kernel
matrices, [48] proposed a scheme for learning data-sparse marginal kernels
defined via Kronecker products and evaluated their technique on the Amazon
baby registry dataset of [49]. While DPP learning schemes are beyond the scope
of this manuscript, the dramatic factorization speedups that can be supplied
by sparse-direct solvers suggests the possibility of learning an entrywise
sparse marginal kernel matrix, perhaps via incorporating entrywise soft-
thresholding into an iterative learning scheme [49].

We therefore implemented tandem sparse-direct factorization and sparse-direct
Hermitian DPP samplers to connect performance results of synthetic sparse DPPs
and discretized Partial Differential Equation solvers. Such a statically-
pivoted sparse-direct $LDL^{H}$ factorization is also the critical
computational kernel of primal-dual Interior Point Methods [50, 51], where the
first-order optimality conditions are symmetric quasi-semidefinite [52, 53].

Performance results for the dynamically scheduled sparse-direct solver on a 3D
Helmholtz equation discretized with trilinear, hexahedral elements Fig. 9 are
shown in Fig. 10; results for a 2D DPP analogue, where the marginal kernel
takes the form of a scaled 2D Laplacian, are shown in Fig. 11.

Figure 9: A 2D slice for a solution of a 3D Helmholtz equation discretized
with trilinear, hexahedra elements and PML boundary conditions for a Gaussian
point source near the boundary and a converging lens near the center of the
domain.

$20$$40$$60$$80$$100$$0$$0.2$$0.4$$0.6$$\mathbf{d}$
($\mathbf{n=(d+1)^{3}}$)TFlop/sfp64 (with symbolic)fp64 (w/o symbolic)

$20$$40$$60$$80$$100$$10^{-1}$$10^{0}$$10^{1}$$\mathbf{d}$
($\mathbf{n=(d+1)^{3}}$)Secondsfp64 (with symbolic)fp64 (w/o symbolic)

Figure 10: Sparse-direct performance on trilinear, hexahedral discretization
of 3D Helmholtz equation with Perfectly Matched Layer boundary conditions.

(a) Sample with log-likelihood of -27472.2.

(b) Greedy maximum-likelihood sample with log-likelihood of -26058.

Figure 11: Samples from $-\sigma\Delta$ over a 200 x 200 grid, where
$\sigma=0.72$.

In the case of the 2D Laplacian sparse-direct DPP sampling over a $200\times
200$ grid, the timings were roughly 0.01 seconds on the 16-core i9-7960x.
Extrapolating from Fig. 5 to a dense matrix size of $40,000$, it is clear that
several more orders of magnitude of efficiency are gained with the sparse-
direct formulation.

We close this section by noting that static pivoting should apply equally well
to non-Hermitian sparse-direct DPP sampling, as the nature of small pivots
being, by definition, rare, suggests a certain degree of stability. As in the
case of dense DPP sampling, probabilistic error analysis is warranted.

##  5 Conclusions

A unified, generic framework for sequentially sampling both Hermitian and non-
Hermitian Determinantal Point Processes directly from their marginal kernel
matrices was presented. The prototype algorithm consisted of a small
modification of an unpivoted LU factorization – in the Hermitian case, it
reduces to a modification of an unpivoted $LDL^{H}$ factorization – and high-
performance, tiled algorithms were presented in the dense case, and an
analogue of CholMod [33] was presented for the sparse case. Further, it was
shown that both greedy, maximum-likelihood sampling and elementary DPP
sampling can be understood and efficiently implemented using small
modifications of traditional matrix factorization techniques.

In addition to generalizing sampling algorithms from Hermitian to non-
Hermitian marginal kernels, the proposed approaches were implemented within
the open source, permissively licensed, Catamari package and shown to lead to
orders of magnitude speedups, even in the dense regime. The Hermitian sparse-
direct DPP sampler leads to further asymptotic speedups. Future work includes
theoretical and empirical exploration of the stability of dynamic pivoting
techniques, both for dense and sparse-direct factorization-based DPP sampling,
including the incorporation of maximum-entropy pivoting for large instances of
Kenyon-formula Aztec domino tilings.

Further exploration of the performance tradeoffs between left-looking and
right-looking elementary DPP sampling, and at which rank it becomes beneficial
to use our generic sampling approach, is warranted – especially on multi-core
and GPU-accelerated architectures. And, lastly, the incorporation of a tiled
extension of [54] for converting a Hermitian L-ensemble kernel into a marginal
kernel via $K=I-(L+I)^{-1}$ [25] is a natural extension of our proposed
techniques and could be expected to require the equivalent of 3 samples worth
of time due to the computational complexity.

But perhaps the most pressing future direction is to investigate the potential
of learning sparse marginal kernels on a benchmark problem, such as the Amazon
baby registry of [49]. Ideally such a scheme would expose a roughly block-
diagonal covariance matrix which captures item clusters and allows for
sufficiently sparse triangular factors.

## Reproducibility

The experiments in this paper were performed using version 0.2.5 of Catamari,
with Intel’s Math Kernel Library BLAS on an Ubuntu 18.04 workstation with 64
GB of RAM and an Intel i9-7960x processor. Catamari can be downloaded from
https://gitlab.com/hodge_star/catamari. Figs. 1 and 2 were generated using the
TIFF output from example driver
`example/uniform_spanning_tree.cc`, while Figs. 3, 4 and 7 similarly used
`example/aztec_diamond.cc`.

Figs. 5 and 6 were generated using data from `example/dense_dpp.cc` and
`example/dense_factorization.cc`; we emphasize that explicitly setting both
`OMP_NUM_THREADS` and `MKL_NUM_THREADS` to 16 on our 16-core machine led to
significant performance improvements (as opposed to the default exploitation
of hyperthreading). And Fig. 8 was produced with
`example/dense_elementary_dpp.cc` – with `OMP_NUM_THREADS` and
`MKL_NUM_THREADS` explicitly set to one. The algorithmic `block_size` was also
set to a value at least as large as the rank so that a left-looking approach
was used throughout.

Figs. 9 and 10 were generated from `example/helmholtz_3d_pml.cc`, with the
2D-slice of the former taken slightly away from the center plane. Fig. 11 was
output from the TIFF support of
`example/dpp_shifted_2d_negative_laplacian.cc`.

## Acknowledgements

The author would like to sincerely thank Guillaume Gautier and Rémi Bardenet
for detailed comments on an early draft of this manuscript, as well as their
concise explanation of the correctness of the factorization-based sampler at
https://dppy.readthedocs.io/en/latest/finite_dpps/exact_sampling.html#id17. I
would also like to thank Jennifer A. Gillenwater for suggesting the
$I-(L+I)^{-1}$ conversion from an L-ensemble and Alex Kulesza for answering
some of my early DPP questions. Lastly, I am thankful for the insightful
comments and suggestions of the two anonymous referees: they significantly
improved the quality of the manuscript.

## References

  * [1] Odile Macchi.  The coincidence approach to stochastic point processes.  Advances in Applied Probability, 7(1):83–122, 1975.
  * [2] Alex Kulesza and Ben Taskar.  Determinantal point processes for machine learning.  Foundations and Trends in Machine Learning, 5(2–3):123–286, 2012\.
  * [3] Claire Launay, Bruno Galerne, and Agnès Desolneux.  Exact Sampling of Determinantal Point Processes without Eigendecomposition.  arXiv e-prints, page arXiv:1802.08429, Feb 2018.
  * [4] Laming Chen and Guoxin Zhang.  Fast greedy map inference for determinantal point process to improve recommendation diversity.  In NeurIPS, 2018.
  * [5] Odile Macchi.  The Fermion Process — A Model of Stochastic Point Process with Repulsive Points, pages 391–398.  Springer Netherlands, Dordrecht, 1977.
  * [6] Madan L. Mehta and Michel Gaudin.  On the density of eigenvalues of a random matrix.  Nuclear Physics, 18:420–427, 1960.
  * [7] Freeman J. Dyson.  Statistical theory of the energy levels of complex systems. i.  Journal of Mathematical Physics, 3(1):140–156, 1962.
  * [8] Freeman J. Dyson.  Statistical theory of the energy levels of complex systems. ii.  Journal of Mathematical Physics, 3(1):157–165, 1962.
  * [9] Freeman J. Dyson.  Statistical theory of the energy levels of complex systems. iii.  Journal of Mathematical Physics, 3(1):166–175, 1962.
  * [10] Jean Ginibre.  Statistical Ensembles of Complex, Quaternion and Real Matrices.  J. Math. Phys., 6:440–449, 1965.
  * [11] Persi Diaconis.  Patterns in eigenvalues: the 70th Josiah Willard Gibbs lecture.  Bull. Amer. Math. Soc., 40:155–178, 2003.
  * [12] Alan Edelman and N. Raj Rao.  Random matrix theory.  Acta Numerica, 14:233–297, 2005.
  * [13] Robert Burton and Robin Pemantle.  Local characteristics, entropy and limit theorems for spanning trees and domino tilings via transfer-impedances.  The Annals of Probability, 21(3):1329–1371, 1993.
  * [14] Itai Benjamini, Russell Lyons, Yuval Peres, and Oded Schramm.  Uniform spanning forests.  Ann. Probab., 29(1):1–65, 02 2001.
  * [15] Kurt Johansson.  Determinantal processes with number variance saturation.  Communications in Mathematical Physics, 252(1):111–148, Dec 2004\.
  * [16] Peter W. Kasteleyn.  The statistics of dimers on a lattice : I. The number of dimer arrangements on a quadratic lattice.  Physica, 27:1209–1225, December 1961.
  * [17] Harold N. V. Temperley and Michael E. Fisher.  Dimer problem in statistical mechanics-an exact result.  The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics, 6(68):1061–1063, 1961.
  * [18] Peter W. Kasteleyn.  Dimer statistics and phase transitions.  Journal of Mathematical Physics, 4(2):287–293, 1963.
  * [19] Sunil Chhita, Kurt Johansson, and Benjamin Young.  Asymptotic domino statistics in the aztec diamond.  Ann. Appl. Probab., 25(3):1232–1278, 06 2015.
  * [20] Persi Diaconis and Steven N. Evans.  Immanants and finite point processes.  Journal of Combinatorial Theory, Series A, 91(1):305–321, 2000\.
  * [21] William Jockusch, James Propp, and Peter Shor.  Random Domino Tilings and the Arctic Circle Theorem.  arXiv Mathematics e-prints, page math/9801068, Jan 1998.
  * [22] Victor-Emmanuel Brunel.  Learning signed determinantal point processes through the principal minor assignment problem.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 7365–7374. Curran Associates, Inc., 2018.
  * [23] Alexander Soshnikov.  Determinantal random point fields.  Russian Mathematical Surveys, 55(5):923–975, oct 2000.
  * [24] J. Ben Hough, Manjunath Krishnapur, Yuval Peres, and Balint Virag.  Determinantal processes and independence.  Probab. Surveys, 3:206–229, 2006.
  * [25] Jennifer A. Gillenwater.  Approximate Inference for Determinantal Point Processes.  PhD thesis, University of Pennsylvania, 2014.
  * [26] Nicholas J. Higham.  Analysis of the Cholesky decomposition of a semi-definite matrix.  In Reliable Numerical Computation, pages 161–185. University Press, 1990.
  * [27] Mike Gartrell, Victor-Emmanuel Brunel, Elvis Dohmatob, and Syrine Krichene.  Learning nonsymmetric determinantal point processes.  CoRR, abs/1905.12962, 2019.
  * [28] Edward Anderson, Zhaojun Bai, Jack J. Dongarra, Anne Greenbaum, Alan McKenney, Jeremy Du Croz, Sven Hammarling, James Demmel, Christian Bischof, and Danny Sorensen.  LAPACK: A portable linear algebra library for high-performance computers.  In Proceedings of the 1990 ACM/IEEE Conference on Supercomputing, Supercomputing ’90, pages 2–11, Los Alamitos, CA, USA, 1990\. IEEE Computer Society Press.
  * [29] Ernie Chan, Field G. Van Zee, Paolo Bientinesi, Enrique S. Quintana-Orti, Gregorio Quintana-Orti, and Robert van de Geijn.  SuperMatrix: A multithreaded runtime scheduling system for algorithms-by-blocks.  In Proceedings of the 13th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, PPoPP ’08, pages 123–132, New York, NY, USA, 2008. ACM.
  * [30] Alfredo Buttari, Julien Langou, Jakub Kurzak, and Jack Dongarra.  A class of parallel tiled linear algebra algorithms for multicore architectures.  Parallel Comput., 35(1):38–53, January 2009.
  * [31] Robert Schreiber.  A new implementation of sparse Gaussian elimination.  ACM Trans. Math. Softw., 8(3):256–276, September 1982.
  * [32] Cleve Ashcraft and Roger Grimes.  The influence of relaxed supernode partitions on the multifrontal method.  ACM Trans. Math. Softw., 15(4):291–309, December 1989.
  * [33] Yanqing Chen, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam.  Algorithm 887: Cholmod, supernodal sparse cholesky factorization and update/downdate.  ACM Trans. Math. Softw., 35(3):22:1–22:14, October 2008.
  * [34] Inderjit S. Dhillon, Beresford N. Parlett, and Christof Vömel.  The design and implementation of the mrrr algorithm.  ACM Trans. Math. Softw., 32(4):533–560, December 2006.
  * [35] Bruce Hendrickson, Elizabeth Jessup, and Christopher Smith.  Toward an efficient parallel eigensolver for dense symmetric matrices.  SIAM J. Sci. Comput., 20(3):1132–1154, January 1999.
  * [36] Christian H. Bischof, Bruno Lang, and Xiaobai Sun.  A framework for symmetric band reduction.  ACM Trans. Math. Softw., 26(4):581–601, December 2000.
  * [37] Azzam Haidar, Hatem Ltaief, and Jack Dongarra.  Parallel reduction to condensed forms for symmetric eigenvalue problems using aggregated fine-grained and memory-aware kernels.  In Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’11, pages 8:1–8:11, New York, NY, USA, 2011. ACM.
  * [38] Grey Ballard, James Demmel, and Nicholas Knight.  Avoiding communication in successive band reduction.  ACM Trans. Parallel Comput., 1(2):11:1–11:37, February 2015.
  * [39] Charles L. Lawson, Richard J. Hanson, Fred T. Krogh, and David R. Kincaid.  Algorithm 539: Basic Linear Algebra Subprograms for Fortran usage [F1].  ACM Trans. Math. Softw., 5(3):324–325, September 1979.
  * [40] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard J. Hanson.  An extended set of FORTRAN Basic Linear Algebra Subprograms.  ACM Trans. Math. Softw., 14(1):1–17, March 1988.
  * [41] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and I. S. Duff.  A set of level 3 Basic Linear Algebra Subprograms.  ACM Trans. Math. Softw., 16(1):1–17, March 1990.
  * [42] OpenMP Architecture Review Board.  OpenMP application program interface version 4.0, 2013.
  * [43] Jack J. Dongarra, Piotr Luszczek, and Antoine Petitet.  The LINPACK benchmark: Past, present, and future, 2002.
  * [44] Guillaume Gautier, Rémi Bardenet, and Michal Valko.  DPPy: Sampling determinantal point processes with Python.  CoRR, abs/1809.07258, 2018.
  * [45] Michal Derezinski, Daniele Calandriello, and Michal Valko.  Exact sampling of determinantal point processes with sublinear time preprocessing.  CoRR, abs/1905.13476, 2019.
  * [46] Iain S. Duff and John K. Reid.  The multifrontal solution of indefinite sparse symmetric linear.  ACM Trans. Math. Softw., 9(3):302–325, September 1983.
  * [47] Jonathan Hogg and Jennifer Scott.  New parallel sparse direct solvers for multicore architectures.  Algorithms, 6(4):702–725, 2013.
  * [48] Zelda Mariet and Suvrit Sra.  Kronecker determinantal point processes.  CoRR, abs/1605.08374, 2016.
  * [49] Jennifer Gillenwater, Alex Kulesza, Emily Fox, and Ben Taskar.  Expectation-maximization for learning determinantal point processes.  In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2, NIPS’14, pages 3149–3157, Cambridge, MA, USA, 2014. MIT Press.
  * [50] Anna Altman and Jacek Gondzio.  HOPDM – a higher order primal-dual method for large scale linear programming.  European Journal of Operational Research, 66(1):158 – 160, 1993\.
  * [51] Jacek Gondzio.  Multiple centrality corrections in a primal-dual method for linear programming.  Computational Optimization and Applications, 6(2):137–156, Sep 1996\.
  * [52] Robert Vanderbei.  Symmetric quasidefinite matrices.  SIAM Journal on Optimization, 5(1):100–113, 1995.
  * [53] Alan George, Khakim Ikramov, and Andrey B. Kucherov.  Some properties of symmetric quasi-definite matrices.  SIAM Journal on Matrix Analysis and Applications, 21(4):1318–1323, 2000.
  * [54] Paolo Bientinesi, Brian Gunter, and Robert A. van de Geijn.  Families of algorithms related to the inversion of a symmetric positive definite matrix.  ACM Trans. Math. Softw., 35(1):3:1–3:22, July 2008.


# Generating Nontrivial Melodies for Music as a Service

###### Abstract

We present a hybrid neural network and rule-based system that generates pop
music. Music produced by pure rule-based systems often sounds mechanical.
Music produced by machine learning sounds better, but still lacks hierarchical
temporal structure. We restore temporal hierarchy by augmenting machine
learning with a temporal production grammar, which generates the music’s
overall structure and chord progressions. A compatible melody is then
generated by a conditional variational recurrent autoencoder.

The autoencoder is trained with eight-measure segments from a corpus of 10,000
MIDI files, each of which has had its melody track and chord progressions
identified heuristically.

The autoencoder maps melody into a multi-dimensional feature space,
conditioned by the underlying chord progression. A melody is then generated by
feeding a random sample from that space to the autoencoder’s decoder, along
with the chord progression generated by the grammar. The autoencoder can make
musically plausible variations on an existing melody, suitable for recurring
motifs. It can also reharmonize a melody to a new chord progression, keeping
the rhythm and contour.

The generated music compares favorably with that generated by other academic
and commercial software designed for the music-as-a-service industry.

##  1 Introduction

Computer-generated music has started to expand from its pure artistic and
academic roots into commerce. Companies such as Jukedeck and Amper offer so-
called music as a service, by analogy with software as a service. However,
their melodies, when present at all, often just arpeggiate the underlying
chord.

We extend this approach by generating music with both chord progressions and
interesting, nontrivial melodies. We expand a song structure such as
$AA^{\prime}BA$ into a harmonic plan, and then add a melody compatible with
this structure and harmony. This compatibility uses a chord-melody
relationship found by applying machine learning to a corpus of MIDI
transcriptions of pop music (Figure 1).

ML training Crawled MIDI files Harmonic analysis Learned modelMelodyChords ML
generation Chord grammarGenerated melody Figure 1: Machine learning (ML)
workflow for generating music from a MIDI corpus.

Prior research is discussed in section 2. Harmonic analysis is detailed in
sections 3 and 4. Hierarchy generation and melody generation are described in
section 5.

##  2 Related Work

Recent approaches to machine composition use neural networks (NNs), hoping to
approximate how humans compose. Chu et al [5] generate a melody with a
hierarchical NN that encodes a composition strategy for pop music, and then
accompany the melody with chords and percussion. However, this music lacks
hierarchical temporal structure. Boulanger-Lewandowski et al [3] investigate
hierarchical temporal dependencies and long-term polyphonic structure.
Inspired by how an opening theme often recurs at a song’s end, they detect
patterns with a recurrent temporal restricted Boltzmann machine (RTRBM). This
can represent more complicated temporal distributions of notes. Similarly,
Huang and Wu [10] generate structured music with a 2-layer Long Short Term
Memory (LSTM) network. Although the resulting music often sounds plausible, it
cannot produce clearly repeated melodic themes, just like a Markov resynthesis
of the text of the famous poem “Jabberwocky” is unlikely to replicate the
identical opening and closing stanzas of the original. Despite the LSTM
network’s theoretical capability of long-term memory, it fails to generalize
to arbitrary time lengths [8], and its generated melodies remain
unimaginative.

In these approaches, tonic chords dominate, and melody is little more than
arpeggiation. To avoid this banality, we work in reverse. We first create
structure and chords, and then fit melody to that. This mimics how classical
western Roman-numeral harmony is taught to beginners: only after one has the
underlying chord sequence, can one explain the melody in terms of chord tones,
passing tones, appoggiaturas, and so on.

##  3 Melody identification

For pop music, a catchy and memorable melody is crucial. To generate melodies
that sound less robotic than those generated by other algorithms, we use
machine learning. To create a learning database, we started with a corpus of
10,000 MIDI files [16], from which we extracted useful training data (melodies
that sound vivid or fun). In particular, the training data was eight-measure
excerpts labelled as melody and chords. We thus had to identify which of a
MIDI file’s several tracks contained the melody. To do so, we assigned each
track the sum of a rubric score and an entropy score. Whichever track scored
highest was declared to be the melody. (Ties between high-scoring tracks were
broken arbitrarily, because they were usually due to several tracks having
identical notes, differing only in which instrument played them.)

###  3.1 Rubric Score

Our rubric considered attributes such as instrumentation, note density, and
pitch range.

We first considered a track’s instrument name (MIDI meta-event FF 04). Certain
instruments are more common for melody, such as violin or flute. Others are
more likely to be applied as accompaniment or long sustained notes, such as
low brass. A third category is likely used as unpitched percussion. The
instrument’s category then adjusted the rubric’s score.

We also considered the track’s note density, how often at least one note is
sounding (between corresponding MIDI note-on and note-off events), as a
fraction of the track’s full duration. A track scored higher if this was
between 0.4 and 0.8, a typical value for pop melodies.

Finally we considered pitch range, because we observed that pop melodies often
lie between C3 and C5. The score was higher for a pitch range between C3 and
C6, to exclude bass tracks from consideration.

The values for these attributes were chosen based on manual inspection of 100
files in the corpus.

###  3.2 Entropy score

We empirically observed that melody tracks often have a greater variety of
pitches than other tracks. Thus, to quantify how varied, complex, and dynamic
a track was, we calculated each track’s entropy

$H(X)=-\sum\nolimits_{i=1}^{12}{P(x_{i})\log P(x_{i})}$ (1)

where $x_{i}$ represents the event that a particular note in the octave is $i$
semitones from the pitch C, and $P(x_{i})$ represents that event’s
probability. Higher entropy corresponds to a greater number of distinct
pitches.

###  3.3 Evaluation

To measure how well this scoring identified melody tracks, we manually tagged
the melody track of 160 randomly selected MIDI files. Comparing the scored
prediction to this ground truth showed that the error rate was 15%.

##  4 Chord Detection

To identify the chords in a MIDI file, we considered three aspects of how pop
music differs from genres like classical music. First, chord inversions (where
the lowest note is not the chord’s root) are rare. When a new chord is
presented, it is often in root position: most pop songs have a clear melody
line and bass line [14], and the onset of a new chord is marked with the
chord’s root in that bass line. Second, chords may contain extensions
(seventh), substitutions (flattened fifth), doublings, drop voicings (changing
which octave a pitch sounds in), and omissions (third or fifth). Although such
modifications complicate the task of functional harmony analysis, this is not
a concern for our application. Third, new chord onsets are often at the start
of a measure; rarely are there more than two chords per measure. Combining
these observations led us to the following chord detection algorithm.

We first partition the song into segments with constant time signatures.
(these are explicitly stated as MIDI meta messages). Then each segment is
evenly divided into bins, where we try to match the entire bin to a chord.
Because chords have different durations, we try different bin lengths: half a
measure, one measure, and two measures. Then for each bin, containing all the
notes sounding during that time interval, we add all these notes to a set that
is matched against a fixed collection of chords, based on how close the
pitches are, with a cost function:

Chord Detection: Cost

1:function BestChordInBin($Pitches$)

2: $Root\leftarrow\text{Lowest note starting before first upbeat}$

3: $Chords\leftarrow\text{All chords, as array of intervals}$

4:return $\operatorname*{argmin}_{C\in
Chords}\\{\textsc{Cost}(Pitches,C,Root)\\}$

5:function Cost($Pitches,Chord,Root$)

6: $PitchCost\leftarrow 0$

7: for $P\in Pitches$ do

8: $interval\leftarrow\text{No. semitones of }P\text{ from }Root$

9: $d\leftarrow\min_{voice\in Chord}\\{dist(interval,voice)\\}$

10: $PitchCost\leftarrow PitchCost+d$

11: $ChordCost\leftarrow 0$

12: for $voice\in Chord$ do

13: $d\leftarrow\min_{P\in Pitches}\\{dist(P-Root,voice)\\}$

14: $ChordCost\leftarrow ChordCost+d$

15:return $PitchCost+ChordCost$

Each chord’s cost is the sum of the distance of the nearest interval in the
chord (from the root) to each interval in the input pitches, and the distance
of the nearest interval in the input pitches (from its “root”) to each
interval in the chord, based on some definition of distance. The cost function
then returns the lowest-cost chord.

Defining the distance in terms of mere pitch difference in semitones would be
simple, but performs poorly. For example, matching the pitch set $[C,E,G]$ to
the chord $[C,E\flat,G\flat]$ would yield a cost of two, which is far too low.
Instead, our distance function reflects how compatible intervals are. The
unison is the most compatible, with distance zero; fourths and fifths are
next, with distance one (Table 1). This conveniently handles omitted-fifth
chords, because the chord’s root matches the omitted fifth with a distance of
only one.

Distance in semitones | Compatibility distance
---|---
0 | 0
1 | 6
2 | 2
3 | 2
4 | 2
5 | 1
6 | 4
7 | 1
Table 1: Interval compatibility.

[options=-Oscores/= -c -s 1.0]scores/fly-me-to-the-moon

Figure 2: Example of chord detection.

Figure 2 demonstrates chord detection on the song Fly me to the moon. The bin
size is half a measure, yielding 8 identified chords. The $A7/aug/\flat 9$
chord resulted from the accompaniment notes $A,G,B\flat(\text{flat
ninth}),C\sharp,E$, and the melody notes $A,G,F(\text{augmented fifth})$.

##  5 Writing Music

To output pieces with audibly hierarchical structure, we start with the
harmonic structure produced by a temporal generative grammar. Then an
autoencoder recurrent neural network (RNN) generates a melody compatible with
this harmonic scaffold. The RNN learns to play using the chord’s notes, with
occasional surprising non-chord tone decorations such as passing tones and
appoggiaturas.

###  5.1 Generating Melody

We first search for a representation of the melody using ML. This is
traditionally done by an autoencoder, a pair of NNs that maps high-dimensional
input data to and from a lower-dimensional space. Although this dimensionality
reduction can eliminate perfect mappings, this turns out not to be a problem
because the subspace of “pleasant” music within all possible musics is
sufficiently small. Thus, the autoencoder can extract the pleasant content and
map only that into the representation space.

It is tempting to feed a random point from the representation space to the
autoencoder’s decoder, and observe how much sense it makes of that point.
However, because one cannot control the shape of the distribution of melody
representations, one cannot guarantee that a given point from the
representation space would be similar to those seen by the decoder during
training. Thus, the vanilla autoencoder architecture [2] is not viable as a
generative model. We propose the following improvements for generating
melodies:

  1. 1.

_Condition the NN on the chord progression._ The chord progression is provided
to the NN at every level, so when reproducing a melody, the decoder has access
to both the representation and the chord progression. This is useful because a
melody has rhythmic information, intervallic content, and contour. The decoder
can ignore the separately provided harmonic information, and use only the
melody’s other aspects. This also lets the representation remain constant
while altering the chord progression, so the NN can adapt a melody to a
changed chord progression, such as what happens when a key changes from minor
to major.

  2. 2.

_Add a stochastic layer._ Autoencoders which learn a stochastic representation
are called variational autoencoders, and perform well in generative modelling
of images [11]. The representation is not deterministic. We assume a
particular (Gaussian) distribution in the representation space, and then train
the NN to transform this distribution to match the distribution of input
melodies in their high dimensional space. This ensures that we can take a
random sampling of the representation space following its associated
probability distribution, then feed it through the decoder and expect a melody
similar to the set of musically sensible melodies.

  3. 3.

_Use recurrent connections._ Pop music has many time-invariant elements,
especially at time scales below a few beats. A recurrent NN shares the same
processing infrastructure for note sequences starting at different times, and
thereby accelerates learning.

  4. 4.

_Normalize all other notes relative to the tonic._ Pop music is also largely
pitch invariant, insofar as a song transposed by a few semitones still sounds
perceptually similar. The NN ignores the song’s key and considers the tonic
pitch to be abstract, as far as pitches in melody and chords are concerned.

####  5.1.1 Implementation

The input melody is quantized to sixteenth notes. Only sections with an
unchanging duple or quadruple meter are kept. The melody is converted to a
series of one-hot vectors, whose slots represent offsets from the tonic in the
range of $-16$ to $16$ semitones, with one more slot representing silence.
There is also an attack channel, where a value of 1 indicates that the note is
being rearticulated at the current time step. The encoding for chords supports
up to two chords per measure, and uses a one-hot vector for scale degrees and
separate boolean channels for chord qualities (Table 2). (Note that because
this encoding uses just seven Roman-numeral symbols, it does not try to
represent chords outside the current mode. Before training, we removed from
the corpus the few songs that contained significant occurrences of this.) We
use the basic triad form for each chord identified using techniques from
section 4, marking compatible chord qualities. For example, G7 is encoded by
marking a 1 in the $Maj$ and $Pwr$ columns. (The chord quality encoding could
be extended to seventh and ninth chords.) The table’s gray rows are data the
network is conditioned on, while the other rows are input data that the
network tries to reproduce. For an 8-measure example, the input and output
vector size is $35\times 8\times 16=4480$, and the conditional vector size is
$8\times 16+5\times 16+8=216$.

{57mm[16x] | -16 | … | 16 | Silent | Attack |  }22.5mm[x8]
---|---|---|---|---|---|---
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| I | … | VI | VII | Silent |
| Pwr | Maj | Min | Dim | Aug |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| Major | Dorian | … | Locrian | Jazz Minor |
Table 2: An encoding of 8 measures (see section 5.1.1).

The network has 24 recurrent layers, 12 each for the encoder and decoder
(Figure 3). Drawing on ideas of deep residual learning from computer vision
[9], we make additional connections from the input to every third hidden
layer. To improve learning, the network accesses both the original melody and
the transformed results from previous layers during processing. The
conditional part (chords and mode) is also provided to the network at every
recurrent layer, as extra incoming connections.

Input vector 600 hidden unit 600 hidden unit 600 hidden unit 600 hidden unit
600 hidden unit 600 hidden unit 600 hidden unit 600 hidden unit 600 hidden
unit 600 hidden unit 600 hidden unit 600 hidden unit chords and mode 8x300 to
1x1200 800 FC 800 FC Latent distribution 600 FC 600 hidden unit 600 hidden
unit 600 hidden unit 600 hidden unit 600 hidden unit 600 hidden unit 600
hidden unit 600 hidden unit 600 hidden unit 600 hidden unit 600 hidden unit
600 hidden unit chords and mode 128x300 to 128x35 Output
vector6006003006006003006006003006006002400800800MeanStandard
DeviationSampling600600300600600300600600300600600300 Figure 3: Network
architecture. Rectangles are bidirectional recurrent neural network cells.
Ellipses are strided time-convolution cells. Rounded rectangles are fully
connected (FC) layers. Numbered arrows indicate a connection’s dimension.

The network is implemented in Tensorflow, a machine learning library for rapid
prototyping and production training [1]. It was trained for four days on an
Nvidia Tesla K80 GPU. We used Gated Recurrent Units [4] to build the
bidirectional recurrent layer and Exponential Linear Units [7] as activation
functions. These significantly accelerate training while simplifying the
network [6, 7]. Figure 4 shows the training error (the sum of model
reproduction errors) and the difference of the latent distribution from a unit
Gaussian distribution, as measured by Kullback-Leibler divergence [12]. The
network’s input data (available at https://goo.gl/VezNNA) is a set of MIDI
songs from various online sources. Our harmonic analysis converted this to
$1.9\times 10^{6}$ measures of melodies and corresponding chords. We
implemented KL warm up, because that is crucial to learning for a variational
autoencoder [15]. But instead of linearly scaling the KL term for this, we
found that a sigmoid reduced the network’s reproduction loss.

$0$$2$$4$$6$$\cdot 10^{5}$$200$$400$$600$Reproduction lossKL divergence Figure
4: Training error and Kullback-Leibler divergence of the NN. The horizontal
axis indicates how many training segments have elapsed ($\times 10^{5}$).
Initial outliers have been removed.

###  5.2 Generating Hierarchy and Chords

Hierarchy and chords are generated simultaneously, using a temporal generative
grammar [13], modified to suit the harmonies of pop music, and extended to
enable repeated motifs with variations. The original temporal generative
grammar has notions of sharing by binding a section to a symbol. For example,
the rule

$\text{let}\ x=\text{I in I }M5(x)\text{ I }M5(x)\text{ I},$ (2)

where $M5$ indicates modulating to the $5^{\text{th}}$ degree, would expand to
five sections, with the second and fourth identical because $x$ is reused. We
extend this by having symbols $x$ carry along a number: $x_{1},x_{2},...$.
Different subscripts of the same symbol still expand to the same chord
progression, but denote slightly different latent representations when
generating corresponding melodies for those sections. The latent
representations corresponding to $x_{i>1}$ are derived from that of $x_{1}$ by
adding random Gaussian perturbations. This yields variations on the original
melody.

###  5.3 Training Examples in the Representation Space

We randomly chose 130 songs from the training set, fed them through the
network, and performed t-SNE analysis on the resulting 130 locations in the
representation space. Although a melody maps to a distribution in the
representation space, Figure 5 plots only each distribution’s mean, for visual
clarity. This t-SNE analysis effectively reduces the 800-dimensional
representation space into a low-dimensional human-readable format [17]. (A
larger interactive visualization of 1,680 songs is at
https://composing.ai/tsne.)

Figure 5: Example melodies in a t-SNE plot of the representation space.

Two songs that are both in the techno genre, Indica by Jushi and Control by
Junk Project, are indeed very near in the t-SNE plot, almost overlapping.
Excerpts from them show that both have a staccato rhythm with notes landing on
the upbeat, and have similar contours (Figure 6).

[options=-Oscores/= -c -s 1.2]scores/indigo [options=-Oscores/= -c -s
1.2]scores/control

Figure 6: Four-bar excerpts from the songs Indica (top) and Control (bottom).

###  5.4 Reharmonizing Melody

We hypothesized that, when building the neural network architecture, providing
the chord progression to both the encoder and the decoder would not preserve
that information in the representation space, thus saving space for rhythmic
nuances and contour. To test this hypothesis, we gave the network songs
disjoint from the training set and collected their representations. We then
fed these representations along with a new chord progression to the network.
We hoped that it would respond by generating a melody that was harmonically
compatible with the new chord progression, while still resembling the original
melody. We demonstrate this with the Chinese folk song Jasmine Flower, in a
genre unfamiliar to the NN (Figure 7). Note that we supplied the chords in
Figure 7 (bottom), for which the NN filled in the melody. The network
flattened the E, A, and B, by observing that the chord progression looked
minor. This is typically how a human would perform the reharmonization,
demonstrating the network’s comprehension of how melody and harmony interact.

Although the NN struggled to reproduce the melody, it provided interesting
modifications. The grace notes in measure 6 could be due to similar ones in
the training set, or due to vacillation between the A$\natural$ from the
representation and the A$\flat$ from the chord conditioning.

[options=-Oscores/= -c -s 1.2]scores/jasmine-flower-original
[options=-Oscores/= -c -s 1.2]scores/jasmine-flower-reharm

Figure 7: The song Jasmine Flower with original chords (top), and adapted to a
new chord progression (bottom).

###  5.5 Examples of generated Melodies

Because an entire multi-section composition cannot fit here, we merely show
excerpts from two shorter examples.

Figure 8 and Figure 9 demonstrate melodies generated from points in the
representation space that are not near any particular previously known melody.
Structure is evident in Figure 8: measures 1–3 present a short phrase, and
measure 4 leads to the next four measures, which recapitulate the first three
measures with elaborate variation. Figure 9 shows an energetic melody where
the grammar only produced C minor chords. Although the final two measures
wander off, the first six have consistent style and natural contour.

[options=-Ogenerated_scores/= -c -s 1.2]generated_scores/sample33-2

Figure 8: Generated melody for a grammar-generated chord progression.

[options=-Ogenerated_scores/= -c -s 1.2]generated_scores/sample33-16

Figure 9: Generated melody for an extended C minor chord.

##  6 Conclusion and future work

We have combined generative grammars for structure and harmony with a NN,
trained on a large corpus, to emit melodies compatible with a given chord
progression. This system generates compositions in a pop music style whose
melody, harmony, motivic development, and hierarchical structure all fit the
genre.

This system is currently limited by assuming that the input data’s chords are
in root position. More sophisticated chord detection would still let it
exploit the relative harmonic rigidity of popular music. Also, by
investigating the representation found by the NN, meaning could be assigned to
some of its 800 dimensions, such as intensity, consonance, and contour. This
would let us boost or attenuate a given melody along those dimensions.

Acknowledgements. The authors are grateful to Mark Hasegawa-Johnson for
overall guidance, and to the anonymous reviewers for insights into both
philosophical issues and minute details.

## References

  * [1] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, and X. Zheng.  Tensorflow: a system for large-scale machine learning.  In Proc. USENIX Conf. Operating Systems Design and Implementation, pages 265–283. USENIX Association, 2016.
  * [2] Bengio.  Learning deep architectures for AI.  Foundations and Trends in Machine Learning, 2(1):1–127, 2009.
  * [3] N. Boulanger-Lewandowski, Y. Bengio, and P. Vincent.  Modeling temporal dependencies in high-dimensional sequences: Application to polyphonic music generation and transcription.  arxiv.org/abs/1206.6392, 2012.
  * [4] K. Cho, B. van Merrienboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arxiv.org/abs/1406.1078, 2014.
  * [5] H. Chu, R. Urtasun, and S. Fidler.  Song from PI: a musically plausible network for pop music generation.  arxiv.org/abs/1611.03477, 2016.
  * [6] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arxiv.org/abs/1412.3555, 2014.
  * [7] D. Clevert, T. Unterthiner, and S. Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arxiv.org/abs/1511.07289, 2015.
  * [8] A. Graves, G. Wayne, and I. Danihelka.  Neural turing machines.  arxiv.org/abs/1410.5401, 2014.
  * [9] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  arxiv.org/abs/1512.03385, 2015.
  * [10] A. Huang and R. Wu.  Deep learning for music.  arxiv.org/abs/1606.04930, 2016.
  * [11] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  arxiv.org/abs/1312.6114, 2013.
  * [12] S. Kullback and R. A. Leibler.  On information and sufficiency.  Annals of Mathematical Statistics, 22(1):79–86, 1951.
  * [13] D. Quick and P. Hudak.  A temporal generative graph grammar for harmonic and metrical structure.  In Proc. International Computer Music Conference, pages 177–184, 2013.
  * [14] M. P. Ryynänen and A. P. Klapuri.  Automatic transcription of melody, bass line, and chords in polyphonic music.  Computer Music Journal, 32(3):72–86, 2008.
  * [15] C. K. Sønderby, T. Raiko, L. Maaløe, S. K. Sønderby, and O. Winther.  Ladder variational autoencoders.  arxiv.org/abs/1602.02282, 2016.
  * [16] Y. Teng and A. Zhao.  Composing.AI.  http://composing.ai/, 2017.
  * [17] L. van der Maaten and G. E. Hinton.  Visualizing high-dimensional data using t-SNE.  Journal of Machine Learning Research, 9:2579–2605, 2008.