diff --git a/gtsam/linear/ConjugateGradientSolver.h b/gtsam/linear/ConjugateGradientSolver.h
index ddb6e26bb..5bf793505 100644
--- a/gtsam/linear/ConjugateGradientSolver.h
+++ b/gtsam/linear/ConjugateGradientSolver.h
@@ -93,6 +93,10 @@ public:
  * System class should support residual(v, g), multiply(v,Av), scal(alpha,v), dot(v,v), axpy(alpha,x,y)
  * leftPrecondition(v, L^{-1}v, rightPrecondition(v, L^{-T}v) where preconditioner M = L*L^T
  * Note that the residual is in the preconditioned domain. Refer to Section 9.2 of Saad's book.
+ *
+ ** REFERENCES:
+ * [1] Y. Saad, "Preconditioned Iterations," in Iterative Methods for Sparse Linear Systems,
+ * 2nd ed. SIAM, 2003, ch. 9, sec. 2, pp.276-281.
  */
 template <class S, class V>
 V preconditionedConjugateGradient(const S &system, const V &initial, const ConjugateGradientParameters &parameters) {
@@ -101,8 +105,8 @@ V preconditionedConjugateGradient(const S &system, const V &initial, const Conju
   estimate = residual = direction = q1 = q2 = initial;
 
   system.residual(estimate, q1);                /* q1 = b-Ax */
-  system.leftPrecondition(q1, residual);        /* r = S^{-T} (b-Ax) */
-  system.rightPrecondition(residual, direction);/* d = S^{-1} r */
+  system.leftPrecondition(q1, residual);        /* r = L^{-1} (b-Ax) */
+  system.rightPrecondition(residual, direction);/* p = L^{-T} r */
 
   double currentGamma = system.dot(residual, residual), prevGamma, alpha, beta;
 
@@ -125,20 +129,20 @@ V preconditionedConjugateGradient(const S &system, const V &initial, const Conju
     if ( k % iReset == 0 ) {
       system.residual(estimate, q1);                      /* q1 = b-Ax */
       system.leftPrecondition(q1, residual);              /* r = L^{-1} (b-Ax) */
-      system.rightPrecondition(residual, direction);      /* d = L^{-T} r */
+      system.rightPrecondition(residual, direction);      /* p = L^{-T} r */
       currentGamma = system.dot(residual, residual);
     }
-    system.multiply(direction, q1);                       /* q1 = A d */
-    alpha = currentGamma / system.dot(direction, q1);     /* alpha = gamma / (d' A d) */
-    system.axpy(alpha, direction, estimate);              /* estimate += alpha * direction */
+    system.multiply(direction, q1);                       /* q1 = A p */
+    alpha = currentGamma / system.dot(direction, q1);     /* alpha = gamma / (p' A p) */
+    system.axpy(alpha, direction, estimate);              /* estimate += alpha * p */
     system.leftPrecondition(q1, q2);                      /* q2 = L^{-1} * q1 */
-    system.axpy(-alpha, q2, residual);                    /* residual -= alpha * q2 */
+    system.axpy(-alpha, q2, residual);                    /* r -= alpha * q2 */
     prevGamma = currentGamma;
-    currentGamma = system.dot(residual, residual);        /* gamma = |residual|^2 */
+    currentGamma = system.dot(residual, residual);        /* gamma = |r|^2 */
     beta = currentGamma / prevGamma;
-    system.rightPrecondition(residual, q1);               /* q1 = L^{-T} residual */
+    system.rightPrecondition(residual, q1);               /* q1 = L^{-T} r */
     system.scal(beta, direction);
-    system.axpy(1.0, q1, direction);                      /* direction = q1 + beta * direction */
+    system.axpy(1.0, q1, direction);                      /* p = q1 + beta * p */
 
     if (parameters.verbosity() >= ConjugateGradientParameters::ERROR )
        std::cout << "[PCG] k = " << k