circuit.angles() seems incorrect

[wangzbo]

It is easy to recreate. Here is the circuit description using openQASM:

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
creg c[2];
h q[0];
h q[1];
s q[0];
t q[1];
cx q[0], q[1];
sx q[0];

I call circuit.angles() and it returns "theta: 0.785398163399, phi: 0, thetaDeg: 45, phiDeg: 0, radius: 0.7071068" for q[0]. The partial trace of q[0] is diag((2+√2)/4, (2-√2)/4), so the correct angle of q[0] should be "theta: 0, phi: 0, thetaDeg: 0, phiDeg: 0, radius: 0.7071068"?

[perak]

Hi @wangzbo thank you for reporting this.

I am trying to find what is wrong, but I am getting results which looks good to me.

Please see this:

import numpy as np


def partial_trace(rho, dims, traced_out):
    """
    Compute the partial trace of a density matrix.
    :param rho: Density matrix (numpy array)
    :param dims: List of subsystem dimensions (e.g., [2, 2] for 2 qubits)
    :param traced_out: Index of the qubit to trace out (0 or 1)
    :return: Reduced density matrix after tracing out the specified qubit
    """
    dimA, dimB = dims
    reshaped_rho = rho.reshape(dimA, dimB, dimA, dimB)
    
    if traced_out == 0:
        reduced_rho = np.trace(reshaped_rho, axis1=0, axis2=2)
    elif traced_out == 1:
        reduced_rho = np.trace(reshaped_rho, axis1=1, axis2=3)
    else:
        raise ValueError("Invalid qubit index to trace out")
    
    return reduced_rho


def bloch_sphere_angles(rho):
    """
    Compute Bloch sphere angles (theta, phi) and radius from a single-qubit density matrix.
    :param rho: 2x2 density matrix
    :return: (theta, phi, radius) for Bloch sphere representation
    """
    pauli_matrices = [
        np.array([[0, 1], [1, 0]]),  # sigma_x
        np.array([[0, -1j], [1j, 0]]),  # sigma_y
        np.array([[1, 0], [0, -1]])   # sigma_z
    ]
    
    # Compute the Bloch vector components
    bloch_vector = np.array([np.trace(rho @ sigma).real for sigma in pauli_matrices])
    
    # Calculate theta and phi
    theta = np.arccos(bloch_vector[2])  # cos(theta) = z-component
    phi = np.arctan2(bloch_vector[1], bloch_vector[0])  # atan2(y, x)
    
    # Calculate the radius (magnitude of the Bloch vector)
    radius = np.linalg.norm(bloch_vector)
    
    return theta, phi, radius

# Define the state vector
state_vector = np.array([ 0.25+0.60355339059327j, -0.10355339059327-0.25j, 0.25+0.60355339059327j, 0.10355339059327+0.25j ])

# Calculate the density matrix
density_matrix = np.outer(state_vector, np.conj(state_vector))

reduced_density_matrix_0 = partial_trace(density_matrix, [2, 2], traced_out=0)
reduced_density_matrix_1 = partial_trace(density_matrix, [2, 2], traced_out=1)


# Print the results
print("Density Matrix:")
print(repr(density_matrix))

print("\n*** Qubit 0 ***")
print("\nReduced Density Matrix:\n")
print(repr(reduced_density_matrix_0))
theta_0, phi_0, radius_0 = bloch_sphere_angles(reduced_density_matrix_0)
print("\nBloch sphere angles:\n")
print(f"Theta: {theta_0}, Phi: {phi_0}, Radius: {radius_0}")

print("\n*** Qubit 1 ***")
print("\nReduced Density Matrix:\n")
print(repr(reduced_density_matrix_1))
theta_1, phi_1, radius_1 = bloch_sphere_angles(reduced_density_matrix_1)
print("\nBloch sphere angles:\n")
print(f"Theta: {theta_1}, Phi: {phi_1}, Radius: {radius_1}")

Prints:

Density Matrix:
array([[ 0.4267767+0.00000000e+00j, -0.1767767+2.66453526e-15j,
         0.4267767+0.00000000e+00j,  0.1767767-2.66453526e-15j],
       [-0.1767767-2.66453526e-15j,  0.0732233+0.00000000e+00j,
        -0.1767767-2.66453526e-15j, -0.0732233+0.00000000e+00j],
       [ 0.4267767+0.00000000e+00j, -0.1767767+2.66453526e-15j,
         0.4267767+0.00000000e+00j,  0.1767767-2.66453526e-15j],
       [ 0.1767767+2.66453526e-15j, -0.0732233+0.00000000e+00j,
         0.1767767+2.66453526e-15j,  0.0732233+0.00000000e+00j]])

*** Qubit 0 ***

Reduced Density Matrix:

array([[0.85355339+0.j, 0.        +0.j],
       [0.        +0.j, 0.14644661+0.j]])

Bloch sphere angles:

Theta: 0.785398163397459, Phi: 0.0, Radius: 0.7071067811865399

*** Qubit 1 ***

Reduced Density Matrix:

array([[0.5       +0.j, 0.35355339+0.j],
       [0.35355339+0.j, 0.5       +0.j]])

Bloch sphere angles:

Theta: 1.5707963267948966, Phi: 0.0, Radius: 0.7071067811865399

OR with Qiskit:

import numpy as np

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.quantum_info import Statevector, DensityMatrix, partial_trace, Pauli

qc = QuantumCircuit()

q = QuantumRegister(2, 'q')
c = ClassicalRegister(2, 'c')

qc.add_register(q)
qc.add_register(c)

qc.h(q[0])
qc.h(q[1])
qc.s(q[0])
qc.t(q[1])
qc.cx(q[0], q[1])
qc.sx(q[0])

# Get the statevector of the circuit
state = Statevector.from_instruction(qc)

# Convert to density matrix
density_matrix = DensityMatrix(state)

reduced_density_matrix_0 = np.round(partial_trace(density_matrix, [1]).data, 15)
reduced_density_matrix_1 = np.round(partial_trace(density_matrix, [0]).data, 15)


r_x_0 = np.real(reduced_density_matrix_0[0, 1] + reduced_density_matrix_0[1, 0])
r_y_0 = np.imag(reduced_density_matrix_0[0, 1] - reduced_density_matrix_0[1, 0])
r_z_0 = np.real(reduced_density_matrix_0[0, 0] - reduced_density_matrix_0[1, 1])

# Calculate the radius, theta, and phi
radius_0 = np.sqrt(r_x_0**2 + r_y_0**2 + r_z_0**2)
theta_0 = np.arccos(r_z_0)
phi_0 = np.arctan2(r_y_0, r_x_0)

r_x_1 = np.real(reduced_density_matrix_1[0, 1] + reduced_density_matrix_1[1, 0])
r_y_1 = np.imag(reduced_density_matrix_1[0, 1] - reduced_density_matrix_1[1, 0])
r_z_1 = np.real(reduced_density_matrix_1[0, 0] - reduced_density_matrix_1[1, 1])

# Calculate the radius, theta, and phi
radius_1 = np.sqrt(r_x_1**2 + r_y_1**2 + r_z_1**2)
theta_1 = np.arccos(r_z_1)
phi_1 = np.arctan2(r_y_1, r_x_1)


# Print the results
print("Density Matrix:")
print(repr(density_matrix))

print("\n*** Qubit 0 ***")
print("\nReduced Density Matrix:\n")
print(repr(reduced_density_matrix_0))
print("\nBloch sphere angles:\n")
print(f"Theta: {theta_0}, Phi: {phi_0}, Radius: {radius_0}")

print("\n*** Qubit 1 ***")
print("\nReduced Density Matrix:\n")
print(repr(reduced_density_matrix_1))
print("\nBloch sphere angles:\n")
print(f"Theta: {theta_1}, Phi: {phi_1}, Radius: {radius_1}")

Prints:

Density Matrix:
DensityMatrix([[ 0.4267767+5.74836793e-18j, -0.1767767+4.62223187e-33j,
                 0.4267767+5.74836793e-18j,  0.1767767-6.93889390e-18j],
               [-0.1767767+1.03494385e-18j,  0.0732233+1.68365798e-18j,
                -0.1767767+1.03494385e-18j, -0.0732233+1.19052598e-18j],
               [ 0.4267767+5.74836793e-18j, -0.1767767+4.62223187e-33j,
                 0.4267767+5.74836793e-18j,  0.1767767-6.93889390e-18j],
               [ 0.1767767+5.90395006e-18j, -0.0732233-1.68365798e-18j,
                 0.1767767+5.90395006e-18j,  0.0732233-1.19052598e-18j]],
              dims=(2, 2))

*** Qubit 0 ***

Reduced Density Matrix:

array([[0.85355339+0.j, 0.        -0.j],
       [0.        +0.j, 0.14644661+0.j]])

Bloch sphere angles:

Theta: 0.7853981633974491, Phi: -0.0, Radius: 0.707106781186547

*** Qubit 1 ***

Reduced Density Matrix:

array([[0.5       +0.j, 0.35355339+0.j],
       [0.35355339+0.j, 0.5       +0.j]])

Bloch sphere angles:

Theta: 1.5707963267948966, Phi: 0.0, Radius: 0.707106781186548

Now... it is possible that I am doing something wrong. Please give me more details on how you calculate density matrix, partial trace and angles.

[wangzbo]

Hi @perak , thank you for your reply. I think all look correct but the step to compute the r, theta and phi.
The partial trace returned from quantum-circuit api partial_trace() is correct. That is how I calculate the partial trace.
the coordinates x, y z are also correct returned with(0, 0, √2/2). This is a mixed state not a pure state(r < 1), so when caculating the r, theta and phi from x, y and z, it should be something like(from wiki):

This is why we get different results for (r, theta, phi).