diff --git a/internal/ceres/program.cc b/internal/ceres/program.cc
index 3a88b13..1d0a157 100644
--- a/internal/ceres/program.cc
+++ b/internal/ceres/program.cc
@@ -261,6 +261,24 @@
   return true;
 }
 
+Program* Program::CreateReducedProgram(vector<double*>* removed_parameter_blocks,
+                                       double* fixed_cost,
+                                       string* error) const {
+  CHECK_NOTNULL(removed_parameter_blocks);
+  CHECK_NOTNULL(fixed_cost);
+  CHECK_NOTNULL(error);
+
+  scoped_ptr<Program> reduced_program(new Program(*this));
+  if (!reduced_program->RemoveFixedBlocks(removed_parameter_blocks,
+                                          fixed_cost,
+                                          error)) {
+    return NULL;
+  }
+
+  reduced_program->SetParameterOffsetsAndIndex();
+  return reduced_program.release();
+}
+
 bool Program::RemoveFixedBlocks(vector<double*>* removed_parameter_blocks,
                                 double* fixed_cost,
                                 string* error) {
diff --git a/internal/ceres/program.h b/internal/ceres/program.h
index 6b9190d..c7b22c4 100644
--- a/internal/ceres/program.h
+++ b/internal/ceres/program.h
@@ -129,16 +129,22 @@
   // Caller owns the result.
   TripletSparseMatrix* CreateJacobianBlockSparsityTranspose() const;
 
-  // Removes constant parameter blocks and residual blocks with no
-  // varying parameter blocks while preserving order.
+  // Create a copy of this program and removes constant parameter
+  // blocks and residual blocks with no varying parameter blocks while
+  // preserving their relative order.
   //
-  // WARNING: It is the caller's responsibility to track ownership of
-  // the parameter and residual blocks removed from the
-  // program. Generally speaking this method should only be called on
-  // a copy of an existing program.
-  bool RemoveFixedBlocks(vector<double*>* removed_parameter_blocks,
-                         double* fixed_cost,
-                         string* message);
+  // removed_parameter_blocks on exit will contain the list of
+  // parameter blocks that were removed.
+  //
+  // fixed_cost will be equal to the sum of the costs of the residual
+  // blocks that were removed.
+  //
+  // If there was a problem, then the function will return a NULL
+  // pointer and error will contain a human readable description of
+  // the problem.
+  Program* CreateReducedProgram(vector<double*>* removed_parameter_blocks,
+                                double* fixed_cost,
+                                string* error) const;
 
   // See problem.h for what these do.
   int NumParameterBlocks() const;
@@ -157,6 +163,21 @@
   string ToString() const;
 
  private:
+  // Remove constant parameter blocks and residual blocks with no
+  // varying parameter blocks while preserving their relative order.
+  //
+  // removed_parameter_blocks on exit will contain the list of
+  // parameter blocks that were removed.
+  //
+  // fixed_cost will be equal to the sum of the costs of the residual
+  // blocks that were removed.
+  //
+  // If there was a problem, then the function will return false and
+  // error will contain a human readable description of the problem.
+  bool RemoveFixedBlocks(vector<double*>* removed_parameter_blocks,
+                         double* fixed_cost,
+                         string* message);
+
   // The Program does not own the ParameterBlock or ResidualBlock objects.
   vector<ParameterBlock*> parameter_blocks_;
   vector<ResidualBlock*> residual_blocks_;
diff --git a/internal/ceres/program_test.cc b/internal/ceres/program_test.cc
index f273872..10bfa12 100644
--- a/internal/ceres/program_test.cc
+++ b/internal/ceres/program_test.cc
@@ -88,15 +88,18 @@
   problem.AddResidualBlock(new BinaryCostFunction(), NULL, &x, &y);
   problem.AddResidualBlock(new TernaryCostFunction(), NULL, &x, &y, &z);
 
-  Program program(problem.program());
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
-  EXPECT_EQ(program.NumParameterBlocks(), 3);
-  EXPECT_EQ(program.NumResidualBlocks(), 3);
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
+
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 3);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 3);
   EXPECT_EQ(removed_parameter_blocks.size(), 0);
   EXPECT_EQ(fixed_cost, 0.0);
 }
@@ -109,15 +112,17 @@
   problem.AddResidualBlock(new UnaryCostFunction(), NULL, &x);
   problem.SetParameterBlockConstant(&x);
 
-  Program program(problem.program());
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
-  EXPECT_EQ(program.NumParameterBlocks(), 0);
-  EXPECT_EQ(program.NumResidualBlocks(), 0);
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 0);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 0);
   EXPECT_EQ(removed_parameter_blocks.size(), 1);
   EXPECT_EQ(removed_parameter_blocks[0], &x);
   EXPECT_EQ(fixed_cost, 9.0);
@@ -134,15 +139,17 @@
   problem.AddParameterBlock(&y, 1);
   problem.AddParameterBlock(&z, 1);
 
-  Program program(problem.program());
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
-  EXPECT_EQ(program.NumParameterBlocks(), 0);
-  EXPECT_EQ(program.NumResidualBlocks(), 0);
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 0);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 0);
   EXPECT_EQ(removed_parameter_blocks.size(), 3);
   EXPECT_EQ(fixed_cost, 0.0);
 }
@@ -161,15 +168,17 @@
   problem.AddResidualBlock(new BinaryCostFunction(), NULL, &x, &y);
   problem.SetParameterBlockConstant(&x);
 
-  Program program(problem.program());
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
-  EXPECT_EQ(program.NumParameterBlocks(), 1);
-  EXPECT_EQ(program.NumResidualBlocks(), 1);
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 1);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 1);
 }
 
 TEST(Program, RemoveFixedBlocksNumEliminateBlocks) {
@@ -186,15 +195,17 @@
   problem.AddResidualBlock(new BinaryCostFunction(), NULL, &x, &y);
   problem.SetParameterBlockConstant(&x);
 
-  Program program(problem.program());
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
-  EXPECT_EQ(program.NumParameterBlocks(), 2);
-  EXPECT_EQ(program.NumResidualBlocks(), 2);
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 2);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 2);
 }
 
 TEST(Program, RemoveFixedBlocksFixedCost) {
@@ -211,27 +222,29 @@
   problem.AddResidualBlock(new BinaryCostFunction(), NULL, &x, &y);
   problem.SetParameterBlockConstant(&x);
 
-  Program program(problem.program());
-
-  double expected_fixed_cost;
-  ResidualBlock *expected_removed_block = program.residual_blocks()[0];
+  ResidualBlock *expected_removed_block = problem.program().residual_blocks()[0];
   scoped_array<double> scratch(
       new double[expected_removed_block->NumScratchDoublesForEvaluate()]);
+  double expected_fixed_cost;
   expected_removed_block->Evaluate(true,
                                    &expected_fixed_cost,
                                    NULL,
                                    NULL,
                                    scratch.get());
 
+
   vector<double*> removed_parameter_blocks;
   double fixed_cost = 0.0;
   string message;
-  EXPECT_TRUE(program.RemoveFixedBlocks(&removed_parameter_blocks,
-                                        &fixed_cost,
-                                        &message));
+  scoped_ptr<Program> reduced_program(
+      CHECK_NOTNULL(problem
+                    .program()
+                    .CreateReducedProgram(&removed_parameter_blocks,
+                                          &fixed_cost,
+                                          &message)));
 
-  EXPECT_EQ(program.NumParameterBlocks(), 2);
-  EXPECT_EQ(program.NumResidualBlocks(), 2);
+  EXPECT_EQ(reduced_program->NumParameterBlocks(), 2);
+  EXPECT_EQ(reduced_program->NumResidualBlocks(), 2);
   EXPECT_DOUBLE_EQ(fixed_cost, expected_fixed_cost);
 }
 
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index 11398b1..dfdfecb 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -664,42 +664,44 @@
                                           string* error) {
   CHECK_NOTNULL(options->linear_solver_ordering.get());
   Program* original_program = problem_impl->mutable_program();
-  scoped_ptr<Program> transformed_program(new Program(*original_program));
+
+  vector<double*> removed_parameter_blocks;
+  scoped_ptr<Program> reduced_program(
+      original_program->CreateReducedProgram(&removed_parameter_blocks,
+                                             fixed_cost,
+                                             error));
+  if (reduced_program.get() == NULL) {
+    return NULL;
+  }
+
+  VLOG(2) << "Reduced problem: "
+          << reduced_program->NumParameterBlocks()
+          << " parameter blocks, "
+          << reduced_program->NumParameters()
+          << " parameters,  "
+          << reduced_program->NumResidualBlocks()
+          << " residual blocks, "
+          << reduced_program->NumResiduals()
+          << " residuals.";
+
+  if (reduced_program->NumParameterBlocks() == 0) {
+    LOG(WARNING) << "No varying parameter blocks to optimize; "
+                 << "bailing early.";
+    return reduced_program.release();
+  }
 
   ParameterBlockOrdering* linear_solver_ordering =
       options->linear_solver_ordering.get();
   const int min_group_id =
       linear_solver_ordering->group_to_elements().begin()->first;
-  vector<double*> removed_parameter_blocks;
-  if (!transformed_program->RemoveFixedBlocks(&removed_parameter_blocks,
-                                              fixed_cost,
-                                              error)) {
-    return NULL;
-  }
-
   linear_solver_ordering->Remove(removed_parameter_blocks);
+
   ParameterBlockOrdering* inner_iteration_ordering =
       options->inner_iteration_ordering.get();
   if (inner_iteration_ordering != NULL) {
     inner_iteration_ordering->Remove(removed_parameter_blocks);
   }
 
-  VLOG(2) << "Reduced problem: "
-          << transformed_program->NumParameterBlocks()
-          << " parameter blocks, "
-          << transformed_program->NumParameters()
-          << " parameters,  "
-          << transformed_program->NumResidualBlocks()
-          << " residual blocks, "
-          << transformed_program->NumResiduals()
-          << " residuals.";
-
-  if (transformed_program->NumParameterBlocks() == 0) {
-    LOG(WARNING) << "No varying parameter blocks to optimize; "
-                 << "bailing early.";
-    return transformed_program.release();
-  }
-
   if (IsSchurType(options->linear_solver_type) &&
       linear_solver_ordering->GroupSize(min_group_id) == 0) {
     // If the user requested the use of a Schur type solver, and
@@ -730,11 +732,11 @@
             options->sparse_linear_algebra_library_type,
             problem_impl->parameter_map(),
             linear_solver_ordering,
-            transformed_program.get(),
+            reduced_program.get(),
             error)) {
       return NULL;
     }
-    return transformed_program.release();
+    return reduced_program.release();
   }
 
   if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY &&
@@ -742,16 +744,16 @@
     if (!ReorderProgramForSparseNormalCholesky(
             options->sparse_linear_algebra_library_type,
             linear_solver_ordering,
-            transformed_program.get(),
+            reduced_program.get(),
             error)) {
       return NULL;
     }
 
-    return transformed_program.release();
+    return reduced_program.release();
   }
 
-  transformed_program->SetParameterOffsetsAndIndex();
-  return transformed_program.release();
+  reduced_program->SetParameterOffsetsAndIndex();
+  return reduced_program.release();
 }
 
 LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options,