Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
Please give the definition of G, otherwise we can't help easily. –  tkott Oct 8 '12 at 14:19
3  
You can get numerical approximations pretty much instantaneously: In[8]:= Eigenvalues[N[NL1]] Out[8]= {1.71461, 1.61191, 1.58333, 1.56951, 1.49703, 1.44858, \ 1.41692, 1.3931, 1.35178, 1.26802, 1.1593, 1.10538, 1., 1., 1., 1., \ 1., 1., 1., 1., 1., 1., 0.906816, 0.864833, 0.822943, 0.770911, \ 0.739958, 0.707208, 0.648993, 0.612231, 0.387313, 0.287049, 0.132272, \ -1.49256*10^-16} Obtaining these as exact values is a very different fish. –  Daniel Lichtblau Oct 8 '12 at 14:38
6  
I once did a pseudo inverse on a 100x200 matrix with exact numbers. Mathematica ran for a week and delivered a matrix taking about 20 MB, as every entry consisted of a rational number with hundreds of digits in numerator and denominator. In cases like this, indeed, a healthy application of N will work wonders. –  Sjoerd C. de Vries Oct 8 '12 at 15:43
1  
...so, did you really need exact eigenvectors, or were you fine with the numerical approximations? –  J. M. Oct 8 '12 at 17:05

2 Answers 2

up vote 10 down vote accepted

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.

share|improve this answer
    
"...by applying #/Norm[#] & to them." - or by using the built-in function Normalize[]. –  J. M. Oct 8 '12 at 23:32
    
"For speeding up Mathematica code, a little analysis of a problem often goes a long way" an excellent point. +1 for that and +1 for the rest (if I could) –  acl Oct 11 '12 at 0:30
    
According to the docs, eigenvectors for a single eigenvalue are chosen to be linearly independent. So, for each eigenspace, I recommend running Orthogonalize on them. This has bitten me several times with regards to numerical eigenvectors, but I don't recall if it applies to analytic ones, too. –  rcollyer Oct 11 '12 at 0:50
    
@rcollyer It all depends on whether you need those eigenvectors to be orthonormal: not all applications do. Moreover, orthonormalizing them using exact arithmetic could take some time, so it's worth considering whether to first convert the exact eigenvectors to machine precision and then orthonormalizing them. –  whuber Oct 11 '12 at 1:13
    
"Thus, if higher accuracy is desired, this approach makes that feasible." Is it generally true that doing an analytical calculation first (which yields complicated expressions) and then converting to floating point numbers gives higher accuracy? I sometimes find myself with expressions which are "ill-conditioned" and give funny (presumably wrong) results when evaluated numerically. Fiddling with them to get a form which is well behaved can often be a pain. Or am I just doing something wrong? –  sebhofer Oct 11 '12 at 17:02

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]]

enter image description here

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.

share|improve this answer
    
I don't know exactly what it does. I suspect some flavor of Krylov solving to get the largest eigenvalues for the inverse matrix. –  Daniel Lichtblau Oct 8 '12 at 20:37
    
The docs seem to indicate that it uses the ARPACK library for the Arnoldi method under the hood. Saad's book has a good discussion on Arnoldi and other Krylov subspace methods. –  J. M. Oct 8 '12 at 23:31

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.