[SciPy-User] Eigenvectors of sparse symmetric matrix

Lutz Maibaum lutz.maibaum@gmail....
Mon Oct 25 19:39:50 CDT 2010


On Oct 25, 2010, at 3:27 PM, Pauli Virtanen wrote:
> Mon, 25 Oct 2010 13:35:21 -0700, Lutz Maibaum wrote:
> [clip]
>> I think I am seeing the same problem as Hao. I have uploaded a sample
>> matrix exhibiting this problem:
>>    http://stanford.edu/~maibaum/matrix.mtx.gz
> 
> Thanks! Not sure if I can reproduce the issue yet (it's been running an 
> hour, and the end is not in sight :)

That's odd. On my machine the 32 bit computation takes about a minute, whereas the 64 bit one takes about 15 minutes before it terminates.

>> I think this is essentially a convergence issue. Calling the eigensolver
>> with an extra argument tol=1e-4 will give reasonable results.
> 
> The wrappers are currently written to print a warning if the iteration 
> does not converge, and I do not see yet why that doesn't happen.

Interesting, I don't get a warning. I see that there is a convergence check in _SymmetricArpackParams.iterate(), but not in _UnsymmetricArpackParams.iterate(). Could that be the reason?

>> 1. The sparse eigensolvers take a "maxiter" argument, which defaults to
>> None. I would have expected this to mean that there is no cap on the
>> amount of iterations. Looking at the code, it seems like it is set to 10
>> times the size of the matrix by default. Is this correct?
> 
> Correct. ARPACK wants some limit, so something has to be given. I don't 
> know if 10*n is the correct choice, though.

Thanks for confirming this. It would be great if that could be added to the docs.

>> 2. The eigensolver also takes an optional tolerance argument, which
>> defaults to 0. This seems like an odd choice for an iterative numerical
>> algorithm. The ARPACK manual mentions that the default tolerance is set
>> to the machine precision, is this what a tolerance of 0 means? This
>> might explain why I don't get reasonable results using 64 bit floats,
>> because the precision is just too high.
> 
> Yes, the default 0 means machine epsilon, same as in ARPACK (where it's 
> also a 'default'). Documentation is lacking, though.
> 
>> 3. Is there any way for the calling function to know whether convergence
>> was achieved, or whether the eigensolver returned because the maxiter
>> iterations were performed? Something like a "NotConverged" exception,
>> maybe?
> 
> It will be changed to raise an exception in Scipy 0.9 on non-convergence, 
> with the obtained partial result stuffed in the exception.

That would very useful. Thanks for all your help!

  Lutz



More information about the SciPy-User mailing list