diff --git a/cpp/iterative-inl.h b/cpp/iterative-inl.h
index e1b0d007c..cd5bef43a 100644
--- a/cpp/iterative-inl.h
+++ b/cpp/iterative-inl.h
@@ -21,7 +21,7 @@ namespace gtsam {
 
 		// Start with g0 = A'*(A*x0-b), d0 = - g0
 		// i.e., first step is in direction of negative gradient
-		V g = gradient(Ab, x);
+		V g = Ab.gradient(x);
 		V d = -g;
 		double dotg0 = dot(g, g), prev_dotg = dotg0;
 		double threshold = epsilon * epsilon * dotg0;
@@ -31,7 +31,7 @@ namespace gtsam {
 				<< threshold << endl;
 
 		// loop maxIterations times
-		for (size_t k = 0; k < maxIterations; k++) {
+		for (size_t k = 1;; k++) {
 
 			// calculate optimal step-size
 			E Ad = Ab * d;
@@ -39,6 +39,7 @@ namespace gtsam {
 
 			// do step in new search direction
 			x = x + alpha * d;
+			if (k==maxIterations) break;
 
 			// update gradient
 			g = g + alpha * (Ab ^ Ad);
diff --git a/cpp/iterative.cpp b/cpp/iterative.cpp
index 18a279ad3..ec2f8ef20 100644
--- a/cpp/iterative.cpp
+++ b/cpp/iterative.cpp
@@ -13,26 +13,6 @@ using namespace std;
 namespace gtsam {
 
 	/* ************************************************************************* */
-
-	/** gradient of objective function 0.5*|Ax-b|^2 at x = A'*(Ax-b) */
-	Vector gradient(const System& Ab, const Vector& x) {
-		const Matrix& A = Ab.first;
-		const Vector& b = Ab.second;
-		return A ^ (A * x - b);
-	}
-
-	/** Apply operator A */
-	Vector operator*(const System& Ab, const Vector& x) {
-		const Matrix& A = Ab.first;
-		return A * x;
-	}
-
-	/** Apply operator A^T */
-	Vector operator^(const System& Ab, const Vector& x) {
-		const Matrix& A = Ab.first;
-		return A ^ x;
-	}
-
 	Vector steepestDescent(const System& Ab, const Vector& x, bool verbose,
 			double epsilon, size_t maxIterations) {
 		return conjugateGradients<System, Vector, Vector> (Ab, x, verbose, epsilon,
@@ -48,23 +28,19 @@ namespace gtsam {
 	/* ************************************************************************* */
 	Vector steepestDescent(const Matrix& A, const Vector& b, const Vector& x,
 			bool verbose, double epsilon, size_t maxIterations) {
-		System Ab = make_pair(A, b);
+		System Ab(A, b);
 		return conjugateGradients<System, Vector, Vector> (Ab, x, verbose, epsilon,
 				maxIterations, true);
 	}
 
 	Vector conjugateGradientDescent(const Matrix& A, const Vector& b,
 			const Vector& x, bool verbose, double epsilon, size_t maxIterations) {
-		System Ab = make_pair(A, b);
+		System Ab(A, b);
 		return conjugateGradients<System, Vector, Vector> (Ab, x, verbose, epsilon,
 				maxIterations);
 	}
 
 	/* ************************************************************************* */
-	VectorConfig gradient(const GaussianFactorGraph& fg, const VectorConfig& x) {
-		return fg.gradient(x);
-	}
-
 	VectorConfig steepestDescent(const GaussianFactorGraph& fg,
 			const VectorConfig& x, bool verbose, double epsilon, size_t maxIterations) {
 		return conjugateGradients<GaussianFactorGraph, VectorConfig, Errors> (fg,
diff --git a/cpp/iterative.h b/cpp/iterative.h
index 46a5046a6..09768a148 100644
--- a/cpp/iterative.h
+++ b/cpp/iterative.h
@@ -11,8 +11,34 @@ namespace gtsam {
 	class GaussianFactorGraph;
 	class VectorConfig;
 
-	/** typedef for combined system |Ax-b|^2 */
-	typedef std::pair<Matrix, Vector> System;
+	/** combined system |Ax-b_|^2 */
+	class System {
+	private:
+		const Matrix& A_;
+		const Vector& b_;
+
+	public:
+
+		System(const Matrix& A, const Vector& b) :
+			A_(A), b_(b) {
+		}
+
+		/** gradient of objective function 0.5*|Ax-b_|^2 at x = A_'*(Ax-b_) */
+		Vector gradient(const Vector& x) const {
+			return A_ ^ (A_ * x - b_);
+		}
+
+		/** Apply operator A_ */
+		inline Vector operator*(const Vector& x) const {
+			return A_ * x;
+		}
+
+		/** Apply operator A_^T */
+		inline Vector operator^(const Vector& e) const {
+			return A_ ^ e;
+		}
+
+	};
 
 	/**
 	 * Method of conjugate gradients (CG) template
diff --git a/cpp/testIterative.cpp b/cpp/testIterative.cpp
index 42bfe52fb..f13de3aeb 100644
--- a/cpp/testIterative.cpp
+++ b/cpp/testIterative.cpp
@@ -51,7 +51,7 @@ TEST( Iterative, conjugateGradientDescent )
 	Vector expectedX = Vector_(6, -0.1, 0.1, -0.1, -0.1, 0.1, -0.2);
 
 	// Do conjugate gradient descent, System version
-	System Ab = make_pair(A, b);
+	System Ab(A, b);
 	Vector actualX = conjugateGradientDescent(Ab, x0);
 	CHECK(assert_equal(expectedX,actualX,1e-9));