[SciPy-User] Eigenvectors of sparse symmetric matrix
Hao Xiong
I am trying to compute the eigenvectors corresponding to the d+1 smallest
eigenvalues of A=W.T*W. I started with W as a dense matrix and then
W = sparse.csr_matrix(W)
A = W.dot(W) # W.T * W
W,V = eigen_symmetric(A,d+1, which='SM')
The biggest problem is that the algorithm fails to converge and I get
all zeros as eigenvectors for a testing dataset. Using dense SVD I
got the expected results.
The second problem is that this sparse version is much slower than
the dense version as
u,s,vh = svd(W)
The testing data only has 1000x1000, while I expect the real data
will have millions by millions of entries. Each row will have only
a dozen to at most dozes of non-zero entries.
Thanks,
Hao
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Lutz Maibaum
> I am trying to compute the eigenvectors corresponding to the d+1 smallest
> eigenvalues of A=W.T*W. I started with W as a dense matrix and then
> W = sparse.csr_matrix(W)
> A = W.dot(W) # W.T * W
> W,V = eigen_symmetric(A,d+1, which='SM')
>
> The biggest problem is that the algorithm fails to converge and I get
> all zeros as eigenvectors for a testing dataset. Using dense SVD I
> got the expected results.
What operating system are you using? The sparse eigensolvers have some issues on 64-bit OS-X (see http://projects.scipy.org/scipy/ticket/1220).
I am currently having similar issues. Oddly enough, I get reasonable results if I convert the matrix from the native float64 to the float32 data type. If you are on a 64 bit system, could you try that and let me know if it changes anything? For example, define A as
A = W.dot(W).astype(numpy.float32)
Best,
Lutz
Pauli Virtanen
> I am trying to compute the eigenvectors corresponding to the d+1
> smallest eigenvalues of A=W.T*W. I started with W as a dense matrix and
> then
> W = sparse.csr_matrix(W)
> A = W.dot(W) # W.T * W
That is W*W and not (W.T)*W
> W,V = eigen_symmetric(A,d+1, which='SM')
>
> The biggest problem is that the algorithm fails to converge and I get
> all zeros as eigenvectors for a testing dataset. Using dense SVD I got
> the expected results.
You can try playing with setting the maxiter parameter to allow ARPACK to
spend more iterations on the problem.
--
Pauli Virtanen
Hao Xiong
Don't know if this matters but the data originally are integers, 0, 1,
and 2.
They are the encoding of genotypes.
Regards,
Hao
Hao Xiong
>> I am trying to compute the eigenvectors corresponding to the d+1
>> smallest eigenvalues of A=W.T*W. I started with W as a dense matrix and
>> then
>> W = sparse.csr_matrix(W)
>> A = W.dot(W) # W.T * W
> That is W*W and not (W.T)*W
Thanks, Pauli. Somehow I convinced myself it was otherwise. I have corrected that.
>> W,V = eigen_symmetric(A,d+1, which='SM')
>>
>> The biggest problem is that the algorithm fails to converge and I get
>> all zeros as eigenvectors for a testing dataset. Using dense SVD I got
>> the expected results.
> You can try playing with setting the maxiter parameter to allow ARPACK to
> spend more iterations on the problem.
I tried maxiter=100000 and still got zero vectors. I must be missing something.
Regards,
Hao
josef...@gmail.com
>
>
>>> I am trying to compute the eigenvectors corresponding to the d+1
>>> smallest eigenvalues of A=W.T*W. I started with W as a dense matrix and
>>> then
>>> W = sparse.csr_matrix(W)
>>> A = W.dot(W) # W.T * W
>
>> That is W*W and not (W.T)*W
>
> Thanks, Pauli. Somehow I convinced myself it was otherwise. I have corrected that.
>
>>> W,V = eigen_symmetric(A,d+1, which='SM')
>>>
>>> The biggest problem is that the algorithm fails to converge and I get
>>> all zeros as eigenvectors for a testing dataset. Using dense SVD I got
>>> the expected results.
>
>> You can try playing with setting the maxiter parameter to allow ARPACK to
>> spend more iterations on the problem.
>
> I tried maxiter=100000 and still got zero vectors. I must be missing something.
just a weird idea, since I have no idea what eigen_symmetric is doing,
and there are no docs that I have seen for the extra options:
Is it possible to run a dense svd on a (random) subset of the data and
then use those as starting values for the sparse decompositions?
Josef
Pauli Virtanen
[clip]
>> You can try playing with setting the maxiter parameter to allow ARPACK
>> to spend more iterations on the problem.
>
> I tried maxiter=100000 and still got zero vectors. I must be missing
> something.
Can you perhaps write a small routine that generates such a fake matrix,
or save it to a file and put one online. This might help in finding out
why the routine doesn't behave as expected.
--
Pauli Virtanen
Lutz Maibaum
> Mon, 25 Oct 2010 11:14:16 -0700, Hao Xiong wrote:
>>> You can try playing with setting the maxiter parameter to allow ARPACK
>>> to spend more iterations on the problem.
>>
>> I tried maxiter=100000 and still got zero vectors. I must be missing
>> something.
>
> Can you perhaps write a small routine that generates such a fake matrix,
> or save it to a file and put one online. This might help in finding out
> why the routine doesn't behave as expected.
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
A simple test to illustate the problem is the following:
In [1]: from scipy.io import mmread
In [2]: import scipy.sparse.linalg.eigen
In [3]: import numpy as np
In [4]: a = mmread("matrix.mtx.gz")
In [5]: scipy.sparse.linalg.eigen(a.astype(np.float32))[0]
Out[5]:
array([ 1.00000119 +0.00000000e+00j, 0.99999952 +9.27230701e-07j,
0.99999952 -9.27230701e-07j, 1.00000107 +0.00000000e+00j,
1.00000715 +0.00000000e+00j, 1.00001323 +0.00000000e+00j])
In [6]: scipy.sparse.linalg.eigen(a.astype(np.float64))[0]
Out[6]:
array([ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])
Note that the 64-bit version takes about 10 times longer than the 32 bit version, and returns only zeros for the eigenvalues. This happens on NumPy 1.5.0, SciPy 0.8.0, Python 2.6.6 on 64-bit OS-X, where float64 is the default type of floats.
I think this is essentially a convergence issue. Calling the eigensolver with an extra argument tol=1e-4 will give reasonable results. While debugging this, I came up with a few questions that the docs don't quite answer, and it would be great if someone could help me out:
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?
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.
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?
Hao, could you try to call the eigensolver with a generous tolerance argument, and see if this fixes your problem? If your issue is unrelated, I apologize for hijacking this thread.
Thanks,
Lutz
Pauli Virtanen
[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 :)
[clip]
> 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.
[clip]
> 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.
> 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.
As it is now, it only raises a warning if it does not convergence. But
no, there is currently no good programmatic way of catching this.
--
Pauli Virtanen
Lutz Maibaum
> 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
Lutz Maibaum
> On Oct 25, 2010, at 3:27 PM, Pauli Virtanen wrote:
>> 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?
The test is actually there, the symmetric case just has an additional test.
However, I noticed that the warning is triggered when the returned "info" status is -1. Is it possible that it should be compared to 1? From http://www.caam.rice.edu/software/ARPACK/UG/node137.html:
c INFO Integer. (INPUT/OUTPUT)
c If INFO .EQ. 0, a randomly initial residual vector is used.
c If INFO .NE. 0, RESID contains the initial residual vector,
c possibly from a previous run.
c Error flag on output.
c = 0: Normal exit.
c = 1: Maximum number of iterations taken.
c All possible eigenvalues of OP has been found. IPARAM(5)
c returns the number of wanted converged Ritz values.
c = 2: No longer an informational error. Deprecated starting
c with release 2 of ARPACK.
c = 3: No shifts could be applied during a cycle of the
c Implicitly restarted Arnoldi iteration. One possibility
c is to increase the size of NCV relative to NEV.
c See remark 4 below.
c = -1: N must be positive.
c = -2: NEV must be positive.
c = -3: NCV-NEV >= 2 and less than or equal to N.
c = -4: The maximum number of Arnoldi update iteration
c must be greater than zero.
c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
c = -6: BMAT must be one of 'I' or 'G'.
c = -7: Length of private work array is not sufficient.
c = -8: Error return from LAPACK eigenvalue calculation;
c = -9: Starting vector is zero.
c = -10: IPARAM(7) must be 1,2,3,4.
c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
c = -12: IPARAM(1) must be equal to 0 or 1.
c = -9999: Could not build an Arnoldi factorization.
c IPARAM(5) returns the size of the current Arnoldi
c factorization.
Thanks,
Anne Archibald
> Sun, 24 Oct 2010 18:42:14 -0700, Hao Xiong wrote:
>> I am trying to compute the eigenvectors corresponding to the d+1
>> smallest eigenvalues of A=W.T*W. I started with W as a dense matrix and
>> then
>> W = sparse.csr_matrix(W)
>> A = W.dot(W) # W.T * W
>
> That is W*W and not (W.T)*W
Just quickly: is your matrix really symmetric? This could lead to all
sorts of bizarre failures to converge.
Anne
Hao Xiong
I think our problems are related but with different symptoms.
First, changing to float32 did not help; I got the same errors.
Second, changing tolerance to 1e-4, 1e-5, 1e-15, all quash warning but
do not solve the problem: all zero eigenvalues and eigenvectors.
I have uploaded a 100x100 symmetric matrix at
http://www.epibiostat.ucsf.edu/biostat/XiongH/sym_sp.zip
On my machine, eigen_symmetric(A,3, which='SM', maxiter=10000)
results in all zero eigenvalues and eigenvectors, if A is the uploaded
matrix read using mmread.
However, converting it to dense and calling np.linalg.eig results in
correct answer.
One way to check this is that the eigenvector corresponding to the least
eigenvalue, which should be 0, is a multiple of a vector of 1's.
I am wondering if a sensible starting value is needed. Iterative methods
often need a reasonable starting value.
Best,
Hao
Hao Xiong
I have a Debian machine running testing version but with
a custom compiled Python-2.7 and NumPy-1.5.1RC1
and SciPy-0.8.0. I specified MKL when I compiled
NumPy and Scipy and I passed -O3 to GCC. NumPy's
basic testing finishes successively, but SciPy's crashes
with a segfault.
On a Gentoo machine with Python-2.6 and NumPy-1.5
and SciPy-0.8.0 linked against Atlas library, I can get
NumPy's basic test finish successively, but SciPy's fails
with one failure.
Hope this helps.
Thanks,
Hao
Lutz Maibaum
> Second, changing tolerance to 1e-4, 1e-5, 1e-15, all quash warning but
> do not solve the problem: all zero eigenvalues and eigenvectors.
Interesting. Did you get a warning about non-convergence before? I'm
not sure what's going on with the 1e-15 tolerance, but the other ones
are probably too large because the smallest eigenvalues seem to be
very close to zero. For your matrix, I get
In [28]: scipy.sparse.linalg.eigen_symmetric(a,3, which='SM',
tol=1e-12, maxiter=10000000)
Out[28]:
(array([ 2.13492457e-17, 9.54397558e-07, 1.77823892e-06]),
array([[-0.1 , 0.05414805, -0.03697394],
[-0.1 , 0.04864488, 0.12934357],
[-0.1 , 0.09785515, -0.02710373],
[-0.1 , 0.12079696, -0.095882 ],
[-0.1 , -0.0923547 , 0.06131896],
[-0.1 , 0.08566995, 0.02677637],
[-0.1 , -0.00999687, -0.03477074],
[-0.1 , 0.10420394, 0.0904605 ],
[-0.1 , -0.01824658, 0.12656016],
[-0.1 , 0.03850215, -0.16266434],
[-0.1 , 0.10411247, 0.12564386],
[-0.1 , 0.0847864 , -0.08006353],
[-0.1 , 0.06881941, -0.03651171],
[-0.1 , 0.02100945, 0.12596802],
[-0.1 , 0.00698142, -0.14227155],
[-0.1 , -0.0994901 , -0.01865521],
[-0.1 , 0.05150979, -0.13437869],
[-0.1 , 0.12983085, 0.10247783],
[-0.1 , 0.20811634, 0.04037155],
[-0.1 , 0.14284242, 0.07440172],
[-0.1 , 0.08759283, 0.00897286],
[-0.1 , 0.11652933, 0.11583934],
[-0.1 , 0.08273911, -0.12928089],
[-0.1 , 0.15103551, 0.08544608],
[-0.1 , -0.10887856, -0.03683742],
[-0.1 , 0.08946787, 0.01810116],
[-0.1 , -0.21466925, 0.08808048],
[-0.1 , 0.01112506, 0.11875543],
[-0.1 , 0.03862264, -0.03816272],
[-0.1 , -0.08819346, 0.0469191 ],
[-0.1 , -0.08715582, -0.10397484],
[-0.1 , 0.09957673, 0.12540574],
[-0.1 , -0.10165562, 0.10154619],
[-0.1 , -0.02138075, 0.06997714],
[-0.1 , -0.02087899, -0.04523328],
[-0.1 , 0.07205966, 0.00801408],
[-0.1 , 0.06474043, 0.00830429],
[-0.1 , 0.08648864, -0.00438077],
[-0.1 , 0.09298343, 0.04886763],
[-0.1 , 0.07158097, 0.0782138 ],
[-0.1 , 0.01239778, -0.15765419],
[-0.1 , -0.05888361, 0.03320853],
[-0.1 , 0.08010641, 0.08588525],
[-0.1 , 0.03127534, -0.15888655],
[-0.1 , 0.15375382, -0.00072328],
[-0.1 , 0.1309185 , 0.01948518],
[-0.1 , -0.21072633, 0.05625481],
[-0.1 , 0.00123581, -0.19868411],
[-0.1 , -0.04948594, -0.1179604 ],
[-0.1 , 0.03724257, -0.18880828],
[-0.1 , -0.05376647, 0.11361879],
[-0.1 , 0.05143578, -0.11411724],
[-0.1 , -0.04570302, 0.13384669],
[-0.1 , -0.05617232, -0.09347502],
[-0.1 , -0.20512585, 0.0484587 ],
[-0.1 , -0.027912 , -0.1848302 ],
[-0.1 , 0.14621125, -0.00988872],
[-0.1 , -0.10030626, -0.09077817],
[-0.1 , -0.0363287 , 0.02784762],
[-0.1 , -0.21623947, 0.06780762],
[-0.1 , -0.06138235, -0.1349 ],
[-0.1 , -0.09814152, -0.04398249],
[-0.1 , 0.12720599, 0.00705402],
[-0.1 , -0.01507454, -0.18508998],
[-0.1 , -0.00798772, -0.23027451],
[-0.1 , 0.0084914 , -0.14105232],
[-0.1 , -0.00326878, 0.1905542 ],
[-0.1 , -0.11332749, -0.01003244],
[-0.1 , 0.16373491, 0.07366324],
[-0.1 , -0.15722344, 0.05073253],
[-0.1 , 0.04282908, 0.05747035],
[-0.1 , -0.11459224, 0.1258188 ],
[-0.1 , -0.03079556, 0.12889243],
[-0.1 , -0.06469642, 0.15025778],
[-0.1 , 0.18106343, 0.04211254],
[-0.1 , -0.11284705, 0.05143415],
[-0.1 , -0.14384552, -0.01344659],
[-0.1 , 0.00723068, -0.19844225],
[-0.1 , -0.05921825, -0.12152038],
[-0.1 , -0.20116698, 0.08917197],
[-0.1 , -0.17052863, 0.05183343],
[-0.1 , -0.01618908, -0.05137175],
[-0.1 , -0.04433511, -0.06585839],
[-0.1 , 0.03211274, 0.1278789 ],
[-0.1 , 0.12588347, -0.08004173],
[-0.1 , -0.08788273, -0.10250587],
[-0.1 , 0.01156012, -0.00283793],
[-0.1 , -0.05788733, -0.08254978],
[-0.1 , 0.08076025, -0.03895826],
[-0.1 , 0.02232925, 0.1221555 ],
[-0.1 , -0.11745394, -0.02053012],
[-0.1 , 0.01355673, 0.12304368],
[-0.1 , -0.17090312, 0.03684983],
[-0.1 , 0.21020815, 0.00314479],
[-0.1 , -0.05367038, -0.13541344],
[-0.1 , -0.11600589, -0.07075897],
[-0.1 , -0.02728704, 0.09827715],
[-0.1 , -0.09539752, 0.04939529],
[-0.1 , 0.0058363 , 0.18280319],
[-0.1 , -0.04509233, -0.0022041 ]]))
which seems in good agreement with the dense solution.
Best,
Pauli Virtanen
> On Oct 25, 2010, at 5:39 PM, Lutz Maibaum wrote:
>> On Oct 25, 2010, at 3:27 PM, Pauli Virtanen wrote:
>>> 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?
>
> The test is actually there, the symmetric case just has an additional
> test.
>
> However, I noticed that the warning is triggered when the returned
> "info" status is -1. Is it possible that it should be compared to 1?
> From http://www.caam.rice.edu/software/ARPACK/UG/node137.html:
[clip]
> c INFO Integer. (INPUT/OUTPUT)
> c If INFO .EQ. 0, a randomly initial residual vector is used.
> c If INFO .NE. 0, RESID contains the initial residual vector,
> c possibly from a previous run.
> c Error flag on output.
> c = 0: Normal exit.
> c = 1: Maximum number of iterations taken.
[clip]
Seems indeed that the check is wrong. Needs fixing.
--
Pauli Virtanen
Pauli Virtanen
> Seems indeed that the check is wrong. Needs fixing.
If you have time, please test:
http://github.com/pv/scipy-work/tree/bug/1313-arpack
Note that the function is there renamed from `eigen` to `eigs`, as the
name shadowed with the package name (plus obvious reasons).
Hao Xiong
I am beginning to suspect my problem is a compilation issue.
On my Gentoo machine after installing arpack package separately,
scipy.sparse.linalg.test() gets
Running unit tests for scipy.sparse.linalg
NumPy version 1.5.0
NumPy is installed in /usr/lib64/python2.6/site-packages/numpy
SciPy version 0.8.0
SciPy is installed in /usr/lib64/python2.6/site-packages/scipy
Python version 2.6.5 (release26-maint, Oct 24 2010, 22:13:01) [GCC 4.4.4]
nose version 0.11.4
..../usr/lib64/python2.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:259:
DeprecationWarning: scipy.sparse.linalg.dsolve.umfpack will be removed,
install scikits.umfpack instead
' install scikits.umfpack instead', DeprecationWarning )
../usr/lib64/python2.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:75:
DeprecationWarning: scipy.sparse.linalg.dsolve.umfpack will be removed,
install scikits.umfpack instead
' install scikits.umfpack instead', DeprecationWarning )
....................K...Segmentation fault
So a few questions:
1. Should I install arpack package separately? I think SciPy source
comes with arpack, so it may not be necessary;
but could it hurt?
2. What is the recommended way of supplying flags to SciPy's build
scripts? I have set my CFLAGS to
"-march=native -mtune=native -pipe -O3 -msse." But if I compile
SciPy directly with this, build
fails. Gentoo filters flags so its build system succeeds.
Thanks,
Hao
Lutz Maibaum
> If you have time, please test:
>
> http://github.com/pv/scipy-work/tree/bug/1313-arpack
I am not sure how to use github, so I just replaced the arpack.py with the new version. The 32 bit version still works, and the 64 bit version now correctly tries to raise an exception. It seems to fail in doing so, however:
Traceback (most recent call last):
File "check.py", line 38, in <module>
compute_eigenvalues(a.astype(np.float64))
File "check.py", line 22, in compute_eigenvalues
print scipy.sparse.linalg.eigs(a)[0]
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 521, in eigs
params.iterate()
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 332, in iterate
self._raise_no_convergence()
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 201, in _raise_no_convergence
ev, vec = self.extract(True)
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 358, in extract
raise ArpackError(ierr, infodict=_NEUPD_ERRORS)
File "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 121, in __init__
RuntimeError.__init__("ARPACK error %d: %s" % (info, msg))
TypeError: descriptor '__init__' requires a 'exceptions.RuntimeError' object but received a 'str'
Again, this might be because I simply replaced the arpack.py file.
Would it be possible to add to the docstrings that maxiter=None is interpreted as maxiter=10*n (size of matrix)?
> Note that the function is there renamed from `eigen` to `eigs`, as the
> name shadowed with the package name (plus obvious reasons).
That's good to know. A while ago there was some discussion about this (http://mail.scipy.org/pipermail/scipy-user/2010-August/026401.html). In SciPy 0.7, I used
scipy.sparse.linalg.eigen.arpack.eigen()
which no longer worked in SciPy 0.8, where I now use
scipy.sparse.linalg.eigen()
Does your change mean that for SciPy 0.9, I will have to use
scipy.sparse.linalg.eigen.arpack.eigs() ?
I help maintain a software package that is supposed to work with all 3 versions of SciPy. What would be a good way to handle this? To support 0.7 and 0.8, I am currently using something like
try:
import scipy.sparse.linalg.eigen.arpack as arpack # this works for 0.7 but fails for 0.8
scipy.sparse.linalg.eigen = arpack.eigen
except:
import scipy.sparse.linalg.eigen # this works for 0.8
Any suggestions would be much appreciated.
Thanks for all your help,
Lutz
David
> On Oct 26, 2010, at 9:33 AM, Pauli Virtanen wrote:
>> If you have time, please test:
>>
>> http://github.com/pv/scipy-work/tree/bug/1313-arpack
>
> I am not sure how to use github
There is a big download button so that you can download the working tree
corresponding to that revision :)
http://github.com/pv/scipy-work/archives/bug/1313-arpack
> I help maintain a software package that is supposed to work with all 3 versions of SciPy. What would be a good way to handle this? To support 0.7 and 0.8, I am currently using something like
>
> try:
> import scipy.sparse.linalg.eigen.arpack as arpack # this works for 0.7 but fails for 0.8
> scipy.sparse.linalg.eigen = arpack.eigen
> except:
> import scipy.sparse.linalg.eigen # this works for 0.8
>
> Any suggestions would be much appreciated.
I think this is a fined method, except that you should use ImportError
in the except clause to avoid hiding arbitrary exception (if you make a
typo in the try section, for example).
An alternative is to explicitly handle the versions (using
scipy.__version__). This has the advantage of being more explicit, but
potentially more fragile (making errors in version handling is
suprisingly easy in my experience),
cheers,
David