Simulation

The simulation module contains the functions used by the simulation backend to construct the transverse Ising Hamiltonian \(H\) and compute the Gibbs state density matrix \(\rho = e^{-\beta H} / \mathcal{Z}\) exactly over the full \(2^{n_{\text{qubits}}}\)-dimensional Hilbert space. Note that the matrices scale exponentially in the number of qubits, so this backend is only feasible for small models. The Hamiltonian formulation is explained in more detail on the Theory page.

compute_H(h, J, A, B, n_qubits, pauli_kron)

Computes the Hamiltonian of the annealer at relative time s.

Parameters:

Name	Type	Description	Default
h	ndarray	Linear Ising terms.	required
J	ndarray	Quadratic Ising terms.	required
A	float	Coefficient of the off-diagonal terms, e.g. A(s).	required
B	float	Coefficient of the diagonal terms, e.g. B(s).	required
n_qubits	int	Number of qubits.	required
pauli_kron	PauliKron	Kronecker product Pauli matrices dict.	required

Returns:

Type	Description
ndarray	Hamiltonian matrix H.

Raises:

Type	Description
ValueError	If the length of h or the shape of J does not match n_qubits.

Source code in src/qbm/simulation/simulation.py

def compute_H(
    h: np.ndarray,
    J: np.ndarray,
    A: float,
    B: float,
    n_qubits: int,
    pauli_kron: PauliKron,
) -> np.ndarray:
    """
    Computes the Hamiltonian of the annealer at relative time s.

    Args:
        h: Linear Ising terms.
        J: Quadratic Ising terms.
        A: Coefficient of the off-diagonal terms, e.g. A(s).
        B: Coefficient of the diagonal terms, e.g. B(s).
        n_qubits: Number of qubits.
        pauli_kron: Kronecker product Pauli matrices dict.

    Returns:
        Hamiltonian matrix H.

    Raises:
        ValueError: If the length of h or the shape of J does not match n_qubits.
    """
    if len(h) != n_qubits:
        raise ValueError(f"h has length {len(h)}, expected n_qubits = {n_qubits}")
    if J.shape != (n_qubits, n_qubits):
        raise ValueError(f"J has shape {J.shape}, expected ({n_qubits}, {n_qubits})")
    # Diagonal terms
    H_diag = np.zeros(2**n_qubits)
    for i in range(n_qubits):
        # Linear terms
        if h[i] != 0:
            H_diag += (B * h[i]) * pauli_kron["z_diag", i]

        # Quadratic terms
        for j in range(i + 1, n_qubits):
            if J[i, j] != 0:
                H_diag += (B * J[i, j]) * pauli_kron["zz_diag", i, j]

    # Return just the diagonal if H is a diagonal matrix
    if A == 0:
        return np.diag(H_diag)

    # Off-diagonal terms
    H = csr_matrix((2**n_qubits, 2**n_qubits), dtype=np.float64)
    for i in range(n_qubits):
        H -= A * pauli_kron["x", i]

    return (H + diags(H_diag, format="csr")).toarray()

compute_rho(H, beta, diagonal=False)

Computes the trace normalized density matrix rho.

Parameters:

Name	Type	Description	Default
H	ndarray	Hamiltonian matrix.	required
beta	float	Inverse temperature beta = 1 / (k_B * T).	required
diagonal	bool	Flag to indicate whether H is a diagonal matrix or not.	False

Returns:

Type	Description
ndarray	Density matrix rho.

Raises:

Type	Description
ValueError	If H is not a square matrix or if beta is not positive.

Source code in src/qbm/simulation/simulation.py

def compute_rho(H: np.ndarray, beta: float, diagonal: bool = False) -> np.ndarray:
    """
    Computes the trace normalized density matrix rho.

    Args:
        H: Hamiltonian matrix.
        beta: Inverse temperature beta = 1 / (k_B * T).
        diagonal: Flag to indicate whether H is a diagonal matrix or not.

    Returns:
        Density matrix rho.

    Raises:
        ValueError: If H is not a square matrix or if beta is not positive.
    """
    if H.ndim != 2 or H.shape[0] != H.shape[1]:
        raise ValueError(f"H must be a square matrix (got shape {H.shape})")
    if beta <= 0:
        raise ValueError(f"beta must be positive (got {beta})")
    # If diagonal then compute directly, else use eigen decomposition
    if diagonal:
        Lambda = H.diagonal()
        exp_beta_Lambda = np.exp(-beta * (Lambda - Lambda.min()))
        return np.diag(exp_beta_Lambda / exp_beta_Lambda.sum())
    else:
        Lambda, S = eigh(H)
        exp_beta_Lambda = np.exp(-beta * (Lambda - Lambda.min()))
        return (S * (exp_beta_Lambda / exp_beta_Lambda.sum())) @ S.T

get_pauli_kron(n_visible, n_hidden)

Computes the necessary Pauli Kronecker product (sparse) matrices for a n_visible + n_hidden qubit problem. Used as an argument to compute_H, e.g. one would instantiate pauli_kron as pauli_kron = get_pauli_kron(n_visible, n_hidden), then pass to compute_H when computing the Hamiltonian.

Parameters:

Name	Type	Description	Default
n_visible	int	Number of visible units.	required
n_hidden	int	Number of hidden units.	required

Returns:

Type	Description
PauliKron	Dictionary of Kronecker product Pauli terms, with keys ("x", i) mapping to
PauliKron	sparse matrices I ⊗ σ_x^(i) ⊗ I, and keys ("z_diag", i) and ("zz_diag", i, j)
PauliKron	mapping to the diagonals of I ⊗ σ_z^(i) ⊗ I and their pairwise products.

Raises:

Type	Description
ValueError	If n_visible or n_hidden is not positive.

Source code in src/qbm/simulation/simulation.py

def get_pauli_kron(n_visible: int, n_hidden: int) -> PauliKron:
    """
    Computes the necessary Pauli Kronecker product (sparse) matrices for a n_visible
    + n_hidden qubit problem. Used as an argument to compute_H, e.g. one would
    instantiate pauli_kron as pauli_kron = get_pauli_kron(n_visible, n_hidden), then
    pass to compute_H when computing the Hamiltonian.

    Args:
        n_visible: Number of visible units.
        n_hidden: Number of hidden units.

    Returns:
        Dictionary of Kronecker product Pauli terms, with keys ("x", i) mapping to
        sparse matrices I ⊗ σ_x^(i) ⊗ I, and keys ("z_diag", i) and ("zz_diag", i, j)
        mapping to the diagonals of I ⊗ σ_z^(i) ⊗ I and their pairwise products.

    Raises:
        ValueError: If n_visible or n_hidden is not positive.
    """
    if n_visible <= 0:
        raise ValueError(f"n_visible must be positive (got {n_visible})")
    if n_hidden <= 0:
        raise ValueError(f"n_hidden must be positive (got {n_hidden})")
    # Set Kronecker product Pauli matrices
    n_qubits = n_visible + n_hidden
    pauli_kron = {}
    for i in range(n_qubits):
        pauli_kron["x", i] = sparse_kron(i, n_qubits, sparse_X)
        pauli_kron["z_diag", i] = sparse_kron(i, n_qubits, sparse_Z).diagonal()
    for i in range(n_qubits):
        for j in range(i + 1, n_qubits):
            pauli_kron["zz_diag", i, j] = (
                pauli_kron["z_diag", i] * pauli_kron["z_diag", j]
            )

    return pauli_kron

sparse_kron(i, n_qubits, A)

Compute I_{2^i} ⊗ A ⊗ I_{2^(n_qubits-i-1)}.

Parameters:

Name	Type	Description	Default
i	int	Index of the "A" matrix.	required
n_qubits	int	Total number of qubits.	required
A	spmatrix	Matrix to tensor with identities.	required

Returns:

Type	Description
Any	I_{2^i} ⊗ A ⊗ I_{2^(n_qubits-i-1)}.

Source code in src/qbm/simulation/simulation.py

def sparse_kron(i: int, n_qubits: int, A: spmatrix) -> Any:
    """
    Compute I_{2^i} ⊗ A ⊗ I_{2^(n_qubits-i-1)}.

    Args:
        i: Index of the "A" matrix.
        n_qubits: Total number of qubits.
        A: Matrix to tensor with identities.

    Returns:
        I_{2^i} ⊗ A ⊗ I_{2^(n_qubits-i-1)}.
    """
    if i != 0 and i != n_qubits - 1:
        return kron(kron(identity(2**i), A), identity(2 ** (n_qubits - i - 1)))
    if i == 0:
        return kron(A, identity(2 ** (n_qubits - 1)))
    if i == n_qubits - 1:
        return kron(identity(2 ** (n_qubits - 1)), A)