diff --git a/matrix.cpp b/matrix.cpp
index cb8166fb..7a992373 100644
--- a/matrix.cpp
+++ b/matrix.cpp
@@ -392,11 +392,6 @@ matrix<S,T>::rotation (T angle, util::vector<3,T> about)
     } };
 }
 
-//-----------------------------------------------------------------------------
-template <size_t S, typename T>
-const matrix<S,T>
-matrix<S,T>::ZEROES { 0 };
-
 
 //-----------------------------------------------------------------------------
 template struct util::matrix<2,float>;
diff --git a/matrix.hpp b/matrix.hpp
index 731ec9a2..2e9d1046 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -78,8 +78,8 @@ namespace util {
         static matrix<4,T> rotation    (T angle, util::vector<3,T> about);
 
         // Constant matrices
-        static const matrix IDENTITY;
-        static const matrix ZEROES;
+        static constexpr matrix identity ();
+        static constexpr matrix zeroes ();
     };
 
 
diff --git a/matrix.ipp b/matrix.ipp
index f23061b8..81ef768d 100644
--- a/matrix.ipp
+++ b/matrix.ipp
@@ -183,6 +183,29 @@ MATRIX_SCALAR_OP(-)
 #undef MATRIX_SCALAR_OP
 
 
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S, typename T>
+constexpr
+util::matrix<S,T>
+util::matrix<S,T>::zeroes (void)
+{
+    return {0};
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <size_t S, typename T>
+constexpr
+util::matrix<S,T>
+util::matrix<S,T>::identity (void)
+{
+    auto m = zeroes ();
+    for (size_t i = 0; i < S; ++i)
+        m[i][i] = 1;
+    return m;
+}
+
+
 ///////////////////////////////////////////////////////////////////////////////
 template <size_t S, typename T>
 constexpr
diff --git a/matrix2.cpp b/matrix2.cpp
index 74386a1f..4bb7f2f4 100644
--- a/matrix2.cpp
+++ b/matrix2.cpp
@@ -33,7 +33,7 @@ template float  util::determinant (const matrix<2,float>&);
 template double util::determinant (const matrix<2,double>&);
 
 
-//-----------------------------------------------------------------------------
+///////////////////////////////////////////////////////////////////////////////
 template <size_t S, typename T>
 matrix<S,T>
 util::inverse (const matrix<S,T> &m)
@@ -46,19 +46,12 @@ util::inverse (const matrix<S,T> &m)
     } / determinant (m);
 }
 
+
+//-----------------------------------------------------------------------------
 template util::matrix<2,float>  util::inverse (const matrix<2,float>&);
 template util::matrix<2,double> util::inverse (const matrix<2,double>&);
 
 
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-const matrix<S,T>
-matrix<S,T>::IDENTITY = { {
-    { 1, 0, },
-    { 0, 1  }
-} };
-
-
 ///////////////////////////////////////////////////////////////////////////////
 template struct util::matrix<2,float>;
 template struct util::matrix<2,double>;
diff --git a/matrix3.cpp b/matrix3.cpp
index 3b8ee1e0..4e9d232c 100644
--- a/matrix3.cpp
+++ b/matrix3.cpp
@@ -61,15 +61,6 @@ template util::matrix<3,float>  util::inverse (const matrix<3,float>&);
 template util::matrix<3,double> util::inverse (const matrix<3,double>&);
 
 
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-const matrix<S,T>
-matrix<S,T>::IDENTITY = { {
-    { 1, 0, 0, },
-    { 0, 1, 0, },
-    { 0, 0, 1  }
-} };
-
 ///////////////////////////////////////////////////////////////////////////////
 template struct util::matrix<3,float>;
 template struct util::matrix<3,double>;
diff --git a/matrix4.cpp b/matrix4.cpp
index d405b369..efa1b9d5 100644
--- a/matrix4.cpp
+++ b/matrix4.cpp
@@ -75,15 +75,6 @@ template util::matrix<4,float> util::inverse  (const matrix<4,float>&);
 template util::matrix<4,double> util::inverse (const matrix<4,double>&);
 
 
-///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-const matrix<S,T>
-matrix<S,T>::IDENTITY = { { { 1, 0, 0, 0 },
-                            { 0, 1, 0, 0 },
-                            { 0, 0, 1, 0 },
-                            { 0, 0, 0, 1 } } };
-
-
 ///////////////////////////////////////////////////////////////////////////////
 template struct util::matrix<4,float>;
 template struct util::matrix<4,double>;
diff --git a/test/matrix.cpp b/test/matrix.cpp
index f9ddd170..844a8ef8 100644
--- a/test/matrix.cpp
+++ b/test/matrix.cpp
@@ -48,7 +48,7 @@ main (void)
     {
         // Identity matrix-vector multiplication
         auto v = util::vector4f { 1.f, 2.f, 3.f, 4.f };
-        auto r = util::matrix4f::IDENTITY * v;
+        auto r = util::matrix4f::identity () * v;
         tap.expect_eq (r, v, "identity matrix-vector multiplication");
     }
 
@@ -101,7 +101,7 @@ main (void)
         bool success = true;
 
         // Ensure identity inverts to identity
-        auto m = util::matrix4f::IDENTITY.inverse ();
+        auto m = util::matrix4f::identity ().inverse ();
         for (size_t r = 0; r < m.rows; ++r)
             for (size_t c = 0; c < m.cols; ++c)
                 if (r == c)
diff --git a/test/quaternion.cpp b/test/quaternion.cpp
index eba9c9a0..63b1948b 100644
--- a/test/quaternion.cpp
+++ b/test/quaternion.cpp
@@ -65,7 +65,7 @@ main (void)
 
     tap.expect_eq (
         quaternionf::identity ().as_matrix (),
-        util::matrix4f::IDENTITY,
+        util::matrix4f::identity (),
         "identity quaternion to matrix"
     );
 
@@ -93,7 +93,7 @@ main (void)
         }
 
         auto q = quaternionf::identity ();
-        auto m = util::matrix4f::IDENTITY;
+        auto m = util::matrix4f::identity ();
 
         for (auto r: ROTATIONS) {
             q = q.angle_axis (r.mag, r.axis) * q;