diff --git a/sophus/se2.hpp b/sophus/se2.hpp
index 60d4654c..f68ca5ef 100644
--- a/sophus/se2.hpp
+++ b/sophus/se2.hpp
@@ -6,7 +6,7 @@
 #include "so2.hpp"
 
 namespace Sophus {
-template <class Scalar_, int Options = 0>
+template <class Scalar_, int Options = 0, bool PreserveInvariant = true>
 class SE2;
 using SE2d = SE2<double>;
 using SE2f = SE2<float>;
@@ -15,27 +15,33 @@ using SE2f = SE2<float>;
 namespace Eigen {
 namespace internal {
 
-template <class Scalar_, int Options>
-struct traits<Sophus::SE2<Scalar_, Options>> {
+template <class Scalar_, int Options_, bool PreserveInvariant_>
+struct traits<Sophus::SE2<Scalar_, Options_, PreserveInvariant_>> {
+  static constexpr int Options = Options_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
   using Scalar = Scalar_;
   using TranslationType = Sophus::Vector2<Scalar, Options>;
-  using SO2Type = Sophus::SO2<Scalar, Options>;
+  using SO2Type = Sophus::SO2<Scalar, Options, PreserveInvariant>;
 };
 
-template <class Scalar_, int Options>
-struct traits<Map<Sophus::SE2<Scalar_>, Options>>
-    : traits<Sophus::SE2<Scalar_, Options>> {
+template <class Scalar_, int SE2Options, bool PreserveInvariant_, int Options_>
+struct traits<
+    Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_>, Options_>>
+    : traits<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_>> {
   using Scalar = Scalar_;
-  using TranslationType = Map<Sophus::Vector2<Scalar>, Options>;
-  using SO2Type = Map<Sophus::SO2<Scalar>, Options>;
+  using TranslationType = Map<Sophus::Vector2<Scalar>, Options_>;
+  using SO2Type =
+      Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_>, Options_>;
 };
 
-template <class Scalar_, int Options>
-struct traits<Map<Sophus::SE2<Scalar_> const, Options>>
-    : traits<Sophus::SE2<Scalar_, Options> const> {
+template <class Scalar_, int SE2Options, bool PreserveInvariant_, int Options_>
+struct traits<
+    Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_> const, Options_>>
+    : traits<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_> const> {
   using Scalar = Scalar_;
-  using TranslationType = Map<Sophus::Vector2<Scalar> const, Options>;
-  using SO2Type = Map<Sophus::SO2<Scalar> const, Options>;
+  using TranslationType = Map<Sophus::Vector2<Scalar> const, Options_>;
+  using SO2Type =
+      Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_> const, Options_>;
 };
 }  // namespace internal
 }  // namespace Eigen
@@ -78,6 +84,8 @@ class SE2Base {
   using Hyperplane = Hyperplane2<Scalar>;
   using Tangent = Vector<Scalar, DoF>;
   using Adjoint = Matrix<Scalar, DoF, DoF>;
+  static constexpr bool PreserveInvariant =
+      Eigen::internal::traits<Derived>::PreserveInvariant;
 
   /// For binary operations the return type is determined with the
   /// ScalarBinaryOpTraits feature of Eigen. This allows mixing concrete and Map
@@ -88,7 +96,9 @@ class SE2Base {
       Scalar, typename OtherDerived::Scalar>::ReturnType;
 
   template <typename OtherDerived>
-  using SE2Product = SE2<ReturnScalar<OtherDerived>>;
+  using SE2Product =
+      SE2<ReturnScalar<OtherDerived>, Eigen::internal::traits<Derived>::Options,
+          PreserveInvariant && OtherDerived::PreserveInvariant>;
 
   template <typename PointDerived>
   using PointProduct = Vector2<ReturnScalar<PointDerived>>;
@@ -115,9 +125,12 @@ class SE2Base {
   /// Returns copy of instance casted to NewScalarType.
   ///
   template <class NewScalarType>
-  SOPHUS_FUNC SE2<NewScalarType> cast() const {
-    return SE2<NewScalarType>(so2().template cast<NewScalarType>(),
-                              translation().template cast<NewScalarType>());
+  SOPHUS_FUNC SE2<NewScalarType, Eigen::internal::traits<Derived>::Options,
+                  PreserveInvariant>
+  cast() const {
+    return SE2<NewScalarType, Eigen::internal::traits<Derived>::Options,
+               PreserveInvariant>(so2().template cast<NewScalarType>(),
+                                  translation().template cast<NewScalarType>());
   }
 
   /// Returns derivative of  this * exp(x)  wrt x at x=0.
@@ -156,9 +169,12 @@ class SE2Base {
 
   /// Returns group inverse.
   ///
-  SOPHUS_FUNC SE2<Scalar> inverse() const {
-    SO2<Scalar> const invR = so2().inverse();
-    return SE2<Scalar>(invR, invR * (translation() * Scalar(-1)));
+  SOPHUS_FUNC
+  SE2<Scalar, Eigen::internal::traits<Derived>::Options, PreserveInvariant>
+  inverse() const {
+    auto const invR = so2().inverse();
+    return SE2<Scalar, Eigen::internal::traits<Derived>::Options,
+               PreserveInvariant>(invR, invR * (translation() * Scalar(-1)));
   }
 
   /// Logarithmic map
@@ -388,12 +404,13 @@ class SE2Base {
 };
 
 /// SE2 using default storage; derived from SE2Base.
-template <class Scalar_, int Options>
-class SE2 : public SE2Base<SE2<Scalar_, Options>> {
+template <class Scalar_, int Options, bool PreserveInvariant_>
+class SE2 : public SE2Base<SE2<Scalar_, Options, PreserveInvariant_>> {
  public:
-  using Base = SE2Base<SE2<Scalar_, Options>>;
+  using Base = SE2Base<SE2<Scalar_, Options, PreserveInvariant_>>;
   static int constexpr DoF = Base::DoF;
   static int constexpr num_parameters = Base::num_parameters;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
 
   using Scalar = Scalar_;
   using Transformation = typename Base::Transformation;
@@ -401,7 +418,7 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   using HomogeneousPoint = typename Base::HomogeneousPoint;
   using Tangent = typename Base::Tangent;
   using Adjoint = typename Base::Adjoint;
-  using SO2Member = SO2<Scalar, Options>;
+  using SO2Member = SO2<Scalar, Options, PreserveInvariant>;
   using TranslationMember = Vector2<Scalar, Options>;
 
   using Base::operator=;
@@ -448,7 +465,7 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   /// of 1.
   ///
   SOPHUS_FUNC
-  SE2(typename SO2<Scalar>::Transformation const& rotation_matrix,
+  SE2(typename SO2Member::Transformation const& rotation_matrix,
       Point const& translation)
       : so2_(rotation_matrix), translation_(translation) {}
 
@@ -606,9 +623,10 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   /// ``expmat(.)`` being the matrix exponential and ``hat(.)`` the hat-operator
   /// of SE(2), see below.
   ///
-  SOPHUS_FUNC static SE2<Scalar> exp(Tangent const& a) {
+  SOPHUS_FUNC static SE2<Scalar, Options, PreserveInvariant> exp(
+      Tangent const& a) {
     Scalar theta = a[2];
-    SO2<Scalar> so2 = SO2<Scalar>::exp(theta);
+    SO2Member so2 = SO2Member::exp(theta);
     Scalar sin_theta_by_theta;
     Scalar one_minus_cos_theta_by_theta;
     using std::abs;
@@ -626,16 +644,18 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
     Vector2<Scalar> trans(
         sin_theta_by_theta * a[0] - one_minus_cos_theta_by_theta * a[1],
         one_minus_cos_theta_by_theta * a[0] + sin_theta_by_theta * a[1]);
-    return SE2<Scalar>(so2, trans);
+    return SE2<Scalar, Options, PreserveInvariant>(so2, trans);
   }
 
   /// Returns closest SE3 given arbitrary 4x4 matrix.
   ///
   template <class S = Scalar>
-  static SOPHUS_FUNC std::enable_if_t<std::is_floating_point<S>::value, SE2>
+  static SOPHUS_FUNC std::enable_if_t<std::is_floating_point<S>::value,
+                                      SE2<Scalar, Options, PreserveInvariant>>
   fitToSE2(Matrix3<Scalar> const& T) {
-    return SE2(SO2<Scalar>::fitToSO2(T.template block<2, 2>(0, 0)),
-               T.template block<2, 1>(0, 2));
+    return SE2<Scalar, Options, PreserveInvariant>(
+        SO2Member::fitToSO2(T.template block<2, 2>(0, 0)),
+        T.template block<2, 1>(0, 2));
   }
 
   /// Returns the ith infinitesimal generators of SE(2).
@@ -709,7 +729,7 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   /// Construct pure rotation.
   ///
   static SOPHUS_FUNC SE2 rot(Scalar const& x) {
-    return SE2(SO2<Scalar>(x), Sophus::Vector2<Scalar>::Zero());
+    return SE2(SO2Member(x), Sophus::Vector2<Scalar>::Zero());
   }
 
   /// Draw uniform sample from SE(2) manifold.
@@ -719,7 +739,7 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   template <class UniformRandomBitGenerator>
   static SE2 sampleUniform(UniformRandomBitGenerator& generator) {
     std::uniform_real_distribution<Scalar> uniform(Scalar(-1), Scalar(1));
-    return SE2(SO2<Scalar>::sampleUniform(generator),
+    return SE2(SO2Member::sampleUniform(generator),
                Vector2<Scalar>(uniform(generator), uniform(generator)));
   }
 
@@ -727,11 +747,11 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   ///
   template <class T0, class T1>
   static SOPHUS_FUNC SE2 trans(T0 const& x, T1 const& y) {
-    return SE2(SO2<Scalar>(), Vector2<Scalar>(x, y));
+    return SE2(SO2Member(), Vector2<Scalar>(x, y));
   }
 
   static SOPHUS_FUNC SE2 trans(Vector2<Scalar> const& xy) {
-    return SE2(SO2<Scalar>(), xy);
+    return SE2(SO2Member(), xy);
   }
 
   /// Construct x-axis translation.
@@ -774,8 +794,8 @@ class SE2 : public SE2Base<SE2<Scalar_, Options>> {
   TranslationMember translation_;
 };
 
-template <class Scalar, int Options>
-SOPHUS_FUNC SE2<Scalar, Options>::SE2()
+template <class Scalar, int Options, bool PreserveInvariant>
+SOPHUS_FUNC SE2<Scalar, Options, PreserveInvariant>::SE2()
     : translation_(TranslationMember::Zero()) {
   static_assert(std::is_standard_layout<SE2>::value,
                 "Assume standard layout for the use of offset of check below.");
@@ -795,11 +815,13 @@ namespace Eigen {
 /// Specialization of Eigen::Map for ``SE2``; derived from SE2Base.
 ///
 /// Allows us to wrap SE2 objects around POD array.
-template <class Scalar_, int Options>
-class Map<Sophus::SE2<Scalar_>, Options>
-    : public Sophus::SE2Base<Map<Sophus::SE2<Scalar_>, Options>> {
+template <class Scalar_, int SE2Options, bool PreserveInvariant_, int Options>
+class Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_>, Options>
+    : public Sophus::SE2Base<
+          Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_>, Options>> {
  public:
-  using Base = Sophus::SE2Base<Map<Sophus::SE2<Scalar_>, Options>>;
+  using Base = Sophus::SE2Base<
+      Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_>, Options>>;
   using Scalar = Scalar_;
   using Transformation = typename Base::Transformation;
   using Point = typename Base::Point;
@@ -814,15 +836,22 @@ class Map<Sophus::SE2<Scalar_>, Options>
   SOPHUS_FUNC
   explicit Map(Scalar* coeffs)
       : so2_(coeffs),
-        translation_(coeffs + Sophus::SO2<Scalar>::num_parameters) {}
+        translation_(coeffs + Sophus::SO2<Scalar, SE2Options,
+                                          PreserveInvariant_>::num_parameters) {
+  }
 
   /// Mutator of SO3
   ///
-  SOPHUS_FUNC Map<Sophus::SO2<Scalar>, Options>& so2() { return so2_; }
+  SOPHUS_FUNC
+  Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_>, Options>& so2() {
+    return so2_;
+  }
 
   /// Accessor of SO3
   ///
-  SOPHUS_FUNC Map<Sophus::SO2<Scalar>, Options> const& so2() const {
+  SOPHUS_FUNC
+  Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_>, Options> const& so2()
+      const {
     return so2_;
   }
 
@@ -839,18 +868,21 @@ class Map<Sophus::SE2<Scalar_>, Options>
   }
 
  protected:
-  Map<Sophus::SO2<Scalar>, Options> so2_;
+  Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_>, Options> so2_;
   Map<Sophus::Vector2<Scalar>, Options> translation_;
 };
 
 /// Specialization of Eigen::Map for ``SE2 const``; derived from SE2Base.
 ///
 /// Allows us to wrap SE2 objects around POD array.
-template <class Scalar_, int Options>
-class Map<Sophus::SE2<Scalar_> const, Options>
-    : public Sophus::SE2Base<Map<Sophus::SE2<Scalar_> const, Options>> {
+template <class Scalar_, int SE2Options, bool PreserveInvariant_, int Options>
+class Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_> const, Options>
+    : public Sophus::SE2Base<
+          Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_> const,
+              Options>> {
  public:
-  using Base = Sophus::SE2Base<Map<Sophus::SE2<Scalar_> const, Options>>;
+  using Base = Sophus::SE2Base<
+      Map<Sophus::SE2<Scalar_, SE2Options, PreserveInvariant_> const, Options>>;
   using Scalar = Scalar_;
   using Transformation = typename Base::Transformation;
   using Point = typename Base::Point;
@@ -863,23 +895,28 @@ class Map<Sophus::SE2<Scalar_> const, Options>
 
   SOPHUS_FUNC explicit Map(Scalar const* coeffs)
       : so2_(coeffs),
-        translation_(coeffs + Sophus::SO2<Scalar>::num_parameters) {}
+        translation_(coeffs + Sophus::SO2<Scalar, SE2Options,
+                                          PreserveInvariant_>::num_parameters) {
+  }
 
   /// Accessor of SO3
   ///
-  SOPHUS_FUNC Map<Sophus::SO2<Scalar> const, Options> const& so2() const {
+  SOPHUS_FUNC
+  Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_> const, Options> const&
+  so2() const {
     return so2_;
   }
 
   /// Accessor of translation vector
   ///
-  SOPHUS_FUNC Map<Sophus::Vector2<Scalar> const, Options> const& translation()
-      const {
+  SOPHUS_FUNC
+  Map<Sophus::Vector2<Scalar> const, Options> const& translation() const {
     return translation_;
   }
 
  protected:
-  Map<Sophus::SO2<Scalar> const, Options> const so2_;
+  Map<Sophus::SO2<Scalar, SE2Options, PreserveInvariant_> const, Options> const
+      so2_;
   Map<Sophus::Vector2<Scalar> const, Options> const translation_;
 };
 }  // namespace Eigen
diff --git a/sophus/so2.hpp b/sophus/so2.hpp
index 8c05916a..8bfd5a80 100644
--- a/sophus/so2.hpp
+++ b/sophus/so2.hpp
@@ -13,7 +13,7 @@
 #include "types.hpp"
 
 namespace Sophus {
-template <class Scalar_, int Options = 0>
+template <class Scalar_, int Options = 0, bool PreserveInvariant = true>
 class SO2;
 using SO2d = SO2<double>;
 using SO2f = SO2<float>;
@@ -22,25 +22,30 @@ using SO2f = SO2<float>;
 namespace Eigen {
 namespace internal {
 
-template <class Scalar_, int Options_>
-struct traits<Sophus::SO2<Scalar_, Options_>> {
+template <class Scalar_, int Options_, bool PreserveInvariant_>
+struct traits<Sophus::SO2<Scalar_, Options_, PreserveInvariant_>> {
   static constexpr int Options = Options_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
   using Scalar = Scalar_;
   using ComplexType = Sophus::Vector2<Scalar, Options>;
 };
 
-template <class Scalar_, int Options_>
-struct traits<Map<Sophus::SO2<Scalar_>, Options_>>
-    : traits<Sophus::SO2<Scalar_, Options_>> {
+template <class Scalar_, int SO2Options, bool PreserveInvariant_, int Options_>
+struct traits<
+    Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>, Options_>>
+    : traits<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>> {
   static constexpr int Options = Options_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
   using Scalar = Scalar_;
   using ComplexType = Map<Sophus::Vector2<Scalar>, Options>;
 };
 
-template <class Scalar_, int Options_>
-struct traits<Map<Sophus::SO2<Scalar_> const, Options_>>
-    : traits<Sophus::SO2<Scalar_, Options_> const> {
+template <class Scalar_, int SO2Options, bool PreserveInvariant_, int Options_>
+struct traits<
+    Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_> const, Options_>>
+    : traits<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_> const> {
   static constexpr int Options = Options_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
   using Scalar = Scalar_;
   using ComplexType = Map<Sophus::Vector2<Scalar> const, Options>;
 };
@@ -98,6 +103,10 @@ class SO2Base {
   using Hyperplane = Hyperplane2<Scalar>;
   using Tangent = Scalar;
   using Adjoint = Scalar;
+  static constexpr bool PreserveInvariant =
+      Eigen::internal::traits<Derived>::PreserveInvariant;
+  static constexpr auto Normalize = std::bool_constant<true>{};
+  static constexpr auto DontNormalize = std::bool_constant<false>{};
 
   /// For binary operations the return type is determined with the
   /// ScalarBinaryOpTraits feature of Eigen. This allows mixing concrete and Map
@@ -108,7 +117,8 @@ class SO2Base {
       Scalar, typename OtherDerived::Scalar>::ReturnType;
 
   template <typename OtherDerived>
-  using SO2Product = SO2<ReturnScalar<OtherDerived>>;
+  using SO2Product = SO2<ReturnScalar<OtherDerived>, Options,
+                         PreserveInvariant && OtherDerived::PreserveInvariant>;
 
   template <typename PointDerived>
   using PointProduct = Vector2<ReturnScalar<PointDerived>>;
@@ -129,8 +139,16 @@ class SO2Base {
   /// Returns copy of instance casted to NewScalarType.
   ///
   template <class NewScalarType>
-  SOPHUS_FUNC SO2<NewScalarType> cast() const {
-    return SO2<NewScalarType>(unit_complex().template cast<NewScalarType>());
+  SOPHUS_FUNC SO2<NewScalarType, Options, PreserveInvariant> cast() const {
+    return SO2<NewScalarType, Options, PreserveInvariant>(
+        unit_complex().template cast<NewScalarType>());
+  }
+
+  /// Cast to SO2 type with different invariant-preserving property
+  ///
+  template <bool OtherPreserveInvariant>
+  SOPHUS_FUNC operator SO2<Scalar, Options, OtherPreserveInvariant>() const {
+    return SO2<Scalar, Options, OtherPreserveInvariant>(unit_complex());
   }
 
   /// This provides unsafe read/write access to internal data. SO(2) is
@@ -146,8 +164,10 @@ class SO2Base {
 
   /// Returns group inverse.
   ///
-  SOPHUS_FUNC SO2<Scalar> inverse() const {
-    return SO2<Scalar>(unit_complex().x(), -unit_complex().y());
+  SOPHUS_FUNC SO2<Scalar, Options, PreserveInvariant> inverse() const {
+    // Conjugation preserves the norm, thus no re-normalization is necessary
+    return SO2<Scalar, Options, PreserveInvariant>(
+        unit_complex().x(), -unit_complex().y(), DontNormalize);
   }
 
   /// Logarithmic map
@@ -209,6 +229,7 @@ class SO2Base {
   template <typename OtherDerived>
   SOPHUS_FUNC SO2Product<OtherDerived> operator*(
       SO2Base<OtherDerived> const& other) const {
+    using Product = SO2Product<OtherDerived>;
     using ResultT = ReturnScalar<OtherDerived>;
     Scalar const lhs_real = unit_complex().x();
     Scalar const lhs_imag = unit_complex().y();
@@ -218,20 +239,30 @@ class SO2Base {
     ResultT const result_real = lhs_real * rhs_real - lhs_imag * rhs_imag;
     ResultT const result_imag = lhs_real * rhs_imag + lhs_imag * rhs_real;
 
-    ResultT const squared_norm =
-        result_real * result_real + result_imag * result_imag;
-    // We can assume that the squared-norm is close to 1 since we deal with a
-    // unit complex number. Due to numerical precision issues, there might
-    // be a small drift after pose concatenation. Hence, we need to renormalizes
-    // the complex number here.
-    // Since squared-norm is close to 1, we do not need to calculate the costly
-    // square-root, but can use an approximation around 1 (see
-    // http://stackoverflow.com/a/12934750 for details).
-    if (squared_norm != ResultT(1.0)) {
-      ResultT const scale = ResultT(2.0) / (ResultT(1.0) + squared_norm);
-      return SO2Product<OtherDerived>(result_real * scale, result_imag * scale);
+    if constexpr (Product::PreserveInvariant) {
+      ResultT const squared_norm =
+          result_real * result_real + result_imag * result_imag;
+      // We can assume that the squared-norm is close to 1 since we deal with
+      // a unit complex number. Due to numerical precision issues, there might
+      // be a small drift after pose concatenation. Hence, we need to
+      // renormalize the complex number here. Since squared-norm is close to
+      // 1, we do not need to calculate the costly square-root, but can use an
+      // approximation around 1 (see http://stackoverflow.com/a/12934750 for
+      // details).
+      if (squared_norm != ResultT(1.0)) {
+        ResultT const scale = ResultT(2.0) / (ResultT(1.0) + squared_norm);
+        ResultT const scaled_real = result_real * scale;
+        ResultT const scaled_imag = result_imag * scale;
+        ResultT const scaled_squared_norm =
+            scaled_real * scaled_real + scaled_imag * scaled_imag;
+        if (scaled_squared_norm != ResultT(1.0)) {
+          return Product(scaled_real, scaled_imag, Normalize);
+        }
+        return Product(scaled_real, scaled_imag, DontNormalize);
+      }
+      return Product(result_real, result_imag, DontNormalize);
     }
-    return SO2Product<OtherDerived>(result_real, result_imag);
+    return Product(result_real, result_imag, DontNormalize);
   }
 
   /// Group action on 2-points.
@@ -334,7 +365,9 @@ class SO2Base {
   ///
   SOPHUS_FUNC void setComplex(Point const& complex) {
     unit_complex_nonconst() = complex;
-    normalize();
+    if constexpr (PreserveInvariant) {
+      normalize();
+    }
   }
 
   /// Accessor of unit quaternion.
@@ -355,12 +388,13 @@ class SO2Base {
 };
 
 /// SO2 using  default storage; derived from SO2Base.
-template <class Scalar_, int Options>
-class SO2 : public SO2Base<SO2<Scalar_, Options>> {
+template <class Scalar_, int Options, bool PreserveInvariant_>
+class SO2 : public SO2Base<SO2<Scalar_, Options, PreserveInvariant_>> {
  public:
-  using Base = SO2Base<SO2<Scalar_, Options>>;
+  using Base = SO2Base<SO2<Scalar_, Options, PreserveInvariant_>>;
   static int constexpr DoF = Base::DoF;
   static int constexpr num_parameters = Base::num_parameters;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
 
   using Scalar = Scalar_;
   using Transformation = typename Base::Transformation;
@@ -371,7 +405,7 @@ class SO2 : public SO2Base<SO2<Scalar_, Options>> {
   using ComplexMember = Vector2<Scalar, Options>;
 
   /// ``Base`` is friend so unit_complex_nonconst can be accessed from ``Base``.
-  friend class SO2Base<SO2<Scalar, Options>>;
+  friend class SO2Base<SO2<Scalar, Options, PreserveInvariant>>;
 
   using Base::operator=;
 
@@ -394,7 +428,11 @@ class SO2 : public SO2Base<SO2<Scalar_, Options>> {
   ///
   template <class OtherDerived>
   SOPHUS_FUNC SO2(SO2Base<OtherDerived> const& other)
-      : unit_complex_(other.unit_complex()) {}
+      : unit_complex_(other.unit_complex()) {
+    if constexpr (PreserveInvariant && !OtherDerived::PreserveInvariant) {
+      Base::normalize();
+    }
+  }
 
   /// Constructor from rotation matrix
   ///
@@ -410,15 +448,39 @@ class SO2 : public SO2Base<SO2<Scalar_, Options>> {
 
   /// Constructor from pair of real and imaginary number.
   ///
+  /// Normalization is performed based on PreserveInvariant
+  ///
   /// Precondition: The pair must not be close to zero.
   ///
   SOPHUS_FUNC SO2(Scalar const& real, Scalar const& imag)
       : unit_complex_(real, imag) {
-    Base::normalize();
+    if (PreserveInvariant) {
+      Base::normalize();
+    }
+  }
+
+  /// Constructor from pair of real and imaginary number with optional
+  /// normalization
+  ///
+  /// Normalization is performed based on normalize flag
+  /// SO2Base::Normalize and SO2Base::DontNormalize constants can be used as a
+  /// shorthand for calling this constructor
+  ///
+  /// Precondition: The pair must not be close to zero.
+  ///
+  template <bool normalize>
+  SO2(Scalar const& real, Scalar const& imag,
+      const std::bool_constant<normalize>)
+      : unit_complex_(real, imag) {
+    if (normalize) {
+      Base::normalize();
+    }
   }
 
   /// Constructor from 2-vector.
   ///
+  /// Normalization is performed based on PreserveInvariant
+  ///
   /// Precondition: The vector must not be close to zero.
   ///
   template <class D>
@@ -426,13 +488,16 @@ class SO2 : public SO2Base<SO2<Scalar_, Options>> {
       : unit_complex_(complex) {
     static_assert(std::is_same<typename D::Scalar, Scalar>::value,
                   "must be same Scalar type");
-    Base::normalize();
+    if constexpr (PreserveInvariant) {
+      Base::normalize();
+    }
   }
 
   /// Constructor from an rotation angle.
   ///
   SOPHUS_FUNC explicit SO2(Scalar theta) {
-    unit_complex_nonconst() = SO2<Scalar>::exp(theta).unit_complex();
+    unit_complex_nonconst() =
+        SO2<Scalar, Options, PreserveInvariant>::exp(theta).unit_complex();
   }
 
   /// Accessor of unit complex number
@@ -450,10 +515,12 @@ class SO2 : public SO2Base<SO2<Scalar_, Options>> {
   /// with ``expmat(.)`` being the matrix exponential and ``hat(.)`` being the
   /// hat()-operator of SO(2).
   ///
-  SOPHUS_FUNC static SO2<Scalar> exp(Tangent const& theta) {
+  SOPHUS_FUNC static SO2<Scalar, Options, PreserveInvariant> exp(
+      Tangent const& theta) {
     using std::cos;
     using std::sin;
-    return SO2<Scalar>(cos(theta), sin(theta));
+    return SO2<Scalar, Options, PreserveInvariant>(cos(theta), sin(theta),
+                                                   Base::DontNormalize);
   }
 
   /// Returns derivative of exp(x) wrt. x.
@@ -578,12 +645,15 @@ namespace Eigen {
 ///
 /// Allows us to wrap SO2 objects around POD array (e.g. external c style
 /// complex number / tuple).
-template <class Scalar_, int Options>
-class Map<Sophus::SO2<Scalar_>, Options>
-    : public Sophus::SO2Base<Map<Sophus::SO2<Scalar_>, Options>> {
+template <class Scalar_, int SO2Options, bool PreserveInvariant_, int Options>
+class Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>, Options>
+    : public Sophus::SO2Base<
+          Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>, Options>> {
  public:
-  using Base = Sophus::SO2Base<Map<Sophus::SO2<Scalar_>, Options>>;
+  using Base = Sophus::SO2Base<
+      Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>, Options>>;
   using Scalar = Scalar_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
 
   using Transformation = typename Base::Transformation;
   using Point = typename Base::Point;
@@ -592,7 +662,8 @@ class Map<Sophus::SO2<Scalar_>, Options>
   using Adjoint = typename Base::Adjoint;
 
   /// ``Base`` is friend so unit_complex_nonconst can be accessed from ``Base``.
-  friend class Sophus::SO2Base<Map<Sophus::SO2<Scalar_>, Options>>;
+  friend class Sophus::SO2Base<
+      Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_>, Options>>;
 
   using Base::operator=;
   using Base::operator*=;
@@ -623,12 +694,16 @@ class Map<Sophus::SO2<Scalar_>, Options>
 ///
 /// Allows us to wrap SO2 objects around POD array (e.g. external c style
 /// complex number / tuple).
-template <class Scalar_, int Options>
-class Map<Sophus::SO2<Scalar_> const, Options>
-    : public Sophus::SO2Base<Map<Sophus::SO2<Scalar_> const, Options>> {
+template <class Scalar_, int SO2Options, bool PreserveInvariant_, int Options>
+class Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_> const, Options>
+    : public Sophus::SO2Base<
+          Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_> const,
+              Options>> {
  public:
-  using Base = Sophus::SO2Base<Map<Sophus::SO2<Scalar_> const, Options>>;
+  using Base = Sophus::SO2Base<
+      Map<Sophus::SO2<Scalar_, SO2Options, PreserveInvariant_> const, Options>>;
   using Scalar = Scalar_;
+  static constexpr bool PreserveInvariant = PreserveInvariant_;
   using Transformation = typename Base::Transformation;
   using Point = typename Base::Point;
   using HomogeneousPoint = typename Base::HomogeneousPoint;
@@ -642,8 +717,8 @@ class Map<Sophus::SO2<Scalar_> const, Options>
 
   /// Accessor of unit complex number.
   ///
-  SOPHUS_FUNC Map<Sophus::Vector2<Scalar> const, Options> const& unit_complex()
-      const {
+  SOPHUS_FUNC
+  Map<Sophus::Vector2<Scalar> const, Options> const& unit_complex() const {
     return unit_complex_;
   }
 
diff --git a/test/core/test_se2.cpp b/test/core/test_se2.cpp
index c0e5d8f6..8d6705e7 100644
--- a/test/core/test_se2.cpp
+++ b/test/core/test_se2.cpp
@@ -9,12 +9,16 @@
 namespace Eigen {
 template class Map<Sophus::SE2<double>>;
 template class Map<Sophus::SE2<double> const>;
+template class Map<Sophus::SE2<double, 0, false>>;
+template class Map<Sophus::SE2<double, 0, false> const>;
 }  // namespace Eigen
 
 namespace Sophus {
 
 template class SE2<double, Eigen::AutoAlign>;
 template class SE2<double, Eigen::DontAlign>;
+template class SE2<double, Eigen::AutoAlign, false>;
+template class SE2<double, Eigen::DontAlign, false>;
 #if SOPHUS_CERES
 template class SE2<ceres::Jet<double, 3>>;
 #endif
@@ -72,6 +76,7 @@ class Tests {
     passed &= testRawDataAcces();
     passed &= testMutatingAccessors();
     passed &= testConstructors();
+    passed &= testNormalizationModes();
     passed &= testFit();
     processTestResult(passed);
   }
@@ -216,6 +221,72 @@ class Tests {
     return passed;
   }
 
+  bool testNormalizationModes() {
+    bool passed = true;
+    // Proactively preserving SO2 norm invariant (default)
+    using SE2N = SE2<Scalar>;
+    // Non-enforcing norm invariant
+    using SE2U = SE2<Scalar, 0, false>;
+    using SO2U = SO2<Scalar, 0, false>;
+
+    SE2U unnormalized(SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize),
+                      Point(Scalar(3), Scalar(4)));
+    SE2N normalized_set_complex;
+    SE2U unnormalized_set_complex;
+    normalized_set_complex.so2().setComplex(Point(Scalar(2), Scalar(0)));
+    unnormalized_set_complex.so2().setComplex(Point(Scalar(2), Scalar(0)));
+    SE2N normalized = unnormalized;
+
+    SOPHUS_TEST_APPROX(passed, normalized.so2().unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_APPROX(passed,
+                       normalized_set_complex.so2().unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_EQUAL(passed, unnormalized_set_complex.so2().unit_complex().x(),
+                      Scalar(2),
+                      "Unnormalized SE2 setComplex preserves magnitude");
+
+    unnormalized.normalize();
+    SOPHUS_TEST_APPROX(passed, unnormalized.so2().unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+
+    SE2U inverse = SE2U(SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize),
+                        Point(Scalar(1), Scalar(0)))
+                       .inverse();
+    SE2N normalized_inverse =
+        SE2N(SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize),
+             Point(Scalar(1), Scalar(0)))
+            .inverse();
+    SOPHUS_TEST_APPROX(passed, normalized_inverse.so2().unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_EQUAL(passed, inverse.so2().unit_complex().x(), Scalar(2),
+                      "Unnormalized inverse preserves rotation magnitude");
+
+    static_assert(
+        std::is_same_v<typename SE2N::template SE2Product<SE2N>, SE2N>);
+    static_assert(
+        std::is_same_v<typename SE2U::template SE2Product<SE2N>, SE2U>);
+    static_assert(
+        std::is_same_v<typename SE2N::template SE2Product<SE2U>, SE2U>);
+    static_assert(
+        std::is_same_v<typename SE2U::template SE2Product<SE2U>, SE2U>);
+
+    SE2N normalized_product =
+        SE2N::trans(Scalar(0), Scalar(0)) *
+        SE2N(SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize),
+             Point(Scalar(0), Scalar(0)));
+    SE2U unnormalized_product =
+        SE2N::trans(Scalar(0), Scalar(0)) *
+        SE2U(SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize),
+             Point(Scalar(0), Scalar(0)));
+    SOPHUS_TEST_APPROX(passed, normalized_product.so2().unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_APPROX(passed, unnormalized_product.so2().unit_complex().norm(),
+                       Scalar(2), Constants<Scalar>::epsilon(), "");
+
+    return passed;
+  }
+
   template <class S = Scalar>
   std::enable_if_t<std::is_floating_point<S>::value, bool> testFit() {
     bool passed = true;
diff --git a/test/core/test_so2.cpp b/test/core/test_so2.cpp
index 9039c3a6..03c92393 100644
--- a/test/core/test_so2.cpp
+++ b/test/core/test_so2.cpp
@@ -8,12 +8,16 @@
 namespace Eigen {
 template class Map<Sophus::SO2<double>>;
 template class Map<Sophus::SO2<double> const>;
+template class Map<Sophus::SO2<double, 0, false>>;
+template class Map<Sophus::SO2<double, 0, false> const>;
 }  // namespace Eigen
 
 namespace Sophus {
 
 template class SO2<double, Eigen::AutoAlign>;
 template class SO2<float, Eigen::DontAlign>;
+template class SO2<double, Eigen::AutoAlign, false>;
+template class SO2<float, Eigen::DontAlign, false>;
 #if SOPHUS_CERES
 template class SO2<ceres::Jet<double, 3>>;
 #endif
@@ -53,6 +57,7 @@ class Tests {
     passed &= testUnity();
     passed &= testRawDataAcces();
     passed &= testConstructors();
+    passed &= testNormalizationModes();
     passed &= testFit();
     processTestResult(passed);
   }
@@ -144,6 +149,70 @@ class Tests {
     return passed;
   }
 
+  bool testNormalizationModes() {
+    bool passed = true;
+    // Proactively preserving norm invariant (default)
+    using SO2N = SO2<Scalar>;
+    // Non-enforcing norm invariant
+    using SO2U = SO2<Scalar, 0, false>;
+
+    Vector2<Scalar> diverged(Scalar(2), Scalar(0));
+    SO2N normalized(diverged);
+    SO2U unnormalized(diverged);
+    SO2N normalized_set_complex;
+    SO2U unnormalized_set_complex;
+    normalized_set_complex.setComplex(diverged);
+    unnormalized_set_complex.setComplex(diverged);
+
+    SOPHUS_TEST_APPROX(passed, normalized.unit_complex().norm(), Scalar(1),
+                       Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_EQUAL(passed, unnormalized.unit_complex().x(), diverged.x(),
+                      "Unnormalized constructor preserves real part");
+    SOPHUS_TEST_EQUAL(passed, unnormalized.unit_complex().y(), diverged.y(),
+                      "Unnormalized constructor preserves imaginary part");
+    SOPHUS_TEST_APPROX(passed, normalized_set_complex.unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_EQUAL(passed, unnormalized_set_complex.unit_complex().x(),
+                      diverged.x(),
+                      "Unnormalized setComplex preserves real part");
+    SOPHUS_TEST_EQUAL(passed, unnormalized_set_complex.unit_complex().y(),
+                      diverged.y(),
+                      "Unnormalized setComplex preserves imaginary part");
+
+    SO2N normalized_from_unnormalized = unnormalized;
+    SOPHUS_TEST_APPROX(passed,
+                       normalized_from_unnormalized.unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(),
+                       "Copy to normalized performs normalization");
+
+    unnormalized.normalize();
+    SOPHUS_TEST_APPROX(passed, unnormalized.unit_complex().norm(), Scalar(1),
+                       Constants<Scalar>::epsilon(), "");
+
+    static_assert(
+        std::is_same_v<typename SO2N::template SO2Product<SO2N>, SO2N>);
+    static_assert(
+        std::is_same_v<typename SO2U::template SO2Product<SO2N>, SO2U>);
+    static_assert(
+        std::is_same_v<typename SO2N::template SO2Product<SO2U>, SO2U>);
+    static_assert(
+        std::is_same_v<typename SO2U::template SO2Product<SO2U>, SO2U>);
+
+    SO2N normalized_product =
+        SO2N(Scalar(2), Scalar(0)) * normalized_from_unnormalized;
+    SO2U unnormalized_product =
+        SO2U(Scalar(2), Scalar(0), SO2U::DontNormalize) *
+        SO2U(Scalar(1), Scalar(0), SO2U::DontNormalize);
+    SOPHUS_TEST_APPROX(passed, normalized_product.unit_complex().norm(),
+                       Scalar(1), Constants<Scalar>::epsilon(), "");
+    SOPHUS_TEST_EQUAL(passed, unnormalized_product.unit_complex().x(),
+                      Scalar(2), "");
+    SOPHUS_TEST_EQUAL(passed, unnormalized_product.unit_complex().y(),
+                      Scalar(0), "");
+
+    return passed;
+  }
+
   template <class S = Scalar>
   std::enable_if_t<std::is_floating_point<S>::value, bool> testFit() {
     bool passed = true;