## 3 Direct Numerical Methods for Strongly Disordered Systems
### Diagonalization Methods  
Diagonalization methods are designed to solve the matrix problem, \(H\psi=E\psi\), directly. A full matrix diagonalization scales with the third power of the matrix size. Therefore, practical calculations are often limited to matrix sizes of the order of ten thousand. For the study of the localization transition, we are usually interested in the states close to the Fermi level. Indeed, most of the numerical studies of the Anderson model are focused on the energy at the band center [41]. Methods have been proposed for calculating the eigenvalues and eigenvectors for sparse matrices in the vicinity of a target eigenvalue, \(\sigma\). Particularly, the Lanczos [140] and Arnoldi [141] methods have been widely used for strongly correlated systems [142; 143; 144]. The feature common to these methods is the Krylov subspace, \(K\), generated by repeatedly multiplying a matrix, \(H\), on an initial trial vector, \(\psi_{t}\),  
\[K^{j}=\{\psi_{t},H\psi_{t},H^{2}\psi_{t},H^{3}\psi_{t},\cdots\,H^{j-1}\psi_{t}\}. \tag{14}\]  
As all the vectors generated converge towards the eigenvector with the lowest eigenvalue, the basis set that is generated is ill-conditioned for large \(j\). The solution is to orthogonalize the basis at each step of the iteration via the Gram-Schmidt process. In essence, the difference between the Lanczos and Arnoldi methods is in the number of vectors in the Gram-Schmidt process. The Arnoldi method uses all the vectors and the Lanczos method only uses the two most recently generated vectors. The original matrix can then be projected into the Krylov subspace of much smaller size, where it may be fully diagonalized [145]. The dominant component of the computation is the matrix-vector multiplication described above. This scales only linearly with the matrix size. For the ground state calculation, matrix sizes of over one billion are routinely done [146]; however, calculating the inner spectrum is somewhat more difficult. The matrix must be shifted and then inverted to transform the target eigenvalue, \(\Lambda\), to the extremal eigenvalue. \[(H-\Lambda\mathbf{I})^{-1}\psi=\frac{1}{E-\Lambda}\psi, \tag{15}\]  
The inverse of the Hamiltonian with a shifted spectrum is generally not known. Then, instead of expanding the basis in the Krylov subspace, the Jacobi-Davidson method (JDM) is often employed [147]. It expands the basis \((\mathbf{u_{0}},\mathbf{u_{1}},\mathbf{u_{2}},\cdots)\) using the Jacobi orthogonal component correction which may be written as  
\[H(\mathbf{u_{j}}+\delta)=(\theta_{j}+\epsilon)(\mathbf{u_{j}}+\delta)\quad \forall\quad\mathbf{u_{j}}\perp\delta, \tag{16}\]  
where \((\mathbf{u_{j}},\theta_{j})\) and \((\mathbf{u_{j}}+\delta,\theta_{j}+\epsilon)\) are the approximate and the exact eigenvector and eigenvalue pairs, respectively. Upon solving the equation for the vector \(\delta\), a new basis vector \(\mathbf{u_{j+1}}=\mathbf{u_{j}}+\delta\) is included in the subspace. Matrix inversion is again involved in solving the equation. Various pre-conditioners are proposed for a quick approximation of the matrix inverse [147]. JADAMILU is a popular package which implements the JDM with an incomplete LU factorization [148; 149] as a pre-conditioner [150]. The scaling of this method seems to be strongly dependent on the Hamiltonian. It tends to be more efficient for matrices which are diagonally dominant, but much less so when off-diagonal matrix elements are large. This is probably due to the difficulty of obtaining a good approximation of the inverse based on the incomplete LU factorization used as a pre-conditioner. Exact diagonalization methods provide an accurate variational approximation for the eigenvalues and eigenvectors of the Hamiltonian, thus allowing the calculation of quantities such as multifractal spectrum and entanglement spectrum which are difficult to obtain from other approaches [151; 152]. On the other hand, Krylov subspace methods are not a good option for calculating the DOS as only one, or a few, eigenstates are targeted at each calculation.