Basic questions regarding slepc4py/petsc4py

Question

I am brand new to slepc4py and petsc4py (and SLEPc and PETSc in general), and I have a few basic questions. I apologize if these have been answered somewhere - I've looked around for a while and haven't had much luck. If you know the answer to one (or more) of these questions, I would very much appreciate your insight. Alternatively, if you could point me to some resource that I should consult, I would appreciate that as well.

I'm trying to implement spectral clustering for a relatively large data set using petsc4py and slepc4py. The steps involved are the following (paraphrasing Ulrike von Luxburg's Tutorial on Spectral Clustering):

1. Supposing that we have a similarity graph represented by a weighted adjacency matrix W, compute the unnormalized graph Laplacian, L = D - W, where D is a diagonal matrix containing the degrees of each node.
2. Compute the first k eigenvectors of L. Let U be the matrix containing these vectors as columns.
3. Each row of U now represents a point in k-dimensional space. Cluster these points using the k-means algorithm.

My questions are as follows:

1. What is the best way to assemble a PETSc sparse matrix? I have a file that contains the nonzero entries of my matrix stored as 3-tuples (i, j, val). I tried modeling the matrix assembly after the example in ex1.py that comes with slepc4py, by doing the following:

A = PETSc.Mat()
A.create()
A.setSizes([max_i, max_j])
A.setUp()
for line in file:
    # split line, etc
    A[i,j] = val
A.assemble()

This took an unreasonably long time to do. I had more success creating a scipy sparse array from the file and then using PETSc.Mat().createAIJ(scipy_mat) to create a PETSc matrix from the scipy matrix. Am I approaching this the right way?

2. What is the best way to compute the graph Laplacian? On a related note, once I have the weighted adjacency matrix loaded in, it is not clear how I should assemble the graph Laplacian. In Matlab or scipy I would iterate over the rows, summing the weights for the edges in a row to get the diagonal element, etc. Is that kind of thinking appropriate with PETSc matrices, or should I be going about this in a different way?

3. How can I create the U matrix from the eigenvectors? Following the examples that come with slepc4py, I was able to compute the first k eigenvalues/eigenvectors for my sample matrix. The example shows how to print out the eigenvalues. It is less clear how to use the eigenvectors after the eigenproblem has been solved. In Matlab or scipy creating a matrix with the eigenvectors as columns would be straightforward. Is it possible to stick a bunch of PETSc Vector objects together to make a PETSc Matrix? I have not found anything to suggest that it is, but I also haven't been able to find any other good way of doing this.

Thanks in advance for any help you might be able to provide!

Cheers, -Rob

Answer 1

When you're setting up a matrix in PETSc, it's best to preallocate where the non-zero entries of the matrix are going to go before you start filling in the values. It looks like you already know the number of non-zero entries, so you can use createAIJ and assemble first to build the matrix structure, then fill the entries.

Internally, PETSc uses a sparse matrix storage scheme, which has the advantage that matrix-vector multiplication is extremely fast, but changing the structure of the matrix once you've already defined it is quite slow. So each time you add a new entry, it has to deallocate the matrix structure, reallocate it to accommodate the new value, then insert it.