2

I have an eigenvalue problem which I would like to solve using SciPy. In this problem it is strongly preferable to use an iterative solver which will return a subset of the (largest/smallest) eigenvalues. My matrix is complex and non-symmetric and so the obvious choice seems to be scipy.sparse.linalg.eigs.

The problem is that the largest/smallest eigenvalues found do not correspond to any of those found by directly computing all eigenvalues, for example with np.linalg.eig or scipy.linalg.eig.

There are two ways I have formulated the eigenproblem, neither of which provide a correct subset of eigenvalues corresponding to some of those found when using scipy.linalg.eigs.

For this example I am trying to calculate the eigenvalues of an 8x8 matrix. There are two different matrices I can use to find the eigenvlaues, T and B, which are complex128 matrices.

print(T)
>>>[[-1.09958765e+00+3.22410005e-07j  1.00537570e+00+7.30703771e-07j
   2.12573887e+00+3.24311380e-07j -1.00543668e+00-1.20938846e-07j
   3.21561757e-12-3.21564175e-14j -7.94387165e-12+7.94365030e-14j
  -8.60430404e-12+8.60436492e-14j  2.75834876e-12-2.75818666e-14j]
 [-4.68667570e-01-7.28085494e-07j  9.92912060e-01+1.75780147e-07j
   4.60413319e-01-5.63575113e-07j  7.00842564e-03+6.19343218e-07j
   7.94387165e-12-7.94365030e-14j -8.71252269e-12+8.71260741e-14j
   2.88892676e-12-2.88876502e-14j  1.97286237e-12-1.97308898e-14j]
 [ 2.08493767e+00+3.15947000e-07j  9.94704212e-01+1.25687994e-07j
  -1.11031077e+00+3.15378798e-07j -9.94644059e-01-7.27211368e-07j
  -8.60430404e-12+8.60436492e-14j -2.88892676e-12+2.88876502e-14j
   3.27314942e-12-3.27317058e-14j  7.85942560e-12-7.85920910e-14j]
 [-4.84795013e-01+5.42421874e-07j -5.89429262e-03+6.14998457e-07j
   4.80114808e-01+7.24892136e-07j  1.00581231e+00+2.04857403e-07j
  -2.75834876e-12+2.75818666e-14j  1.97286237e-12-1.97308898e-14j
  -7.85942560e-12+7.85920910e-14j -8.77410575e-12+8.77417811e-14j]
 [ 3.87718452e+11+3.87709389e+09j  1.49301876e+08+1.31818744e+06j
  -3.92737411e+11-3.92752141e+09j -1.37929035e+08-1.31818533e+06j
  -1.09958765e+00+3.22410005e-07j  4.68667570e-01+7.28085494e-07j
   2.08493767e+00+3.15947000e-07j  4.84795013e-01-5.42421874e-07j]
 [-1.49301876e+08-1.31818744e+06j  9.03206411e+07+8.26900912e+05j
   3.97988227e+08+4.04443606e+06j -6.67845152e+07-8.26891445e+05j
  -1.00537570e+00-7.30703771e-07j  9.92912060e-01+1.75780147e-07j
  -9.94704212e-01-1.25687994e-07j -5.89429262e-03+6.14998457e-07j]
 [-3.92737411e+11-3.92752141e+09j -3.97988227e+08-4.04443606e+06j
   3.97324134e+11+3.97314939e+09j  3.86769360e+08+4.04443403e+06j
   2.12573887e+00+3.24311380e-07j -4.60413319e-01+5.63575113e-07j
  -1.11031077e+00+3.15378798e-07j -4.80114808e-01-7.24892136e-07j]
 [ 1.37929035e+08+1.31818533e+06j -6.67845152e+07-8.26891445e+05j
  -3.86769360e+08-4.04443403e+06j  9.09527203e+07+8.26901241e+05j
   1.00543668e+00+1.20938846e-07j  7.00842564e-03+6.19343218e-07j
   9.94644059e-01+7.27211368e-07j  1.00581231e+00+2.04857403e-07j]]

print(B)
>>>[[ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   1.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j]
 [ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  1.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j]
 [ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   1.        +0.00000000e+00j  0.        +0.00000000e+00j]
 [ 0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  0.        +0.00000000e+00j
   0.        +0.00000000e+00j  1.        +0.00000000e+00j]
 [-0.79928622+1.87502871e-07j -1.43762007+1.02122070e-07j
   1.14981416+2.91511693e-07j -0.78754995+5.79168236e-08j
  -3.32019757+3.58806586e-07j  2.76830601+1.40564911e-06j
   4.06144987+4.60131633e-07j -0.54325794-3.46158151e-07j]
 [ 0.05767153-3.25797920e-07j -0.77301871+1.90361027e-07j
   1.29667889-5.59059891e-08j -1.12859667-2.86954366e-07j
  -2.19932049-1.39749785e-06j  3.10159289+6.37002086e-07j
   0.87977116-7.93318428e-07j -0.20013653+1.04983798e-06j]
 [ 1.13144449+2.88366134e-07j  0.75604402-5.03575308e-08j
  -0.81383949+1.79831306e-07j  1.43905136-1.01142516e-07j
   3.97388643+4.53177472e-07j  0.58049841+3.50436705e-07j
  -3.38007844+3.37107848e-07j -2.77547348-1.40198341e-06j]
 [-1.29860856+4.74789429e-08j -1.15176544-2.92841003e-07j
  -0.01185552+3.30116757e-07j -0.84010699+1.76973150e-07j
  -0.94093613+7.63127061e-07j -0.18462657+1.05164339e-06j
   2.29381822+1.40564792e-06j  3.17633502+7.03936184e-07j]]

Now I can get values through direct solution

vals, vecs = scipy.linalg.eig(T)
print(vals)

>>>[-6.26066826+2.77160294e-07j -0.15972736-7.07114310e-09j
  1.22206543+2.81763275e-01j  1.22222891-2.81656470e-01j
  0.77691915+1.79037089e-01j  0.77698287-1.79143629e-01j
  0.99998715+1.22009182e-02j  0.999864  -1.21994156e-02j]

With the same result for scipy.linalg.eig(B), numpy.linalg.eig(T), and numpy.linalg.eig(B)

Using scipy.sparse.linalg.eigs, however, produces

vals, vecs = scipy.sparse.linalg.eigs(T, k=2)
print(vals)
>>>[  93.65239733+1529.24383301j -100.07923265-1529.25018015j]

and

vals, vecs = scipy.sparse.linalg.eigs(B, k=2)
print(vals)
>>>[-6.26066826+2.77160294e-07j  1.22222891-2.81656470e-01j]

What is source of this discrepancy?

astnstn
  • 21
  • 1
  • 1
    Your `T` matrix is very numerically "difficult". You have values spanning many orders of magnitude. Of course this is not a good excuse and the condition number is not terrible but it seems to be the cause of the numerical instabilities. For example dividing `T` by the max abs value makes sparse and regular `eigs` agree to within the tolerance. Of course when multiplying back out the tolerance blows up.... – Simon Tartakovksy May 31 '23 at 18:59

0 Answers0