# Computing signature

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

-
Sure, Sylvester inertia law en.wikipedia.org/wiki/Sylvester's_law_of_inertia –  Fernando Muro Sep 19 '12 at 17:03
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? –  Terry Loring Sep 19 '12 at 21:14
@Terry: There is a chicken and egg aspect here; the hidden agenda of my question is how big matrices can I deal with... –  Igor Rivin Sep 20 '12 at 0:53
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. –  Terry Loring Sep 20 '12 at 6:54

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.

-

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.

-
Maybe I am missing something, but you cannot just diagonalize a matrix merely by orthogonal transformations (at best one gets a tridiagonal form) --- –  Suvrit Sep 19 '12 at 19:51
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. –  Denis Serre Sep 19 '12 at 20:17
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.) –  David Speyer Sep 19 '12 at 20:37
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. –  Suvrit Sep 19 '12 at 21:03
@ 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. –  Ian Agol Sep 19 '12 at 21:41

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.

-
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? –  David Speyer Sep 19 '12 at 20:38
@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. –  Denis Serre Sep 20 '12 at 15:37
@Denis: For polynomials with all roots real, Descartes' rule gives the exact number of roots. –  Markus Schweighofer Nov 17 '12 at 9:42

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.

-

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.

-

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

-