# Inverse of a sparse 21x21 symbolic matrix

I am trying to compute the inverse of the following symmetric sparse(density 0.37) 21x21 matrix involving symbolic expressions in the entries(the variables a1,a2,a3,a4,a5,a6 shall correspond to positive half-integers). While there have been a number of similar questions floating around, I have, as yet, not managed to glean enough information to solve the problem. If one of the existing threads does indeed provide an answer which I am not recognising then feel free to call my attention to it, and sorry for the duplicate.

M = SparseArray[{
{1/2 (-2 a4 - 2 a5), 0, 0, a4/2, (I a4)/2, 0, a5/2, 0, (I a5)/2, a5/2, 0,
-((I a5)/2), a4/2, -((I a4)/2), 0, 0, 0, 0, 0, 0, 0},
{0, 1/2 (-2 a1 - 2 a4), 0, -((I a4)/2), a4/2, 0, 0, 0, 0, 0, 0, 0, (I a4)/2,
a4/2, 0, 0, 0, 0, 0, a1/2, (I a1)/2},
{0, 0,1/2 (-2 a1 - 2 a5), 0, 0, 0, ((I a5)/2), 0, a5/2, (I a5)/2, 0, a5/2,
0, 0, 0, 0, 0, 0, 0, -((I a1)/2), a1/2},
{a4/2, -((I a4)/2), 0,1/2 (-2 a2 - 2 a4), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
a4/2, (I a4)/2, 0, a2/2, 0, -((I a2)/2)},
{(I a4)/2, a4/2, 0, 0, 1/2 (-2 a4 - 2 a6), 0, 0, a6/2, (I a6)/2, 0, a6/2,
-((I a6)/2), 0, 0, 0, -((I a4)/2), a4/2, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1/2 (-2 a2 - 2 a6), 0, -((I a6)/2), a6/2, 0, (I a6)/2, a6/2,
0, 0, 0, 0, 0, 0, (I a2)/2, 0, a2/2},
{a5/2, 0, -((I a5)/2), 0, 0, 0, 1/2 (-2 a3 - 2 a5), 0, 0, 0, 0, 0, 0, 0, 0,
a5/2, 0, (I a5)/2, a3/2, (I a3)/2, 0},
{0, 0, 0, 0, a6/2, -((I a6)/2), 0, 1/2 (-2 a3 - 2 a6), 0, 0, 0, 0, 0, a6/2,
(I a6)/2, 0, 0, 0, -((I a3)/2), a3/2, 0},
{(I a5)/2, 0, a5/2, 0, (I a6)/2, a6/2, 0, 0, 1/2 (-2 a5 - 2 a6), 0, 0, 0,
0,-((I a6)/2), a6/2, -((I a5)/2), 0, a5/2, 0, 0, 0},
{a5/2, 0, (I a5)/2, 0, 0, 0, 0, 0, 0, 1/2 (-2 a3 - 2 a5), 0, 0, 0, 0, 0,
a5/2, 0, -((I a5)/2), a3/ 2, -((I a3)/2), 0},
{0, 0, 0, 0, a6/2, (I a6)/2, 0, 0, 0, 0, 1/2 (-2 a3 - 2 a6), 0, 0, a6/2,
((I a6)/2), 0, 0, 0, (I a3)/2, a3/2, 0},
{-((I a5)/2), 0, a5/2, 0, -((I a6)/2), a6/2, 0, 0, 0, 0, 0,
1/2 (-2 a5 - 2 a6), 0, (I a6)/2, a6/2, (I a5)/2, 0, a5/2, 0, 0, 0},
{a4/2, (I a4)/2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/2 (-2 a2 - 2 a4), 0, 0,
a4/2, -((I a4)/2), 0, a2/2, 0, (I a2)/2},
{-((I a4)/2), a4/2, 0, 0, 0, 0, 0, a6/2, -((I a6)/2), 0, a6/2,
(I a6)/2, 0, 1/2 (-2 a4 - 2 a6), 0, (I a4)/2, a4/2, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, (I a6)/2, a6/2, 0, -((I a6)/2), a6/2, 0, 0,
1/2 (-2 a2 - 2 a6), 0, 0, 0, -((I a2)/2), 0, a2/2}, {0, 0, 0, a4/2,
-((I a4)/2), 0, a5/2, 0, -((I a5)/2), a5/2, 0, (I a5)/2, a4/2, (
I a4)/2, 0, 1/2 (-2 a4 - 2 a5), 0, 0, 0, 0, 0},
{0, 0, 0, (I a4)/2, a4/2, 0, 0, 0, 0, 0, 0, 0, -((I a4)/2), a4/2, 0, 0,
1/2 (-2 a1 - 2 a4), 0, 0, a1/2, -((I a1)/2)},
{0, 0, 0, 0, 0, 0, (I a5)/2, 0, a5/2, -((I a5)/2), 0, a5/2, 0, 0, 0, 0, 0,
1/2 (-2 a1 - 2 a5), 0, (I a1)/2, a1/2},
{0, 0, 0, a2/2, 0, (I a2)/2, a3/2, -((I a3)/2), 0, a3/2, (I a3)/2, 0, a2/2,
0, -((I a2)/2), 0, 0, 0, 1/2 (-2 a2 - 2 a3), 0, 0},
{0, a1/2, -((I a1)/2), 0, 0, 0, (I a3)/2, a3/2, 0, -((I a3)/2), a3/2, 0, 0,
0, 0, 0, a1/2, (I*a1)/2, 0, 1/2 (-2 a1 - 2 a3), 0},
{0, (I a1)/2, a1/2, -((I a2)/2), 0, a2/2, 0, 0, 0, 0, 0, 0, (I a2)/2, 0,
a2/2, 0, -((I a1)/2), a1/2, 0, 0, 1/2 (-2 a1 - 2 a2)}
}]


Naively trying Inverse[M] results in exhausting memory pretty quickly. Since it has been possible to compute the determinant, my best shot has been trying to employ the Cayley-Hamilton method with the use of complete Bell polynomials (here on wikipedia). The inverse of M is given by the product of 1/Det[M] with a whose essential components are powers of M and Bell polynomials of traces of powers of M. Implementing the formula in the link:

t[l_] := -(l - 1)! Tr[M^l] // Simplify
tlist = Table[t[i], {i, 1, 21-1}];
B[x_List] := Sum[BellY[Length[x], k, x], {k, 1, Length[x]}]; (*build
complete Bell Polynomial from partial Bell Polynomials*)

S = IdentityMatrix[21] (-1)^(21 - 1)/(21 - 1)! B[Take[tlist, 21 - 1]] +
Sum[M^(s - 1) (-1)^(21 - 1)/(21 - s)! B[Take[tlist, 21 - s]], {s, 2, 21-1}]
(*the first summand is just manually doing the s=1 part of the whole sum
where M^0=IdentityMatrix[21]*)


The inverse of M is then Sum/Det[M]. S does indeed compute but it gives me such a large output that Mathematica won't display it without crashing let alone run Simplify[] on it(it saves to a 100kB .mx file). Of course, this is is a result but I am guessing the final result should be somewhat more compact.

I'm not at all sure whether the Cayley-Hamilton method is the best way to go, but I'm happy for any improvements and indeed other suggestions as to how to compute the inverse of this matrix.

• For what do you need the inverse of the matrix? What's the original problem you are solving? – Thies Heidecke Aug 21 '17 at 10:52
• If you're multiplying the inverse with a vector or a matrix, as opposed to the needing the bare inverse, then you don't need the inverse. Use LinearSolve[]. – J. M. is away Aug 21 '17 at 13:05
• (1) I agree with @J.M. that LinearSolve is almost certainly the way to go. (2) There is a typo of a missing space, so that a "variable" is Ia1. (3) I'm trying this with Method->"OneStepRowReduction". On the forward elimination it seems to get at least to row 18 (in a developer build with debugging information available). But it is bogging down. Not sure if it completes this part what will happen on the backward elimination. – Daniel Lichtblau Aug 21 '17 at 16:39
• Here is a result of sorts. In[306]:= AbsoluteTiming[ mInv = Inverse[mmat, Method -> "OneStepRowReduction"];] Out[306]= {1929.418725, Null} In[307]:= LeafCount[mInv] Out[307]= 642060435 In[309]:= LeafCount[mInv[[1, 1]]] Out[309]= 1300467 In[310]:= t11 = Together[mInv[[1, 1]]]; In[311]:= LeafCount[t11] Out[311]= 133802 So it may be possible to get a factor of 10 or so compression by using Together on all entries. It will be slow though. – Daniel Lichtblau Aug 21 '17 at 17:13
• Further information. Together on the first row takes around 23 minutes on my desktop machine and produces a total leaf count of slightly under a meg. The original leaf count of that row is around 41.5 meg. – Daniel Lichtblau Aug 21 '17 at 18:58