Bugfixes and efficiency improvements for basis conversion, inner products, and computing Frobenius distance

pygsti/algorithms/contract.py

@@ -358,7 +358,7 @@ def _contract_to_valid_spam(model, verbosity=0):
 
     # ** assumption: only the first vector element of pauli vectors has nonzero trace
     dummyVec = _np.zeros((model.dim, 1), 'd'); dummyVec[0, 0] = 1.0
-    firstElTrace = _np.real(_tools.trace(_tools.ppvec_to_stdmx(dummyVec)))  # == sqrt(2)**nQubits
+    firstElTrace = _np.real(_np.trace(_tools.ppvec_to_stdmx(dummyVec)))  # == sqrt(2)**nQubits
     diff = 0
 
     # rhoVec must be positive semidefinite and trace = 1

pygsti/baseobjs/basis.py

@@ -547,8 +547,7 @@ def is_normalized(self):
         """
         if self.elndim == 2:
             for i, mx in enumerate(self.elements):
-                t = _np.trace(_np.dot(mx, mx))
-                t = _np.real(t)
+                t = _np.linalg.norm(mx.ravel(), 2) # == sqrt(tr(mx mx))
                 if not _np.isclose(t, 1.0): return False
             return True
         elif self.elndim == 1:

pygsti/modelmembers/operations/__init__.py

@@ -475,18 +475,16 @@ def optimize_operation(op_to_optimize, target_op):
         return
 
     from pygsti import optimize as _opt
-    from pygsti.tools import matrixtools as _mt
     assert(target_op.dim == op_to_optimize.dim)  # operations must have the same overall dimension
     targetMatrix = target_op.to_dense() if isinstance(target_op, LinearOperator) else target_op
 
     def _objective_func(param_vec):
         op_to_optimize.from_vector(param_vec)
-        return _mt.frobeniusnorm(op_to_optimize.to_dense() - targetMatrix)
+        return _np.linalg.norm((op_to_optimize.to_dense() - targetMatrix).ravel())
 
     x0 = op_to_optimize.to_vector()
     minSol = _opt.minimize(_objective_func, x0, method='BFGS', maxiter=10000, maxfev=10000,
                            tol=1e-6, callback=None)
 
     op_to_optimize.from_vector(minSol.x)
-    #print("DEBUG: optimized operation to min frobenius distance %g" %
-    #      _mt.frobeniusnorm(op_to_optimize-targetMatrix))
+    return

pygsti/modelmembers/operations/fullcptpop.py

@@ -93,7 +93,7 @@ def _set_params_from_choi_mx(self, choi_mx, truncate):
             Lmx = _np.linalg.cholesky(choi_mx)
 
         #check TP condition: that diagonal els of Lmx squared add to 1.0
-        Lmx_norm = _np.trace(_np.dot(Lmx.T.conjugate(), Lmx))  # sum of magnitude^2 of all els
+        Lmx_norm = _np.linalg.norm(Lmx.ravel())  # = sqrt(tr(Lmx' Lmx))
         assert(_np.isclose(Lmx_norm, 1.0)), "Cholesky decomp didn't preserve trace=1!"
 
         self.params = _np.empty(dim**2, 'd')

pygsti/modelmembers/states/__init__.py

@@ -426,17 +426,16 @@ def optimize_state(vec_to_optimize, target_vec):
         return
 
     from pygsti import optimize as _opt
-    from pygsti.tools import matrixtools as _mt
     assert(target_vec.dim == vec_to_optimize.dim)  # vectors must have the same overall dimension
     targetVector = target_vec.to_dense() if isinstance(target_vec, State) else target_vec
 
     def _objective_func(param_vec):
         vec_to_optimize.from_vector(param_vec)
-        return _mt.frobeniusnorm(vec_to_optimize.to_dense() - targetVector)
+        return _np.linalg.norm((vec_to_optimize.to_dense() - targetVector).ravel())
 
     x0 = vec_to_optimize.to_vector()
     minSol = _opt.minimize(_objective_func, x0, method='BFGS', maxiter=10000, maxfev=10000,
                            tol=1e-6, callback=None)
 
     vec_to_optimize.from_vector(minSol.x)
-    #print("DEBUG: optimized vector to min frobenius distance %g" % _mt.frobeniusnorm(vec_to_optimize-targetVector))
+    return

pygsti/modelmembers/states/cptpstate.py

@@ -150,7 +150,7 @@ def _set_params_from_vector(self, vector, truncate):
             Lmx = _np.linalg.cholesky(density_mx)
 
         #check TP condition: that diagonal els of Lmx squared add to 1.0
-        Lmx_norm = _np.trace(_np.dot(Lmx.T.conjugate(), Lmx))  # sum of magnitude^2 of all els
+        Lmx_norm = _np.linalg.norm(Lmx.ravel())  # = sqrt(tr(Lmx' Lmx))
         assert(_np.isclose(Lmx_norm, 1.0)), \
             "Cholesky decomp didn't preserve trace=1!"
 
@@ -180,7 +180,7 @@ def _construct_vector(self):
             for j in range(i):
                 self.Lmx[i, j] = (self.params[i * dmDim + j] + 1j * self.params[j * dmDim + i]) / paramNorm
 
-        Lmx_norm = _np.trace(_np.dot(self.Lmx.T.conjugate(), self.Lmx))  # sum of magnitude^2 of all els
+        Lmx_norm = _np.linalg.norm(self.Lmx.ravel())  # = sqrt(tr(Lmx' Lmx))
         assert(_np.isclose(Lmx_norm, 1.0)), "Violated trace=1 condition!"
 
         #The (complex, Hermitian) density matrix is build by
@@ -192,7 +192,7 @@ def _construct_vector(self):
         # write density matrix in given basis: = sum_i alpha_i B_i
         # ASSUME that basis is orthogonal, i.e. Tr(Bi^dag*Bj) = delta_ij
         basis_mxs = _np.rollaxis(self.basis_mxs, 2)  # shape (dmDim, dmDim, len(vec))
-        vec = _np.array([_np.trace(_np.dot(M.T.conjugate(), density_mx)) for M in basis_mxs])
+        vec = _np.array([_np.vdot(M, density_mx) for M in basis_mxs])
 
         #for now, assume Liouville vector should always be real (TODO: add 'real' flag later?)
         assert(_np.linalg.norm(_np.imag(vec)) < IMAG_TOL)

pygsti/report/reportables.py

@@ -537,8 +537,13 @@ def evaluate_nearby(self, nearby_model):
             JAstd = self.d * _tools.fast_jamiolkowski_iso_std(
                 nearby_model.sim.product(self.circuit), mxBasis)
             JBstd = self.d * _tools.fast_jamiolkowski_iso_std(self.B, mxBasis)
-            Jt = (JBstd - JAstd).T
-            return 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag))
+            J = JBstd - JAstd
+            # Old code. Keep for now, to make it easier to see correctness of new code.
+            #   Jt = J.T
+            #   val = 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag))
+            # New code.
+            val = 0.5 * (_np.vdot(J.real, self.W.real) + _np.vdot(J.imag, self.W.imag))
+            return val
 
     #def circuit_half_diamond_norm(model_a, model_b, circuit):
     #    A = model_a.sim.product(circuit) # "gate"
@@ -1238,12 +1243,17 @@ def evaluate_nearby(self, nearby_model):
             -------
             float
             """
-            gl = self.oplabel; mxBasis = nearby_model.basis
-            JAstd = self.d * _tools.fast_jamiolkowski_iso_std(
-                nearby_model.operations[gl].to_dense(on_space='HilbertSchmidt'), mxBasis)
+            mxBasis = nearby_model.basis
+            A = nearby_model.operations[self.oplabel].to_dense(on_space='HilbertSchmidt')
+            JAstd = self.d * _tools.fast_jamiolkowski_iso_std(A, mxBasis)
             JBstd = self.d * _tools.fast_jamiolkowski_iso_std(self.B, mxBasis)
-            Jt = (JBstd - JAstd).T
-            return 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag))
+            J = JBstd - JAstd
+            # Old code. Keep for now, to make it easier to see correctness of new code.
+            #   Jt = J.T
+            #   val = 0.5 * _np.trace(_np.dot(Jt.real, self.W.real) + _np.dot(Jt.imag, self.W.imag))
+            # New code.
+            val = 0.5 * (_np.vdot(J.real, self.W.real) + _np.vdot(J.imag, self.W.imag))
+            return val
 
     def half_diamond_norm(a, b, mx_basis):
         """
@@ -2021,11 +2031,26 @@ def error_generator_jacobian(opstr):
                 noise = first_order_noise(opstr, errOnGate, gl)
                 jac[:, i * nSuperOps + k] = [_np.vdot(errOut.flatten(), noise.flatten()) for errOut in error_superops]
 
-                #DEBUG CHECK
-                check = [_np.trace(_np.dot(
-                    _tools.jamiolkowski_iso(errOut, mxBasis, mxBasis).conj().T,
-                    _tools.jamiolkowski_iso(noise, mxBasis, mxBasis))) * 4  # for 1-qubit...
-                    for errOut in error_superops]
+                # DEBUG CHECK  # keep old code for now, to show correctness of new code.
+                #
+                # ------------------------------ original code ----------------------------
+                # check = [_np.trace(_np.dot(
+                #     _tools.jamiolkowski_iso(errOut, mxBasis, mxBasis).conj().T,
+                #     _tools.jamiolkowski_iso(noise,  mxBasis, mxBasis))) * 4  # for 1-qubit...
+                #     for errOut in error_superops]
+                # --------------------------- more readable code --------------------------
+                # check = []
+                # for errOut in error_superops:
+                #     arg1 = _tools.jamiolkowski_iso(errOut, mxBasis, mxBasis).conj().T
+                #     arg2 = _tools.jamiolkowski_iso(noise,  mxBasis, mxBasis)
+                #     check.append(_np.trace(_np.dot(arg1, arg2)) * 4)
+                # ---------------------------- efficient code -----------------------------
+                check = []
+                for errOut in error_superops:
+                    arg1 = _tools.jamiolkowski_iso(errOut, mxBasis, mxBasis)
+                    arg2 = _tools.jamiolkowski_iso(noise,  mxBasis, mxBasis)
+                    check.append(_np.vdot(arg1, arg2) * 4)
+
                 assert(_np.allclose(jac[:, i * nSuperOps + k], check))
 
         assert(_np.linalg.norm(jac.imag) < 1e-6), "error generator jacobian should be real!"

pygsti/tools/basistools.py

@@ -549,9 +549,9 @@ def stdmx_to_vec(m, basis):
     v = _np.empty((basis.size, 1))
     for i, mx in enumerate(basis.elements):
         if basis.real:
-            v[i, 0] = _np.real(_mt.trace(_np.dot(mx, m)))
+            v[i, 0] = _np.real(_np.vdot(mx, m))
         else:
-            v[i, 0] = _np.real_if_close(_mt.trace(_np.dot(mx, m)))
+            v[i, 0] = _np.real_if_close(_np.vdot(mx, m))
     return v
 
 

pygsti/tools/jamiolkowski.py

@@ -117,9 +117,14 @@ def jamiolkowski_iso(operation_mx, op_mx_basis='pp', choi_mx_basis='pp'):
     for i in range(M):
         for j in range(M):
             BiBj = _np.kron(BVec[i], _np.conjugate(BVec[j]))
-            BiBj_dag = _np.transpose(_np.conjugate(BiBj))
-            choiMx[i, j] = _mt.trace(_np.dot(opMxInStdBasis, BiBj_dag)) \
-                / _mt.trace(_np.dot(BiBj, BiBj_dag))
+            # BiBj_dag = _np.transpose(_np.conjugate(BiBj))
+            # num = _np.trace(_np.dot(opMxInStdBasis, BiBj_dag))  # original code
+            #     = _np.trace(_np.dot(BiBj_dag, opMxInStdBasis))  # cycle the trace
+            #     = _np.vdot(BiBj, opMxInStdBasis)                # efficient version
+            # den = _np.trace(_np.dot(BiBj, BiBj_dag))
+            num = _np.vdot(BiBj, opMxInStdBasis)
+            den = _np.linalg.norm(BiBj.ravel()) ** 2
+            choiMx[i, j] = num / den 
 
     # This construction results in a Jmx with trace == dim(H) = sqrt(operation_mx.shape[0])
     #  (dimension of density matrix) but we'd like a Jmx with trace == 1, so normalize:

pygsti/tools/lindbladtools.py

@@ -13,7 +13,6 @@
 import numpy as _np
 import scipy.sparse as _sps
 
-from pygsti.tools import matrixtools as _mt
 from pygsti.tools.basistools import basis_matrices
 
 

pygsti/tools/matrixtools.py

@@ -31,37 +31,6 @@
 EXPM_DEFAULT_TOL = 2**-53  # Scipy default
 
 
-def trace(m):  # memory leak in numpy causes repeated trace calls to eat up all memory --TODO: Cython this
-    """
-    The trace of a matrix, sum_i m[i,i].
-
-    A memory leak in some version of numpy can cause repeated calls to numpy's
-    trace function to eat up all available system memory, and this function
-    does not have this problem.
-
-    Parameters
-    ----------
-    m : numpy array
-        the matrix (any object that can be double-indexed)
-
-    Returns
-    -------
-    element type of m
-        The trace of m.
-    """
-    return sum([m[i, i] for i in range(m.shape[0])])
-#    with warnings.catch_warnings():
-#        warnings.filterwarnings('error')
-#        try:
-#            ret =
-#        except Warning:
-#            print "BAD trace from:\n"
-#            for i in range(M.shape[0]):
-#                print M[i,i]
-#            raise ValueError("STOP")
-#    return ret
-
-
 def is_hermitian(mx, tol=1e-9):
     """
     Test whether mx is a hermitian matrix.
@@ -125,47 +94,7 @@ def is_valid_density_mx(mx, tol=1e-9):
     bool
         True if mx is a valid density matrix, otherwise False.
     """
-    return is_hermitian(mx, tol) and is_pos_def(mx, tol) and abs(trace(mx) - 1.0) < tol
-
-
-def frobeniusnorm(ar):
-    """
-    Compute the frobenius norm of an array (or matrix),
-
-    sqrt( sum( each_element_of_a^2 ) )
-
-    Parameters
-    ----------
-    ar : numpy array
-        What to compute the frobenius norm of.  Note that ar can be any shape
-        or number of dimenions.
-
-    Returns
-    -------
-    float or complex
-        depending on the element type of ar.
-    """
-    return _np.sqrt(_np.sum(ar**2))
-
-
-def frobeniusnorm_squared(ar):
-    """
-    Compute the squared frobenius norm of an array (or matrix),
-
-    sum( each_element_of_a^2 ) )
-
-    Parameters
-    ----------
-    ar : numpy array
-        What to compute the squared frobenius norm of.  Note that ar can be any
-        shape or number of dimenions.
-
-    Returns
-    -------
-    float or complex
-        depending on the element type of ar.
-    """
-    return _np.sum(ar**2)
+    return is_hermitian(mx, tol) and is_pos_def(mx, tol) and abs(_np.trace(mx) - 1.0) < tol
 
 
 def nullspace(m, tol=1e-7):

pygsti/tools/optools.py

@@ -105,17 +105,15 @@ def fidelity(a, b):
         evals = _np.linalg.eigvals(a)
         _warnings.warn(("sqrtm(a) failure when computing fidelity - beware result. "
                         "Maybe due to rank defficiency - eigenvalues of a are: %s") % evals)
-    F = (_mt.trace(_hack_sqrtm(_np.dot(sqrtA, _np.dot(b, sqrtA)))).real)**2  # Tr( sqrt{ sqrt(a) * b * sqrt(a) } )^2
+    F = (_np.trace(_hack_sqrtm(_np.dot(sqrtA, _np.dot(b, sqrtA)))).real)**2  # Tr( sqrt{ sqrt(a) * b * sqrt(a) } )^2
     return float(F)
 
 
 def frobeniusdist(a, b):
     """
-    Returns the frobenius distance between gate or density matrices.
+    Returns the frobenius distance between arrays: ||a - b||_Fro.
 
-    This is given by :
-
-      `sqrt( sum( (a_ij-b_ij)^2 ) )`
+    This could be inlined, but we're keeping it for API consistency with other distance functions.
 
     Parameters
     ----------
@@ -130,16 +128,14 @@ def frobeniusdist(a, b):
     float
         The resulting frobenius distance.
     """
-    return _mt.frobeniusnorm(a - b)
+    return _np.linalg.norm((a - b).ravel())
 
 
 def frobeniusdist_squared(a, b):
     """
-    Returns the square of the frobenius distance between gate or density matrices.
-
-    This is given by :
+    Returns the square of the frobenius distance between arrays: (||a - b||_Fro)^2.
 
-      `sum( (A_ij-B_ij)^2 )`
+    This could be inlined, but we're keeping it for API consistency with other distance functions.
 
     Parameters
     ----------
@@ -154,7 +150,7 @@ def frobeniusdist_squared(a, b):
     float
         The resulting frobenius distance.
     """
-    return _mt.frobeniusnorm_squared(a - b)
+    return frobeniusdist(a, b)**2
 
 
 def residuals(a, b):
@@ -742,10 +738,7 @@ def unitarity(a, mx_basis="gm"):
         B = _bt.change_basis(a, mx_basis, "gm")  # everything should be able to be put in the "gm" basis
 
     unital = B[1:d**2, 1:d**2]
-    #old version
-    #u = _np.trace(_np.dot(_np.conj(_np.transpose(unital)), unital)) / (d**2 - 1)
-    #new version
-    u= _np.einsum('ij,ji->', unital.conjugate().T, unital ) / (d**2 - 1)
+    u = _np.linalg.norm(unital.ravel())**2 / (d**2 - 1)
     return u
 
 
@@ -778,7 +771,7 @@ def fidelity_upper_bound(operation_mx):
     # # gives same result:
     # closestUnitaryJmx = _np.dot(choi_evecs, _np.dot( _np.diag(new_evals), _np.linalg.inv(choi_evecs) ) )
     closestJmx = _np.kron(closestVec, _np.transpose(_np.conjugate(closestVec)))  # closest rank-1 Jmx
-    closestJmx /= _mt.trace(closestJmx)  # normalize so trace of Jmx == 1.0
+    closestJmx /= _np.trace(closestJmx)  # normalize so trace of Jmx == 1.0
 
     maxF = fidelity(choi, closestJmx)
 
@@ -791,7 +784,7 @@ def fidelity_upper_bound(operation_mx):
         #    print "DEBUG choi_evals = ",choi_evals, " iMax = ",iMax
         #    #print "DEBUG: J = \n", closestUnitaryJmx
         #    print "DEBUG: eigvals(J) = ", _np.linalg.eigvals(closestJmx)
-        #    print "DEBUG: trace(J) = ", _mt.trace(closestJmx)
+        #    print "DEBUG: trace(J) = ", _np.trace(closestJmx)
         #    print "DEBUG: maxF = %f,  maxF_direct = %f" % (maxF, maxF_direct)
         #    raise ValueError("ERROR: maxF - maxF_direct = %f" % (maxF -maxF_direct))
         assert(abs(maxF - maxF_direct) < 1e-6)