Fast method to calculating signature of a matrix

Question

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?

Answer 1

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] Statistics`Library`VectorToSymmetricMatrix[
     #[[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.