det simplification

Question

I have a big matrix with 8 by 8 and each of element is a polynomial function with regard to w. I need to calculate det[Matrix-identitymatrix], but since the large calculation, Mathematica cannot handle it, but actually, I only need the first coefficient of w in det[Matrix-identitymatrix], so I get rid of high power: Det[QM - IdentityMatrix] // Function[y, Normal[y + O[w]^2]]. but still cannot get the result after half hour, do I have any other option? Thanks, guys!

Answer 1

Could just interpolate. Alternatively, to get just constant and linear terms could do as follows.

(1) Evaluate determinant at w->0 for constant term.

(2) Successively take derivative of each row to obtain 8 new matrices, evaluate determinant of each with w->0, and take sum. This gives the coefficient of the linear term.

Random example:

randomPoly[deg_, max_, x_] := 
 RandomInteger[{-max, max}, deg + 1].x^Range[0, deg]
randomMatrix[m_, n_, deg_, max_, x_] := 
 Table[randomPoly[deg, max, x], {m}, {n}]
SeedRandom[1111];
n = 4;
mat = randomMatrix[n, n, 3, 10, x]

(* Out[62]= {{-8 + 5 x - 3 x^2 - 8 x^3, -4 + 2 x + 4 x^2 + 2 x^3, -2 + 
   7 x - x^2 - 10 x^3, 
  9 - 6 x - 3 x^2 - 2 x^3}, {5 - 2 x + 4 x^2 + 6 x^3, 
  4 + 5 x^2 - 8 x^3, -5 + 4 x + 5 x^2 - 3 x^3, -6 - 6 x - 7 x^2 - 
   9 x^3}, {8 + 9 x + 2 x^2 + 10 x^3, 1 + 5 x + 8 x^2 - 8 x^3, 
  4 + 5 x - 2 x^2 + x^3, -9 - 10 x - 4 x^2 - 3 x^3}, {3 + 10 x, -2 + 
   4 x + x^2 + 7 x^3, 
  10 - 5 x - 3 x^2 - 10 x^3, -6 - 8 x - 6 x^2 + 7 x^3}} *)

First find the determinant the hard way.

Det[mat - IdentityMatrix[n]]

(* Out[64]= -588 + 608 x + 4765 x^2 + 2367 x^3 + 3916 x^4 - 812 x^5 - 
 21041 x^6 + 255 x^7 - 34513 x^8 + 13674 x^9 - 18497 x^10 - 
 2682 x^11 + 4336 x^12 *)

Constant term:

const = Det[mat - IdentityMatrix[n] /. x -> 0]

(* Out[68]= -588 *)

Linear term coefficient:

dmats = 
  Table[MapAt[D[#, x] &, mat - IdentityMatrix[n], j], {j, n}];
lincoeff = Total[Map[Det, dmats] /. x -> 0]

(* ut[67]= 608 *)

Answer 2

One can make use of matrix differentials, e.g. formula (42) of The Matrix Cookbook:

$$\partial \left(\operatorname{det} \left(X\right) \right) = \operatorname{det} \left(X\right) \operatorname{Tr}\left[X^{-1}.\partial X\right]$$

Using @DanielLichtblau's example:

m[x_] = mat - IdentityMatrix[4];
Det[m[x]]

-588 + 608 x + 4765 x^2 + 2367 x^3 + 3916 x^4 - 812 x^5 - 21041 x^6 + 255 x^7 - 34513 x^8 + 13674 x^9 - 18497 x^10 - 2682 x^11 + 4336 x^12

we get:

Det[m[0]]
Det[m[0]] Tr[Inverse[m[0]] . m'[0]]

-588

608

in agreement with his answer.

Addendum

It is possible to define differentiation rules to allow one to compute higher order terms. The first step is to teach Mathematica special rules for derivatives of Det, Tr and Inverse. In order to do so, we need to prevent the normal derivative rules from applying, and this can be done by including these symbols in a system differentiate option for "ExcludedFunctions". The following code does this, and also defines the needed special rules for derivatives of these symbols:

MatrixD[expr_, x__] := With[
    {old = OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"]},

    Internal`WithLocalSettings[
        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->Join[old, {Det, Inverse, Tr}]];

        Unprotect[D];
        (* handle list derivatives *)
        D[h:((Det|Tr|Inverse)[m_]), {z_, n_Integer}] := Nest[D[#, Replace[z, _List :> {z}]]&, h, n];
        D[h:((Det|Tr|Inverse)[m_]), {z_List}] := D[h, #]& /@ z;
        D[h:((Det|Tr|Inverse)[m_]), z_, y___] := D[D[h, z], y];

        (* define derivatives for Det, Tr, and Inverse *)
        D[Det[m_], z:Except[_List]] := Det[m] Tr[Inverse[m] . D[m,z]];
        D[Tr[m_], z:Except[_List]] := Tr[D[m,z]];
        D[Inverse[m_], z:Except[_List]] := -Inverse[m] . D[m, z] . Inverse[m],

        D[expr, x],

        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->old];
        Clear[D];
        Protect[D]
    ]
]

Let's use the same example as before. We need to use mm[x] instead of m[x] so that Det doesn't evaluate to a polynomial prematurely.

Table[
    MatrixD[Det[mm[x]], {x, n}]/n! /. 
        x->0 /. 
        Inverse[mm[0]]->Inverse[m[0]] /. 
        mm->m,
    {n, 0, 6}
]

CoefficientList[Det[m[x]], x][[;;7]]

{-588, 608, 4765, 2367, 3916, -812, -21041}

{-588, 608, 4765, 2367, 3916, -812, -21041}