# Determinant of large matrices: there's GOTTA be a faster way

WARNING this is a very long report and is likely going to cause boredom. Be warned!!

I've heard of the determinant of small matrices, such as:

$$\det \begin{pmatrix} a&b\\ c&d\\ \end{pmatrix} = ad-bc$$

case in point:

$$\det \begin{pmatrix} 57&48\\ 79&102\\ \end{pmatrix} = 57\times 102-48\times 79 =5814-3792 =2022$$

This is a pretty hefty example i found in one of my books on vectors and matrices. And there are much more complex examples. for instance, to find the determinant of a matrix of order 3, you do this:

\begin{align} &\det \begin{pmatrix} a&b&c\\ d&e&f\\ g&h&i\\ \end{pmatrix}\\ &=a\times \det \begin{bmatrix} e&f\\ h&i\\ \end{bmatrix}\\ &-b\times \det \begin{bmatrix} d&f\\ g&i\\ \end{bmatrix}\\ &+c\times \det \begin{bmatrix} d&e\\ g&h\\ \end{bmatrix} \end{align}

This sequence looks a bit simple, but in reality it blows up(becoimes increasingly large) after a while. for instance, with a $$5\times 5$$ matrix someone asked me to model, this is how my 'fun time' went:

\begin{align} &\det \begin{Bmatrix} a&b&c&d&e\\ f&g&h&i&j\\ k&l&m&n&o\\ p&q&r&s&t\\ u&v&w&x&y\\ \end{Bmatrix}\\ &=a\times \det \begin{Bmatrix} g&h&i&j\\ l&m&n&o\\ q&r&s&t\\ v&w&x&y\\ \end{Bmatrix} -b\times \det \begin{Bmatrix} f&h&i&j\\ k&m&n&o\\ p&r&s&t\\ u&w&x&y\\ \end{Bmatrix} +c\times \det \begin{Bmatrix} f&g&i&j\\ k&l&n&o\\ p&q&s&t\\ u&v&x&y\\ \end{Bmatrix}\\ &-d\times \det \begin{Bmatrix} f&g&h&j\\ k&l&m&o\\ p&q&r&t\\ u&v&w&y\\ \end{Bmatrix} +e\times \det \begin{Bmatrix} f&g&h&i\\ k&l&m&n\\ p&q&r&s\\ u&v&w&x\\ \end{Bmatrix} \end{align}

This is a complex wad of calculations for me to completely do. so I'll break it down into the 5 conponents: A, B, C, D, and E, respectively.

$$A=a\times \det \begin{Bmatrix} g&h&i&j\\ l&m&n&o\\ q&r&s&t\\ v&w&x&y\\ \end{Bmatrix} \\ =a\left( g\times \det \begin{Bmatrix} m&n&o\\ r&s&t\\ w&x&y\\ \end{Bmatrix} -h\times \det \begin{Bmatrix} l&n&o\\ q&s&t\\ v&x&y\\ \end{Bmatrix} +i\times \det \begin{Bmatrix} l&m&o\\ q&r&t\\ v&w&y\\ \end{Bmatrix} -j\times \det \begin{Bmatrix} l&m&n\\ q&r&s\\ v&w&x\\ \end{Bmatrix} \right)\\ =a\left( g\left( m\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} \right)\\ -h\left( l\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} \right)\\ +i\left( l\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} \right) -j\left( l\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} \right) \right)\\ =a\left( g\left(m(sy-xt)-n(ry-wt)+o(rx-ws)\right)\\ -h\left(l(sy-xt)-n(qy-vt)+o(qx-vs)\right)\\ +i\left(l(ry-wt)-m(qy-vt)+o(qw-vr)\right)\\ -j\left(l(rx-ws)-m(qx-vs)+n(qw-vr)\right) \right)$$

(If you want to see this behemoth in code form, go to this page, but i'm not $$100$$% sure that it will work.)

$$B= -b\times \det \begin{Bmatrix} f&h&i&j\\ k&m&n&o\\ p&r&s&t\\ u&w&x&y\\ \end{Bmatrix}\\ -b\left( f\times \det \begin{Bmatrix} m&n&o\\ r&s&t\\ w&x&y\\ \end{Bmatrix} -h\times \det \begin{Bmatrix} k&n&o\\ p&s&t\\ u&x&y\\ \end{Bmatrix} +i\times \det \begin{Bmatrix} k&m&o\\ p&r&t\\ u&w&y\\ \end{Bmatrix} -j\times \det \begin{Bmatrix} k&m&n\\ p&r&s\\ u&w&x\\ \end{Bmatrix} \right)\\ =-b\left( f\left( m\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} \right)\\ -h\left( k\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} \right)\\ +i\left( k\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} \right) -j\left( k\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} \right) \right)\\ =-b\left( f\left(m(sy-xt)-n(ry-wt)+o(rx-ws)\right)\\ -h\left(k(sy-xt)-n(py-ut)+o(px-us)\right)\\ +i\left(k(ry-wt)-m(py-ut)+o(pw-ur)\right)\\ -j\left(k(rx-ws)-m(px-us)+n(pw-ur)\right) \right)$$

and that is part b! this is a grueling amount of code for me to place. $$\frac{3}{5}$$ way to go...

$$C=c\times \det \begin{Bmatrix} f&g&i&j\\ k&l&n&o\\ p&q&s&t\\ u&v&x&y\\ \end{Bmatrix}\\ =c\left( f\times \det \begin{Bmatrix} l&n&o\\ q&s&t\\ v&x&y\\ \end{Bmatrix} -g\times \det \begin{Bmatrix} k&n&o\\ p&s&t\\ u&x&y\\ \end{Bmatrix} +i\times \det \begin{Bmatrix} k&l&o\\ p&q&t\\ u&v&y\\ \end{Bmatrix} -j\times \det \begin{Bmatrix} k&l&n\\ p&q&s\\ u&v&x\\ \end{Bmatrix} \right)\\ =c\left( f\left( l\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} \right)\\ -g\left( k\times \det \begin{Bmatrix} s&t\\ x&y\\ \end{Bmatrix} -n\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} \right)\\ +i\left( k\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right)\\ -j\left( k\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right) \right)\\ =c\left( f\left(l(sy-xt)-n(qy-vt)+o(qx-vs)\right)\\ -g\left(k(sy-xt)-n(py-ut)+o(px-us)\right)\\ +i\left(k(qy-vt)-l(py-ut)+o(pv-uq)\right)\\ -j\left(k(qx-vs)-l(px-us)+n(pv-uq)\right) \right)$$

That's the C-section. now to get to the D-section...

$$D=-d\times \det \begin{Bmatrix} f&g&h&j\\ k&l&m&o\\ p&q&r&t\\ u&v&w&y\\ \end{Bmatrix}\\ =-d\left( f\times \det \begin{Bmatrix} l&m&o\\ q&r&t\\ v&w&y\\ \end{Bmatrix} -g\times \det \begin{Bmatrix} k&m&o\\ p&r&t\\ u&w&y\\ \end{Bmatrix} +h\times \det \begin{Bmatrix} k&l&o\\ p&q&t\\ u&v&y\\ \end{Bmatrix} -j\times \det \begin{Bmatrix} k&l&m\\ p&q&r\\ u&v&w\\ \end{Bmatrix} \right)\\ =-d\left( f\left( l\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} \right)\\ -g\left( k\times \det \begin{Bmatrix} r&t\\ w&y\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} \right)\\ +h\left( k\times \det \begin{Bmatrix} q&t\\ v&y\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&t\\ u&y\\ \end{Bmatrix} +o\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right)\\ -j\left( k\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} +m\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right) \right)\\ =-d\left( f\left(l(ry-wt)-m(qy-vt)+o(qw-vr)\right)\\ -g\left(k(ry-wt)-m(py-ut)+o(pw-ur)\right)\\ +h\left(k(qy-vt)-l(py-ut)+o(pv-uq)\right)\\ -j\left(k(qw-vr)-l(pw-ur)+m(pv-uq)\right) \right)$$

Are you bored yet? I am. Luckily, I got one more section to go...

$$E=e\times \det \begin{Bmatrix} f&g&h&i\\ k&l&m&n\\ p&q&r&s\\ u&v&w&x\\ \end{Bmatrix} =e\left( f\times \det \begin{Bmatrix} l&m&n\\ q&r&s\\ v&w&x\\ \end{Bmatrix} -g\times \det \begin{Bmatrix} k&m&n\\ p&r&s\\ u&w&x\\ \end{Bmatrix} +h\times \det \begin{Bmatrix} k&l&n\\ p&q&s\\ u&v&x\\ \end{Bmatrix} -i\times \det \begin{Bmatrix} k&l&m\\ p&q&r\\ u&v&w\\ \end{Bmatrix} \right)\\ =e\left( f\left( l\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} \right)\\ -g\left( k\times \det \begin{Bmatrix} r&s\\ w&x\\ \end{Bmatrix} -m\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} \right)\\ +h\left( k\times \det \begin{Bmatrix} q&s\\ v&x\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&s\\ u&x\\ \end{Bmatrix} +n\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right)\\ -i\left( k\times \det \begin{Bmatrix} q&r\\ v&w\\ \end{Bmatrix} -l\times \det \begin{Bmatrix} p&r\\ u&w\\ \end{Bmatrix} +m\times \det \begin{Bmatrix} p&q\\ u&v\\ \end{Bmatrix} \right) \right)\\ =e\left( f\left(l(rx-ws)-m(qx-vs)+n(qw-vr)\right)\\ -g\left(k(rx-ws)-m(px-us)+n(pw-ur)\right)\\ +h\left(k(qx-vs)-l(px-us)+n(pv-uq)\right)\\ -i\left(k(qw-vr)-l(pw-ur)+m(pv-uq)\right) \right)$$

ZZZZZZZZZZZZ... GAH! okay... to recap:

$$\det \begin{Bmatrix} a&b&c&d&e\\ f&g&h&i&j\\ k&l&m&n&o\\ p&q&r&s&t\\ u&v&w&x&y\\ \end{Bmatrix}\\ =a\left( g\left(m(sy-xt)-n(ry-wt)+o(rx-ws)\right)\\ -h\left(l(sy-xt)-n(qy-vt)+o(qx-vs)\right)\\ +i\left(l(ry-wt)-m(qy-vt)+o(qw-vr)\right)\\ -j\left(l(rx-ws)-m(qx-vs)+n(qw-vr)\right) \right)\\ -b\left( f\left(m(sy-xt)-n(ry-wt)+o(rx-ws)\right)\\ -h\left(k(sy-xt)-n(py-ut)+o(px-us)\right)\\ +i\left(k(ry-wt)-m(py-ut)+o(pw-ur)\right)\\ -j\left(k(rx-ws)-m(px-us)+n(pw-ur)\right) \right)\\ +c\left( f\left(l(sy-xt)-n(qy-vt)+o(qx-vs)\right)\\ -g\left(k(sy-xt)-n(py-ut)+o(px-us)\right)\\ +i\left(k(qy-vt)-l(py-ut)+o(pv-uq)\right)\\ -j\left(k(qx-vs)-l(px-us)+n(pv-uq)\right) \right)\\ -d\left( f\left(l(ry-wt)-m(qy-vt)+o(qw-vr)\right)\\ -g\left(k(ry-wt)-m(py-ut)+o(pw-ur)\right)\\ +h\left(k(qy-vt)-l(py-ut)+o(pv-uq)\right)\\ -j\left(k(qw-vr)-l(pw-ur)+m(pv-uq)\right) \right)\\ +e\left( f\left(l(rx-ws)-m(qx-vs)+n(qw-vr)\right)\\ -g\left(k(rx-ws)-m(px-us)+n(pw-ur)\right)\\ +h\left(k(qx-vs)-l(px-us)+n(pv-uq)\right)\\ -i\left(k(qw-vr)-l(pw-ur)+m(pv-uq)\right) \right)$$

Now that THAT'S over (STOP SCROLLING!!), I must mention that I pretty much blew my friend's mind showing him this. NOW he wants me to figure out a matrix of order 10. AURRRRRRRRUUUUUUUUUUUUGGGGGGGGGGHHHHHHHHHH!!!!!!! I DONT HAVE THE TIME!!!! Therefore, I am wondering if there is a faster way to calculate the determinant of a HUGE matrix. hope there is. Thanks in advance.

EDIT i was conversating with my friend, explaining how timewasting calculating a matrix of order 10 is, and i convinced him to drop the 'do by hand' idea, and instead do it on the computer.

• ... use a computer? – Bobson Dugnutt May 7 '17 at 0:58
• @Lovsovs any other ideas? – Alexander Day May 7 '17 at 0:59
• BTW is my answer to the determinant of order 5 matrix correct? i may have a few errors from here to there. – Alexander Day May 7 '17 at 1:01
• Computer is the way to go. There is no possible reason to ever use the formula you just wrote out... seeing it alone is awful, nevermind trying to actually use it – Brevan Ellefsen May 7 '17 at 1:05
• just curious, why the $#$%% WOULD YOU WRITE ALL OF THAT IN LATEX – Saketh Malyala Jun 4 '17 at 0:12

No, this is not the way that any (sane) person would compute a determinant. This is not even the way a computer would calculate a determinant! It requires a sum over $n!$ terms, which quickly becomes infeasible even for a computer, around $n = 15$ or so. An elementary way to compute a determinant quickly is by using Gaussian elimination.

We know a few facts about the determinant:

1. Adding a scalar multiple of one row to another does not change the determinant.
2. Interchanging two rows negates the determinant.
3. Scaling a row by a constant multiplies the determinant by that constant.

So, now take the matrix

$$A = \begin{bmatrix}-4 & 3 &3 \\ 8 & 7 & 3 \\ 4 & 3 & 3\end{bmatrix}$$

By fact (1) above, I can add twice the top row to the middle row, and also the top row to the bottom row, without affecting the determinant. So:

$$\det A = \det \begin{bmatrix}-4 & 3 &3 \\ 0 & 13 & 9 \\ 0 & 6 & 6\end{bmatrix}$$

Now, I can interchange the bottom two rows, and and scale the row with only $6$'s, at a cost of $-6$:

$$\det A = - \det \begin{bmatrix}-4 & 3 &3 \\ 0 & 6 & 6 \\ 0 & 13 & 9 \end{bmatrix} = - 6 \det \begin{bmatrix}-4 & 3 &3 \\ 0 & 1 & 1 \\ 0 & 13 & 9 \end{bmatrix}$$

Now I can subtract 13 lots of the middle row from the bottom row:

$$\det A = - 6 \det \begin{bmatrix}-4 & 3 &3 \\ 0 & 1 & 1 \\ 0 & 13 & 9 \end{bmatrix} = - 6 \det \begin{bmatrix}-4 & 3 &3 \\ 0 & 1 & 1 \\ 0 & 0 & -4 \end{bmatrix}$$

Now the matrix is upper-triangular, and so the determinant is just the product of the diagonal entries. So we have

$$\det A = -6 (-4 \times 1 \times -4) = -96$$

So there you have it: computing a determinant is as easy as finding row-echelon form.

• I'm agree with you. Computing the Gauss reduction and then multiplying the elements of the diagonal is possibly one of the best methods to compute the determinant for me. (Obviously you have to keep in mind that when you perform a permutation the symbol of the original determinant changes) – Rubén Ballester May 7 '17 at 12:15
• @Joppy interesting... this will be useful in the future!! – Alexander Day May 7 '17 at 12:42
• @AlexanderDay: Yes, it's quite easy to calculate determinants by hand this way. However, if I had to do a lot of them, I would rather use a computer. – Joppy May 7 '17 at 15:19

In general, determinants of large matrices are not computed by cofactor expansion but rather by factoring the matrix into factors whose determinants are easy to compute.

For example, you can factor an $n$ by $n$ matrix $A$ as

$A=P^{T}LU$

where $P$ is a permutation matrix, $L$ is lower triangular, and $U$ is upper triangular. This computation takes $O(n^{3})$ time for an $n$ by $n$ matrix $A$.

Since

$\det(A)=\det(P^{T})\det(L)\det(U)$

and the determinants of $P^{T}$, $L$, and $U$ are easy to compute (the determinant of a lower or upper triangular matrix is the product of the diagonal elements and you can easily do cofactor expansion on a permutation matrix), you can quickly find the determinant of $A$.

If you want to try a computational experiment, test MATLAB's det() function on randomly generated matrices of size $n$ by $n$ for $n=1000, 2000,\ldots, 10000$ and use tic/toc to see how long the computation takes.

If the entries of your matrix belong to a field, then you can compute the determinant easily using either LPU decomposition or PLU decomposition. These algorithms take $O \left(n^3\right)$ time, where $n$ is the size of the matrix.

If the entries of your matrix belong to an arbitrary commutative ring, then there are still $O \left(n^4\right)$-time algorithms to compute the determinant. See Günter Rote, Division-free algorithms for the determinant and the Pfaffian: Algebraic and Combinatorial Approaches, §2. (If I understand correctly, the rough idea of at least one of the algorithms is to replace the matrix $A \in R^{n\times n}$ by the matrix $1 - At$ over the power series ring $R \left[\left[t\right]\right]$, then compute the determinant of the latter via LU decomposition (which always exists in the power series ring), and then obtain $\det A$ as a coefficient of this polynomial. Of course, power series per se are uncomputable, but here only the first few coefficients need to be taken care of.)

Of course, the algorithms cannot do magic. The running time estimates of $O \left(n^3\right)$ and $O \left(n^4\right)$ assume that the fundamental operations of the base ring ($+$, $\cdot$, $-$ and taking inverses in the case of a field) require constant time and the overhead of storing and copying matrix entries is negligible. This assumption is justified if the base ring is a finite field or (to some extent) if the base "ring" is the floating-point reals (although these don't really form a ring, so you might end up with completely wrong results due to numerical instability), but not if the base ring is the integers (because integers get harder to work with the larger they become), the rational numbers or a polynomial ring. When the entries of your matrix are algebraically independent indeterminates, then you should not expect anything too fast, at least if you want the result in expanded form; after all, the result will simply be the general formula for the determinant of an $n \times n$-matrix, which "contains" a list of all $n!$ permutations of $\left\{1,2,\ldots,n\right\}$, and clearly such list requires at least $O \left(n!\right)$ time to write down! There might be some faster algorithms that result in non-expanded versions (similarly to Horner's scheme for polynomial evaluation), but I wouldn't expect anything with polynomial running time unless you allow the algorithm to return a recursion instead of an explicit sum-of-products-sums-of-products-of-etc..

Despite appearances this exercise was not a waste of time. Denote the matrix in quesation by $A_{n \times n}$. Since the determinant is a rotational invariant, for every $n\ge 2$, it is possible to express it in a closed form as a function of all other rotational invariants, i.e traces of powers of this matrix $\left( Tr[A^p] \right)_{p=1}^n$. The formula in question reads: $$\det(A) = \sum\limits_{m=1}^n \frac{(-1)^{m+n}}{m!} \cdot \sum\limits_{\begin{array}{r} p_1+p_2+\cdots+p_m=n \\ p_1\ge 1,\cdots,p_m\ge 1\end{array}} \prod\limits_{q=1}^m \frac{Tr[A^{p_q}]}{p_q}$$

Here the second sum on the right hand side runs over all partitions of $n$ into strictly positive integers.

See Calculate a multiple sum of inverse integers. for the derivation.

Denote ${\mathcal S}_n := n! \det(A)$ and $a(p) := Tr[A^p]$ for $p\ge 1$ then we have: \begin{eqnarray} {\mathcal S}_2 &=& a(1)^2-a(2)\\ {\mathcal S}_3 &=& a(1)^3-3 a(2) a(1)+2 a(3)\\ {\mathcal S}_4 &=& a(1)^4-6 a(2) a(1)^2+8 a(3) a(1)+3 a(2)^2-6 a(4)\\ {\mathcal S}_5 &=& a(1)^5-10 a(2) a(1)^3+20 a(3) a(1)^2+15 a(2)^2 a(1)-30 a(4) a(1)-20 a(2) a(3)+24 a(5)\\ {\mathcal S}_6 &=& a(1)^6-15 a(2) a(1)^4+40 a(3) a(1)^3+45 a(2)^2 a(1)^2-90 a(4) a(1)^2-120 a(2) a(3) a(1)+144 a(5) a(1)-15 a(2)^3+40 a(3)^2+90 a(2) a(4)-120 a(6)\\ {\mathcal S}_7 &=& a(1)^7-21 a(2) a(1)^5+70 a(3) a(1)^4+105 a(2)^2 a(1)^3-210 a(4) a(1)^3-420 a(2) a(3) a(1)^2+504 a(5) a(1)^2-105 a(2)^3 a(1)+280 a(3)^2 a(1)+630 a(2) a(4) a(1)-840 a(6) a(1)+210 a(2)^2 a(3)-420 a(3) a(4)-504 a(2) a(5)+720 a(7) \end{eqnarray}

It is not hard to see that in each of those formulas the coefficients in front of the product of powers of traces have the property that they sum up to zero and that their absolute values sum up to $n!$. In the table below I show the number of different terms on the right hand side $\nu(n)$ as a function of the dimension $n$. We have: \begin{array}{rr} n & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16\\ \nu(n) & 2 & 3 & 5 & 7 &11 & 15 & 22 & 30 & 42 & 56 & 77 & 101 & 135 & 176 & 231 \end{array}

Therefore it is completely not true that with this method we need to carry out $n!$ calculations. As a matter of fact for $n=49$ there are only $173525$ terms to sum up (see https://oeis.org/A000041/list).

A famous formula says determinant is product of eigenvalues. If you can find eigenvalues faster than doing such an expansion, then it could be useful

$$\det({\bf A}) = \prod_{i=1}^N \lambda_i({\bf A})$$

One way to find eigenvalues numerically is power method. You will get good estimate quickly. In terms of complexity $O(n^2)$ for matrix-vector multiplication, $n$ eigenvalues to find, that gives $O(n^3)$

All in all you will get away with a small polynomial in complexity compared to the horrendous combinatorial $O(n!)$ complexity of doing what you are doing right now.

Obviously a 2x2 matrix is easy to take a determinant of. There is a quick way to do a 3x3 as well.

$\newcommand\RED{\color{red}} \newcommand\BLUE{\color{blue}} \newcommand\GREEN{\color{green}} A= \begin{bmatrix} \ a & b & c\\ d & e & f\\ g & h & i\\ \end{bmatrix}$

Start along the top row and multiply the three diagonals shown in red, blue and green and sum them.

$\newcommand\RED{\color{red}} \newcommand\BLUE{\color{blue}} \newcommand\GREEN{\color{green}} A= \begin{bmatrix} \ \RED a & \BLUE b & \GREEN c\\ \GREEN d & \RED e & \BLUE f\\ \BLUE g & \GREEN h & \RED i\\ \end{bmatrix}$

Now ,start along the bottom row and multiply the second set of three diagonals as shown below in red, green and blue, and subtract them.

$\newcommand\RED{\color{red}} \newcommand\BLUE{\color{blue}} \newcommand\GREEN{\color{green}} A= \begin{bmatrix} \ \BLUE a & \GREEN b & \RED c\\ \GREEN d & \RED e & \BLUE f\\ \RED g & \BLUE h & \GREEN i\\ \end{bmatrix}$

The result is $det(A)=aei+bfg+cdh-gec-hfa-idb$

Essentially, this lets you just write down the answer to the determinant of any 3x3 matrix you see with minimal effort. If you know Cramers Rule it lets you immediately write down the answers to a any 3x3 system of linear equations.

Unfortunately this is a mathematical coincidence. It is NOT the case that the determinant of a square matrix is just a sum and difference of all the products of the diagonals. For a 4x4 matrix, you expand across the first column by co-factors, then take the determinant of the resulting 3x3 matrices as above. Again, if you know Cramers Rule, you can write down the answers to a 4x4 system of linear equations with only slightly more effort.

For anything larger though, it becomes absurdly complex.(use a computer, that is what they are for) There are two terms when calculating the determinant of a 2x2 matrix. There are six terms for a 3x3 matrix. For a 4x4 matrix there are 24 terms. For a 5x5 matrix, there are 120 terms. (expand by co-factors, then expand each of the 5 resulting 4x4 matrices by co-factors and then take the determinant of the resulting 3x3 matrices by diagonals. If your really careful you might be able to get through a 5x5 matrix, but a 6x6 matrix will give 720 terms and a 7x7 matrix yields 5040 terms. For an n x n, it requires n! terms! The computer uses much more efficient algorithms.