Laplacian Eigenmaps for Neural Circuit Analysis

Algorithm

In Belkin & Niyogi, 2002, the authors provide an algorithm for computing an embedding for a dataset. The general algorithm proceeds in 4 main steps: (1) computing a nearest neighbors’ graph (2) computing the graph Laplacian, (3) computing all or the top $r$ eigenvalues and eigenvectors (eigenpairs) of the returned Laplacian graph (4) using the smallest $r$ eigenpairs to embed the graph into a lower dimensional space. Where the last 3 steps comprise the Laplacian eigenmap.

Algorithm: Laplacian Eigenmap

1. Construct the Laplacian matrix $L = D - A$, where $D$ is the diagonal matrix of degree sums.
2. Solve the generalized eigenvector problem: $L\mathbf{y} = \lambda D\mathbf{y}$.
3. Use the eigenvectors corresponding to the smallest non-zero eigenvalues to define the embedding.

Nearest Neighbors Graph

The nearest neighbors graph is constructed based on a similarity function. For $k$ points $x_1,\dots, x_k$ in a dataset $D$, the edges are drawn based on a similarity function. The proposed method in Belkin & Niyogi, 2002 is an $\varepsilon$-neighborhood where points are connected if they are less than $\varepsilon$ away from each other under the Euclidean metric. Selection of the $\varepsilon$ parameter is crucial but will not be discussed here.

After the edges of the graph have been realized, weights can be assigned based on domain knowledge. For example, the number of synapses between two neurons can serve as the weights of an undirected graph of neurons. Alternatively, a heuristic can be used like the proposed heat kernel shown below, which penalizes points being far away from each other¹.

In both cases, the goal is to construct a graph whose Laplacian is a positive semidefinite matrix.

Laplacian Eigenmap

Once the graph has been constructed, the graph Laplacian $L$ can be computed as the difference between the degree matrix $D$ and the weight matrix $W$. The computational complexity of this step is $O(n^2)$ for subtraction in each entry.

Computing Eigenvalues and Eigenvectors

The method proposed by Belkin and Niyogi reduces the problem to a generalized eigenvalue problem. Different solvers have different implementations to compute the eigenvalues and eigenvectors of the nearest neighbor’s graph. For example, the SciPy implementation solves for these using QR-iteration Mohanty & Gopalan, 2012 whereas MATLAB’s offers the use of the QZ-algorithm Moler & Stewart, 1973. Both algorithms first convert the matrices into upper Hessenberg form which has computational complexity $O(n^2)$ Mohanty & Gopalan, 2012.

Overall, both algorithms have computational complexity of $O(n^3)$ although the QZ-algorithm has better numeric stability Moler & Stewart, 1973. The sparse algorithms in SciPy, which is wrapper around ARPACK Lehoucq et al., 1997 allows for the computational complexity of the eigenpair computation to be reduced to $O(n^2)$ depending on the sparsity of the graph. For practical implementations, for $A\in\mathbb{R}^{5\times 5}$, the feasibility of computing the eigenpairs becomes less reasonable for non-iterative methods.

Embedding

Framed as a minimization problem, given the graph Laplacian $L$, the goal is to find a set of vectors $\mathbf{y}_1, \mathbf{y}_2, \ldots, \mathbf{y}_k$ that minimize the objective function:

subject to the constraint:

Since the embedding step is the solution to the minimization problem and this problem is solved by computing the eigenvectors, there is no additional computational complexity from this embedding.

Testing

We implement the algorithm exactly as we describe above with a subset of the fly_dataset representing a core component of the T4a neural circuit. In contrast to the above, we have a graph that is symmetric but also sparse. In practice, there are algorithms that can leverage the sparsity of the graph structure to reduce the computational expense.

We will examine the speed of the algorithm during the Laplacian eigenmap step as this is both the step that produces the embedding of the graph and also the most computationally expensive step. All tests will be run on a DELL XPS 15 9570 with 16 gigabytes of RAM, an Intel Core™ i7-8750H CPU @ 2.20GHz×12 running Ubuntu 22.04 LTS with a 64-bit architecture. All code will be written in Python 3.10.12 with the duration of 200 reruns of the algorithm averaged and the final time reported.

The adjacency graph is an $N=2391$ square 32-bit sparse float array. Computing the Laplacian of given the adjacency graph using the Laplacian function from SciPy’s sparse graph library takes $\sim 5.00 \times 10^{-4}$ while a similar computation using the NumPy library takes $\sim 9.93 \times 10^{-3}$ seconds with 32-bit floats.

Next, the eigenvalues and eigenvectors of the Laplacian are computed. Similarly, we use the SciPy sparse implementation and the NumPy dense implementation.

However, the sparse implementation for sparse matrices requires that only the first $r$ eigenvector-eigenvalue pairs are returned. This is because the SciPy implementation uses Arnoldi iteration Arnoldi, 1951. Since we would like to embed the graph into a 2D or 3D graph, we set $r=2$. The sparse implementation takes $\sim 1.52 \times 10^{-1}$ seconds while the NumPy implementation takes over 468 minutes (about 8 hours) and 39 seconds.

We can further explore the difference between the sparse and non-sparse implementations on a synthetic dataset. We can create sparse and dense graphs by assigning a probability $p$ of nodes being connected. This is easy to do with the networkx library in Python by its erdos_renyi_graph function.

For our dense graphs, we set $p=0.9$ and for the sparse graphs we set $p=0.1$. On the same device and with 10 loops per node count ranging from 10 to 3000 nodes. The graphs for the dense and sparse implementations are as follows:

The two figures show the comparison of average computation times for obtaining the graph Laplacian using SciPy’s sparse and Numpy dense matrix implementations over varying node sizes. The dense matrix implementation is surprisingly faster than the sparse implementation when compared across varying node sizes.

Next, we compare SciPy’s handling of sparse and dense graphs on another synthetic dataset where the number of nodes range from 100 to 2000, increasing in increments of 100. As opposed to the previous tests, all tests here are ran in SciPy because its sparse linalg generalized eigenvalue solver handles both dense and sparse matrices.

In contrast to the previous experiment, the dense implementation takes longer than the sparse one despite both using the same generalized eigenvalue eigenvector (GEEV) iterative solver Lehoucq et al., 1997. Indeed, both plots show the average computation time for obtaining the smallest two eigenvalues using SciPy’s GEEV. They also show that the computational complexity is related to the number of nodes being plotted. However, when the sparse graphs are compared to the dense graphs, we can see that the sparse graphs show lower complexity than the dense graphs.

References

Arnoldi, W. E. (1951). The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quarterly of Applied Mathematics, 9(1), 17-29.

Belkin, M., & Niyogi, P. (2002). Laplacian Eigenmaps and Spectral Techniques for Embedding and Clustering. In T. G. Dietterich, S. Becker, & Z. Ghahramani (Eds.), Advances in Neural Information Processing Systems 14 (pp. 585-592). The MIT Press. https://doi.org/10.7551/mitpress/1120.003.0080

Borst, A., Drews, M., & Meier, M. (2020). The Neural Network behind the Eyes of a Fly. Current Opinion in Physiology, 16, 33-42. https://doi.org/10.1016/j.cophys.2020.05.004

Chung, F. R. K. (1995). Eigenvalues of graphs. In Proceedings of the International Congress of Mathematicians: August 3-11, 1994 Zürich, Switzerland (pp. 1333-1342). Springer.

Ding, C. (2004). A tutorial on spectral clustering. Talk presented at ICML.

Dorkenwald, S., McKellar, C. E., Macrina, T., et al. (2022). FlyWire: online community for whole-brain connectomics. Nature Methods, 19(1), 119-128. https://doi.org/10.1038/s41592-021-01330-0

FlyWire (2024). Visual Columns Challenge. Retrieved from https://codex.flywire.ai/app/visual_columns_challenge

Francis, J. G. F. (1961). The QR Transformation A Unitary Analogue to the LR Transformation—Part 1. The Computer Journal, 4(3), 265-271. https://doi.org/10.1093/comjnl/4.3.265

Shinomiya, K., Huang, G., Lu, Z., et al. (2019). Comparisons between the ON- and OFF-edge motion pathways in the Drosophila brain. eLife, 8, e40025. https://doi.org/10.7554/eLife.40025

Lehoucq, R. B., Sorensen, D. C., & Yang, C. (1997). ARPACK: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods. Available from netlib@ornl.gov.

Mohanty, S. K., & Gopalan, S. (2012). I/O efficient QR and QZ algorithms. In 2012 19th International Conference on High Performance Computing (pp. 1-9). https://doi.org/10.1109/HiPC.2012.6507492

Mohar, B. (1991). Eigenvalues, diameter, and mean distance in graphs. Graphs and Combinatorics, 7(1), 53-64.

Moler, C. B., & Stewart, G. W. (1973). An Algorithm for Generalized Matrix Eigenvalue Problems. SIAM Journal on Numerical Analysis, 10(2), 241-256. https://doi.org/10.1137/0710024

Nériec, N., & Desplan, C. (2016). From The Eye To The Brain: Development Of The Drosophila Visual System. Current Topics in Developmental Biology, 116, 247-271. https://doi.org/10.1016/bs.ctdb.2015.11.032

Scheffer, L. K., Xu, C. S., Januszewski, M., et al. (2020). A connectome and analysis of the adult Drosophila central brain. eLife, 9, e57443. https://doi.org/10.7554/eLife.57443

Von Luxburg, U. (2007). A tutorial on spectral clustering. Statistics and Computing, 17, 395-416.

Wikipedia (2023). Rand index. Retrieved from https://en.wikipedia.org/wiki/Rand_index