diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 5281cad13..ccbd80848 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -265,8 +265,9 @@ def opcoeffs_to_fogiv_components_array(self, op_coeffs):
         errorgen_vec = _np.zeros(self.errorgen_space_dim, 'd')
         for i, (op_label, elem_lbl) in enumerate(self.errorgen_space_op_elem_labels):
             errorgen_vec[i] += op_coeffs[op_label].get(elem_lbl, 0.0)
-        return self.errorgen_vec_to_fogi_components_array(errorgen_vec), \
-            self.errorgen_vec_to_fogv_components_array(errorgen_vec)
+        out1 = self.errorgen_vec_to_fogi_components_array(errorgen_vec)
+        out2 = self.errorgen_vec_to_fogv_components_array(errorgen_vec)
+        return out1, out2
 
     def fogi_components_array_to_errorgen_vec(self, fogi_components):
         assert(self._dependent_fogi_action == 'drop'), \
diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py
index e486b559a..8bcb20fd6 100644
--- a/pygsti/tools/matrixtools.py
+++ b/pygsti/tools/matrixtools.py
@@ -206,7 +206,6 @@ def nullspace_qr(m, tol=1e-7):
     return q[:, rank:]
 
 
-#TODO: remove the orthogonalize argument (requires changing functions that call this one)
 def nice_nullspace(m, tol=1e-7, orthogonalize=False):
     """
     Computes the nullspace of a matrix, and tries to return a "nice" basis for it.
@@ -229,21 +228,19 @@ def nice_nullspace(m, tol=1e-7, orthogonalize=False):
     -------
     An matrix of shape (M,K) whose columns contain nullspace basis vectors.
     """
-
-    #
-    # nullsp = nullspace(m, tol)
-    # dim_ker = nullsp.shape[1]
-    # _, _, p = _spl.qr(nullsp.T.conj(), mode='raw', pivoting=True)
-    # ret = nullsp @ (nullsp.T[:, p[dim_ker]]).conj()
-    #
-    ## ^ Equivalent to, but faster than the following
-    ##
-    ##     nullsp_projector = nullsp @ nullsp.T.conj()
-    ##     ret = nullsp_projector[:, p[:dim_ker]]
-    ##
-    #
-
-    ret = nullspace(m, tol)
+    nullsp = nullspace(m, tol)
+    dim_ker = nullsp.shape[1]
+    if dim_ker == 0:
+        return nullsp  # empty 0-by-N array
+    _, _, p = _spl.qr(nullsp.T.conj(), mode='raw', pivoting=True)
+    ret = nullsp @ (nullsp.T[:, p[:dim_ker]]).conj()
+    # ^ That's equivalent to, but faster than:
+    #     nullsp_projector = nullsp @ nullsp.T.conj()
+    #     _, _, p = _spl.qr(nullsp_projector mode='raw', pivoting=True)
+    #     ret = nullsp_projector[:, p[:dim_ker]]
+
+    if orthogonalize:
+        ret, _ = _spl.qr(ret, mode='economic')
     for j in range(ret.shape[1]):  # normalize columns so largest element is +1.0
         imax = _np.argmax(_np.abs(ret[:, j]))
         if abs(ret[imax, j]) > 1e-6:
diff --git a/test/unit/objects/test_fogi.py b/test/unit/objects/test_fogi.py
index 5676d853f..d55314aa2 100644
--- a/test/unit/objects/test_fogi.py
+++ b/test/unit/objects/test_fogi.py
@@ -219,8 +219,8 @@ def test_cloud_crosstalk_fogi(self):
         nprefix = mdl.num_params - nfogi  # reparameterization *prefixes* FOGI params with "unused" params
 
         self.assertEqual(nprefix, 0)  # because include_spam=True above
-        self.assertArraysAlmostEqual(mdl.fogi_errorgen_components_array(include_fogv=False, normalized_elem_gens=True),
-                                     mdl.to_vector()[nprefix:])
+        temp = mdl.fogi_errorgen_components_array(include_fogv=False, normalized_elem_gens=True)
+        self.assertArraysAlmostEqual(temp, mdl.to_vector()[nprefix:])
 
         v = mdl.to_vector()  # just test this works