Specifying initial vector for finding Eigenvectors using Arnoldi method

Question

I am trying to speedup the calculation of eigenvalues, given that I have good guesses for the eigenvectors. From what I know of Arnoldi/Lanczos, my good guesses should be helpful. Unfortunately, I am unable to see any speedup when I give these guesses to Mathematica.

Here is a simple and very artificial example.

I first create a random vector:

testvec = SparseArray[Table[RandomInteger[2000] -> Random[], {500}],{2000}]

I then create a matrix for which this is an eigenvector:

testmat = KroneckerProduct[testvec, testvec]

By construction the eigenvalue is

ev=(testvec.testmat.testvec)/(testvec.testvec)

On my computer (a macbook air running Mathematica 10), the following takes about 0.24 seconds:

AbsoluteTiming[Eigenvalues[testmat, 1, Method -> "Arnoldi"]]

Indeed it finds the right eigenvalue. I now feed it information about the eigenvector:

AbsoluteTiming[Eigenvalues[testmat, 1,Method -> {"Arnoldi", "StartingVector" -> testvec}]]

This once again takes 0.24 seconds.

1. Why does this not give me a speedup?
2. What can I do differently to get a speedup?

Note: I realize that I can get about a factor of 4 speedup by using the "Shift" option -- but it seems I should be able to get a couple orders of magnitude better by using a good starting vector.

Answer 1

Since the test matrix you use has rank 1, many iterative schemes will converge very rapidly. Starting with a more general real symmetric, positive-definite matrix, we compare 3 ways of computing the largest eigenvalue (and corresponding eigenvector)

a = RandomReal[NormalDistribution[0, 1], {2000, 2000}];
testmat = Transpose[a].a;

A generic approach

AbsoluteTiming[{val, vec} = Eigensystem[testmat, 1];]
(* {0.893087, Null} *)

The uninitialised Arnoldi method is 5 x faster

AbsoluteTiming[{val1, vec1} = 
   Eigensystem[testmat, 1, Method -> "Arnoldi"];]
(* {0.162465, Null} *)

The initialised Arnoldi method is almost 4 x faster again

AbsoluteTiming[{val2, vec2} = 
   Eigensystem[testmat, 1, 
    Method -> {"Arnoldi", "StartingVector" -> First[vec]}];]
(* {0.046821, Null} *)

All give the same eigenvalue

Flatten[{val, val1, val2}]
(* {7935.59, 7935.59, 7935.59} *)

EDIT

I also note that your examples evaluate much faster if converted to full, rather than sparse arrays.