Calculating the trace with Haar random pure states

Maximally mixed state

A ubiquitous quantum state is the maximally mixed state,

$$\tau = \frac{I}{d},$$

where $d$ is the total Hilbert space size and $I$ is the identity operator. It is a mixed quantum state where every state is given equal weight. You can arrive at a maximally mixed state in a variety of ways, like taking the infinite temperature limit of a Gibbs state, or looking at the reduced density matrix of a maximally entangled pure quantum state. To better understand the maximally mixed state, we want to take expectation values of a general operator expression, $\hat{A}$,

$$\langle \hat{A} \rangle = \frac{1}{d} \text{tr} \left( \hat{A} \right)$$

Numerically calculating the trace

where $\hat{A}$ could be a Hamiltonian, a time evolved observable, a quantum circuit or any other operator expression for which the trace is defined. Calculating the trace is numerically expensive. For this tutorial let's assume we want to calculate the trace exactly (or within a well controlled error), and that there isn't anything special about $\hat{A}$ that would allow us to symbolically or analytically reduce the complexity of calculating the trace. The first, most memory intensive way to calculate the trace could be casting our operator expression $\hat{A}$ to a matrix,

/*
    Build A as some expression, maybe an operator_prod
    or an operator_sum
*/
var mat_A = matrix(A)

print(trace(mat_A))

There are a few problems with this. If we are discussing a spin $\frac{1}{2}$ system with $N$ qubits, then we are storing $2^{2N}$ real or complex numbers, restricting us to $N=16$ or so spins. If we don't want to pay this cost of memory we could instead move our cost to complexity,

var psi = zeros([2**N], as_matrix);
var sum = 0.0;
for(var m = 0; m<2**N; ++m)
{
    psi[m] = 1;
    sum+= real(expval(A,psi));
    psi[m] = 0;
}
print(sum);

If $\hat{A}$ was a summation of non-branching operators you could also use Qbit and the >> in place operation to completely avoid expensive state vector storage and simply act operators in place on a Qbit. Regardless, this strategy is extremely expensive in terms of complexity. I encourage you to try it.

Approximating traces with Haar random states

We can exploit something called quantum typicality to approximate trace expressions with haar random states,

$$\frac{1}{d} \text{tr} \left( \hat{A} \right) = \langle \psi | \hat{A} | \psi \rangle + O\left( \frac{1}{\sqrt{d}} \frac{1}{d} \text{tr} \left( \hat{A}^\dagger \hat{A} \right) \right),$$

where $|\psi\rangle$ is a Haar random pure state. For sufficiently large Hilbert space dimension typically one vector is sufficient to approximate the trace to high precision. It is crucial to work in contexts where $\frac{1}{d} \text{tr} \left( \hat{A}^\dagger \hat{A} \right)$ is $O(1)$. This is a standard result in quantum information theory, for more information you can consult typicality results like those referenced in section 3 of this article.

In aleph spawning a Haar random vector is easy,

var N = 18;
var dim = 2**N;
// the vector can be any size you want
// it doesn't need to be a power of 2. 
var my_vec = haar_vector(dim);

Calculating correlation functions

To give us a concrete example to use typicality let's first write out an example Hamiltonian. For this I will choose a non-integrable Ising chain written below,

$$\hat{H} = \sum_{m=0}^{L-1} J_1 \gamma_1 Z_m Z_{m+1} + J_2+ \gamma_2 Z_m Z_{m+2} + \lambda_z Z_m + \lambda_x X_m.$$

And for simplicity we will calculate the autocorrelation function,

$$C(t) = \frac{1}{d} \text{ tr } \left( Z_0(t) Z_0(0) \right)$$

which is an expectation value of a correlation function taken with the maximally mixed state. Alternatively you could identify this with the infinite temperature ensemble. The point here really isn't the physical context, it is accurately approximately $\frac{1}{d} \text{tr} \left( \cdot \right)$ using quantum typicality (a Haar random pure state). In this example we will utilize Krylov time evolution which was covered here.

// make our model with open boundaries
// open boundaries are often better to allow information
// to propagate without encountering finite size effects
def make_hamiltonian(real J1, real J2, real lambdax, real lambdaz, integer L)
{
    var ham = operator_sum(as_complex)
    var s1; var s2;
    var s3;
    // we are doing open boundaries
    for(m : [0..L-1])
    {
        s1 = m;
        s2 = (m+1)%L
        s3 = (m+2)%L
        ham += J1*ZZ(s1,s2) + J2*ZZ(s1,s3) + lambdax*X(s1) + lambdaz*Z(s1)
    } 
    return ham;
}

// chosen constants for the example
var J1 = 1.0;
var J2 = -0.5;
var lambdax = 0.5;
var  lambdaz = -0.5;
var L = 20;

var ham = make_hamiltonian(J1,J2,lambdax,lambdaz,L);

// making our symbolic unitary
var dt = 0.1;
var Nkrylov = 10;
var expiHt = fun[dt](x) { return exp(-1i * dt * x) }
var U = operator_function(ham, expiHt, ["krylov_dimension" : Nkrylov]);

// we want to do the calculation as, 
// <psi(t) | Z_0 | zpsi(t)> 
var psi = haar_vector(2**L);
var zpsi = Z(0)*psi; 
var psi_temp = Z(0)*zpsi;
print(dot(psi,psi_temp));
// doing 10 steps to showcase
for(steps:[0..10])
{
    psi = U * psi;
    zpsi = U * zpsi;
    psi_temp = Z(0)*zpsi;
    print(dot(psi,psi_temp));

}

What's next?

• Write a brute force exact diagonalization time evolution verifying that this example does approximate our autocorrelation function
• Can you make this script faster?
• Is our Krylov dimension too small? Too big?
• What other trace expressions would you use this for?

Documentation Contributors

Jonathon Riddell