# How to define a matrix product based on a different scalar product

I want to define a product of 3x3 Matrices where the coefficient behave in a different way as usual. I defined the multiplication between coefficient as the function MyScalProd

MyScalProd[x_, y_] := (x^2*y + y^2*x)/2; (*scalar product*)


and now I want to use it to multiply coefficients involved in the matrix product. I tried as I could this way

ProdMat[z_, a_, b_] := (
z = Array[cc, {3, 3}];
cc[1, 1] =
MyScalProd[a[[1, 1]], b[[1, 1]]] +
MyScalProd[a[[1, 2]], b[[2, 1]]] +
MyScalProd[a[[1, 3]], b[[3, 1]]];
cc[1, 2] =
MyScalProd[a[[1, 1]], b[[1, 2]]] +
MyScalProd[a[[1, 2]], b[[2, 2]]] +
MyScalProd[a[[1, 3]], b[[3, 2]]];
cc[1, 3] =
MyScalProd[a[[1, 1]], b[[1, 3]]] +
MyScalProd[a[[1, 2]], b[[2, 3]]] +
MyScalProd[a[[1, 3]], b[[3, 3]]];
cc[2, 1] =
MyScalProd[a[[2, 1]], b[[1, 1]]] +
MyScalProd[a[[2, 2]], b[[2, 1]]] +
MyScalProd[a[[2, 3]], b[[3, 1]]];
cc[2, 2] =
MyScalProd[a[[2, 1]], b[[1, 2]]] +
MyScalProd[a[[2, 2]], b[[2, 2]]] +
MyScalProd[a[[2, 3]], b[[3, 2]]];
cc[2, 3] =
MyScalProd[a[[2, 1]], b[[1, 3]]] +
MyScalProd[a[[2, 2]], b[[2, 3]]] +
MyScalProd[a[[2, 3]], b[[3, 3]]];
cc[3, 1] =
MyScalProd[a[[3, 1]], b[[1, 1]]] +
MyScalProd[a[[3, 2]], b[[2, 1]]] +
MyScalProd[a[[3, 3]], b[[3, 1]]];
cc[3, 2] =
MyScalProd[a[[3, 1]], b[[1, 2]]] +
MyScalProd[a[[3, 2]], b[[2, 2]]] +
MyScalProd[a[[3, 3]], b[[3, 2]]];
cc[3, 3] =
MyScalProd[a[[3, 1]], b[[1, 3]]] +
MyScalProd[a[[3, 2]], b[[2, 3]]] +
MyScalProd[a[[3, 3]], b[[3, 3]]];
)(*MatrixProduct based on MyScalProd*)

1. I would like to have a function so that given two matrices I can write "Z=ProdMat[A,B]"
2. the above code seem to work only the first time I call the procedure, while the second time gives problems involving the reuse of the old variables. Let me know if you have solutions
• While off-hand I see several issues, it is important to ask: what is the syntax with which you call this function, such that you have success on the first call, but subsequent calls result in errors, as you mentioned? Commented Mar 11, 2022 at 6:25
• Is something like ProdMat[A1, X, Y]; ProdMat[A2, Y, X] ... the result is something like this gyazo.com/b571084a443d77c4b1eacdd417f5e5bb
– Dac0
Commented Mar 11, 2022 at 6:37

The reason why this works only once is that

• you define the matrix z symbolically in terms of the symbols c[i,j]
• then you define the symbols c[i,j].
• but the symbol c is not scoped.

When you call ProdMat the second time then c[i,j] is already defined; so z = Array[cc, {3, 3}]; writes directly to the variable z. This is why later changes to c cannot affect z.

Also you use z as input variable. Mathematica does not allow to modify the inputs unless ProdMat has one of the attributes HoldFirst or HoldAll.

A short and more flexible approach would be something like

ProdMat[a_, b_] := Table[
Sum[MyScalProd[a[[i, k]], b[[k, j]]],{k,1,Dimensions[a][[2]]}]
,{i,1,Dimensions[a][[1]]}
,{j,1,Dimensions[b][[2]]}
]


A repaired version of your original code look like this:

ProdMatRepaired[a_, b_] := Module[{z},
z = ConstantArray[0, {Dimensions[a][[1]], Dimensions[b][[2]]}];
Do[
z[[i, j]] = Sum[
MyScalProd[a[[i, k]], b[[k, j]]]
, {k, 1, Dimensions[a][[2]]}
]
, {i, 1, Dimensions[a][[1]]}, {j, 1, Dimensions[b][[2]]}];
z
]


Note how I used Module to scope the variable z and how I set the return value of the function to z by placing z at the end of the Module without a semicolon.

Now you can do

a = RandomInteger[{1, 10}, {3, 3}];
b = RandomInteger[{1, 10}, {3, 3}];
res1 = ProdMat[a, b];
res2 = ProdMatFixed[a, b];


Both lead to the same result:

res1 == res2


True

• just to be sure I understood the problem: if I had put "z = Array[cc, {3, 3}];" after the definition of the symbols c[i,j] it would have worked?
– Dac0
Commented Mar 11, 2022 at 10:41
• Yes, it would. But it would have been more robust (and probably more efficent) to create an empty matrix z = ConstantArra[0,{3,3}] and then write z[[i,j]] = ... instead of cc[2, 3] = .... Commented Mar 11, 2022 at 10:43
• and so I would have to write in the syntax: ProdMat[A,B]; res=z; t have the result stored in matrix "res"?
– Dac0
Commented Mar 11, 2022 at 10:45
• No. I would return z as the result of ProdMat so that one can simply write res = ProdMat[A,B]. Like I did in my answer. Commented Mar 11, 2022 at 10:46
• Sorry I didn't understand... but it would be very helpful for me. Could you please rewrite in your answer the fix you would have done to the original code? Thank you very much
– Dac0
Commented Mar 11, 2022 at 10:48

The function Inner allows you to define an inner product, where you can specify the functions for multiplication and addition. So you can define for example a new operator

CircleDot[x_List,y_List]:=Inner[MyScalProd,x,y,Plus];
CircleDot[x_List,y_List,z__]:=CircleDot[CircleDot[x,y],z];


which leads the same result as the implementation of Henrik Schumacher

a ⊙ b === ProdMat[a,b]
(* True *)

b ⊙ a ⊙ b
(* {{4799469, 4613164, 4067140}, {824411, 803114, 681023}, {3772584,  4714759, 1968481}} *)

• This is a neat implementation. I don't follow why matrix multiplication is an inner product though?
– Jojo
Commented Mar 11, 2022 at 17:54
• @Joe Yes you are right, it's not an inner product for matrices. What I really meant is that Inner is a generalization of Dot, which acts an inner product for vectors. Commented Mar 11, 2022 at 18:39