diff --git a/include/rotgen/impl/matrix_impl64.hpp b/include/rotgen/impl/matrix_impl64.hpp
index e3345a7..ff9d349 100644
--- a/include/rotgen/impl/matrix_impl64.hpp
+++ b/include/rotgen/impl/matrix_impl64.hpp
@@ -43,6 +43,8 @@ namespace rotgen
     double const& operator()(std::size_t index) const;
 
     matrix_impl64& operator+=(matrix_impl64 const& rhs);
+    matrix_impl64& operator-=(matrix_impl64 const& rhs);
+    matrix_impl64& operator-();
     matrix_impl64& operator*=(matrix_impl64 const& rhs);
     matrix_impl64& operator*=(double  d);
 
diff --git a/include/rotgen/matrix.hpp b/include/rotgen/matrix.hpp
index a481c23..b7ef164 100644
--- a/include/rotgen/matrix.hpp
+++ b/include/rotgen/matrix.hpp
@@ -69,6 +69,17 @@ namespace rotgen
       return *this;
     }
 
+    matrix& operator-=(matrix const& rhs)
+    {
+      static_cast<parent&>(*this) -= static_cast<parent const&>(rhs);
+      return *this;
+    }
+
+    matrix& operator-()
+    {
+      return static_cast<parent const&>(*this).operator-();
+    }
+
     matrix& operator*=(matrix const& rhs)
     {
       static_cast<parent&>(*this) *= static_cast<parent const&>(rhs);
@@ -89,6 +100,13 @@ namespace rotgen
     return that += rhs;
   }
 
+  template<typename S, int R, int C, int O, int MR, int MC>
+  matrix<S,R,C,O,MR,MC> operator-(matrix<S,R,C,O,MR,MC> const& lhs, matrix<S,R,C,O,MR,MC> const& rhs)
+  {
+    matrix<S,R,C,O,MR,MC> that(lhs);
+    return that -= rhs;
+  }
+
   template<typename S, int R, int C, int O, int MR, int MC>
   matrix<S,R,C,O,MR,MC> operator*(matrix<S,R,C,O,MR,MC> const& lhs, matrix<S,R,C,O,MR,MC> const& rhs)
   {
diff --git a/src/matrix_impl64.cpp b/src/matrix_impl64.cpp
index e1fec52..515659a 100644
--- a/src/matrix_impl64.cpp
+++ b/src/matrix_impl64.cpp
@@ -95,6 +95,19 @@ namespace rotgen
     return *this;
   }
 
+  matrix_impl64& matrix_impl64::operator-=(matrix_impl64 const& rhs)
+  {
+    storage_->data -= rhs.storage_->data;
+    return *this;
+  }
+
+  matrix_impl64& matrix_impl64::operator-()
+  {
+    matrix_impl64 result(*this);
+    result.storage_->data = -result.storage_->data;
+    return result;
+  }
+
   matrix_impl64& matrix_impl64::operator*=(matrix_impl64 const& rhs)
   {
     storage_->data *= rhs.storage_->data;