efficient way to compute the inversion of the following matrix

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!

-
What are you looking for? Good theoretical complexity or fast answer in practice? These are two different questions. –  Thierry Zell Mar 21 '11 at 6:21
I am just wondering whether there would be some effciient way to optimize the ${E-a \cdot X}^{-1}$ computation? –  Nancy Mar 21 '11 at 6:47
There's also the question of whether you really want the inverse matrix. (For example, if you are trying to solve a system of linear equations, then computing the inverse is not the best approach.) The matrix $E-aX$ is diagonally dominant, which is an important hypothesis for analyzing iterative algorithms, such as the Jacobi method for solving systems of equations. –  Henry Cohn Mar 21 '11 at 6:50
@Nancy, I think what Thierry is saying is that what's efficient in theory may not be efficient in practice, and vice versa. Also, what's efficient when $n=10$ may not be when $n=10,000$, and also the answer may depend on just how sparse is sparse. –  Gerry Myerson Mar 21 '11 at 6:52
Yeah, you can alteratively think of this problem as solve a linear equation in terms of $Y$ such that $(E-aX)Y=E$. –  Nancy Mar 21 '11 at 9:30

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.

-
Matrix multiplication can be done faster than $O(n^3)$ even for dense matrices; for example there are implementations of Strassen multiplication faster than naive one for $n>200$. On the other hand, you're multiplying (possibly) dense matrix by a sparse matrix. So there may be faster multiplication algorithm still. –  Vít Tuček Mar 22 '11 at 0:32
Strassen's algorithm (and other similar "faster than O(n^3)" algorithms simply haven't proven to be useful in practice, even for fairly large (e.g. 10000 by 10000 matrix times 10000 by 10000 matrix) matrix-matrix multiplies. This is mostly because memory bandwidth and access times rather than flops are the real driving factor in matrix multiplication times. In this case, you'd definitely want to take advantage of the sparsity of X in doing the matrix multiplications. –  Brian Borchers Mar 22 '11 at 2:24

${(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.

-

$(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.

-