Neural Nonnegative Matrix Factorization for Hierarchical Multilayer Topic Modeling

TL;DR

Neural NMF is a novel layered approach based on nonnegative matrix factorization that uncovers hierarchical structures in data, improving interpretability and performance over existing methods across various datasets.

Contribution

This paper introduces Neural NMF, a neural network-based recursive NMF method for hierarchical topic modeling, enhancing structure discovery and interpretability.

Findings

Neural NMF outperforms existing hierarchical NMF methods.

It provides better hierarchical structure learning.

It improves interpretability of topics.

Abstract

We introduce a new method based on nonnegative matrix factorization, Neural NMF, for detecting latent hierarchical structure in data. Datasets with hierarchical structure arise in a wide variety of fields, such as document classification, image processing, and bioinformatics. Neural NMF recursively applies NMF in layers to discover overarching topics encompassing the lower-level features. We derive a backpropagation optimization scheme that allows us to frame hierarchical NMF as a neural network. We test Neural NMF on a synthetic hierarchical dataset, the 20 Newsgroups dataset, and the MyLymeData symptoms dataset. Numerical results demonstrate that Neural NMF outperforms other hierarchical NMF methods on these data sets and offers better learned hierarchical structure and interpretability of topics.

Figures

Tables

Table 1. Table 1 . Reconstruction error / classification accuracy for various supervision levels and layer structures on the synthetic dataset. Entry marked with ∗ * corresponds to the experiment represented in Figure 2 , and entry marked with ∗ ⁣ ∗ ** corresponds to the experiment represented in Figure 3 .

	Layers	Hier. NMF	Deep NMF	Neural NMF
Unsuper.	1	0.053 / 0.111	0.031 / 0.111	0.029 / 1
	2	0.399 / 0.222	0.414 / 0.222	0.310 / 0.995 ∗
	3	0.860 / 0.356	0.838 / 0.356	0.492 / 1
Semisuper.	1	0.049 / 0.933	0.031 / 0.947	0.042 / 1
	2	0.374 / 0.926	0.394 / 0.911	0.305 / 1
	3	0.676 / 0.930	0.733 / 0.930	0.496 / 0.990 ∗∗
Supervised	1	0.052 / 0.960	0.042 / 0.962	0.042 / 1
	2	0.311 / 0.984	0.310 / 0.984	0.307 / 1
	3	0.495 / 1	0.494 / 1	0.498 / 1

Table 2. Table 2 . Classification accuracies of each layer given for a two layer unsupervised experiment and a two layer semisupervised experiment on the subsampled 20 Newsgroups dataset.

	Layer	Hier. NMF	Deep NMF	Neural NMF
Unsuper.	1	0.593	0.638	0.604
	2	0.507	0.444	0.532
Semisuper.	1	0.593	0.642	0.690
	2	0.546	0.536	0.654

Table 3. Table 3 . Topic keywords for layer 1 of the subsampled 20 Newsgroups dataset produced by Neural NMF.

	Topic 1	Topic 2	Topic 3	Topic 4	Topic 5	Topic 6	Topic 7	Topic 8	Topic 9	Topic 10
1	drive	new	sale	bike	msg	geb	space	humanist	people	jesus
2	apple	dj	offer	dod	food	pitt	launch	article	gun	god
3	video	following	drive	motorcycle	used	banks	nasa	politics	israel	people
4	mhz	st	mb	helmet	object	gordon	shuttle	omran	fbi	bible
5	card	computer	best	games	dietz	n3jxp	moon	bedouin	government	christians
6	graphics	fm	color	stadium	disease	chastity	ideas	backcountry	us	israel
7	mac	mower	software	got	responses	cadre	orbit	speaking	say	christian
8	sound	sold	disks	baseball	patients	skepticism	lunar	liar	arab	religion
9	powerbook	battery	game	year	diet	dsl	centaur	absood	jews	order
10	know	remains	shipping	players	epilepsy	shameful	medical	john	killed	rosicrucian

Table 4. Table 4 . Topic keywords for layer 2 of the subsampled 20 Newsgroups dataset produced by Neural NMF.

	Topic 1	Topic 2	Topic 3	Topic 4	Topic 5	Topic 6
1	drive	sale	bike	space	people	jesus
2	apple	new	dod	know	gun	god
3	video	computer	motorcycle	msg	israel	people
4	graphics	offer	helmet	launch	article	bible
5	mhz	drive	stadium	diet	fbi	israel
6	card	shipping	got	cost	government	christians
7	sound	mb	games	disease	guns	christian
8	mac	sell	baseball	heard	us	religion
9	powerbook	color	players	used	mr	jews
10	projector	best	uhhhh	centaur	dear	rosicrucian

Table 5. Table 5 . Topic keywords for layer 1 of the subsampled 20 Newsgroups dataset produced by HNMF.

	Topic 1	Topic 2	Topic 3	Topic 4	Topic 5	Topic 6	Topic 7	Topic 8	Topic 9	Topic 10
1	mail	mb	sale	team	geb	know	fbi	people	israel	jesus
2	modem	drive	games	runs	pitt	files	space	gun	arab	god
3	mac	card	offer	win	gordon	program	koresh	hudson	israeli	christians
4	keyboard	ram	shipping	games	banks	motorcycle	government	us	lebanon	bible
5	internal	vga	game	year	dsl	bike	fire	say	soldiers	christian
6	apple	hard	power	pitching	chastity	postscript	time	guns	peace	christ
7	computer	video	best	game	cadre	format	launch	moral	arabs	life
8	use	floppy	system	last	njxp	file	handheld	morality	jews	law
9	cable	color	sound	fans	skepticism	question	jmd	way	lebanese	love
10	software	mhz	super	rbi	shameful	dod	com	data	occupied	jews

Table 6. Table 6 . Topic keywords for layer 2 of the subsampled 20 Newsgroups dataset produced by Hierarchical.

	Topic 1	Topic 2	Topic 3	Topic 4	Topic 5	Topic 6
1	mb	sale	team	geb	israel	jesus
2	drive	games	runs	pitt	arab	god
3	mail	offer	win	gordon	israeli	christians
4	modem	shipping	year	banks	people	bible
5	mac	game	games	dsl	government	christian
6	know	power	pitching	chastity	lebanon	christ
7	keyboard	system	game	cadre	fbi	life
8	computer	best	last	njxp	soldiers	law
9	internal	disks	fans	skepticism	peace	love
10	apple	sound	rbi	shameful	arabs	jews

Equations

$$\underset{\mathbf{A}\ge 0,\mathbf{S}\ge 0}{\min }||\mathbf{X}-\mathbf{A}\mathbf{S}|{|}^2.$$

$$\mathrm{rank}(\mathbf{A}\mathbf{S})\le \min (\mathrm{rank}(\mathbf{A}),\mathrm{rank}(\mathbf{S}))\le k.$$

$${\mathbf{W}}_{n,m}=\{\begin{array}{ll}1& {\mathbf{X}}_{n,m}\text{ is known}\\ 0& \text{otherwise,}\end{array}$$

$${\mathbf{Z}}_{:,m}=\{\begin{array}{ll}{1}_P& \text{label for object }m\text{ is known}\\ 0& \text{otherwise}.\end{array}$$

$$\underset{\mathbf{A}\ge 0,\mathbf{S}\ge 0,\mathbf{B}\ge 0}{\min }\underset{\text{Reconstruction Error}}{\underbrace{\Vert \mathbf{W}\odot (\mathbf{X}-\mathbf{A}\mathbf{S}){\Vert }_F^2}}+\lambda \underset{\text{Classification Error}}{\underbrace{\Vert \mathbf{Z}\odot (\mathbf{Y}-\mathbf{B}\mathbf{S}){\Vert }_F^2}}.$$

$$\begin{array}{rl}\mathbf{X}& \approx {\mathbf{A}}^{(0)}{\mathbf{S}}^{(0)},\\ \mathbf{X}& \approx {\mathbf{A}}^{(0)}{\mathbf{A}}^{(1)}{\mathbf{S}}^{(1)},\\ & ⋮\\ \mathbf{X}& \approx {\mathbf{A}}^{(0)}{\mathbf{A}}^{(1)}\cdots {\mathbf{A}}^{(\mathcal{L})}{\mathbf{S}}^{(\mathcal{L})}.\end{array}$$

$$\Vert \mathbf{X}-{\mathbf{A}}^{(0)}{\mathbf{A}}^{(1)}\cdots {\mathbf{A}}^{(\mathcal{L})}{\mathbf{S}}^{(\mathcal{L})}{\Vert }_F$$

$${\mathbf{S}}^{(\ell )}$$

$$q(\mathbf{A},\mathbf{X})=\underset{\mathbf{S}\ge 0}{\mathrm{argmin}}\Vert \mathbf{X}-\mathbf{A}\mathbf{S}\Vert .$$

$${\mathbf{S}}^{(0)}$$

$${\mathbf{S}}^{(\ell )}$$

$${\left(\frac{\partial q}{\partial \mathbf{x}}(\mathbf{A},\mathbf{x})\right)}_{T,:}={\mathbf{A}}_{:,T}^{†},\text{ and }{\left(\frac{\partial q}{\partial \mathbf{x}}(\mathbf{A},\mathbf{x})\right)}_{T^c,:}=0.$$

$${\left(\frac{\partial q}{\partial {\mathbf{A}}_{i,:}}(\mathbf{A},\mathbf{x})\right)}_{T,T}=-{\left({\mathbf{A}}_{:,T}^{†}\right)}_{:,i}{\left({\mathbf{A}}_{:,T}^{†}\mathbf{x}\right)}^{\top }+{\left(\left(\mathbf{I}-{\mathbf{A}}_{:,T}{\mathbf{A}}_{:,T}^{†}\right)\mathbf{x}\right)}_i{\mathbf{A}}_{:,T}^{†}{\left({\mathbf{A}}_{:,T}^{†}\right)}^{\top },\text{ and }{\left(\frac{\partial q}{\partial {\mathbf{A}}_{i,:}}(\mathbf{A},\mathbf{x})\right)}_{(T\times T{)}^c}=0.$$

$${\mathbf{S}}^{(0)}=q({\mathbf{A}}^{(0)},\mathbf{X}),\text{ and }{\mathbf{S}}^{(\ell )}=q({\mathbf{A}}^{(\ell )},{\mathbf{S}}^{(\ell -1)}),$$

$${\mathbf{\Phi }}^{({\ell }_1,{\ell }_2),m}={\left({{\mathbf{A}}_{:,T_m^{({\ell }_2)}}^{({\ell }_2)}}^{†}\right)}_{:,T_m^{({\ell }_2-1)}}{\left({{\mathbf{A}}_{:,T_m^{({\ell }_2-1)}}^{({\ell }_2-1)}}^{†}\right)}_{:,T_m^{({\ell }_2-2)}}\dots {\left({{\mathbf{A}}_{:,T_m^{({\ell }_1+1)}}^{({\ell }_1+1)}}^{†}\right)}_{:,T_m^{({\ell }_1)}}\left({{\mathbf{A}}_{:,T_m^{({\ell }_1)}}^{({\ell }_1)}}^{†}\right).$$

$${\mathbf{d}}^{({\ell }_1,{\ell }_2),m}={\left({\mathbf{\Phi }}^{({\ell }_1,{\ell }_2),m}\right)}^{\top }{\left(\frac{\partial L}{\partial {\mathbf{S}}^{({\ell }_2)}}\right)}_{T_m^{({\ell }_2)},m}^{\text{*}},$$

$${\mathbf{U}}_{:,T_m^{({\ell }_1)}}^{({\ell }_1,{\ell }_2),m}=-{\mathbf{d}}^{({\ell }_1,{\ell }_2),m}{\left({\mathbf{S}}_{T_m^{({\ell }_1)},m}^{({\ell }_1)}\right)}^{\top }+{\left({\mathbf{S}}^{({\ell }_1-1)}-{\mathbf{A}}^{({\ell }_1)}{\mathbf{S}}^{({\ell }_1)}\right)}_{:,m}{\left({\mathbf{d}}^{({\ell }_1,{\ell }_2,),m}\right)}^{\top }{\left({{\mathbf{A}}_{:,T_m^{({\ell }_1)}}^{({\ell }_1)}}^{†}\right)}^{\top },\text{ and }{\mathbf{U}}_{:,{T_m^{({\ell }_1)}}^c}^{({\ell }_1,{\ell }_2),m}=0.$$

$${\mathbf{S}}^{(0)}=q({\mathbf{A}}^{(0)},\mathbf{X}),\text{ and }{\mathbf{S}}^{(\ell )}=q({\mathbf{A}}^{(\ell )},{\mathbf{S}}^{(\ell -1)}),$$

$$\frac{\partial L}{\partial {\mathbf{A}}^{({\ell }_1)}}={\left(\frac{\partial L}{\partial {\mathbf{A}}^{({\ell }_1)}}\right)}^{\mathbf{S}}+\sum _{\begin{array}{c}{\ell }_1\le {\ell }_2\le \mathcal{L}\\ 1\le m\le M\end{array}}{\mathbf{U}}^{({\ell }_1,{\ell }_2),m}.$$

$$q_{:,m}(\mathbf{A},\mathbf{X})=q(\mathbf{A},{\mathbf{X}}_{:,m}).$$

$$q_T(\mathbf{A},\mathbf{x})={\mathbf{A}}_{:,T}^{†}\mathbf{x},q_{T^c}(\mathbf{A},\mathbf{x})=0.$$

$$\Vert \mathbf{x}-\mathbf{A}\mathbf{s}\Vert \ge \Vert \mathbf{x}-\mathbf{A}q(\mathbf{A},\mathbf{x})\Vert .$$

$$\Vert \mathbf{x}-\mathbf{A}q(\mathbf{A},\mathbf{x})\Vert =\Vert \mathbf{x}-{\mathbf{A}}_{:,T}{q}_T(\mathbf{A},\mathbf{x})-{\mathbf{A}}_{:,T^c}{q}_{T^c}(\mathbf{A},\mathbf{x})\Vert =\Vert \mathbf{x}-{\mathbf{A}}_{:,T}{q}_T(\mathbf{A},\mathbf{x})\Vert ,$$

$$\Vert \mathbf{x}-\mathbf{A}\mathbf{s}\Vert \ge \Vert \mathbf{x}-{\mathbf{A}}_{:,T}{q}_T(\mathbf{A},\mathbf{x})\Vert$$

$$\Vert \mathbf{x}-{\mathbf{A}}_{:,T}\mathbf{w}\Vert =\Vert \mathbf{x}-({\mathbf{A}}_{:,T}{\mathbf{v}}_T+{\mathbf{A}}_{:,T^c}{\mathbf{v}}_{T^c})\Vert =\Vert \mathbf{x}-\mathbf{A}\mathbf{v}\Vert \ge \Vert \mathbf{x}-{\mathbf{A}}_{:,T}{q}_T(A,x)\Vert .$$

$$q_T(\mathbf{A},\mathbf{x})\in \underset{\mathbf{w}\in (0,\infty {)}^{|T|}}{\mathrm{argmin}}\Vert \mathbf{x}-{\mathbf{A}}_{:,T}\mathbf{w}\Vert .$$

$$\mathbf{w}\mapsto \Vert \mathbf{x}-{\mathbf{A}}_{:,T}\mathbf{w}\Vert .$$

$$q_T(\mathbf{A},\mathbf{x})={\mathbf{A}}_{:,T}^{†}\mathbf{x}.$$

$${\mathbf{A}}^{\prime }\mapsto \Vert {{\mathbf{A}}_{:,T}^{\prime }}^{†}\Vert$$

$$\Vert {{\mathbf{A}}_{:,T}^T}^{†}\Vert \ge \Vert {{\mathbf{A}}_{:,T}^{\prime }}^{†}\Vert$$

$$J_1=\underset{\begin{array}{c}T\subset [s]\\ T\ne \varnothing \end{array}}{\max }\Vert {{\mathbf{A}}_{:,T}^T}^{†}\Vert .$$

Full text

DN, JH, ES, JV, and DM are grateful to and were partially supported by NSF CAREER DMS #1348721 and NSF BIGDATA #1740325. This work is based upon work completed at the UCLA CAM REU during Summer 2018 which was funded by NSF DMS #1659676. JH is additionally grateful to and was partially supported by NSF DMS #2211318.

1. Introduction

As the size of available data continues to grow, scalable approaches for extracting meaningful latent trends within large-scale data, or reducing redundant information within the data have become active and important focuses of research. Within this flourishing area of research, topic modeling approaches have received particular interest; topic modeling is a popular class of machine learning techniques that cluster and classify observations to reveal latent themes in a dataset. Algorithms for topic modeling often find application in the domain of document classification and clustering [5, 16, 38, 41, 3, 35], but more recently have found use in applications such as image processing  [20, 24, 29, 12], financial data mining [10], audio processing [7, 17], genetics [32], and bioinformatics [33]. An approach related to, but distinct from, topic modeling is feature extraction. While topic modeling seeks to abstract the dataset and represent data points by these topics, feature extraction aims to find a few key, representative features that best represent the data set for the task at hand (e.g., classification) [21]. Many of the most popular models for these tasks are built upon, or are simply, dimension-reduction techniques which aim to reduce the dimension of the representation of the data; see e.g., [4].

Nonnegative matrix factorization (NMF) is a popular method in machine learning because it is able to both perform feature extraction and generate topic models [6, 9]. Users of NMF specify the number of topics believed to be in the dataset; the model then identifies representative dictionary elements (topics) and coefficients which represent each element of the data set in terms of the topics (thereby lowering the representative dimension of the data). Users often explore topics at different resolutions (number of topics) and seek relationships between the topics learned at various levels of granularity. NMF, however, does not inherently discover the relationships between topics learned at these different levels. A popular alternative, known as hierarchical NMF (HNMF), is to sequentially apply NMF to learn the relationship between NMF topics at different levels of granularity.

Borrowing techniques from neural networks, we seek to modify HNMF to

illustrate the relationships between topics learned at differing levels of granularity, and to specifically provide a hierarchical representation of how topics at finer granularity relate to topics at courser granularity, while avoiding the often high approximation error of naive application of HNMF. An advantage of methods which illustrate hierarchical relationships among the topics over classical NMF is that practitioners can examine the results for multiple numbers of topics without running NMF multiple times. Additionally, while the recovery error increases as the number of layers in a hierarchical model increases, these models have the desirable property that they immediately illustrate the hierarchical structure of the latent topics. In applications in which the learned sub-topics represent known data clusters, the learned hierarchy could provide unknown cluster-level relationships within the data.

We propose a new hierarchical NMF method, which we denote Neural NMF, and illustrate its promise on a synthetic hierarchical data set, the 20 Newsgroups dataset, and MyLymeData, a real dataset containing survey data of Lyme disease patients collected by Lymedisease.org.111An extended abstract of this work appeared in Proc. Interational Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP) 2019 [15]. On the synthtic data set, Neural NMF outerperforms other hierarchical NMF methods in both reconstruction error and classification accuracy. On the 20 Newsgroups dataset, Neural NMF outperforms other methods in classification accuracy and offers a better hierarchical structure and interpretability of the topics. We also illustrate the ability of Neural NMF to identify meaningful, and even surprising, hierarchical topic structure on real world data with the MyLymeData set. These initial results indicate that the Neural NMF method can vastly improve both the reconstruction of the overall hierarchical NMF model as well as improving the learned inter-layer hierarchical structure and interpretability of topics at each layer.

1.1. Organization

The remainder of our paper is organized as follows. In the remainder of the introduction, we briefly introduce notations and conventions in Section 1.2; introduce the foundational models which inspired our own: NMF in Section 1.3, semi-supervised NMF in Section 1.4, hierarchical NMF in Section 1.5, and Deep NMF in Section 1.6; and finally review briefly some further related work in Section 1.7. In Section 2, we introduce our approach which consists of a forward-propogation process and a backward-propogation process (detailed in Section 2.1; here we also include the statements of our main theoretical results that derive the necessary gradient information for backpropagation). In Section 3, we empirically test Neural NMF on a small, synthetic dataset, the 20 Newsgroups dataset, and the MyLymeData dataset. Finally, in Section 4 we offer some conclusions and future directions. The proofs of our main theoretical results are left to Appendix References so as to not distract from the main focus of the paper: the promise of Neural NMF in applications.

1.2. Notation

We distinguish matrices and vectors from scalar quantities using bold font. For a matrix ${\mathbf{F}}$, the notations ${\mathbf{F}}_{i,:}$ and ${\mathbf{F}}_{:,j}$ denote row $i$ and column $j$, respectively. For sets of indices $T$ and $S$, we take ${\mathbf{F}}_{T,:}$ and ${\mathbf{F}}_{:,S}$ to mean the matrix obtained by restricting to the rows of ${\mathbf{F}}$ in $T$ or the columns of ${\mathbf{F}}$ in $S$, respectively. By extension, ${\mathbf{v}}_{T}$ is the vector ${\mathbf{v}}$ restricted to the entries with indices in T. We denote the Moore-Penrose pseudoinverse of ${\mathbf{F}}$ as ${\mathbf{F}}^{\dagger}$. Entrywise (Hadamard) multiplication and division between ${\mathbf{F}}$ and ${\mathbf{G}}$ are denoted by ${\mathbf{F}}\odot{\mathbf{G}}$ and $\frac{{\mathbf{F}}}{{\mathbf{G}}}$, respectively. The vector of length $k$ with all entries equal to one is denoted ${\mathbf{1}}_{k}$, while ${\mathbf{I}}$ indicates the identity matrix of compatible dimensions in all circumstances. We perform subscript (indicial) operations before superscript (e.g., pseudoinversion, transposition) operations whenever applicable. The set $[0,\infty)^{k}$ is denoted $\mathbb{R}_{+}^{k}$. In methods with ${\mathcal{L}}$ layers, we use ${\mathbf{F}}^{(\ell)}$ to mean the value of matrix ${\mathbf{F}}$ at layer $\ell$. We similarly use $k^{(\ell)}$ as the number of topics at layer $\ell$. In any supervised setting, we use $P$ as the number of total classes. Finally, we denote the set of integers $\{1,2,\dots,m\}$ as $[m]$. We let $\|\cdot\|$ denote the Frobenius norm ($\ell_{2}$ norm in the case of vector input) throughout, unless otherwise noted.

1.3. NMF

For a given data matrix ${\mathbf{X}}\in\mathbb{R}_{+}^{N\times M}$ and model rank $k$, nonnegative matrix factorization (NMF) seeks ${\mathbf{A}}\in\mathbb{R}_{+}^{N\times k}$ and ${\mathbf{S}}\in\mathbb{R}_{+}^{k\times M}$ such that ${\mathbf{X}}\approx{\mathbf{A}}{\mathbf{S}}$. To find ${\mathbf{A}}$ and ${\mathbf{S}}$, we wish to solve the optimization problem

[TABLE]

The nonnegativity restriction on ${\mathbf{X}}$, ${\mathbf{A}}$, and ${\mathbf{S}}$ differentiates NMF from other topic modeling and feature extraction techniques, including principal component analysis (PCA) and autoencoding. Since all values are nonnegative, only additive combinations are allowed to approximate every data point in ${\mathbf{X}}$. This allows for a natural and intuitive ‘parts-based’ representation [29]; there are no topics that contribute in a negative manner. We can view the columns of the ${\mathbf{A}}$ matrix as vectors of important features of ${\mathbf{X}}$, or from a topic modeling perspective, as $k$ hidden themes in our data. In this view, the ${\mathbf{S}}$ matrix provides the coefficients to represent each data point of ${\mathbf{X}}$ as a linear combination of the hidden themes (the columns of ${\mathbf{A}}$). Moreover, we can view the product ${\mathbf{A}}{\mathbf{S}}$ as a low-rank approximation of ${\mathbf{X}}$, since

[TABLE]

The problem (1.1) is convex in ${\mathbf{A}}$ and ${\mathbf{S}}$ separately (holding the other fixed), but nonconvex in both. We therefore cannot expect to consistently find a global minimum as a solution. Several techniques exist for finding local minima; most schemes alternate fixing ${\mathbf{A}}$ and ${\mathbf{S}}$ while updating the other variable in order to iteratively decrease (1.1). One of the most popular of these methods is the multiplicative updates methods [30]. Other algorithms for approximating the solution to (1.1) include the alternating least squares and rank-one residue iteration (RRI) methods [23]. Although these offer potential advantages over a multiplicative updates scheme, the latter can be readily extended to deal with more complex cost functions, such as that of semisupervised NMF, which is discussed in the next section.

1.4. Semisupervised NMF (SSNMF)

A natural extension of NMF is to take advantage of any known label information in the factorization [31]. Suppose ${\mathbf{Y}}\in\mathbb{R}^{P\times M}$ is a matrix containing label information for $M$ objects in $P$ classes. Let ${\mathbf{W}}\in\mathbb{R}^{N\times M}$ be the binary indicator matrix for the data, that is

[TABLE]

and let ${\mathbf{Z}}\in\mathbb{R}^{P\times M}$ be the binary indicator matrix for the labels, with

[TABLE]

We can incorporate label information and manage missing data by adjusting our problem to

[TABLE]

The resulting ${\mathbf{B}}\in\mathbb{R}^{P\times k}$ is a linear classification matrix, with ${\mathbf{B}}_{i,:}$ indicating the association between class $i$ and each of the learned topics. The error (1.3) balances the reconstruction error on the data that are known with the classification error on the labels that are known. The relative importance of the classification error is controlled by the user-defined hyperparameter $\lambda$. The multiplicative update minimization strategy discussed in Section 1.3 can be extended to handle this setting [31]. Generalizations of this SSNMF model have been recently proposed in [2].

1.5. Hierarchical NMF (HNMF)

A further extension of both NMF and semisupervised NMF seeks to illuminate hierarchical structure by recursively factorizing the ${\mathbf{S}}$ matrices. By performing NMF with $k=k^{(0)}$, we reveal $k^{(0)}$ hidden topics in the data; by repeating the factorization on the ${\mathbf{S}}$ matrix with $k=k^{(1)}$, we further collect the data into $k^{(1)}$ supertopics among the $k^{(0)}$ topics. This process for ${\mathcal{L}}$ layers approximately factors the data matrix as

[TABLE]

The ${\mathbf{A}}^{(i)}$ matrix represents the how the subtopics at layer $i$ collect into the supertopics at layer $i+1$.

Note that as ${\mathcal{L}}$ increases, the error

[TABLE]

necessarily increases as error propagates with each step. As a result, significant error is introduced when ${\mathcal{L}}$ is large. Choosing $k^{(0)},k^{(1)},...k^{({\mathcal{L}})}$ in practice proves difficult when the number of topics at each layer is unknown, as the number of possibilities grow combinatorially. Additionally, large differences between the number of topics for adjacent layers introduces large error into the factorization. To alleviate the error propagation between layers, HNMF can be endowed with the structure of a neural network. The next sections focus on developing those extensions, inspired by [13].

Finally, as the NMF problem is ill-posed and has an infinite number of global minima (e.g., one may appropriately rescale each of the factors), this ill-posedness is exacerbated in HNMF. There are unicity results for NMF when the matrix ${\mathbf{X}}$ satisfies specific constraints (see e.g., [11, 25, 18, 14, 27]), but we do not know of such results for the HNMF model; this is an important theoretical direction for future investigation.

1.6. Deep NMF (DNMF)

In [13], the authors make a first step toward bridging the gap between HNMF and neural networks. Their method achieves a considerable performance improvement over standard NMF in classification. The forward process for DNMF is HNMF with pooling operator, $p$, applied after each layer of decomposition to introduce nonlinearity and minimize overfitting. Without the pooling operation, the DNMF model is identical to HNMF. The other major contribution of [13] is a proposed backpropagation algorithm meant to refine the result obtained from the forward process. However, the backpropagation technique introduced in [13] differs from backpropagation techniques in neural network settings, as it only propagates one layer at a time and uses multiplicative updates instead of gradient descent to update the values of ${\mathbf{A}}$ and ${\mathbf{S}}$.

1.7. Other Related Work

Similar ideas were explored in [40], [28], and [39]. In [40], the authors develop a hierarchical model in which some of the nonnegativity constraints are relaxed; however, this lacks our proposed backpropagation algorithm for training the model. In [28], the authors propose a NMF backpropagation algorithm using an “unfolding" approach; however, their method does not allow for hierarchy. Finally, a method similar to ours was developed in [39], but differs in that it lacks the nonnegativity constraints that makes our method applicable in topic modeling and feature extraction.

2. Proposed Method: Neural NMF

The proposed backpropagation of DNMF [13] behaves differently than the backpropagation technique typically employed by neural networks. Traditional backpropagation determines the gradient of a cost function with respect to all the weights in the network, so that all the weights in the network may be updated at once. In [13], the update for ${\mathbf{S}}^{(\ell)}$ only depends on ${\mathbf{A}}^{(\ell+1)}$ and ${\mathbf{S}}^{(\ell+1)}$. One of the barriers to formulating a proper backpropagation step for DNMF is that in optimization methods like multiplicative updates [30] or alternating least squares [34], the ${\mathbf{A}}$ and ${\mathbf{S}}$ matrices take turns acting as the independent and dependent variables in the updates. This is in contrast to neural networks, where the weights connecting neurons between layers are always the independent variables, while the activations of the neurons are dependent on the weights. This separation of independent and dependent variables allows for calculation of the derivatives of the dependent variables with respect to the independent variables in a relatively simple way.

To proceed, we choose to regard the ${\mathbf{A}}$ matrices as the independent variables in our model. This is natural since the ${\mathbf{S}}$ matrix at one layer is “passed” to the next layer for factorization, analogous to the neurons’ activations in a neural network being passed to the next layer. In this analogy, the entries of the ${\mathbf{A}}$ matrix now act as the weights between the neurons.

Since we choose to regard the ${\mathbf{A}}$ matrices as the independent variables, we need to determine the ${\mathbf{S}}$ matrices from the ${\mathbf{A}}$ matrices. The natural way to do this is to always require the ${\mathbf{S}}$ matrices to solve the nonnegative least squares problem (2.1). Suppose ${\mathbf{A}}^{(0)},\dots,{\mathbf{A}}^{({\mathcal{L}})}$ are given and define ${\mathbf{S}}^{(-1)}={\mathbf{X}}$. Then we let

[TABLE]

We define the forward-propagation function $q({\mathbf{A}},{\mathbf{X}})$, for any nonnegative matrices ${\mathbf{X}}$ and ${\mathbf{A}}$ with the same number of rows, as the solution to the nonnegative least-squares problem between ${\mathbf{X}}$ and ${\mathbf{A}}$; that is,

[TABLE]

We note that this problem is ill-posed if ${\mathbf{A}}$ does not have full column rank. In what follows, we make the assumption that every ${\mathbf{A}}$ matrix has full column rank. This assumption is reasonable, as the ${\mathbf{A}}$ matrices are always tall matrices (i.e., have more rows than columns), and this condition would only be violated in pathological circumstances. Using this notation, we can define the ${\mathbf{S}}$ matrices as

[TABLE]

These equations show that ${\mathbf{S}}^{(\ell_{j})}$ depends on ${\mathbf{A}}^{(\ell_{i})}$ for $\ell_{i}\leq\ell_{j}$, but not for $\ell_{i}>\ell_{j}$ and they form the forward-propogation stage of Neural NMF. In Algorithms 1 and 2 we display the pseudocode for our proposed method. The partial derivatives $\frac{\partial L}{\partial{\mathbf{A}}^{(i)}}$ will be derived in Theorem 2.6.

In the case that partial label information is provided for the last layer, we can perform a semi-supervised HNMF by additionally calculating the supervision matrix ${\mathbf{B}}$ as ${\mathbf{B}}=({\mathbf{Z}}\odot{\mathbf{Y}})({\mathbf{S}}^{({\mathcal{L}}})^{\dagger}$ where ${\mathbf{Y}}$ and ${\mathbf{Z}}$ are the label and indicator matrices, respectively, defined in Section 1.4. One can include a term in loss function $L$ encouraging matrix similarity between ${\mathbf{Y}}$ and ${\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})}$, such as $\|{\mathbf{Z}}\odot({\mathbf{Y}}-{\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})})\|$, which will influence the learned ${\mathbf{A}}$ and ${\mathbf{S}}$ matrices via backpropagation.

2.1. Backpropagation

In order to derive a backpropagation update for the matrices ${\mathbf{A}}$, we differentiate a cost function which depends on both the ${\mathbf{A}}$ and ${\mathbf{S}}$ matrices, $L=f({\mathbf{X}},{\mathbf{S}}^{(0)},\ldots,{\mathbf{S}}^{({\mathcal{L}})},{\mathbf{A}}^{(0)},\ldots,{\mathbf{A}}^{({\mathcal{L}})})$, with respect to the ${\mathbf{A}}$ matrices. A natural choice for the cost $L$ is the HNMF error (1.5), however our method can be employed for other choices of cost functions. To differentiate, we employ the chain rule, differentiating $q({\mathbf{A}},{\mathbf{X}})$ with respect to both ${\mathbf{A}}$ and ${\mathbf{X}}$. See Figure 1 for a visualization of the computational graph for the forward propagation.

We begin by computing formulas for the derivatives of $q({\mathbf{A}},{\mathbf{X}})$ with respect to ${\mathbf{A}}$ and ${\mathbf{X}}$. Lemma A.1 (presented in the appendix), shows that $\frac{\partial q}{\partial X}({\mathbf{A}},{\mathbf{X}})$ can be computed columnwise (with respect to columns of ${\mathbf{X}}$), so we need only determine derivative formulas for $q({\mathbf{A}},{\mathbf{x}})$, where ${\mathbf{x}}$ is a column of matrix ${\mathbf{X}}$. We define a space in which differentiation of $q({\mathbf{A}},{\mathbf{X}})$ is relatively simple below.

Definition 2.1**.**

Let $U_{r,s}\subset\mathbb{R}_{+}^{r\times s}\times\mathbb{R}_{+}^{r}$ ($r\geq s$) be the set of all matrix-vector pairs $({\mathbf{A}},{\mathbf{x}})$ such that ${\mathbf{A}}$ has full (column) rank and the support of $q$ does not change in some neighborhood of $({\mathbf{A}},{\mathbf{x}})$.

Lemma 2.2**.**

The set defined in Definition 2.1, $U_{r,s}$ is dense, open, and has full (Lebesgue) measure.

Roughly, this lemma says that most pairs ${\mathbf{A}}$ and ${\mathbf{x}}$ encountered in the course of our method will be such that the support of the computed column of $S$ will not change; that is, we can use supports learned during forward propagation for backpropagation gradients. We delay the somewhat technical proof of this lemma to Appendix A. Now we utilize Lemma A.2 to differentiate $q({\mathbf{A}},{\mathbf{x}})$.

Theorem 2.3**.**

Let $U_{r,s}$ be the set defined in Definition 2.1, and suppose $({\mathbf{A}},{\mathbf{x}})\in U_{r,s}$. Let $T=\operatorname*{supp}q({\mathbf{A}},{\mathbf{x}})$. The partial derivative of $q$ with respect to ${\mathbf{x}}$ is given entry-wise as

[TABLE]

The partial derivative of $q$ with respect to row $i$ of ${\mathbf{A}}$ for $i\in[r]$ is given entry-wise as

[TABLE]

These formulas for the derivatives of $q$ are sufficient for implementing a backwards propagation algorithm in a machine learning library, such as PyTorch or TensorFlow. We summarize in Theorem 2.6 how these partial derivatives are applied via the chain rule for an arbitrary differentiable cost function $L$. Assume there are given values for the variables ${\mathbf{A}}^{(0)},\cdots,{\mathbf{A}}^{({\mathcal{L}})}$ and the dependent (matrix) variables ${\mathbf{S}}^{(\ell)}$ are defined as functions of the independent (matrix) variables ${\mathbf{A}}^{(\ell)}$ recursively as

[TABLE]

for $\ell\in[\mathcal{L}]$, where $\mathcal{L}$ is the number of layers.

This calculation requires one to collect derivatives along all paths between $L$ and the variable of interest in the computational graph in Figure 1. We achieve this via the use of several auxiliary variables which collect partial derivatives along different segments of the paths. First, the variable $\Phi$ collects (column-wise) derivatives along the central path between ${\mathbf{S}}^{(\ell_{2})}$ and ${\mathbf{S}}^{(\ell_{1})}$.

Definition 2.4**.**

Let $T_{m}^{(\ell)}=\operatorname*{supp}{\mathbf{S}}^{(\ell)}_{:,m}$ where $m\in[k^{(\ell-1)}]$ is a column index of ${\mathbf{S}}^{(\ell)}$ which has $k^{(\ell-1)}$ columns. Define ${\mathbf{\Phi}}^{(\ell_{1},\ell_{2}),m}$ for ${\mathcal{L}}\geq\ell_{2}\geq\ell_{1}\geq 0$ by

[TABLE]

The variable $\mathbf{d}$ utilizes $\Phi$ to collect (column-wise) derivatives along the path that follows the edge from $L$ back to ${\mathbf{S}}^{(\ell_{2})}$, and along the central path between ${\mathbf{S}}^{(\ell_{2})}$ to ${\mathbf{S}}^{(\ell_{1})}$. Finally, the variable $\mathbf{U}$ adds to this path the edge between ${\mathbf{S}}^{(\ell_{1})}$ and ${\mathbf{A}}^{(\ell_{1})}$, which is the variable of interest.

Definition 2.5**.**

Suppose $L$ is a cost function depending on all the variables ${\mathbf{S}}^{(\ell)}$ and ${\mathbf{A}}^{(\ell)}$. Let $\left(\frac{\partial L}{\partial{\mathbf{S}}^{(\ell_{2})}}\right)^{\textbf{*}}$ denote the derivative of $L$ with respect to ${\mathbf{S}}^{(\ell_{2})}$ holding ${\mathbf{S}}^{(\ell_{2}+1)}$, $\ldots$ , ${\mathbf{S}}^{({\mathcal{L}})}$ constant. We define

[TABLE]

and ${\mathbf{U}}^{(\ell_{1},\ell_{2}),m}$ with

[TABLE]

The final desired derivative, $\frac{\partial L}{\partial{\mathbf{A}}^{(\ell_{1})}}$, then sums these collected derivatives over all paths between $L$ and ${\mathbf{A}}^{(\ell_{1})}$, including the direct edge.

Theorem 2.6**.**

Let the dependent (matrix) variables ${\mathbf{S}}^{(\ell)}$ be defined as functions of the independent (matrix) variables ${\mathbf{A}}^{(\ell)}$ recursively as

[TABLE]

for $\ell\in[\mathcal{L}]$, where $\mathcal{L}$ is the number of layers. Let $T_{m}^{(\ell)}$ and ${\mathbf{\Phi}}^{(\ell_{1},\ell_{2}),m}$ for ${\mathcal{L}}\geq\ell_{2}\geq\ell_{1}\geq 0$ be as defined in Definition 2.4. Fixing a point in the space of ${\mathbf{A}}$-matrices, suppose that for all $0\leq\ell\leq\mathcal{L}$ and $m\in[k^{(\ell-1)}]$, $({\mathbf{A}}^{(\ell)},{\mathbf{S}}^{(\ell-1)}_{:,m})$ is in the set $U_{k^{(\ell-1)},k^{(\ell)}}$ defined in Definition 2.1. Suppose $L$ is a cost function depending on all the variables ${\mathbf{S}}^{(\ell)}$ and ${\mathbf{A}}^{(\ell)}$. Let $\left(\frac{\partial L}{\partial{\mathbf{S}}^{(\ell_{2})}}\right)^{\textbf{*}}$, ${\mathbf{d}}^{(\ell_{1},\ell_{2}),m}$, and ${\mathbf{U}}^{(\ell_{1},\ell_{2}),m}$ be as defined in Definition 2.5. Then, if we let $\left(\frac{\partial L}{\partial{\mathbf{A}}^{(\ell_{1})}}\right)^{\mathbf{S}}$ be the derivative of $L$ with respect to ${\mathbf{A}}^{(\ell_{1})}$, holding the ${\mathbf{S}}$ matrices constant, we have

[TABLE]

That is, if the matrices ${\mathbf{A}}$ and ${\mathbf{S}}$ are such that $({\mathbf{A}}^{(\ell)},{\mathbf{S}}^{(\ell-1)}_{:,m})\in U_{k^{(\ell-1)},k^{(\ell)}}$ of Definition 2.1 then we can compute the derivative of differentiable cost function $L$ with respect to the independent variables ${\mathbf{A}}$ using simple matrix operations and the support information calculated during forward propagation.

3. Experimental Results

We test Neural NMF on three datasets: a small, synthetic dataset, the 20 Newsgroups dataset, and the MyLymeData dataset. The synthetic dataset is a small block matrix, with three different levels of hierarchy in the blocks. The 20 Newsgroups dataset is a common benchmark dataset in which hierarchy between the topics of the dataset is known. Meanwhile, the MyLymeData dataset represents the symptoms experienced by a group of Lyme disease patients, and hierarchy in this dataset is not known a priori. Our implementation is available in the indexed Python package NeuralNMF and the code for experiments is provided on Github [22].

On the synthetic dataset, we compare Neural NMF to HNMF and Deep NMF in the unsupervised, semi-supervised, and fully supervised settings with 1, 2, and 3 layers, and report classification accuracy and reconstruction loss. We see that in each of these settings, Neural NMF outperforms HNMF and Deep NMF in both classification and reconstruction, and forms a better low rank representation of the data set that appears to preserve more of the coarser block structure.

On the 20 Newsgroups dataset, we compare Neural NMF to HNMF in the unsupervised and semi-supervised settings for 2 layer experiments, and report classification accuracy at each layer. We focus on classification and qualitative analysis rather than reconstruction because on this data set, one cannot ask to produce a highly accurate low-rank reconstruction of the data but instead seek to form meaningful and class-discriminatory topics that form a hierarchical structure. Our experimental results show that for both the unsupervised and supervised settings, Neural NMF attains a higher classification accuracy than HNMF at each layer, and the topics formed by Neural NMF have significantly better interpretability and hierarchical structure. We also see that despite only being provided partial label information at the second layer, Neural NMF is able to improve the classification accuracy at the first layer with the additional of this label information, demonstrating that unlike HNMF, Neural NMF is able to propagate information provided at the last layer to earlier layers.

On the MyLymeData set, we see the potential for Neural NMF to produce hierarchical topic structure on a real-world large-scale survey dataset. While the ground truth hierarchical structure is unknown for this real, messy dataset, we note that the results produced by Neural NMF yield interpretable results that reflect both what is well-known and unknown about Lyme disease patients and the manifestation of their symptoms.

3.1. Synthetic Data

We first test the reconstruction and classification ability of Neural NMF in an idealized setting: on a $90\times 87$ noisy toy dataset with a clear three-layer hierarchical structure. Starting with two large blocks, we overlay increasingly smaller and more intense asymmetric regions along the diagonal of a matrix, and finally add a $\text{uniform}(0,1)$ noise to the entire matrix; see the left plot of Figure 2. We know the optimal model rank sequence to be $k^{(0)}=9,k^{(1)}=4,k^{(2)}=2$ a priori. We test HNMF, Deep NMF, and Neural NMF with one, two, or three layer structure, and various levels of supervision. The labels representing to which of the nine classes the data points belong (grouped by the highest intensity blocks) are given for 40% of the data (semisupervised) or 100% of the data (supervised). For each level of supervision and depth, the results are averaged over 25 trials. We present the recovery error and classification accuracy measuring the discrepancy between ${\mathbf{Y}}$ and the computed matrices ${\mathbf{B}}$ and ${\mathbf{S}}$ for these experiments in Table 1. The reconstruction error is computed relative to the norm of the original matrix as $\|{\mathbf{X}}-{\mathbf{A}}^{(0)}{\mathbf{A}}^{(1)}\cdots{\mathbf{A}}^{({\mathcal{L}})}{\mathbf{S}}^{({\mathcal{L}})}\|/\|{\mathbf{X}}\|$. Object $m$ is predicted to have label $p$ if $({\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})})_{pm}=\max$ $({\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})})_{:m}$, and classification accuracy is the proportion of predicted labels that match the true labels. Noteworthy improvements of Neural NMF over HNMF and Deep NMF are bolded. We comment that we expect the advantage Neural NMF enjoys over HNMF and Deep NMF is due to the backpropagation method allowing it to avoid suboptimal local minima found by HNMF and Deep NMF. Additionally, in this example the classification accuracy is inherently tied to the reconstruction error (as the labels are generated by entries of the matrix), so Neural NMF is able to achieve good accuracy even in the unsupervised setting.

While the approximations produced at the $\ell$th layer have rank $k^{(\ell)}$, the rank of the final approximation produced by the hierarchical model will be $k^{(\mathcal{L})}$. For this reason, we can only seek the $k^{(\mathcal{L})}$ most representative features (blocks) when qualitatively evaluating the final approximations produced by hierarchical models.

In Figure 2, we visualize the reconstructions produced by each method with no supervision and two layer structure ($k^{(0)}=9$ and $k^{(1)}=4$). We cannot hope to resolve the highest-intensity features from the original data in our reconstructions, as the NMF approximations have lower than necessary rank. Instead, we consider how accurately the methods reconstruct the two-layer block structure. Although each method is able to capture some of the structure, it is clear that Neural NMF outperforms HNMF and Deep NMF, resolving sharper blocks.

In Figure 3, we see that Neural NMF similarly outperforms HNMF and Deep NMF in a semisupervised three-layer trial ($k^{(0)}=9,k^{(1)}=4,k^{(2)}=2$); note that HNMF and Deep NMF produce many columns with extemely low intensity (HNMF entirely misses the second of the two coarsest blocks), while Neural NMF produces a rank-2 approximation which correctly reconstructs the coarsest two-block structure of least intensity (lightest blue). We note that adding label information is expected to lead to worse reconstruction because the optimization task will focus on improving classification, often at the detriment of higher reconstruction loss.

3.2. 20 Newsgroups Data

The 20 Newsgroups dataset is a collection of approximately 20,000 text documents containing the text of messages from 20 different newsgroups on the distributed discussion system Usenet which functioned similarly to current internet discussion forums. The documents are partitioned nearly evenly across the 20 newsgroups which can be further classified into six supergroups (computers, for sale, sports/recreation, politics, science, religion) [26]. This clear hierarchical topic structure makes this an appropriate testing ground for Neural NMF.

In our experiments, we subsample 100 documents from each of 10 subtopics from the 20 Newsgroups dataset. We encode the text data in a word frequency vector representation of the bag-of-words model. We perform experiments for both the semi-supervised and unsupervised tasks. For semi-supervision, we provide labels for 75% of the documents and compute classification accuracy on the 25% of documents without provided labels. We present the keywords from each experiment (those words which are represented with largest magnitude in each topic) and the classification accuracy measuring the discrepancy between ${\mathbf{Y}}$ and the computed matrices ${\mathbf{B}}$ and ${\mathbf{S}}$. The reconstruction error is computed as $\|{\mathbf{X}}-{\mathbf{A}}^{(0)}{\mathbf{A}}^{(1)}\cdots{\mathbf{A}}^{({\mathcal{L}})}{\mathbf{S}}^{({\mathcal{L}})}\|/\|{\mathbf{X}}\|$. Object $m$ is predicted to have label $p$ if $({\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})})_{pm}=\max$ $({\mathbf{B}}{\mathbf{S}}^{({\mathcal{L}})})_{:m}$, and classification accuracy is the proportion of predicted labels that match the true labels.

In Table 2, we display the classifications accuracies for Neural NMF, Deep NMF, and HNMF for the first and second layers of an unsupervised and semi-supervised 2 layer experiment, where supervision labels are provided only for the last layer. Each experiments was run for ten trials and we report the average of the trials. We see that Neural NMF outperforms HMF in all setting, and outerperforms Deep NMF when supervision information is provided. Comparing the unsupervised and supervised experiments, we see that even though label information is provided only at the second layer, the first layer of Neural NMF gains substantial improvement when supervision is added, suggesting that the classification information at the second layer successfully propagated to the first layer. This is not possible for HNMF, where each layer’s factorization is computed separately. We also see that for Deep NMF the first layer does not have significant improvement when supervision information is added.

In Tables 3 and 4, we display the keywords learned by Neural NMF at the first and second layers, respectively, for the semi-supervised two layer experiment. We see that the words are meaningful and representative of the topics within the 20 Newsgroups dataset, and we note that at rank 6, each of the six topics related directly to one of the six super topics of this dataset. We are also able to see clear hierarchical structure, such as a medical topic and space topic at the first layer (topics 5 and 7 in Table 3) that combine into a science topic in the second layer (topic 4 in Table 4). This hierarchical relationship is also evident from Figure 4, where we display a heat map of the $A_{2}$ matrix for the 2 layer semi-supervised Neural NMF experiment. We see a clear relationship between topics at rank 10 and rank 6, which agrees with the known hierarchy. The topic labels at rank 10 were chosen qualitatively based on the keywords seen in Tables 3, and the topic labels at rank 6 were determined by the presence of 20 Newsgroups data set classes in the rows of the ${\mathbf{S}}^{(1)}$ matrix.222At first glance, the keywords in Topic 6 of layer 1 do not appear to correspond to any 20 Newsgroups document topic, but all 10 keywords come from the email signature of a prolific user within medicine newsgroups with the name Gordon Banks, who includes his radio call sign N3JXP, his email geb@cadre.dsl.pitt.edu, and the quote “Skepticism is the chastity of the intellect, and it is shameful to surrender it too soon” (see e.g., [8, Pg. 259]).

In Tables 5 and 6 we display the keywords learned by HNMF at the first and second layers, respectively, for the semi-supervised two-layer experiment, and in Figure 4 we display a heat map of the ${\mathbf{A}}_{2}$ matrix for this experiment which shows the relationship between the topics at each layer. We see that while most of the topics provide salient topic modelling information corresponding to the 20 Newsgroups topics, some topics are fairly unclear (e.g., Topic 6, 7 of rank 10) and the learned hierarchical structure does not adhere well to the expected structure (e.g., baseball and motorcycles do not collect into the recreation super-topic, medicine and space do not collect into the science super-topic).

3.3. Lyme Data

We conclude with a real world example using Lyme disease data. The MyLymeData dataset [1] used in this analysis consists of survey responses of approximately 4000 current and former Lyme disease patients (which has since grown to over 17,000). The questions cover demographics, symptoms at various stages of the disease, medical procedures, and more. Responses may take binary, categorical, or scalar values. Each patient is self-identified as ‘well’ or ‘unwell’ at the time of the survey. We center our analysis on a subset of the dataset concerning patient symptom information shortly after an initial tick bite and at the time of diagnosis. These questions yield binary data indicating whether or not, for example, a patient observed a ‘bulls-eye rash’ somewhere on his or her body at the time of possible diagnosis. As we do not know the hierarchy of this data a priori, we simply run exploratory experiments in this section. We highlight potential advantages of Neural NMF illustrated in these results.

Figure 6 is an example of the ${\mathbf{S}}$ matrix when NMF is performed with model rank $k=5$ alongside the ${\mathbf{S}}$ matrices at the second layer for HNMF and Neural NMF with full network structure $k^{(0)}=6,k^{(1)}=5$, and $k^{(2)}=4$. First, note in Figure 6 that at the second layer, Neural NMF reveals topics that are extremely similar to those produced by simple NMF. The Neural NMF topic structure reveals in particular an interesting observation about the appearance of a bulls eye rash, previously thought to be a critical element to a Lyme disease diagnosis. However, in the Neural NMF topic structure, this symptom actually plays very little role in any of the topics, versus appearing fairly strongly in the other two structures. Further investigation here is critical, as physicians are now beginning to agree that the prevalence (and thus importance) of the bulls eye rash symptom may be seriously less than previously believed.

To further hone this observation, we present the results of Neural NMF with $k^{(0)}=6$ and $k^{(1)}=5$ on a subset of the data corresponding to symptoms experienced by patients at diagnosis but for unwell and well patients separately. This is shown in Figure 7, where we see drastically different positioning of the bulls eye rash symptom. In well patients, it forms a very strong topic indicating a prevalence in that patient group. However, in unwell patients, it is mildly represented among a topic seeming to indicate general malaise. This warrants further investigation again, as this rash may play a role in whether a patient becomes well or not. A final interesting observation is that twitching, a physical symptom, appears in the same topic as cognitive symptoms such as sleep impairment and memory loss. These are all symptoms of a so-called neurological Lyme disease manifestation that is still not yet understood. In summary, studying an accurate hierarchical topic structure via Neural NMF leads to important and interesting directions of further study both mathematically and medically.

4. Conclusion

We present a novel method for hierarchical multilayer nonnegative matrix factorization which incorporates the backwards propagation technique from deep learning to minimize error accumulation. Preliminary tests on toy datasets show this method outperforms existing multilayer NMF algorithms. The forward and backwards propagation steps of Neural NMF may offer decreased reconstruction and classification error over single-application HNMF. Additionally, it seems that Neural NMF often better resolves data points into a single topic.

Future directions include to further compare Neural NMF and others on various datasets in order to find precise regimes in which it offers substantial improvement. Furthermore, theoretically analyzing the convergence of Neural NMF on ideal datasets (those containing latent hierarchy) is an important future direction of work.

5. Acknowledgements

The authors would like to thank LymeDisease.org for the use of data derived from MyLymeData to conduct our experiments and the patients for their contributions to MyLymeData. They additionally thank Dr. Anna Ma for her instruction on the MyLymeData dataset, LymeDisease.org CEO Lorraine Johnson for her collaboration, Dr. Blake Hunter for proposing the project, and Professor Andrea Bertozzi, Director of Applied Math at UCLA, for organizing the REU program through which this research was conducted.

Appendix A Proofs of Lemma A.3 and Proposition 2.2

First, we state two lemmas about $q$, which will lead to the derivative formula. Our first lemma demonstrates that $q$ actually acts column-wise on ${\mathbf{X}}$.

Lemma A.1**.**

Suppose ${\mathbf{A}}\in\mathbb{R}_{+}^{N\times k}$ and ${\mathbf{X}}\in\mathbb{R}_{+}^{N\times M}$ are nonnegative matrices (with the same number of rows), and ${\mathbf{A}}$ has full column rank. Then, for all $m$,

[TABLE]

This lemma is a consequence of the column-wise separability of the Frobenius norm objective and simplifies the problem of finding the derivatives of $q({\mathbf{A}},{\mathbf{X}})$ for general ${\mathbf{A}}$ and ${\mathbf{X}}$ to the problem of finding the derivatives of $q({\mathbf{A}},{\mathbf{x}})$ where ${\mathbf{x}}$ is a column vector of size $N$. The following lemma now gives us a formula for $q({\mathbf{A}},{\mathbf{x}}),$ which we will be able to differentiate.

Lemma A.2**.**

Suppose ${\mathbf{A}}\in\mathbb{R}_{+}^{r,s}$ has full column rank, and ${\mathbf{x}}\in\mathbb{R}_{+}^{r}$. Let $T$ be the support of $q({\mathbf{A}},{\mathbf{x}})$, i.e., the set of indices on which the vector is strictly positive. Then $q({\mathbf{A}},{\mathbf{x}})$ is given by

[TABLE]

Proof.

Since $q({\mathbf{A}},{\mathbf{x}})=\operatorname*{argmin}_{{\mathbf{s}}\geq 0}\|{\mathbf{x}}-{\mathbf{A}}{\mathbf{s}}\|$, for any nonnegative vector ${\mathbf{s}}$,

[TABLE]

Then, using that $q({\mathbf{A}},{\mathbf{x}})$ has support $T$ (and so is zero outside that support), we get

[TABLE]

so

[TABLE]

for any nonnegative vector ${\mathbf{s}}$. Suppose ${\mathbf{w}}$ is in $(0,\infty)^{|T|}$, i.e., it is a positive vector of length $|T|$. Let ${\mathbf{v}}$ be the vector in $\mathbb{R}_{+}^{r}$ with ${\mathbf{v}}_{T}={\mathbf{w}}$ and ${\mathbf{v}}_{T^{c}}=0$. Then

[TABLE]

This holds for any ${\mathbf{w}}$ in $(0,\infty)^{|T|}$.

Since $q({\mathbf{A}},{\mathbf{x}})$ has support equal to $T$, all of the entries of $q_{T}({\mathbf{A}},{\mathbf{x}})$ are positive, so it is in $(0,\infty)^{|T|}$. So, by what we just showed,

[TABLE]

Thus, since $(0,\infty)^{|T|}$ is an open set, we have that $q_{T}({\mathbf{A}},{\mathbf{x}})$ is a local minimum of the function

[TABLE]

Since ${\mathbf{A}}$ has full column rank, ${\mathbf{A}}_{:,T}$, a matrix whose columns are a subset of those of ${\mathbf{A}}$, must also have full column rank. But then we know that this function has only one local minimum, equal to its global minimum, and is given by ${\mathbf{A}}_{:,T}^{\dagger}{\mathbf{x}}$. Thus, we have that

[TABLE]

Using the fact that, by the definition of support, $q_{T^{c}}({\mathbf{A}},{\mathbf{x}})=0$, we arrive at the desired result. ∎

This result cannot be used to calculate $q({\mathbf{A}},{\mathbf{x}})$ since it requires knowledge a priori of the answer’s support. However, if the support of $q({\mathbf{A}},{\mathbf{x}})$ is locally constant in $({\mathbf{A}},{\mathbf{x}})$-space, then once we have calculated the support of $q({\mathbf{A}},{\mathbf{x}})$ we can use this formula to differentiate $q$ with respect to ${\mathbf{A}}$ and ${\mathbf{x}}$. Lemma 2.2 guarantees that we can do this in almost all circumstances. The following lemma will be used in the proof of Lemma 2.2.

Lemma A.3**.**

Let $V_{r,s}$ be the set of matrices in $\mathbb{R}_{+}^{r\times s}$ ($r\geq s$) with full column rank. Then $q$ is continuous on $V_{r,s}\times\mathbb{R}_{+}^{n}$.

Remark**.**

In the following proof, all matrix norms are taken to be spectral norms, rather than Frobenius norms.

Proof.

Fix a point $({\mathbf{A}},{\mathbf{x}})$ in $V_{r,s}\times\mathbb{R}_{+}^{r}$; we will show $q$ is continuous there. Let $N$ be a bounded neighborhood of ${\mathbf{A}}$ such that $\overline{N}\subset V_{r,s}$. Let $p_{T}:\mathbb{R}^{r\times s}\rightarrow\mathbb{R}^{r\times|T|}$ be the projection which maps ${\mathbf{A}}\mapsto{\mathbf{A}}_{:,T}$ for $T\subset[s]$. As a projection map, $p_{T}$ is clearly continuous, and so $p_{T}\left(\overline{N}\right)$ is compact, since $\overline{N}$ is compact. Since ${\mathbf{A}}^{\prime}_{:,T}$ must have full column rank if ${\mathbf{A}}^{\prime}$ does, and since $\overline{N}\subset V_{r,s}$, we must have that $p_{T}\left(\overline{N}\right)\subset V_{r,|T|}$. Since the pseudoinverse operation is continuous on sets of matrices with constant rank, it is continuous on $p_{T}\left(\overline{N}\right)$. Additionally, the spectral norm is continuous everywhere, so we get that the map $n_{T}:\overline{N}\rightarrow\mathbb{R}$ given by

[TABLE]

is continuous for each nonempty $T\subset[s]$. By the extreme value theorem, for each nonempty $T\subset[s]$, there is a matrix ${\mathbf{A}}^{T}$ in $\overline{N}$ such that

[TABLE]

for all ${\mathbf{A}}^{\prime}$ in $\overline{N}$. Then we let

[TABLE]

Thus, for all ${\mathbf{A}}^{\prime}$ in $\overline{N}$, we must have that $\left\|{{\mathbf{A}}^{\prime}_{:,T}}^{\dagger}\right\|\leq J_{1}$ for any set of indices $T\subset[s]$. In particular, for any $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N\times\mathbb{R}_{+}^{r}$, by Lemma A.2, we have that for $T=\text{supp}(q({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime}))$, $q_{T}({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})={{\mathbf{A}}^{\prime}_{:,T}}^{\dagger}{\mathbf{x}}^{\prime}$ and $q_{T^{c}}({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})=0$, so

[TABLE]

where the last inequality follows from the definition of $J_{1}$. Let $N^{\prime}$ be a bounded neighborhood of ${\mathbf{x}}$ in $\mathbb{R}_{+}^{r}$. Let $J_{2}=\sup_{x^{\prime}\in N^{\prime}}\|x^{\prime}\|$. This finally gives that for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N\times N^{\prime}$,

[TABLE]

Now, suppose $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ is in $N\times N^{\prime}$. By the definition of $q$, we have that

[TABLE]

and

[TABLE]

We claim that the map $\|{\mathbf{x}}-{\mathbf{A}}q(\cdot)\|$ is continuous at $({\mathbf{A}},{\mathbf{x}})$, i.e., for each $\epsilon>0$, there is an open set $W\subset\mathbb{R}_{+}^{r\times s}\times\mathbb{R}_{+}^{r}$ containing $({\mathbf{A}},{\mathbf{x}})$ such that for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $W$ we have

[TABLE]

To see this, fix $\epsilon>0$. We have that for $({\mathbf{A}}^{\prime},\bm{x}^{\prime})$ and $({\mathbf{A}}^{\prime\prime},{\mathbf{x}}^{\prime\prime})$ in $N\times N^{\prime}$, by the reverse triangle inequality,

[TABLE]

In particular, the above holds for $({\mathbf{A}}^{\prime\prime},{\mathbf{x}}^{\prime\prime})=({\mathbf{A}},{\mathbf{x}})$ and $({\mathbf{A}}^{\prime\prime},{\mathbf{x}}^{\prime\prime})=({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$, so

[TABLE]

[TABLE]

From these, we get

[TABLE]

[TABLE]

Now, using equation (A.1), we get

[TABLE]

Now, using equations (A.3) and (A.4), we get

[TABLE]

Finally, applying equation (A.2), we have

[TABLE]

Putting this all together gives

[TABLE]

for any $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N\times N^{\prime}$. Let

[TABLE]

Then let $B_{\delta}({\mathbf{x}})$ be the open $\ell_{2}$ ball centered at ${\mathbf{x}}$ in $\mathbb{R}^{r}$ of radius $\delta$, and $B_{\delta}({\mathbf{A}})$ be the open spectral norm ball centered at ${\mathbf{A}}$ in $\mathbb{R}^{r\times s}$ of radius $\delta$. Then we have that for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $(B_{\delta}({\mathbf{x}})\cap N)\times(B_{\delta}({\mathbf{A}})\cap N^{\prime})$,

[TABLE]

Thus, setting $W=(B_{\delta}({\mathbf{x}})\cap N)\times(B_{\delta}({\mathbf{A}})\cap N^{\prime})$ gives what we wanted.

Now we show that $q$ is continuous at $({\mathbf{A}},{\mathbf{x}})$. We fix $\epsilon>0$ and attempt to find an open set $W$ of $({\mathbf{A}},{\mathbf{x}})$ so that $\|q({\mathbf{A}},{\mathbf{x}})-q({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})\|<\epsilon$ for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $W$. Since $q({\mathbf{A}},{\mathbf{x}})$ is the unique global minimum (since ${\mathbf{A}}$ has full column rank) of the convex function

[TABLE]

on $\mathbb{R}_{+}^{r}$, we can apply [37, Theorem 27.2]. Thus, there is a $\delta>0$ such that, if ${\mathbf{s}}$ is a nonnegative vector, then

[TABLE]

Then since we proved that the function $\|{\mathbf{x}}-{\mathbf{A}}q(\cdot)\|$ is continuous at $({\mathbf{A}},{\mathbf{x}})$, there is a neighborhood $W$ of $({\mathbf{A}},{\mathbf{x}})$ in $\mathbb{R}_{+}^{r\times s}\times\mathbb{R}_{+}^{r}$ such that for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $W$,

[TABLE]

Then clearly for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $W$,

[TABLE]

so by (A.5), we get

[TABLE]

for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $W$, making $W$ our desired neighborhood. Thus, $q$ is continuous at $({\mathbf{A}},{\mathbf{x}})$. Since $({\mathbf{A}},{\mathbf{x}})$ was chosen arbitrarily in $V_{r,s}\times\mathbb{R}_{+}^{r}$, we immediately have that $q$ is continuous on this entire set, as desired.

∎

We may now prove Lemma 2.2.

Proof of Lemma 2.2.

Since in the regular Lebesgue measure, we have that full-measure sets are dense, it suffices to show that $U_{r,s}$ is open and has full-measure. First we show that it is open. We again adopt the notation $V_{r,s}$ to denote the set of matrices in $\mathbb{R}_{+}^{r\times s}$ ($r\geq s$) with full (column) rank. By its definition, since $q$ is only defined on $V_{r,s}\times\mathbb{R}_{+}^{n}$, we have that $({\mathbf{A}},{\mathbf{x}})\in U_{r,s}$ exactly when there is a neighborhood $N$ of $({\mathbf{A}},{\mathbf{x}})$ contained in $V_{r,s}\times\mathbb{R}_{+}^{r}$ such that $\operatorname*{supp}q({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ is constant on $N$. For each $({\mathbf{A}},{\mathbf{x}})\in U_{r,s}$, let $N_{{\mathbf{A}},{\mathbf{x}}}$ be such a neighborhood. Then for all $({\mathbf{A}},{\mathbf{x}})$ in $U_{r,s}$, $N_{{\mathbf{A}},{\mathbf{x}}}\subset U_{r,s}$, since for each $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N_{{\mathbf{A}},{\mathbf{x}}}$, $N_{{\mathbf{A}},{\mathbf{x}}}$ is a neighborhood of $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $V_{r,s}\times\mathbb{R}_{+}^{r}$ on which the support of $q$ is constant, so $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})\in U_{r,s}$. Thus, we must have

[TABLE]

which demonstrates that $U_{r,s}$ is open.

We now show that $U^{c}_{r,s}$ has zero measure. By the definition of $U_{r,s}$, we have

[TABLE]

Since $V_{r,s}^{c}$ has zero measure, so does $V_{r,s}^{c}\times\mathbb{R}_{+}^{r}$, meaning that it suffices to show that

[TABLE]

has zero measure. Suppose that $({\mathbf{A}},{\mathbf{x}})$ is in $W$. Then there is a sequence $({\mathbf{A}}^{k},{\mathbf{x}}^{k})$ in $V_{r,s}\times\mathbb{R}_{+}^{r}$ such that $({\mathbf{A}}^{k},{\mathbf{x}}^{k})\rightarrow({\mathbf{A}},{\mathbf{x}})$ but $\operatorname*{supp}q({\mathbf{A}}^{k},{\mathbf{x}}^{k})\neq\operatorname*{supp}q({\mathbf{A}},{\mathbf{x}})$. Since there are only finitely many possible supports, at least one support must be represented infinitely often in the sequence $\operatorname*{supp}q({\mathbf{A}}^{k},{\mathbf{x}}^{k})$, so, by possibly restricting to a subsequence, we may assume without loss of generality that $\operatorname*{supp}q({\mathbf{A}}^{k},{\mathbf{x}}^{k})$ is constant. Let $T_{1}=\operatorname*{supp}({\mathbf{A}}^{k},{\mathbf{x}}^{k})$ and let $T_{0}=\operatorname*{supp}({\mathbf{A}},{\mathbf{x}})$. By hypothesis, $T_{0}\neq T_{1}$.

We claim that $T_{0}\subset T_{1}$. To see this, suppose it were not true. Then there would be an index $i$ in $T_{0}$ which is not in $T_{1}$. Then, since $i$ is in $T_{0}$, $q_{i}({\mathbf{A}},{\mathbf{x}})>0$. On the other hand, $i$ is not in $T_{1}$, so $q_{i}({\mathbf{A}}^{k},{\mathbf{x}}^{k})=0$. Since $({\mathbf{A}}^{k},{\mathbf{x}}^{k})\rightarrow({\mathbf{A}},{\mathbf{x}})$, and $q$ is continuous by Lemma A.3, we must have that $q_{i}({\mathbf{A}}^{k},{\mathbf{x}}^{k})\rightarrow q_{i}({\mathbf{A}},{\mathbf{x}})$, but this would imply that the zero sequence converges to a positive number, which is false. This contradiction shows that $T_{0}\subset T_{1}$.

Since $T_{0}\subset T_{1}$ and $T_{0}\neq T_{1}$, there must be an index $i$ in $T_{1}$ which is not in $T_{0}$. Let $T_{1}(i)$ be the index of $i$ in $T_{1}$, i.e., so that ${\mathbf{F}}_{:,i}$ and $\left({\mathbf{F}}_{:,T_{1}}\right)_{:,T_{1}(i)}$ are the same vector. Then we have that, by Lemma A.2,

[TABLE]

whereas $q_{i}({\mathbf{A}}^{k},{\mathbf{x}}^{k})=0$ since $i$ is not in $T_{0}$. Again, since $q$ is continuous, we must have that $q_{i}({\mathbf{A}}^{k},{\mathbf{x}}^{k})\rightarrow q_{i}({\mathbf{A}},{\mathbf{x}})$, so

[TABLE]

as $k\rightarrow\infty$. Since $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})\mapsto\left({{\mathbf{A}}^{\prime}_{:,T_{1}}}^{\dagger}x^{\prime}\right)_{T_{1}(i)}$ is a continuous map on $V_{r,s}\times\mathbb{R}_{+}^{r}$ (as the pseudoinversion, projection, and matrix multiplication operations are all continuous on this set), and $({\mathbf{A}}^{k},{\mathbf{x}}^{k})\rightarrow({\mathbf{A}},{\mathbf{x}})$, we must have

[TABLE]

Thus,

[TABLE]

Since $({\mathbf{A}},{\mathbf{x}})$ was chosen arbitrarily in $W$, we have shown that for every $({\mathbf{A}},{\mathbf{x}})$ in $W$, there is a set of indices $T$ and an index $i$ in $T$ such that

[TABLE]

Let $f_{T}:V_{r,s}\times\mathbb{R}_{+}^{r}\rightarrow\mathbb{R}^{|T|}$ be the map given by

[TABLE]

Let $Z_{d}$ be the set of all vectors in $\mathbb{R}^{d}$ which have at least one zero entry. Then we have shown that

[TABLE]

Since the union above is finite, it suffices to show that $f_{T}^{-1}(Z_{|T|})$ has zero measure for each $T\subset[s]$. So we are done if we can show that the preimage of measure zero sets under $f_{T}$ have measure zero, as $Z_{d}$ is a measure zero subset of $\mathbb{R}^{d}$ (it is just the union of the $k$ codimension-1 coordinate planes, each of which has measure zero).

By [36, Theorem 2], since $f_{T}$ is a map from a higher dimensional space to a lower dimensional space, this will be true if we can show that $f_{T}$ is continuous everywhere, differentiable almost everywhere, and has a full (row) rank derivative almost everywhere. Since pseudoinversion is continuous and differentiable on the full rank matrices, and column projection and matrix multiplication are differentiable everywhere, we immediately have that $f_{T}$ is continuous and differentiable everywhere on $V_{r,s}\times\mathbb{R}_{+}^{r}$. Thus, we are done if we can show that $f_{T}$ has a full rank derivative almost everywhere. We write the derivative of $f_{T}$ in block form,

[TABLE]

where we have linearized the indices in ${\mathbf{A}}$ (i.e., $\frac{\partial f_{T}}{\partial{\mathbf{A}}}$ is a $|T|\times rs$ dimensional matrix). This will have full (row) rank if either $\frac{\partial f_{T}}{\partial{\mathbf{A}}}$ or $\frac{\partial f_{T}}{\partial{\mathbf{x}}}$ has full (row) rank. Looking at the formula (A.6), we can immediately see that

[TABLE]

Since every ${\mathbf{A}}$ in $V_{r,s}$ has full (column) rank, ${\mathbf{A}}_{:,T}$ will also always have full (column) rank, so ${\mathbf{A}}_{:,T}^{\dagger}$ always has full (row) rank. Thus, $\frac{\partial f_{T}}{\partial{\mathbf{x}}}$ has full row-rank at every point in $V_{r,s}\times\mathbb{R}_{+}^{r}$, meaning $df_{T}$ does as well, so we are done.

∎

We may now prove Theorem 2.3.

Proof of Theorem 2.3.

Since $({\mathbf{A}},{\mathbf{x}})$ is in $U_{r,s}$, there is a neighborhood $N$ of $({\mathbf{A}},{\mathbf{x}})$ such that for every $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N$, $q({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ also has support $T$. Then for any $({\mathbf{A}}^{\prime},\bm{x}^{\prime})$ in $N$, by Lemma A.2, we have

[TABLE]

From this, we can immediately see that differentiating $q$ with respect to ${\mathbf{x}}^{\prime}$ and plugging in $({\mathbf{A}},{\mathbf{x}})$ gives equation (2.3).

On the other hand, by the second equation in (A.7), if $\alpha\notin T$, then $q_{\alpha}({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})=0$ for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N$, so

[TABLE]

for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N$, and therefore for $({\mathbf{A}},{\mathbf{x}})$ in particular. Furthermore, if $\beta\notin T$, then neither ${{\mathbf{A}}^{\prime}_{:,T}}^{\dagger}{\mathbf{x}}^{\prime}$ nor [math] depends on ${\mathbf{A}}^{\prime}_{i,\beta}$, so

[TABLE]

for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})\in N$, and therefore for $({\mathbf{A}},{\mathbf{x}})$ in particular. Thus, if either $\alpha\notin T$ or $\beta\notin T$, then

[TABLE]

Equivalently, the above equation holds if $(\alpha,\beta)\in(T\times T)^{c}$. This gives the second equation in (2.4).

So all that is left is to show is that

[TABLE]

We have that

[TABLE]

for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N$ by the first equation (A.7). In order to compute the last derivative above, we appeal to the formula for the derivative of the pseudoinverse operation. Specifically, if $B$ is a matrix, then we have that

[TABLE]

where ${\mathbf{E}}^{i,\alpha}$ is a matrix of the same size as ${\mathbf{F}}$, all of whose entries are [math], except for the $(i,\alpha)^{th}$ entry, which is $1$. See Theorem 4.3 in [19] for a proof of this formula. Eventually, we will plug in ${\mathbf{A}}^{\prime}_{:,T}$ in for ${\mathbf{F}}$ in the formula above, and so since ${\mathbf{A}}^{\prime}$ has full column rank (as $({\mathbf{A}},{\mathbf{x}})\in U_{r,s}$), so does ${\mathbf{A}}^{\prime}_{:,T}$, and so we may assume that ${\mathbf{F}}$ has full column rank. In particular, this means that ${\mathbf{I}}-{\mathbf{F}}^{\dagger}{\mathbf{F}}=0$, so the third term of the equation above drops out, giving

[TABLE]

Note that, by our definition of ${\mathbf{E}}^{i,\alpha}$, we have that for any matrices ${\mathbf{F}}$ and ${\mathbf{G}}$, we have

[TABLE]

and, by taking transposes,

[TABLE]

Thus,

[TABLE]

Then we get

[TABLE]

Now, the matrix $({\mathbf{F}}^{\dagger})_{:,i}({\mathbf{F}}^{\dagger}{\mathbf{x}}^{\prime})^{\top}$ has $(\beta,\alpha)^{th}$ entry equal to $({\mathbf{F}}^{\dagger})_{\beta,i}({\mathbf{F}}^{\dagger}{\mathbf{x}}^{\prime})_{\alpha}$, and so we get that

[TABLE]

Setting ${\mathbf{F}}={\mathbf{A}}^{\prime}_{:,T}$, this gives

[TABLE]

which holds for all $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ in $N$. Plugging in $({\mathbf{A}},{\mathbf{x}})$ for $({\mathbf{A}}^{\prime},{\mathbf{x}}^{\prime})$ gives the desired result.

∎

We finally can prove Theorem 2.6, which gives the required derivatives to apply backpropagation.

Proof of Theorem 2.6.

The chain rule tells us that

[TABLE]

since ${\mathbf{S}}^{(\ell_{2})}_{\alpha,m}$ for $\ell_{1}\leq\ell_{2}\leq{\mathcal{L}}$ are the only variables which $C$ depends on, which themselves depend on ${\mathbf{A}}^{(\ell_{1})}$. Therefore we are done if we can show that

[TABLE]

Now, we apply Lemma A.1 to the second equation in (2.5) to get

[TABLE]

By repeated application of the chain rule to this equation, we get that

[TABLE]

Applying Theorem 2.3 to equation A.8 gives

[TABLE]

Now we use the following fact: if ${\mathbf{F}}$ and ${\mathbf{G}}$ are compatible matrices, $R$ is a subset of the row indices of ${\mathbf{G}}$, and ${\mathbf{G}}_{R^{c},:}=0$, then

[TABLE]

Combining equation (A.9) and the second equation in (A.10) with this fact gives

[TABLE]

And now we can apply the first equation in (A.10) to get

[TABLE]

Now, we apply the chain rule again to A.8, giving

[TABLE]

On the other hand,

[TABLE]

where $\left(\frac{\partial L}{\partial{\mathbf{S}}^{(\ell_{2})}_{:,m}}\right)$ is taken to be a column vector. Thus,

[TABLE]

We can apply equation (A.11) again, using both parts of A.12, giving

[TABLE]

By Theorem 2.3, the rows of $\left(\frac{\partial{\mathbf{S}}^{(\ell_{1})}_{:,m}}{\partial{\mathbf{A}}^{(\ell_{1})}_{i,:}}\right)$ with indices outside $T_{m}^{(\ell_{1})}$, are zero, so we can apply our fact again to get

[TABLE]

Now, Theorem 2.3 again tells us that the columns of $\left(\frac{\partial{\mathbf{S}}^{(\ell_{1})}_{:,m}}{\partial{\mathbf{A}}^{(\ell_{1})}_{r,:}}\right)$ outside $T_{m}^{(\ell_{1})}$ are zero, so we have

[TABLE]

which is what we wanted, since ${\mathbf{U}}^{(\ell_{1},\ell_{2}),m}_{:,{T_{m}^{(\ell_{1})}}^{c}}=0$. Thus, all that is left to show is that

[TABLE]

i.e.,

[TABLE]

Applying Theorem 2.3 once more gives that

[TABLE]

Now, using equation (A.8) and Lemma A.2, we get that

[TABLE]

Thus, we have,

[TABLE]

From the definition of ${\mathbf{\Phi}}$, we see that

[TABLE]

so

[TABLE]

Since ${\mathbf{d}}^{(\ell_{1},\ell_{2}),m}$ is a column vector, we get

[TABLE]

Plugging this into A.17, we get

[TABLE]

Thus, from equations (A.15) and (A.13), it suffices to show that

[TABLE]

First, we note that, by equation (A.16), we have

[TABLE]

Since ${\mathbf{S}}^{(\ell_{1})}_{{T_{m}^{(\ell_{1})}}^{c},m}=0$, we apply equation (A.11) once more to get that

[TABLE]

Thus,

[TABLE]

This gives

[TABLE]

Using equation (A.18) and the definition of ${\mathbf{d}}$, we get that

[TABLE]

Thus,

[TABLE]

which is what we wanted.

∎