Expectation values of Pauli strings?

[BoltzmannEntropy]

Hello, For a tutorial, I am trying to do something very simple, calculating the expectation value of: Instead of:

hamiltonian = random_pauli_str_generator(N, terms=1)

How can i specify my own Hamiltonian to compare the value with the theoretical one? Thanks

[LeiZhang-116-4]

One naive example for the expectance calculation could be

from paddle_quantum import Hamiltonian
from paddle_quantum.state import random_state

psi = random_state(2)
pauli_str = [(1.0, 'Z0, Z1'), (2.0, 'X0, Y1')]
H = Hamiltonian(pauli_str)
print("The expectation value of the observable ZZ + 2XY is", psi.expec_val(H))

Above example cannot be used in a QNN, in which case you should try

from paddle_quantum.loss import ExpecVal

loss_fcn = ExpecVal(H)
print("The loss value is", loss_fcn(psi))

[BoltzmannEntropy]

Thank you very much. For teaching purposes, is this the preferred method?

%reset -f
import numpy as np
from numpy import pi as PI
import paddle 
import paddle_quantum
from paddle import matmul
from paddle_quantum.ansatz import Circuit
from paddle_quantum.qinfo import random_pauli_str_generator, pauli_str_to_matrix
from paddle_quantum.linalg import dagger
from paddle_quantum import Hamiltonian
from paddle_quantum.state import random_state

x0=np.array([1/np.sqrt(2), 1/np.sqrt(2)])
x1=np.array([1,0])
x2=np.array([0,1])

state_lst=(x0,x1,x2)

pauli_str_lst = ([(1.0, 'X0')],[(1.0, 'Y0')],[(1.0, 'Z0')])

for gate in pauli_str_lst:    
#     print (gate)
    H = Hamiltonian(gate)
    for state in state_lst:
        psi=paddle_quantum.State(paddle.to_tensor(state,dtype='complex64'), dtype='complex64')
        print("Gate {}, State {}, Exp-val {}".format(gate, psi, psi.expec_val(H).numpy()))

Gate [(1.0, 'X0')], State [0.70710677+0.j 0.70710677+0.j], Exp-val [0.99999994]
Gate [(1.0, 'X0')], State [1.+0.j 0.+0.j], Exp-val [0.]
Gate [(1.0, 'X0')], State [0.+0.j 1.+0.j], Exp-val [0.]
Gate [(1.0, 'Y0')], State [0.70710677+0.j 0.70710677+0.j], Exp-val [0.]
Gate [(1.0, 'Y0')], State [1.+0.j 0.+0.j], Exp-val [0.]
Gate [(1.0, 'Y0')], State [0.+0.j 1.+0.j], Exp-val [0.]
Gate [(1.0, 'Z0')], State [0.70710677+0.j 0.70710677+0.j], Exp-val [0.]
Gate [(1.0, 'Z0')], State [1.+0.j 0.+0.j], Exp-val [1.]
Gate [(1.0, 'Z0')], State [0.+0.j 1.+0.j], Exp-val [-1.]

Thanks

[LeiZhang-116-4]

Some codes in terms of Hamiltonian and State could be further simplified, as shown below

import numpy as np
import paddle_quantum as pq
from paddle_quantum import Hamiltonian
from paddle_quantum.state import State

pq.set_backend('state_vector')
pq.set_dtype('complex128')

x0 = np.array([1, 1]) / np.sqrt(2)
x1 = np.array([1, 0])
x2 = np.array([0, 1])
state_list = (x0,x1,x2)
pauli_str_list = ([(1.0, 'X0')],[(1.0, 'Y0')],[(1.0, 'Z0')])

for pauli_str in pauli_str_list:
    H = Hamiltonian(pauli_str)
    for vec in state_list:
        psi = State(vec)
        print(f"Observable {H.pauli_words}, State {psi.numpy()}, Exp val {psi.expec_val(H)}")

You can refer to the API documents of Hamiltonian and State for more information.

[BoltzmannEntropy]

Here,

from sympy import Matrix, symbols, sqrt, init_printing
from sympy.physics.quantum import TensorProduct
from IPython.display import display_pretty

init_printing(use_latex=True)


U_I = Matrix([[1,0],
              [0,1]])
U_H = 1/sqrt(2)*Matrix([[1, 1],
                        [1,-1]])

U_Z=Matrix(2,2,[1,0,0,-1])
# U_Z=Matrix(2,2,[1,0,0,-1])

U_ZZ = TensorProduct(U_Z,U_Z)
phi_plus = Matrix(1,4,np.array([1, 0, 0, 1])/np.sqrt(2)).T

# phi_plus*(U_ZZ*phi_plus)
U_ZZ, phi_plus, U_ZZ*phi_plus, U_Z, U_ZZ,  phi_plus.T*(U_ZZ*phi_plus)

I get expectation value of 1 from phi_plus.T*(U_ZZ*phi_plus).

But in paddle:

# matmul(matmul(U, self.rho), dagger(U))
# paddle.kron(proj, X) 
U_Z=paddle.to_tensor([[1,0],[0,-1]], dtype='complex128')
phi_plus = paddle.Tensor(np.array([1, 0, 0, 1])/np.sqrt(2))# | Phi^+ >
phi_plus=paddle.unsqueeze(phi_plus, axis=1)
# U_ZZ = paddle.kron(U_Z,pq.linalg.dagger(U_Z)).numpy()
U_ZZ = paddle.kron(U_Z,U_Z)
exp_val=(phi_plus.T*(U_ZZ*phi_plus)).numpy()
U_ZZ, phi_plus, U_ZZ*phi_plus, exp_val 

I get a matrix instead of a scaler, what am I doing wrong?

In numpy:

U_Z=[[1,0],[0,-1]]
phi_plus = (np.matrix([1, 0, 0, 1])/np.sqrt(2)).T# | Phi^+ >
phi_plus.conj().T*(np.kron(U_Z, U_Z)*phi_plus)
matrix([[1.]])

[BoltzmannEntropy]

?

[LeiZhang-116-4]

The operation * for matrix-like classes (such as sympy.Matrix, numpy.matrix) means matrix multiplication, whereas the operation * for array-like classes (such as numpy.ndarray, paddle.Tensor) means entry-wise multiplication. Please use the symbol @ or the function paddle.matmul instead.

[imppresser]

One naive example for the expectance calculation could be

from paddle_quantum import Hamiltonian
from paddle_quantum.state import random_state

psi = random_state(2)
pauli_str = [(1.0, 'Z0, Z1'), (2.0, 'X0, Y1')]
H = Hamiltonian(pauli_str)
print("The expectation value of the observable ZZ + 2XY is", psi.expec_val(H))

Above example cannot be used in a QNN, in which case you should try

from paddle_quantum.loss import ExpecVal

loss_fcn = ExpecVal(H)
print("The loss value is", loss_fcn(psi))

Hey, thank you, but what if I need to compute <ψ|H^2|ψ>, how to use the psi.expec_val(H)?

[BoltzmannEntropy]

@imppresser see the example here: https://faculty.washington.edu/seattle/physics541/11solved.pdf

[LeiZhang-116-4]

Hi there. The function State.expec_val() only accepts inputs from class Hamiltonian, to be compatible with the situation when the backend is switched from simulators to real quantum devices.

To calculate $\langle \psi | H^2 |\psi\rangle$, you can use the bra-ket notation and matrix multiplications to complete the job.

from paddle_quantum import Hamiltonian
from paddle_quantum.state import random_state
from paddle_quantum.qinfo import random_pauli_str_generator

num_qubits = 3
psi = random_state(num_qubits)
H = Hamiltonian(random_pauli_str_generator(num_qubits))
H_matrix = paddle.to_tensor(H.construct_h_matrix())

expect_value = psi.bra @ H_matrix @ H_matrix @ psi.ket
print("The pauli string of this Hamiltonian is \n", H.pauli_str)
print("The expectation value is", expect_value.item())

The pauli string of this Hamiltonian is 
 [[0.6976072696164013, 'Z1,Z0'], [0.6145107870116349, 'Z1,X0'], [-0.5089016842519825, 'Y1,Y2,Y0']]
The expectation value is (1.177162766456604-2.9802322387695312e-08j)

[imppresser]

@imppresser see the example here: https://faculty.washington.edu/seattle/physics541/11solved.pdf

Thank you for your detailed lecture note, caculating the H^2 is surely the way. However, in my case need more general way of implementing code of H^2, because the variety of H. Very happy to find the course, it will help a lot!

[imppresser]

Hi there. The function State.expec_val() only accepts inputs from class Hamiltonian, to be compatible with the situation when the backend is switched from simulators to real quantum devices.

To calculate ⟨ψ|H2|ψ⟩, you can use the bra-ket notation and matrix multiplications to complete the job.

from paddle_quantum import Hamiltonian
from paddle_quantum.state import random_state
from paddle_quantum.qinfo import random_pauli_str_generator

num_qubits = 3
psi = random_state(num_qubits)
H = Hamiltonian(random_pauli_str_generator(num_qubits))
H_matrix = paddle.to_tensor(H.construct_h_matrix())

expect_value = psi.bra @ H_matrix @ H_matrix @ psi.ket
print("The pauli string of this Hamiltonian is \n", H.pauli_str)
print("The expectation value is", expect_value.item())

The pauli string of this Hamiltonian is 
 [[0.6976072696164013, 'Z1,Z0'], [0.6145107870116349, 'Z1,X0'], [-0.5089016842519825, 'Y1,Y2,Y0']]
The expectation value is (1.177162766456604-2.9802322387695312e-08j)

Thanks for the detailed feedback. I noticed the "expect_val()" function changed (for compatibility). The solution you give will help me out!

[imppresser]

When N(as qubits) turns big like 18, or above, it shows errors for it cannot allocate this big memory. But using "ExpecVal(H)" would not cause error even as large as 20+. If it is possible to run big tasks while expect_value(H^2)?

[BoltzmannEntropy]

Are you using a GPU? Can I see the a full code snippet?

[imppresser]

Are you using a GPU? Can I see the a full code snippet?

Not using GPU. The problem is to construct h_matrix, it is natural cannot build 2^18 x 2^18 complex64 matrix on PC, it will need 500+GB.

[LeiZhang-116-4]

You have reached another motivation to use the class Hamiltonian. To be specific, the matrices of physical Hamiltonians are usually enriched with zero entries, i.e. the sparse matrix. The linear operations for such matrices can be sped up by tensor contractions (see paddle.einsum used in Paddle Quantum), saving a significant amout of computational resources.

Regarding your question, multiplication for the Hamiltonian class is not yet implemented. But, it has been included for future discussions. For now one naive and possibly expensive approach is to extract the Pauli strings of $H$, calculate the Pauli strings of $H^2$ using the identity of Pauli basis $$\sigma_j \ \sigma_k = \delta_{jk} \ I + i \epsilon_{jkl} \ \sigma_{l},$$ and then feed the result back into ExpecVal or State.expec_val.

[imppresser]

You have reached another motivation to use the class Hamiltonian. To be specific, the matrices of physical Hamiltonians are usually enriched with zero entries, i.e. the sparse matrix. The linear operations for such matrices can be sped up by tensor contractions (see paddle.einsum used in Paddle Quantum), saving a significant amout of computational resources.

Regarding your question, multiplication for the Hamiltonian class is not yet implemented. But, it has been included for future discussions. For now one naive and possibly expensive approach is to extract the Pauli strings of H, calculate the Pauli strings of H2 using the identity of Pauli basis σj σk=δjk I+iϵjkl σl, and then feed the result back into ExpecVal or State.expec_val.

I will try this, thank you. Here I have two questions.

1. I use HVA Ansatz(VQE) to get the ground state energy in 1d Transverse-field Ising model, and can get energy below the real ground state energy, when N=4 with periodic boundry condition and certain transverse strength g=1, the ground state energy is -5.22625186 (Exact Diagonalization). The Loss gives me -5.226253032684326 which belows the real ground state energy. If that is possible, why?
2. I use "state.bra @ H_matrix @ state.ket" to caculate the $\langle \psi | H |\psi\rangle$, the result is slightly different from the "ExpecVal or State.expec_val". How could that be? regE: -5.226253509521484 loss: -5.226253032684326 loss_expec_val: -5.226253032684326 regE is the result of state.bra @ H_matrix @ state.ket; loss is the result of ExpecVal(Hamiltonian); loss_expec_val is the result of State.expec_val(Hamiltonian) Many thanks to your support!

[imppresser]

You have reached another motivation to use the class Hamiltonian. To be specific, the matrices of physical Hamiltonians are usually enriched with zero entries, i.e. the sparse matrix. The linear operations for such matrices can be sped up by tensor contractions (see paddle.einsum used in Paddle Quantum), saving a significant amout of computational resources.

Regarding your question, multiplication for the Hamiltonian class is not yet implemented. But, it has been included for future discussions. For now one naive and possibly expensive approach is to extract the Pauli strings of H, calculate the Pauli strings of H2 using the identity of Pauli basis σj σk=δjk I+iϵjkl σl, and then feed the result back into ExpecVal or State.expec_val.

I tried to construct pauli string using σj σk=δjk I+iϵjkl σl, one simple example like Z0X0=iY0, the coefficient will need to multiply i. An original pauli string [[1.0,'Z0, X0']] (cannot use directly, will cause error AssertionError: each Pauli operator should act on different qubit), converted to [[1.0j,'Y0']], It will cause error, because it only accept coefficient as float type.

[LeiZhang-116-4]

1. I use HVA Ansatz(VQE) to get the ground state energy in 1d Transverse-field Ising model, and can get energy below the real ground state energy, when N=4 with periodic boundry condition and certain transverse strength g=1, the ground state energy is -5.22625186 (Exact Diagonalization). The Loss gives me -5.226253032684326 which belows the real ground state energy. If that is possible, why?
2. I use "state.bra @ H_matrix @ state.ket" to caculate the ⟨ψ|H|ψ⟩, the result is slightly different from the "ExpecVal or State.expec_val". How could that be?

The variations among State.expec_val, state.bra @ H_matrix @ state.ket and the theorectical value are due to the inevitable precision problem, which increases exponentially as the number of qubits gets larger (and hence motivates the development of quantum computers). Particularly, the precision differences between State.expec_val and state.bra @ H_matrix @ state.ket are essentially the difference between tensor contractions and matrix multiplications, where tensor contraction is slightly better since it deals with smaller matrices.

To relieve the precision problem when the qubit size is not too large, you can switch the data type (dtype) from complex64 to complex128, improving the precision from $\mathcal{O}\left( \ 10^{-8}\ \right)$ to $\mathcal{O}\left( \ 10^{-16}\ \right)$ in average, as demonstrated in the following example:

import paddle_quantum as pq
from paddle_quantum import Hamiltonian
from paddle_quantum.state import ghz_state
from paddle_quantum.qinfo import random_pauli_str_generator

num_qubits = 4
pauli_str = random_pauli_str_generator(num_qubits, terms=10)

# Compute error from imaginary part of expectation value
def expect_error() -> float:
    psi = ghz_state(num_qubits)
    H = Hamiltonian(pauli_str)
    H_matrix = paddle.to_tensor(H.construct_h_matrix())
    value = psi.bra @ H_matrix @ H_matrix @ psi.ket
    return abs(value.imag().item())

# Calculate the expectance value under complex64
pq.set_dtype('complex64')
error_64 = expect_error()

# Calculate the expectance value under complex128
pq.set_dtype('complex128')
error_128 = expect_error()

print("The error for complex64 is", error_64)
print("The error for complex128 is", error_128)

The error for complex64 is 1.1195256277574117e-08
The error for complex128 is 9.813077866773593e-18

Note that the function paddle_quantum.set_dype should be excecuted at the very beginning of the program, so that every operation in Paddle Quantum can be computed under dtype complex128. Also, as a trade-off, computations in complex128 are slower than those in complex64.

[LeiZhang-116-4]

I tried to construct pauli string using σj σk=δjk I+iϵjkl σl, one simple example like Z0X0=iY0, the coefficient will need to multiply i. An original pauli string [[1.0,'Z0, X0']] (cannot use directly, will cause error AssertionError: each Pauli operator should act on different qubit), converted to [[1.0j,'Y0']], It will cause error, because it only accept coefficient as float type.

The reason why the class Hamiltonian only accepts real coefficients, is that Hamiltonians are Hermitian matrices and thus have only real eigenvalues. For example ZX cannot be a Hamiltonian since its eigenvalues are i and -i.

For the same reason, I believe every imaginary coefficient occured in one term of $H^2$ can be vanished somewhere for meeting its complex conjugates. For example, consider $H = X + Z$. Then

$$H^2 = X^2 + XZ + ZX + Z^2 = I - iY + iY + I = 2I.$$

Such logic may need to be considered during the construction of Pauli strings of $H^2$.

[imppresser]

I tried to construct pauli string using σj σk=δjk I+iϵjkl σl, one simple example like Z0X0=iY0, the coefficient will need to multiply i. An original pauli string [[1.0,'Z0, X0']] (cannot use directly, will cause error AssertionError: each Pauli operator should act on different qubit), converted to [[1.0j,'Y0']], It will cause error, because it only accept coefficient as float type.

The reason why the class Hamiltonian only accepts real coefficients, is that Hamiltonians are Hermitian matrices and thus have only real eigenvalues. For example ZX cannot be a Hamiltonian since its eigenvalues are i and -i.

For the same reason, I believe every imaginary coefficient occured in one term of H2 can be vanished somewhere for meeting its complex conjugates. For example, consider H=X+Z. Then

H2=X2+XZ+ZX+Z2=I−iY+iY+I=2I.

Such logic may need to be considered during the construction of Pauli strings of H2.

Yes, indeed. I retain every term when I trying to build the H^2 pauli string function. And I should care the imaginary coefficient can be vanished, thanks!

[imppresser]

1. I use HVA Ansatz(VQE) to get the ground state energy in 1d Transverse-field Ising model, and can get energy below the real ground state energy, when N=4 with periodic boundry condition and certain transverse strength g=1, the ground state energy is -5.22625186 (Exact Diagonalization). The Loss gives me -5.226253032684326 which belows the real ground state energy. If that is possible, why?
2. I use "state.bra @ H_matrix @ state.ket" to caculate the ⟨ψ|H|ψ⟩, the result is slightly different from the "ExpecVal or State.expec_val". How could that be?

The variations among State.expec_val, state.bra @ H_matrix @ state.ket and the theorectical value are due to the inevitable precision problem, which increases exponentially as the number of qubits gets larger (and hence motivates the development of quantum computers). Particularly, the precision differences between State.expec_val and state.bra @ H_matrix @ state.ket are essentially the difference between tensor contractions and matrix multiplications, where tensor contraction is slightly better since it deals with smaller matrices.

To relieve the precision problem when the qubit size is not too large, you can switch the data type (dtype) from complex64 to complex128, improving the precision from O( 10−8 ) to O( 10−16 ) in average, as demonstrated in the following example:

import paddle_quantum as pq
from paddle_quantum import Hamiltonian
from paddle_quantum.state import ghz_state
from paddle_quantum.qinfo import random_pauli_str_generator

num_qubits = 4
pauli_str = random_pauli_str_generator(num_qubits, terms=10)

# Compute error from imaginary part of expectation value
def expect_error() -> float:
    psi = ghz_state(num_qubits)
    H = Hamiltonian(pauli_str)
    H_matrix = paddle.to_tensor(H.construct_h_matrix())
    value = psi.bra @ H_matrix @ H_matrix @ psi.ket
    return abs(value.imag().item())

# Calculate the expectance value under complex64
pq.set_dtype('complex64')
error_64 = expect_error()

# Calculate the expectance value under complex128
pq.set_dtype('complex128')
error_128 = expect_error()

print("The error for complex64 is", error_64)
print("The error for complex128 is", error_128)

The error for complex64 is 1.1195256277574117e-08
The error for complex128 is 9.813077866773593e-18

Note that the function paddle_quantum.set_dype should be excecuted at the very beginning of the program, so that every operation in Paddle Quantum can be computed under dtype complex128. Also, as a trade-off, computations in complex128 are slower than those in complex64.

Nice solution, I will try, thanks 👍