Inverting matrix of unknown size

I have a matrix that has a variable number of columns and rows (but it is always square). The matrix is:

$$\Gamma = \begin{bmatrix} \Theta & 0 & 0 & 0 & T & 0 & R_j \\ u \Theta & \rho & 0 & 0 & u T & 0 & R_j \\ v \Theta & 0 & \rho & 0 & v T & 0 & R_j \\ w \Theta & 0 & 0 & \rho & w T & 0 & R_j \\ H \Theta -1 & u \rho & v \rho & w \rho & \Omega \rho +H T & \frac{5 \rho }{3} & R_j \\ k \Theta & 0 & 0 & 0 & k T & \rho & k R_j \\ Y_i \Theta & 0 & 0 & 0 & Y_i T & 0 & Y_i R_j+\rho \delta_{ij} \\ \end{bmatrix}$$

where the indices $i,j$ are the row and column indices respectively. For example, if $i,j={1,2}$ the matrix would look like:

$$\Gamma = \begin{bmatrix} \Theta & 0 & 0 & 0 & T & 0 & R_1 & R_2 \\ u \Theta & \rho & 0 & 0 & u T & 0 & R_1 & R_2 \\ v \Theta & 0 & \rho & 0 & v T & 0 & R_1 & R_2 \\ w \Theta & 0 & 0 & \rho & w T & 0 & R_1 & R_2 \\ H \Theta -1 & u \rho & v \rho & w \rho & \Omega \rho +H T & \frac{5 \rho }{3} & R_1 & R_2 \\ k \Theta & 0 & 0 & 0 & k T & \rho & k R_1 & k R_2 \\ Y_1 \Theta & 0 & 0 & 0 & Y_1 T & 0 & Y_1 R_1+\rho & Y_1 R_2 \\ Y_2 \Theta & 0 & 0 & 0 & Y_2 T & 0 & Y_2 R_1 & Y_2 R_2 + \rho \\ \end{bmatrix}$$

The values for $i$ and $j$ can be arbitrary.

When I put this in as is (just treating the subscripts as notation), I get terms divided by the Kronecker delta only which is impossible. Is there a way to invert this and keep it general in terms of $i,j$? The code to enter the matrix is:

G =
{
{Θ, 0, 0, 0, T, 0, Rj},
{Θ*u, ρ, 0, 0, T*u, 0, Rj},
{Θ*v, 0, ρ, 0, T*v, 0, Rj},
{Θ*w, 0, 0, ρ, T*w, 0, Rj},
{Θ*H - 1, ρ*u, ρ*v, ρ*w, T*H + ρ*Ω, 5/3*ρ, Rj},
{Θ*k, 0, 0, 0, T*k, ρ, Rj*k},
{Θ*Yi, 0, 0, 0, T*Yi, 0, Rj*Yi + ρ*δ}
};
Ginv = FullSimplify[Inverse[G]];

• Can you explain what is wrong with this? Ginv.G gives the right thing. Also, what does "terms divided by the Kronecker delta only" mean, and why is it impossible? – acl Jul 10 '14 at 22:51
• Mathematica does not know what these symbols mean mathematically. there are just symbols. If you want actual dirac delta, you should use DiracDelta[] btw, you might want to add space here "[Rho][Delta]" to make it "[Rho] [Delta]" just in case. – Nasser Jul 10 '14 at 22:53
• @acl If I look at the Ginv, there are terms like $Y_i/\delta_{ij}$ which will be $Y_i/0$ for all but 1 column, which is ill-defined. – tpg2114 Jul 10 '14 at 22:54
• @Nasser But how can I tell Mathematica that $i,j$ are for columns and rows? Without that information, the DiracDelta[] wouldn't know what to do either right? – tpg2114 Jul 10 '14 at 22:55
• I do not understand what you have there, so can't help. I do not know what Y_i supposed to mean. And do not know what R_j suppose to mean. Where did these come from? What is "i" there? what is "j" ? your question is fully described well to answer it (for me at least) – Nasser Jul 10 '14 at 23:03

This seems to work although I have not tested it thoroughly. Also, it's a rather brute force approach, so I hope someone else will post a more elegant solution.

g[0] =
{{Θ, 0, 0, 0, T, 0},
{u Θ, ρ, 0, 0, T u, 0},
{v Θ, 0, ρ, 0, T v, 0},
{w Θ, 0, 0, ρ, T w, 0},
{-1 + H Θ, u ρ, v ρ, w ρ, H T + ρ Ω, (5 ρ)/3},
{k Θ, 0,0, 0, k T, ρ}}

upperRight[n_] := Transpose@Table[Join[ConstantArray[R[j], 5], {k R[j]}], {j, n}]

lower[n_] :=
Table[Join[{Θ*Y[i], 0, 0, 0, T*Y[i], 0}, Table[Y[i] R[j], {j, n}]], {i, n}] /.
p : R[x_] Y[x_] :> p + ρ

g[n_] := Join[MapThread[Join, {g[0], upperRight[n]}, 1], lower[n]]

g2 = g[2]

{{Θ, 0, 0, 0, T, 0, R[1], R[2]},
{u Θ, ρ, 0, 0, T u, 0, R[1], R[2]},
{v Θ, 0, ρ, 0, T v, 0, R[1], R[2]},
{w Θ, 0, 0, ρ, T w, 0, R[1], R[2]},
{-1 + H Θ, u ρ, v ρ, w ρ, H T + ρ Ω, (5 ρ)/3, R[1], R[2]},
{k Θ, 0, 0, 0, k T, ρ, k R[1], k R[2]},
{Θ Y[1], 0, 0, 0, T Y[1], 0, ρ + R[1] Y[1], R[2] Y[1]},
{Θ Y[2], 0, 0, 0, T Y[2], 0, R[1] Y[2], ρ + R[2] Y[2]}}


The above is invertible although the result is messy.

Short[Inverse[g2], 4]


As it turns out, by employing block matrix inversion, one can derive a general formula for the inverse of $\mathbf \Gamma$. I'll leave the rigorous algebra/justification for other people, but here is how to apply block inversion:

emat = {{Θ, 0, 0, 0, T, 0},
{Θ u, ρ, 0, 0, T u, 0},
{Θ v, 0, ρ, 0, T v, 0},
{Θ w, 0, 0, ρ, T w, 0},
{Θ H - 1, ρ u, ρ v, ρ w, T H + ρ Ω, 5 ρ/3},
{Θ k, 0, 0, 0, T k, ρ}};

n = 7; (* size of lower right block, change as needed *)

(* upper right *)
ff = DiagonalMatrix[{1, 1, 1, 1, 1, k}].
Transpose[Table[ConstantArray[Subscript[R, j], 6], {j, n}]];

(* lower left *)
gg = Table[{Subscript[Y, i] Θ, 0, 0, 0, Subscript[Y, i] T, 0}, {i, n}];

(* lower right *)
hh = Outer[Times, Table[Subscript[Y, i], {i, n}], Table[Subscript[R, j], {j, n}]] +
ρ IdentityMatrix[n];


Generate the full matrix $\mathbf \Gamma$:

gmat = ArrayFlatten[{{emat, ff}, {gg, hh}}];


Now, assemble the inverse in blocks:

einv = Simplify[Inverse[emat]]; (* invert upper left block *)

vec = Simplify[Total[einv.DiagonalMatrix[{1, 1, 1, 1, 1, k}], {2}]];

dp = Table[Subscript[R, j], {j, n}].Table[Subscript[Y, j], {j, n}];


Build $\mathbf \Gamma^{-1}$:

ginv = ArrayFlatten[{{einv + Transpose[PadRight[{vec dp/ρ}, {6, 6}]],
-Outer[Times, vec, Table[Subscript[R, j], {j, n}]]/ρ},
IdentityMatrix[n]/ρ}}];


Check:

gmat.ginv == IdentityMatrix[6 + n] // Simplify
True


In symbolic format, the inverse can be expressed as

$$\mathbf\Gamma^{-1}=\small \begin{pmatrix} \frac{3 H T-5 k T-3 \left(u^2+v^2+w^2\right)T+3\rho\Omega}{3(T+\Theta\rho\Omega)}+\frac{\left(T\left(H-u^2+u-v^2+v-w^2+w-1\right)+\rho\Omega\right)\sum_k R_k Y_k}{\rho(\Theta\rho\Omega+T)}&\frac{T u}{T+\Theta\rho\Omega} &\frac{T v}{T+\Theta\rho\Omega}&\frac{T w}{T+\Theta\rho \Omega}&-\frac{T}{T+\Theta\rho\Omega}&\frac{5T}{3(T+\Theta\rho\Omega)}&-\frac{T\left(H-u^2+u-v^2+v-w^2+w-1\right)+\rho\Omega}{\rho(\Theta\rho\Omega+T)}R_j\\ \frac{(1-u)\sum_k R_k Y_k}{\rho^2}-\frac{u}{\rho}&\frac{1}{\rho} & 0 & 0 & 0 & 0 &\frac{u-1}{\rho^2}R_j\\ \frac{(1-v)\sum_k R_k Y_k}{\rho^2}-\frac{v}{\rho} & 0 & \frac{1}{\rho} & 0 & 0 & 0 &\frac{v-1}{\rho^2}R_j\\ \frac{(1-w)\sum_k R_k Y_k}{\rho^2}-\frac{w}{\rho} & 0 & 0 & \frac{1}{\rho} & 0 & 0 &\frac{w-1}{\rho^2}R_j\\ \frac{-3 H\Theta+5 k\Theta+3\left(u^2+v^2+w^2\right)\Theta+3}{3(T+\Theta\rho\Omega)}+\frac{\left(1-\Theta\left(H-u^2+u-v^2+v-w^2+w-1\right)\right)\sum_k R_k Y_k}{\rho(\Theta\rho\Omega +T)}& -\frac{u\Theta}{T+\Theta\rho\Omega}& -\frac{v\Theta}{T+\Theta\rho\Omega} & -\frac{w\Theta}{T+\Theta\rho\Omega} & \frac{\Theta}{T+\Theta\rho\Omega}& -\frac{5\Theta}{3 (T+\Theta\rho\Omega)}&-\frac{1-\Theta\left(H-u^2+u-v^2+v-w^2+w-1\right)}{\rho(\Theta\rho\Omega +T)}R_j\\ -\frac{k}{\rho} & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 \\ -\frac{Y_i}{\rho} & 0 & 0 & 0 & 0 & 0 & \frac{\delta_{ij}}{\rho} \end{pmatrix}$$

The block inversion worked as well as it did because of fortuitous simplifications. For instance, here are some useful identities I determined:

einv.ff == Outer[Times, vec, Table[Subscript[R, j], {j, n}]] // Simplify
True

gg.einv == PadRight[Table[{Subscript[Y, i]}, {i, n}], {n, 6}] // Simplify
True

Inverse[hh - gg.einv.ff] == IdentityMatrix[n]/ρ // Simplify
True


It was by using these identities that I was able to construct the explicit expression for $\mathbf \Gamma^{-1}$.