diff --git a/gext/fitting.py b/gext/fitting.py
index 41fc250b5fa617723da8b9cf25259bd948b20635..dac86d64782226a53ecd3d21a67d95e92d65a383 100644
--- a/gext/fitting.py
+++ b/gext/fitting.py
@@ -76,22 +76,13 @@ class DiffFitting(AbstractFitting):
         if self.options["regularization"] > 0.0:
             a += np.identity(len(b))*self.options["regularization"]
         coefficients = np.linalg.solve(a, b)
-        return np.array(coefficients, dtype=np.float64)
 
-    def linear_combination(self, vectors: List[np.ndarray],
-            coefficients: np. ndarray) -> np.ndarray:
-        """Given a set of vectors (or matrices) and the corresponding
-        coefficients, build their linear combination."""
-        if len(vectors) == 1:
-            raise ValueError("DiffFitting does not work for one vector")
-        result = np.zeros(vectors[0].shape, dtype=np.float64)
-        diff_vectors = []
-        for i in range(1,len(vectors)):
-              diff_vectors.append(vectors[i-1]-vectors[-1])
-        for coeff, vector in zip(coefficients, diff_vectors):
-            result += vector*coeff
-        result=result+vectors[-1]
-        return result
+        # Convert diff coefficients to normal coefficients (like in
+        # least square fitting)
+        coefficients = np.concatenate((coefficients, [1.0 - np.sum(coefficients)]))
+
+        return coefficients
+
 
 class LeastSquare(AbstractFitting):
 
@@ -158,6 +149,8 @@ class QuasiTimeReversible(AbstractFitting):
             a += np.identity(len(b))*self.options["regularization"]
         coefficients = np.linalg.solve(a, b)
 
+        # Convert quasi time reversible coefficients to normal
+        # coefficients (like in least square fitting)
         if q == 1:
             full_coefficients = np.concatenate(([-1.0], coefficients))
         elif q%2 == 0: