efficient way to compute the inversion of the following matrix

Question

Hi, there

I have looked it up in the current textbook. The conventional numerical method to compute the inversion of an $n \times n$ matrix requires $O(n^3)$. However, for the following special matrix ${(E-a \cdot X)}^{-1}$,

where $E$ is an $n \times n$ identity matrix, $ a \in (0,1)$ is a scalar, and $X$ is a sparse stochastic matrix (the sum of each row is 1, and all its entries are between 0 and 1),

do you have some ideas to compute ${(E-a \cdot X)}^{-1}$ as fast as you can? (to reduce its $O(n^3)$ complexity is preferable, or approximate solution is also acceptable) ? Thanks in advance!

Answer 1

As others have already indicated, there's a good chance that you don't actually need $(E-aX)^{-1}$, but rather $(E-aX)^{-1}v$ for some vector $v$. In that case, you might be better off using an iterative method to solve $(E-aX)y=v$.

Another option to consider is the power series:

$ (E-aX)^{-1}=E+aX+a^2X^2+a^3X^3+... $

which is convergent for $|a|<1$, and can converge quite quickly if $a$ is small. If $a$ is small enough, then perhaps using a few terms of this power series will be sufficient. Computing powers of $X$ also takes $O(n^3)$ time in the worst case, but if $X$ is extremely sparse it might be much faster in practice than computing the inverse of $(E-aX)$.

In computing a partial sum of this series, you can use Horner's rule and get successive partial sums:

$ M^{1}=E+aX $

$ M^{2}=E+aX(M^{1})=E+aX+a^2X^2 $

$ M^{3}=E+aX(M^{2})=E+aX+a^2X^2+a^3X^3 $

$ \ldots $

Note that $(E-aX)^{-1}$ will typically be fully dense, but the sum of the first few terms of this series could be relatively sparse.

Answer 2

${(E-a \cdot X)}^{-1}=E+a \cdot X+\cdots+(a \cdot X)^k + O(a^{k+1}).$ General matrix multiplication can be $O(n^{\log_2 7})$ (I'm not sure how bad the constant is though) and since you say that $X$ is sparse there may be faster methods.

Answer 3

$(E-aX)^{-1}$ can be solved iteratively by representing $X$ in rank-1 decomposition. Since we are not using the special structure of X, absorb $a$ into $X$ for simplification. Assume $X$ has rank $m<n$. Let $u, v$ vectors of size n. Assuming you have access to rank-1 decomposition, $ X = \sum_{i=1}^{m} u_i v'_i $

Define $X_k=X_{k-1}+u_kv'_k,$ with $X_1 = u_1v'_1, $ and $X=X_m.$ Now let $ M_k = (E + X_k)^{-1},$ with $M_0=E$, also define $w'_k=v'_kM_{k-1}.$ We can write

\begin{eqnarray*} M_k = (E + X_k)^{-1} & = & (M_{k-1}^{-1} + u_kv'_k)^{-1} \\ & = & (M_{k-1}-\frac{M_{k-1}u_kw'_k}{1+w'_ku_k}) \end{eqnarray*} Note that $1+w'_ku_k$ is a scalar quantity and no matrix-matrix multiplication is required.

Answer 4

I have a similar issue with inverting sparse matrices (inverting Hessian for Newton-like optimization). As was pointed out to me in this blog, inverse of sparse symmetric matrices not only does not have to be sparse, it cannot be even approximated with sparse matrices: an inverse of a tridiagonal matrix may have dense entries that are far from zero. Therefore, by inverting your matrix you may be loosing the performance advantage of a sparse matrix.

Many operations on sparse matrices are O(n^2) even though they are O(n^3) (or O(n^2.8)-ish) for general matrices. Besides that, each of those O(n^2) operations would be much faster on a typical computer because less amount of data passed back and forth would lead to better cache performance.

Therefore, as pointed out in other answers, you'll be much better of if you don't inverse the matrix, but use other means to solve the end problem you are working on. Something like a decomposition that would preserve sparsity would serve your (supposed) end goals much better.