diff --git a/gtsam/linear/GaussianFactor.h b/gtsam/linear/GaussianFactor.h
index c02f59b89..92bc246e1 100644
--- a/gtsam/linear/GaussianFactor.h
+++ b/gtsam/linear/GaussianFactor.h
@@ -107,7 +107,7 @@ namespace gtsam {
     virtual GaussianFactor::shared_ptr negate() const = 0;
 
     /// y += alpha * A'*A*x
-    virtual void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y) = 0;
+    virtual void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y) const = 0;
 
     /// A'*b for Jacobian, eta for Hessian
     virtual VectorValues gradientAtZero() const = 0;
diff --git a/gtsam/linear/HessianFactor.cpp b/gtsam/linear/HessianFactor.cpp
index 5ee9b8bb5..970051faf 100644
--- a/gtsam/linear/HessianFactor.cpp
+++ b/gtsam/linear/HessianFactor.cpp
@@ -501,7 +501,7 @@ GaussianFactor::shared_ptr HessianFactor::negate() const
 
 /* ************************************************************************* */
 void HessianFactor::multiplyHessianAdd(double alpha, const VectorValues& x,
-    VectorValues& y) {
+    VectorValues& y) const {
 
   for (size_t i = 0; i < size(); ++i) {
     pair<VectorValues::iterator, bool> it = y.tryInsert(keys_[i], Vector());
diff --git a/gtsam/linear/HessianFactor.h b/gtsam/linear/HessianFactor.h
index 6b524c3d1..39c4009d2 100644
--- a/gtsam/linear/HessianFactor.h
+++ b/gtsam/linear/HessianFactor.h
@@ -357,7 +357,7 @@ namespace gtsam {
     void updateATA(const HessianFactor& update, const Scatter& scatter);
 
     /** y += alpha * A'*A*x */
-    void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y);
+    void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y) const;
 
     /// eta for Hessian
     VectorValues gradientAtZero() const;
diff --git a/gtsam/linear/JacobianFactor.cpp b/gtsam/linear/JacobianFactor.cpp
index d0512f794..bd06c5f8d 100644
--- a/gtsam/linear/JacobianFactor.cpp
+++ b/gtsam/linear/JacobianFactor.cpp
@@ -434,7 +434,7 @@ namespace gtsam {
 
   /* ************************************************************************* */
   void JacobianFactor::multiplyHessianAdd(double alpha, const VectorValues& x,
-      VectorValues& y) {
+      VectorValues& y) const {
     Vector Ax = (*this)*x;
     transposeMultiplyAdd(alpha,Ax,y);
   }
diff --git a/gtsam/linear/JacobianFactor.h b/gtsam/linear/JacobianFactor.h
index 608560ba8..d6e937769 100644
--- a/gtsam/linear/JacobianFactor.h
+++ b/gtsam/linear/JacobianFactor.h
@@ -260,7 +260,7 @@ namespace gtsam {
     Vector operator*(const VectorValues& x) const;
 
     /** y += alpha * A'*A*x */
-    void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y);
+    void multiplyHessianAdd(double alpha, const VectorValues& x, VectorValues& y) const;
 
     /** x += A'*e.  If x is initially missing any values, they are created and assumed to start as
      *  zero vectors. */