Solve eigenvalue problem using finite differences without vectorization

Question

I am interested in solving the problem $-A u = \lambda u$ on a finite differences grid on a square. In my case, the operator $A$ is of the type $-\Delta + \mu I$, where $\mu$ is there to impose some constraints. The method I used until now was to consider a finite differences grid and get $$ 4u_{i,j} - u_{i,j+1}-u_{i,j-1}-u_{i+1,j}-u_{i-1,j} + \mu_{i,j} u_{i,j}=\lambda u_{i,j}$$ Then we transform $u_{i,j}$ in a vector and write the above problem in matrix form $Au = \lambda u$. The problem I have with this is that the size of the matrix $A$ is $n \times n$, where $n$ is the number of points in the grid. Therefore, for a grid with $N^2$ points the dimension of the matrix is $N^2 \times N^2$. Therefore, the size of the matrices grows very fast with respect to the discretization.

Is it possible to solve this kind of problem without vectorization, i.e. using only matrices of size $N \times N$?

My idea is the following: If we denote $U=(u_{i,j})$ then the discrete PDE can be written as $$ (4I -S_{0,1}-S_{0,-1}-S_{1,0}-S_{-1,0})U+\mu \odot U=\lambda U,$$ where $S_{i,j}$ is the shift of the matrix with $i$ on horizontal and $j$ on vertical directions. The $\odot$ is "pointwise" product. Can this lead to a solution?

Answer 1

For sure it is possible with iterative methods, because you can apply your operator $u\mapsto Au$ cheaply while keeping $u$ in "matrix form". It is simply a double nested for cycle. Check out ARPACK, or more in general Arnoldi-(or Krylov-)type methods: in most implementations, you can use as input a 'black-box' routine that computes $Au$ from $u$, rather than an explicit matrix.

These methods are best suited to computing only a small section of the spectrum (for instance, the largest or smallest eigenvalues and their corresponding eigenvectors), and nowadays they are the standard method for solving large and sparse eigenproblems numerically.