# Numerically stable Lanczos process? I need to compute Elements of inverse in sparsity pattern of A

I have a large sparse symmetric positive definite matrix NxN matrix $A$. Let $s$ be the average number of non-zeros per row (i.e. $sN$ total non-zero elements). I would like to compute the elements of the inverse of A only at locations where A is non-zero as quickly as possible (this arises when computing expected pairwise marginal statistics of Gaussian Markov Random Fields).

My idea is as follows: compute an Orthogonal Tridiagonalization of A with the Lanczos process (that is, $A=Q^TTQ$ where $Q$ is orthogonal and $T$ is tridiagonal NxN matrices). The full tridiagonalization in theory will arise after N iterations, each iteration requiring a single matrix vector product, which costs $O(sN)$ time, yielding an $O(sN^2)$ algorithm.

Let $q_i$ be the $i^{th}$ column of $Q$. If $A_{ij}$ is non-zero, we can compute $A^{-1}_{ij} = e_i^TA^{-1}e_j = q_i^T T^{-1} q_j$. This requires a tridiagonal system solve and a matrix vector product, for a total of $O(N)$ operations. Computing the desired $sN$ elements thus requires $O(sN^2)$ operations, yielding a full algorithm complexity of $O(sN^2)$.

From what I understand, the "textbook" Lanczos process is not numerically stable, Q-vectors lose orthogonality and must be re-orthogonalized with respect to the already computed basis, thus destroying the $O(sN^2)$ complexity.

Is there any numerically stable way to compute the tri-diagonalization in $O(sN^2)$ (or a way to compute the desired elements of the inverse via another method in $O(sN^2)$ time)? If not, it seems like the only option is full inversion which would cost $O(N^3)$, a MUCH larger complexity. Could we get away with the numerically unstable algorithm by using extended precision arithmetic or is this a bad idea?