Convert from qubit operation to axis-angle Bloch sphere rotation

Question

Given the 2x2 unitary matrix representation of an operation to apply to a single qubit, how do I figure out the rotation it corresponds to on the Bloch sphere?

For example, the Hadamard matrix is a 180 degree rotation around the X+Z axis. How do I get from [[1,1],[1,-1]]*sqrt(0.5) to (X+Z, 180 deg)?

Answer 1

Single-qubit operations are basically just unit quaternions, but with an extra phase factor. The similarity is because the Pauli matrices, times sqrt(-1), satisfy the i^2=j^2=k^2=ijk=-1 relation that defines quaternions.

As a result, the hard part of the conversion method is already taken care of by any "quaternion to axis angle" code. Just pull out the phased quaternion components, figure out the phase factor, then apply the quaternion-to-angle-axis method.

import math
import cmath

def toBlochAngleAxis(matrix):
    """
    Breaksdown a matrix U into axis, angle, and phase_angle components satisfying
    U = exp(i phase_angle) (I cos(angle/2) - axis sigma i sin(angle/2))

    :param matrix: The 2x2 unitary matrix U
    :return: The breakdown (axis(x, y, z), angle, phase_angle)
    """
    [[a, b], [c, d]] = matrix

    # --- Part 1: convert to a quaternion ---

    # Phased components of quaternion.
    wp = (a + d) / 2.0
    xp = (b + c) / 2.0j
    yp = (b - c) / 2.0
    zp = (a - d) / 2.0j

    # Arbitrarily use largest value to determine the global phase factor.
    phase = max([wp, xp, yp, zp], key=abs)
    phase /= abs(phase)

    # Cancel global phase factor, recovering quaternion components.
    w = complex(wp / phase).real
    x = complex(xp / phase).real
    y = complex(yp / phase).real
    z = complex(zp / phase).real

    # --- Part 2: convert from quaternion to angle-axis ---

    # Floating point error may have pushed w outside of [-1, +1]. Fix that.
    w = min(max(w, -1), +1)

    # Recover angle.
    angle = -2*math.acos(w)

    # Normalize axis.
    n = math.sqrt(x*x + y*y + z*z);
    if n < 0.000001:
        # There's an axis singularity near angle=0.
        # Just default to no rotation around the Z axis in this case.
        angle = 0
        x = 0
        y = 0
        z = 1
        n = 1
    x /= n
    y /= n
    z /= n

    # --- Part 3: (optional) canonicalize ---

    # Prefer angle in [-pi, pi]
    if angle <= -math.pi:
        angle += 2*math.pi
        phase *= -1

    # Prefer axes that point positive-ward.
    if x + y + z < 0:
        x *= -1
        y *= -1
        z *= -1
        angle *= -1

    phase_angle = cmath.polar(phase)[1]

    return (x, y, z), angle, phase_angle

Testing it out:

print(toBlochAngleAxis([[1, 0], [0, 1]])) # Identity
# ([0, 0, 1], 0, 0.0)

print(toBlochAngleAxis([[0, 1], [1, 0]])) # Pauli X, 180 deg around X
# ([1.0, -0.0, -0.0], 3.141592653589793, 1.5707963267948966)

print(toBlochAngleAxis([[0, -1j], [1j, 0]])) # Pauli Y, 180 deg around Y
# ([-0.0, 1.0, -0.0], 3.141592653589793, 1.5707963267948966)

print(toBlochAngleAxis([[1, 0], [0, -1]])) # Pauli Z, 180 deg around Z
# ([-0.0, -0.0, 1.0], 3.141592653589793, 1.5707963267948966)

s = math.sqrt(0.5)
print(toBlochAngleAxis([[s, s], [s, -s]])) # Hadamard, 180 deg around X+Z
# ([0.7071067811865476, -0.0, 0.7071067811865476], 3.141592653589793, 1.5707963267948966)

print(toBlochAngleAxis([[s, s*1j], [s*1j, s]])) # -90 deg X axis, no phase
# ((1.0, 0.0, 0.0), -1.5707963267948966, 0.0)