diff --git a/gtsam/linear/JacobianFactor.cpp b/gtsam/linear/JacobianFactor.cpp
index 2bd985048..6130a7b06 100644
--- a/gtsam/linear/JacobianFactor.cpp
+++ b/gtsam/linear/JacobianFactor.cpp
@@ -657,7 +657,6 @@ std::pair<boost::shared_ptr<GaussianConditional>,
   // Combine and sort variable blocks in elimination order
   JacobianFactor::shared_ptr jointFactor;
   try {
-    cout << "JacobianFactor make_shared" << endl;
     jointFactor = boost::make_shared<JacobianFactor>(factors, keys);
   } catch (std::invalid_argument&) {
     throw InvalidDenseElimination(
@@ -665,8 +664,6 @@ std::pair<boost::shared_ptr<GaussianConditional>,
             "involved in the provided factors.");
   }
 
-  jointFactor->print("JointFactor0:");
-
   // Do dense elimination
   if (jointFactor->model_)
     jointFactor->model_ = jointFactor->model_->QR(jointFactor->Ab_.matrix());
@@ -676,13 +673,9 @@ std::pair<boost::shared_ptr<GaussianConditional>,
   // Zero below the diagonal
   jointFactor->Ab_.matrix().triangularView<Eigen::StrictlyLower>().setZero();
 
-  factors.print("Factors to eliminate: ");
-  jointFactor->print("JointFactor1:");
-
   // Split elimination result into conditional and remaining factor
   GaussianConditional::shared_ptr conditional = jointFactor->splitConditional(
       keys.size());
-  cout << "JacobianFactor split conditoinal ok!" << endl;
 
   return make_pair(conditional, jointFactor);
 }
diff --git a/gtsam/linear/NoiseModel.cpp b/gtsam/linear/NoiseModel.cpp
index 989d04c66..1a0f170f4 100644
--- a/gtsam/linear/NoiseModel.cpp
+++ b/gtsam/linear/NoiseModel.cpp
@@ -415,7 +415,7 @@ Constrained::shared_ptr Constrained::unit(size_t augmentedDim) const {
 // It is Gram-Schmidt orthogonalization rather than Householder
 // Previously Diagonal::QR
 SharedDiagonal Constrained::QR(Matrix& Ab) const {
-  bool verbose = true;
+  bool verbose = false;
   if (verbose) cout << "\nStarting Constrained::QR" << endl;
 
   // get size(A) and maxRank
@@ -428,34 +428,29 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
 
   Vector pseudo(m); // allocate storage for pseudo-inverse
   Vector invsigmas = reciprocal(sigmas_);
-  gtsam::print(invsigmas, "invsigmas: ");
   Vector weights = emul(invsigmas,invsigmas); // calculate weights once
 
   // Denote weights of inequality constraints with the negative sign
   // TODO: might be slow!
   for (size_t s = 0; s<sigmas_.size(); ++s)
     weights[s] = (sigmas_[s]<0)?-weights[s]:weights[s];
-  gtsam::print(weights, "weights: ");
+  if (verbose) gtsam::print(weights, "weights: ");
 
   // We loop over all columns, because the columns that can be eliminated
   // are not necessarily contiguous. For each one, estimate the corresponding
   // scalar variable x as d-rS, with S the separator (remaining columns).
   // Then update A and b by substituting x with d-rS, zero-ing out x's column.
-  gtsam::print(Ab, "Ab = ");
-  cout << " n = " << n << endl;
+  if (verbose) gtsam::print(Ab, "Ab = ");
   for (size_t j=0; j<n; ++j) {
-    cout << "--------------------" << endl;
-    cout << "j: " << j << endl;
+    if (verbose) cout << "j: " << j << endl;
     // extract the first column of A
     Vector a = Ab.col(j);
-    gtsam::print(a, "a = ");
 
     // Calculate weighted pseudo-inverse and corresponding precision
     gttic(constrained_QR_weightedPseudoinverse);
     double precision = weightedPseudoinverse(a, weights, pseudo);
     gttoc(constrained_QR_weightedPseudoinverse);
-    cout << "precision: " << precision << endl;
-    gtsam::print(pseudo, "pseudo: ");
+    if (verbose) gtsam::print(pseudo, "pseudo: ");
 
     // If precision is zero, no information on this column
     // This is actually not limited to constraints, could happen in Gaussian::QR
@@ -468,14 +463,12 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
     rd(j)=1.0; // put 1 on diagonal
     for (size_t j2=j+1; j2<n+1; ++j2) { // and fill in remainder with dot-products
       Vector Abj2 = Ab.col(j2);
-      gtsam::print(Abj2, "Ab.col(j2): ");
       rd(j2) = pseudo.dot(Ab.col(j2));
     }
     gttoc(constrained_QR_create_rd);
 
     // construct solution (r, d, sigma)
     Rd.push_back(boost::make_tuple(j, rd, precision));
-    gtsam::print(rd, "rd = ");
 
     // exit after rank exhausted
     if (Rd.size()>=maxRank) break;
@@ -484,16 +477,12 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
     gttic(constrained_QR_update_Ab);
     Ab.middleCols(j+1,n-j) -= a * rd.segment(j+1, n-j).transpose();
     gttoc(constrained_QR_update_Ab);
-    gtsam::print(Ab, "Updated Ab = ");
-
   }
 
-  cout << "Create storage for precisions" << endl;
   // Create storage for precisions
   Vector precisions(Rd.size());
 
   gttic(constrained_QR_write_back_into_Ab);
-  cout << "write back result" << endl;
   // Write back result in Ab, imperative as we are
   // TODO: test that is correct if a column was skipped !!!!
   size_t i = 0; // start with first row
@@ -511,8 +500,6 @@ SharedDiagonal Constrained::QR(Matrix& Ab) const {
   }
   gttoc(constrained_QR_write_back_into_Ab);
 
-  cout << "return " << (int) mixed << endl;
-  gtsam::print(precisions, "precisions:");
   // Must include mu, as the defaults might be higher, resulting in non-convergence
   return mixed ? Constrained::MixedPrecisions(mu_, precisions) : Diagonal::Precisions(precisions);
 }
diff --git a/gtsam/linear/QPSolver.cpp b/gtsam/linear/QPSolver.cpp
index eeebab7e7..58e8c119e 100644
--- a/gtsam/linear/QPSolver.cpp
+++ b/gtsam/linear/QPSolver.cpp
@@ -210,7 +210,7 @@ bool QPSolver::updateWorkingSetInplace(GaussianFactorGraph& workingGraph,
 /* ************************************************************************* */
 boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGraph& workingGraph,
     const VectorValues& xk, const VectorValues& p) const {
-  static bool debug = true;
+  static bool debug = false;
 
   double minAlpha = 1.0;
   int closestFactorIx = -1, closestSigmaIx = -1;
@@ -228,9 +228,7 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
           Vector aj = jacobian->getA(xj).row(s);
           ajTp += aj.dot(pj);
         }
-        if (debug) {
-          cout << "s, ajTp: " << s << " " << ajTp << endl;
-        }
+        if (debug) cout << "s, ajTp: " << s << " " << ajTp << endl;
 
         // Check if  aj'*p >0. Don't care if it's not.
         if (ajTp<=0) continue;
@@ -242,15 +240,11 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
           Vector aj = jacobian->getA(xj).row(s);
           ajTx += aj.dot(xkj);
         }
-        if (debug) {
-          cout << "b[s], ajTx: " << b[s] << " " << ajTx << " " << ajTp << endl;
-        }
+        if (debug) cout << "b[s], ajTx: " << b[s] << " " << ajTx << " " << ajTp << endl;
 
         // alpha = (bj - aj'*xk) / (aj'*p)
         double alpha = (b[s] - ajTx)/ajTp;
-        if (debug) {
-          cout << "alpha: " << alpha << endl;
-        }
+        if (debug) cout << "alpha: " << alpha << endl;
 
         // We want the minimum of all those max alphas
         if (alpha < minAlpha) {
@@ -266,7 +260,7 @@ boost::tuple<double, int, int> QPSolver::computeStepSize(const GaussianFactorGra
 
 /* ************************************************************************* */
 bool QPSolver::iterateInPlace(GaussianFactorGraph& workingGraph, VectorValues& currentSolution) const {
-  static bool debug = true;
+  static bool debug = false;
   // Obtain the solution from the current working graph
   VectorValues newSolution = workingGraph.optimize();
   if (debug) newSolution.print("New solution:");
@@ -294,9 +288,7 @@ bool QPSolver::iterateInPlace(GaussianFactorGraph& workingGraph, VectorValues& c
     int factorIx, sigmaIx;
     VectorValues p = newSolution - currentSolution;
     boost::tie(alpha, factorIx, sigmaIx) = computeStepSize(workingGraph, currentSolution, p);
-    if (debug) {
-      cout << "alpha, factorIx, sigmaIx: " << alpha << " " << factorIx << " " << sigmaIx << endl;
-    }
+    if (debug) cout << "alpha, factorIx, sigmaIx: " << alpha << " " << factorIx << " " << sigmaIx << endl;
     // also add to the working set the one that complains the most
     updateWorkingSetInplace(workingGraph, factorIx, sigmaIx, 0.0);
     // step!