# Low precision from function returning matrix

Got a function depending on single parameter "t" defined below which returns a 2 by 2 matrix

H[t_] := PauliMatrix[1] (3 - Cos[Pi t]) + PauliMatrix[3];
U[t_] := Total[Map[With[{x = Normalize[#]}, TensorProduct[x, x /. t -> 0]] &, Eigenvectors[H[t]]], 1];


When I call the function the result gets rounded to the identity

 U[0.2]
{{1., 0.}, {0., 1.}}


But if I use a replacement I get the answer to a suitable precision

 U[t] /. t -> 0.2
{{0.999843, -0.0177345}, {0.0177345, 0.999843}}


How do I change this behaviour?

This is not a precision problem. You can check that with

U[2/10] - U[0.2] // N


{{1.11022*10^-16, 5.55112*10^-17}, {5.55112*10^-17, 1.66533*10^-16}}

Eigenvectors simply does not (ortho-)normalize eigenvectors for symbolic input. Mostly because this would blow-up the complexity of the returned expressions even more.

• Im sorry I don't see what your example illustrates. Why does U[2/10] give anything different to U[0.2]? I know Eigenvectors doesn't orthonormalise for symbolic input, which is why I'm forcing it to normalise the symbolic expression? If I evaluate the RHS side with "t' left symbolic, and then define the function from that expression it behaves as I want it to. How do you suggest I write the function? – oweydd Nov 21 '18 at 13:10
• Up to machine precision, U[2/10] does not give anything different from U[0.2]. – Henrik Schumacher Nov 21 '18 at 13:24
• I don't see how you orthonormalize the eigenvectors. Did it occur to you that the replacement rule t -> 0 does nothing if you call U with anything else but t? Better use Orthogonalize. – Henrik Schumacher Nov 21 '18 at 13:24
• The eigenvectors are automatically orthogonal, so I don't need to worry about that. The replacement rule does in fact work if I call U with any symbolic input. None of this answers my question. – oweydd Nov 21 '18 at 13:37
• Yes you are right: The replacement rule does something also for other symbols. But it does not what you expect for numerical input. Here is something to mull over: f[t_] := {1, 2, 3, t} /. t -> 0; f /@ Range[3] f /@ N@Range[3]. – Henrik Schumacher Nov 21 '18 at 13:48

The problem with the original code is an evaluation problem. @Henrik points this out in a comment.

The rule t -> 0 becomes 0.2 -> 0 in the call U[0.2]:

Hold[U[0.2]] /. DownValues@U
(*
Hold[Total[(With[{x$$= Normalize[#1]}, x$$ \[TensorProduct] (x\$ /. 0.2 -> 0)] &) /@ Eigenvectors[H[0.2]], 1]]
*)


However the rule is applied to the eigenvectors of H[0.2], which do not contain a 0.2 to replace:

Eigenvectors[H[0.2]]
(*  {{0.540734, -0.841193}, {-0.841193, -0.540734}}  *)


It's different in the call U[t], in which the rule is t -> 0 and the eigenvectors of H[t] contain a t to replace with zero.

One way around the difficulty is to use a different symbol for t in H[t]. Formal symbols cannot normally be assigned a value, so \[FormalT] is a good candidate:

UU[t_] := With[{v = Normalize /@ Eigenvectors[H[\[FormalT]]]},
TensorProduct, {v /. \[FormalT] -> t, v /. \[FormalT] -> 0}], 1]
];


Since Eigenvectors[H[t]] is always the same, one might want to precompute the formulas for U[t] and use them to define the function U. It's a good idea to protect t with Block, in case it has been assigned a value. Also, it turns out formal symbols can be assigned a value by a user with Block[{\[FormalT] = 5}, <code>], so it would be even safer to use Module to create a temporary symbol for H[t]:

Block[{t},
UU[t_] = Module[{$$t, v}, v = Normalize /@ Eigenvectors[H[$$t]];
Total[MapThread[TensorProduct, {v /. $$t -> t, v /.$$t -> 0}], 1]
]
];


Defining the function as

 U[t_] := Evaluate[Total[Map[With[{x = Normalize[#]}, TensorProduct[x,x /. t -> 0]] &, Eigenvectors[H[t]]], 1]];


gives the behaviour I want. But I still don't understand why the original definition did't work.

• If this resolves your problem, then you should learn about the difference between Set (=) and SetDelayed (:=)... – Henrik Schumacher Nov 21 '18 at 13:45