Merge pull request #397 from NumSim-Stack/389-390-gram-and-sympart

#389 #390: derive symmetric(+PSD) for Gram forms and the symmetric part

Author: petlenz

include/numsim_cas/tensor/tensor_operators.h

@@ -83,6 +83,19 @@ inline expression_holder<tensor_expression> tag_invoke(add_fn, L &&lhs,
         result.data()->set_space({Skew{}, AnyTraceTag{}});
       return result;
     }
+    // trans(A) + A and A + trans(A) are symmetric in any dimension: the
+    // symmetric part 2*sym(A). Mirror of the skew branch above (and of the
+    // sub operator's skew rule). The skew branch matched trans(A) + (-A)
+    // above, so a match here means a genuine plus. Closes #390.
+    if (is_trans_of(lhs, rhs) || is_trans_of(rhs, lhs)) {
+      auto &_lhs{lhs.template get<tensor_visitable_t>()};
+      simplifier::tensor_detail::add_base visitor(std::forward<L>(lhs),
+                                                  std::forward<R>(rhs));
+      auto result = _lhs.accept(visitor);
+      if (result.is_valid())
+        result.data()->set_space({Symmetric{}, AnyTraceTag{}});
+      return result;
+    }
   }
 
   auto &_lhs{lhs.template get<tensor_visitable_t>()};
@@ -166,6 +179,22 @@ tag_invoke(mul_fn, L &&lhs, [[maybe_unused]] R &&rhs) {
       return make_expression<identity_tensor>(lhs.get().dim(), std::size_t{2});
     if (is_trans_of(lhs, rhs) && is_orthogonal(rhs))
       return make_expression<identity_tensor>(rhs.get().dim(), std::size_t{2});
+    // Gram form: trans(X)*X and X*trans(X) are symmetric and positive
+    // semidefinite for any X (x . (F^T F) . x = ||F x||^2 >= 0). Orthogonal X
+    // already returned I above; X invertibility is unknown, so PSD only (not
+    // PD). Closes #389.
+    if (is_trans_of(lhs, rhs) || is_trans_of(rhs, lhs)) {
+      auto &_lhs{lhs.template get<tensor_visitable_t>()};
+      tensor_detail::simplifier::mul_base visitor(std::forward<L>(lhs),
+                                                  std::forward<R>(rhs));
+      auto result = _lhs.accept(visitor);
+      if (result.is_valid()) {
+        result.data()->set_space({Symmetric{}, AnyTraceTag{}});
+        result.data()->tensor_algebra_assumptions().insert(
+            positive_semidefinite{});
+      }
+      return result;
+    }
   }
   auto &_lhs{lhs.template get<tensor_visitable_t>()};
   tensor_detail::simplifier::mul_base visitor(std::forward<L>(lhs),

tests/TensorSpacePropagationTest.h

@@ -424,6 +424,75 @@ TEST_F(TensorSpacePropagationTest, DiffAddSymmetric) {
   EXPECT_PRINT(d, "2*P_sym{4}");
 }
 
+// ═══════════════════════════════════════════════════════════════════════════════
+// Gram forms (#389): trans(X)*X and X*trans(X) are symmetric + PSD for any X
+// ═══════════════════════════════════════════════════════════════════════════════
+
+TEST_F(TensorSpacePropagationTest, GramTransXtimesXIsSymmetricPSD) {
+  // trans(X)*X = right Cauchy-Green; symmetric and PSD regardless of X.
+  auto g = trans(X) * X;
+  EXPECT_TRUE(is_symmetric(g)) << "trans(X)*X must be symmetric";
+  EXPECT_TRUE(is_positive_semidefinite(g)) << "trans(X)*X must be PSD";
+  EXPECT_FALSE(is_skew(g));
+  EXPECT_FALSE(is_positive_definite(g)) << "X invertibility unknown → PSD only";
+}
+
+TEST_F(TensorSpacePropagationTest, GramXtimesTransXIsSymmetricPSD) {
+  // X*trans(X) = left Cauchy-Green / Finger tensor.
+  auto g = X * trans(X);
+  EXPECT_TRUE(is_symmetric(g));
+  EXPECT_TRUE(is_positive_semidefinite(g));
+}
+
+TEST_F(TensorSpacePropagationTest, GramNormalizesUnderSymProjector) {
+  // The derived symmetric annotation must let sym() short-circuit.
+  EXPECT_PRINT(sym(trans(X) * X), ::testcas::S(trans(X) * X));
+}
+
+TEST_F(TensorSpacePropagationTest, OrthogonalTransXtimesXStillFoldsToIdentity) {
+  // Safety: for orthogonal Q the trans(Q)*Q -> I fold must still win over the
+  // Gram annotation branch (I is stronger than "some symmetric PSD tensor").
+  auto Q = std::get<0>(make_tensor_variable(std::tuple{"Q", dim, 2}));
+  assume_orthogonal(Q);
+  EXPECT_TRUE(is_same<identity_tensor>(trans(Q) * Q));
+  EXPECT_TRUE(is_same<identity_tensor>(Q * trans(Q)));
+}
+
+TEST_F(TensorSpacePropagationTest, DistinctSymmetricProductStaysNonSymmetric) {
+  // Safety: A*B of two DISTINCT symmetric tensors is symmetric only if they
+  // commute — the Gram rule must not over-generalize to any product.
+  auto B = std::get<0>(make_tensor_variable(std::tuple{"B", dim, 2}));
+  assume_symmetric(B);
+  EXPECT_FALSE(is_symmetric(C * B));
+  EXPECT_FALSE(is_positive_semidefinite(C * B));
+}
+
+// ═══════════════════════════════════════════════════════════════════════════════
+// Symmetric part (#390): trans(X)+X and X+trans(X) are symmetric for any X
+// ═══════════════════════════════════════════════════════════════════════════════
+
+TEST_F(TensorSpacePropagationTest, SymmetricPartXPlusTransXIsSymmetric) {
+  // X + trans(X) = 2*sym(X); symmetric in any dimension.
+  EXPECT_TRUE(is_symmetric(X + trans(X)));
+  EXPECT_TRUE(is_symmetric(trans(X) + X));
+  EXPECT_FALSE(is_skew(X + trans(X)));
+}
+
+TEST_F(TensorSpacePropagationTest, SymmetricPartDoesNotStealTheSkewSpelling) {
+  // Safety: trans(X)-X and trans(X)+(-X) must remain skew — the skew branch
+  // is checked before the symmetric-part branch.
+  EXPECT_TRUE(is_skew(trans(X) - X));
+  EXPECT_TRUE(is_skew(trans(X) + (-X)));
+  EXPECT_TRUE(is_skew(X - trans(X)));
+}
+
+TEST_F(TensorSpacePropagationTest, GeneralSumStaysUnannotated) {
+  // Safety: a generic X+Y (Y not trans(X)) carries no space.
+  auto Y = std::get<0>(make_tensor_variable(std::tuple{"Y", dim, 2}));
+  EXPECT_FALSE(is_symmetric(X + Y));
+  EXPECT_FALSE(is_skew(X + Y));
+}
+
 } // namespace numsim::cas
 
 #endif // TENSORSPACEPROPAGATIONTEST_H