Computing signature

Question

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

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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

Answer 6

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.