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!
2 Answers
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 *)
-
$\begingroup$ As a beginner, I'm quite confused with seed random, and could you explain it in a simple way, please? Thanks! $\endgroup$– Ivy GaoDec 17, 2017 at 13:08
-
$\begingroup$
SeedRandom
allows one to get "Random" results in a way that can be repeated. $\endgroup$ Dec 17, 2017 at 17:01
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}
-
1$\begingroup$ +1.
LinearSolve[m[0], m'[0]]
might be a more efficient (and stable on floating-point) alternative toInverse[m[0]].m'[0]
. $\endgroup$ Dec 16, 2017 at 2:13
w
the only parameter that is not numerical in the matrix? What are the degrees of the polynomials? $\endgroup$w
before the calculation of determinant. TrySeries
function. $\endgroup$Det[Map[Normal@Series[#, {w, 0, 1}] &, QM - IdentityMatrix[8], {2}]]
? $\endgroup$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$