Faster Inference of Cell Complexes from Flows via Matrix Factorization

TL;DR

This paper introduces a fast, matrix-factorization heuristic for inferring cell complexes from flow data on graphs, improving efficiency and performance in noisy scenarios compared to previous methods.

Contribution

It proposes a novel matrix-factorization-based heuristic to efficiently infer cell complexes from flow data, addressing NP-hardness and outperforming prior heuristics in speed and accuracy.

Findings

Significantly reduces computational time compared to prior heuristics.

Achieves comparable solution quality in most settings.

Outperforms previous methods in noisy data scenarios.

Abstract

We consider the following inference problem: Given a set of edge-flow signals observed on a graph, lift the graph to a cell complex, such that the observed edge-flow signals can be represented as a sparse combination of gradient and curl flows on the cell complex. Specifically, we aim to augment the observed graph by a set of 2-cells (polygons encircled by closed, non-intersecting paths), such that the eigenvectors of the Hodge Laplacian of the associated cell complex provide a sparse, interpretable representation of the observed edge flows on the graph. As it has been shown that the general problem is NP-hard in prior work, we here develop a novel matrix-factorization-based heuristic to solve the problem. Using computational experiments, we demonstrate that our new approach is significantly less computationally expensive than prior heuristics, while achieving only marginally worse performance in most settings. In fact, we find that for specifically noisy settings, our new approach outperforms the previous state of the art in both solution quality and computational speed.

Equations

$$C_0$$

$$\mathcal{L}(\mathcal{C},\mathbf{F})=\Vert {\text{harm}}_{\mathcal{C}}(\mathbf{F}){\Vert }_F={\left(\sum _{i=1}^s{\left|\left|{\text{harm}}_{\mathcal{C}}({\mathbf{f}}_i)\right|\right|}_2^2\right)}^{1/2}$$

$$\underset{\mathcal{C}}{\min }\mathcal{L}(\mathcal{C},\mathbf{F})\text{s.t.}\mathcal{C}\text{ has G as 1-skeleton}\text{ and }\left|{\mathcal{C}}_2\right|\le k$$

$${\mathcal{C}}^{(i)}\leftarrow {\mathcal{C}}^{(i-1)}\cup \{{\theta }_j^{(i)}:j\in [1,\dots ,l^{\prime }]\}$$

$$\underset{{\mathbf{B}}_2,\mathbf{C}}{\min }\Vert \mathbf{F}-{\mathbf{B}}_2\mathbf{C}\Vert \text{s.t.}{\mathbf{B}}_2\in {\mathcal{B}}_2^{(k^{\prime })},\mathbf{C}\in {\mathbb{R}}^{k^{\prime }\times s}$$

$${\mathbf{H}}^{(i)}\leftarrow {\mathbf{H}}^{(i-1)}-\left[{\hat{\mathbf{b}}}_1^{(i)}\cdots {\hat{\mathbf{b}}}_{l^{\prime }}^{(i)}\right]\cdot {\left[{\hat{\mathbf{b}}}_1^{(i)}\cdots {\hat{\mathbf{b}}}_{l^{\prime }}^{(i)}\right]}^{†}\cdot {\mathbf{B}}^{(i)}{\mathbf{C}}^{(i)}$$

Full text

I Introduction

Graphs have become a prevalent abstraction in data science due to their ability to model many real-life systems [1]. Graph signal processing (GSP) [2, 3] enables signal processing for signals that are defined on such graphs, such as temperatures measured at different locations or neuron activity in different parts of the brain. However, in many applications, the recorded data represents flows, e.g. of people [4], traffic [5], or money [6]. Such flows are often more naturally represented with an (oriented) signal on the edges of a graph. To process such edge signals, recent extensions of GSP towards topological signal processing (TSP) [7, 8, 9, 10, 11, 12] utilize simplicial complexes or cell complexes, and their associated Hodge Laplacians.

These Hodge Laplacian operators, a generalization of the graph Laplacian to higher dimensions, are a central pillar of TSP, and can be used in lieu of the graph Laplacian as shift operators in signal processing tasks for signals defined on the edges or higher-dimensional cells of a complex. Importantly, the 1-Hodge Laplacian induces a decomposition of the space of edge flows into gradient, curl, and harmonic flows [13, 14]. Similarly to the graph Laplacian in GSP, the eigenvalues and corresponding eigenspaces of this Hodge Laplacian enable the construction of low- and high-pass filters, denoising, signal compression, and other classical signal processing tasks [9].

However, the higher-order cell structure of the complexes used to construct these Hodge-Laplacians is typically not readily available. Specifically, when considering edge-flow data, commonly only the underlying graph structure and the associated edge flows are observable. Therefore, analogously to inferring a graph structure from data on the nodes [3], different methods to infer simplicial complex from edge data on a graph have been proposed [8, 15]. The general problem of inferring a cell complex from such data, i.e., finding an optimal set of $2$-cells (polygons) such that the edge-flows can be represented by a sparse set of eigenvectors of the Hodge Laplacian, was introduced in [16]. It was shown that this problem is NP-hard, and thus [16] introduced a heuristic algorithm to solve it.

Contribution. In this paper, we introduce an alternative approach to solving the cell inference problem using matrix decompositions. As we show, our novel approach is significantly faster and incurs only a marginally higher approximation error than the originally proposed heuristic. In our experiments, this trade-off is particularly advantageous on larger networks, when inferring a large absolute number of 2-cells, and in the presence of high noise.

II Background

This introduction to $2$-dimensional cell complexes is inspired from [17].

We consider $2$-dimensional complexes constructed by augmenting simple undirected graphs as follows. Consider a graph ${\mathcal{G}=(\mathcal{V},\mathcal{E})}$ with a set of vertices $\mathcal{V}=\{v_{1},\ldots,v_{n}\}$ and a set of edges by $\mathcal{E}=\{e_{1},\ldots,e_{m}\}$. Each edge $e_{k}$ consists of a pair of vertices $e=(v_{i},v_{j})$, with arbitrary but fixed order. We call $v_{i}$ the source and $v_{j}$ the target of $e_{k}$. We encode the structure of the graph in an incidence matrix $\mathbf{B}_{1}\in\{0,\pm 1\}^{n\times m}$. For every edge $e_{k}$, oriented from $v_{i}$ to $v_{j}$, we have $(\mathbf{B}_{1})_{i,k}=1$ and $(\mathbf{B}_{1})_{j,k}=-1$, and $(\mathbf{B}_{1})_{-,k}=0$ otherwise.

Definition 1** (Simple cell complexes of dimension $2$).**

An abstract regular cell complex $\mathcal{C}$ of dimension $2$ consists of a graph $\mathcal{G}$, and a non-empty ordered set of $2$-cells (polygons) $\mathcal{C}_{2}$, encoded via a boundary matrix $\mathbf{B}_{2}\in\{0,\pm 1\}^{m\times|\mathcal{C}_{2}|}$. Specifically for every $2$-cell $\theta\in\mathcal{C}_{2}$, we have a set of edges $\{e_{1},\ldots,e_{m}\}$ corresponding to a simple cycle in $\mathcal{G}$ which forms the boundary of $\theta$. Each column of $\mathbf{B}_{2}$ corresponds to such a cycle and the entries of the associated edges are set to $\pm 1$ such that $\mathbf{B}_{1}\mathbf{B}_{2}=0$.

Remark 2* (Orientation).*

Note that the above construction amounts to equipping each edge and each polygon with a reference orientation. This built-in notion of orientation allows for a natural representation of physical data like flows via positive or negative values (with or opposite to the orientation).

Oriented signals on CCs can be represented as chains:

Definition 3** (Signal (or chain) space of a cell complex).**

Given a cell complex $\mathcal{C}$, we denote by $C_{k}=\mathbb{R}^{|\mathcal{C}_{k}|}$ the $k$-th signal space of $\mathcal{C}$ and obtain an associated sequence of signal spaces

[TABLE]

The boundary matrices $\mathbf{B}_{k}$ are linear maps between the signal spaces $C_{k}$. In this paper, we analyze flows from the signal space on edges $C_{1}$.

Hodge Laplacian and Hodge decomposition. The Hodge Laplacian $\mathbf{L}_{k}=\mathbf{B}_{k}^{\top}\mathbf{B}_{k}+\mathbf{B}_{k+1}\mathbf{B}_{k+1}^{\top}$ is a generalization of the Graph Laplacian [13, 14, 18]. The Hodge decomposition divides the space of edge flows $C_{1}$ into three eigenspaces associate to $\mathbf{L}_{1}$: The gradient space $\operatorname{Im}\mathbf{B}_{1}^{\top}$, the curl space $\operatorname{Im}\mathbf{B}_{2}$, and the harmonic space $\ker\mathbf{B}_{1}^{\top}\cap\ker\mathbf{B}_{2}$. Intuitively, the gradient space consists of flows based on the difference between potentials on the nodes and the curl space consists of flows around the boundaries of $2$-cells.

For a flow $\mathbf{f}\in C_{1}$ on a cell complex $\mathcal{C}$, we define the gradient flow $\operatorname{grad}_{\mathcal{C}}(\mathbf{f})=\mathbf{B}_{1}^{\top}(\mathbf{B}_{1}^{\top})^{\dagger}\mathbf{f}$, curl flow $\operatorname{curl}_{\mathcal{C}}(\mathbf{f})=\mathbf{B}_{2}(\mathbf{B}_{2})^{\dagger}\mathbf{f}$, and harmonic flow ${\operatorname{harm}_{\mathcal{C}}(\mathbf{f})=(\mathbf{I}-\mathbf{L}_{1}(\mathbf{L}_{1})^{\dagger})\mathbf{f}}$. $(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudoinverse.

III The Cell Inference Problem

Given a graph and a set of observed flows, the cell inference problem is to find a sparse set of $2$-cells that minimize the projection of the flows into the harmonic space when added. We denote the given graph by $\mathcal{G}=(\mathcal{V},\mathcal{E})$ and the $s\in\mathbb{N}$ sampled edge flows by $\mathbf{f}_{i}\in C_{1}$ for $i\in\{1,\ldots,s\}$. In the following, we use $\mathbf{F}$ to denote the matrix of all flows $\mathbf{F}=\left[\mathbf{f}_{1},\ldots,\mathbf{f}_{s}\right]$.

More formally, the problem can be modeled as an optimization problem: find a $2$-dimensional cell complex $\mathcal{C}$ which minimizes the harmonic projection of the flows $\mathbf{F}$ onto the harmonic space of $\mathcal{C}$.:

[TABLE]

As minimizing this loss without constraints on the number of cells would result in a $\mathcal{C}$ with an empty harmonic space, we formulate the minimization problem such that the number of 2-cells is constrained:

[TABLE]

Without loss of generality, we will in the following assume that the flows $\mathbf{F}$ are gradient free, as any gradient component will not alter the optimal $\mathcal{C}$ in our optimization problem.

IV Methods

As proposed in [16], we can proceed in an iterative fashion to solve the cell inference problem. We start with (i) a cell complex $\mathcal{C}^{(0)}$ equivalent to $\mathcal{G}$, (ii) the number of cells to add per iteration $l^{\prime}$, and (iii) the number $l$ of cell candidates to consider, with $l^{\prime}\leq l$. In each iteration $i$, we add $l^{\prime}$ new 2-cells $\theta^{(i)}_{1},...,\theta^{(i)}_{l^{\prime}}$ until the desired amount of 2-cells are inferred:

[TABLE]

To operationalize this algorithmic idea, we have to address two questions. First, how do we obtain possible cell candidates $\theta^{(i)}_{1},...,\theta^{(i)}_{l}$ in each iteration? Second, how do we evaluate the candidates in each iteration and select the best?

In the original work [16], a Spanning Tree Heuristic (SPH) was proposed that constructs candidate 2-cells from one or multiple constructed spanning trees. The candidate cells that (when added to the cell complex) maximally decrease the loss are then selected greedily. This requires computing several harmonic projections of an (augmented) cell complex per iteration, which is the computationally most expensive part of the algorithm (even though efficient algorithms like LSMR [19] can be used to compute the projections).

There is also a more subtle issue associated with the SPH approach, that leads to many iterations required to infer appropriate cells. When using spanning trees to generate candidate cells, each candidate cell contains exactly one edge outside the spanning tree, and at least two edges from the tree. Importantly, candidate cells that significantly lower the loss are likely to share edges in the spanning tree. Adding one of those cells to the complex will account for the flow on all the shared spanning tree edges. The other candidates are thus unlikely to decrease the loss significantly, as they share many edges with the already added cell. Hence, in [16], only one cell is added in each iteration.

Matrix-factorization-based approach

In the following we present a novel, modular framework called Matrix-Factorization-Cell-Inference (MFCI) (cf. Figure 1, Algorithm 1) to solve the cell inference problem. For this, we interpret Equation 5 as a problem to find a matrix $\mathbf{B}_{2}$ that minimizes:

[TABLE]

where $\mathcal{B}_{2}^{(k^{\prime})}$ is the set of valid edge-to-cell boundary matrices of cell complexes with $k^{\prime}\leq k$ 2-cells and a 1-skeleton given by $\mathcal{G}$; $\mathbf{C}$ is a matrix of real coefficients. Given a matrix $\mathbf{B}_{2}$, an optimal $\mathbf{C}$ can be obtained by solving a least squares problem.

Our key idea now is to first find a low-rank approximation of the flows $\mathbf{F}\approx\mathbf{B}\cdot\mathbf{C}$ using a matrix factorization, without considering the constraints on $\mathbf{B}$; and then discretize the columns of $\mathbf{B}$ to correspond to valid boundary vectors of 2-cells. This idea can be iteratively applied several times, adding multiple cells in each iteration. To create a practical algorithm from this idea we need to choose a suitable matrix factorization method, a heuristic to find cell candidates from this factorization, and a procedure to evaluate and select the best candidates. As a secondary consideration, we also consider an efficient approximation for calculating the harmonic flows.

Calculating (approximate) harmonic flows. In each step $i$, SPH calculates the harmonic flows ${\mathbf{H}^{(i)}:=\operatorname{harm}_{\mathcal{C}^{(i)}}(\mathbf{F})}$. To increase speed, MFCI can (optionally) use an approximate update of the harmonic flow instead. Consider that, in sparse configurations, there are few overlaps between $2$-cells from different iterations. Thus, the change in harmonic flow depends mostly on the boundaries $\hat{\mathbf{b}}^{(i)}_{1},\ldots,\hat{\mathbf{b}}^{(i)}_{l^{\prime}}$ of the chosen cells. With $\mathbf{H}^{(0)}=\mathbf{F}$, we can approximate the harmonic flow as:

[TABLE]

Matrix factorization. To find a suitable matrix factorization $\mathbf{F}\approx\mathbf{B}\cdot\mathbf{C}$ for our application, let us first collect some desirable properties for this factorization. First, the approximation error $\left\|\mathbf{F-B\cdot C}\right\|$ should be low. To make the discretization of the columns of $\mathbf{B}$ to 2-cells easy, the entries of $\mathbf{B}$ should also almost fulfill the properties of a boundary matrix already: First, the entries of $\mathbf{B}$ should be [math] or $\pm 1$. Second, the boundary of $\mathbf{B}$ should be zero: $\mathbf{B}_{1}\mathbf{B}=0$.

The first possible factorization we consider is a (truncated) Singular Value Decomposition (SVD) of $\mathbf{F}$, which is the optimal low-rank approximation. However, the entries are typically not close to [math] or $\pm 1$. Hence, we also consider Independent Component Analysis (ICA) [20], which separates a matrix into statistically independent components that are not necessarily orthogonal. In practice, this is often closer to the desired decomposition into a boundary matrix $\mathbf{B}_{2}$ and associated circulations around $2$-cells.

Obtaining 2-cell candidates. We apply the chosen matrix factorization strategy iteratively to find 2-cell candidates. In each iteration, we calculate a low-rank matrix approximation $\mathbf{B}^{(i)}\cdot\mathbf{C}^{(i)}\approx\mathbf{H}^{(i-1)}$ (e.g. using the SVD). In each iteration, we first calculate the approximation error of each column of the matrix $\mathbf{B}^{(i)}$ as $\|\mathbf{H}^{(i-1)}-\mathbf{B}^{(i)}_{-,j}\mathbf{C}^{(i)}_{j,-}\|_{1}$, where $\|A\|_{1}=\sum_{i,j}|A_{i,j}|$ denotes the $L_{1}$ norm. We keep the $l$ columns $\mathbf{b}^{(i)}_{1},\ldots,\mathbf{b}^{(i)}_{l}$ of $\mathbf{B}^{(i)}$ with the lowest approximation error and discretize them to valid 2-cell candidates using one of two heuristics:

• •

The deterministic heuristic creates a candidate by adding edges in decreasing order of their absolute values in $\mathbf{b}^{(i)}_{j}$ to an empty graph until a (unique) cycle is formed.

• •

The random-walk based heuristic simulates a random walker, using the corresponding values from chosen vectors of the matrix factorization $\mathbf{b}^{(i)}_{j}$ as weights for the probability distribution. Once the random walk forms a cycle, this cycle is returned as the candidate. Each cycle identifies a 2-cell $\theta^{\prime(i)}_{j}$.

We add the $l^{\prime}$ best candidates to the cell complex (see Equation 6); for $l=l^{\prime}$, we skip the evaluation.

Our approach infers a 2-cell candidate for each column of the approximation. The number of 2-cell additions in each iteration thus depends on the number of chosen columns, and the rank of the matrix approximation. Since the rank of $\mathbf{F}$ is bounded by $s$, multiple cells can only be added for $s>1$. As with a larger number of samples $s$, the low-rank approximation is more likely to eliminate noise in the observed flow data, we expect MFCI to be most useful when many flows are observed.

V Experiments

For the empirical evaluation, we compare both the approximation error and the compute time of our new approach to the state-of-the-art method SPH111Unless stated otherwise, we use the similarity SPH heuristic with 11 clusters and 11 cell candidates like in the original paper. In our experiments, a larger number did not significantly improve accuracy. [16] and the non-discretized SVD, which is the optimal continuous solution. We choose the rank of the matrix factorization to be equal to the number of candidates $l$.

We generate synthetic random CCs using the algorithm from [21]. We sample flow signals $\mathbf{f}_{i}=\mathbf{B}_{2}c_{i}+\mathbf{f}^{\prime}_{i}$, where $c_{i}\in C_{2}$ is a signal obtained from an i.i.d. gaussian distribution on all $2$-cells and $\mathbf{f}^{\prime}_{i}\in C_{1}$ is i.i.d. gaussian noise on the edges.

We also consider real-world data based on taxi trips in New York City [22, 23]. The transitions between neighborhoods in a driver’s trajectory are represented as a directed flow. These flows are then aggregated into $s$ flows via summation.

Synthetic Data. Figure 2 shows the approximation errors and compute times on synthetic data. As the ICA-based factorization approach outperforms the SVD-based approach in both approximation error and compute time, we only show the former. For the most significant speed-up, we skip evaluation of candidates and use the approximate update of the harmonic flow. With this configuration (denoted “8oo-1” in Figure 2), no calls to LSMR are required. Thus, the computation is fast while the approximation error is comparable to SPH. By evaluating candidates and only adding the best in each iteration (like SPH does), MFCI outperforms SPH in terms of approximation error, with a similar (if slightly larger) computational requirement. In fact, the approximation error is close to the ground truth of the sampled cells.

Real-World Data. In Figure 3, we compare MFCI to SPH on the taxi dataset. The maximum spanning tree heuristic is configured to only evaluate one candidate per iteration, which requires only one call to LSMR.

In contrast to the synthetic data, the SPH heuristics are more accurate than MFCI. Furthermore, in MFCI, using the SVD significantly outperforms the ICA, which is just slightly better than the random baseline. However, it is also notable that, after 40 added 2-cells, the difference in error between SPH and MFCI stays approximately constant. The computational times behave similarly to the synthetic data, with MFCI (without candidate evaluation) being significantly faster than SPH.

Harmonic Flow Approximation. We evaluate the usefulness of the proposed approximate calculation of the harmonic flow using the pseudoinverse. Figure 4 compares the configurations of MFCI from Figure 2 to the equivalent approximation variants. For MFCI with candidate evaluation, we see a small increase in approximation error with a negligible performance improvement when using the approximation. However, for the faster variant without candidate evaluation, the error does not increase, but the computation time is less than half of the original. Overall, this suggests that the approximation is a good trade-off for the faster variant of MFCI.

Robustness to Noise. Figure 5a shows the relative performance of MFCI compared to SPH for different noise levels. We observe that for low noise levels, SPH is more accurate than MFCI, but the opposite is true for higher noise levels. We hypothesize that the low-rank matrix factorization filters out noise in the data, leading to better candidate cells than in SPH.

VI Conclusion

We presented a new framework, MFCI, that approximately solves the general cell inference problem by using matrix factorization. Overall, our experiments show that compared to the previous state-of-the-art, MFCI achieves a significant speed-up at the expense of a small increase in the approximation error in certain settings. Furthermore, in noisy configurations, MFCI achieves a better approximation error.

Importantly, we saw qualitatively different results on synthetic and real-world data. The synthetic data was based on the same assumptions made in the development of MFCI, namely, that there are underlying $2$-cells that generate flow independently of each other. It is likely that real-world data has strong correlations, leading to a less useful matrix decomposition. This hypothesis is consistent with the observation that ICA, which employs similar inherent assumptions, performs poorly on real-world data. Nevertheless, MFCI with SVD is useful for applications where its significant speed-up (and thus increased scalability) outweighs its relatively small increase in approximation error. The approximation error is based on signal compression and likely translates well to signal processing tasks, but it is unclear how it affects the performance of other downstream tasks built on the inferred cells.

There is a large space of possible configurations for MFCI that exceeds the scope of this paper. First, any component of MFCI could be replaced, mainly the matrix factorization method and the heuristic to obtain candidates. Furthermore, the projection could use a hybrid approach where the approximate method is used, but the projection is done explicitly in some iterations. Second, due to the iterative nature of both MFCI and SPH, it is also possible to combine the two methods. As we have seen, SPH performs particularly well in the first iterations, where it is also comparatively fast. Therefore, it may be advantageous to design a hybrid approach, where the first few iterations are performed with SPH and the remaining iterations with MFCI. Ideally, this could combine the accuracy of SPH in the early iterations with the better computational performance of MFCI in the later, more computationally expensive, iterations.