
\documentclass[a4paper, 12pt]{article}
% \usepackage{bbold}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{mathrsfs}
\usepackage{listings}
\usepackage{graphicx}
% \usepackage{natbib}

% \usepackage[margin=0.5in]{geometry}

\DeclareMathOperator{\Ker}{Ker}
\DeclareMathOperator{\Img}{Im}
\DeclareMathOperator{\Dom}{Dom}
\DeclareMathOperator{\supp}{Supp}
\DeclareMathOperator{\Id}{Id}
\DeclareMathOperator{\prox}{prox}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator{\sign}{sign}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{claim}{Claim}
\newtheorem*{case}{Case}
\newtheorem*{notation}{Notation}
\newtheorem*{prop}{Proposition}
\newtheorem*{problem}{Problem}
\newtheorem*{solution}{Solution}
\newtheorem*{definition}{Definition}


\begin{document}
\title{Matrix Organization Problems via Bi-Folderings}
\date{December 18, 2017}
\author{Julien Clancy}

\maketitle

\section{Introduction} % (fold)
\label{sec:introduction}

In a series of papers \cite{GAVISH2012354},\cite{Coifman2011},\cite{coifman2017} Coifman et al recently formulated a generalization of multiscale analysis to spaces of points equipped with a set of functions, {\it assuming} smoothness of the functions to build a metric on the points. Such a set of functions can be conveniently organized into a matrix, if the functions are taken to be the columns and the points are the rows. This induced metric on the rows can then be used to build a metric on the columns, this time considering the {\it rows} as functions on the point set of the {\it columns}. Alternating in this fashion, we obtain a pair of coupled geometries on the rows and columns, which hopefully reveal some structure of the matrix. The inventors' work focuses on ``real data'', e.g. \cite{coifman:7548324}, starting with a kind of earth-mover's distance between the rows or columns. We hope to apply these methods to solve problems relating to recovering low-rank blocks in matrices. The aim of this project is to use these methods on ``planted low-rank block'' problems to gauge their performance, and combine them with techniques more traditionally used for low-rank matrices. We obtain empirical solutions to open problems in the area, and cheap and simple solutions to existing problems.

% section introduction (end)


\section{The Problems} % (fold)
\label{sec:the_problems}

\subsection{Subspace Clustering} % (fold)
\label{sec:subspace_clustering}

The most elementary version of the problem at hand is subspace clustering (see \cite{candes2012} and references therein). We are given a $n \times m$ matrix $X$ with the property that there are subspaces $V_i$ such that all $X_j$ (columns of $X$) are in some $V_i$, and $\sum \dim V_i \ll m$; the goal is to recover $V_i$ and the assignment of each $X_j$ to some $V_i$ (or, in the case where there are also outlying points that belong to no subspace, just the subspaces $V_i$, from which the non-outlier membership can be easily recovered). Equivalently, we are guaranteed that after a permutation the columns of $X$ are in contiguous low-rank blocks, and we must recover the permutation. A matrix (with the correct column ordering) generated according to this problem is shown in figure 1, with outlines for relief.

\begin{figure}
	\centering
	\includegraphics[width=32em]{subspace_clustering_3.png}
	\caption{\small An instance of the subspace clustering problem on three subspaces, with the correct permutation. Subspaces and points in them were chosen randomly. Notice the repetitive rectangular patterns in the blocks, which are consistent within but inconsistent across blocks.}
\end{figure}

What \cite{candes2012} analyzed was the ``sparse subspace clustering'' (SSC) algorithm, where we take each point $X_i$ and try to write it as a linear combination of other points $X_j$ in the sparsest way possible. Specifically, they solve
\begin{equation*}
	\min \, \lVert C \rVert_1 \text{ s.t. } XC = X, C_{i,i} = 0
\end{equation*}
and use the sparsity pattern of $C$ as a similarity matrix for spectral clustering. They prove that with overwhelming probability over a natural generative model (random subspaces with random vectors, $\rho$ vectors per subspace) $C$ only contains nonzeros where two vectors truly are in the same subspace. The generalization to the noisy case, given in \cite{Candes2013}, is the Lasso-SSC program
\begin{equation*}
	\min \, \frac{1}{2} \lVert X - XC \rVert_F^2 + \lambda \lVert C \rVert_1 \text{ s.t. } C_{i,i} = 0
\end{equation*}
One of the goals of this project (and our ongoing research) is to integrate this penalty into tree-based clustering to improve performance.

% subsection subspace_clustering (end)

\subsection{Subspace Clustering with Irrelevant Features} % (fold)
\label{sub:subspace_clustering_with_irrelevant_features}

In \cite{NIPS2015_5707} the authors introduced and solved a more complex version of subspace clustering, where some one the features (rows) are meaningless, or in the model, replaced with noise. Otherwise the task is the same as above. Looking at each point in the Lasso-SSC program individually, expanding the Frobenius norm, and assuming that $\lVert x_i \rVert = 1$, we see that it is equivalent to
\begin{equation*}
	\min \,  \frac{1}{2} \lVert X_{\sim i} c_i \rVert_2^2 - \langle x_i, X_{\sim i} c_i \rangle + \lambda \lVert c_i \rVert_1
\end{equation*}
Their trick is to replace inner products $\langle a, b \rangle = \sum_{i=1}^n a_i b_i$ by a robust counterpart, where instead of summing all $a_i b_i$ we take only the $k$ smallest entries, for some $k < n$. This allows us to redefine the norm and inner product above and run the minimization. What the authors show is that indeed this works.

\begin{figure}
	\centering
	\includegraphics[width=32em]{ssc_irrelvant_features.png}
	\caption{\small Subspace clustering with irrelevant features (shown in ``planted'' form, without permutations on the rows and columns).}
\end{figure}

% subsection subspace_clustering_with_irrelevant_features (end)

\subsection{Subspace-Dependent Irrelevant Features} % (fold)
\label{sub:subspace_dependent_irrelevant_features}

In \cite{NIPS2015_5707} the authors also bring up a related problem, where which features are relevant depends on the subspace. This setup is shown in figure 3. The major goal of this project was to solve this problem, hopefully integrating the methods of \cite{candes2012} with those of \cite{Coifman2011}.

\begin{figure}
	\centering
	\includegraphics[width=32em]{ssc_irrelevant_dependent.png}
	\caption{\small Recovering planted low-rank blocks, where the blocks do not intersect across columns. As before this matrix is shown without permutations for clarity.}
\end{figure}

% subsection subspace_dependent_irrelevant_features (end)

% section the_problems (end)

\section{Bi-Foldering} % (fold)
\label{sec:bi_foldering}

Here we describe the simplified version of the algorithm in \cite{Coifman2011} for matrix organization. We take a matrix to be a function $M \colon X \times Y \to \mathbb{R}$, where $X$ is the set of rows and $Y$ is the set of columns. A tree $\mathcal{T}$ is a hierarchy of partitions on a point space which has the property that each ``folder'' in a partition intersects only one folder in the partition immediately above it --- think of a binary tree where the first level is just the whole space $X$, the second is a bipartition of $X$, the third is a division of $X$ into four disjoint subsets, and so on, with all partitions nested. Given a tree $\mathcal{T}$ we say its $\ell$th level is $\mathcal{T}^\ell$ and the $i$th folder on that level is $\mathcal{T}^\ell_i$. Notice that a tree naturally induces a geometry on the space that its folders partition, via the tree metric.\footnote{This is the distance of two points to their least common ancestor in the tree. One might think that if two people live in the same house their tree distance is, say, $2$, while if they live only in the same city it is $3$, and if they live only in the same state it is $4$.} Denoting by $d_\mathcal{T}$ this metric, let $W_\mathcal{T}(x_1,x_2) = d_\mathcal{T}(x_1,x_2)$ be its affinity matrix. It is well-known that the top eigenvectors of affinity matrices embed their metrics in low-dimensional space, so if we do this, we could hopefully recover the tree via recursive $2$-means clustering. It is not particularly useful to recover a tree from its embedding, but it is useful to point out that trees and affinity matrices are in some sense one and the same.

Given a partition tree $\mathcal{T}$ on a set $X$ and a collection of functions $f_j \colon X \to \mathbb{R}$, we can define new affinity matrix on $Y = \{f_j\}$ by
\begin{equation*}
	W_Y(f_j,f_k) = \sum_{\ell} \sum_{i} \frac{1}{|\mathcal{T}^\ell_i|^\alpha} \left|\langle f_j,f_k \rangle_{\mathcal{T}^\ell_i}\right|
\end{equation*}
where $\langle f, g \rangle_S = \sum_{i \in S} f_i g_i$ and $\alpha$ is some parameter that represents scale-sensitivity of the metric --- we might wish to undercount small folders because their contents are noisy or simply because there are so many of them, corresponding to $\alpha \leq 1$, or we might want to treat every folder independently of its size by choosing $\alpha = 0$, or with the central limit theorem in mind we might take $\alpha = 1/2$. Throughout all of the experiments recorded here we used $\alpha = 1$, but changing to $\alpha = 1/2$ did not substantively change the results. Using the procedure outlined above we embed the set $\{f_j\}$ in in Euclidean space using the eigenvectors of $W_Y$, then construct a binary tree $\mathcal{S}$ using recursive $2$-means. Since the functions $f_j$ are actually the columns of $M$, we then use $\mathcal{S}$ to organize the functions given by the {\it rows} of $M$, and alternate in this fashion until convergence (which usually means $5$ to $10$ iterations).

% section bi_foldering (end)

\section{Experiments} % (fold)
\label{sec:experiments}

As a sanity check we ran this algorithm on the first example problem above, subspace clustering. Its performance was competitive with the state of the art. We also ran experiments on the problem of subspace clustering with irrelevant features, and this time achieved significantly better performance than the state of the art --- while the method in \cite{NIPS2015_5707} can reportedly identify the subspaces with up to half of the features corrupted, our method works at least down to 90\% corruptions, and possily more. On the problem of finding low-rank blocks our method performed extremely well, outperforming reasonable expectations. For example, in our experiments we started with a $256 \times 1024$ matrix of noise, and planted eight $64 \times 128$ blocks containing random samples from $8$-dimensional subspaces. We ensured that the statistics of noise and non-noise patches matched by scaling the noise levels, and ran the simple algorithm above. Our metric for recovery of a block was if one level in the column partition tree very closely matched with which columns were in some block, since given this information we can easily recover the block itself. Indeed we found that every subspace had a level in the final partition tree that contained its and only its columns, modulo an error of no more than 5\%.

We also implemented the Lasso-SSC algorithm using FISTA \cite{Beck:2009:FIS:1658360.1658364} to construct affinity matrices at each level. While more experiments are necessary to calibrate scaling constants relating to different levels in the tree and the sparsity-fit tradeoff in the coefficient matrix $C$, we found that it performed slightly worse than the tree-based method and required much more computation time (since the affinity solution step is itself a convex solve, not a matrix multiply). In the future we plan to implement d'Aspremont et al's acceleration methods in \cite{NIPS2016_6267} or the restart schemes described in \cite{ODonoghue2015} to further speed up this solve. We believe that a SSC-type method has the potential to be very successful in this area, but further work is needed.

If we are dealing with very small numbers of subspaces then looking at the permuted matrix is enough to recognize whether the algorithm worked. However, especially in the case of irrelevant features that depend on the subspace, it is difficult to identify optically whether the algorithm actually succeeded. A better way is to allow for the user to explore the tree organization of the columns of the matrix. While Julia's plotting is quite good, we found that the most intuitive solution involves user interaction. This is much simpler to do in JavaScript using the D3.js library, so we installed the HTTPServer.jl package to interface the algorithms we implemented with the D3 visualization. The tree visalization in the lower-left is largely based off of the examples of Mike Bostock (the primary D3 developer). The code operates on a randomly permuted version of the example matrix (generated from one of the models above), but the interface shows the matrix in the correct order, so that contiguous blocks of columns being in the same folder corresponds to algorithmic correctness. The three-dimensional scatterplot shows the columns plotted according to the top three\footnote{Actually, the top four minus the first, which is almost constant beacause affinity matrices are typically close to being Markov.} eigenvectors of the affinity matrix generated by the learned tree. To navigate the data, single-click or double-click nodes in the tree to highlight the corresponding columns or expand a node's children.

% section experiments (end)

\section{Code Structure} % (fold)
\label{sec:code_structure}

The server is currently configured to generate and analyze an instance of the third, hardest model. The visualization is shown at 127.0.0.1:8000.

The algorithm itself operates by traversing partition trees, so the workhorses of our code are types for Folders, Partitions, and PartitionTrees. The functions for generating the tree-based similarity matrix and deriving a dual geometry from it proceed as explained above; see \cite{Coifman2011} for more details. The sparse subspace clustering code uses an $\ell^1$ minimization solver based on the accelerated gradient descent algorithm derived in\cite{Beck:2009:FIS:1658360.1658364}. We also implemented a fast top eigenvector solver desribed in \cite{NIPS2015_5735}, but found it to be slightly inferior in runtime compared to Julia's built-in svds function at a given accuracy level (benchmarked against an unlimited runtime full svd). We believe memory is the key culprit here, which is an unavoidable feature of the algorithm --- theoretically it is only a log factor worse than Lanczos iteration, but the constants are probably unfavorable. We also implemented unit tests for the portions of our code amenable to it.\\

Our implementation of the main algorithm is based off of one in MATLAB due to Ronald Coifman and his collaborators, which we found difficult to modify without reimplementation.

% section code_structure (end)

\section{Future Work} % (fold)
\label{sec:future_work}

The feeling of the author is that subspace clustering should be present somehow in any solution of the low-rank blocks problem --- this is the most prominent avenue for future research. We found a naive integration with tree geometries to be unsuccessful, but it could be that it would work better without a multiscale framework, say with selective masking of entries.

The models considered here are in some sense toys. There are many real-world examples where low-rank blocks appear continuously rather than with well-defined boundaries, for instance in hierarchical matrices arirising from potential operators \cite{Hackbusch:2015:HMA:2988144}. Matrices with arbitrary low-rank block decompositions do not necessarily fall on a perfect tensor grid, and so the bi-tree setup is probably unsuitabe to analyze them. This is the other direction in which we intend to continue research.

% section future_work (end)

\small
\bibliography{bib} 
\bibliographystyle{ieeetr}

\end{document}


