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. 



For positive definite $A,B$ why does $AB+BA$ tend to be positive definite?
$\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.


