python - Using ARPACK solving eigenvalueproblem, but getting inconsistent results with Matlab -

python - Using ARPACK solving eigenvalueproblem, but getting inconsistent results with Matlab -

i'm new arpack, downloaded script following

import time import numpy np scipy.linalg import eigh scipy.sparse.linalg import eigs np.set_printoptions(suppress=true) n=30 rstart=0 rend=n a=np.zeros(shape=(n,n)) # first row if rstart == 0: a[0, :2] = [2, -1] rstart += 1 # lastly row if rend == n: a[n-1, -2:] = [-1, 2] rend -= 1 # other rows in range(rstart, rend): a[i, i-1:i+2] = [-1, 2, -1] a[0,8]=30 start_time = time.time() evals_large, evecs_large = eigs(a, 10, sigma=3.6766133, which='lm') print evals_large end_time=time.time()-start_time print(" elapsed time: %12f seconds " % end_time)

it solves simple eigenvalue problem (the matrix a there not symmetric, set a[0,8] 30). 3 eigenvalues closest 3.6766133 (sigma=3.6766133 in setting) according arpack results are

[ 3.68402411+0.j 3.82005897+0.j 3.51120293+0.j]

then go matlab, , solve same eigenvalue problem, , results are

4.144524409923138 + 0.000000000000000i 3.642801014184622 + 0.497479798520641i 3.642801014184622 - 0.497479798520641i 2.372392770347609 + 0.762183281789166i 2.372392770347609 - 0.762183281789166i 3.979221766266502 + 0.000000000000000i 3.918541441830947 + 0.000000000000000i 3.820058967057387 + 0.000000000000000i 3.684024113506185 + 0.000000000000000i 3.511202932803536 + 0.000000000000000i 3.307439963195127 + 0.000000000000000i 3.080265978640102 + 0.000000000000000i 2.832849552917550 + 0.000000000000000i 2.565972630556613 + 0.000000000000000i 2.283744793210587 + 0.000000000000000i 1.996972474451519 + 0.000000000000000i 0.927737801889518 + 0.670252740725955i 0.927737801889518 - 0.670252740725955i 1.714561796881689 + 0.000000000000000i -0.015193770830045 + 0.264703483268519i -0.015193770830045 - 0.264703483268519i 1.438919271663752 + 0.000000000000000i 0.019951101383019 + 0.000000000000000i 0.080534338862828 + 0.000000000000000i 0.181591307101504 + 0.000000000000000i 0.318955140475174 + 0.000000000000000i 0.488231021129767 + 0.000000000000000i 0.688030188040126 + 0.000000000000000i 1.171318650526539 + 0.000000000000000i 0.917612528393044 + 0.000000000000000i

apparently, sec mode 3.642801014184622 + 0.497479798520641i more close sigma=3.6766133, arpack didn't pick out.

what problem? help me figure out? lot.

a few things first matlab functions:

eigenvalues returned eig not sorted. in [v,d] = eig(a) guaranteed columns of v corresponding right eigenvectors eigenvalues in d(i,i). on other hand, svd returns singular values sorted in decreasing order.

d = eigs(a,k) homecoming k largest-magnitude eigenvalues. intended big , sparse matrices, , not substitute for:

class="lang-matlab prettyprint-override">d = eig(full(a)); d = sort(d, 'descend'); d = d(1:k);

(eigs based on arpack, while eig uses lapack routines).

there no natural ordering of complex numbers. convention sort function sorts complex elements first magnitude (i.e. abs(x)), phase angle on [-pi,pi] interval (i.e. angle(x)) if magnitudes equal.

matlab

with in mind, consider next matlab code:

class="lang-matlab prettyprint-override">% create same banded matrix you're using n = 30; = spdiags(ones(n,1)*[-1,2,-1], [-1 0 1], n, n); a(1,9) = 30; %a = full(a); % k eigenvalues closest sigma k = 10; sigma = 3.6766133; d = eigs(a, k, sigma); % lets check indeed sorted distance sigma dist = abs(d-sigma); issorted(dist)

i get:

>> d d = 3.684024113506185 + 0.000000000000000i 3.820058967057386 + 0.000000000000000i 3.511202932803535 + 0.000000000000000i 3.918541441830945 + 0.000000000000000i 3.979221766266508 + 0.000000000000000i 3.307439963195125 + 0.000000000000000i 4.144524409923134 + 0.000000000000000i 3.642801014184618 + 0.497479798520640i 3.642801014184618 - 0.497479798520640i 3.080265978640096 + 0.000000000000000i >> dist dist = 0.007410813506185 0.143445667057386 0.165410367196465 0.241928141830945 0.302608466266508 0.369173336804875 0.467911109923134 0.498627536953383 0.498627536953383 0.596347321359904

you can seek similar results using dense eig:

class="lang-matlab prettyprint-override">% closest k eigenvalues sigma ev = eig(full(a)); [~,idx] = sort(ev - sigma); ev = ev(idx(1:k)) % compare against eigs norm(d - ev)

the difference acceptably little (close machine epsilon):

>> norm(ev-d) ans = 1.257079405021441e-14 python

similarly in python:

class="lang-py prettyprint-override">import numpy np scipy.sparse import spdiags scipy.sparse.linalg import eigs # create banded matrix n = 30 = spdiags((np.ones((n,1))*[-1,2,-1]).t, [-1,0,1], n, n).todense() a[0,8] = 30 # eigs: k closest eigenvalues sigma k = 10 sigma = 3.6766133 d = eigs(a, k, sigma=sigma, which='lm', return_eigenvectors=false) d = d[::-1] x in d: print '{:.16f}'.format(x) # eig ev,_ = np.linalg.eig(a) idx = np.argsort(np.abs(ev - sigma)) ev = ev[idx[:k]] x in ev: print '{:.16f}'.format(x)

with similar results:

# eigs 3.6840241135061853+0.0000000000000000j 3.8200589670573866+0.0000000000000000j 3.5112029328035343+0.0000000000000000j 3.9185414418309441+0.0000000000000000j 3.9792217662665070+0.0000000000000000j 3.3074399631951246+0.0000000000000000j 4.1445244099231351+0.0000000000000000j 3.6428010141846170+0.4974797985206380j 3.6428010141846170-0.4974797985206380j 3.0802659786400950+0.0000000000000000j # eig 3.6840241135061880+0.0000000000000000j 3.8200589670573906+0.0000000000000000j 3.5112029328035339+0.0000000000000000j 3.9185414418309468+0.0000000000000000j 3.9792217662665008+0.0000000000000000j 3.3074399631951201+0.0000000000000000j 4.1445244099231271+0.0000000000000000j 3.6428010141846201+0.4974797985206384j 3.6428010141846201-0.4974797985206384j 3.0802659786400906+0.0000000000000000j

python matlab scipy eigenvalue arpack