diff --git a/gtsam/linear/ConjugateGradientSolver.h b/gtsam/linear/ConjugateGradientSolver.h
index 4d4623d05..ddb6e26bb 100644
--- a/gtsam/linear/ConjugateGradientSolver.h
+++ b/gtsam/linear/ConjugateGradientSolver.h
@@ -90,8 +90,8 @@ public:
 /*********************************************************************************************/
 /*
  * A template of linear preconditioned conjugate gradient method.
- * System class should support residual(v, g), multiply(v,Av), leftPrecondition(v, S^{-t}v,
- * rightPrecondition(v, S^{-1}v), scal(alpha,v), dot(v,v), axpy(alpha,x,y)
+ * 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.
  */
 template <class S, class V>
@@ -124,19 +124,19 @@ 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 = 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);      /* d = 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.leftPrecondition(q1, q2);                      /* q2 = S^{-T} * q1 */
+    system.leftPrecondition(q1, q2);                      /* q2 = L^{-1} * q1 */
     system.axpy(-alpha, q2, residual);                    /* residual -= alpha * q2 */
     prevGamma = currentGamma;
     currentGamma = system.dot(residual, residual);        /* gamma = |residual|^2 */
     beta = currentGamma / prevGamma;
-    system.rightPrecondition(residual, q1);               /* q1 = S^{-1} residual */
+    system.rightPrecondition(residual, q1);               /* q1 = L^{-T} residual */
     system.scal(beta, direction);
     system.axpy(1.0, q1, direction);                      /* direction = q1 + beta * direction */