# This code is part of Qiskit.
#
# (C) Copyright IBM 2021.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Constrained convex least-squares tomography fitter.
"""
from typing import Optional, Dict, Tuple, Union
import time
import numpy as np
from qiskit_experiments.library.tomography.basis import (
MeasurementBasis,
PreparationBasis,
)
from . import cvxpy_utils
from .cvxpy_utils import cvxpy
from . import lstsq_utils
from .fitter_data import _basis_dimensions
[docs]
@cvxpy_utils.requires_cvxpy
def cvxpy_linear_lstsq(
outcome_data: np.ndarray,
shot_data: np.ndarray,
measurement_data: np.ndarray,
preparation_data: np.ndarray,
measurement_basis: Optional[MeasurementBasis] = None,
preparation_basis: Optional[PreparationBasis] = None,
measurement_qubits: Optional[Tuple[int, ...]] = None,
preparation_qubits: Optional[Tuple[int, ...]] = None,
conditional_measurement_indices: Optional[Tuple[int, ...]] = None,
conditional_preparation_indices: Optional[Tuple[int, ...]] = None,
trace: Union[None, float, str] = "auto",
psd: bool = True,
trace_preserving: Union[None, bool, str] = "auto",
partial_trace: Optional[np.ndarray] = None,
weights: Optional[np.ndarray] = None,
**kwargs,
) -> Tuple[np.ndarray, Dict]:
r"""Constrained weighted linear least-squares tomography fitter.
Overview
This fitter reconstructs the maximum-likelihood estimate by using
``cvxpy`` to minimize the constrained least-squares negative log
likelihood function
.. math::
\hat{\rho}
&= -\mbox{argmin }\log\mathcal{L}{\rho} \\
&= \mbox{argmin }\sum_i w_i^2(\mbox{Tr}[E_j\rho] - \hat{p}_i)^2 \\
&= \mbox{argmin }\|W(Ax - y) \|_2^2
subject to
- *Positive-semidefinite* (``psd=True``): :math:`\rho \gg 0` is constrained
to be a positive-semidefinite matrix.
- *Trace* (``trace=t``): :math:`\mbox{Tr}(\rho) = t` is constrained to have
the specified trace.
- *Trace preserving* (``trace_preserving=True``): When performing process
tomography the Choi-state :math:`\rho` represents is constrained to be
trace preserving.
where
- :math:`A` is the matrix of measurement operators
:math:`A = \sum_i |i\rangle\!\langle\!\langle M_i|`
- :math:`y` is the vector of expectation value data for each projector
corresponding to estimates of :math:`b_i = Tr[M_i \cdot x]`.
- :math:`x` is the vectorized density matrix (or Choi-matrix) to be fitted
:math:`x = |\rho\rangle\\!\rangle`.
.. note:
Various solvers can be called in CVXPY using the `solver` keyword
argument. When ``psd=True`` the optimization problem is a case of a
*semidefinite program* (SDP) and requires a SDP compatible solver
for CVXPY. CVXPY includes an SDP compatible solver `SCS`` but it
is recommended to install the the open-source ``CVXOPT`` solver
or one of the supported commercial solvers. See the `CVXPY
documentation
<https://www.cvxpy.org/tutorial/advanced/index.html#solve-method-options>`_
for more information on solvers.
.. note::
Linear least-squares constructs the full basis matrix :math:`A` as a dense
numpy array so should not be used for than 5 or 6 qubits. For larger number
of qubits try the
:func:`~qiskit_experiments.library.tomography.fitters.linear_inversion`
fitter function.
Args:
outcome_data: measurement outcome frequency data.
shot_data: basis measurement total shot data.
measurement_data: measurement basis index data.
preparation_data: preparation basis index data.
measurement_basis: Optional, measurement matrix basis.
preparation_basis: Optional, preparation matrix basis.
measurement_qubits: Optional, the physical qubits that were measured.
If None they are assumed to be ``[0, ..., M-1]`` for M measured qubits.
preparation_qubits: Optional, the physical qubits that were prepared.
If None they are assumed to be ``[0, ..., N-1]`` for N prepared qubits.
conditional_measurement_indices: Optional, conditional measurement data
indices. If set this will return a list of fitted states conditioned
on a fixed basis measurement of these qubits.
conditional_preparation_indices: Optional, conditional preparation data
indices. If set this will return a list of fitted states conditioned
on a fixed basis preparation of these qubits.
trace: trace constraint for the fitted matrix. If "auto" this will be set
to 1 for QST or the input dimension for QST (default: "auto").
psd: If True rescale the eigenvalues of fitted matrix to be positive
semidefinite (default: True)
trace_preserving: Enforce the fitted matrix to be trace preserving when
fitting a Choi-matrix in quantum process
tomography. If "auto" this will be set to True for
QPT and False for QST (default: "auto").
partial_trace: Enforce conditional fitted Choi matrices to partial
trace to POVM matrices.
weights: Optional array of weights for least squares objective. Weights
should be the same shape as the outcome_data.
kwargs: kwargs for cvxpy solver.
Raises:
QiskitError: If CVXPy is not installed on the current system.
AnalysisError: If analysis fails.
Returns:
The fitted matrix rho that maximizes the least-squares likelihood function.
"""
t_start = time.time()
if measurement_basis and measurement_qubits is None:
measurement_qubits = tuple(range(measurement_data.shape[1]))
if preparation_basis and preparation_qubits is None:
preparation_qubits = tuple(range(preparation_data.shape[1]))
input_dims = _basis_dimensions(
basis=preparation_basis,
qubits=preparation_qubits,
conditional_indices=conditional_preparation_indices,
)
output_dims = _basis_dimensions(
basis=measurement_basis,
qubits=measurement_qubits,
conditional_indices=conditional_measurement_indices,
)
if trace_preserving == "auto" and preparation_data.shape[1] > 0:
trace_preserving = True
if trace == "auto" and output_dims is not None:
trace = np.prod(output_dims)
metadata = {
"fitter": "cvxpy_linear_lstsq",
"input_dims": input_dims,
"output_dims": output_dims,
}
# Conditional measurement indices
if conditional_measurement_indices:
cond_meas_data = measurement_data[:, np.array(conditional_measurement_indices, dtype=int)]
cond_meas_indices = np.unique(cond_meas_data, axis=0)
num_meas_cond = len(cond_meas_indices)
metadata["conditional_measurement_index"] = []
metadata["conditional_measurement_outcome"] = []
else:
num_meas_cond = 0
cond_meas_data = np.zeros((measurement_data.shape[0], 0), dtype=int)
cond_meas_indices = np.zeros((1, 0), dtype=int)
if conditional_preparation_indices:
cond_prep_data = preparation_data[:, np.array(conditional_preparation_indices, dtype=int)]
cond_prep_indices = np.unique(cond_prep_data, axis=0)
num_prep_cond = len(cond_meas_indices)
metadata["conditional_preparation_index"] = []
else:
num_prep_cond = 0
cond_prep_data = np.zeros((preparation_data.shape[0], 0), dtype=int)
cond_prep_indices = np.zeros((1, 0), dtype=int)
if outcome_data.shape[0] > 1:
metadata["conditional_circuit_outcome"] = []
fits = []
for cond_prep_idx in cond_prep_indices:
for cond_meas_idx in cond_meas_indices:
# Mask for specified conditional indices
cond_mask = np.all(cond_meas_data == cond_meas_idx, axis=1) & np.all(
cond_prep_data == cond_prep_idx, axis=1
)
if weights is None:
cond_weights = None
else:
cond_weights = weights[:, cond_mask]
basis_matrix, probability_data, probability_weights = lstsq_utils.lstsq_data(
outcome_data[:, cond_mask],
shot_data[cond_mask],
measurement_data[cond_mask],
preparation_data[cond_mask],
measurement_basis=measurement_basis,
preparation_basis=preparation_basis,
measurement_qubits=measurement_qubits,
preparation_qubits=preparation_qubits,
weights=cond_weights,
conditional_measurement_indices=conditional_measurement_indices,
conditional_preparation_indices=conditional_preparation_indices,
)
# Since CVXPY only works with real variables we must specify the real
# and imaginary parts of matrices separately: rho = rho_r + 1j * rho_i
num_circ_components, num_tomo_components, _ = probability_data.shape
dim = int(np.sqrt(basis_matrix.shape[1]))
# Generate list of conditional components for block diagonal matrix
# rho = sum_k |k><k| \otimes rho(k)
rhos_r = []
rhos_i = []
cons = []
for i in range(num_circ_components):
for j in range(num_tomo_components):
rho_r, rho_i, cons_i = cvxpy_utils.complex_matrix_variable(
dim, hermitian=True, psd=psd
)
rhos_r.append(rho_r)
rhos_i.append(rho_i)
cons.append(cons_i)
if num_circ_components > 1:
metadata["conditional_circuit_outcome"].append(i)
if num_meas_cond:
metadata["conditional_measurement_index"].append(tuple(cond_meas_idx))
metadata["conditional_measurement_outcome"].append(j)
if num_prep_cond:
metadata["conditional_preparation_index"].append(tuple(cond_prep_idx))
# Partial trace when fitting Choi-matrices for quantum process tomography.
# This applied to the sum of conditional components
# Note that this adds an implicitly
# trace preserving is a specific partial trace constraint ptrace(rho) = I
# Note: partial trace constraints implicitly define a trace constraint,
# so if a different trace constraint is specified it will be ignored
joint_cons = None
if partial_trace is not None:
for rho_r, rho_i, povm in zip(rhos_r, rhos_i, partial_trace):
joint_cons = cvxpy_utils.partial_trace_constaint(rho_r, rho_i, povm)
elif trace_preserving:
input_dim = np.prod(input_dims)
joint_cons = cvxpy_utils.trace_preserving_constaint(
rhos_r,
rhos_i,
input_dim=input_dim,
hermitian=True,
)
elif trace is not None:
joint_cons = cvxpy_utils.trace_constraint(
rhos_r, rhos_i, trace=trace, hermitian=True
)
# OBJECTIVE FUNCTION
# The function we wish to minimize is || arg ||_2 where
# arg = bm * vec(rho) - data
# Since we are working with real matrices in CVXPY we expand this as
# bm * vec(rho) = (bm_r + 1j * bm_i) * vec(rho_r + 1j * rho_i)
# = bm_r * vec(rho_r) - bm_i * vec(rho_i)
# + 1j * (bm_r * vec(rho_i) + bm_i * vec(rho_r))
# = bm_r * vec(rho_r) - bm_i * vec(rho_i)
# where we drop the imaginary part since the expectation value is real
# Construct block diagonal fit variable from conditional components
# Construct objective function
if probability_weights is not None:
probability_data = probability_weights * probability_data
bms_r = []
bms_i = []
for i in range(num_circ_components):
for j in range(num_tomo_components):
weighted_mat = probability_weights[i, j][:, None] * basis_matrix
bms_r.append(np.real(weighted_mat))
bms_i.append(np.imag(weighted_mat))
else:
bm_r = np.real(basis_matrix)
bm_i = np.imag(basis_matrix)
bms_r = [bm_r] * num_circ_components * num_tomo_components
bms_i = [bm_i] * num_circ_components * num_tomo_components
# Stack lstsq objective from sum of components
args = []
idx = 0
for i in range(num_circ_components):
for j in range(num_tomo_components):
model = bms_r[idx] @ cvxpy.vec(rhos_r[idx], order="F") - bms_i[idx] @ cvxpy.vec(
rhos_i[idx], order="F"
)
data = probability_data[i, j]
args.append(model - data)
idx += 1
# Combine all variables and constraints into a joint optimization problem
# if there is a joint constraint
if joint_cons:
args = [cvxpy.hstack(args)]
for cons_i in cons:
joint_cons += cons_i
cons = [joint_cons]
# Solve each component separately
metadata["cvxpy_solver"] = None
metadata["cvxpy_status"] = []
for arg, con in zip(args, cons):
# Optimization problem
obj = cvxpy.Minimize(cvxpy.norm(arg, p=2))
prob = cvxpy.Problem(obj, con)
# Solve SDP
cvxpy_utils.solve_iteratively(prob, 5000, **kwargs)
# Return optimal values and problem metadata
metadata["cvxpy_solver"] = prob.solver_stats.solver_name
metadata["cvxpy_status"].append(prob.status)
fits += [rho_r.value + 1j * rho_i.value for rho_r, rho_i in zip(rhos_r, rhos_i)]
# Add additional metadata
if psd:
metadata["psd_constraint"] = True
if partial_trace is not None:
metadata["partial_trace"] = partial_trace
elif trace_preserving:
metadata["trace_preserving"] = True
elif trace is not None:
metadata["trace"] = trace
t_stop = time.time()
metadata["fitter_time"] = t_stop - t_start
if len(fits) == 1:
return fits[0], metadata
return fits, metadata
[docs]
@cvxpy_utils.requires_cvxpy
def cvxpy_gaussian_lstsq(
outcome_data: np.ndarray,
shot_data: np.ndarray,
measurement_data: np.ndarray,
preparation_data: np.ndarray,
measurement_basis: Optional[MeasurementBasis] = None,
preparation_basis: Optional[PreparationBasis] = None,
measurement_qubits: Optional[Tuple[int, ...]] = None,
preparation_qubits: Optional[Tuple[int, ...]] = None,
conditional_measurement_indices: Optional[Tuple[int, ...]] = None,
trace: Union[None, float, str] = "auto",
psd: bool = True,
trace_preserving: Union[None, bool, str] = "auto",
partial_trace: Optional[np.ndarray] = None,
outcome_prior: Union[np.ndarray, int] = 0.5,
max_weight: float = 1e10,
**kwargs,
) -> Dict:
r"""Constrained Gaussian linear least-squares tomography fitter.
.. note::
This function calls :func:`cvxpy_linear_lstsq` with a Gaussian weights
vector. Refer to its documentation for additional details.
Overview
This fitter reconstructs the maximum-likelihood estimate by using
``cvxpy`` to minimize the constrained least-squares negative log
likelihood function
.. math::
\hat{\rho}
&= \mbox{argmin} (-\log\mathcal{L}{\rho}) \\
&= \mbox{argmin }\|W(Ax - y) \|_2^2 \\
-\log\mathcal{L}(\rho)
&= |W(Ax -y) \|_2^2 \\
&= \sum_i \frac{1}{\sigma_i^2}(\mbox{Tr}[E_j\rho] - \hat{p}_i)^2
Additional Details
The Gaussian weights are estimated from the observed frequency and shot data
via a Bayesian update of a Dirichlet distribution with observed outcome data
frequencies :math:`f_i(s)`, and Dirichlet prior :math:`\alpha_i(s)` for
tomography basis index `i` and measurement outcome `s`.
The mean posterior probabilities are computed as
.. math:
p_i(s) &= \frac{f_i(s) + \alpha_i(s)}{\bar{\alpha}_i + N_i} \\
Var[p_i(s)] &= \frac{p_i(s)(1-p_i(s))}{\bar{\alpha}_i + N_i + 1}
w_i(s) = \sqrt{Var[p_i(s)]}^{-1}
where :math:`N_i = \sum_s f_i(s)` is the total number of shots, and
:math:`\bar{\alpha}_i = \sum_s \alpha_i(s)` is the norm of the prior.
Args:
outcome_data: measurement outcome frequency data.
shot_data: basis measurement total shot data.
measurement_data: measurement basis index data.
preparation_data: preparation basis index data.
measurement_basis: Optional, measurement matrix basis.
preparation_basis: Optional, preparation matrix basis.
measurement_qubits: Optional, the physical qubits that were measured.
If None they are assumed to be ``[0, ..., M-1]`` for M measured qubits.
preparation_qubits: Optional, the physical qubits that were prepared.
If None they are assumed to be ``[0, ..., N-1]`` for N prepared qubits.
conditional_measurement_indices: Optional, conditional measurement data
indices. If set this will return a list of conditional fitted states
conditioned on a fixed basis measurement of these qubits.
trace: trace constraint for the fitted matrix. If "auto" this will be set
to 1 for QST or the input dimension for QST (default: "auto").
psd: If True rescale the eigenvalues of fitted matrix to be positive
semidefinite (default: True)
trace_preserving: Enforce the fitted matrix to be trace preserving when
fitting a Choi-matrix in quantum process
tomography. If "auto" this will be set to True for
QPT and False for QST (default: "auto").
partial_trace: Enforce conditional fitted Choi matrices to partial
trace to POVM matrices.
outcome_prior: The Bayesian prior :math:`\alpha` to use computing Gaussian
weights. See additional information.
max_weight: Set the maximum value allowed for weights vector computed from
tomography data variance.
kwargs: kwargs for cvxpy solver.
Raises:
QiskitError: If CVXPY is not installed on the current system.
AnalysisError: If analysis fails.
Returns:
The fitted matrix rho that maximizes the least-squares likelihood function.
"""
t_start = time.time()
weights = lstsq_utils.binomial_weights(
outcome_data,
shot_data=shot_data,
outcome_prior=outcome_prior,
max_weight=max_weight,
)
fits, metadata = cvxpy_linear_lstsq(
outcome_data,
shot_data,
measurement_data,
preparation_data,
measurement_basis=measurement_basis,
preparation_basis=preparation_basis,
measurement_qubits=measurement_qubits,
preparation_qubits=preparation_qubits,
conditional_measurement_indices=conditional_measurement_indices,
trace=trace,
psd=psd,
trace_preserving=trace_preserving,
partial_trace=partial_trace,
weights=weights,
**kwargs,
)
t_stop = time.time()
# Update metadata
metadata["fitter"] = "cvxpy_gaussian_lstsq"
metadata["fitter_time"] = t_stop - t_start
return fits, metadata