Accuracy limitations of singular value decomposition?

Question

in the process of working on a physics problem I have found the need to use the singular value decomposition function built into Mathematica. I have encountered what seem to be limitations to the accuracy of the SVD when the elements of the matrix we are decomposing become large. Let me be more explicit:

The matrix to decompose is defined in terms of the tensor which can be constructed by running the following code:

T = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, 3}];

Part[Part[Part[T, 1], 1], 1] = e^(3 (KK + h));
Part[Part[Part[T, 1], 3], 3] = 3 e^(-KK + h);
Part[Part[Part[T, 2], 2], 2] = e^(3 (KK - h));
Part[Part[Part[T, 2], 3], 3] = 3 e^-(KK + h);
T = Normal[Symmetrize[T]];

The matrix that we want to use the SVD on is then

Partition[Flatten[TensorContract[TensorProduct[T,T],{3,4}]],9]

Now we just have to choose some values for the parameters $KK, h$ and let the SVD do its magic. Since it might be relevant to a diagnosis, here is my code for running the SVD:

singout = SingularValueDecomposition[
N[Partition[Flatten[TensorContract[TensorProduct[T, T], {3, 4}]], 
9]], Tolerance -> 0]

In the context of my physics problem, it's reasonable to take $h=0$, so let's do that. If we run the SVD with $KK=1$, I find that the code works exactly as intended: the input matrix is reproduced almost exactly by the matrix product

u.w.Transpose[v]

where $u, w$, and $v$ are as defined in the SVD documentation. There is a bit of numerical noise, but nothing goes dramatically wrong.

In contrast, we can take $KK=10$, which is, again, a reasonable idea in the context of my physics problem. I'm not sure how to format matrices here, so I've attached a picture of the output. The first matrix is the exact input, and the second is the SVD output, chopped at $10^{-10}$ to avoid clutter. There are some small discrepancies between the two matrices, but by far the most important error is the bottom-left entry of the SVD output: it is returning zero instead of a number on the order of $10^8$. I need to understand why this is happening so that I can control the approximation from the SVD appropriately. Any suggestions for understanding this would be very helpful.

I should perhaps also note that this seems to be a problem in general for large values of $KK$, but I'm not sure why that is the case. Thanks for any help you can offer.

Answer 1

ClearAll["Global`*"]
T = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, 3}];

Part[Part[Part[T, 1], 1], 1] = E^(3 (KK + h));
Part[Part[Part[T, 1], 3], 3] = 3 E^(-KK + h);
Part[Part[Part[T, 2], 2], 2] = E^(3 (KK - h));
Part[Part[Part[T, 2], 3], 3] = 3 E^-(KK + h);
T = Normal[Symmetrize[T]];

h = 0;
KK = 10;
im = Partition[Flatten[TensorContract[TensorProduct[T, T], {3, 4}]], 
   9];
singout = SingularValueDecomposition[N[im, 20], Tolerance -> 0];

{u, w, v} = singout;
res = (u.w.Transpose[v]);

N@MinMax[res - im]

{-4.12231*10^-9, 3.54693*10^-41}

As you can see, very small difference... and as noted, chopping and comparing visually it matches your image.