# Computing determinants of matrices of linear forms

Suppose we have three $n \times n$ matrices $A$, $B$, $C$ with floating point entries. We would like to compute the polynomial $\det (xA+yB+zC)$. At least in Mathematica, and I think in all computer algebra systems, this will take $n!$ steps; Mathematica chokes around $n=15$.

What's the smart way to do this? I have some ideas, but they are all complicated enough that I don't want to implement them before hearing from others.

This question is on the border between MO, cstheory and Mathematica, but I suspect that this is the right place to start.

ADDED IN RESPONSE TO QUESTIONS BELOW The trouble with Gaussian elimination is that you have to divide by polynomials in $(x,y,z)$, and the expressions soon get huge. Interpolation, probably at roots of unity so that the interpolation matrix will be unitary, is my best idea, but I wanted to see if there was a better one before I tried it.

-
What goes wrong with Gaussian elimination? Is it numerically unstable? – Qiaochu Yuan Aug 13 '12 at 18:44
The problem is that the matrix is not a numeric matrix, so you will wind up with horribly huge rational functions (with floating point coefficients). – Igor Rivin Aug 13 '12 at 19:06
Are the floating point entries rational, real algebraic, or other? – Will Jagy Aug 13 '12 at 20:19
@Will: What exactly do you mean by your question? – Igor Rivin Aug 13 '12 at 21:40
@Igor, mostly if they are rational and a common denominator can be found, I get a different sense of this. You can do Gaussian elimination over the integers, and over polynomials with integer coefficients, by keeping track of a steadily growing common divisor. – Will Jagy Aug 14 '12 at 2:41

The smart way of doing this is by interpolation (this, by the way, is how Mathematica computes the characteristic polynomial), and in Mathematica you can use InterpolatingPolynomial[] to do this in a couple-of-line program.

Remark In fact, to compute the characteristic polynomial of an integer matrix, Mathematica uses interpolation twice, in effect. That's because to compute the ordinary determinant, it uses Chinese remaindering....

-

(1) Look up Gauss-Bareiss and Dodgson condensation - the determinant can be computed fraction-free over a PID.

(2) Try interpolation. I've used this in anger for the two-variable case, it does look less pleasant with three variables.

-
First, it is Dodgson CONDENSATION. Second, what is so unpleasant about interpolation? I would conjecture that it would be MUCH faster than the other algorithms you suggest, and the implementation (at least in mathematica) is completely trivial. – Igor Rivin Aug 13 '12 at 21:49
Fixed the terminology. – Chris Godsil Aug 13 '12 at 22:00
Confused about suggestion (1). $\mathbb{R}[x,y,z]$ is not a PID (though it is a UFD) and every version of Dodgson condensation I know requires division. – David Speyer Aug 13 '12 at 23:22
But it requires no 'true' division; all division are a priori known to be exact. And PID is not really relevant. What one need to be able to do is: given A and B where it is known that A divides B in the ring, compute (an approximation to) the co-divior. I am not completely sure now, but this seems feasible also for multi var poly (though perhaps it is too expensive to keep the algo competitive relative to approximation). – user9072 Aug 14 '12 at 0:07
@David Speyer: what quid says. (Dodgson will work over a UFD, but PID was as much as I felt safe asserting from memory.) I suspect Igor is right about the relative speed, but both algorithms are easy to implement. My impression is that if the number of variables increases, Dodgson may become more attractive. – Chris Godsil Aug 14 '12 at 0:31