Subtituting kronecker deltas for products in simplification

I have a problem where I have some vector, $$a = (a_1, a_2, \ldots, a_n)^T$$ that I want to have a Kronecker delta multiplication rule, $$a[i]a[j] = \delta_{ij}$$ (these are basically basis vectors, but for my purposes it's useful to just treat them as symbols that have that product relationship). Then I have some $$n \times n$$ matrices $$B$$ and $$C$$ and I want to compute $$(Ba)^T(Ca)$$. Note that this is just a polynomial in the components of the $$a$$ vector, and I'd like to simplify it according to my multiplication rule. For instance,

avec = {a[1], a[2], a[3]}
Bmat = {{B[1, 1], B[1, 2], B[1, 3]}, {B[2, 1], B[2, 2],
B[2, 3]}, {B[3, 1], B[3, 2], B[3, 3]}}
Cmat = {{C[1, 1], C[1, 2], C[1, 3]}, {C[2, 1], C[2, 2],
C[2, 3]}, {C[3, 1], C[3, 2], C[3, 3]}}

Transpose[(Bmat . avec)] . (Cmat . avec) // Expand


gives

a[1]^2 B[1, 1] C[1, 1] + a[1] a[2] B[1, 2] C[1, 1] +
a[1] a[3] B[1, 3] C[1, 1] + a[1] a[2] B[1, 1] C[1, 2] +
a[2]^2 B[1, 2] C[1, 2] + a[2] a[3] B[1, 3] C[1, 2] +
a[1] a[3] B[1, 1] C[1, 3] + a[2] a[3] B[1, 2] C[1, 3] +
a[3]^2 B[1, 3] C[1, 3] + a[1]^2 B[2, 1] C[2, 1] +
a[1] a[2] B[2, 2] C[2, 1] + a[1] a[3] B[2, 3] C[2, 1] +
a[1] a[2] B[2, 1] C[2, 2] + a[2]^2 B[2, 2] C[2, 2] +
a[2] a[3] B[2, 3] C[2, 2] + a[1] a[3] B[2, 1] C[2, 3] +
a[2] a[3] B[2, 2] C[2, 3] + a[3]^2 B[2, 3] C[2, 3] +
a[1]^2 B[3, 1] C[3, 1] + a[1] a[2] B[3, 2] C[3, 1] +
a[1] a[3] B[3, 3] C[3, 1] + a[1] a[2] B[3, 1] C[3, 2] +
a[2]^2 B[3, 2] C[3, 2] + a[2] a[3] B[3, 3] C[3, 2] +
a[1] a[3] B[3, 1] C[3, 3] + a[2] a[3] B[3, 2] C[3, 3] +
a[3]^2 B[3, 3] C[3, 3]


and I'd like it to just automatically set all those cross terms to 0.

I just literally don't even know where to begin. The only think I can figure out how to do is add an assumption that captures the non-crossed terms,

FullSimplify[Transpose[(Bmat . avec)] . (Cmat . avec) // Expand,
a[_]^2 == 1]


but that leaves all the cross terms. I've tried the only real thing I can think of,

FullSimplify[Transpose[(Bmat . avec)] . (Cmat . avec) // Expand,
a[_] a[_] == KroneckerDelta[_, _]]


but that doesn't seem to work. Any help is appreciated!

You can use ReplaceAll. Also, avoid capitalizing your symbols.

aVec = Array[a, 3];
bMat = Array[b, {3, 3}];
cMat = Array[c, {3, 3}];

rules = {a[i_] a[j_] :> KroneckerDelta[i, j], a[_]^2 :> 1};

res = (Transpose[(bMat . aVec)] . (cMat . aVec) // Expand) /. rules


b[1, 1] c[1, 1] + b[1, 2] c[1, 2] + b[1, 3] c[1, 3] + b[2, 1] c[2, 1] + b[2, 2] c[2, 2] + b[2, 3] c[2, 3] + b[3, 1] c[3, 1] + b[3, 2] c[3, 2] + b[3, 3] c[3, 3]

• That's correct (+1), @Domen. Commented Mar 10, 2023 at 22:31

A variation of @Domen's answer is as follows:

expr = Transpose[(Bmat . avec)] . (Cmat . avec) // Expand;

rule1 = a[i_]^2 :> HoldForm[a[i]*a[i]];
rule2 = a[i_] *a[j_] :> HoldForm[KroneckerDelta[i, j]];

expr2 = expr /. rule1 /. rule2


Then using ReleaseHold:

Apply[ReleaseHold@*ReleaseHold, {expr2}]


Remark: This way of applying the rules allows us to observe that the replacements are made correctly.

Or equivalently:

expr /. {a[i_]*a[j_] :> KroneckerDelta[i, j], a[i_]*a[i_] :> KroneckerDelta[i, i]}

• This is great and seems to work. Could you provide some insight into what's going on? Commented Mar 10, 2023 at 22:06
• Does this really work? Am I completely confused or does it just removes all occurences of a ... Commented Mar 10, 2023 at 22:11
• @Domen True, it is a bit confusing. Maybe it would be better to set a criterion to identify the pattern using If and MatchQ. Commented Mar 10, 2023 at 22:15
• Thanks for correcting me, @Domen :-) Commented Mar 10, 2023 at 22:46