diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 4f754c9..5bc5382 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -1537,7 +1537,15 @@
 
 .. code-block:: c++
 
-   ProductManifold se3(new QuaternionManifold, new EuclideanManifold<3>);
+   ProductManifold<QuaternionManifold, EuclideanManifold<3>> se3{QuaternionManifold{}, EuclideanManifold<3>{}};
+
+
+In C++17, the template parameters can be left out as they are automatically
+deduced making the initialization much simpler:
+
+.. code-block:: c++
+
+   ProductManifold se3{QuaternionManifold{}, EuclideanManifold<3>{}};
 
 
 :class:`QuaternionManifold`
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index b3590e9..d42a0d7 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -306,7 +306,8 @@
   if (CERES_GET_FLAG(FLAGS_use_quaternions) &&
       CERES_GET_FLAG(FLAGS_use_manifolds)) {
     Manifold* camera_manifold =
-        new ProductManifold(new QuaternionManifold, new EuclideanManifold<6>);
+        new ProductManifold<QuaternionManifold, EuclideanManifold<6>>{
+            QuaternionManifold{}, EuclideanManifold<6>{}};
     for (int i = 0; i < bal_problem->num_cameras(); ++i) {
       problem->SetManifold(cameras + camera_block_size * i, camera_manifold);
     }
diff --git a/include/ceres/ceres.h b/include/ceres/ceres.h
index b67b5e6..c32477d 100644
--- a/include/ceres/ceres.h
+++ b/include/ceres/ceres.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2019 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -64,6 +64,7 @@
 #include "ceres/numeric_diff_options.h"
 #include "ceres/ordered_groups.h"
 #include "ceres/problem.h"
+#include "ceres/product_manifold.h"
 #include "ceres/sized_cost_function.h"
 #include "ceres/solver.h"
 #include "ceres/sphere_manifold.h"
diff --git a/include/ceres/manifold.h b/include/ceres/manifold.h
index d50de7f..24eac23 100644
--- a/include/ceres/manifold.h
+++ b/include/ceres/manifold.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2021 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -35,6 +35,7 @@
 #include <algorithm>
 #include <array>
 #include <memory>
+#include <utility>
 #include <vector>
 
 #include "ceres/internal/disable_warnings.h"
@@ -337,72 +338,6 @@
   std::vector<bool> constancy_mask_;
 };
 
-// Construct a manifold by taking the Cartesian product of a number of other
-// manifolds. This is useful, when a parameter block is the cartesian product
-// of two or more manifolds. For example the parameters of a camera consist of
-// a rotation and a translation, i.e., SO(3) x R^3.
-//
-// Example usage:
-//
-// ProductParameterization se3(new Quaternion(), new EuclideanManifold(3));
-//
-// is the manifold for a rigid transformation, where the rotation is
-// represented using a quaternion.
-class CERES_EXPORT ProductManifold final : public Manifold {
- public:
-  ProductManifold(const ProductManifold&) = delete;
-  ProductManifold& operator=(const ProductManifold&) = delete;
-  // NOTE: Do not remove the trivial destructor as this will cause linker
-  // errors in MSVC builds.
-  ~ProductManifold() override;
-
-  // NOTE: The constructor takes ownership of the input
-  // manifolds.
-  //
-  template <typename... Manifolds>
-  explicit ProductManifold(Manifolds*... manifolds)
-      : manifolds_(sizeof...(Manifolds)) {
-    constexpr int kNumManifolds = sizeof...(Manifolds);
-    static_assert(kNumManifolds >= 2,
-                  "At least two manifolds must be specified.");
-
-    using ManifoldPtr = std::unique_ptr<Manifold>;
-
-    // Wrap all raw pointers into std::unique_ptr for exception safety.
-    std::array<ManifoldPtr, kNumManifolds> manifolds_array{
-        ManifoldPtr(manifolds)...};
-
-    // Initialize internal state.
-    for (int i = 0; i < kNumManifolds; ++i) {
-      ManifoldPtr& manifold = manifolds_[i];
-      manifold = std::move(manifolds_array[i]);
-
-      buffer_size_ = (std::max)(
-          buffer_size_, manifold->TangentSize() * manifold->AmbientSize());
-      ambient_size_ += manifold->AmbientSize();
-      tangent_size_ += manifold->TangentSize();
-    }
-  }
-
-  int AmbientSize() const override;
-  int TangentSize() const override;
-
-  bool Plus(const double* x,
-            const double* delta,
-            double* x_plus_delta) const override;
-  bool PlusJacobian(const double* x, double* jacobian) const override;
-  bool Minus(const double* y,
-             const double* x,
-             double* y_minus_x) const override;
-  bool MinusJacobian(const double* x, double* jacobian) const override;
-
- private:
-  std::vector<std::unique_ptr<Manifold>> manifolds_;
-  int ambient_size_ = 0;
-  int tangent_size_ = 0;
-  int buffer_size_ = 0;
-};
-
 // Implements the manifold for a Hamilton quaternion as defined in
 // https://en.wikipedia.org/wiki/Quaternion. Quaternions are represented as
 // unit norm 4-vectors, i.e.
diff --git a/include/ceres/product_manifold.h b/include/ceres/product_manifold.h
new file mode 100644
index 0000000..bd6131f
--- /dev/null
+++ b/include/ceres/product_manifold.h
@@ -0,0 +1,292 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2022 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//         sergiu.deitsch@gmail.com (Sergiu Deitsch)
+//
+
+#ifndef CERES_PUBLIC_PRODUCT_MANIFOLD_H_
+#define CERES_PUBLIC_PRODUCT_MANIFOLD_H_
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <numeric>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/fixed_array.h"
+#include "ceres/internal/port.h"
+#include "ceres/manifold.h"
+
+namespace ceres {
+
+// Construct a manifold by taking the Cartesian product of a number of other
+// manifolds. This is useful, when a parameter block is the Cartesian product
+// of two or more manifolds. For example the parameters of a camera consist of
+// a rotation and a translation, i.e., SO(3) x R^3.
+//
+// Example usage:
+//
+// ProductManifold<QuaternionManifold, EuclideanManifold<3>>
+//      se3(QuaternionManifold(), EuclideanManifold<3>());
+//
+// is the manifold for a rigid transformation, where the rotation is
+// represented using a quaternion.
+//
+// In C++17, the template parameters can be left out as they are automatically
+// deduced making the initialization much simpler:
+//
+// ProductManifold se3{QuaternionManifold{}, EuclideanManifold<3>{}};
+//
+// The manifold implementations must be either copyable or moveable to be usable
+// in a ProductManifold.
+template <typename... Manifolds>
+class ProductManifold final : public Manifold {
+  struct InitTag {};
+
+ public:
+  template <
+      typename Manifold0,
+      typename Manifold1,
+      typename... ManifoldN,
+      std::enable_if_t<std::is_constructible<std::tuple<Manifolds...>,
+                                             Manifold0,
+                                             Manifold1,
+                                             ManifoldN...>::value>* = nullptr>
+  ProductManifold(Manifold0&& manifold0,
+                  Manifold1&& manifold1,
+                  ManifoldN&&... manifolds)
+      : ProductManifold{
+            InitTag{},  // Use InitTag to disambiguate public/private
+                        // constructors, which both use variadic arguments
+            std::forward<Manifold0>(manifold0),
+            std::forward<Manifold1>(manifold1),
+            std::forward<ManifoldN>(manifolds)...} {}
+
+  int AmbientSize() const override { return ambient_size_; }
+  int TangentSize() const override { return tangent_size_; }
+
+  bool Plus(const double* x,
+            const double* delta,
+            double* x_plus_delta) const override {
+    return PlusImpl(
+        x, delta, x_plus_delta, std::make_index_sequence<kNumManifolds>{});
+  }
+
+  bool Minus(const double* y,
+             const double* x,
+             double* y_minus_x) const override {
+    return MinusImpl(
+        y, x, y_minus_x, std::make_index_sequence<kNumManifolds>{});
+  }
+
+  bool PlusJacobian(const double* x, double* jacobian_ptr) const override {
+    MatrixRef jacobian(jacobian_ptr, AmbientSize(), TangentSize());
+    jacobian.setZero();
+    internal::FixedArray<double> buffer(buffer_size_);
+
+    return PlusJacobianImpl(
+        x, jacobian, buffer, std::make_index_sequence<kNumManifolds>{});
+  }
+
+  bool MinusJacobian(const double* x, double* jacobian_ptr) const override {
+    MatrixRef jacobian(jacobian_ptr, TangentSize(), AmbientSize());
+    jacobian.setZero();
+    internal::FixedArray<double> buffer(buffer_size_);
+
+    return MinusJacobianImpl(
+        x, jacobian, buffer, std::make_index_sequence<kNumManifolds>{});
+  }
+
+ private:
+  static constexpr std::size_t kNumManifolds = sizeof...(Manifolds);
+
+  // In the public constructor, we enforce at least two parameters. Then, the
+  // public constructor delegates the arguments to this private constructor to
+  // be able to construct-initialize the manifolds tuple. The delegation is
+  // necessary to avoid requiring the manifolds to be default constructible.
+  // InitTag is used as a token to disambiguate between both variadic parameter
+  // constructors.
+  template <typename... Args>
+  explicit ProductManifold(InitTag, Args&&... manifolds)
+      : buffer_size_{(std::max)(
+            {(manifolds.TangentSize() * manifolds.AmbientSize())...})},
+        ambient_sizes_{manifolds.AmbientSize()...},
+        tangent_sizes_{manifolds.TangentSize()...},
+        ambient_offsets_{ExclusiveScan(ambient_sizes_)},
+        tangent_offsets_{ExclusiveScan(tangent_sizes_)},
+        ambient_size_{
+            std::accumulate(ambient_sizes_.begin(), ambient_sizes_.end(), 0)},
+        tangent_size_{
+            std::accumulate(tangent_sizes_.begin(), tangent_sizes_.end(), 0)},
+        manifolds_{std::forward<Args>(manifolds)...} {}
+
+  template <std::size_t Index0, std::size_t... Indices>
+  bool PlusImpl(const double* x,
+                const double* delta,
+                double* x_plus_delta,
+                std::index_sequence<Index0, Indices...>) const {
+    if (!std::get<Index0>(manifolds_)
+             .Plus(x + ambient_offsets_[Index0],
+                   delta + tangent_offsets_[Index0],
+                   x_plus_delta + ambient_offsets_[Index0])) {
+      return false;
+    }
+
+    return PlusImpl(x, delta, x_plus_delta, std::index_sequence<Indices...>{});
+  }
+
+  static constexpr bool PlusImpl(const double* /*x*/,
+                                 const double* /*delta*/,
+                                 double* /*x_plus_delta*/,
+                                 std::index_sequence<>) noexcept {
+    return true;
+  }
+
+  template <std::size_t Index0, std::size_t... Indices>
+  bool MinusImpl(const double* y,
+                 const double* x,
+                 double* y_minus_x,
+                 std::index_sequence<Index0, Indices...>) const {
+    if (!std::get<Index0>(manifolds_)
+             .Minus(y + ambient_offsets_[Index0],
+                    x + ambient_offsets_[Index0],
+                    y_minus_x + tangent_offsets_[Index0])) {
+      return false;
+    }
+
+    return MinusImpl(y, x, y_minus_x, std::index_sequence<Indices...>{});
+  }
+
+  static constexpr bool MinusImpl(const double* /*y*/,
+                                  const double* /*x*/,
+                                  double* /*y_minus_x*/,
+                                  std::index_sequence<>) noexcept {
+    return true;
+  }
+
+  template <std::size_t Index0, std::size_t... Indices>
+  bool PlusJacobianImpl(const double* x,
+                        MatrixRef& jacobian,
+                        internal::FixedArray<double>& buffer,
+                        std::index_sequence<Index0, Indices...>) const {
+    if (!std::get<Index0>(manifolds_)
+             .PlusJacobian(x + ambient_offsets_[Index0], buffer.data())) {
+      return false;
+    }
+
+    jacobian.block(ambient_offsets_[Index0],
+                   tangent_offsets_[Index0],
+                   ambient_sizes_[Index0],
+                   tangent_sizes_[Index0]) =
+        MatrixRef(
+            buffer.data(), ambient_sizes_[Index0], tangent_sizes_[Index0]);
+
+    return PlusJacobianImpl(
+        x, jacobian, buffer, std::index_sequence<Indices...>{});
+  }
+
+  static constexpr bool PlusJacobianImpl(
+      const double* /*x*/,
+      MatrixRef& /*jacobian*/,
+      internal::FixedArray<double>& /*buffer*/,
+      std::index_sequence<>) noexcept {
+    return true;
+  }
+
+  template <std::size_t Index0, std::size_t... Indices>
+  bool MinusJacobianImpl(const double* x,
+                         MatrixRef& jacobian,
+                         internal::FixedArray<double>& buffer,
+                         std::index_sequence<Index0, Indices...>) const {
+    if (!std::get<Index0>(manifolds_)
+             .MinusJacobian(x + ambient_offsets_[Index0], buffer.data())) {
+      return false;
+    }
+
+    jacobian.block(tangent_offsets_[Index0],
+                   ambient_offsets_[Index0],
+                   tangent_sizes_[Index0],
+                   ambient_sizes_[Index0]) =
+        MatrixRef(
+            buffer.data(), tangent_sizes_[Index0], ambient_sizes_[Index0]);
+
+    return MinusJacobianImpl(
+        x, jacobian, buffer, std::index_sequence<Indices...>{});
+  }
+
+  static constexpr bool MinusJacobianImpl(
+      const double* /*x*/,
+      MatrixRef& /*jacobian*/,
+      internal::FixedArray<double>& /*buffer*/,
+      std::index_sequence<>) noexcept {
+    return true;
+  }
+
+  template <typename T, std::size_t N>
+  static std::array<T, N> ExclusiveScan(const std::array<T, N>& values) {
+    std::array<T, N> result;
+    T init = 0;
+
+    // TODO Replace by std::exclusive_scan once C++17 is available
+    for (std::size_t i = 0; i != N; ++i) {
+      result[i] = init;
+      init += values[i];
+    }
+
+    return result;
+  }
+
+  int buffer_size_;
+  std::array<int, kNumManifolds> ambient_sizes_;
+  std::array<int, kNumManifolds> tangent_sizes_;
+  std::array<int, kNumManifolds> ambient_offsets_;
+  std::array<int, kNumManifolds> tangent_offsets_;
+  int ambient_size_;
+  int tangent_size_;
+  std::tuple<Manifolds...> manifolds_;
+};
+
+#ifdef CERES_HAS_CPP17
+// C++17 deduction guide that allows the user to avoid explicitly specifying
+// the template parameters of ProductManifold. The class can instead be
+// instantiated as follows:
+//
+//   ProductManifold manifold{QuaternionManifold{}, EuclideanManifold<3>{}};
+//
+template <typename Manifold0, typename Manifold1, typename... Manifolds>
+ProductManifold(Manifold0&&, Manifold1&&, Manifolds&&...)
+    -> ProductManifold<Manifold0, Manifold1, Manifolds...>;
+#endif
+
+}  // namespace ceres
+
+#endif  // CERES_PUBLIC_PRODUCT_MANIFOLD_H_
diff --git a/internal/ceres/manifold.cc b/internal/ceres/manifold.cc
index b2ff953..f412a79 100644
--- a/internal/ceres/manifold.cc
+++ b/internal/ceres/manifold.cc
@@ -267,112 +267,6 @@
   return true;
 }
 
-// While the trivial destructor is not necessary, MSVC (1930) complains in
-// shared library builds by issuing a linker error:
-//
-// bundle_adjuster.obj : error LNK2019: unresolved external symbol
-// "__declspec(dllimport) const ceres::ProductManifold::`vftable'"
-// (__imp_??_7ProductManifold@ceres@@6B@) referenced in function "public:
-// __cdecl ceres::ProductManifold::ProductManifold<class
-// ceres::QuaternionManifold,class ceres::EuclideanManifold>(class
-// ceres::QuaternionManifold *,class ceres::EuclideanManifold *)"
-// (??$?0VQuaternionManifold@ceres@@VEuclideanManifold@1@@ProductManifold@ceres@@QEAA@PEAVQuaternionManifold@1@PEAVEuclideanManifold@1@@Z)
-//
-// Providing a trivial implementation here resolves the issue.
-ProductManifold::~ProductManifold() = default;
-
-int ProductManifold::AmbientSize() const { return ambient_size_; }
-
-int ProductManifold::TangentSize() const { return tangent_size_; }
-
-bool ProductManifold::Plus(const double* x,
-                           const double* delta,
-                           double* x_plus_delta) const {
-  int ambient_cursor = 0;
-  int tangent_cursor = 0;
-  for (const auto& m : manifolds_) {
-    if (!m->Plus(x + ambient_cursor,
-                 delta + tangent_cursor,
-                 x_plus_delta + ambient_cursor)) {
-      return false;
-    }
-    tangent_cursor += m->TangentSize();
-    ambient_cursor += m->AmbientSize();
-  }
-
-  return true;
-}
-
-bool ProductManifold::PlusJacobian(const double* x,
-                                   double* jacobian_ptr) const {
-  MatrixRef jacobian(jacobian_ptr, AmbientSize(), TangentSize());
-  jacobian.setZero();
-  internal::FixedArray<double> buffer(buffer_size_);
-
-  int ambient_cursor = 0;
-  int tangent_cursor = 0;
-  for (const auto& m : manifolds_) {
-    const int ambient_size = m->AmbientSize();
-    const int tangent_size = m->TangentSize();
-
-    if (!m->PlusJacobian(x + ambient_cursor, buffer.data())) {
-      return false;
-    }
-
-    jacobian.block(ambient_cursor, tangent_cursor, ambient_size, tangent_size) =
-        MatrixRef(buffer.data(), ambient_size, tangent_size);
-
-    ambient_cursor += ambient_size;
-    tangent_cursor += tangent_size;
-  }
-
-  return true;
-}
-
-bool ProductManifold::Minus(const double* y,
-                            const double* x,
-                            double* y_minus_x) const {
-  int ambient_cursor = 0;
-  int tangent_cursor = 0;
-  for (const auto& m : manifolds_) {
-    if (!m->Minus(y + ambient_cursor,
-                  x + ambient_cursor,
-                  y_minus_x + tangent_cursor)) {
-      return false;
-    }
-    tangent_cursor += m->TangentSize();
-    ambient_cursor += m->AmbientSize();
-  }
-
-  return true;
-}
-
-bool ProductManifold::MinusJacobian(const double* x,
-                                    double* jacobian_ptr) const {
-  MatrixRef jacobian(jacobian_ptr, TangentSize(), AmbientSize());
-  jacobian.setZero();
-  internal::FixedArray<double> buffer(buffer_size_);
-
-  int ambient_cursor = 0;
-  int tangent_cursor = 0;
-  for (const auto& m : manifolds_) {
-    const int ambient_size = m->AmbientSize();
-    const int tangent_size = m->TangentSize();
-
-    if (!m->MinusJacobian(x + ambient_cursor, buffer.data())) {
-      return false;
-    }
-
-    jacobian.block(tangent_cursor, ambient_cursor, tangent_size, ambient_size) =
-        MatrixRef(buffer.data(), tangent_size, ambient_size);
-
-    ambient_cursor += ambient_size;
-    tangent_cursor += tangent_size;
-  }
-
-  return true;
-}
-
 bool QuaternionManifold::Plus(const double* x,
                               const double* delta,
                               double* x_plus_delta) const {
diff --git a/internal/ceres/manifold_test.cc b/internal/ceres/manifold_test.cc
index b697d02..60f98d2 100644
--- a/internal/ceres/manifold_test.cc
+++ b/internal/ceres/manifold_test.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2021 Google Inc. All rights reserved.
+// Copyright 2022 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -41,6 +41,7 @@
 #include "ceres/line_manifold.h"
 #include "ceres/manifold_test_utils.h"
 #include "ceres/numeric_diff_options.h"
+#include "ceres/product_manifold.h"
 #include "ceres/rotation.h"
 #include "ceres/sphere_manifold.h"
 #include "ceres/types.h"
@@ -158,54 +159,64 @@
 }
 
 TEST(ProductManifold, Size2) {
-  Manifold* manifold1 = new SubsetManifold(5, {2});
-  Manifold* manifold2 = new SubsetManifold(3, {0, 1});
-  ProductManifold manifold(manifold1, manifold2);
+  SubsetManifold manifold1(5, {2});
+  SubsetManifold manifold2(3, {0, 1});
+  ProductManifold<SubsetManifold, SubsetManifold> manifold(manifold1,
+                                                           manifold2);
 
   EXPECT_EQ(manifold.AmbientSize(),
-            manifold1->AmbientSize() + manifold2->AmbientSize());
+            manifold1.AmbientSize() + manifold2.AmbientSize());
   EXPECT_EQ(manifold.TangentSize(),
-            manifold1->TangentSize() + manifold2->TangentSize());
+            manifold1.TangentSize() + manifold2.TangentSize());
 }
 
 TEST(ProductManifold, Size3) {
-  Manifold* manifold1 = new SubsetManifold(5, {2});
-  Manifold* manifold2 = new SubsetManifold(3, {0, 1});
-  Manifold* manifold3 = new SubsetManifold(4, {1});
+  SubsetManifold manifold1(5, {2});
+  SubsetManifold manifold2(3, {0, 1});
+  SubsetManifold manifold3(4, {1});
 
-  ProductManifold manifold(manifold1, manifold2, manifold3);
+  ProductManifold<SubsetManifold, SubsetManifold, SubsetManifold> manifold(
+      manifold1, manifold2, manifold3);
 
   EXPECT_EQ(manifold.AmbientSize(),
-            manifold1->AmbientSize() + manifold2->AmbientSize() +
-                manifold3->AmbientSize());
+            manifold1.AmbientSize() + manifold2.AmbientSize() +
+                manifold3.AmbientSize());
   EXPECT_EQ(manifold.TangentSize(),
-            manifold1->TangentSize() + manifold2->TangentSize() +
-                manifold3->TangentSize());
+            manifold1.TangentSize() + manifold2.TangentSize() +
+                manifold3.TangentSize());
 }
 
 TEST(ProductManifold, Size4) {
-  Manifold* manifold1 = new SubsetManifold(5, {2});
-  Manifold* manifold2 = new SubsetManifold(3, {0, 1});
-  Manifold* manifold3 = new SubsetManifold(4, {1});
-  Manifold* manifold4 = new SubsetManifold(2, {0});
+  SubsetManifold manifold1(5, {2});
+  SubsetManifold manifold2(3, {0, 1});
+  SubsetManifold manifold3(4, {1});
+  SubsetManifold manifold4(2, {0});
 
-  ProductManifold manifold(manifold1, manifold2, manifold3, manifold4);
+  ProductManifold<SubsetManifold,
+                  SubsetManifold,
+                  SubsetManifold,
+                  SubsetManifold>
+      manifold(manifold1, manifold2, manifold3, manifold4);
 
   EXPECT_EQ(manifold.AmbientSize(),
-            manifold1->AmbientSize() + manifold2->AmbientSize() +
-                manifold3->AmbientSize() + manifold4->AmbientSize());
+            manifold1.AmbientSize() + manifold2.AmbientSize() +
+                manifold3.AmbientSize() + manifold4.AmbientSize());
   EXPECT_EQ(manifold.TangentSize(),
-            manifold1->TangentSize() + manifold2->TangentSize() +
-                manifold3->TangentSize() + manifold4->TangentSize());
+            manifold1.TangentSize() + manifold2.TangentSize() +
+                manifold3.TangentSize() + manifold4.TangentSize());
 }
 
 TEST(ProductManifold, NormalFunctionTest) {
-  Manifold* manifold1 = new SubsetManifold(5, {2});
-  Manifold* manifold2 = new SubsetManifold(3, {0, 1});
-  Manifold* manifold3 = new SubsetManifold(4, {1});
-  Manifold* manifold4 = new SubsetManifold(2, {0});
+  SubsetManifold manifold1(5, {2});
+  SubsetManifold manifold2(3, {0, 1});
+  SubsetManifold manifold3(4, {1});
+  SubsetManifold manifold4(2, {0});
 
-  ProductManifold manifold(manifold1, manifold2, manifold3, manifold4);
+  ProductManifold<SubsetManifold,
+                  SubsetManifold,
+                  SubsetManifold,
+                  SubsetManifold>
+      manifold(manifold1, manifold2, manifold3, manifold4);
 
   for (int trial = 0; trial < kNumTrials; ++trial) {
     const Vector x = Vector::Random(manifold.AmbientSize());
@@ -218,29 +229,29 @@
     int ambient_cursor = 0;
     int tangent_cursor = 0;
 
-    EXPECT_TRUE(manifold1->Plus(&x[ambient_cursor],
-                                &delta[tangent_cursor],
-                                &x_plus_delta_expected[ambient_cursor]));
-    ambient_cursor += manifold1->AmbientSize();
-    tangent_cursor += manifold1->TangentSize();
+    EXPECT_TRUE(manifold1.Plus(&x[ambient_cursor],
+                               &delta[tangent_cursor],
+                               &x_plus_delta_expected[ambient_cursor]));
+    ambient_cursor += manifold1.AmbientSize();
+    tangent_cursor += manifold1.TangentSize();
 
-    EXPECT_TRUE(manifold2->Plus(&x[ambient_cursor],
-                                &delta[tangent_cursor],
-                                &x_plus_delta_expected[ambient_cursor]));
-    ambient_cursor += manifold2->AmbientSize();
-    tangent_cursor += manifold2->TangentSize();
+    EXPECT_TRUE(manifold2.Plus(&x[ambient_cursor],
+                               &delta[tangent_cursor],
+                               &x_plus_delta_expected[ambient_cursor]));
+    ambient_cursor += manifold2.AmbientSize();
+    tangent_cursor += manifold2.TangentSize();
 
-    EXPECT_TRUE(manifold3->Plus(&x[ambient_cursor],
-                                &delta[tangent_cursor],
-                                &x_plus_delta_expected[ambient_cursor]));
-    ambient_cursor += manifold3->AmbientSize();
-    tangent_cursor += manifold3->TangentSize();
+    EXPECT_TRUE(manifold3.Plus(&x[ambient_cursor],
+                               &delta[tangent_cursor],
+                               &x_plus_delta_expected[ambient_cursor]));
+    ambient_cursor += manifold3.AmbientSize();
+    tangent_cursor += manifold3.TangentSize();
 
-    EXPECT_TRUE(manifold4->Plus(&x[ambient_cursor],
-                                &delta[tangent_cursor],
-                                &x_plus_delta_expected[ambient_cursor]));
-    ambient_cursor += manifold4->AmbientSize();
-    tangent_cursor += manifold4->TangentSize();
+    EXPECT_TRUE(manifold4.Plus(&x[ambient_cursor],
+                               &delta[tangent_cursor],
+                               &x_plus_delta_expected[ambient_cursor]));
+    ambient_cursor += manifold4.AmbientSize();
+    tangent_cursor += manifold4.TangentSize();
 
     for (int i = 0; i < x.size(); ++i) {
       EXPECT_EQ(x_plus_delta[i], x_plus_delta_expected[i]);
@@ -252,9 +263,10 @@
 }
 
 TEST(ProductManifold, ZeroTangentSizeAndEuclidean) {
-  Manifold* subset_manifold = new SubsetManifold(1, {0});
-  Manifold* euclidean_manifold = new EuclideanManifold<2>;
-  ProductManifold manifold(subset_manifold, euclidean_manifold);
+  SubsetManifold subset_manifold(1, {0});
+  EuclideanManifold<2> euclidean_manifold;
+  ProductManifold<SubsetManifold, EuclideanManifold<2>> manifold(
+      subset_manifold, euclidean_manifold);
   EXPECT_EQ(manifold.AmbientSize(), 3);
   EXPECT_EQ(manifold.TangentSize(), 2);
 
@@ -276,9 +288,10 @@
 }
 
 TEST(ProductManifold, EuclideanAndZeroTangentSize) {
-  Manifold* subset_manifold = new SubsetManifold(1, {0});
-  Manifold* euclidean_manifold = new EuclideanManifold<2>;
-  ProductManifold manifold(euclidean_manifold, subset_manifold);
+  SubsetManifold subset_manifold(1, {0});
+  EuclideanManifold<2> euclidean_manifold;
+  ProductManifold<EuclideanManifold<2>, SubsetManifold> manifold(
+      euclidean_manifold, subset_manifold);
   EXPECT_EQ(manifold.AmbientSize(), 3);
   EXPECT_EQ(manifold.TangentSize(), 2);
 
@@ -297,6 +310,124 @@
   }
 }
 
+struct CopyableManifold : ceres::Manifold {
+  CopyableManifold() = default;
+  CopyableManifold(const CopyableManifold&) = default;
+  // Do not care about copy-assignment
+  CopyableManifold& operator=(const CopyableManifold&) = delete;
+  // Not moveable
+  CopyableManifold(CopyableManifold&&) = delete;
+  CopyableManifold& operator=(CopyableManifold&&) = delete;
+
+  int AmbientSize() const override { return 3; }
+  int TangentSize() const override { return 2; }
+
+  bool Plus(const double* x,
+            const double* delta,
+            double* x_plus_delta) const override {
+    return true;
+  }
+
+  bool PlusJacobian(const double* x, double* jacobian) const override {
+    return true;
+  }
+
+  bool RightMultiplyByPlusJacobian(const double* x,
+                                   const int num_rows,
+                                   const double* ambient_matrix,
+                                   double* tangent_matrix) const override {
+    return true;
+  }
+
+  bool Minus(const double* y,
+             const double* x,
+             double* y_minus_x) const override {
+    return true;
+  }
+
+  bool MinusJacobian(const double* x, double* jacobian) const override {
+    return true;
+  }
+};
+
+struct MoveableManifold : ceres::Manifold {
+  MoveableManifold() = default;
+  MoveableManifold(MoveableManifold&&) = default;
+  // Do not care about move-assignment
+  MoveableManifold& operator=(MoveableManifold&&) = delete;
+  // Not copyable
+  MoveableManifold(const MoveableManifold&) = delete;
+  MoveableManifold& operator=(const MoveableManifold&) = delete;
+
+  int AmbientSize() const override { return 3; }
+  int TangentSize() const override { return 2; }
+
+  bool Plus(const double* x,
+            const double* delta,
+            double* x_plus_delta) const override {
+    return true;
+  }
+
+  bool PlusJacobian(const double* x, double* jacobian) const override {
+    return true;
+  }
+
+  bool RightMultiplyByPlusJacobian(const double* x,
+                                   const int num_rows,
+                                   const double* ambient_matrix,
+                                   double* tangent_matrix) const override {
+    return true;
+  }
+
+  bool Minus(const double* y,
+             const double* x,
+             double* y_minus_x) const override {
+    return true;
+  }
+
+  bool MinusJacobian(const double* x, double* jacobian) const override {
+    return true;
+  }
+};
+
+TEST(ProductManifold, CopyableOnly) {
+  ProductManifold<CopyableManifold, EuclideanManifold<3>> manifold1{
+      CopyableManifold{}, EuclideanManifold<3>{}};
+
+  CopyableManifold inner2;
+  ProductManifold<CopyableManifold, EuclideanManifold<3>> manifold2{
+      inner2, EuclideanManifold<3>{}};
+
+  EXPECT_EQ(manifold1.AmbientSize(), manifold2.AmbientSize());
+  EXPECT_EQ(manifold1.TangentSize(), manifold2.TangentSize());
+}
+
+TEST(ProductManifold, MoveableOnly) {
+  ProductManifold<MoveableManifold, EuclideanManifold<3>> manifold1{
+      MoveableManifold{}, EuclideanManifold<3>{}};
+
+  MoveableManifold inner2;
+  ProductManifold<MoveableManifold, EuclideanManifold<3>> manifold2{
+      std::move(inner2), EuclideanManifold<3>{}};
+
+  EXPECT_EQ(manifold1.AmbientSize(), manifold2.AmbientSize());
+  EXPECT_EQ(manifold1.TangentSize(), manifold2.TangentSize());
+}
+
+TEST(ProductManifold, CopyableOrMoveable) {
+  const CopyableManifold inner12{};
+  ProductManifold<MoveableManifold, CopyableManifold> manifold1{
+      MoveableManifold{}, inner12};
+
+  MoveableManifold inner21;
+  CopyableManifold inner22;
+  ProductManifold<MoveableManifold, CopyableManifold> manifold2{
+      std::move(inner21), inner22};
+
+  EXPECT_EQ(manifold1.AmbientSize(), manifold2.AmbientSize());
+  EXPECT_EQ(manifold1.TangentSize(), manifold2.TangentSize());
+}
+
 TEST(QuaternionManifold, PlusPiBy2) {
   QuaternionManifold manifold;
   Vector x = Vector::Zero(4);