1

I have two matrices I want to add, and one of the matrices is a tensor product of two vectors. I've used a SetDelayed to define the summed matrix, because I want to evaluate it for different values of these vectors. However, when I try to evaluate it using a replacement, I get a new array with the wrong dimensionality. Here's a minimal(ish?) example:

mat := DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec]

Then using a replace gives me something horrible:

mat /. vec -> {2, 2, 2, 2}

{{{{5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 5}}, {{4, 4, 4,
 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}}, {{4, 4, 4, 4}, {4,
 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}}, {{4, 4, 4, 4}, {4, 4, 4, 
4}, {4, 4, 4, 4}, {4, 4, 4, 4}}}, {{{4, 4, 4, 4}, {4, 4, 4, 
4}, {4, 4, 4, 4}, {4, 4, 4, 4}}, {{5, 5, 5, 5}, {5, 5, 5, 5}, {5, 
5, 5, 5}, {5, 5, 5, 5}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}, {4, 4, 4, 4}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 
4, 4, 4}}}, {{{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}, {{5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 
5}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}}, {{{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}, {{4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 4}, {4, 4, 4, 
4}}, {{5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 5}, {5, 5, 5, 5}}}}

But a block (or setting vec={2,2,2,2}) gives the correct behavior:

Block[{vec = {2, 2, 2, 2}}, mat]

{{5, 4, 4, 4}, {4, 5, 4, 4}, {4, 4, 5, 4}, {4, 4, 4, 5}}

(Sorry, not sure if I typed in the output correctly here.)

What's going on?

Adam
  • 709
  • 6
  • 11
  • The way to debug this type of problem is to evaluate your expression (DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec]) symbolically. Try that and you'll see exactly what's going wrong. – Szabolcs Jun 02 '14 at 15:35

1 Answers1

4

The way to debug this type of problem is to evaluate your expression (DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec]) symbolically. Try that and you'll see exactly what's going wrong.

Addition auto-threads inside of arrays, i.e. {1,2,3} + a auto-evaluates to {1 + a, 2 + a, 3 + a}. This is what happens in your case, with the diagonal matrix standing in place of {1,2,3}, and the still symbolic TensorProduct[vec, vec] in place of a.

You could use

DiagonalMatrix[diag] + TensorProduct[vec, vec]  /. {diag -> {1,1,1,1}, vec -> {2,2,2,2}}

or

Hold[DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec]] /. vec -> {2,2,2,2} // ReleaseHold

or

mat[vec_] := DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec]
Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
  • Thanks @Szabolcs! This makes perfect sense. What I've ended up doing is defining the added matrix as above, but then assigning a value to vec before ever evaluating the addition. Is there anything especially wrong or problematic about that solution? – Adam Jun 02 '14 at 15:52
  • @Adam Do you mean doing mat := DiagonalMatrix[{1, 1, 1, 1}] + TensorProduct[vec, vec] then vec = ...; mat? If so, it works, but I'm a bit uncomfortable as it's error-prone. If you set the value of vec once only, then I'd rather do vec = ...; mat = ... (note = vs :=). If you keep changing the value of vec, then I'd do mat[vec_] := ... to make it a parameter. Sometimes I end up with several parameters that originated from global variables. In that case I sometimes use withRules described here. – Szabolcs Jun 02 '14 at 16:23
  • That is essentially an equivalent of the Block solution you show with the difference that it is easy to store parameter values in lists as {vec -> ...} and re-use them later. – Szabolcs Jun 02 '14 at 16:23
  • The real problem is that it only works when you actually can do the replacements before. but sometimes (when generating compilation bodies e.g.) you cannot. I find this a design omittance. See also remarks [here]. (http://chat.stackexchange.com/transcript/message/15676082#15676082) and here It is always a hassle to work around the missing abstract datatype features of Mathematica (and even the Wolfram Language). – Rolf Mertig Jun 02 '14 at 16:25