Physical Process GST#
This tutorial shows how to use the InterpolatedDenseOp and InterpolatedOpFactory to create quick-to-evaluate operations by interpolating between the discrete points at quick a more computationally-intensive process is performed. Often the computationally intensive process simulates the physics of a qubit gate, and would not practially work as a custom model operation because of the time required to evaluate it.
In order to turn such physical processes into gates, you should implement a custom PhysicalProcess object and then use the InterpolatedDenseOp or InterpolatedOpFactory class to interpolate the values of the custom process on a set of pre-defined points. All the physics simulation is then done at the time of creating the interpolated operation (or factory), after which the object can be saved for later use. An InterpolatedDenseOp or InterpolatedOpFactory object can be evaluated at any parameter-space point within the ranges over which the initial interpolation was performed.
All of this functionality is currently provided within the pygsti.extras.interpygate sub-package. This tutorial demonstrates how to setup a custom physical process and create an interpolated gate and factory object from it.
We’ll begin by some standard imports and by importing the interpygate sub-package. We get a MPI communicator if we can, as usually the physical simulation is performed using multiple processors.
import numpy as np
from scipy.linalg import expm
import pygsti
import pygsti.extras.interpygate as interp
try:
# See if MPI is available.
# For robustness, force the Libfabric/OFI provider to fallback to standard `sockets`.
# This prevents a hard C-level process abort on some platforms. One should remove
# this environment-variable assignment unless actually needed on your machine.
import os
os.environ['FI_PROVIDER'] = 'sockets'
from mpi4py import MPI
comm = MPI.COMM_WORLD
except ImportError:
comm = None
Defining a physical process#
We create a physical process simulator by deriving from the PhysicalProcess class and implementing its create_process_matrix function. This is the computationally intensive method that generates a process matrix based on some set of parameters. Every physical process has a fixed number of parameters that define the space that will be interpolated over. The generated process matrix is expected to be in whatever basis the ultimate Model operations will be in - usually the Pauli-product basis specified by "pp" - and have a fixed shape. This shape, given by process_shape below, is almost always a square matrix of dimension \(4^n\) where \(n\) is the number of qubits. Specifying an auxiliary information shape (aux_shape below) and implementing the create_aux_info will allow additional (floating point) values that describe the process to be interpolated.
Below we create a physical process that evolves a quantum state for some time (also a parameter) using a parameterized Lindbladian. Process tomography is used to construct a process matrix from the state evolution. The process has 6 parameters.
class ExampleProcess(interp.PhysicalProcess):
def __init__(self):
self.Hx = np.array([[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, -1],
[0, 0, 1, 0]], dtype='float')
self.Hy = np.array([[0, 0, 0, 0],
[0, 0, 0, 1],
[0, 0, 0, 0],
[0, -1, 0, 0]], dtype='float')
self.Hz = np.array([[0, 0, 0, 0],
[0, 0, -1, 0],
[0, 1, 0, 0],
[0, 0, 0, 0]], dtype='float')
self.dephasing_generator = np.diag([0, -1, -1, 0])
self.decoherence_generator = np.diag([0, -1, -1, -1])
num_params = 6 # omega (0), phase (1), detuning (2), dephasing (3), decoherence (4), time (5)
process_shape = (4, 4)
super().__init__(num_params, process_shape,
aux_shape=()) # our auxiliary information is a single float (None means no info)
def advance(self, state, v):
""" Evolves `state` in time """
state = np.array(state, dtype='complex')
omega, phase, detuning, dephasing, decoherence, t = v #Here are all our parameters
H = (omega * np.cos(phase) * self.Hx + omega * np.sin(phase) * self.Hy + detuning * self.Hz)
L = dephasing * self.dephasing_generator + decoherence * self.decoherence_generator
process = pygsti.tools.change_basis(expm((H + L) * t),'pp', 'col')
dim = state.size
vec_density_in = np.outer(state, state.conj()).ravel(order='F')
vec_density_out = process @ vec_density_in
state = vec_density_out.reshape((dim, dim), order='F')
return state
def create_process_matrix(self, v, comm=None):
def state_to_process_mxs(state):
return self.advance(state, v)
processes = interp.run_process_tomography(state_to_process_mxs, n_qubits=1,
basis='pp', comm=comm) # returns None on all but root processor
return np.array(processes) if (processes is not None) else None
def create_aux_info(self, v, comm=None):
omega, phase, detuning, dephasing, decoherence, t = v
return t*omega # matches aux_shape=() above
We can call create_process_matrix to generate a process matrix at a given set of parameters. Below we compute the ideal “target” operation by choosing the parameters corresponding to no errors.
example_process = ExampleProcess()
target_mx = example_process.create_process_matrix(np.array([1.0, 0.0, 0.0, 0.0, 0.0, np.pi/2]), comm=comm)
target_op = pygsti.modelmembers.operations.StaticArbitraryOp(target_mx)
print(target_op)
Making things more efficient#
We note that since our physical process is just an evolution in time, process matrices corresponding to different values of (just) the time parameter are especially easy to compute - a single evolution could compute, in one shot, the process matrices for an entire range of times.
The PhysicalProcess class contains support for such “easy-to-compute” parameters via the num_params_evaluated_as_group argument to its constructor. This argument defaults to 0, and specifies how many of the parameters, starting with the last one and working backward, should be evaluated within the same function call. If num_params_evaluated_as_group is set higher than 0, the derived class must implement the create_process_matrices and (optionally) create_aux_infos methods instead of create_process_matrix and create_aux_info. These methods take an additional grouped_v argument that contains arrays of values for the final num_params_evaluated_as_group parameters, and are expected return arrays of process matrices with corresponding shape (i.e., there is a leading index in the retured values for each “grouped” parameter).
We demonstrate this more complex usage below, where values for our final time argument are handled all at once.
class ExampleProcess_GroupTime(interp.PhysicalProcess):
def __init__(self):
self.Hx = np.array([[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, -1],
[0, 0, 1, 0]], dtype='float')
self.Hy = np.array([[0, 0, 0, 0],
[0, 0, 0, 1],
[0, 0, 0, 0],
[0, -1, 0, 0]], dtype='float')
self.Hz = np.array([[0, 0, 0, 0],
[0, 0, -1, 0],
[0, 1, 0, 0],
[0, 0, 0, 0]], dtype='float')
self.dephasing_generator = np.diag([0, -1, -1, 0])
self.decoherence_generator = np.diag([0, -1, -1, -1])
num_params = 6 # omega (0), phase (1), detuning (2), dephasing (3), decoherence (4), time (5)
process_shape = (4, 4)
super().__init__(num_params, process_shape,
aux_shape=(), # a single float
num_params_evaluated_as_group=1) # time values can be evaluated all at once
def advance(self, state, v, times):
state = np.array(state, dtype='complex')
omega, phase, detuning, dephasing, decoherence = v
H = (omega * np.cos(phase) * self.Hx + omega * np.sin(phase) * self.Hy + detuning * self.Hz)
L = dephasing * self.dephasing_generator + decoherence * self.decoherence_generator
processes = [pygsti.tools.change_basis(expm((H + L) * t),'pp', 'col') for t in times]
vec_density_in = np.outer(state, state.conj()).ravel(order='F')
dim = state.size
states = []
for process in processes:
vec_density_out = process @ vec_density_in
state = vec_density_out.reshape((dim, dim), order='F')
states.append(state)
return states
def create_process_matrices(self, v, grouped_v, comm=None):
assert(len(grouped_v) == 1) # we expect a single "grouped" parameter
times = grouped_v[0]
def state_to_process_mxs(state):
return self.advance(state, v, times)
processes = interp.run_process_tomography(state_to_process_mxs, n_qubits=1,
basis='pp', time_dependent=True, comm=comm)
return np.array(processes) if (processes is not None) else None
def create_aux_infos(self, v, grouped_v, comm=None):
omega, phase, detuning, dephasing, decoherence = v
times = grouped_v[0]
return np.array([t*omega for t in times], 'd')
We can similarly create a target operation from this physical process, but now we must specify a list of times.
example_process = ExampleProcess_GroupTime()
target_mx = example_process.create_process_matrices(np.array([1.0, 0.0, 0.0, 0.0, 0.0]), [[np.pi/2]], comm=comm)[0]
target_op = pygsti.modelmembers.operations.StaticArbitraryOp(target_mx)
print(target_op)
Creating an interpolated operation (gate)#
Now that we’ve done the hard work of creating the physical process, it’s easy to create an operator that evaluates the physical process on a grid of points and interpolates between them. The resulting InterpolatedDenseOp can be evaluated (i.e. from_vector can be invoked) at any point within the range being interpolated.
The parameters of the resulting InterpolatedDenseOp are the same as those of the underlying PhysicalProcess, and ranges are specified using either a (min, max, num_points) tuple or an array of values. Below we use only 2 points in most directions so it doesn’t take too long to run.
Creating the object also requires a target operation, for which we use target_op as defined above. This is required because internally it is the error generator rather than the process matrix itself that is interpolated. The target operation can be parameterized by any contiguous subset of the physical process’s parameters, starting with the first one. In our example, target_op is a StaticArbitraryOp and so takes 0 parameters. This should be interpreted as the “first 0 parameters of our example process”.
param_ranges = ([(0.9, 1.1, 2), # omega
(-.1, .1, 2), # phase
(-.1, .1, 2), # detuning
(0, 0.1, 2), # dephasing
(0, 0.1, 2), # decoherence
np.linspace(np.pi / 2, np.pi / 2 + .5, 10) # time
])
interp_op = interp.InterpolatedDenseOp.create_by_interpolating_physical_process(
target_op, example_process, param_ranges, comm=comm)
The created interp_op can then be evaluated (quickly) at points in parameter space.
interp_op.from_vector([1.1, 0.01, 0.01, 0.055, 0.055, 1.59])
interp_op.to_dense()
The auxiliary information can be retrieved from any interpolated operator via its aux_info attribute.
interp_op.aux_info
Creating an interpolated operation factory#
Operation factories in pyGSTi take “arguments” provided by in-circuit labels and produce operations. For example, the value of the rotation angle might be specified over a continuous interval by the algorithm being run, rather than being noise parameter that is fit to data when a model is optimized (e.g. in GST).
The InterpolatedOpFactory object interpolates a physical process, similar to InterpolatedDenseOp, but allows the user to divide the parameters of the physical process into factory arguments and operation parameters. The first group is meant to range over different intended (target) operations, and the latter group is meant to be unkonwn quantities determined by fitting a model to data. To create an InterpolatedOpFactory, we must first create a custom factory class that creates the target operation corresponding to a given set of arguments. As in the case of InterpolatedDenseOp, the target operations can be parameterized by any contiguous subset of the factory’s parameters, starting with the first one.
We choose to make a factory that takes as arguments the time and omega physical process parameters.
class TargetOpFactory(pygsti.modelmembers.operations.OpFactory):
def __init__(self):
self.process = ExampleProcess_GroupTime()
pygsti.modelmembers.operations.OpFactory.__init__(self, state_space=1, evotype="densitymx")
def create_object(self, args=None, sslbls=None):
assert(sslbls is None)
assert(len(args) == 2) # t (time), omega
t, omega = args
mx = self.process.create_process_matrices(np.array([omega, 0.0, 0.0, 0.0, 0.0]), [[t]], comm=None)[0]
#mx = self.process.create_process_matrix(np.array([omega, 0.0, 0.0, 0.0, 0.0, t]), comm=None) # Use this if using our initial ExampleProcess above.
return pygsti.modelmembers.operations.StaticArbitraryOp(mx)
We can then create an InterpolatedOpFactory similarly to how we created an InterpolatedDenseOp except now we separately specify factory argument and optimization parameter ranges, and specify which of the underlying physical process’s parameters are turned into factory arguments (arg_indices below).
arg_ranges = [np.linspace(np.pi / 2, np.pi / 2 + .5, 10), # time
(0.9, 1.1, 2) # omega
]
param_ranges = [(-.1, .1, 2), # phase
(-.1, .1, 2), # detuning
(0, 0.1, 2), # dephasing
(0, 0.1, 2) # decoherence
]
arg_indices = [5, 0] #indices for time and omega within ExampleProcess_GroupTime's parameters
opfactory = interp.InterpolatedOpFactory.create_by_interpolating_physical_process(
TargetOpFactory(), example_process, arg_ranges, param_ranges, arg_indices, comm=comm)
Note that the factory has only 4 parameters (whereas the physical process and the interpolated operator we made above have 6). This is because 2 of the physical process parameters have been turned into factory arguments.
print(opfactory.num_params)
print(interp_op.num_params)
print(example_process.num_params)
We can use the factory to create an InterpolatedDenseOp operation at a given time and omega pair:
opfactory.from_vector(np.array([0.01, 0.01, 0.055, 0.055]))
op = opfactory.create_op((1.59, 1.1))
op.to_dense()
op.aux_info