0

I need to invert a large number (currently 1e6, could maybe be optimized to 3e3) of symmetric complex tridiagonal matrices. The matrices M are all very similar. In fact M is a function of two parameters M(a, b), and I need to calculate the inverse on a grid.

What I currently do is 

import numpy as np
from numpy.linalg import inv

N = 40
mat_indices = range(0, N)

a = np.linspace(-5, 5, dtype=np.complex256, num=500)                      
b = np.linspace(-5, 5, dtype=np.complex256, num=500)                    
w += 5j * max((w[1] - w[0], eps[1] - eps[0]))                                   

M = np.zeros(a.shape + b.shape + (N, N), dtype=np.complex)           
for lay in layers:                                                              
    M[..., lay, lay] = fa(a[:, np.newaxis]) + fb(b[np.newaxis, :])
for lay in layers[:-1]:                                                         
    G_latt_inv[..., lay, lay+1] = fsub(a, b)                            
    G_latt_inv[..., lay+1, lay] = fsub(a, b)                                    

M_inv = inv(M)

Obviously this is a terrible solution, I neither use symmetry nor sparsity of the matrix. Furthermore, an iterative method should be advantageous, as my matrices are all very similar.

I know there is scipy.sparse, probably I could use a solve method for b=Unity. But I still don't use any of my many advantageous properties.

I do not need an optimal solution, I rather prefer simplicity because I don't have a lot of time to test and verify my results. But currently it is way too inefficient.

DerWeh
  • 1,721
  • 1
  • 15
  • 26
  • A suggestion for using symmetry is [here](https://stackoverflow.com/a/40710036). Unfortunately, there is no magic "go faster" button and I don't have a lot of time to test and verify your results, either. –  Dec 22 '17 at 03:31
  • 1
    Of course there is no magic. But plenty of linear algebra you can use to speed things up. I tried also solving it with `scipy.linalg.solveh_banded`. The tradeoff is that I need a double loop to use it. depending on the dimensions one or the other is more suitable. But for the relevant ranges the difference is only a factor of 2. – DerWeh Dec 22 '17 at 13:13
  • I just discovered that [scipy.linalg.solve](https://docs.scipy.org/doc/scipy-1.0.0/reference/generated/scipy.linalg.solve.html) allows to use `assume_a=sym` to assume a symmetric matrix. Can this maybe also done for the `solve_banded`? – DerWeh Dec 27 '17 at 09:58
  • [Tridiagonal solver numpy](https://gist.github.com/ofan666/1875903) should work right off for complex too. ([Tridiagonal solver cython](https://github.com/cpcloud/PyTDMA) ? [complex-numbers-in-cython](https://stackoverflow.com/questions/30054019/complex-numbers-in-cython) look messy, at least gcc.) – denis Jan 23 '18 at 14:27

0 Answers0