Eigenvalue / Eigenvector Calculation

Question

I'm currently trying to compute the three smallest eigenvectors for a 34 by 34 matrix. While I was expecting this to take some time, Mathematica has been running for the past 3 hours, which seems maybe slightly disproportional to the task at hand. The matrix is a laplacian matrix, so it is relatively sparse. Here is my input (G is the graph itself):

A1 = AdjacencyMatrix[G];
D1 = DiagonalMatrix[VertexDegree[G]];
NL1 = IdentityMatrix[34] - MatrixPower[D1, -1/2].A1.MatrixPower[D1, -1/2];
EvecsNL1 = Eigenvectors[NL1,-3]

Any help is much appreciated.

Edit: Sorry about that:

G = ExampleData[{"NetworkGraph", "ZacharysKarateClub"}]

I'm trying to get accustomed to working with graph laplacians in mathematica, so I'm messing around with some included datasets.

Edit 2: Wow, I feel foolish. Thank you for your help.

Answer 1

For speeding up Mathematica code, a little analysis of a problem often goes a long way.

Analysis

Writing $\mathbb{G}$ for the diagonal matrix of vertex degrees and $\mathbb{A}$ for the adjacency matrix, this question seeks eigenvectors of $\mathbb{1} - \mathbb{G}^{-1/2} \mathbb{A} \mathbb{G}^{-1/2}$. That is, it looks for solutions $(\mathbf{z}, \lambda)$ to

$$\mathbf{z} \left(\mathbb{1} - \mathbb{G}^{-1/2} \mathbb{A} \mathbb{G}^{-1/2} \right) = \lambda \mathbf{z}.$$

The square roots ought to give us pause: they are likely a significant complication in performing the arithmetic.

It is straightforward to rewrite this system of equations, though, as

$$\left(\mathbf{z} \mathbb{G}^{1/2}\right)\left(\mathbb{G^{-1} A}\right) = \left(1-\lambda\right)\left(\mathbf{z}\mathbb{G}^{1/2}\right).$$

Whence, if we solve

$$\mathbf{w}\left(\mathbb{G^{-1} A}\right) = \mu \mathbf{w}$$

then $\lambda = 1-\mu$ and $\mathbf{z} = \mathbf{w} \mathbb{G}^{-1/2}$. The square roots come in only at the end as easy matrix multiplications, while $\mathbb{GA}$ is a rational matrix. Hopefully, calculations involving it will be faster.

Solution

Here is an implementation:

f[g_, a_] := Module[{e, v},
   {e, v} = Eigensystem[DiagonalMatrix[1/g] . a];
   {1 - e, v.DiagonalMatrix[Sqrt[g]]}
   ];

Let's carry it out on the example data:

g = ExampleData[{"NetworkGraph", "ZacharysKarateClub"}];
{l, z} = f[VertexDegree[g], AdjacencyMatrix[g]]; // AbsoluteTiming

{31.8538219, Null}

---

Remarks

A half minute is reasonable (the output is huge and complicated). The results agree with the (much faster) numerical solution to about five significant figures--but no more. Thus, if higher accuracy is desired, this approach makes that feasible. Typically, you would find the exact eigensystem and then immediately convert that to floating-point representation for all the (much less numerically involved) post-processing operations.

This implementation did not bother to assure that the eigenvectors are still of unit length. If desired, do this after the fact by applying #/Norm[#]& to them. Also, it returns all the eigenvectors and eigenvalues, rather than the three corresponding to the smallest absolute values of the eigenvalues. But finding those three requires only a simple post-check of the list of eigenvalues. This is made necessary because we are really looking for the eigenvectors for which $\mu$ is as close to $1$ as possible and there's no savings associated with computing just them.

The output of this solution and that of a direct attack on the numerical version of the original problem may differ in two ways. First, the eigenvalues will occur in a different order. Second, even when the eigenvectors are normalized, they are determined only up to sign. These differences are of no mathematical or practical consequence.

Finally, Eigensystem arranges to return a set of mutually orthonormal eigenvectors for any eigenvalue of multiplicity greater than one. The modified approach described here may ruin their orthonormality. That property can be restored, if desired, by applying a Gram-Schmidt process to each degenerate eigenspace.

Answer 2

Note that your A1 is a SparseArray, while your NL1 is a full list of lists, since you've formed it using IdentityMatrix and DiagonalMatrix which return full lists of lists. This could have a huge impact on memory footprint as well as speed and accuracy of your computations, particularly if you're working with very large, sparse graphs. I would take this into account when forming the Laplacian matrix.

As I understand it, the Laplacian matrix of a graph is its adacency matrix minus its vertex degree matrix. We could implement this as a sparse matrix like so

g = GraphData[{"Kneser", {14, 6}}];
vertexDegrees = VertexDegree[g];
diagonalRules = Table[{i, i} -> -vertexDegrees[[i]],
  {i, 1, Length[vertexDegrees]}];
adjacencyRules = ArrayRules[AdjacencyMatrix[g]];
laplacian = SparseArray[Join[diagonalRules, adjacencyRules]];
Dimensions[laplacian]

(* Out: {3003, 3003} *)

It's a large matrix; let's look at just a piece of it.

MatrixForm[Normal[laplacian[[1282 ;; 1296, 1282 ;; 1296]]] /. 
  0 -> Style[0, LightGray]]

I guess your version of the command to compute the 10 smallest eigenvalues would be

Eigenvalues[N[laplacian], -10] // AbsoluteTiming

(* Out: {2.454553, {-22.1063, -22., -22., -21.8781, -13.0058, 
   -13., -13., -13., -12.9948, 4.5206*10^-15}} *)

You might try something like so:

Eigenvalues[N[laplacian], 10, Method -> {"Arnoldi",
  "Shift" -> 0, "BasisSize" -> 50}] // AbsoluteTiming

(* Out: {1.691656, {-13., -13., -13., -13., -13., -13., -13.,
   -13., -12.6622, -6.90978*10^-16}} *)

You can use the Eigensystem command to check for accuracy.

{vals, vecs} = Eigensystem[N[laplacian], 10, Method -> {"Arnoldi",
  "Shift" -> 0, "BasisSize" -> 50}];
Norm[laplacian.vecs[[1]] - vals[[1]]*vecs[[1]]]

(* Out: 2.30196*10^-9 *)

{vals, vecs} = Eigensystem[N[laplacian], -10];
Norm[laplacian.vecs[[1]] - vals[[1]]*vecs[[1]]]

(* 0.284252 *)

It appears that the Arnoldi specified version is faster and more accurate. I don't know a lot of the details behind the computations, particularly the second syntax with the -10. Perhaps @DanielLichtblau or @ruebenko could comment on that.