4
$\begingroup$

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!

$\endgroup$
  • 2
    $\begingroup$ Can you provide (code that generates) a typical matrix that we can work with? Is w the only parameter that is not numerical in the matrix? What are the degrees of the polynomials? $\endgroup$ – Marius Ladegård Meyer Dec 14 '17 at 9:32
  • $\begingroup$ I guess you need in the first place to expand the 64 entries of your "big matrix" and keep up to the linear terms in w before the calculation of determinant. Try Series function. $\endgroup$ – Αλέξανδρος Ζεγγ Dec 14 '17 at 9:52
  • $\begingroup$ Just to expand on @AlexanderZeng's comment, have you tried Det[Map[Normal@Series[#, {w, 0, 1}] &, QM - IdentityMatrix[8], {2}]]? $\endgroup$ – aardvark2012 Dec 14 '17 at 10:17
  • $\begingroup$ I tried the method you mentioned, it becomes a little bit better. and I'm sorry I cannot give you the code as it's a very long code, and each matrix elements comes from a system and each of them may contain thousand item and 5 other variables. in fact, I have no idea what's it looks like, I can only guess. but thanks for the idea, it's very convenient and help. $\endgroup$ – Ivy Gao Dec 14 '17 at 13:45
  • 1
    $\begingroup$ If they are so complicated, and depend on other parameters than w, what do you expect to gain from your calculation? What will you learn from a super complicated symbolic coefficient? What are you going to use it for? Maybe this is an XY-problem... $\endgroup$ – Marius Ladegård Meyer Dec 14 '17 at 20:25
5
$\begingroup$

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 *)
$\endgroup$
  • $\begingroup$ As a beginner, I'm quite confused with seed random, and could you explain it in a simple way, please? Thanks! $\endgroup$ – Ivy Gao Dec 17 '17 at 13:08
  • $\begingroup$ SeedRandom allows one to get "Random" results in a way that can be repeated. $\endgroup$ – Daniel Lichtblau Dec 17 '17 at 17:01
6
$\begingroup$

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}

$\endgroup$
  • 1
    $\begingroup$ +1. LinearSolve[m[0], m'[0]] might be a more efficient (and stable on floating-point) alternative to Inverse[m[0]].m'[0]. $\endgroup$ – Michael E2 Dec 16 '17 at 2:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.