17
$\begingroup$

I have a feeling that this might have already been asked, but can't find the question. Anyway, the question is: given a symmetric $n\times n$ matrix, is there a faster way to compute its signature than computing the full eigen-decomposition? (the gt tag is because my application is in computing signatures of four-manifolds, and some others might want to do the same computation...)

$\endgroup$
4
  • 6
    $\begingroup$ Sure, Sylvester inertia law en.wikipedia.org/wiki/Sylvester's_law_of_inertia $\endgroup$ Sep 19, 2012 at 17:03
  • 2
    $\begingroup$ Can you say something about the matrix? When I needed signatures of matrices back in the 1990s, the matrices were on the order of 400 by 400. Today this is trivially small for eigensolvers. Also, if you know how big a spectral gap you have, you know how much error you can tolerate in the computations. Finally, did you mean compute as in on a von Neumann computing machine, or compute as in paper and pencil? $\endgroup$ Sep 19, 2012 at 21:14
  • $\begingroup$ @Terry: There is a chicken and egg aspect here; the hidden agenda of my question is how big matrices can I deal with... $\endgroup$
    – Igor Rivin
    Sep 20, 2012 at 0:53
  • 1
    $\begingroup$ If you have a reasonable gap in the spectrum at $0$ then I think 1000 by 1000 matrices will be no trouble, working in double floating point, even in Matlab on a basic recent desktop. If you have access to a parallelized code of inversion then the Newton's method I mention should get you up to 3000 by 3000 matrices in under an hour, on an appropriate machine with lots of RAM and cores. This is for dense matrices. For sparse matrices, I don't know. $\endgroup$ Sep 20, 2012 at 6:54

6 Answers 6

8
$\begingroup$

In the software Regina what I do is apply Descartes Rule of Signs on the characteristic polynomial. The algorithm for computing the determinant we use has complexity of the order of $n^4$, but it seems to work quite effectively for the signature computations we've been doing so far.

The primary reason we use these techniques is they're fairly easy to implement. Also, we're worried about rigour so we want to avoid floating-point data types that could have round-off errors. We use arbitrary precision integers for the computations, via the GNU Multi-Precision libraries.

$\endgroup$
6
$\begingroup$

As Fernando Muro points out in the comments, Sylvester's law of inertia is probably the easiest way to determine the signature. You diagonalize the symmetric matrix by the Gram-Schmidt process. This is essentially as easy as Gaussian elimination.

A symmetric matrix is equivalent to a bilinear form, with basis $\{v_1, \ldots, v_n\}$, and matrix entries inner products $\langle v_i,v_j\rangle, i,j\leq n$. If $\langle v_1,v_1\rangle \neq 0$, then replace $v_i$ with $v_i-v_1 \langle v_1,v_i\rangle/\langle v_1,v_1\rangle$, and proceed by induction (this has the effect of performing the same row and column operations, so as to keep the matrix symmetric). The only variation on the standard Gram-Schmidt is what to do when $\langle v_1,v_1\rangle=0$, which clearly means the form is indefinite. If $\langle v_1, v_i\rangle=0$ for all $i\leq n$, then you can proceed by induction. Otherwise, by permuting indices, you may assume that $\langle v_1,v_2\rangle \neq 0$. Replace $v_1$ with $v_2+tv_1$, such that $\langle v_2,v_2\rangle+2t\langle v_1,v_2\rangle\neq 0$, and proceed as with the usual Gram-Schmidt. Once you have diagonalized in this way, Sylvester's law implies that the signature of the diagonal matrix is the same as that of the original one.

Note added: Since you claim you're interested in computing the signatures of 4-manifolds, there's some special points one might be able to take advantage of. Certainly then the form will be integral unimodular, so presumably your symmetric matrix has integer entries and determinant $\pm 1$. If the 4-manifolds are smooth, then by Donaldson's theorem, if the intersection form is (positive) definite, then it is integrally equivalent to the diagonal form $1\oplus 1\oplus \cdots\oplus 1=1^n$. Otherwise, if the form is indefinite, then it is either odd and integrally equivalent to a diagonal form $1^k\oplus (-1)^{n-k}$, or it is even and there is a unique equivalence class for a given $n$ and signature. One can then perform the diagonalization process over $\mathbb{Z}$ using some version of the Euclidean algorithm, except in the even case (in which case there is a standard tridiagonal integral form). Actually, in the even case one could take a sum $\oplus \pm 1$ to make the form odd, and then diagonalize over the integers (a similar trick taking a sum with a hyperbolic form makes things integrally diagonalizable in the general case if you are also considering non-smooth manifolds which might have non-diagonalizable forms). See Milnor-Husemoller or Conway-Sloane for more details.

$\endgroup$
7
  • 1
    $\begingroup$ Maybe I am missing something, but you cannot just diagonalize a matrix merely by orthogonal transformations (at best one gets a tridiagonal form) --- $\endgroup$
    – Suvrit
    Sep 19, 2012 at 19:51
  • $\begingroup$ I confirm Suvrit's comment. Gram-Schmidt is not a mean for diagonalization. It is used to perform the QR factorization, which is a piece of a step of the QR algorithm. More generally, Abel's Theorem states that solving a general polynomial equation $P(x)=0$ is impossible in finitely many operations. Because you can associate a symmetric matrix with an arbitrary polynomial with real roots, you cannot diagonalize a symmetric real matrix by finitely many operations. $\endgroup$ Sep 19, 2012 at 20:17
  • 5
    $\begingroup$ But Denis, we're acting on this matrix by $A \mapsto S A S^T$, not $A \mapsto S A S^{-1}$. With the $A \mapsto S A S^T$ transformation, you can get to a diagonal matrix with only rational operations, exactly as Agol says. (I'm not sure whether or not this what you call Gram-Schmidt, although I always thought of it as Gram-Schmidt.) $\endgroup$ Sep 19, 2012 at 20:37
  • $\begingroup$ Thanks David, that was the clarification I sought; The name Gram-Schmidt invokes the idea of orthogonal transformations, whereas Agol is invoking the (constructive) a general congruence transform. $\endgroup$
    – Suvrit
    Sep 19, 2012 at 21:03
  • $\begingroup$ @ Suvrit and Denis: you're right, I shouldn't have used the term "orthogonal transformation", which I've elminated from the description. I guess I meant "orthogonalization" or something, but as David points out, it's just rational operations. $\endgroup$
    – Ian Agol
    Sep 19, 2012 at 21:41
5
$\begingroup$

If $n$ is not too big, compute the characteristic polynomial $P$ by the method of Leverrier, modified by Preparata & Sarwate; this takes $O(n^{3.5})$ operations, see the 2nd edition of my book Matrices (GTM 216, Springer-Verlag). Then apply the Routh-Hurwitz algorithm to find the number of ots positive roots. This second step, and variants, are very well described by Gantmacher in the 2nd volume of his book on matrices.

Notice that this takes finitely many operations. Actually the whole book of Gantmacher is motivated by stability issues in dynamical systems. Thus its main purpose is to check whether a given matrix has its spectrum in a given half of the complex plane.

$\endgroup$
3
  • $\begingroup$ Since we know that all the roots of the characteristic polynomial will be real, why not just use Descartes' Rule of signs to compute how many of them are positive? $\endgroup$ Sep 19, 2012 at 20:38
  • $\begingroup$ @David. Descartes' rule does not give the exact number of positive roots of a polynomial, but only an upper bound, and an information mod 2. $\endgroup$ Sep 20, 2012 at 15:37
  • $\begingroup$ @Denis: For polynomials with all roots real, Descartes' rule gives the exact number of roots. $\endgroup$ Nov 17, 2012 at 9:42
4
$\begingroup$

Here is a textbook solution that runs in time $O(n^3)$ (assuming that Gaussian Elimination can be done in that time, which might not be immediate unless one takes steps to rule out intermediate coefficient growth). The first part of the statement below is merely LU decomposition. The last sentence is what you're looking for.

Theorem (Stewart) Let $A$ be symmetric and suppose that the leading principal submatrices of $A$ are nonsingular. Then, there exists a unit lower triangular matrix $L$ and upper triangular matrix $U$ such that $A=LU$. Moreover, if $D_U=\text{diag}(u_{11},\ldots,u_{nn})$, then $\text{Inertia}(D_U)=\text{Inertia}(A)$.

Unfortunately, if one pivots (which in practice one must), the above rule breaks. In this case, I would recommend first tridiagonalizing the matrix using orthgonal transformations---this will cost $O(n^3)$. Then, you can compute the eigenvalues of the tridiagonal matrix in time $O(n^2)$ (pseudo-polynomial time, though), and then you'll have the inertia.

$\endgroup$
3
$\begingroup$

This seems to be an $O(n^3)$ problem, so if your matrix is smaller than 100 by 100 you might as well use an eigensolver.

For bigger matrices you can use Newton's method (see Higham's papers or book) to get a real almost orthogonal matrix and take the trace and round off. To get near the orthogonal in the polar decomposition iteratively replace $A$ by the average of $A$ with the transpose of the inverse of $A$.

That is, let $A_0 = A$ and $A_{n+1} = \left( \frac{1}{2}A + \frac{1}{2}\left( A^{\mathrm{T}} \right)^{-1} \right) $ and then the signature is $$ \lim_{n \rightarrow \infty} \mathrm{Trace} \left( A_N \right ) $$

If you have a priori information on the norm of $A$ and the size of the spectral gap at zero, you can figure how many iterations to do easily. If not many, this is very fast.

Nicholas J. Higham and Pythagoras Papadimitriou, A parallel algorithm for computing the polar decomposition, Elsevier, Parallel Computing, 20, (1994).

$\endgroup$
2
$\begingroup$

May I attract your attention on the following formula, where $B$ is a $n\times n$ non-singular symmetric real-valued matrix $$ \int_{\mathbb R^n}e^{-2i\pi x\cdot \xi} e^{i\pi Bx\cdot x}dx=e^{\frac{i\pi}{4}\text{signature } B}\vert\det B\vert^{-1/2}e^{-i\pi B^{-1}\xi\cdot \xi}, $$ so that $ e^{\frac{i\pi}{4}\text{signature } B}=\vert\det B\vert^{1/2} \int_{\mathbb R^n} e^{i\pi Bx\cdot x}dx. $ Of course the Fresnel-type integral on the rhs of last formula is not absolutely converging but, as an oscillatory integral, will become so after integration by parts.

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