[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