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) # == 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)
 
     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)  # = 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)
 
     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)  # = 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)  # = 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,9 @@ 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
+            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 +1239,13 @@ 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
+            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 +2023,13 @@ 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
+                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,9 @@ 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))
+            num = _np.vdot(BiBj, opMxInStdBasis)
+            den = _np.linalg.norm(BiBj) ** 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):