codenamev said:
Homework Statement
Hi, I'm trying to calculate the inverse of a really big matrix (4096x4096) using matlab. The inversion process takes too long on my pc, so i want to find a faster way.
The matrix has a strange property that i think could be useful to solve the problem:
splitting it in 64 blocks (8 rows and 8 columns of 256x256 submatrices), each submatrice is diagonal. Inverting a diagonal matrix is really simple, so i am asking you if this property can be used to calculate the inversion in a faster way. Homework Equations
The Attempt at a Solution
i found something about block diagonal matrix inversion, but my matrix is a bit different from a block diagonal matrix, because i have no zero submatrices... each block is diagonal...
i'm trying to use the matrix inversion in block form but i feel that there must be a smarter solution...
thanks
You are correct: the special structure allow you to reduce the computations from many billions of multiplications and additions, to a few hundred thousand. It gets complicated (but straightforward in principle), so I will just illustrate it for two cases: (I) the 2n x 2n case (4 nxn diagonal blocks); and (II) the 3n x 3n case (9 nxn diagonal blocks).
(I) M = \pmatrix{A&B\\D&E}, \text{ where } A = \text{diag}(a_1, a_2,\ldots,a_n), \text{ etc.} We want to solve the equations MZ = R, where
R = \pmatrix{u\\v}, u, \, v = \text{given n-dimensional vectors} and
Z = \pmatrix{x\\y}, \; x,\, y = \text{unknown n-dimensional vectors.}
These become
Ax + By = u \longrightarrow x = A^{-1}u - A^{-1}B y\\<br />
Dx+Ey = v \longrightarrow D(A^{-1}u - A^{-1}B y) + Ey = v,\\<br />
\text{so } E'y = v', \text{where } E' = E - D A^{-1}B, \; v' = v - D A^{-1} u.<br />
So, we can get y, then can get x from a previous formula.
Note: this is simplified by using the following simple facts: (1) The inverse of a diagonal matrix is diagonal (if it exists), that is: (\text{diag}(a_1,\ldots,a_n))^{-1} = \text{diag}(1/a_1, \ldots, 1/a_n). (2) The product of two diagonal matrices is diagonal; that is:
\text{diag}(a_1,\ldots,a_n) \cdot \text{diag}(b_1,\ldots,b_n)<br />
= \text{diag}(a_1 b_1, \ldots, a_n b_n)..
(3) The product of a scalar and a diagonal matrix is diagonal; that is
c\: \text{diag}(a_1,\ldots,a_n) = \text{diag}(ca_1,\ldots, ca_n), for any constant c.
(4) The sum of diagonal matrices is diagonal; that is
\text{diag}(a_1,\ldots,a_n) + \text{diag}(b_1,\ldots,n_n) <br />
= \text{diag}(a_1 + b_1, \ldots, a_n + b_n).
Therefore, the matrix E' above is diagonal, with (i,i)-element equal to
E'(i,i) \equiv e'_i = e_i - d_i b_i/a_i, so no actual matrix multiplications are needed. Further, solving E'y = v' for y and Ax = u-By for x is easy.
(II) Solve
\pmatrix{A&B&C\\D&E&F\\G&H&K} \pmatrix{x\\y\\z} = \pmatrix{u\\v\\w},
where each of the listed sub-matrices are nxn diagonal.
We can write the first two equations as
Ax + By = U, Dx + Ey = V, where U = u-Cz,\; V = v - Fz, so from above we have
E' y = V' = V - DA^{-1}U = v - DA^{-1}u + (DA^{-1}C -F)z,
which we can write as E' y = v' + F' z.
We can substitute this into the equation for x to get Ax = u' + C'z, where we get u' and C' from messy but straightforward algebra. So now we have x and y solved in terms of z, and can substitute these expressions into the third equation
Gx + Hy + Kz = w,
to get an equation of the form K' z = w'. All the matrices E', F', C', K', etc., are diagonal, so all the computations are very easy.
You can extend the method to 64 blocks in the same way; it just needs lots of labels to keep track of the different matrices throughout. All of the matrices encountered are diagonal, so it is enough to store their diagonals, and no acutal matrix multiplications or inversions are needed.
RGV
