How-to use advanced system modelling functionality¶
The systems module contains tools for building and representing models of systems on tensor
product spaces. The high level usage of this module is demonstrated in the Building and solving
models of quantum systems tutorial, in which pre-built models are used
to simulate transmon systems. This userguide walks through how to use the underlying operator
construction tools to build new models, or to modify models in non-trivial ways: e.g. restriction to
a subspace of interest.
This user-guide walks through the following tasks:
How-to define operators acting on tensor product spaces.
How-to define an operator acting on a subspace.
How-to define an operator acting on the logical subspace of a multi-transmon model.
How-to restrict an operator to a low energy subspace of a 3-transmon model.
1. How-to build operators acting on tensor product spaces¶
Here we will walk-through the construction of operators acting on a tri-partite system. First, we define three 2-dimensional subsystems:
from qiskit_dynamics.systems import Subsystem
Q1 = Subsystem("Q1", dim=2)
Q2 = Subsystem("Q2", dim=2)
Q3 = Subsystem("Q3", dim=2)
Define several operators acting on the individual subsystems.
from qiskit_dynamics.systems import I, X, Y, Z, N
# X acting on Q1
X1 = X(Q1)
# Y acting on Q2
Y2 = Y(Q2)
# Z acting on Q3
Z3 = Z(Q3)
print(X1)
X(Q1)
Note that while these operators are defined only on individual subsystems, if they are used in the
context of multi-subsystem models, the operators implicitly act as the identity on the subsystems
lying outside the operator’s definition. E.g. X(Q1) is defined on Q1 only, but when thought
of as an operator on \(Q_1 \otimes Q_2\), it represents the operator \(X \otimes I\).
Once an operator is construct, the matrix() method returns its matrix form. By default, the
matrix is constructed only on the subsystems the operator explicitly acts on:
X1.matrix()
array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]])
To construct the matrix corresponding to X1 when viewed as an operator acting on the combined
tensor product space \(Q_1 \otimes Q_2\), explicitly pass the list of subsystems to the
matrix() method:
X1.matrix([Q1, Q2])
array([[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j]])
We can also reverse the ordering of the tensor product representation by reversing the order of the subsystem list:
X1.matrix([Q2, Q1])
array([[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]])
More complicated operators can be built through algebraic operations, e.g. addition:
X1 + Y2
X(Q1) + Y(Q2)
This new composite operator acts on the combined set of subsystems that X1 and Y2 act on:
(X1 + Y2).subsystems
[Subsystem(name=Q1, dim=2), Subsystem(name=Q2, dim=2)]
When calling the matrix() method, a matrix will be constructed on the tensor product system
in the above order.
(X1 + Y2).matrix()
array([[0.+0.j, 1.+0.j, 0.-1.j, 0.+0.j],
[1.+0.j, 0.+0.j, 0.+0.j, 0.-1.j],
[0.+1.j, 0.+0.j, 0.+0.j, 1.+0.j],
[0.+0.j, 0.+1.j, 1.+0.j, 0.+0.j]])
Similarly, we can build the matrix for X1 + Y2 when viewed as an operator on the tripartite
system \(Q_1 \otimes Q_2 \otimes Q_3\).
(X1 + Y2).matrix([Q1, Q2, Q3])
array([[0.+0.j, 1.+0.j, 0.-1.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[1.+0.j, 0.+0.j, 0.+0.j, 0.-1.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+1.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+1.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.-1.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.-1.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+1.j, 0.+0.j, 0.+0.j, 1.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+1.j, 1.+0.j, 0.+0.j]])
Matrix multiplication can also be performed:
X1 @ Z3
X(Q1) @ Z(Q3)
In the above case, as X1 acts on Q1 and Z3 acts on Q3, X1 @ Z3 represents the
operator \(X \otimes Z\) acting on the space \(Q_1 \otimes Q_3\).
Lastly, we can multiply and add scalars to operators. Scalars under addition are treated as multiples of the identity.
1 + 2 * X1
2 * X(Q1) + 1
2. How-to define an operator acting on a subspace¶
It is common in quantum information to define operators on subspaces, e.g. the computational subspace of a physical system. Here we walk through constructing an \(X\) operator on the first two levels of a \(4\)-dimensional system. In mathematical notation, we want to construct the operator \(X \oplus 0\), where both \(X\) and \(0\) are \(2 \times 2\) matrices.
First, define Subsystem instances representing both the \(2\)-dimensional subspace and
the full \(4\)-dimensional space.
# define subsystem for the subspace
C2 = Subsystem("C2", dim=2)
# define the higher dimensional space
C4 = Subsystem("C4", dim=4)
Next, define a ONBasis instance for the subspace of C4 representing how the standard
basis elements of C2 are mapped into C4. Here we will use the first two standard basis
elements of C4, representing the first two levels of C4.
from qiskit_dynamics.systems import ONBasis
import numpy as np
basis = ONBasis(
basis_vectors=np.eye(4, 2),
subsystems=[C4]
)
# view the basis vectors
basis.basis_vectors
array([[1., 0.],
[0., 1.],
[0., 0.],
[0., 0.]])
Using the basis, we explicitly construct the injection of C2 into C4 using the
SubsystemMapping class. This class represents a linear map between vector spaces
\(V \rightarrow W\) of the form \(v \mapsto Av\) for an operator \(A\). Acting on an
operator \(X\), a SubsystemMapping will perform the transformation
\(X \mapsto A X A^\dagger\).
Here, the operator \(A\) is defined as a matrix given by the first two basis elements of C4:
from qiskit_dynamics.systems import SubsystemMapping
injection = SubsystemMapping(
matrix=basis.basis_vectors,
in_subsystems=[C2],
out_subsystems=[C4]
)
\(X\) acting on the first two levels of C4 is then constructed as:
subspace_X = injection(X(C2))
subspace_X
MappedOperator(X(C2), [Subsystem(name=C2, dim=2)] ->
[Subsystem(name=C4, dim=4)])
Observe the desired matrix:
subspace_X.matrix()
array([[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])
3. How-to define an operator acting on the logical subspace of a multi-transmon model¶
In this section we work through a more advanced version of the previous example. Here, we consider the problem of constructing the operator “\(X\) acting on the computional subspace of the first qubit in a two-transmon system”. Mathematically, this means the matrix \(A(X \otimes I)A^\dagger\), where \(X\) and \(I\) are \(2 \times 2\) matrices, and \(A\) is the isometry mapping the two qubit computational subspace (the first 4 energy levels) into the two transmon physical space.
For this, we walk through the following steps:
Define subsystems for both the logical/computational spaces, and the physical spaces.
Construct the standard static Hamiltonian for a 2 transmon model, and compute the dressed basis (the basis of energy eigenstates).
Construct a basis for the computational subspace within the physical space.
Define the operator \(X\) acting on the logical qubit \(0\).
“Expand” this operator into the full physical space, creating the desired operator \(A(X \otimes I)A^\dagger\)
First, construct the Subsystem instances we will work with:
# logical subsystems
L0 = Subsystem("L0", dim=2)
L1 = Subsystem("L1", dim=2)
# physical subsystems
Q0 = Subsystem("Q0", dim=3)
Q1 = Subsystem("Q1", dim=3)
Define the 2 transmon Hamiltonian.
# define a standard Hamiltonian
H = (2 * np.pi * 5. * N(Q0) +(- 0.33) * np.pi * N(Q0) @ (N(Q0) + (-1 * I(Q0))) +
2 * np.pi * 5.5 * N(Q1) +(- 0.33) * np.pi * N(Q1) @ (N(Q1) + (-1 * I(Q1))) +
2 * np.pi * 0.002 * X(Q0) @ X(Q1))
Compute the dressed basis, and retrieve the ONBasis instance corresponding to the
computational states.
from qiskit_dynamics.systems import DressedBasis
# Get the dressed basis with an explicit tensor product ordering
dressed_basis = DressedBasis.from_hamiltonian(H, [Q0, Q1])
# retrieve the computational states
computational_states = dressed_basis.computational_states
Define the mapping of the logical space \(L_0 \otimes L_1\) into the computational subspace of the physical space \(Q_0 \otimes Q_1\) specified by the matrix of basis vectors for the computational subspace.
injection = SubsystemMapping(
matrix=computational_states.basis_vectors,
in_subsystems=[L0, L1],
out_subsystems=[Q0, Q1]
)
Finally, define X acting on L0, and inject it into the full two-transmon physical space
using injection. Note that as the injection acts on the combined \(L_0 \otimes L_1\) system,
X(L0) will be treated as X(L0) @ I(L1) when performing the injection (i.e. with implicit
identity on \(L_1\)).
op = X(L0)
injected_X0 = injection(op)
injected_X0
MappedOperator(X(L0), [Subsystem(name=L0, dim=2), Subsystem(name=L1, dim=2)] ->
[Subsystem(name=Q0, dim=3), Subsystem(name=Q1, dim=3)])
4. How-to restrict an operator to a low energy subspace of a 3-transmon model¶
Similarly to defining an operator on a subspace and expanding it into the full space, we may want to restrict on operator or model to a subspace. For example, restricting a model to a low energy subspace is a common technique to reduce the dimension of a model.
Here, we walk through the problem of restricting an operator to a low energy subspace of a 3 transmon system with the following steps:
Build the static Hamiltonian of a 3 transmon system.
Restrict it to a subspace with bounded energy.
Restrict the X operator acting on one of the transmons to the same subspace.
Define a 3 transmon static Hamiltonian:
# physical subsystems
Q0 = Subsystem("Q0", dim=3)
Q1 = Subsystem("Q1", dim=3)
Q2 = Subsystem("Q2", dim=3)
# define a standard Hamiltonian
H = (
2 * np.pi * 5. * N(Q0) +(- 0.33) * np.pi * N(Q0) @ (N(Q0) + (-1 * I(Q0))) +
2 * np.pi * 5.5 * N(Q1) +(- 0.33) * np.pi * N(Q1) @ (N(Q1) + (-1 * I(Q1))) +
2 * np.pi * 5.3 * N(Q2) +(- 0.33) * np.pi * N(Q2) @ (N(Q2) + (-1 * I(Q2))) +
2 * np.pi * 0.002 * X(Q0) @ X(Q1) +
2 * np.pi * 0.002 * X(Q1) @ X(Q2)
)
Construct the dressed basis and view eigenvalues.
from qiskit_dynamics.systems import DressedBasis
# Get the dressed basis
dressed_basis = DressedBasis.from_hamiltonian(H, [Q0, Q1, Q2])
dressed_basis.evals
array([-4.72070179e-06, 3.14158690e+01, 6.07583390e+01, 3.45576854e+01,
6.59733237e+01, 9.53158459e+01, 6.70414964e+01, 9.84573332e+01,
1.27799613e+02, 3.33007493e+01, 6.47166230e+01, 9.40590930e+01,
6.78589204e+01, 9.92745594e+01, 1.28617081e+02, 1.00343272e+02,
1.31759110e+02, 1.61101389e+02, 6.45282158e+01, 9.59440896e+01,
1.25286560e+02, 9.90853795e+01, 1.30501019e+02, 1.59843541e+02,
1.31570205e+02, 1.62986043e+02, 1.92328322e+02])
Construct a basis for a low energy subspace below a given cutoff.
low_energy_states = dressed_basis.low_energy_states(cutoff_energy=70.)
low_energy_states.evals
array([-4.72070179e-06, 3.14158690e+01, 6.07583390e+01, 3.45576854e+01,
6.59733237e+01, 6.70414964e+01, 3.33007493e+01, 6.47166230e+01,
6.78589204e+01, 6.45282158e+01])
Observe the standard basis labelling, and note that this energy cutoff happens to correspond to the subspace with at most 2-excitations in the full system.
low_energy_states.labels
[{'index': (0, 0, 0), 'eval': np.float64(-4.72070178672793e-06)},
{'index': (1, 0, 0), 'eval': np.float64(31.415869000063307)},
{'index': (2, 0, 0), 'eval': np.float64(60.75833903196221)},
{'index': (0, 1, 0), 'eval': np.float64(34.55768535500057)},
{'index': (1, 1, 0), 'eval': np.float64(65.97332367727418)},
{'index': (0, 2, 0), 'eval': np.float64(67.04149638508592)},
{'index': (0, 0, 1), 'eval': np.float64(33.300749290809065)},
{'index': (1, 0, 1), 'eval': np.float64(64.71662302225069)},
{'index': (0, 1, 1), 'eval': np.float64(67.85892038225253)},
{'index': (0, 0, 2), 'eval': np.float64(64.52821584309763)}]
Restrict the Hamiltonian to this low energy space. Note that we first need to define a
Subsystem instance representing this subspace in isolation.
LESpace = Subsystem("LES", dim=len(low_energy_states))
restriction = SubsystemMapping(
matrix=low_energy_states.basis_vectors_adj,
in_subsystems=[Q0, Q1, Q2],
out_subsystems=[LESpace]
)
low_energy_H = restriction(H)
Looking at the diagonal of low_energy_H, we can confirm that the entries are the eigenvalues
below the cutoff.
np.diag(low_energy_H.matrix())
array([-4.72070179e-06+0.00000000e+00j, 3.14158690e+01+0.00000000e+00j,
6.07583390e+01+0.00000000e+00j, 3.45576854e+01+0.00000000e+00j,
6.59733237e+01+0.00000000e+00j, 6.70414964e+01+0.00000000e+00j,
3.33007493e+01-7.88860905e-31j, 6.47166230e+01+0.00000000e+00j,
6.78589204e+01+0.00000000e+00j, 6.45282158e+01+0.00000000e+00j])
With this restriction mapping, we can also restrict other operators to this subspace, e.g. the
\(X\) operator acting on the physical Q1 system. Note that when the restriction map is
applied to X(Q1), the operator is interpreted as I(Q0) @ X(Q1) @ I(Q2), i.e. \(X\)
acting on Q1, and the identity on the remaining subsystems in the input space of restriction
that X(Q1) does not explicitly act on.
drive_op = restriction(X(Q1))
drive_op
MappedOperator(X(Q1), [Subsystem(name=Q0, dim=3), Subsystem(name=Q1, dim=3), Subsystem(name=Q2, dim=3)] ->
[Subsystem(name=LES, dim=10)])