Merge pull request #4346 from jwpeterson/is_hpd

Add TypeTensor check for Hermitian positive-definiteness

Author: jwpeterson

include/numerics/type_tensor.h

@@ -436,6 +436,13 @@ class TypeTensor
   void write_unformatted (std::ostream & out_stream,
                           const bool newline = true) const;
 
+  /**
+   * Returns true if the TypeTensor is Hermitian (or symmetric, when
+   * T==Real) and positive-definite to within the provided relative
+   * tolerance, false otherwise.
+   */
+  bool is_hpd(Real rel_tol = TOLERANCE*TOLERANCE) const;
+
 protected:
 
   /**

src/numerics/type_tensor.C

@@ -63,6 +63,61 @@ void TypeTensor<T>::write_unformatted (std::ostream & out_stream,
     out_stream << '\n';
 }
 
+template <typename T>
+bool TypeTensor<T>::is_hpd(Real rel_tol) const
+{
+  // Convenient reference to this object, to be used for matrix
+  // operations.
+  const auto & A = *this;
+
+  // The norm of this matrix, needed for relative tolerance checks.
+  auto A_norm = A.norm();
+
+  // The zero matrix is only positive semi-definite
+  if (A_norm == 0)
+    return false;
+
+  // An absolute tolerance (based on the user's provided relative
+  // tolerance) to be used in the floating point comparisons below.
+  auto abs_tol = rel_tol * A_norm;
+
+  // Form the complex conjugate transpose of A
+  TypeTensor<T> A_conjugate_transpose;
+  for (unsigned int i=0; i<LIBMESH_DIM; i++)
+    for (unsigned int j=0; j<LIBMESH_DIM; j++)
+      A_conjugate_transpose(i,j) = libmesh_conj(A(j,i));
+
+  // Check if Hermitian
+  if ((A - A_conjugate_transpose).norm() > abs_tol)
+    return false;
+
+  // If we made it here, then we are Hermitian, so now we just need to
+  // check if we are positive-definite by checking that all principal
+  // minors are positive. Note: the determinant of a Hermitian matrix
+  // and all principal minors are real-valued. Since we already
+  // checked that the matrix is Hermitian above, we don't bother to
+  // check the complex parts are also zero now.
+
+  // For 3x3 and 2x2, check the 1x1 determinant
+#if LIBMESH_DIM > 1
+  if (std::real(A(0,0)) < -abs_tol)
+    return false;
+#endif
+
+  // For 3x3, check the upper 2x2 determinant
+#if LIBMESH_DIM > 2
+  if (std::real(A(0,0)*A(1,1) - A(0,1)*A(1,0)) < -abs_tol)
+    return false;
+#endif
+
+  // Finally, check the full determinant
+  if (std::real(A.det()) < -abs_tol)
+    return false;
+
+  // If we made it here, then the matrix is Hermitian
+  // positive-definite.
+  return true;
+}
 
 
 template <>

tests/numerics/type_tensor_test.C

@@ -24,6 +24,7 @@ public:
 #endif
   CPPUNIT_TEST(testRowCol);
   CPPUNIT_TEST(testIsZero);
+  CPPUNIT_TEST(testIsHPD);
 #ifdef LIBMESH_HAVE_METAPHYSICL
   CPPUNIT_TEST(testReplaceAlgebraicType);
 #endif
@@ -104,6 +105,87 @@ private:
     }
   }
 
+  void testIsHPD()
+  {
+    LOG_UNIT_TEST;
+
+    {
+#if LIBMESH_DIM == 3
+      {
+        TensorValue<double> tensor(2., 1., 0.,
+                                   1., 2., 1.,
+                                   0., 1., 2.);
+
+        CPPUNIT_ASSERT(tensor.is_hpd(/*rel_tol=*/0.));
+      }
+      {
+        // Symmetric but not positive-definite, because the upper 2x2
+        // principal minor is zero.
+        TensorValue<double> tensor(1,  0., 0.,
+                                   0., 0., 1.,
+                                   0., 1., 0.);
+
+        CPPUNIT_ASSERT(!tensor.is_hpd());
+      }
+      {
+        // Random matrix that is SPD by construction
+
+        auto get_random = [&]()
+        {
+          return static_cast<Real>(std::rand()) / RAND_MAX; // in range [0,1]
+        };
+
+        TensorValue<double> tensor(get_random(), get_random(), get_random(),
+                                   get_random(), get_random(), get_random(),
+                                   get_random(), get_random(), get_random());
+
+        // Make symmetric
+        tensor = 0.5*(tensor + tensor.transpose());
+
+        // Make positive definite
+        for (auto i : make_range(LIBMESH_DIM))
+          tensor(i,i) += 3.;
+
+        CPPUNIT_ASSERT(tensor.is_hpd());
+      }
+
+#ifdef LIBMESH_USE_COMPLEX_NUMBERS
+
+      auto i = Complex(0, 1);
+      auto one = Complex(1, 0);
+      {
+        // Symmetric but not Hermitian, because Hermitian matrices must
+        // have real-valued entries on the diagonal.
+        TensorValue<Complex> tensor(i,  0., 0.,
+                                    0., 1., 0.,
+                                    0., 0., 1.);
+
+        CPPUNIT_ASSERT(!tensor.is_hpd());
+      }
+      {
+        // Symmetric but not Hermitian, because the off-diagonal
+        // entries must be complex conjugates of one another.
+        TensorValue<Complex> tensor(2,  i,  0.,
+                                    i,  2., i,
+                                    0., i,  2.);
+
+        CPPUNIT_ASSERT(!tensor.is_hpd());
+      }
+      {
+        // Both Hermitian and positive-definite
+        TensorValue<Complex> tensor(2.,     one+i,  0.,
+                                    one-i,  2.,     one+i,
+                                    0.,     one-i,  2.);
+
+        CPPUNIT_ASSERT(tensor.is_hpd());
+      }
+
+#endif
+
+#endif // LIBMESH_DIM == 3
+    }
+  }
+
   void testRotation()
   {
     LOG_UNIT_TEST;