I think you have two main alternatives:

1) Use a MatrixSymbol object with explicit symmetry in its two levels:

In[1]:= A = MatrixSymbol["A", {n, n}, Symmetric[{1, 2}]]

In[2]:= A + Transpose[A] // TensorReduce
Out[2]= 2 MatrixSymbol["A", {n, n}, Complexes, Symmetric[{1, 2}]]

In[3]:= A - Transpose[A] // TensorReduce
Out[3]= SymbolicZerosArray[{n, n}]

In[4]:= Inverse[A] . Transpose[A] // TensorReduce
Out[4]= MatrixPower[MatrixSymbol["A", {n, n}, Complexes, Symmetric[{1, 2}]], 0]

2) or declare assumptions for a symbol B (I use a separate symbol here to avoid Element resolving to True immediately):

In[5]:= $Assumptions = Element[B, Arrays[{n, n}, Symmetric[{1, 2}]]];

In[6]:= B + Transpose[B] // TensorReduce
Out[6]= 2 B

In[7]:= B - Transpose[B] // TensorReduce
Out[7]= SymbolicZerosArray[{n, n}]

In[8]:= Inverse[B] . Transpose[B] // TensorReduce
Out[8]= MatrixPower[B, 0]

Though correct, MatrixPower[squarematrix, 0] should be converted into IdentityMatrix[n].

I think this type of computation requires guidance, mainly to specify what needs to be replaced by what. For example, I'd do something like this:

$Assumptions = {
   Element[A, Matrices[{n, n}, Symmetric[{1, 2}]]],
   Element[x | y | b | c, Vectors[{n}]]
};

rules = {
   x -> Inverse[A] . b,
   y -> Inverse[A] . c
};

In[3]:= b . y == c . x /. rules // TensorReduce
Out[3]= True

Note that the traditional notation that transposes the vector u on the LHS of a two-vector inner product uT . v can be represented more easily in WL as Dot[u, v].

Thanks Jose, I appreciate it. While this is very useful, the fact that I need to know the result before verifying it is a limiting factor.

I think it would be great if the WL (Full)Simplify and Eliminate commands could be extended to transform matrix statements. Perhaps LLMs can help? Just for the heck, I entered the following prompt in Claude.ai:

Prompt: Consider an nxn matrix A, and n-vectors b, c, x and y. Assume that A.x = b, A.y = c, and Transpose(A) = A. Develop an alternative expression for Transpose[b].y that does not involve b or y, and is as simple as possible, avoiding A if at all possible. *** Response: Let's approach this step-by-step:

We're given that A.x = b This means b = A.x We're also given that A.y = c This means y = A^(-1).c (assuming A is invertible) We're told that Transpose(A) = A This means A is symmetric We need to find an expression for Transpose[b].y Let's start with Transpose[b].y: Transpose[b].y = (A.x)^T.y (substituting from step 1) = x^T.A^T.y (property of matrix transpose) = x^T.A.y (because A is symmetric) = x^T.c (substituting from step 2)

Therefore, Transpose[b].y = x^T.c

Perhaps WL can call an LLM under the hood to produce candidate results that next are (automatically) verified using the 14.1 approach you described?