diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 28bbcb8e3..8e4947ad1 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -20,6 +20,10 @@ from pygsti import tools as _tools from pygsti.tools import mpitools as _mpit from pygsti.tools import slicetools as _slct +from pygsti.models import ( + ExplicitOpModel as _ExplicitOpModel +) +from pygsti.models.model import OpModel as _OpModel from pygsti.models.gaugegroup import ( TrivialGaugeGroupElement as _TrivialGaugeGroupElement, GaugeGroupElement as _GaugeGroupElement @@ -27,29 +31,70 @@ from pygsti.models import ExplicitOpModel from pygsti.modelmembers.operations import LinearOperator as _LinearOperator -from typing import Callable, Union, Optional +from typing import Callable, Union, Optional, Any -def gaugeopt_to_target(model, target_model, item_weights=None, - cptp_penalty_factor=0, spam_penalty_factor=0, - gates_metric="frobenius", spam_metric="frobenius", - gauge_group=None, method='auto', maxiter=100000, - maxfev=None, tol=1e-8, oob_check_interval=0, - convert_model_to=None, return_all=False, comm=None, - verbosity=0, check_jac=False, n_leak=0): +class GaugeoptToTargetArgs: + """ + This class is basically a namespace. It was introduced to strip out tons of complexity + in gaugeopt_to_target(...) without breaking old code that might call gaugeopt_to_target. """ - Optimize the gauge degrees of freedom of a model to that of a target. - - Parameters - ---------- - model : Model - The model to gauge-optimize - - target_model : Model - The model to optimize to. The metric used for comparing models - is given by `gates_metric` and `spam_metric`. - item_weights : dict, optional + old_trailing_positional_args = ( + 'item_weights', 'cptp_penalty_factor', 'spam_penalty_factor', + 'gates_metric', 'spam_metric', 'gauge_group', 'method', + 'maxiter', 'maxfev', 'tol', 'oob_check_interval', + 'convert_model_to', 'return_all', 'comm', 'verbosity', + 'check_jac' + ) + + @staticmethod + def parsed_model(model: Union[_OpModel, _ExplicitOpModel], convert_model_to: Optional[Any]) -> _OpModel: + if convert_model_to is None: + return model + + if not isinstance(model, _ExplicitOpModel): + raise ValueError('Gauge optimization only supports model conversion for ExplicitOpModels.') + conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) + model_out = model.copy() # don't alter the original model's parameterization (this would be unexpected) + for args in conversion_args: + if isinstance(args, str): + model_out.convert_members_inplace(args, set_default_gauge_group=True) + elif isinstance(args, dict): + model_out.convert_members_inplace(**args) + else: + raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) + return model_out + + @staticmethod + def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: str, spam_metric: str, n_leak: int) -> str: + ls_mode_allowed = bool(target_model is not None + and gates_metric.startswith("frobenius") + and spam_metric.startswith("frobenius") + and n_leak == 0) + # and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits + # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... + if method == "ls" and not ls_mode_allowed: + raise ValueError("Least-squares method is not allowed! Target" + " model must be non-None and frobenius metrics" + " must be used.") + elif method == "auto": + return 'ls' if ls_mode_allowed else 'L-BFGS-B' + else: + return method + + # The static dicts of default values are substituted in gaugeopt_to_target's kwargs. + # This is a safe thing to do because no invocation of "gaugeopt_to_target" within pyGSTi + # used positional arguments past target_model. + + create_objective_passthrough_kwargs : dict[str,Any] = dict( + item_weights=None, # gets cast to a dict. + cptp_penalty_factor=0, + spam_penalty_factor=0, + check_jac=False, + ) + """ + item_weights : dict Dictionary of weighting factors for gates and spam operators. Keys can be gate, state preparation, or POVM effect, as well as the special values "spam" or "gates" which apply the given weighting to *all* spam operators @@ -58,16 +103,56 @@ def gaugeopt_to_target(model, target_model, item_weights=None, "spam" values. The precise use of these weights depends on the model metric(s) being used. - cptp_penalty_factor : float, optional + cptp_penalty_factor : float If greater than zero, the objective function also contains CPTP penalty terms which penalize non-CPTP-ness of the gates being optimized. This factor multiplies these CPTP penalty terms. - spam_penalty_factor : float, optional + spam_penalty_factor : float If greater than zero, the objective function also contains SPAM penalty terms which penalize non-positive-ness of the state preps being optimized. This factor multiplies these SPAM penalty terms. + check_jac : bool + When True, check least squares analytic jacobian against finite differences. + """ + + gaugeopt_custom_passthrough_kwargs : dict[str, Any] = dict( + maxiter=100000, + maxfev=None, + tol=1e-8, + oob_check_interval=0, + verbosity=0, + ) + """ + maxiter : int + Maximum number of iterations for the gauge optimization. + + maxfev : int + Maximum number of function evaluations for the gauge optimization. + Defaults to maxiter. + + tol : float + The tolerance for the gauge optimization. + + oob_check_interval : int + If greater than zero, gauge transformations are allowed to fail (by raising + any exception) to indicate an out-of-bounds condition that the gauge optimizer + will avoid. If zero, then any gauge-transform failures just terminate the + optimization. + """ + + other_kwargs : dict[str, Any] = dict( + gates_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + spam_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + gauge_group=None, # used iff convert_model_to is not None, then passed to _create_objective_fn and gaugeopt_custom + method='auto', # validated, then passed to _create_objective_fn and gaugeopt_custom + return_all=False, # used in the function body (only for branching after passing to gaugeopt_custom) + convert_model_to=None, # passthrough to parsed_model. + comm=None, # passthrough to _create_objective_fn and gaugeopt_custom + n_leak=0 # passthrough to parsed_objective and _create_objective_fn + ) + """ gates_metric : {"frobenius", "fidelity", "tracedist"}, optional The metric used to compare gates within models. "frobenius" computes the normalized sqrt(sum-of-squared-differences), with weights @@ -100,22 +185,6 @@ def gaugeopt_to_target(model, target_model, item_weights=None, - 'evolve' -- evolutionary global optimization algorithm using DEAP - 'brute' -- Experimental: scipy.optimize.brute using 4 points along each dimensions - maxiter : int, optional - Maximum number of iterations for the gauge optimization. - - maxfev : int, optional - Maximum number of function evaluations for the gauge optimization. - Defaults to maxiter. - - tol : float, optional - The tolerance for the gauge optimization. - - oob_check_interval : int, optional - If greater than zero, gauge transformations are allowed to fail (by raising - any exception) to indicate an out-of-bounds condition that the gauge optimizer - will avoid. If zero, then any gauge-transform failures just terminate the - optimization. - convert_model_to : str, dict, list, optional For use when `model` is an `ExplicitOpModel`. When not `None`, calls `model.convert_members_inplace(convert_model_to, set_default_gauge_group=False)` if @@ -132,58 +201,103 @@ def gaugeopt_to_target(model, target_model, item_weights=None, When not None, an MPI communicator for distributing the computation across multiple processors. - verbosity : int, optional - How much detail to send to stdout. + n_leak : int + Used in leakage modeling. If positive, this specifies how modify defintiions + of gate and SPAM metrics to reflect the fact that target gates do not have + well-defined actions outside the computational subspace. + """ - check_jac : bool - When True, check least squares analytic jacobian against finite differences + @staticmethod + def create_full_kwargs(args: tuple[Any,...], kwargs: dict[str, Any]): + full_kwargs = GaugeoptToTargetArgs.create_objective_passthrough_kwargs.copy() + full_kwargs.update(GaugeoptToTargetArgs.gaugeopt_custom_passthrough_kwargs) + full_kwargs.update(GaugeoptToTargetArgs.other_kwargs) + full_kwargs.update(kwargs) + + if extra_posargs := len(args) > 0: + msg = \ + f""" + Recieved {extra_posargs} positional arguments past `model` and `target_model`. + These should be passed as appropriate keyword arguments instead. This version + of pyGSTi will infer intended keyword arguments based on the legacy argument + positions. Future versions of pyGSTi will raise an error. + """ + _warnings.warn(msg) + for k,v in zip(GaugeoptToTargetArgs.old_trailing_positional_args, args): + full_kwargs[k] = v + + if full_kwargs['item_weights'] is None: + full_kwargs['item_weights'] = dict() + + return full_kwargs + + +def gaugeopt_to_target(model, target_model, *args, **kwargs): + """ + Legacy function to optimize the gauge degrees of freedom of a model to that of a target. + + Use of more than two positional arguments is deprecated; keyword arguments should be used + instead. See the GaugeoptToTargetArgs class for available keyword arguments and their + default values. Returns ------- - model : if return_all == False + model : if kwargs.get('return_all', False) == False + + (goodnessMin, gaugeMx, model) : if kwargs.get('return_all', False) == True - (goodnessMin, gaugeMx, model) : if return_all == True Where goodnessMin is the minimum value of the goodness function (the best 'goodness') found, gaugeMx is the gauge matrix used to transform the model, and model is the final gauge-transformed model. """ - - if item_weights is None: item_weights = {} - - ls_mode_allowed = bool(target_model is not None - and gates_metric.startswith("frobenius") - and spam_metric.startswith("frobenius")) - #and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits - # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... - if method == "ls" and not ls_mode_allowed: - raise ValueError("Least-squares method is not allowed! Target" - " model must be non-None and frobenius metrics" - " must be used.") - if method == "auto": - method = 'ls' if ls_mode_allowed else 'L-BFGS-B' + full_kwargs = GaugeoptToTargetArgs.create_full_kwargs(args, kwargs) + """ + This function handles a strange situation where `target_model` can be None. - if convert_model_to is not None: - conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) - model = model.copy() # don't alter the original model's parameterization (this would be unexpected) - for args in conversion_args: - if isinstance(args, str): - model.convert_members_inplace(args, set_default_gauge_group=True) - elif isinstance(args, dict): - model.convert_members_inplace(**args) - else: - raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) + In this case, the objective function will only depend on `cptp_penality_factor` + and `spam_penalty_factor`; it's forbidden to use method == 'ls'; and the + returned OpModel might not have its basis set. + """ + # arg parsing: validating `method` + gates_metric = full_kwargs['gates_metric'] + spam_metric = full_kwargs['spam_metric'] + n_leak = full_kwargs['n_leak'] + method = full_kwargs['method'] + method = GaugeoptToTargetArgs.parsed_method( + target_model, method, gates_metric, spam_metric, n_leak + ) + + # arg parsing: (possibly) converting `model` + convert_model_to = full_kwargs['convert_model_to'] + model = GaugeoptToTargetArgs.parsed_model(model, convert_model_to) + + # actual work: constructing objective_fn and jacobian_fn + item_weights = full_kwargs['item_weights'] + cptp_penalty = full_kwargs['cptp_penalty_factor'] + spam_penalty = full_kwargs['spam_penalty_factor'] + comm = full_kwargs['comm'] + check_jac = full_kwargs['check_jac'] objective_fn, jacobian_fn = _create_objective_fn( - model, target_model, item_weights, - cptp_penalty_factor, spam_penalty_factor, - gates_metric, spam_metric, method, comm, check_jac, n_leak) - - result = gaugeopt_custom(model, objective_fn, gauge_group, method, - maxiter, maxfev, tol, oob_check_interval, - return_all, jacobian_fn, comm, verbosity) - - #If we've gauge optimized to a target model, declare that the + model, target_model, item_weights, cptp_penalty, spam_penalty, + gates_metric, spam_metric, method, comm, check_jac, n_leak + ) + + # actual work: calling the (wrapper of the wrapper of the ...) optimizer + gauge_group = full_kwargs['gauge_group'] + maxiter = full_kwargs['maxiter'] + maxfev = full_kwargs['maxfev'] + tol = full_kwargs['tol'] + oob_check = full_kwargs['oob_check_interval'] + return_all = full_kwargs['return_all'] + verbosity = full_kwargs['verbosity'] + result = gaugeopt_custom( + model, objective_fn, gauge_group, method, maxiter, maxfev, + tol, oob_check, return_all, jacobian_fn, comm, verbosity + ) + + # If we've gauge optimized to a target model, declare that the # resulting model is now in the same basis as the target. if target_model is not None: newModel = result[-1] if return_all else result @@ -192,9 +306,14 @@ def gaugeopt_to_target(model, target_model, item_weights=None, return result -def gaugeopt_custom(model, objective_fn, gauge_group=None, +GGElObjective = Callable[[_GaugeGroupElement, bool], Union[_np.floating, _np.ndarray]] + +GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] + + +def gaugeopt_custom(model, objective_fn: GGElObjective, gauge_group=None, method='L-BFGS-B', maxiter=100000, maxfev=None, tol=1e-8, - oob_check_interval=0, return_all=False, jacobian_fn=None, + oob_check_interval=0, return_all=False, jacobian_fn: GGElJacobian=None, comm=None, verbosity=0): """ Optimize the gauge of a model using a custom objective function. @@ -205,8 +324,9 @@ def gaugeopt_custom(model, objective_fn, gauge_group=None, The model to gauge-optimize objective_fn : function - The function to be minimized. The function must take a single `Model` - argument and return a float. + The function to be minimized. The function must take a GaugeGroupElement + and a bool. If method == 'ls' then objective_fn must return an ndarray; if + method != 'ls' then objective_fn must return a float. gauge_group : GaugeGroup, optional The gauge group which defines which gauge trasformations are optimized @@ -358,17 +478,70 @@ def _call_jacobian_fn(gauge_group_el_vec): return newModel -GGElObjective = Callable[[_GaugeGroupElement,bool], Union[float, _np.ndarray]] +def _transform_with_oob_check(mdl, gauge_group_el, oob_check): + """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ + mdl = mdl.copy() + if oob_check: + try: + mdl.transform_inplace(gauge_group_el) + except Exception as e: + raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition + else: + mdl.transform_inplace(gauge_group_el) + return mdl -GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] -LabelLike = Union[str, _baseobjs.Label] +def gates_with_instruments(model) -> dict[_baseobjs.Label, _LinearOperator]: + gates = {k:v for (k,v) in model.operations.items()} + for lbl, inst in model.instruments.items(): + gates.update(inst.simplify_operations(lbl)) + return gates -def _create_objective_fn(model, target_model, item_weights: Optional[dict[LabelLike, float]]=None, - cptp_penalty_factor: float=0.0, spam_penalty_factor: float=0.0, - gates_metric="frobenius", spam_metric="frobenius", - method=None, comm=None, check_jac=False, n_leak=0) -> tuple[GGElObjective, GGElJacobian]: +def _gate_fidelity_targets(model, target_model): + from pygsti.report.reportables import eigenvalue_entanglement_infidelity + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + + mdl_ops = gates_with_instruments(model) + tgt_ops = gates_with_instruments(target_model) + + for lbl in tgt_ops: + G_target = tgt_ops[lbl].to_dense() + G_curest = mdl_ops[lbl].to_dense() + t = 1 - eigenvalue_entanglement_infidelity(G_curest, G_target, model.basis) + t = _np.clip(t, a_min=0.0, a_max=1.0) + gate_fidelity_targets[lbl] = (G_target, t) + return gate_fidelity_targets + + +def _prep_fidelity_targets(model, target_model): + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + for preplbl in target_model.preps: + rho_target = _tools.vec_to_stdmx( target_model.preps[preplbl].to_dense(), model.basis ) + rho_curest = _tools.vec_to_stdmx( model.preps[preplbl].to_dense(), model.basis ) + t = _tools.eigenvalue_fidelity(rho_curest, rho_target) + t = _np.clip(t, a_min=0.0, a_max=1.0) + prep_fidelity_targets[preplbl] = (rho_target, t) + return prep_fidelity_targets + + +def _povm_effect_fidelity_targets(model, target_model): + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + for povmlbl in target_model.povms: + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_target.items() } + Es_curest = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_curest.items() } + ts = {elbl : _tools.eigenvalue_fidelity(Es_target[elbl], Es_curest[elbl]) for elbl in M_target } + ts = {elbl : _np.clip(t, a_min=0.0, a_max=1.0) for (elbl, t) in ts.items()} + povm_fidelity_targets[povmlbl] = (Es_target, ts) + return povm_fidelity_targets + + +def _legacy_create_scalar_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, n_leak: int + ) -> tuple[GGElObjective, GGElJacobian]: """ Creates the objective function and jacobian (if available) for gaugeopt_to_target @@ -384,353 +557,392 @@ def _create_objective_fn(model, target_model, item_weights: Optional[dict[LabelL if mxBasis.name == "unknown" and target_model is not None: mxBasis = target_model.basis - def _transform_with_oob_check(mdl, gauge_group_el, oob_check): - """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ - mdl = mdl.copy() - if oob_check: - try: mdl.transform_inplace(gauge_group_el) - except Exception as e: - raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition + assert gates_metric != "frobeniustt" + assert spam_metric != "frobeniustt" + # ^ PR #410 removed support for Frobenius transform-target metrics in this codepath. + + dim = int(_np.sqrt(mxBasis.dim)) + I = _tools.IdentityOperator() + if n_leak > 0: + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + P = B @ B.T.conj() + if _np.linalg.norm(P.imag) > 1e-12: + msg = f"Attempting to run leakage-aware gauge optimization with basis {mxBasis}\n" + msg += "is resulting an orthogonal projector onto the computational subspace that\n" + msg += "is not real-valued. Try again with a different basis, like 'l2p1' or 'gm'." + raise ValueError(msg) else: - mdl.transform_inplace(gauge_group_el) - return mdl + P = P.real + else: + P = I + transform_mx_arg = (P, I) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. + + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + if gates_metric == 'fidelity': + gate_fidelity_targets.update(_gate_fidelity_targets(model, target_model)) + + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + if spam_metric == 'fidelity': + prep_fidelity_targets.update(_prep_fidelity_targets(model, target_model)) + povm_fidelity_targets.update(_povm_effect_fidelity_targets(model, target_model)) + + tgt_ops : dict[_baseobjs.Label, _LinearOperator] = dict() + if target_model is not None: + tgt_ops.update(gates_with_instruments(target_model)) - if method == "ls": - # least-squares case where objective function returns an array of - # the before-they're-squared difference terms and there's an analytic jacobian + def _objective_fn(gauge_group_el, oob_check): + mdl : ExplicitOpModel = _transform_with_oob_check(model, gauge_group_el, oob_check) + mdl_ops : dict[_baseobjs.Label, _LinearOperator] = gates_with_instruments(mdl) + ret: _np.floating = _np.float64(0.0) - assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ - "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" - assert(gates_metric == spam_metric) - frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" + if cptp_penalty_factor > 0: + mdl.basis = mxBasis # set basis for jamiolkowski iso + cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) + ret += _np.sum(cpPenaltyVec) + + if spam_penalty_factor > 0: + mdl.basis = mxBasis + spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) + ret += _np.sum(spamPenaltyVec) + + if target_model is None: + return ret + + if "frobenius" in gates_metric: + if spam_metric == gates_metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + else: + wts = item_weights.copy() + wts['spam'] = 0.0 + for k in wts: + if k in mdl.preps or k in mdl.povms: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in gates_metric: + val = val ** 2 + ret += val + + elif gates_metric == "fidelity": + # Leakage-aware metrics NOT available + for opLbl in mdl_ops: + wt = item_weights.get(opLbl, opWeight) + mop = mdl_ops[opLbl].to_dense() + top, t = gate_fidelity_targets[opLbl] + v = _tools.entanglement_fidelity(mop, top, mxBasis) + z = _np.abs(t - v) + ret += wt * z + + elif gates_metric == "tracedist": + # If n_leak==0, then subspace_jtracedist is just jtracedist. + for opLbl in mdl_ops: + wt = item_weights.get(opLbl, opWeight) + top = tgt_ops[opLbl].to_dense() + mop = mdl_ops[opLbl].to_dense() + ret += wt * _tools.subspace_jtracedist(top, mop, mxBasis, n_leak) - if frobenius_transform_target: - full_target_model = target_model.copy() - full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. else: - full_target_model = None # in case it get's referenced by mistake + raise ValueError("Invalid gates_metric: %s" % gates_metric) - def _objective_fn(gauge_group_el, oob_check): + if "frobenius" in spam_metric and gates_metric == spam_metric: + # We already handled SPAM error in this case. Just return. + return ret - if frobenius_transform_target: - transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) - other = model - else: - transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) - other = target_model + if "frobenius" in spam_metric: + # SPAM and gates can have different choices for squared vs non-squared. + wts = item_weights.copy(); wts['gates'] = 0.0 + for k in wts: + if k in mdl_ops or k in mdl.instruments: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in spam_metric: + val = val ** 2 + ret += val + + elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available + val_prep = 0.0 + for preplbl in target_model.preps: + wt_prep = item_weights.get(preplbl, spamWeight) + rho_curest = _tools.vec_to_stdmx(model.preps[preplbl].to_dense(), mxBasis) + rho_target, t = prep_fidelity_targets[preplbl] + v = _tools.fidelity(rho_curest, rho_target) + val_prep += wt_prep * abs(t - v) + val_povm = 0.0 + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target, ts = povm_fidelity_targets[povmlbl] + for elbl in M_target: + t_e = ts[elbl] + E = _tools.vec_to_stdmx(M_curest[elbl].to_dense(), mxBasis) + v_e = _tools.fidelity(E, Es_target[elbl]) + val_povm += wt_povm * abs(t_e - v_e) + ret += (val_prep + val_povm) # type: ignore + + elif spam_metric == "tracedist": + # Leakage-aware metrics NOT available. + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + else: + raise ValueError("Invalid spam_metric: %s" % spam_metric) - residuals, _ = transformed.residuals(other, None, item_weights) + return ret - # We still the non-target model to be transformed and checked for these penalties - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: - transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + return _objective_fn, None - if cptp_penalty_factor > 0: - transformed.basis = mxBasis - cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) - else: cpPenaltyVec = [] # so concatenate ignores - if spam_penalty_factor > 0: - transformed.basis = mxBasis - spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) - else: spamPenaltyVec = [] # so concatenate ignores +def _legacy_create_least_squares_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, comm: Optional[Any], check_jac: bool + ) -> tuple[GGElObjective, GGElJacobian]: - return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) - else: - return residuals + opWeight = item_weights.get('gates', 1.0) + spamWeight = item_weights.get('spam', 1.0) + mxBasis = model.basis + + #Use the target model's basis if model's is unknown + # (as it can often be if it's just come from an logl opt, + # since from_vector will clear any basis info) + if mxBasis.name == "unknown" and target_model is not None: + mxBasis = target_model.basis + + assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ + "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" + assert(gates_metric == spam_metric) + frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" + + if frobenius_transform_target: + full_target_model = target_model.copy() + full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. + else: + full_target_model = None # in case it get's referenced by mistake - def _jacobian_fn(gauge_group_el): + def _objective_fn(gauge_group_el: _GaugeGroupElement, oob_check: bool) -> _np.ndarray: - #Penalty terms below always act on the transformed non-target model. - original_gauge_group_el = gauge_group_el + if frobenius_transform_target: + transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) + other = model + else: + transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + other = target_model + residuals, _ = transformed.residuals(other, None, item_weights) + + # We still the non-target model to be transformed and checked for these penalties + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: if frobenius_transform_target: - gauge_group_el = gauge_group_el.inverse() - mdl_pre = full_target_model.copy() - mdl_post = mdl_pre.copy() - else: - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - # Indices: Jacobian output matrix has shape (L, N) - start = 0 - d = mdl_pre.dim - N = gauge_group_el.num_params - L = mdl_pre.num_elements - - #Compute "extra" (i.e. beyond the model-element) rows of jacobian - if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) - if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) - - #Set basis for pentaly term calculation - if cptp_penalty_factor != 0 or spam_penalty_factor != 0: - mdl_pre.basis = mxBasis - mdl_post.basis = mxBasis - - jacMx = _np.zeros((L, N)) - - #Overview of terms: - # objective: op_term = (S_inv * gate * S - target_op) - # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) - # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) - - # objective: rho_term = (S_inv * rho - target_rho) - # jac: d(rho_term) = d (S_inv) * rho - # d(rho_term) = -(S_inv * dS * S_inv) * rho - - # objective: ET_term = (E.T * S - target_E.T) - # jac: d(ET_term) = E.T * dS - - #Overview of terms when frobenius_transform_target == True). Note that the objective - #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. - - # objective: op_term = (gate - S * target_op * S_inv) - # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) - # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) - - # objective: rho_term = (rho - S * target_rho) - # jac: d(rho_term) = - dS * target_rho - - # objective: ET_term = (E.T - target_E.T * S_inv) - # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) - # d(ET_term) = target_E.T * (S_inv * dS * S_inv) - - #Distribute computation across processors - allDerivColSlice = slice(0, N) - derivSlices, myDerivColSlice, derivOwners, mySubComm = \ - _mpit.distribute_slice(allDerivColSlice, comm) - if mySubComm is not None: - _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() - + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover - - n = _slct.length(myDerivColSlice) - wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None - my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for - - # S, and S_inv are shape (d,d) - #S = gauge_group_el.transform_matrix - S_inv = gauge_group_el.transform_matrix_inverse - dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n - dS.shape = (d, d, n) # call it (d1,d2,n) - dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) - assert(dS.shape == (n, d, d)) - - # --- NOTE: ordering here, with running `start` index MUST - # correspond to those in Model.residuals, which in turn - # must correspond to those in ForwardSimulator.residuals - which - # currently orders as: gates, simplified_ops, preps, effects. - - # -- LinearOperator terms - # ------------------------- - for lbl, G in mdl_pre.operations.items(): - # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) - # Note: (S_inv * G * S) is G' (transformed G) - wt = item_weights.get(lbl, opWeight) - left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) + transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + + if cptp_penalty_factor > 0: + transformed.basis = mxBasis + cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) + else: cpPenaltyVec = [] # so concatenate ignores + + if spam_penalty_factor > 0: + transformed.basis = mxBasis + spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) + else: spamPenaltyVec = [] # so concatenate ignores + + return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) + else: + return residuals + + def _jacobian_fn(gauge_group_el: _GaugeGroupElement) -> _np.ndarray: + + #Penalty terms below always act on the transformed non-target model. + original_gauge_group_el = gauge_group_el + + if frobenius_transform_target: + gauge_group_el = gauge_group_el.inverse() + mdl_pre = full_target_model.copy() + mdl_post = mdl_pre.copy() + else: + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) + + # Indices: Jacobian output matrix has shape (L, N) + start = 0 + d = mdl_pre.dim + N = gauge_group_el.num_params + L = mdl_pre.num_elements + + #Compute "extra" (i.e. beyond the model-element) rows of jacobian + if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) + if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) + + #Set basis for pentaly term calculation + if cptp_penalty_factor != 0 or spam_penalty_factor != 0: + mdl_pre.basis = mxBasis + mdl_post.basis = mxBasis + + jacMx = _np.zeros((L, N)) + + #Overview of terms: + # objective: op_term = (S_inv * gate * S - target_op) + # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) + # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) + + # objective: rho_term = (S_inv * rho - target_rho) + # jac: d(rho_term) = d (S_inv) * rho + # d(rho_term) = -(S_inv * dS * S_inv) * rho + + # objective: ET_term = (E.T * S - target_E.T) + # jac: d(ET_term) = E.T * dS + + #Overview of terms when frobenius_transform_target == True). Note that the objective + #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. + + # objective: op_term = (gate - S * target_op * S_inv) + # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) + # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) + + # objective: rho_term = (rho - S * target_rho) + # jac: d(rho_term) = - dS * target_rho + + # objective: ET_term = (E.T - target_E.T * S_inv) + # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) + # d(ET_term) = target_E.T * (S_inv * dS * S_inv) + + #Distribute computation across processors + allDerivColSlice = slice(0, N) + derivSlices, myDerivColSlice, derivOwners, mySubComm = \ + _mpit.distribute_slice(allDerivColSlice, comm) + if mySubComm is not None: + _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() + + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover + + n = _slct.length(myDerivColSlice) + wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None + my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for + + # S, and S_inv are shape (d,d) + #S = gauge_group_el.transform_matrix + S_inv = gauge_group_el.transform_matrix_inverse + dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n + dS.shape = (d, d, n) # call it (d1,d2,n) + dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) + assert(dS.shape == (n, d, d)) + + # --- NOTE: ordering here, with running `start` index MUST + # correspond to those in Model.residuals, which in turn + # must correspond to those in ForwardSimulator.residuals - which + # currently orders as: gates, simplified_ops, preps, effects. + + # -- LinearOperator terms + # ------------------------- + for lbl, G in mdl_pre.operations.items(): + # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) + # Note: (S_inv * G * S) is G' (transformed G) + wt = item_weights.get(lbl, opWeight) + left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) + result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) + result = result.reshape((d**2, n)) # must copy b/c non-contiguous + my_jacMx[start:start + d**2] = wt * result + start += d**2 + + # -- Instrument terms + # ------------------------- + for ilbl, Inst in mdl_pre.instruments.items(): + wt = item_weights.get(ilbl, opWeight) + for lbl, G in Inst.items(): + # same calculation as for operation terms + left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) result = result.reshape((d**2, n)) # must copy b/c non-contiguous my_jacMx[start:start + d**2] = wt * result start += d**2 - # -- Instrument terms - # ------------------------- - for ilbl, Inst in mdl_pre.instruments.items(): - wt = item_weights.get(ilbl, opWeight) - for lbl, G in Inst.items(): - # same calculation as for operation terms - left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) - result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) - result = result.reshape((d**2, n)) # must copy b/c non-contiguous - my_jacMx[start:start + d**2] = wt * result - start += d**2 - - # -- prep terms - # ------------------------- - for lbl, rho in mdl_post.preps.items(): - # d(rho_term) = -(S_inv * dS * S_inv) * rho - # Note: (S_inv * rho) is transformed rho - wt = item_weights.get(lbl, spamWeight) - Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) - result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) - my_jacMx[start:start + d] = wt * result + # -- prep terms + # ------------------------- + for lbl, rho in mdl_post.preps.items(): + # d(rho_term) = -(S_inv * dS * S_inv) * rho + # Note: (S_inv * rho) is transformed rho + wt = item_weights.get(lbl, spamWeight) + Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) + result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) + my_jacMx[start:start + d] = wt * result + start += d + + # -- effect terms + # ------------------------- + for povmlbl, povm in mdl_pre.povms.items(): + for lbl, E in povm.items(): + # d(ET_term) = E.T * dS + wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) + result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) + my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) start += d - # -- effect terms - # ------------------------- - for povmlbl, povm in mdl_pre.povms.items(): - for lbl, E in povm.items(): - # d(ET_term) = E.T * dS - wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) - result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) - my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) - start += d - - # -- penalty terms -- Note: still use original gauge transform applied to `model` - # ------------------------- - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: # reset back to non-target-tranform "mode" - gauge_group_el = original_gauge_group_el - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - if cptp_penalty_factor > 0: - start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, cptp_penalty_factor, - mdl_pre.basis, wrtIndices) - - if spam_penalty_factor > 0: - start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, spam_penalty_factor, - mdl_pre.basis, wrtIndices) - - #At this point, each proc has filled the portions (columns) of jacMx that - # it's responsible for, and so now we gather them together. - _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) - #Note jacMx is completely filled (on all procs) - - if check_jac and (comm is None or comm.Get_rank() == 0): - def _mock_objective_fn(v): - return _objective_fn(gauge_group_el, False) - vec = gauge_group_el.to_vector() - _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', - verbosity=1) - - return jacMx - - else: - # non-least-squares case where objective function returns a single float - # and (currently) there's no analytic jacobian - - assert gates_metric != "frobeniustt" - assert spam_metric != "frobeniustt" - # ^ PR #410 removed support for Frobenius transform-target metrics in this codepath. - - dim = int(_np.sqrt(mxBasis.dim)) - if n_leak > 0: - B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) - P = B @ B.T.conj() - if _np.linalg.norm(P.imag) > 1e-12: - msg = f"Attempting to run leakage-aware gauge optimization with basis {mxBasis}\n" - msg += "is resulting an orthogonal projector onto the computational subspace that\n" - msg += "is not real-valued. Try again with a different basis, like 'l2p1' or 'gm'." - raise ValueError(msg) - else: - P = P.real - transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) - # ^ The semantics of this tuple are defined by the frobeniusdist function - # in the ExplicitOpModelCalc class. - else: - transform_mx_arg = None - # ^ It would be equivalent to set this to a pair of IdentityOperator objects. - - def _objective_fn(gauge_group_el, oob_check): - mdl : ExplicitOpModel = _transform_with_oob_check(model, gauge_group_el, oob_check) - mdl_ops : dict[_baseobjs.Label, _LinearOperator] = mdl._excalc().operations - tgt_ops : dict[_baseobjs.Label, _LinearOperator] = dict() - if target_model is not None: - tgt_ops.update(target_model._excalc().operations) - # ^ Use these dicts instead of mdl.operations and target_model.operations, - # since these dicts are updated to include instruments. - ret = 0.0 + # -- penalty terms -- Note: still use original gauge transform applied to `model` + # ------------------------- + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: + if frobenius_transform_target: # reset back to non-target-tranform "mode" + gauge_group_el = original_gauge_group_el + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) if cptp_penalty_factor > 0: - mdl.basis = mxBasis # set basis for jamiolkowski iso - cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) - ret += _np.sum(cpPenaltyVec) + start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, cptp_penalty_factor, + mdl_pre.basis, wrtIndices) if spam_penalty_factor > 0: - mdl.basis = mxBasis - spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) - ret += _np.sum(spamPenaltyVec) - - if target_model is None: - return ret - - if "frobenius" in gates_metric: - if spam_metric == gates_metric: - val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) - else: - non_gates = {'spam'}.union(set(mdl.preps)).union(set(mdl.povms)) - temp_wts = {k: (0.0 if k in non_gates else v) for (k, v) in item_weights.items()} - val = mdl.frobeniusdist(target_model, transform_mx_arg, temp_wts) - if "squared" in gates_metric: - val = val ** 2 - ret += val - - elif gates_metric == "fidelity": - # If n_leak==0, then subspace_entanglement_fidelity is just entanglement_fidelity - for opLbl in mdl_ops: - wt = item_weights.get(opLbl, opWeight) - top = tgt_ops[opLbl].to_dense() - mop = mdl_ops[opLbl].to_dense() - ret += wt * (1.0 - _tools.subspace_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 - - elif gates_metric == "tracedist": - # If n_leak==0, then subspace_jtracedist is just jtracedist. - for opLbl in mdl_ops: - wt = item_weights.get(opLbl, opWeight) - top = tgt_ops[opLbl].to_dense() - mop = mdl_ops[opLbl].to_dense() - ret += wt * _tools.subspace_jtracedist(top, mop, mxBasis, n_leak) + start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, spam_penalty_factor, + mdl_pre.basis, wrtIndices) - else: - raise ValueError("Invalid gates_metric: %s" % gates_metric) - - if "frobenius" in spam_metric and gates_metric == spam_metric: - # We already handled SPAM error in this case. Just return. - return ret - - if "frobenius" in spam_metric: - # SPAM and gates can have different choices for squared vs non-squared. - non_spam = {'gates'}.union(set(mdl_ops)) # use mdl_ops, not mdl.operations! - temp_wts = {k: (0.0 if k in non_spam else v) for (k, v) in item_weights.items()} - val = mdl.frobeniusdist(target_model, transform_mx_arg, temp_wts) - if "squared" in spam_metric: - val = val ** 2 - ret += val - - elif spam_metric == "fidelity": - # Leakage-aware metrics NOT available - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) - ret += wt * (1.0 - fidelity)**2 - - elif spam_metric == "tracedist": - # Leakage-aware metrics NOT available. - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + #At this point, each proc has filled the portions (columns) of jacMx that + # it's responsible for, and so now we gather them together. + _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) + #Note jacMx is completely filled (on all procs) - else: - raise ValueError("Invalid spam_metric: %s" % spam_metric) + if check_jac and (comm is None or comm.Get_rank() == 0): + def _mock_objective_fn(v): + return _objective_fn(gauge_group_el, False) + vec = gauge_group_el.to_vector() + _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', + verbosity=1) - return ret + return jacMx - _jacobian_fn = None return _objective_fn, _jacobian_fn +def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,float]]=None, + cptp_penalty_factor: float=0.0, spam_penalty_factor: float=0.0, + gates_metric="frobenius", spam_metric="frobenius", + method=None, comm=None, check_jac=False, n_leak: int=0) -> tuple[GGElObjective, GGElJacobian]: + if item_weights is None: + item_weights = dict() + if method == 'ls': + return _legacy_create_least_squares_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, comm, check_jac + ) + else: + return _legacy_create_scalar_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, n_leak + ) + + def _cptp_penalty_size(mdl): """ Helper function - *same* as that in core.py. diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index 5d276d9e8..7e3ae7c30 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import collections as _collections import itertools as _itertools import uuid as _uuid @@ -1635,6 +1637,10 @@ def add_availability(opkey, op): availability, instrument_names=list(self.instruments.keys()), nonstd_instruments=self.instruments) + def copy(self) -> ExplicitOpModel: + c = super().copy() # <-- that is indeed an ExplicitOpModel + return c # type: ignore + def create_modelmember_graph(self): return _MMGraph({ 'preps': self.preps, diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 5f6f8e940..562118a50 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import bisect as _bisect import copy as _copy import itertools as _itertools @@ -326,7 +328,7 @@ def _post_copy(self, copy_into, memo): """ pass - def copy(self): + def copy(self) -> Model: """ Copy this model. @@ -2309,7 +2311,7 @@ def _post_copy(self, copy_into, memo): copy_into._reinit_opcaches() super(OpModel, self)._post_copy(copy_into, memo) - def copy(self): + def copy(self) -> OpModel: """ Copy this model. @@ -2320,7 +2322,7 @@ def copy(self): """ self._clean_paramvec() # ensure _paramvec is rebuilt if needed if OpModel._pcheck: self._check_paramvec() - ret = Model.copy(self) + ret = Model.copy(self) # <-- don't be fooled. That's an OpModel! if self._param_bounds is not None and self.parameter_labels is not None: ret._clean_paramvec() # will *always* rebuild paramvec; do now so we can preserve param bounds assert _np.all(self.parameter_labels == ret.parameter_labels) # ensure ordering is the same diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index f92dc8d36..e6e4af76c 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -17,6 +17,7 @@ """ import importlib import warnings as _warnings +from typing import Union import numpy as _np import scipy.linalg as _spl @@ -40,6 +41,8 @@ FINITE_DIFF_EPS = 1e-7 +BasisLike = Union[_Basis, str] + def _null_fn(*arg): return None @@ -1577,32 +1580,67 @@ def eigenvalue_nonunitary_avg_gate_infidelity(a, b, mx_basis): # init args == (model1, model2, op_label) -def eigenvalue_entanglement_infidelity(a, b, mx_basis): + +def eigenvalue_entanglement_infidelity(a: _np.ndarray, b: _np.ndarray, mx_basis: BasisLike, is_tp=None, is_unitary=None, tol=1e-8) -> _np.floating: """ - Eigenvalue entanglement infidelity between a and b + Eigenvalue entanglement infidelity is the infidelity between certain diagonal matrices that + contain [see Note 1] the eigenvalues of J(a) and J(b), where J(.) is the Jamiolkowski + isomorphism map. This is equivalent to computing infidelity of (J(a), J(b)) while pretending + that they share an eigenbasis. Parameters ---------- - a : numpy.ndarray - The first process (transfer) matrix. + is_tp : bool, optional (default None) + Flag indicating that a and b are TP in their provided basis. If None (the default), + an explicit check is performed up to numerical tolerance `tol`. If True/False, then + the check is skipped and this is used as the result of the check as though it were + performed. - b : numpy.ndarray - The second process (transfer) matrix. + is_unitary : bool, optional (default None) + Flag indicating that b is unitary. If None (the default) an explicit check is performed + up to tolerance `tol`. If True/False, then the check is skipped and this is used as the + result of the check as though it were performed. - mx_basis : Basis or {'pp', 'gm', 'std'} - the basis that `a` and `b` are in. + Notes + ----- + [1] The eigenvalues are ordered using a min-weight matching algorithm. - Returns - ------- - float + [2] If a and b are trace-preserving (TP) and b is unitary, then we compute eigenvalue + entanglement fidelity by leveraging three facts: + + 1. If (x,y) share an eigenbasis and have (consistently ordered) eigenvalues (v(x),v(y)), + then Tr(x y^{\\dagger}) is the inner product of v(x) and v(y). + + 2. J(a) and J(b) share an eigenbasis of if `a` and `b` share an eigenbasis. + + 3. If `a` and `b` are TP and `b` is unitary, then their entanglement fidelity can be + expressed as |Tr( a.b^{\\dagger} )| / d^2. + + Then together, we see that if (a,b) are TP and `b` is unitary, then their eigenvalue + entanglement fidelity is expressible as the inner product of v(a) and v(b), where + v(a) and v(b) vectors holding are suitably-ordered eigenvalues of a and b. + """ d2 = a.shape[0] - evA = _np.linalg.eigvals(a) - evB = _np.linalg.eigvals(b) - _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), - return_pairs=True) # just to get pairing - mlPl = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) - return 1.0 - mlPl / float(d2) + + if is_unitary is None: + is_unitary = _np.allclose(_np.eye(d2), b @ b.T.conj(), atol=tol, rtol=tol) + if is_tp is None: + is_tp = _tools.is_trace_preserving(a, mx_basis, tol) and _tools.is_trace_preserving(b, mx_basis, tol) + + if is_unitary and is_tp: + evA = _np.linalg.eigvals(a) + evB = _np.linalg.eigvals(b) + _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), + return_pairs=True) # just to get pairing + fid = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) / d2 + else: + Ja = _tools.fast_jamiolkowski_iso_std(a, mx_basis) + Jb = _tools.fast_jamiolkowski_iso_std(b, mx_basis) + from pygsti.tools.optools import eigenvalue_fidelity + fid = eigenvalue_fidelity(Ja, Jb, gauge_invariant=True) + return 1.0 - fid + Eigenvalue_entanglement_infidelity = _modf.opsfn_factory(eigenvalue_entanglement_infidelity) diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index e9c51c659..233a0a260 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -90,6 +90,16 @@ def is_hermitian(mx, tol=1e-9): return _np.all(_np.abs(mx - mx.T.conj()) <= tol) +def assert_hermitian(mat : _np.ndarray, tol: _np.floating) -> None: + hermiticity_error = _np.abs(mat - mat.T.conj()) + if _np.any(hermiticity_error > tol): + message = f""" + Input matrix 'mat' is not Hermitian, up to tolerance {tol}. + The absolute values of entries in (mat - mat^H) are \n{hermiticity_error}. + """ + raise ValueError(message) + + def is_pos_def(mx, tol=1e-9, attempt_cholesky=False): """ Test whether mx is a positive-definite matrix. diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 51a1c61a2..ae4f4bbbb 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -184,6 +184,67 @@ def psd_square_root(mat): return f +def eigenvalue_fidelity(x: _np.ndarray, y: _np.ndarray, gauge_invariant:bool=True) -> _np.floating: + """ + For positive semidefinite matrices (x, y), define M(x,y) = √(√x · y · √x). + + The fidelity of (x, y) is defined as F(x,y) = Tr(M(x,y))^2. + + Let v(x) and v(y) denote vectors containing the eigenvalues of x and y, + sorted in decreasing order. + + If gauge_invariant is True, then this function returns ⟨√v(x), √v(y)⟩^2, + which is always an upper bound for F(x,y). + + If gauge_invariant is False, then this function uses a heuristic to "match" + eigenvalues of x and y based on similarity of their corresponding eigenvectors. + + Proof + ----- + Let A = √x·√y, so that M(x,y) = √(A A^†). The fact that M(x,y) is a square-root + of a Gram matrix of A implies (via the polar decomposition) the existence of a + unitary matrix U where A = M(x,y)U. This lets us write fidelity as + + F(x,y) = Tr( √x · √y U^† )^2. (1) + + Next, let s(z) denote the vector of singular values of a matrix z, in descending + order. Assuming the claim from Problem III.6.12 of Bhatia's Matrix Analysis, we + have the inequality + + Tr( √x · √y U^† ) ≤ ⟨ s(√x), s(√y U^†) ⟩. (2) + + Using the fact that the singular values of √x are the same as its eigenvalues, + and the singular values of √y U^† are the eigenvalues of √y, we have + + ⟨ s(√x), s(√y U^†) ⟩ = ⟨ √v(x), √v(y) ⟩. (3) + + Combining (1), (2), and (3) proves our claim. + """ + tol = _np.finfo(x.dtype).eps ** 0.75 + _mt.assert_hermitian(x, tol) + _mt.assert_hermitian(y, tol) + if gauge_invariant: + valsX = _np.sort(_spl.eigvalsh(x)) + valsY = _np.sort(_spl.eigvalsh(y)) + pairs = [(i,i) for i in range(valsX.size)] + else: + from pygsti.tools import minweight_match + valsX, vecsX = _spl.eigh(x) + valsY, vecsY = _spl.eigh(y) + dissimilarity = lambda vec_x, vec_y : abs(1 - abs(vec_x @ vec_y)) + _, pairs = minweight_match(vecsX.T.conj(), vecsY.T.conj(), dissimilarity, return_pairs=True) + ind_x, ind_y = zip(*pairs) + arg_x = _np.maximum(valsX[list(ind_x)], 0) + arg_y = _np.maximum(valsY[list(ind_y)], 0) + f = (arg_x**0.5 @ arg_y**0.5)**2 + return f # type: ignore + + +def eigenvalue_infidelity(a, b, gauge_invariant=True) -> _np.floating: + return 1 - eigenvalue_fidelity(a, b, gauge_invariant) + + + def frobeniusdist(a, b) -> _np.floating[Any]: """ Returns the frobenius distance between arrays: ||a - b||_Fro. @@ -377,6 +438,26 @@ def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J return tracedist(JA, JB) +def is_trace_preserving(a: _np.ndarray, mx_basis: BasisLike='pp', tol=1e-8) -> bool: + dim = int(_np.round( a.size**0.5 )) + udim = int(_np.round( dim**0.5 )) + basis = _Basis.cast(mx_basis, dim=dim) + if basis.first_element_is_identity: + checkone = _np.isclose(a[0,0], 1.0, rtol=tol, atol=tol) + checktwo = _np.allclose(a[1:], 0.0, atol=tol) + return checkone and checktwo + # else, check that the adjoint of a is unital. + I_mat = _np.eye(udim) + I_vec = _bt.stdmx_to_vec(I_mat, basis) + if _np.isrealobj(a): + expect_I_vec = a.T @ I_vec + else: + expect_I_vec = a.T.conj() @ I_vec + atol = tol * udim + check = float(_np.linalg.norm(I_vec - expect_I_vec)) <= atol + return check + + def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. @@ -439,10 +520,7 @@ def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary #if the tp flag isn't set we'll calculate whether it is true here if is_tp is None: - def is_tp_fn(x): - return _np.isclose(x[0, 0], 1.0) and _np.allclose(x[0, 1:d2], 0) - - is_tp= (is_tp_fn(a) and is_tp_fn(b)) + is_tp = is_trace_preserving(a, mx_basis) and is_trace_preserving(b, mx_basis) #if the unitary flag isn't set we'll calculate whether it is true here if is_unitary is None: diff --git a/test/unit/algorithms/test_gaugeopt_correctness.py b/test/unit/algorithms/test_gaugeopt_correctness.py index 824517554..5942d9050 100644 --- a/test/unit/algorithms/test_gaugeopt_correctness.py +++ b/test/unit/algorithms/test_gaugeopt_correctness.py @@ -270,7 +270,7 @@ def test_lbfgs_frodist(self): def test_lbfgs_tracedist(self): times = [] - for seed in range(1, 3): + for seed in [1]: dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='tracedist') times.append(dt) print(f'L-BFGS GaugeOpt over UnitaryGaugeGroup w.r.t. trace dist: {times}.')