Skip to content

Quantum Information API Reference

This page documents the classes and functions available in the skq.quantum_info module. This module provides tools for working with quantum states, density matrices, quantum channels, Hamiltonians, and other quantum information concepts.

State Representations

Statevector

The Statevector class provides a representation for pure quantum states.

from skq.quantum_info import Statevector

# Create a state vector
state = Statevector([1/np.sqrt(2), 1/np.sqrt(2)])
print(state.num_qubits)  # 1
print(state.is_normalized())  # True

skq.quantum_info.Statevector

Bases: ndarray

Statevector representation for a quantum (qubit) state. NOTE: Input is assumed to be in big-endian format. Little-endian -> Least significant qubit (LSB) is on the right. Like |q1 q0> where q0 is the LSB. Big-endian -> Least significant qubit (LSB) is on the left. Like |q0 q1> where q0 is the LSB.

Source code in src/skq/quantum_info/state.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
class Statevector(np.ndarray):
    """
    Statevector representation for a quantum (qubit) state.
    NOTE: Input is assumed to be in big-endian format.
    Little-endian -> Least significant qubit (LSB) is on the right. Like |q1 q0> where q0 is the LSB.
    Big-endian -> Least significant qubit (LSB) is on the left. Like |q0 q1> where q0 is the LSB.
    """

    def __new__(cls, input_array):
        arr = np.asarray(input_array, dtype=complex)
        obj = arr.view(cls)
        assert obj.is_1d(), "State vector must be 1D."
        assert obj.is_normalized(), "State vector must be normalized."
        assert obj.is_power_of_two_len(), "State vector length must be a least length 2 and a power of 2."
        return obj

    @property
    def num_qubits(self) -> int:
        """Number of qubits in the state vector."""
        return int(np.log2(len(self)))

    def is_1d(self) -> bool:
        """State vector is 1D."""
        return self.ndim == 1

    def is_normalized(self) -> bool:
        """State vector is normalized: ||ψ|| = 1."""
        return np.isclose(np.linalg.norm(self), 1)

    def is_power_of_two_len(self) -> bool:
        """Check if a number is a power of two"""
        n = len(self)
        return n >= 2 and (n & (n - 1)) == 0

    def is_multi_qubit(self) -> bool:
        """State vector represents a multi-qubit state."""
        return self.num_qubits > 1

    def is_indistinguishable(self, other: "Statevector") -> bool:
        """Two state vectors are indistinguishable if their density matrices are the same."""
        return np.allclose(self.density_matrix(), other.density_matrix())

    def magnitude(self) -> float:
        """Magnitude (or norm) of the state vector. sqrt(<ψ|ψ>)"""
        np.linalg.norm(self)

    def expectation(self, operator: np.array) -> float:
        """
        Expectation value of an observable with respect to the quantum state.
        computes <ψ|O|ψ>
        :param operator: The observable (Hermitian matrix) as a 2D numpy array.
        :return: Expectation value.
        """
        assert len(operator.shape) == 2, "The operator must be a 2D matrix."
        assert np.allclose(operator, operator.T.conj()), "The operator must be Hermitian."
        return float(np.real(self.conjugate_transpose() @ operator @ self))

    def density_matrix(self) -> DensityMatrix:
        """Return the density matrix representation of the state vector."""
        return DensityMatrix(np.outer(self, self.conj()))

    def probabilities(self) -> np.ndarray:
        """Probabilities of all possible states."""
        return np.abs(self) ** 2

    def measure_index(self) -> int:
        """
        Simulate a measurement of the state and get a sampled index.
        :return: Index of the measured state.
        For example:
        |0> -> 0
        |00> -> 0
        |11> -> 3
        """
        return np.random.choice(len(self), p=self.probabilities())

    def measure_bitstring(self) -> str:
        """
        Simulate a measurement of the state vector and get a bitstring representing the sample.
        :return: Bitstring representation of the measured state.
        For example:
        |0> -> "0"
        |00> -> "00"
        |11> -> "11"
        """
        return bin(self.measure_index())[2:].zfill(self.num_qubits)

    def reverse(self) -> "Statevector":
        """Reverse the order of the state vector to account for endianness.
        For example Qiskit uses little endian convention.
        Little-endian -> Least significant qubit (LSB) is on the right. Like |q1 q0> where q0 is the LSB.
        Big-endian -> Least significant qubit (LSB) is on the left. Like |q0 q1> where q0 is the LSB.
        """
        return Statevector(self[::-1])

    def bloch_vector(self) -> np.ndarray:
        """Bloch vector representation of the quantum state."""
        return self.density_matrix().bloch_vector()

    def conjugate_transpose(self) -> np.ndarray:
        """Conjugate transpose (Hermitian adjoint) of the state vector."""
        return self.T.conj()

    def orthonormal_basis(self) -> np.ndarray:
        """
        Orthonormal basis using the Gram-Schmidt process.
        :return: 2D array representing the orthonormal basis.
        """
        return np.array([self / np.linalg.norm(self)]).T

    def schmidt_decomposition(self) -> tuple[np.array, np.array, np.array]:
        """
        Perform Schmidt decomposition on a quantum state.
        :return: Tuple of Schmidt coefficients, Basis A and Basis B
        """
        # Infer dimensions
        N = len(self)
        dim_A = int(np.sqrt(N))
        dim_B = N // dim_A

        # Singular Value Decomposition
        state_matrix = self.reshape(dim_A, dim_B)
        U, S, Vh = svd(state_matrix)

        # Coefficients (S), Basis A (U) and Basis B (Vh^dagger)
        return S, U, Vh.conj().T

    def to_qiskit(self) -> qiskit.quantum_info.Statevector:
        """
        Convert the state vector to a Qiskit QuantumCircuit object.
        Qiskit uses little-endian convention.
        :return: Qiskit StateVector object
        """
        return qiskit.quantum_info.Statevector(self.reverse())

    @staticmethod
    def from_qiskit(statevector: qiskit.quantum_info.Statevector) -> "Statevector":
        """
        Convert a Qiskit StateVector object to a scikit-q StateVector.
        Qiskit uses little-endian convention.
        :param statevector: Qiskit StateVector object
        :return: scikit-q StateVector object
        """
        return Statevector(statevector.data).reverse()

    def to_pyquil(self):
        """
        Convert the state vector to a PyQuil object.
        PyQuil uses little-endian convention.
        :return: PyQuil object
        """
        raise NotImplementedError("Conversion to PyQuil is not implemented.")

    @staticmethod
    def from_pyquil(statevector) -> "Statevector":
        """
        Convert a PyQuil object to a scikit-q StateVector object.
        PyQuil uses little-endian convention.
        :return: scikit-q StateVector object
        """
        raise NotImplementedError("Conversion from PyQuil is not implemented.")

num_qubits property

Number of qubits in the state vector.

bloch_vector()

Bloch vector representation of the quantum state.

Source code in src/skq/quantum_info/state.py 
103
104
105
def bloch_vector(self) -> np.ndarray:
    """Bloch vector representation of the quantum state."""
    return self.density_matrix().bloch_vector()

conjugate_transpose()

Conjugate transpose (Hermitian adjoint) of the state vector.

Source code in src/skq/quantum_info/state.py 
107
108
109
def conjugate_transpose(self) -> np.ndarray:
    """Conjugate transpose (Hermitian adjoint) of the state vector."""
    return self.T.conj()

density_matrix()

Return the density matrix representation of the state vector.

Source code in src/skq/quantum_info/state.py 
65
66
67
def density_matrix(self) -> DensityMatrix:
    """Return the density matrix representation of the state vector."""
    return DensityMatrix(np.outer(self, self.conj()))

expectation(operator)

Expectation value of an observable with respect to the quantum state. computes <ψ|O|ψ> :param operator: The observable (Hermitian matrix) as a 2D numpy array. :return: Expectation value.

Source code in src/skq/quantum_info/state.py 
54
55
56
57
58
59
60
61
62
63
def expectation(self, operator: np.array) -> float:
    """
    Expectation value of an observable with respect to the quantum state.
    computes <ψ|O|ψ>
    :param operator: The observable (Hermitian matrix) as a 2D numpy array.
    :return: Expectation value.
    """
    assert len(operator.shape) == 2, "The operator must be a 2D matrix."
    assert np.allclose(operator, operator.T.conj()), "The operator must be Hermitian."
    return float(np.real(self.conjugate_transpose() @ operator @ self))

from_pyquil(statevector) staticmethod

Convert a PyQuil object to a scikit-q StateVector object. PyQuil uses little-endian convention. :return: scikit-q StateVector object

Source code in src/skq/quantum_info/state.py 
161
162
163
164
165
166
167
168
@staticmethod
def from_pyquil(statevector) -> "Statevector":
    """
    Convert a PyQuil object to a scikit-q StateVector object.
    PyQuil uses little-endian convention.
    :return: scikit-q StateVector object
    """
    raise NotImplementedError("Conversion from PyQuil is not implemented.")

from_qiskit(statevector) staticmethod

Convert a Qiskit StateVector object to a scikit-q StateVector. Qiskit uses little-endian convention. :param statevector: Qiskit StateVector object :return: scikit-q StateVector object

Source code in src/skq/quantum_info/state.py 
143
144
145
146
147
148
149
150
151
@staticmethod
def from_qiskit(statevector: qiskit.quantum_info.Statevector) -> "Statevector":
    """
    Convert a Qiskit StateVector object to a scikit-q StateVector.
    Qiskit uses little-endian convention.
    :param statevector: Qiskit StateVector object
    :return: scikit-q StateVector object
    """
    return Statevector(statevector.data).reverse()

is_1d()

State vector is 1D.

Source code in src/skq/quantum_info/state.py 
29
30
31
def is_1d(self) -> bool:
    """State vector is 1D."""
    return self.ndim == 1

is_indistinguishable(other)

Two state vectors are indistinguishable if their density matrices are the same.

Source code in src/skq/quantum_info/state.py 
46
47
48
def is_indistinguishable(self, other: "Statevector") -> bool:
    """Two state vectors are indistinguishable if their density matrices are the same."""
    return np.allclose(self.density_matrix(), other.density_matrix())

is_multi_qubit()

State vector represents a multi-qubit state.

Source code in src/skq/quantum_info/state.py 
42
43
44
def is_multi_qubit(self) -> bool:
    """State vector represents a multi-qubit state."""
    return self.num_qubits > 1

is_normalized()

State vector is normalized: ||ψ|| = 1.

Source code in src/skq/quantum_info/state.py 
33
34
35
def is_normalized(self) -> bool:
    """State vector is normalized: ||ψ|| = 1."""
    return np.isclose(np.linalg.norm(self), 1)

is_power_of_two_len()

Check if a number is a power of two

Source code in src/skq/quantum_info/state.py 
37
38
39
40
def is_power_of_two_len(self) -> bool:
    """Check if a number is a power of two"""
    n = len(self)
    return n >= 2 and (n & (n - 1)) == 0

magnitude()

Magnitude (or norm) of the state vector. sqrt(<ψ|ψ>)

Source code in src/skq/quantum_info/state.py 
50
51
52
def magnitude(self) -> float:
    """Magnitude (or norm) of the state vector. sqrt(<ψ|ψ>)"""
    np.linalg.norm(self)

measure_bitstring()

Simulate a measurement of the state vector and get a bitstring representing the sample. :return: Bitstring representation of the measured state. For example: |0> -> "0" |00> -> "00" |11> -> "11"

Source code in src/skq/quantum_info/state.py 
84
85
86
87
88
89
90
91
92
93
def measure_bitstring(self) -> str:
    """
    Simulate a measurement of the state vector and get a bitstring representing the sample.
    :return: Bitstring representation of the measured state.
    For example:
    |0> -> "0"
    |00> -> "00"
    |11> -> "11"
    """
    return bin(self.measure_index())[2:].zfill(self.num_qubits)

measure_index()

Simulate a measurement of the state and get a sampled index. :return: Index of the measured state. For example: |0> -> 0 |00> -> 0 |11> -> 3

Source code in src/skq/quantum_info/state.py 
73
74
75
76
77
78
79
80
81
82
def measure_index(self) -> int:
    """
    Simulate a measurement of the state and get a sampled index.
    :return: Index of the measured state.
    For example:
    |0> -> 0
    |00> -> 0
    |11> -> 3
    """
    return np.random.choice(len(self), p=self.probabilities())

orthonormal_basis()

Orthonormal basis using the Gram-Schmidt process. :return: 2D array representing the orthonormal basis.

Source code in src/skq/quantum_info/state.py 
111
112
113
114
115
116
def orthonormal_basis(self) -> np.ndarray:
    """
    Orthonormal basis using the Gram-Schmidt process.
    :return: 2D array representing the orthonormal basis.
    """
    return np.array([self / np.linalg.norm(self)]).T

probabilities()

Probabilities of all possible states.

Source code in src/skq/quantum_info/state.py 
69
70
71
def probabilities(self) -> np.ndarray:
    """Probabilities of all possible states."""
    return np.abs(self) ** 2

reverse()

Reverse the order of the state vector to account for endianness. For example Qiskit uses little endian convention. Little-endian -> Least significant qubit (LSB) is on the right. Like |q1 q0> where q0 is the LSB. Big-endian -> Least significant qubit (LSB) is on the left. Like |q0 q1> where q0 is the LSB.

Source code in src/skq/quantum_info/state.py 
 95
 96
 97
 98
 99
100
101
def reverse(self) -> "Statevector":
    """Reverse the order of the state vector to account for endianness.
    For example Qiskit uses little endian convention.
    Little-endian -> Least significant qubit (LSB) is on the right. Like |q1 q0> where q0 is the LSB.
    Big-endian -> Least significant qubit (LSB) is on the left. Like |q0 q1> where q0 is the LSB.
    """
    return Statevector(self[::-1])

schmidt_decomposition()

Perform Schmidt decomposition on a quantum state. :return: Tuple of Schmidt coefficients, Basis A and Basis B

Source code in src/skq/quantum_info/state.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def schmidt_decomposition(self) -> tuple[np.array, np.array, np.array]:
    """
    Perform Schmidt decomposition on a quantum state.
    :return: Tuple of Schmidt coefficients, Basis A and Basis B
    """
    # Infer dimensions
    N = len(self)
    dim_A = int(np.sqrt(N))
    dim_B = N // dim_A

    # Singular Value Decomposition
    state_matrix = self.reshape(dim_A, dim_B)
    U, S, Vh = svd(state_matrix)

    # Coefficients (S), Basis A (U) and Basis B (Vh^dagger)
    return S, U, Vh.conj().T

to_pyquil()

Convert the state vector to a PyQuil object. PyQuil uses little-endian convention. :return: PyQuil object

Source code in src/skq/quantum_info/state.py 
153
154
155
156
157
158
159
def to_pyquil(self):
    """
    Convert the state vector to a PyQuil object.
    PyQuil uses little-endian convention.
    :return: PyQuil object
    """
    raise NotImplementedError("Conversion to PyQuil is not implemented.")

to_qiskit()

Convert the state vector to a Qiskit QuantumCircuit object. Qiskit uses little-endian convention. :return: Qiskit StateVector object

Source code in src/skq/quantum_info/state.py 
135
136
137
138
139
140
141
def to_qiskit(self) -> qiskit.quantum_info.Statevector:
    """
    Convert the state vector to a Qiskit QuantumCircuit object.
    Qiskit uses little-endian convention.
    :return: Qiskit StateVector object
    """
    return qiskit.quantum_info.Statevector(self.reverse())

Predefined States

SKQ provides several predefined quantum states:

skq.quantum_info.ZeroState

Bases: Statevector

Zero state |0...0>

Source code in src/skq/quantum_info/state.py 
171
172
173
174
175
class ZeroState(Statevector):
    """Zero state |0...0>"""

    def __new__(cls, num_qubits: int):
        return super().__new__(cls, [1] + [0] * (2**num_qubits - 1))

skq.quantum_info.OneState

Bases: Statevector

One state |1...1>

Source code in src/skq/quantum_info/state.py 
178
179
180
181
182
class OneState(Statevector):
    """One state |1...1>"""

    def __new__(cls, num_qubits: int):
        return super().__new__(cls, [0] * (2**num_qubits - 1) + [1])

skq.quantum_info.PlusState

Bases: Statevector

Single Qubit |+> superposition state

Source code in src/skq/quantum_info/state.py 
185
186
187
188
189
class PlusState(Statevector):
    """Single Qubit |+> superposition state"""

    def __new__(cls):
        return super().__new__(cls, [1, 1] / np.sqrt(2))

skq.quantum_info.MinusState

Bases: Statevector

Single Qubit |-> superposition state

Source code in src/skq/quantum_info/state.py 
192
193
194
195
196
class MinusState(Statevector):
    """Single Qubit |-> superposition state"""

    def __new__(cls):
        return super().__new__(cls, [1, -1] / np.sqrt(2))

Bell States

skq.quantum_info.PhiPlusState

Bases: Statevector

Bell state |Φ+>

Source code in src/skq/quantum_info/state.py 
199
200
201
202
203
class PhiPlusState(Statevector):
    """Bell state |Φ+>"""

    def __new__(cls):
        return super().__new__(cls, [1, 0, 0, 1] / np.sqrt(2))

skq.quantum_info.PhiMinusState

Bases: Statevector

Bell state |Φ->

Source code in src/skq/quantum_info/state.py 
206
207
208
209
210
class PhiMinusState(Statevector):
    """Bell state |Φ->"""

    def __new__(cls):
        return super().__new__(cls, [1, 0, 0, -1] / np.sqrt(2))

skq.quantum_info.PsiPlusState

Bases: Statevector

Bell state |Ψ+>

Source code in src/skq/quantum_info/state.py 
213
214
215
216
217
class PsiPlusState(Statevector):
    """Bell state |Ψ+>"""

    def __new__(cls):
        return super().__new__(cls, [0, 1, 1, 0] / np.sqrt(2))

skq.quantum_info.PsiMinusState

Bases: Statevector

Bell state |Ψ->

Source code in src/skq/quantum_info/state.py 
220
221
222
223
224
class PsiMinusState(Statevector):
    """Bell state |Ψ->"""

    def __new__(cls):
        return super().__new__(cls, [0, 1, -1, 0] / np.sqrt(2))

Multi-qubit States

skq.quantum_info.GHZState

Bases: Statevector

GHZ state |0...0> + |1...1>

Source code in src/skq/quantum_info/state.py 
227
228
229
230
231
232
233
234
235
class GHZState(Statevector):
    """GHZ state |0...0> + |1...1>"""

    def __new__(cls, num_qubits: int):
        assert num_qubits >= 3, "GHZ state requires at least 3 qubits."
        state = np.zeros(2**num_qubits)
        state[0] = 1 / np.sqrt(2)
        state[-1] = 1 / np.sqrt(2)
        return super().__new__(cls, state)

skq.quantum_info.WState

Bases: Statevector

W state |001> + |010> + |100>

Source code in src/skq/quantum_info/state.py 
238
239
240
241
242
243
244
245
246
class WState(Statevector):
    """W state |001> + |010> + |100>"""

    def __new__(cls, num_qubits: int):
        assert num_qubits >= 3, "W state requires at least 3 qubits."
        state = np.zeros(2**num_qubits)
        for i in range(num_qubits):
            state[2**i] = 1 / np.sqrt(num_qubits)
        return super().__new__(cls, state)

Density Matrices

Density matrices represent both pure and mixed quantum states.

from skq.quantum_info import DensityMatrix

# Create a density matrix
rho = DensityMatrix(np.array([[0.5, 0], [0, 0.5]]))
print(rho.is_pure())  # False
print(rho.is_mixed())  # True

skq.quantum_info.DensityMatrix

Bases: HermitianOperator

Density matrix representation of a qubit state.

Source code in src/skq/quantum_info/density.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
class DensityMatrix(HermitianOperator):
    """Density matrix representation of a qubit state."""

    def __new__(cls, input_array):
        obj = super().__new__(cls, input_array)
        assert obj.is_positive_semidefinite(), "Density matrix must be positive semidefinite (All eigenvalues >= 0)."
        assert obj.trace_equal_to_one(), "Density matrix must have trace equal to one. Normalize to unit trace if you want to use this matrix as a DensityMatrix."
        return obj

    @property
    def num_qubits(self) -> int:
        """Number of qubits in the density matrix."""
        return int(np.log2(len(self)))

    def is_positive_semidefinite(self):
        """Matrix is positive semidefinite."""
        eigenvalues = np.linalg.eigvalsh(self)
        return np.all(eigenvalues >= 0)

    def is_pure(self) -> bool:
        """Check if density matrix is pure."""
        return np.isclose(np.trace(self @ self), 1)

    def is_mixed(self) -> bool:
        """Check if density matrix is mixed."""
        return not self.is_pure()

    def trace_equal_to_one(self) -> bool:
        """Trace of density matrix is equal to one."""
        return np.isclose(np.trace(self), 1)

    def probabilities(self) -> float:
        """Probabilities of all possible state measurements."""
        return np.diag(self).real

    def is_multi_qubit(self) -> bool:
        """Check if the density matrix represents a multi-qubit state."""
        return self.num_qubits > 1

    def trace_norm(self) -> float:
        """Trace norm of the density matrix."""
        return np.trace(np.sqrt(self.conjugate_transpose() @ self))

    def distance(self, other: "DensityMatrix") -> float:
        """Trace norm distance between two density matrices."""
        assert isinstance(other, DensityMatrix), "'other' argument must be a valid DensityMatrix object."
        return self.trace_norm(self - other)

    def bloch_vector(self) -> np.ndarray:
        """Bloch vector of the density matrix."""
        if self.num_qubits > 1:
            raise NotImplementedError("Bloch vector is not yet implemented for multi-qubit states.")

        # Bloch vector components
        bx = np.trace(np.dot(self, X())).real
        by = np.trace(np.dot(self, Y())).real
        bz = np.trace(np.dot(self, Z())).real
        return np.array([bx, by, bz])

    def kron(self, other: "DensityMatrix") -> "DensityMatrix":
        """
        Kronecker (tensor) product of two density matrices.
        This can be used to create so-called "product states" that represent
        the independence between two quantum systems.
        :param other: DensityMatrix object
        :return: Kronecker product of the two density matrices
        """
        return DensityMatrix(np.kron(self, other))

    def entropy(self):
        """von Neumann entropy."""
        eigenvals = self.eigenvalues()
        # Only consider non-zero eigenvalues
        nonzero_eigenvalues = eigenvals[eigenvals > 0]
        return -np.sum(nonzero_eigenvalues * np.log(nonzero_eigenvalues))

    def internal_energy(self, hamiltonian: np.array) -> float:
        """
        Expected value of the Hamiltonian.
        :param hamiltonian: Hamiltonian matrix
        :return: Expected value of the Hamiltonian
        """
        return np.trace(self @ hamiltonian)

    def to_qiskit(self) -> qiskit.quantum_info.DensityMatrix:
        """
        Convert the density matrix to a Qiskit DensityMatrix object.
        :return: Qiskit DensityMatrix object
        """
        return qiskit.quantum_info.DensityMatrix(self)

    @staticmethod
    def from_qiskit(density_matrix: qiskit.quantum_info.DensityMatrix) -> "DensityMatrix":
        """
        Create a DensityMatrix object from a Qiskit DensityMatrix object.
        :param density_matrix: Qiskit DensityMatrix object
        :return: DensityMatrix object
        """
        return DensityMatrix(density_matrix.data)

    def to_pennylane(self, wires: list[int] | int = None) -> qml.QubitDensityMatrix:
        """
        Convert the density matrix to a PennyLane QubitDensityMatrix.
        :param wires: List of wires to apply the density matrix to
        :return: PennyLane QubitDensityMatrix object
        """
        wires = wires if wires is not None else range(self.num_qubits)
        return qml.QubitDensityMatrix(self, wires=wires)

    @staticmethod
    def from_pennylane(density_matrix: qml.QubitDensityMatrix) -> "DensityMatrix":
        """
        Convert a PennyLane QubitDensityMarix object to a scikit-q StateVector.
        :param density_matrix: PennyLane QubitDensityMatrix object
        :return: scikit-q StateVector object
        """
        return DensityMatrix(density_matrix.data[0])

    @staticmethod
    def from_probabilities(probabilities: np.array) -> "DensityMatrix":
        """
        Create a density matrix from a list of probabilities.
        :param probabilities: A 1D array of probabilities
        :return: Density matrix
        """
        assert np.isclose(np.sum(probabilities), 1), f"Probabilities must sum to one. Got sum: {np.sum(probabilities)}"
        assert len(probabilities.shape) == 1, f"Probabilities must be a 1D array. Got shape: {probabilities.shape}"
        return DensityMatrix(np.diag(probabilities))

num_qubits property

Number of qubits in the density matrix.

bloch_vector()

Bloch vector of the density matrix.

Source code in src/skq/quantum_info/density.py 
59
60
61
62
63
64
65
66
67
68
def bloch_vector(self) -> np.ndarray:
    """Bloch vector of the density matrix."""
    if self.num_qubits > 1:
        raise NotImplementedError("Bloch vector is not yet implemented for multi-qubit states.")

    # Bloch vector components
    bx = np.trace(np.dot(self, X())).real
    by = np.trace(np.dot(self, Y())).real
    bz = np.trace(np.dot(self, Z())).real
    return np.array([bx, by, bz])

distance(other)

Trace norm distance between two density matrices.

Source code in src/skq/quantum_info/density.py 
54
55
56
57
def distance(self, other: "DensityMatrix") -> float:
    """Trace norm distance between two density matrices."""
    assert isinstance(other, DensityMatrix), "'other' argument must be a valid DensityMatrix object."
    return self.trace_norm(self - other)

entropy()

von Neumann entropy.

Source code in src/skq/quantum_info/density.py 
80
81
82
83
84
85
def entropy(self):
    """von Neumann entropy."""
    eigenvals = self.eigenvalues()
    # Only consider non-zero eigenvalues
    nonzero_eigenvalues = eigenvals[eigenvals > 0]
    return -np.sum(nonzero_eigenvalues * np.log(nonzero_eigenvalues))

from_pennylane(density_matrix) staticmethod

Convert a PennyLane QubitDensityMarix object to a scikit-q StateVector. :param density_matrix: PennyLane QubitDensityMatrix object :return: scikit-q StateVector object

Source code in src/skq/quantum_info/density.py 
120
121
122
123
124
125
126
127
@staticmethod
def from_pennylane(density_matrix: qml.QubitDensityMatrix) -> "DensityMatrix":
    """
    Convert a PennyLane QubitDensityMarix object to a scikit-q StateVector.
    :param density_matrix: PennyLane QubitDensityMatrix object
    :return: scikit-q StateVector object
    """
    return DensityMatrix(density_matrix.data[0])

from_probabilities(probabilities) staticmethod

Create a density matrix from a list of probabilities. :param probabilities: A 1D array of probabilities :return: Density matrix

Source code in src/skq/quantum_info/density.py 
129
130
131
132
133
134
135
136
137
138
@staticmethod
def from_probabilities(probabilities: np.array) -> "DensityMatrix":
    """
    Create a density matrix from a list of probabilities.
    :param probabilities: A 1D array of probabilities
    :return: Density matrix
    """
    assert np.isclose(np.sum(probabilities), 1), f"Probabilities must sum to one. Got sum: {np.sum(probabilities)}"
    assert len(probabilities.shape) == 1, f"Probabilities must be a 1D array. Got shape: {probabilities.shape}"
    return DensityMatrix(np.diag(probabilities))

from_qiskit(density_matrix) staticmethod

Create a DensityMatrix object from a Qiskit DensityMatrix object. :param density_matrix: Qiskit DensityMatrix object :return: DensityMatrix object

Source code in src/skq/quantum_info/density.py 
102
103
104
105
106
107
108
109
@staticmethod
def from_qiskit(density_matrix: qiskit.quantum_info.DensityMatrix) -> "DensityMatrix":
    """
    Create a DensityMatrix object from a Qiskit DensityMatrix object.
    :param density_matrix: Qiskit DensityMatrix object
    :return: DensityMatrix object
    """
    return DensityMatrix(density_matrix.data)

internal_energy(hamiltonian)

Expected value of the Hamiltonian. :param hamiltonian: Hamiltonian matrix :return: Expected value of the Hamiltonian

Source code in src/skq/quantum_info/density.py 
87
88
89
90
91
92
93
def internal_energy(self, hamiltonian: np.array) -> float:
    """
    Expected value of the Hamiltonian.
    :param hamiltonian: Hamiltonian matrix
    :return: Expected value of the Hamiltonian
    """
    return np.trace(self @ hamiltonian)

is_mixed()

Check if density matrix is mixed.

Source code in src/skq/quantum_info/density.py 
34
35
36
def is_mixed(self) -> bool:
    """Check if density matrix is mixed."""
    return not self.is_pure()

is_multi_qubit()

Check if the density matrix represents a multi-qubit state.

Source code in src/skq/quantum_info/density.py 
46
47
48
def is_multi_qubit(self) -> bool:
    """Check if the density matrix represents a multi-qubit state."""
    return self.num_qubits > 1

is_positive_semidefinite()

Matrix is positive semidefinite.

Source code in src/skq/quantum_info/density.py 
25
26
27
28
def is_positive_semidefinite(self):
    """Matrix is positive semidefinite."""
    eigenvalues = np.linalg.eigvalsh(self)
    return np.all(eigenvalues >= 0)

is_pure()

Check if density matrix is pure.

Source code in src/skq/quantum_info/density.py 
30
31
32
def is_pure(self) -> bool:
    """Check if density matrix is pure."""
    return np.isclose(np.trace(self @ self), 1)

kron(other)

Kronecker (tensor) product of two density matrices. This can be used to create so-called "product states" that represent the independence between two quantum systems. :param other: DensityMatrix object :return: Kronecker product of the two density matrices

Source code in src/skq/quantum_info/density.py 
70
71
72
73
74
75
76
77
78
def kron(self, other: "DensityMatrix") -> "DensityMatrix":
    """
    Kronecker (tensor) product of two density matrices.
    This can be used to create so-called "product states" that represent
    the independence between two quantum systems.
    :param other: DensityMatrix object
    :return: Kronecker product of the two density matrices
    """
    return DensityMatrix(np.kron(self, other))

probabilities()

Probabilities of all possible state measurements.

Source code in src/skq/quantum_info/density.py 
42
43
44
def probabilities(self) -> float:
    """Probabilities of all possible state measurements."""
    return np.diag(self).real

to_pennylane(wires=None)

Convert the density matrix to a PennyLane QubitDensityMatrix. :param wires: List of wires to apply the density matrix to :return: PennyLane QubitDensityMatrix object

Source code in src/skq/quantum_info/density.py 
111
112
113
114
115
116
117
118
def to_pennylane(self, wires: list[int] | int = None) -> qml.QubitDensityMatrix:
    """
    Convert the density matrix to a PennyLane QubitDensityMatrix.
    :param wires: List of wires to apply the density matrix to
    :return: PennyLane QubitDensityMatrix object
    """
    wires = wires if wires is not None else range(self.num_qubits)
    return qml.QubitDensityMatrix(self, wires=wires)

to_qiskit()

Convert the density matrix to a Qiskit DensityMatrix object. :return: Qiskit DensityMatrix object

Source code in src/skq/quantum_info/density.py 
 95
 96
 97
 98
 99
100
def to_qiskit(self) -> qiskit.quantum_info.DensityMatrix:
    """
    Convert the density matrix to a Qiskit DensityMatrix object.
    :return: Qiskit DensityMatrix object
    """
    return qiskit.quantum_info.DensityMatrix(self)

trace_equal_to_one()

Trace of density matrix is equal to one.

Source code in src/skq/quantum_info/density.py 
38
39
40
def trace_equal_to_one(self) -> bool:
    """Trace of density matrix is equal to one."""
    return np.isclose(np.trace(self), 1)

trace_norm()

Trace norm of the density matrix.

Source code in src/skq/quantum_info/density.py 
50
51
52
def trace_norm(self) -> float:
    """Trace norm of the density matrix."""
    return np.trace(np.sqrt(self.conjugate_transpose() @ self))

Thermal States

skq.quantum_info.GibbsState

Bases: DensityMatrix

Gibbs (mixed) state representation of a quantum state in thermal equilibrium. :param hamiltonian: Hamiltonian matrix of the system :param temperature: Temperature of the system in Kelvin

Source code in src/skq/quantum_info/density.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
class GibbsState(DensityMatrix):
    """
    Gibbs (mixed) state representation of a quantum state in thermal equilibrium.
    :param hamiltonian: Hamiltonian matrix of the system
    :param temperature: Temperature of the system in Kelvin
    """

    def __new__(cls, hamiltonian: np.array, temperature: float):
        cls.hamiltonian = hamiltonian
        cls.temperature = temperature
        cls.beta = 1 / (BOLTZMANN_CONSTANT * temperature)
        cls.exp_neg_beta_H = expm(-cls.beta * hamiltonian)
        cls.partition_function = np.trace(cls.exp_neg_beta_H)
        density_matrix = cls.exp_neg_beta_H / cls.partition_function
        return super().__new__(cls, density_matrix)

    def free_energy(self) -> float:
        """Helmholtz free energy."""
        return -BOLTZMANN_CONSTANT * self.temperature * np.log(self.partition_function)

    def heat_capacity(self) -> float:
        """Calculate the heat capacity."""
        beta = 1 / (BOLTZMANN_CONSTANT * self.temperature)
        energy_squared = np.trace(self @ self.hamiltonian @ self.hamiltonian)
        energy_mean = self.internal_energy(self.hamiltonian) ** 2
        return beta**2 * (energy_squared - energy_mean)

free_energy()

Helmholtz free energy.

Source code in src/skq/quantum_info/density.py 
157
158
159
def free_energy(self) -> float:
    """Helmholtz free energy."""
    return -BOLTZMANN_CONSTANT * self.temperature * np.log(self.partition_function)

heat_capacity()

Calculate the heat capacity.

Source code in src/skq/quantum_info/density.py 
161
162
163
164
165
166
def heat_capacity(self) -> float:
    """Calculate the heat capacity."""
    beta = 1 / (BOLTZMANN_CONSTANT * self.temperature)
    energy_squared = np.trace(self @ self.hamiltonian @ self.hamiltonian)
    energy_mean = self.internal_energy(self.hamiltonian) ** 2
    return beta**2 * (energy_squared - energy_mean)

Quantum Channels

Quantum channels represent physical operations on quantum states.

from skq.quantum_info import DepolarizingChannel

# Create a depolarizing channel with 10% noise
channel = DepolarizingChannel(0.1)

skq.quantum_info.QuantumChannel

Bases: SuperOperator

Quantum Channel representation in choi, stinespring, or kraus form.

Source code in src/skq/quantum_info/channel.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
class QuantumChannel(SuperOperator):
    """Quantum Channel representation in choi, stinespring, or kraus form."""

    def __new__(cls, input_array: np.array, representation="choi"):
        cls.representation = representation
        obj = super().__new__(cls, input_array)
        if cls.representation == "choi":
            obj._validate_choi()
        elif cls.representation == "stinespring":
            obj._validate_stinespring()
        elif cls.representation == "kraus":
            obj._validate_kraus()
        else:
            raise ValueError("Invalid representation. Choose from 'choi', 'stinespring', or 'kraus'.")
        return obj

    def _validate_choi(self):
        """Validate if the input matrix is a valid choi matrix."""
        assert self.is_nd(2), "Choi matrix must be a 2D matrix."
        assert self.is_square(), "Choi matrix must be a square matrix."
        assert self.is_positive_semidefinite(), "Choi matrix must be positive semidefinite."
        assert self.is_trace_preserving(), "Choi matrix must be trace-preserving."

    def _validate_stinespring(self):
        """Validate if the input matrix is a valid stinespring matrix."""
        assert self.is_nd(2), "Stringspring matrix must be a 2D matrix."
        assert self.is_isometry(), "Stinespring matrix must be an isometry."
        assert self.shape[0] > self.shape[1], "Stinespring matrix must have more rows than columns."

    def _validate_kraus(self):
        """Validate if the input list of matrices is a valid kraus representation."""
        assert self.is_nd(3), "Kraus operators must be 3D matrices."
        for kraus_op in self:
            assert kraus_op.shape[0] == kraus_op.shape[1], "Each Kraus operator must be a square matrix."
        assert self.is_trace_preserving(), "Sum of outer products of kraus operators must be identity."

    def is_nd(self, n: int) -> bool:
        """Channel is an n-dimensional matrix."""
        return self.ndim == n

    def is_square(self) -> bool:
        """Check if the superoperator is represented by a square matrix."""
        return self.shape[0] == self.shape[1]

    def is_positive_semidefinite(self) -> bool:
        """Check if channel is positive semidefinite (i.e. all eigenvalues are non-negative).
        Requirement for choi matrices."""
        eigenvalues = np.linalg.eigvalsh(self)
        return np.all(eigenvalues >= 0)

    def is_trace_preserving(self) -> bool:
        """Check if the quantum channel is trace-preserving."""
        if self.representation == "kraus":
            d = self[0].shape[0]
            kraus_sum = sum([np.dot(k.conjugate_transpose(), k) for k in self])
            return np.allclose(kraus_sum, np.eye(d))
        elif self.representation == "choi":
            d = int(np.sqrt(self.shape[0]))
            choi_reshaped = self.reshape(d, d, d, d)
            partial_trace = np.trace(choi_reshaped, axis1=1, axis2=3)
            return np.allclose(partial_trace, np.eye(d)) or np.isclose(np.trace(partial_trace), 1.0)
        else:
            raise ValueError(f"Trace-preserving check is not implemented for representation '{self.representation}'.")

    def is_isometry(self) -> bool:
        """
        Check if the quantum channel is an isometry.
        An isometry is a linear transformation that preserves distances (V^dagger V = I)
        """
        return np.allclose(self.conjugate_transpose() @ self, np.eye(self.shape[1]))

    def compose(self, other: "QuantumChannel") -> "QuantumChannel":
        """
        Compose quantum channels.
        :param other: QuantumChannel object to compose with.
        :return: Composed QuantumChannel object in the Kraus representation.
        """
        assert isinstance(other, QuantumChannel), "Input must be a QuantumChannel."
        self_kraus = self.to_kraus()
        other_kraus = other.to_kraus()
        composed_kraus = []

        for k1 in self_kraus:
            for k2 in other_kraus:
                composed_kraus.append(np.dot(k1, k2))

        composed_kraus = np.array(composed_kraus, dtype=complex)
        kraus_sum = sum(np.dot(k.conj().T, k) for k in composed_kraus)
        normalization_factor = np.linalg.inv(sqrtm(kraus_sum).astype(np.complex128))
        normalized_kraus = [np.dot(normalization_factor, k) for k in composed_kraus]
        return QuantumChannel(np.array(normalized_kraus), representation="kraus")

    def tensor(self, other: "QuantumChannel") -> "QuantumChannel":
        """
        Tensor product with another channel.
        :param other: QuantumChannel object to tensor with.
        :return: QuantumChannel object in the Choi representation
        """
        assert isinstance(other, QuantumChannel), "The other object must be a QuantumChannel."
        return QuantumChannel(np.kron(self.to_choi(), other.to_choi()), representation="choi")

    def fidelity(self, other: "QuantumChannel") -> float:
        """
        Fidelity between two quantum channels.
        :param other: QuantumChannel object to calculate fidelity with.
        :return: Fidelity in range [0...1].
        """
        assert isinstance(other, QuantumChannel), "The other object must be a QuantumChannel."
        choi_self = self.to_choi()
        choi_other = other.to_choi()
        norm_choi_self = choi_self / np.trace(choi_self)
        norm_choi_other = choi_other / np.trace(choi_other)
        sqrt_choi_self = sqrtm(norm_choi_self)

        product_matrix = sqrt_choi_self @ norm_choi_other @ sqrt_choi_self
        fidelity_value = np.trace(sqrtm(product_matrix)) ** 2
        fidelity_value = np.real(fidelity_value)
        fidelity_value = min(max(fidelity_value, 0), 1)
        return fidelity_value

    def to_choi(self):
        """Convert the channel to the choi matrix representation."""
        if self.representation == "choi":
            return self
        elif self.representation == "stinespring":
            return self._stinespring_to_choi()
        elif self.representation == "kraus":
            return self._kraus_to_choi()

    def to_stinespring(self):
        """Convert the channel to the stinespring representation."""
        if self.representation == "stinespring":
            return self
        elif self.representation == "choi":
            return self._choi_to_stinespring()
        elif self.representation == "kraus":
            return self._kraus_to_stinespring()

    def to_kraus(self):
        """Convert the channel to the kraus representation."""
        if self.representation == "kraus":
            return self
        elif self.representation == "choi":
            return self._choi_to_kraus()
        elif self.representation == "stinespring":
            return self._stinespring_to_kraus()

    def _stinespring_to_choi(self) -> "QuantumChannel":
        """Convert stinespring representation to choi matrix."""
        return QuantumChannel(np.dot(self, self.conjugate_transpose()), representation="choi")

    def _choi_to_stinespring(self) -> "QuantumChannel":
        """Convert choi matrix to stinespring representation."""
        d = int(np.sqrt(self.shape[0]))
        w, v = np.linalg.eigh(self)
        # Filter out small eigenvalues to avoid numerical instability
        non_zero_indices = np.where(w > 1e-10)[0]
        sqrt_eigenvals = np.sqrt(w[non_zero_indices])
        v = v[:, non_zero_indices]
        r = len(non_zero_indices)
        stinespring_matrix = np.zeros((d * r, d), dtype=complex)
        for i in range(r):
            stinespring_matrix[i * d : (i + 1) * d, :] = sqrt_eigenvals[i] * v[:, i].reshape(d, d)
        return QuantumChannel(stinespring_matrix, representation="stinespring")

    def _kraus_to_choi(self) -> "QuantumChannel":
        """Convert Kraus representation to Choi matrix using vectorization."""
        d = self[0].shape[0]
        choi_matrix = np.zeros((d * d, d * d), dtype=complex)
        for k in self:
            choi_matrix += np.kron(k, k.conj())
        return QuantumChannel(choi_matrix / d, representation="choi")

    def _choi_to_kraus(self) -> "QuantumChannel":
        """Convert choi matrix to kraus operators."""
        d = int(np.sqrt(self.shape[0]))
        kraus_operators = []
        w, v = np.linalg.eigh(self)
        for i in range(len(w)):
            if np.isclose(w[i], 0):
                continue
            kraus_operators.append(np.sqrt(w[i]) * v[:, i].reshape(d, d))
        return QuantumChannel(np.array(kraus_operators, dtype=complex), representation="kraus")

    def _kraus_to_stinespring(self) -> "QuantumChannel":
        """Convert kraus representation to stinespring representation."""
        d = self[0].shape[0]
        num_kraus = len(self)
        stinespring_matrix = np.zeros((d * num_kraus, d), dtype=complex)
        for i, k in enumerate(self):
            stinespring_matrix[i * d : (i + 1) * d, :] = k
        return QuantumChannel(stinespring_matrix, representation="stinespring")

    def _stinespring_to_kraus(self) -> "QuantumChannel":
        """Convert stinespring representation to kraus operators."""
        d = self.shape[1]
        num_kraus = self.shape[0] // d
        kraus_operators = []
        for i in range(num_kraus):
            kraus_operators.append(self[i * d : (i + 1) * d, :])
        return QuantumChannel(np.array(kraus_operators, dtype=complex), representation="kraus")

    def __call__(self, density_matrix: DensityMatrix) -> DensityMatrix:
        rho_out = np.zeros_like(density_matrix, dtype=complex)
        for K in self.to_kraus():
            rho_out += K @ density_matrix @ K.conj().T
        return DensityMatrix(rho_out)

    def to_qiskit(self) -> qiskit.quantum_info.operators.channel.quantum_channel.QuantumChannel:
        """Convert the channel to a Qiskit object."""
        if self.representation == "kraus":
            # Input must be a list of numpy arrays
            kraus_list = [self[i] for i in range(self.shape[0])]
            return qiskit.quantum_info.Kraus(kraus_list)
        elif self.representation == "stinespring":
            return qiskit.quantum_info.Stinespring(self)
        elif self.representation == "choi":
            return qiskit.quantum_info.Choi(self)

    def from_qiskit(self, channel: qiskit.quantum_info.operators.channel.quantum_channel.QuantumChannel) -> "QuantumChannel":
        """Convert a Qiskit channel to a QuantumChannel object."""
        if isinstance(channel, qiskit.quantum_info.Kraus):
            return QuantumChannel(channel.data, representation="kraus")
        elif isinstance(channel, qiskit.quantum_info.Stinespring):
            return QuantumChannel(channel.data, representation="stinespring")
        elif isinstance(channel, qiskit.quantum_info.Choi):
            return QuantumChannel(channel.data, representation="choi")
        else:
            raise ValueError(f"Invalid Qiskit channel type '{channel}'.")

_choi_to_kraus()

Convert choi matrix to kraus operators.

Source code in src/skq/quantum_info/channel.py 
183
184
185
186
187
188
189
190
191
192
def _choi_to_kraus(self) -> "QuantumChannel":
    """Convert choi matrix to kraus operators."""
    d = int(np.sqrt(self.shape[0]))
    kraus_operators = []
    w, v = np.linalg.eigh(self)
    for i in range(len(w)):
        if np.isclose(w[i], 0):
            continue
        kraus_operators.append(np.sqrt(w[i]) * v[:, i].reshape(d, d))
    return QuantumChannel(np.array(kraus_operators, dtype=complex), representation="kraus")

_choi_to_stinespring()

Convert choi matrix to stinespring representation.

Source code in src/skq/quantum_info/channel.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
def _choi_to_stinespring(self) -> "QuantumChannel":
    """Convert choi matrix to stinespring representation."""
    d = int(np.sqrt(self.shape[0]))
    w, v = np.linalg.eigh(self)
    # Filter out small eigenvalues to avoid numerical instability
    non_zero_indices = np.where(w > 1e-10)[0]
    sqrt_eigenvals = np.sqrt(w[non_zero_indices])
    v = v[:, non_zero_indices]
    r = len(non_zero_indices)
    stinespring_matrix = np.zeros((d * r, d), dtype=complex)
    for i in range(r):
        stinespring_matrix[i * d : (i + 1) * d, :] = sqrt_eigenvals[i] * v[:, i].reshape(d, d)
    return QuantumChannel(stinespring_matrix, representation="stinespring")

_kraus_to_choi()

Convert Kraus representation to Choi matrix using vectorization.

Source code in src/skq/quantum_info/channel.py 
175
176
177
178
179
180
181
def _kraus_to_choi(self) -> "QuantumChannel":
    """Convert Kraus representation to Choi matrix using vectorization."""
    d = self[0].shape[0]
    choi_matrix = np.zeros((d * d, d * d), dtype=complex)
    for k in self:
        choi_matrix += np.kron(k, k.conj())
    return QuantumChannel(choi_matrix / d, representation="choi")

_kraus_to_stinespring()

Convert kraus representation to stinespring representation.

Source code in src/skq/quantum_info/channel.py 
194
195
196
197
198
199
200
201
def _kraus_to_stinespring(self) -> "QuantumChannel":
    """Convert kraus representation to stinespring representation."""
    d = self[0].shape[0]
    num_kraus = len(self)
    stinespring_matrix = np.zeros((d * num_kraus, d), dtype=complex)
    for i, k in enumerate(self):
        stinespring_matrix[i * d : (i + 1) * d, :] = k
    return QuantumChannel(stinespring_matrix, representation="stinespring")

_stinespring_to_choi()

Convert stinespring representation to choi matrix.

Source code in src/skq/quantum_info/channel.py 
157
158
159
def _stinespring_to_choi(self) -> "QuantumChannel":
    """Convert stinespring representation to choi matrix."""
    return QuantumChannel(np.dot(self, self.conjugate_transpose()), representation="choi")

_stinespring_to_kraus()

Convert stinespring representation to kraus operators.

Source code in src/skq/quantum_info/channel.py 
203
204
205
206
207
208
209
210
def _stinespring_to_kraus(self) -> "QuantumChannel":
    """Convert stinespring representation to kraus operators."""
    d = self.shape[1]
    num_kraus = self.shape[0] // d
    kraus_operators = []
    for i in range(num_kraus):
        kraus_operators.append(self[i * d : (i + 1) * d, :])
    return QuantumChannel(np.array(kraus_operators, dtype=complex), representation="kraus")

_validate_choi()

Validate if the input matrix is a valid choi matrix.

Source code in src/skq/quantum_info/channel.py 
26
27
28
29
30
31
def _validate_choi(self):
    """Validate if the input matrix is a valid choi matrix."""
    assert self.is_nd(2), "Choi matrix must be a 2D matrix."
    assert self.is_square(), "Choi matrix must be a square matrix."
    assert self.is_positive_semidefinite(), "Choi matrix must be positive semidefinite."
    assert self.is_trace_preserving(), "Choi matrix must be trace-preserving."

_validate_kraus()

Validate if the input list of matrices is a valid kraus representation.

Source code in src/skq/quantum_info/channel.py 
39
40
41
42
43
44
def _validate_kraus(self):
    """Validate if the input list of matrices is a valid kraus representation."""
    assert self.is_nd(3), "Kraus operators must be 3D matrices."
    for kraus_op in self:
        assert kraus_op.shape[0] == kraus_op.shape[1], "Each Kraus operator must be a square matrix."
    assert self.is_trace_preserving(), "Sum of outer products of kraus operators must be identity."

_validate_stinespring()

Validate if the input matrix is a valid stinespring matrix.

Source code in src/skq/quantum_info/channel.py 
33
34
35
36
37
def _validate_stinespring(self):
    """Validate if the input matrix is a valid stinespring matrix."""
    assert self.is_nd(2), "Stringspring matrix must be a 2D matrix."
    assert self.is_isometry(), "Stinespring matrix must be an isometry."
    assert self.shape[0] > self.shape[1], "Stinespring matrix must have more rows than columns."

compose(other)

Compose quantum channels. :param other: QuantumChannel object to compose with. :return: Composed QuantumChannel object in the Kraus representation.

Source code in src/skq/quantum_info/channel.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def compose(self, other: "QuantumChannel") -> "QuantumChannel":
    """
    Compose quantum channels.
    :param other: QuantumChannel object to compose with.
    :return: Composed QuantumChannel object in the Kraus representation.
    """
    assert isinstance(other, QuantumChannel), "Input must be a QuantumChannel."
    self_kraus = self.to_kraus()
    other_kraus = other.to_kraus()
    composed_kraus = []

    for k1 in self_kraus:
        for k2 in other_kraus:
            composed_kraus.append(np.dot(k1, k2))

    composed_kraus = np.array(composed_kraus, dtype=complex)
    kraus_sum = sum(np.dot(k.conj().T, k) for k in composed_kraus)
    normalization_factor = np.linalg.inv(sqrtm(kraus_sum).astype(np.complex128))
    normalized_kraus = [np.dot(normalization_factor, k) for k in composed_kraus]
    return QuantumChannel(np.array(normalized_kraus), representation="kraus")

fidelity(other)

Fidelity between two quantum channels. :param other: QuantumChannel object to calculate fidelity with. :return: Fidelity in range [0...1].

Source code in src/skq/quantum_info/channel.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def fidelity(self, other: "QuantumChannel") -> float:
    """
    Fidelity between two quantum channels.
    :param other: QuantumChannel object to calculate fidelity with.
    :return: Fidelity in range [0...1].
    """
    assert isinstance(other, QuantumChannel), "The other object must be a QuantumChannel."
    choi_self = self.to_choi()
    choi_other = other.to_choi()
    norm_choi_self = choi_self / np.trace(choi_self)
    norm_choi_other = choi_other / np.trace(choi_other)
    sqrt_choi_self = sqrtm(norm_choi_self)

    product_matrix = sqrt_choi_self @ norm_choi_other @ sqrt_choi_self
    fidelity_value = np.trace(sqrtm(product_matrix)) ** 2
    fidelity_value = np.real(fidelity_value)
    fidelity_value = min(max(fidelity_value, 0), 1)
    return fidelity_value

from_qiskit(channel)

Convert a Qiskit channel to a QuantumChannel object.

Source code in src/skq/quantum_info/channel.py 
229
230
231
232
233
234
235
236
237
238
def from_qiskit(self, channel: qiskit.quantum_info.operators.channel.quantum_channel.QuantumChannel) -> "QuantumChannel":
    """Convert a Qiskit channel to a QuantumChannel object."""
    if isinstance(channel, qiskit.quantum_info.Kraus):
        return QuantumChannel(channel.data, representation="kraus")
    elif isinstance(channel, qiskit.quantum_info.Stinespring):
        return QuantumChannel(channel.data, representation="stinespring")
    elif isinstance(channel, qiskit.quantum_info.Choi):
        return QuantumChannel(channel.data, representation="choi")
    else:
        raise ValueError(f"Invalid Qiskit channel type '{channel}'.")

is_isometry()

Check if the quantum channel is an isometry. An isometry is a linear transformation that preserves distances (V^dagger V = I)

Source code in src/skq/quantum_info/channel.py 
74
75
76
77
78
79
def is_isometry(self) -> bool:
    """
    Check if the quantum channel is an isometry.
    An isometry is a linear transformation that preserves distances (V^dagger V = I)
    """
    return np.allclose(self.conjugate_transpose() @ self, np.eye(self.shape[1]))

is_nd(n)

Channel is an n-dimensional matrix.

Source code in src/skq/quantum_info/channel.py 
46
47
48
def is_nd(self, n: int) -> bool:
    """Channel is an n-dimensional matrix."""
    return self.ndim == n

is_positive_semidefinite()

Check if channel is positive semidefinite (i.e. all eigenvalues are non-negative). Requirement for choi matrices.

Source code in src/skq/quantum_info/channel.py 
54
55
56
57
58
def is_positive_semidefinite(self) -> bool:
    """Check if channel is positive semidefinite (i.e. all eigenvalues are non-negative).
    Requirement for choi matrices."""
    eigenvalues = np.linalg.eigvalsh(self)
    return np.all(eigenvalues >= 0)

is_square()

Check if the superoperator is represented by a square matrix.

Source code in src/skq/quantum_info/channel.py 
50
51
52
def is_square(self) -> bool:
    """Check if the superoperator is represented by a square matrix."""
    return self.shape[0] == self.shape[1]

is_trace_preserving()

Check if the quantum channel is trace-preserving.

Source code in src/skq/quantum_info/channel.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
def is_trace_preserving(self) -> bool:
    """Check if the quantum channel is trace-preserving."""
    if self.representation == "kraus":
        d = self[0].shape[0]
        kraus_sum = sum([np.dot(k.conjugate_transpose(), k) for k in self])
        return np.allclose(kraus_sum, np.eye(d))
    elif self.representation == "choi":
        d = int(np.sqrt(self.shape[0]))
        choi_reshaped = self.reshape(d, d, d, d)
        partial_trace = np.trace(choi_reshaped, axis1=1, axis2=3)
        return np.allclose(partial_trace, np.eye(d)) or np.isclose(np.trace(partial_trace), 1.0)
    else:
        raise ValueError(f"Trace-preserving check is not implemented for representation '{self.representation}'.")

tensor(other)

Tensor product with another channel. :param other: QuantumChannel object to tensor with. :return: QuantumChannel object in the Choi representation

Source code in src/skq/quantum_info/channel.py 
102
103
104
105
106
107
108
109
def tensor(self, other: "QuantumChannel") -> "QuantumChannel":
    """
    Tensor product with another channel.
    :param other: QuantumChannel object to tensor with.
    :return: QuantumChannel object in the Choi representation
    """
    assert isinstance(other, QuantumChannel), "The other object must be a QuantumChannel."
    return QuantumChannel(np.kron(self.to_choi(), other.to_choi()), representation="choi")

to_choi()

Convert the channel to the choi matrix representation.

Source code in src/skq/quantum_info/channel.py 
130
131
132
133
134
135
136
137
def to_choi(self):
    """Convert the channel to the choi matrix representation."""
    if self.representation == "choi":
        return self
    elif self.representation == "stinespring":
        return self._stinespring_to_choi()
    elif self.representation == "kraus":
        return self._kraus_to_choi()

to_kraus()

Convert the channel to the kraus representation.

Source code in src/skq/quantum_info/channel.py 
148
149
150
151
152
153
154
155
def to_kraus(self):
    """Convert the channel to the kraus representation."""
    if self.representation == "kraus":
        return self
    elif self.representation == "choi":
        return self._choi_to_kraus()
    elif self.representation == "stinespring":
        return self._stinespring_to_kraus()

to_qiskit()

Convert the channel to a Qiskit object.

Source code in src/skq/quantum_info/channel.py 
218
219
220
221
222
223
224
225
226
227
def to_qiskit(self) -> qiskit.quantum_info.operators.channel.quantum_channel.QuantumChannel:
    """Convert the channel to a Qiskit object."""
    if self.representation == "kraus":
        # Input must be a list of numpy arrays
        kraus_list = [self[i] for i in range(self.shape[0])]
        return qiskit.quantum_info.Kraus(kraus_list)
    elif self.representation == "stinespring":
        return qiskit.quantum_info.Stinespring(self)
    elif self.representation == "choi":
        return qiskit.quantum_info.Choi(self)

to_stinespring()

Convert the channel to the stinespring representation.

Source code in src/skq/quantum_info/channel.py 
139
140
141
142
143
144
145
146
def to_stinespring(self):
    """Convert the channel to the stinespring representation."""
    if self.representation == "stinespring":
        return self
    elif self.representation == "choi":
        return self._choi_to_stinespring()
    elif self.representation == "kraus":
        return self._kraus_to_stinespring()

Predefined Channels

skq.quantum_info.QubitResetChannel

Bases: QuantumChannel

Reset channel for a qubit system. Initialized in the Kraus representation.

Source code in src/skq/quantum_info/channel.py 
241
242
243
244
245
246
247
248
249
class QubitResetChannel(QuantumChannel):
    """
    Reset channel for a qubit system.
    Initialized in the Kraus representation.
    """

    def __new__(cls):
        kraus_ops = np.array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
        return super().__new__(cls, kraus_ops, representation="kraus")

skq.quantum_info.DepolarizingChannel

Bases: QuantumChannel

Depolarizing noise (qubit) channel as Kraus representation. Special case of PauliNoiseChannel for p_x = p_y = p_z. :param p: Probability of depolarization.

Source code in src/skq/quantum_info/channel.py 
252
253
254
255
256
257
258
259
260
261
262
263
class DepolarizingChannel(QuantumChannel):
    """
    Depolarizing noise (qubit) channel as Kraus representation.
    Special case of PauliNoiseChannel for p_x = p_y = p_z.
    :param p: Probability of depolarization.
    """

    def __new__(cls, p: float):
        assert 0 <= p <= 1, "Depolarization probability must be in range [0...1]."
        cls.p = p
        kraus_ops = np.array([np.sqrt(1 - p) * I(), np.sqrt(p / 3) * X(), np.sqrt(p / 3) * Y(), np.sqrt(p / 3) * Z()])
        return super().__new__(cls, kraus_ops, representation="kraus")

skq.quantum_info.PauliNoiseChannel

Bases: QuantumChannel

Pauli noise channel as Kraus representation. :param p_x: Probability of X error. :param p_y: Probability of Y error. :param p_z: Probability of Z error.

Source code in src/skq/quantum_info/channel.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
class PauliNoiseChannel(QuantumChannel):
    """
    Pauli noise channel as Kraus representation.
    :param p_x: Probability of X error.
    :param p_y: Probability of Y error.
    :param p_z: Probability of Z error.
    """

    def __new__(cls, p_x: float, p_y: float, p_z: float):
        # Ensure the probabilities are within [0, 1] and the total does not exceed 1
        assert 0 <= p_x <= 1, "Probability p_x must be in range [0...1]."
        assert 0 <= p_y <= 1, "Probability p_y must be in range [0...1]."
        assert 0 <= p_z <= 1, "Probability p_z must be in range [0...1]."
        total_p = p_x + p_y + p_z
        assert total_p <= 1, "The sum of probabilities p_x, p_y, and p_z must not exceed 1."

        p_i = 1 - total_p
        kraus_ops = np.array([np.sqrt(p_i) * I(), np.sqrt(p_x) * X(), np.sqrt(p_y) * Y(), np.sqrt(p_z) * Z()])
        return super().__new__(cls, kraus_ops, representation="kraus")

skq.quantum_info.CompletelyDephasingChannel

Bases: QuantumChannel

Dephases a quantum (qubit) state in the computational basis. Initialized in the Kraus representation.

Source code in src/skq/quantum_info/channel.py 
287
288
289
290
291
292
293
294
295
class CompletelyDephasingChannel(QuantumChannel):
    """
    Dephases a quantum (qubit) state in the computational basis.
    Initialized in the Kraus representation.
    """

    def __new__(cls):
        kraus_ops = np.array([np.diag([1 if i == j else 0 for j in range(2)]) for i in range(2)])
        return super().__new__(cls, kraus_ops, representation="kraus")

skq.quantum_info.AmplitudeDampingChannel

Bases: QuantumChannel

Amplitude Damping Channel for a quantum (qubit) system. :param gamma: Damping parameter in range [0...1].

Source code in src/skq/quantum_info/channel.py 
298
299
300
301
302
303
304
305
306
307
308
class AmplitudeDampingChannel(QuantumChannel):
    """
    Amplitude Damping Channel for a quantum (qubit) system.
    :param gamma: Damping parameter in range [0...1].
    """

    def __new__(cls, gamma: float):
        assert 0 <= gamma <= 1, "Gamma must be in range [0...1]."
        cls.gamma = gamma
        kraus_ops = np.array([[[1, 0], [0, np.sqrt(1 - gamma)]], [[0, np.sqrt(gamma)], [0, 0]]])
        return super().__new__(cls, kraus_ops, representation="kraus")

skq.quantum_info.PhaseFlipChannel

Bases: QuantumChannel

Phase flip channel for a qubit system. Initialized in the Kraus representation. :param p: Probability of phase flip.

Source code in src/skq/quantum_info/channel.py 
311
312
313
314
315
316
317
318
319
320
321
322
class PhaseFlipChannel(QuantumChannel):
    """
    Phase flip channel for a qubit system.
    Initialized in the Kraus representation.
    :param p: Probability of phase flip.
    """

    def __new__(cls, p: float):
        assert 0 <= p <= 1, "Depolarization probability must be in range [0...1]."
        cls.p = p
        kraus_ops = np.array([np.sqrt(1 - p) * I(), np.sqrt(p) * Z()])
        return super().__new__(cls, kraus_ops, representation="kraus")

Hamiltonians

Hamiltonians represent the energy of quantum systems.

from skq.quantum_info import Hamiltonian, IsingHamiltonian

# Create an Ising model Hamiltonian for 3 qubits
H = IsingHamiltonian(num_qubits=3, J=1.0, h=0.5)
print(H.ground_state_energy())

skq.quantum_info.Hamiltonian

Bases: HermitianOperator

Class representing a Hamiltonian in quantum computing.

:param input_array: The input array representing the Hamiltonian. Will be converted to a complex numpy array. :param hbar: The reduced Planck constant. Default is 1.0 (natural units). If you want to use the actual physical value, set hbar to 1.0545718e-34.

Source code in src/skq/quantum_info/hamiltonian.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
class Hamiltonian(HermitianOperator):
    """
    Class representing a Hamiltonian in quantum computing.

    :param input_array: The input array representing the Hamiltonian. Will be converted to a complex numpy array.
    :param hbar: The reduced Planck constant. Default is 1.0 (natural units).
    If you want to use the actual physical value, set hbar to 1.0545718e-34.
    """

    def __new__(cls, input_array, hbar: float = 1.0):
        assert hbar > 0, "The reduced Planck constant must be greater than zero."
        obj = super().__new__(cls, input_array)
        obj.hbar = hbar
        return obj

    @property
    def num_qubits(self) -> int:
        """Number of qubits in the Hamiltonian."""
        return int(np.log2(self.shape[0]))

    def is_multi_qubit(self) -> bool:
        """Check if Hamiltonian involves multiple qubits."""
        return self.num_qubits > 1

    def time_evolution_operator(self, t: float) -> np.ndarray:
        """Time evolution operator U(t) = exp(-iHt/hbar)."""
        return scipy.linalg.expm(-1j * self * t / self.hbar)

    def ground_state_energy(self) -> float:
        """Ground state energy. i.e. the smallest eigenvalue."""
        eigenvalues = self.eigenvalues()
        return eigenvalues[0]

    def ground_state(self) -> np.ndarray:
        """The eigenvector corresponding to the smallest eigenvalue."""
        _, eigenvectors = np.linalg.eigh(self)
        return eigenvectors[:, 0]

    def convert_endianness(self) -> "Hamiltonian":
        """Hamiltonian from big-endian to little-endian and vice versa."""
        perm = np.argsort([int(bin(i)[2:].zfill(self.num_qubits)[::-1], 2) for i in range(2**self.num_qubits)])
        return self[np.ix_(perm, perm)]

    def add_noise(self, noise_strength: float, noise_operator: np.array) -> "Hamiltonian":
        """
        Add noise to the Hamiltonian.
        :param noise_strength: Strength of the noise. Generally in range [0.0001, 0.5].
        :param noise_operator: Operator representing the noise. For example Pauli matrices.
        :return: A new Hamiltonian with noise added.
        """
        assert noise_operator.shape == self.shape, "Noise operator must have the same dimensions as the Hamiltonian."
        return Hamiltonian(self + noise_strength * noise_operator, hbar=self.hbar)

    def to_qiskit(self) -> qiskit.quantum_info.Operator:
        """
        Convert the scikit-q Hamiltonian to a Qiskit Operator object.
        Qiskit using little endian convention, so we permute the order of the qubits.
        :return: Qiskit Operator object
        """
        return qiskit.quantum_info.Operator(self.convert_endianness())

    @staticmethod
    def from_qiskit(operator: qiskit.quantum_info.Operator) -> "Hamiltonian":
        """
        Create a scikit-q Hamiltonian object from a Qiskit Operator object.
        Qiskit using little endian convention, so we permute the order of the qubits.
        :param operator: Qiskit Operator object
        :return: Hamiltonian object
        """
        return Hamiltonian(operator.data).convert_endianness()

    def to_pennylane(self, wires: list[int] | int = None, **kwargs) -> "qml.Hamiltonian":
        """
        Convert the scikit-q Hamiltonian to a PennyLane Hamiltonian.
        :param wires: List of wires to apply the Hamiltonian to
        kwargs are passed to the PennyLane Hamiltonian constructor.
        :return: PennyLane Hamiltonian object
        """
        coefficients = [1.0]
        wires = wires if wires is not None else list(range(self.num_qubits))
        observables = [qml.Hermitian(self, wires=wires)]
        return qml.Hamiltonian(coefficients, observables, **kwargs)

    @staticmethod
    def from_pennylane(hamiltonian: qml.Hamiltonian) -> "Hamiltonian":
        """
        Convert a PennyLane Hamiltonian object to a scikit-q Hamiltonian object.
        :param hamiltonian: PennyLane Hamiltonian object
        :return: Hamiltonian object
        """
        assert len(hamiltonian.ops) == 1, "Only single-term Hamiltonians are supported."
        return Hamiltonian(hamiltonian.ops[0].matrix())

num_qubits property

Number of qubits in the Hamiltonian.

add_noise(noise_strength, noise_operator)

Add noise to the Hamiltonian. :param noise_strength: Strength of the noise. Generally in range [0.0001, 0.5]. :param noise_operator: Operator representing the noise. For example Pauli matrices. :return: A new Hamiltonian with noise added.

Source code in src/skq/quantum_info/hamiltonian.py 
53
54
55
56
57
58
59
60
61
def add_noise(self, noise_strength: float, noise_operator: np.array) -> "Hamiltonian":
    """
    Add noise to the Hamiltonian.
    :param noise_strength: Strength of the noise. Generally in range [0.0001, 0.5].
    :param noise_operator: Operator representing the noise. For example Pauli matrices.
    :return: A new Hamiltonian with noise added.
    """
    assert noise_operator.shape == self.shape, "Noise operator must have the same dimensions as the Hamiltonian."
    return Hamiltonian(self + noise_strength * noise_operator, hbar=self.hbar)

convert_endianness()

Hamiltonian from big-endian to little-endian and vice versa.

Source code in src/skq/quantum_info/hamiltonian.py 
48
49
50
51
def convert_endianness(self) -> "Hamiltonian":
    """Hamiltonian from big-endian to little-endian and vice versa."""
    perm = np.argsort([int(bin(i)[2:].zfill(self.num_qubits)[::-1], 2) for i in range(2**self.num_qubits)])
    return self[np.ix_(perm, perm)]

from_pennylane(hamiltonian) staticmethod

Convert a PennyLane Hamiltonian object to a scikit-q Hamiltonian object. :param hamiltonian: PennyLane Hamiltonian object :return: Hamiltonian object

Source code in src/skq/quantum_info/hamiltonian.py 
 93
 94
 95
 96
 97
 98
 99
100
101
@staticmethod
def from_pennylane(hamiltonian: qml.Hamiltonian) -> "Hamiltonian":
    """
    Convert a PennyLane Hamiltonian object to a scikit-q Hamiltonian object.
    :param hamiltonian: PennyLane Hamiltonian object
    :return: Hamiltonian object
    """
    assert len(hamiltonian.ops) == 1, "Only single-term Hamiltonians are supported."
    return Hamiltonian(hamiltonian.ops[0].matrix())

from_qiskit(operator) staticmethod

Create a scikit-q Hamiltonian object from a Qiskit Operator object. Qiskit using little endian convention, so we permute the order of the qubits. :param operator: Qiskit Operator object :return: Hamiltonian object

Source code in src/skq/quantum_info/hamiltonian.py 
71
72
73
74
75
76
77
78
79
@staticmethod
def from_qiskit(operator: qiskit.quantum_info.Operator) -> "Hamiltonian":
    """
    Create a scikit-q Hamiltonian object from a Qiskit Operator object.
    Qiskit using little endian convention, so we permute the order of the qubits.
    :param operator: Qiskit Operator object
    :return: Hamiltonian object
    """
    return Hamiltonian(operator.data).convert_endianness()

ground_state()

The eigenvector corresponding to the smallest eigenvalue.

Source code in src/skq/quantum_info/hamiltonian.py 
43
44
45
46
def ground_state(self) -> np.ndarray:
    """The eigenvector corresponding to the smallest eigenvalue."""
    _, eigenvectors = np.linalg.eigh(self)
    return eigenvectors[:, 0]

ground_state_energy()

Ground state energy. i.e. the smallest eigenvalue.

Source code in src/skq/quantum_info/hamiltonian.py 
38
39
40
41
def ground_state_energy(self) -> float:
    """Ground state energy. i.e. the smallest eigenvalue."""
    eigenvalues = self.eigenvalues()
    return eigenvalues[0]

is_multi_qubit()

Check if Hamiltonian involves multiple qubits.

Source code in src/skq/quantum_info/hamiltonian.py 
30
31
32
def is_multi_qubit(self) -> bool:
    """Check if Hamiltonian involves multiple qubits."""
    return self.num_qubits > 1

time_evolution_operator(t)

Time evolution operator U(t) = exp(-iHt/hbar).

Source code in src/skq/quantum_info/hamiltonian.py 
34
35
36
def time_evolution_operator(self, t: float) -> np.ndarray:
    """Time evolution operator U(t) = exp(-iHt/hbar)."""
    return scipy.linalg.expm(-1j * self * t / self.hbar)

to_pennylane(wires=None, **kwargs)

Convert the scikit-q Hamiltonian to a PennyLane Hamiltonian. :param wires: List of wires to apply the Hamiltonian to kwargs are passed to the PennyLane Hamiltonian constructor. :return: PennyLane Hamiltonian object

Source code in src/skq/quantum_info/hamiltonian.py 
81
82
83
84
85
86
87
88
89
90
91
def to_pennylane(self, wires: list[int] | int = None, **kwargs) -> "qml.Hamiltonian":
    """
    Convert the scikit-q Hamiltonian to a PennyLane Hamiltonian.
    :param wires: List of wires to apply the Hamiltonian to
    kwargs are passed to the PennyLane Hamiltonian constructor.
    :return: PennyLane Hamiltonian object
    """
    coefficients = [1.0]
    wires = wires if wires is not None else list(range(self.num_qubits))
    observables = [qml.Hermitian(self, wires=wires)]
    return qml.Hamiltonian(coefficients, observables, **kwargs)

to_qiskit()

Convert the scikit-q Hamiltonian to a Qiskit Operator object. Qiskit using little endian convention, so we permute the order of the qubits. :return: Qiskit Operator object

Source code in src/skq/quantum_info/hamiltonian.py 
63
64
65
66
67
68
69
def to_qiskit(self) -> qiskit.quantum_info.Operator:
    """
    Convert the scikit-q Hamiltonian to a Qiskit Operator object.
    Qiskit using little endian convention, so we permute the order of the qubits.
    :return: Qiskit Operator object
    """
    return qiskit.quantum_info.Operator(self.convert_endianness())

Predefined Hamiltonians

skq.quantum_info.IsingHamiltonian

Bases: Hamiltonian

Hamiltonian for the Ising model. :param num_qubits: Number of qubits in the system. :param J: Interaction strength between qubits. :param h: Transverse field strength. :param hbar: The reduced Planck constant. Default is 1.0 (natural units).

Source code in src/skq/quantum_info/hamiltonian.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
class IsingHamiltonian(Hamiltonian):
    """
    Hamiltonian for the Ising model.
    :param num_qubits: Number of qubits in the system.
    :param J: Interaction strength between qubits.
    :param h: Transverse field strength.
    :param hbar: The reduced Planck constant. Default is 1.0 (natural units).
    """

    def __new__(cls, num_qubits: int, J: float, h: float, hbar: float = 1.0):
        size = 2**num_qubits
        H = np.zeros((size, size), dtype=complex)
        # Interaction term: -J * sum(Z_i Z_{i+1})
        for i in range(num_qubits - 1):
            term = np.eye(1)
            for j in range(num_qubits):
                if j == i or j == i + 1:
                    term = np.kron(term, Z())
                else:
                    term = np.kron(term, I())
            H -= J * term

        # Transverse field term: -h * sum(X_i)
        for i in range(num_qubits):
            term = np.eye(1)
            for j in range(num_qubits):
                if j == i:
                    term = np.kron(term, X())
                else:
                    term = np.kron(term, I())
            H -= h * term

        return super().__new__(cls, H, hbar)

skq.quantum_info.HeisenbergHamiltonian

Bases: Hamiltonian

Hamiltonian for the Heisenberg model. :param num_qubits: Number of qubits in the system. :param J: Interaction strength between qubits. :param hbar: The reduced Planck constant. Default is 1.0 (natural units).

Source code in src/skq/quantum_info/hamiltonian.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
class HeisenbergHamiltonian(Hamiltonian):
    """
    Hamiltonian for the Heisenberg model.
    :param num_qubits: Number of qubits in the system.
    :param J: Interaction strength between qubits.
    :param hbar: The reduced Planck constant. Default is 1.0 (natural units).
    """

    def __new__(cls, num_qubits: int, J: float, hbar: float = 1.0):
        size = 2**num_qubits
        H = np.zeros((size, size), dtype=complex)

        # Interaction term: J * sum(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1})
        for i in range(num_qubits - 1):
            X_i = np.array([[0, 1], [1, 0]])
            Y_i = np.array([[0, -1j], [1j, 0]])
            Z_i = np.array([[1, 0], [0, -1]])
            H += J * (np.kron(np.kron(np.eye(2**i), np.kron(X_i, X_i)), np.eye(2 ** (num_qubits - i - 2))) + np.kron(np.kron(np.eye(2**i), np.kron(Y_i, Y_i)), np.eye(2 ** (num_qubits - i - 2))) + np.kron(np.kron(np.eye(2**i), np.kron(Z_i, Z_i)), np.eye(2 ** (num_qubits - i - 2))))

        return super().__new__(cls, H, hbar)

Hadamard Matrices

Hadamard matrices are useful in quantum information theory and quantum algorithms.

from skq.quantum_info import generate_hadamard_matrix

# Generate a Hadamard matrix of order 4
H = generate_hadamard_matrix(4)

skq.quantum_info.HadamardMatrix

Bases: Operator

Hadamard matrix representation.

Source code in src/skq/quantum_info/hadamard.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
class HadamardMatrix(Operator):
    """Hadamard matrix representation."""

    def __new__(cls, input_array: np.array):
        obj = super().__new__(cls, input_array)
        obj = np.real(obj).astype(int)
        assert obj.is_binary(), "Hadamard matrix must have entries of +1 or -1."
        assert obj.is_orthogonal(), "Hadamard matrix must have orthogonal rows and columns."
        assert obj.is_hadamard_order(), "The order of a Hadamard matrix must be 1, 2, or a multiple of 4."
        return obj

    def is_binary(self) -> bool:
        """All elements are +1 or -1."""
        return np.all(np.isin(self, [1, -1]))

    def is_orthogonal(self) -> bool:
        """Rows and columns are orthogonal."""
        n = self.shape[0]
        return np.allclose(self @ self.T, n * np.eye(n))

    def is_hadamard_order(self) -> bool:
        """Order of the Hadamard matrix is 1, 2, or a multiple of 4."""
        n = self.shape[0]
        return n == 1 or n == 2 or (n % 4 == 0)

    def determinant(self) -> float:
        """Determinant of the Hadamard matrix."""
        return np.linalg.det(self)

    def equivalence(self, other: "HadamardMatrix") -> bool:
        """Hadamard matrices are equivalent if row/column permutations and sign flips are equal."""
        permuted_matrices = self.permutations_and_sign_flips()
        return any(np.array_equal(matrix, other) for matrix in permuted_matrices)

    def spectral_norm(self) -> float:
        """Spectral norm of the Hadamard matrix.
        :return: largest singular value (i.e. spectral norm) of the Hadamard matrix.
        """
        return np.linalg.norm(self, ord=2)

    def permutations_and_sign_flips(self) -> list[np.array]:
        """Generate all matrix permutations by permuting rows/columns and flipping signs."""
        n = self.shape[0]
        matrices = []
        for row_perm in permutations(range(n)):
            for col_perm in permutations(range(n)):
                permuted_matrix = self[row_perm, :][:, col_perm]
                for signs in product([-1, 1], repeat=n):
                    sign_matrix = np.diag(signs)
                    matrices.append(sign_matrix @ permuted_matrix @ sign_matrix)
        return matrices

determinant()

Determinant of the Hadamard matrix.

Source code in src/skq/quantum_info/hadamard.py 
32
33
34
def determinant(self) -> float:
    """Determinant of the Hadamard matrix."""
    return np.linalg.det(self)

equivalence(other)

Hadamard matrices are equivalent if row/column permutations and sign flips are equal.

Source code in src/skq/quantum_info/hadamard.py 
36
37
38
39
def equivalence(self, other: "HadamardMatrix") -> bool:
    """Hadamard matrices are equivalent if row/column permutations and sign flips are equal."""
    permuted_matrices = self.permutations_and_sign_flips()
    return any(np.array_equal(matrix, other) for matrix in permuted_matrices)

is_binary()

All elements are +1 or -1.

Source code in src/skq/quantum_info/hadamard.py 
18
19
20
def is_binary(self) -> bool:
    """All elements are +1 or -1."""
    return np.all(np.isin(self, [1, -1]))

is_hadamard_order()

Order of the Hadamard matrix is 1, 2, or a multiple of 4.

Source code in src/skq/quantum_info/hadamard.py 
27
28
29
30
def is_hadamard_order(self) -> bool:
    """Order of the Hadamard matrix is 1, 2, or a multiple of 4."""
    n = self.shape[0]
    return n == 1 or n == 2 or (n % 4 == 0)

is_orthogonal()

Rows and columns are orthogonal.

Source code in src/skq/quantum_info/hadamard.py 
22
23
24
25
def is_orthogonal(self) -> bool:
    """Rows and columns are orthogonal."""
    n = self.shape[0]
    return np.allclose(self @ self.T, n * np.eye(n))

permutations_and_sign_flips()

Generate all matrix permutations by permuting rows/columns and flipping signs.

Source code in src/skq/quantum_info/hadamard.py 
47
48
49
50
51
52
53
54
55
56
57
def permutations_and_sign_flips(self) -> list[np.array]:
    """Generate all matrix permutations by permuting rows/columns and flipping signs."""
    n = self.shape[0]
    matrices = []
    for row_perm in permutations(range(n)):
        for col_perm in permutations(range(n)):
            permuted_matrix = self[row_perm, :][:, col_perm]
            for signs in product([-1, 1], repeat=n):
                sign_matrix = np.diag(signs)
                matrices.append(sign_matrix @ permuted_matrix @ sign_matrix)
    return matrices

spectral_norm()

Spectral norm of the Hadamard matrix. :return: largest singular value (i.e. spectral norm) of the Hadamard matrix.

Source code in src/skq/quantum_info/hadamard.py 
41
42
43
44
45
def spectral_norm(self) -> float:
    """Spectral norm of the Hadamard matrix.
    :return: largest singular value (i.e. spectral norm) of the Hadamard matrix.
    """
    return np.linalg.norm(self, ord=2)

skq.quantum_info.generate_hadamard_matrix(order)

Hadamard matrix of the given order using Sylvester's construction method. :param order: The order of the Hadamard matrix. Must be a power of 2. :return: HadamardMatrix of the specified order.

Source code in src/skq/quantum_info/hadamard.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def generate_hadamard_matrix(order: int) -> HadamardMatrix:
    """
    Hadamard matrix of the given order using Sylvester's construction method.
    :param order: The order of the Hadamard matrix. Must be a power of 2.
    :return: HadamardMatrix of the specified order.
    """
    assert order >= 1 and (order & (order - 1)) == 0, "Order must be a power of 2."

    # Recurse block matrices to construct Hadamard matrix
    if order == 2:
        matrix = np.array([[1, 1], [1, -1]])
    else:
        H_prev = generate_hadamard_matrix(order // 2)
        matrix = np.block([[H_prev, H_prev], [H_prev, -H_prev]])
    return HadamardMatrix(matrix)

Superoperators

Superoperators are linear maps that transform operators to operators.

skq.quantum_info.SuperOperator

Bases: ndarray

Base class for quantum superoperators. For example quantum channels.

Source code in src/skq/quantum_info/superoperator.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
class SuperOperator(np.ndarray):
    """Base class for quantum superoperators. For example quantum channels."""

    def __new__(cls, input_array: np.array):
        arr = np.array(input_array, dtype=complex)
        obj = arr.view(cls)
        assert obj.is_at_least_nxn(n=1), "Superoperator must be at least a 1x1 matrix."
        assert obj.is_power_of_n_shape(n=2), "Superoperator must have a power of 2 shape."
        return obj

    def is_at_least_nxn(self, n: int) -> bool:
        """Superoperator is at least an n x n matrix."""
        return self.shape[0] >= n and self.shape[1] >= n

    def is_power_of_n_shape(self, n: int) -> bool:
        """
        Superoperator shape is a power of n.
        :param n: Number to check for power of n shape.
        """

        def _is_power_of_n(x, n):
            if x < 1:
                return False
            while x % n == 0:
                x //= n
            return x == 1

        return _is_power_of_n(self.shape[0], n) and _is_power_of_n(self.shape[1], n)

    def conjugate_transpose(self) -> np.ndarray:
        """
        Conjugate transpose (i.e. Hermitian adjoint or 'dagger operation') of the superoperator.
        1. Transpose the matrix
        2. Take the complex conjugate of each element (Flip the sign of the imaginary part)
        """
        return self.T.conj()

conjugate_transpose()

Conjugate transpose (i.e. Hermitian adjoint or 'dagger operation') of the superoperator. 1. Transpose the matrix 2. Take the complex conjugate of each element (Flip the sign of the imaginary part)

Source code in src/skq/quantum_info/superoperator.py 
33
34
35
36
37
38
39
def conjugate_transpose(self) -> np.ndarray:
    """
    Conjugate transpose (i.e. Hermitian adjoint or 'dagger operation') of the superoperator.
    1. Transpose the matrix
    2. Take the complex conjugate of each element (Flip the sign of the imaginary part)
    """
    return self.T.conj()

is_at_least_nxn(n)

Superoperator is at least an n x n matrix.

Source code in src/skq/quantum_info/superoperator.py 
14
15
16
def is_at_least_nxn(self, n: int) -> bool:
    """Superoperator is at least an n x n matrix."""
    return self.shape[0] >= n and self.shape[1] >= n

is_power_of_n_shape(n)

Superoperator shape is a power of n. :param n: Number to check for power of n shape.

Source code in src/skq/quantum_info/superoperator.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def is_power_of_n_shape(self, n: int) -> bool:
    """
    Superoperator shape is a power of n.
    :param n: Number to check for power of n shape.
    """

    def _is_power_of_n(x, n):
        if x < 1:
            return False
        while x % n == 0:
            x //= n
        return x == 1

    return _is_power_of_n(self.shape[0], n) and _is_power_of_n(self.shape[1], n)