# Fast method to calculating signature of a matrix

For calculation of the topological properties of a hamiltonian, sometimes we need the signature of that matrix. This means we only need number of positive eigenvalues. One simple way is to first calculate the eigenvalues and then find difference of number of positive and negative eigenvalues. for example,

H = RandomReal[{-1, 1}, {100, 100}] // (# + #\[ConjugateTranspose]) &;
Total[If[# > 0, 1, -1] & /@ Eigenvalues[H]]


I am interested to know i) is there any function in Mathematica calculate signature ii) is there any fast algorithm within Mathematica to do so fast?

• You mean Sign? This function is Listable i.e. Sign[{-3, 4, 0}] works BTW. – xzczd Jun 30 '20 at 3:44
• Total@Sign@Eigenvalues[H] – SneezeFor16Min Jun 30 '20 at 4:49
• The problem is not in sign. The problem is how to get rids of the eigenvalue calculation. – Rasoul-Ghadimi Jun 30 '20 at 5:37
• @Rasoul-Ghadimi First of all, there's no function in Mathematica that directly computes the signature of a matrix. Furthermore, all the definitions for signature I can find lead to eigenvalues. Unless there're definitions that avoid eigenvalues (e.g. ask on Mathematics SE), I don't think there's any faster way to do the calculation. – SneezeFor16Min Jun 30 '20 at 6:52
• – Michael E2 Jun 30 '20 at 17:19

SeedRandom[1];
H = RandomReal[{-1, 1}, {100, 100}] // (# + #\[ConjugateTranspose]) &;
Short[H, 3]

{{1.269557961,-0.6813805661,<<96>>,-0.3068229219,-0.9273213366},<<98>>,{-<<19>>,<<99>>}}


# Comparison

Eigenvalues is the fastest:

Total@Sign@Eigenvalues[H] // RepeatedTiming

{0.00044, 2}


Descartes' Sign Rule on the characteristic polynomial:

Total[1 - Ratios[Sign@CoefficientList[
CharacteristicPolynomial[H, x], x]
]] - MatrixRank[H] // RepeatedTiming

Total[RealAbs@Differences@Sign@CoefficientList[
CharacteristicPolynomial[H, x], x]] -
MatrixRank[H] // RepeatedTiming

2 (Length[Split@Sign@CoefficientList[
CharacteristicPolynomial[H, x], x]
] - 1) - MatrixRank[H] // RepeatedTiming

{0.0030, 2}
{0.0031, 2}
{0.0032, 2}


$$LDL^\top$$ decomposition (idea, code from @J.M.'stechnicaldifficulties):

Total@Sign@LDLT[H][[2]] // RepeatedTiming

{0.0036, 2}


## Benchmark

Needs["GeneralUtilities"];
BenchmarkPlot[{sig1Eigen, sig2Poly1Ratio, sig2Poly2Diff,
sig2Poly3Split, sig3LDLT},
n \[Function] StatisticsLibraryVectorToSymmetricMatrix[
#[[n + 1 ;;]], #[[;; n]], n
] &@RandomReal[{-1, 1}, Binomial[n + 1, 2]],
"IncludeFits" -> True, TimeConstraint -> 100]


Definitions are the same as above. (Here real symmetric matrices are generated with an undocumented function seen here, which is a bit faster.) Result:

It seems that all of these methods are in $$\mathcal{O}(n^3)$$ (Ratios one should be the same), but Eigenvalues one has the smallest coefficient.

## Note: Bug?

One strange thing is that CountRoots doesn't give correct answers here, nor does Reduce or Solve. Is this a bug?

CountRoots[CharacteristicPolynomial[H, x], {x, 0, \[Infinity]}]

21 (* Should be 51 *)

Reduce[CharacteristicPolynomial[H, x] == 0 && x > 0, x] // Length (* or Solve *)

15 (* Should be 51 *)

Solve[{CharacteristicPolynomial[H, x] == 0, x > 0}, x,
Complexes] // Length

51 (* Correct *)

Reduce[CharacteristicPolynomial[H, x] == 0, x] // Length (* or Solve *)

100 (* Correct *)


Otherwise, 2 CountRoots[CharacteristicPolynomial[H, x], {x, 0, \[Infinity]}] - MatrixRank[H] can be used.

### Update

This must be a bug -- Some real-valued roots are considered by CountRoots to have "invisible" imaginary parts, with some of them "$$| Im(x) | \gt 1$$"!

CountRoots[
CharacteristicPolynomial[H, x], {x, -I, 100 + I}] // AbsoluteTiming

{44.9898317, 33} (* Still incorrect *)

CountRoots[
CharacteristicPolynomial[H, x], {x, -5 I,
100 + 5 I}] // AbsoluteTiming

{127.967137, 51} (* Correct, but very slow *)


Solve has slighter peoblem. I think it's due to the machine precision:

Solve[CharacteristicPolynomial[H, x] == 0 && Re[x] > 0, x] //
Length // AbsoluteTiming

{0.0206228, 51} (* Correct *)


### Update 2

According to @MichaelE2, this is because the machine precision is not enough for the large coefficients and high degree of the characteristic polynomial.

• It would be interesting to know how the methods scale with the matrix size $n$. What method is the asymptotic winner for $n\rightarrow\infty$? – yarchik Jun 30 '20 at 21:25
• @yarchik, I've updated the answer. – SneezeFor16Min Jul 1 '20 at 5:54
• Nice analysis. I interpret it as the rank-determination is the bottleneck with $\mathcal{O}(n^\omega)$. – yarchik Jul 1 '20 at 6:54
• @yarchik Yes, and also the characteristic polynomial, maybe with $\mathcal{O}\left( n^\omega \log n \right)$. – SneezeFor16Min Jul 1 '20 at 8:02
• Probably not a bug. Machine-precision, high-degree polynomials are known to be perfidious. Examine Plot[Evaluate@CharacteristicPolynomial[H, x], {x, 0, 10}]`. – Michael E2 Jul 1 '20 at 21:52