# Block matrix definition and inversion

I would like to define the following block matrix

$$A=\begin{bmatrix} H-G_1 & -G_1 & \ldots & -G_1 \\ -G_2 & H-G_2 & \ldots & -G_2 \\ \vdots & \ldots & \vdots & \vdots \\ -G_n &\ldots & \ldots & H-G_n \end{bmatrix}$$ where

$H=[H(k,l)]_{m\times m} , G_i=[g(k,l)]_{m\times m}, a_i=(a_{i,1},\ldots,a_{i,m})^T, F_i=(f_{i,1},\ldots,f_{i,m})^T$

I would also like to generate the inverse of matrix $A$:

$$\begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} =A^{-1}\begin{bmatrix} F_1 \\ F_2 \\ \vdots \\ F_n \end{bmatrix}~~~(*)$$

I'm not sure how to tell Mathematica to create the matrix and compute it's inverse. The following (broken) code is a first attempt:

H = Table[H[k, l], {k, 2 M}, {l, 2 M}];

For[i = 1, i <= n, i++, {
G[i] = Table[g[i, k, l], {k, 2 M}, {l, 2 M}];
F[i] = Table[f[i, x[l]], {l, 2 M}];
}];

A = Table[if[i = j, H - G[i], -G[i]], {i, n}, {j, n}];

For[i = 1, i <= n, i++, {a[i] = Table[a[i, k], {k, 2 M}];}];

F = Table[F[i], {i, n}];

a = Table[a[i], {i, n}];

a = Inverse[A]*F;

Print[MatrixForm[a]]

• What you are asking to do is mathematically impossible. The expression size of the inverse of a symbolic matrix increases super-exponentially if I recall correctly, and anything bigger than a $5\times 5$ matrix is completely impossible to work with. Since you're using block symbolic matrices, the size of the whole problem is even worse. Are you sure that you need to calculate the inverse, though? Is there some reason why a floating-point approximation wouldn't work? And in any case, solving linear systems by direct inversion is not good practice, as factorizations are more stable. Nov 23, 2014 at 14:46
• @DumpsterDoofus Matrix elements are real numbers. Nov 23, 2014 at 15:03

There are a number of errors in your code, so hopefully going through them one by one will help you in your task. First, I'll define numbers to insert into the matrices:

h[k_, l_] := k Sin[l];
g[i_, k_, l_] := 10/(1 + i^2 + k^2 + l^2) + RandomReal[];
f[i_, l_] := i l;
x[l_] := 1;
M = 10;
n = 4;


The first line of code,

H = Table[H[k, l], {k, 2 M}, {l, 2 M}];


has a problem: H is defined to be an array and a function simultaneously, which makes no sense, and understandably triggers a rather unwelcome infinite loop. Instead, use a different name h (see my definitions above) for the function that generates the matrix elements, like this:

H = Table[h[k, l], {k, 2 M}, {l, 2 M}];


or even this:

H = Array[h, {2 M, 2 M}];


The second line of code,

For[i = 1, i <= n, i++, {
G[i] = Table[g[i, k, l], {k, 2 M}, {l, 2 M}];
F[i] = Table[f[i, x[l]], {l, 2 M}];
}];


works, but is strongly discouraged. Instead of For loops, you can just delay-set the definitions, like this:

G[i_] := Table[g[i, k, l], {k, 2 M}, {l, 2 M}];
F[i_] := Table[f[i, x[l]], {l, 2 M}];


The third line of code,

A = Table[if[i = j, H - G[i], -G[i]], {i, n}, {j, n}];


contains a syntax error and a conceptual error. All built-in functions are capitalized in Mathematica, as they are considered English proper nouns, and so if should be replaced with If. Second, this yields a $4\times4\times20\times20$ rank-4 tensor, so you want to apply ArrayFlatten to make the block-structure manifest, like this:

A = ArrayFlatten[Table[If[i == j, H - G[i], -G[i]], {i, n}, {j, n}]];


The fourth, sixth, and seventh lines of code,

For[i = 1, i <= n, i++, {
a[i] = Table[a[i, k], {k, 2 M}];
}];
a = Table[a[i], {i, n}];
a = Inverse[A]*F;


is again workable, but discouraged, and also makes no sense; your problem is to solve for a, so why are you generating an array a and then overwriting it? To solve this problem, simply delete your fourth and sixth lines of code.

The fifth line of code,

F = Table[F[i], {i, n}];


generates a $4\times20$ matrix, and also has the same infinite loop error I mentioned in the first line of code, and presumably you want to flatten it, so you can instead do this:

FF = Flatten@Table[F[i], {i, n}];


Finally, the sixth line of code has a syntax error: matrix multiplication in Mathematica is done with ., rather than juxtaposition or the * operator, which is elementwise multiplication (which incidentally results in a size mismatch error). Additionally, the previous definitions of matrix elements were all symbolic, so to numericize them, use the N operator, like this:

a = Inverse[N@A].FF;


In summary, here's the patched-up version:

h[k_, l_] := k Sin[l];
g[i_, k_, l_] := 10/(1 + i^2 + k^2 + l^2) + RandomReal[];
f[i_, l_] := i l;
x[l_] := 1;
M = 10;
n = 4;

H = Table[h[k, l], {k, 2 M}, {l, 2 M}];
G[i_] := Table[g[i, k, l], {k, 2 M}, {l, 2 M}];
F[i_] := Table[f[i, x[l]], {l, 2 M}];
A = ArrayFlatten[Table[If[i == j, H - G[i], -G[i]], {i, n}, {j, n}]];
FF = Flatten@Table[F[i], {i, n}];
a = Inverse[N@A].FF


In addition, as I mentioned in the comment, matrix inversion is bad practice; instead, you can use a more stable factorization-based approach:

factA = LinearSolve[N@A];
a = factA[FF]

• @ DumpsterDoofus, Thanks a lot. I edited my code. Please see it. Nov 24, 2014 at 15:42
• @Selena: What is the purpose of the edit? Is there some aspect of my explanation that you need help understanding? You asked how to invert a matrix, and to do that, you can use Inverse. Nov 24, 2014 at 16:11
• @Selena As @DumpsterDoofus already pointed out something about your code that is discouraged, here information from the Wolfram page: There is a convention that built‐in Wolfram Language objects always have names starting with uppercase (capital) letters. To avoid confusion, you should always choose names for your own variables that start with lowercase letters. As you are already using F and G as well as M you are quite close to E and N, which are no good choices. Nov 28, 2014 at 13:10

On the basis of $A$ being defined as posted:

adef[g_List, h_] := Module[{n = Length[g]},
gs = ArrayFlatten[Table[g, {n}]];
hs = ArrayFlatten[
ReplacePart[ConstantArray[0, {n, n}], {i_, i_} :> h]];
hs - gs]


To illustrate:

adefex[g_List, h_] := Module[{n = Length[g]},
gs = ArrayFlatten[Table[g, {n}]];
hs = ArrayFlatten[
ReplacePart[ConstantArray[0, {n, n}], {i_, i_} :> h]];
{hs, gs, hs - gs}]
testg = RandomInteger[10, {4, 2, 2}];
testh = RandomInteger[10, {2, 2}];
Column[{MatrixForm /@ testg, MatrixForm@testh, The first row is $\{G_i\}$, the second row $H$, the third row block matrices with the last in the row $A$ and the last row the inverse of $A$.