TL;DR - I computed the one-body density matrix (OBDM) stochastically via a method in a paper listed below, and it generates non-physical occupation numbers that 1) either has some negative values or 2) has some values greater than 1. I'm just wondering what's the best course of action to resolve this issue. Any help is greatly appreciated.


I recently read a paper by P. Lowdin (Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural Spin-Orbitals, and Convergence Problems in the Method of Configurational Interaction ), which states a stochastic method to generate the one-body density matrix for an arbitrary wavefunction.

I've used this method to generate the one-body density matrix (OBDM), and it appears to resemble the one-body density matrix from other methods (e.g. FCI), which is good, however, I've noted that when I compute the occupation numbers of the OBDM I get negative values, which make no sense.

I did check if the OBDM is PSD, which it turns out it isn't and this makes sense from the negative eigenvalues, which I assume emerges as an artifact from calculating the OBDM stochastically.

This led me to symmetrize the OBDM via replacing, $M \to M \cdot M^\intercal$, where $M$ is the OBDM of size $N \times N$ and $\intercal$ is simply the transpose operator. This should ensure that the OBDM is symmetric and PSD. I calculated occupation numbers for this 'symmetrized OBDM' and all occupation numbers are positive, but some occupation numbers are greater than 1, which does not make sense physically at all as the occupation numbers are bound to be $0 \leq \lambda \leq 1$ (at least physically).

The question I have is, how can I resolve this issue with non-physical occupation numbers for the OBDM? Were the original occupation numbers valid (even though some were negative?), because in that case all values were bounded below 1. However, when I compute something like the von Neumann entropy, I get an invalid value that's above the maximum entropy for the system, which I've computed as $\ln(N)$ with $N$ being the number of particles for my system in 1D-space.

Any help on this would be greatly appreciated!

  • $\begingroup$ Yes, it's not the best notation. I'll make an edit to it. In short, I'm multiplying it by its transpose, then re-defining that to itself. $\endgroup$ Aug 31, 2023 at 16:58
  • 1
    $\begingroup$ A stochastic method to calculate a matrix will give statistical errors. These will give a bias for the computed eigenvalues. You should estimate the statistical errors and the bias. Similar problems are addressed here: D. M. Ceperley; B. Bernu, J. Chem. Phys. 89, 6316–6328 (1988), The calculation of excited state quantum Monte Carlo, doi.org/10.1063/1.455398 $\endgroup$
    – user200143
    Aug 31, 2023 at 17:14
  • $\begingroup$ Thanks for sharing this paper @user200143. I've done a quick estimate of the error of the occupation numbers, and the negative values do have an error. I've estimated the error as the root mean squared error of T values for the occupation numbers. As I don't have an 'exact' occupation number, would it suffice if I took the standard deviation of the T values and divided it by sqrt(T)? As in, compute an error of the mean of those samples. to define an error for the occupation numbers? $\endgroup$ Aug 31, 2023 at 20:17

1 Answer 1


This was too long for a comment, so I wrote it as an answer. Yes, repeating the calculation T times and calculating the average and standard deviation and using the central limit theorem like you say will give a statistical error estimate. However, the result will be biased. You can understand this from a 2x2 matrix with zero off diagonal elements with statistical errors. Either sign of the fluctuation will split the levels, so the statistical errors in the off diagonal make the average of the lower eigenvalue smaller and the average of the higher eigenvalue bigger than the correct values. Fortunately bias drops off as $T^{-1}$, (i.e. keep more samples in the matrix element averages) so it's usually possible to make it smaller than the statistical errors by increasing the number of samples.

Edited to add:

By splitting I meant for the 2x2 case where for equal diagonal elements $d$, and off diagonals $o$, the eigenvalues are $d \pm o$, so the original degenerate eigenvalues are split.

For the bias, look at calculating a function (like an eigenvalue) of the average of a set of stochastic variables. For the simplest one variable case, you have $f(\bar x +\Delta x)$ where $\bar x$ is the correct average and $\Delta x$ is a fluctuating error. Taking the average and assuming $\Delta x$ is small you get \begin{equation} \langle f(\bar x +\Delta x)\rangle = \left \langle f(\bar x) + \Delta x f'(\bar x) + \frac{1}{2}\Delta x^2 f''(\bar x)+... \right \rangle \end{equation} Since $\langle \Delta x \rangle = 0$, the bias is approximately $\left \langle \frac{1}{2}\Delta x^2 f''(\bar x) \right \rangle$ which is propo rtional to the variance which drops off like $T^{-1}$. For more variables, the correlations between them also contribute but the lowest order terms are always the average of two fluctuating quantities with the same $T^{-1}$ behavior.

The statistical do fall off like $T^{-1/2}$, i.e. they are given by the square root of the variance divided by $T$, where the variance is $\langle [f(x+\Delta x)]^2\rangle-\langle f(x+\Delta x)\rangle^2$ so usually bias is easier to control by taking more samples of $x$ before using it in whatever $f(x)$ you are calculating.

  • $\begingroup$ Thanks for the explanation @user200143, I assume the error emerging comes from the off-diagonal elements because the on-diagonal matrix elements are strictly positive (so no sign problem) whereas the off-diagonal elements are positive/negative and hence have a sign problem. When you say "spilt the levels" are you referring to the eigenvalues of the OBDM being separated by their sign? Additionally, is there a way to estimate the statistical noise of this method? I'd assume it's the RMSE of the matrix elements of the error the $T$ OBDMs? Also why does it decay with $T^{-1}$ and not $T^{-1/2}$? $\endgroup$ Sep 1, 2023 at 17:05
  • 1
    $\begingroup$ I edited the answer since again my comment would not fit. $\endgroup$
    – user200143
    Sep 3, 2023 at 3:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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