8
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ You mean Sign? This function is Listable i.e. Sign[{-3, 4, 0}] works BTW. $\endgroup$ – xzczd Jun 30 at 3:44
  • 1
    $\begingroup$ Total@Sign@Eigenvalues[H] $\endgroup$ – SneezeFor16Min Jun 30 at 4:49
  • $\begingroup$ The problem is not in sign. The problem is how to get rids of the eigenvalue calculation. $\endgroup$ – Rasoul-Ghadimi Jun 30 at 5:37
  • 1
    $\begingroup$ @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. $\endgroup$ – SneezeFor16Min Jun 30 at 6:52
  • 2
    $\begingroup$ Related: mathoverflow.net/questions/107593/computing-signature $\endgroup$ – Michael E2 Jun 30 at 17:19
12
$\begingroup$
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:

Benchmark for Matrix Signatures

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.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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$? $\endgroup$ – yarchik Jun 30 at 21:25
  • 1
    $\begingroup$ @yarchik, I've updated the answer. $\endgroup$ – SneezeFor16Min Jul 1 at 5:54
  • 1
    $\begingroup$ Nice analysis. I interpret it as the rank-determination is the bottleneck with $\mathcal{O}(n^\omega)$. $\endgroup$ – yarchik Jul 1 at 6:54
  • 1
    $\begingroup$ @yarchik Yes, and also the characteristic polynomial, maybe with $\mathcal{O}\left( n^\omega \log n \right)$. $\endgroup$ – SneezeFor16Min Jul 1 at 8:02
  • 1
    $\begingroup$ Probably not a bug. Machine-precision, high-degree polynomials are known to be perfidious. Examine Plot[Evaluate@CharacteristicPolynomial[H, x], {x, 0, 10}]. $\endgroup$ – Michael E2 Jul 1 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.