Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 2 at 15:35

1 Answer 1

up vote 4 down vote accepted

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]
share|improve this answer
    
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 2 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 2 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 2 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]. (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 2 at 16:25

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.