8
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ What goes wrong with Gaussian elimination? Is it numerically unstable? $\endgroup$ Aug 13, 2012 at 18:44
  • 1
    $\begingroup$ 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). $\endgroup$
    – Igor Rivin
    Aug 13, 2012 at 19:06
  • $\begingroup$ Are the floating point entries rational, real algebraic, or other? $\endgroup$
    – Will Jagy
    Aug 13, 2012 at 20:19
  • $\begingroup$ @Will: What exactly do you mean by your question? $\endgroup$
    – Igor Rivin
    Aug 13, 2012 at 21:40
  • $\begingroup$ @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. $\endgroup$
    – Will Jagy
    Aug 14, 2012 at 2:41

3 Answers 3

4
$\begingroup$

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....

$\endgroup$
3
$\begingroup$

(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.

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Igor Rivin
    Aug 13, 2012 at 21:49
  • $\begingroup$ Fixed the terminology. $\endgroup$ Aug 13, 2012 at 22:00
  • 1
    $\begingroup$ 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. $\endgroup$ Aug 13, 2012 at 23:22
  • $\begingroup$ 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). $\endgroup$
    – user9072
    Aug 14, 2012 at 0:07
  • $\begingroup$ @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. $\endgroup$ Aug 14, 2012 at 0:31
3
$\begingroup$

To just add to the answers already here, it is possible to extract the individual coefficients of this polynomial exactly using the following method. Let

\begin{equation} h(x) = \sum_{j = 0}^d c_j x^j \end{equation}

be a polynomial of degree $d$, then

\begin{equation} c_j = \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}}) \end{equation}

since

\begin{align} \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}}) &= \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} \left(\sum_{l = 0}^d c_l e^{\frac{2 \pi i k l}{d + 1}} \right) \\ &= \sum_{l = 0}^d \left(\frac{1}{d + 1} \sum_{k = 0}^d e^{\frac{2 \pi i (l - j) k}{d + 1}} \right) c_l \\ &= \sum_{l = 0}^d \delta_{jl} c_l \\ \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}}) &= c_j \end{align}

where $\delta_{jl}$ is the Kronecker delta and $i$ is the imaginary unit. Letting

\begin{equation} f(x, y, z) = \det(x A + y B + z C) \end{equation}

we can expand $f$ in powers of $x, y, z$ using multilinearity of the rows or columns of the matrix inside the determinant as

\begin{equation} f(x, y, z) = \sum_{j = 0}^n x^j g_j(y, z) = \sum_{j = 0}^n x^j \left(\sum_{k = 0}^{n - j} a_{j, k} y^k z^{n - j - k}\right) \end{equation}

Note that, since one choice of $x, y, z$ must be made for every row or column in the expansion, the powers of these variables must sum to $n$ in every term. This polynomial is completely specified by the coefficients $\{a_{j, k}\}$, and to extract them, it is sufficient to only evaluate $f$ at $z = 1$

\begin{equation} f(x, y, 1) = \sum_{j = 0}^n x^j g_j(y, 1) = \sum_{j = 0}^n \sum_{k = 0}^{n - j} a_{j, k} x^j y^k \end{equation}

One query of the polynomial

\begin{equation} g_j(y, 1) = \sum_{k = 0}^{n - j} a_{j, k} y^k \end{equation}

is given by the above formula

\begin{equation} g_j(y, 1) = \frac{1}{n + 1} \sum_{p = 0}^n e^{-\frac{2 \pi i j p}{n + 1}} f(e^{\frac{2 \pi i p}{n + 1}}, y, 1) \end{equation}

Similarly,

\begin{equation} a_{j,k} = \frac{1}{n - j + 1} \sum_{q = 0}^{n - j} e^{-\frac{2 \pi i k q}{n - j + 1}} g_j(e^{\frac{2 \pi i q}{n - j + 1}}, 1) \end{equation}

Therefore

\begin{equation} \det(x A + y B + z C) = \sum_{j = 0}^n \sum_{k = 0}^{n - j} a_{j,k} x^j y^k z^{n - j - k} \end{equation}

where

\begin{equation}\ a_{j, k} = \frac{1}{(n + 1)(n - j + 1)}\sum_{p = 0}^{n} \sum_{q = 0}^{n - j} e^{-2 \pi i \left(\frac{jp}{n + 1} + \frac{kq}{n - j + 1} \right)} \det(e^{\frac{2 \pi i p}{n + 1}} A + e^{\frac{2 \pi i q}{n - j + 1}} B + C) \end{equation}

requiring the evaluation of only $O(n^4)$ determinants to completely characterize the polynomial.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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