During the last week of SGI, Tiago De Souza Fernandes, Miles Silberling-Cook, and I studied computing all pairs of geodesic distances on a manifold mentored by Nicholas Sharp, a former postdoc at the University of Toronto and a senior research scientist at NVIDIA. These distances are stored in a matrix A, with each row and column corresponding to a point on the manifold and each entry Ai,j representing the distance between the points i and j. This matrix is often very large and its memory usage grows quadratically with the addition of points–A manifold with 1,000 points of interest would need a matrix with 1,000,000 entries to store all geodesic distances.
The issue of memory use is accompanied by the issue of computation time. The most common algorithm for computing an all-pairs geodesics distance matrix consists of running the algorithm to compute single-source geodesic distances once for each point. To address these difficulties, we developed a method to estimate the all-pairs geodesics matrix using stochastic descent, taking advantage of the fact that the matrix is real and symmetric.
Every real symmetric matrix has real eigenvalues and real orthogonal eigenvectors. This means we can decompose an n × n symmetric matrix A into the product A = QΛQT where Q is the matrix whose columns are eigenvectors of A, is the matrix whose diagonal entries are the eigenvalues of A and all other entries are 0, and QT denotes the transpose of Q. Larger eigenvalues have a greater influence on the product QΛQT, so we can choose a positive integer k and store only the k largest eigenvalues and their corresponding eigenvectors to approximate the all-pairs matrix. Storing the matrix this way, memory usage increases linearly in k and n as opposed to quadratically in n.
Rather than first computing the full matrix A and then decomposing it, we can instead compute only a subset of A and attempt to construct the truncated decomposition of the full matrix from only that subset. In particular, computing the single-source geodesic distance from marked points of interest to all other points yields a subset of the rows and columns of A. Each iteration of the single-source geodesic algorithm provides less information than the last, so it seemed reasonable to end this process early and try to estimate the matrix using this information. The proposed algorithm to do this is as follows:
- Compute the single-source geodesic distances for a subset X of the points of interest and store them as in an all-pairs matrix A.
- Randomly initialize an n × n matrix Q (where n is the number of points of interest) and a set λ of size n. Λ will be the matrix whose only nonzero entries are the elements of λ along the diagonal. (If we save memory as described above, the last n – k columns of Q and the corresponding diagonal entries of Λ are set to 0. In this case, the elements of λ must be in descending order of magnitude.)
- Until convergence:
- Randomly sample m entries Ai,j from X and the diagonal of A (these entries are known to be 0).
- For each sampled entry, compute the error |Ai,j – Σl QiλlQjT| where Qi denotes the ith column of Q and λl denotes the lth element of λ.
- Take a gradient step to update Q and λ.
The hope is that this will converge to a reasonable approximation of the all-pairs matrix. I enjoyed studying this problem and look forward to implementing and testing this algorithm on some of the meshes we have been using throughout SGI.