For positive definite $A,B$ why does $AB+BA$ tend to be positive definite?

Question

Let $A$ and $B$ be two positive definite $n \times n$ matrices. It is, of course, not true that $AB+BA$ is necessarily positive definite.

Consider, though, the results of the following numerical experiment. I generated $A$ by letting its eigenvalues be random in $[0,1]$, and selecting its eigenvectors by generating a random matrix of standard Gaussians and applying Gram-Schmidt to it. The matrix $B$ is generated in the same way.

I did this 1000 times and checked what proportion of times the matrix $AB+BA$ has at least one negative eigenvalue [1]. Here are the results for different dimensions $n$:

• $n=2, ~~~~94.8 \%$
• $n=3, ~~~~89.4 \%$.
• $n=4, ~~~~78 \%$.
• $n=5, ~~~~72.7 \%$.
• $n=10, ~~~40.3 \%$.
• $n=15, ~~~20.1 \%$.
• $n=20, ~~~11.4 \%$.
• $n=50, ~~~0.3\%$.
• $n=100, ~~0 \%$.

This suggests that, as a function of $n$, examples with $AB+BA$ not psd tend to get rarer and rarer. Is it possible to give a proof of this?

It may be more natural to consider a different random model of randomly generated psd matrices; I only generated them in the way I described above because it seemed easiest.

[1] Actually, I checked if there is an eigenvalue less then $-1 \cdot 10^{-5}$ in MATLAB to account for rounding error.

Answer 1

$\text{tr}(AB+BA) = 2 \operatorname{tr}(A^{1/2} B A^{1/2}) > 0$, so that may produce some bias toward positive eigenvalues. In particular if you generate your "random" matrices in such a way that the eigenvalues of $AB+BA$ will tend to be concentrated very close together, this may produce the results you observed.

But I tried a different experiment: $A = X^T X$ and $B = Y^T Y$ where $X$ and $Y$ are random $n \times n$ matrices with integer entries in $[-100,100]$.
For the case $n=10$, I found that it was very rare (0 occurrences in 3000 trials) for $AB + BA$ to be positive definite.

Answer 2

Your question appears to be based on a false premise. In fact $AB+BA$ does not tend to be positive definite as $n$ increases, even within the particular distribution you happen to be using.

To demonstrate this, here is a simple piece of Mathematica code that implements precisely the numerical experiment you described:

randMat[n_] := With[{
    eigval = RandomReal[{0, 1}, n],
    eigvec = Orthogonalize @ RandomVariate[NormalDistribution[], {n, n}]},
  Transpose[eigvec] . DiagonalMatrix[eigval] . eigvec];

prob[n_] := Table[With[{A = randMat[n], B = randMat[n]},
  PositiveDefiniteMatrixQ[A.B + B.A]], {10000}];

The results obtained (from running 10,000 trials for each value of $n$) are as follows:

+-----+-------------------------------+-----------------------------------+
| n   | probability that AB+BA is pd  | probability that AB+BA is NOT pd  |
+-----+-------------------------------+-----------------------------------+
| 2   | 94.67%                        | 05.33%                            |
+-----+-------------------------------+-----------------------------------+
| 3   | 87.09%                        | 12.91%                            |
+-----+-------------------------------+-----------------------------------+
| 4   | 79.57%                        | 20.43%                            |
+-----+-------------------------------+-----------------------------------+
| 5   | 71.26%                        | 28.74%                            |
+-----+-------------------------------+-----------------------------------+
| 10  | 39.94%                        | 60.06%                            |
+-----+-------------------------------+-----------------------------------+
| 15  | 21.78%                        | 78.22%                            |
+-----+-------------------------------+-----------------------------------+
| 20  | 10.31%                        | 89.69%                            |
+-----+-------------------------------+-----------------------------------+
| 50  | 00.16%                        | 99.84%                            |
+-----+-------------------------------+-----------------------------------+
| 100 | 00.00%                        | 100.00%                           |
+-----+-------------------------------+-----------------------------------+

From this it is clear that you must have accidentally switched the two columns. As $n$ increases, it is becomes rarer for $AB + BA$ to be positive definite.

As an aside, I find it surprising that this question and its currently accepted answer have received a combined 17 upvotes, when seemingly nobody has even tried to replicate the OP's trivial numerical experiment.

Answer 3

It is in theory possible to use free probability to describe the limit eigenvalue distribution (as the dimension tends to $+\infty$) of the anticommutator $AB+BA$ when $A$ and $B$ are independent random matrices which are unitarily invariant (which is the case for the description used in the OP). The same is true actually for any self-adjoint polynomial in $A$, $B$.

The drawback of this method it that the limit distribution will be obtained in an indirect way. If one does this for the uniform measure on $[0,1]$ as asked in the OP, it may be hard to decide from this approch whether the limit distribution is supported on $\mathbf{R}^+$ or not.

However, such computations are possible in simple cases. Propositon 6.11 from Fevrier, Maxime; Mastnak, Mitja; Nica, Alexandru; Szpojankowski, Kamil, Using Boolean cumulants to study multiplication and anti-commutators of free random variables says the following: if $E$ and $F$ are subspaces of dimension $n/2$ which are chosen independently according to the $O(n)$-invariant measure on the Grassmann manifold, and $A$ and $B$ are the corresponding orthogonal projections, then as $n \to \infty$ the eigenvalue distribution of $AB+BA$ is described by an explicit density which is supported on $[-1/4,2]$ (see also Figure 2 in the same paper). In particular, the probability that $AB+BA$ is positive semidefinite tends to $0$ as $n \to \infty$. For this toy model, this gives a rigorous proof of the fact that the anticommutator of positive matrices is typically not positive. (This is probably overkill, ane can presumably analyze everything in terms of principal angles between the subspaces $E$ and $F$.)