# Accuracy limitations of singular value decomposition?

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.

• Change the precision in N to an appropriate value, e.g. 20 gets 4.85165*10^8 in LL corner... – ciao Apr 21 '16 at 8:47
• Really? I just tried to reproduce that on my machine and there was no change. I increased the working/machine precision and purposely left the second argument for N empty to avoid that problem, but putting it in manually changes nothing for me. – miggle Apr 21 '16 at 8:52
• Beats me then - try restarting your kernel - looking at the result with appropriate N setting matches your top (input matrix) exactly when chopped... – ciao Apr 21 '16 at 9:14
• Hmm, that's bizarre. I closed everything completely and didn't change any internal precision settings, just adjusted the argument of N, and that didn't help me at all. What is happening... – miggle Apr 21 '16 at 9:25
• I'll post my run as answer (too big for comment). – ciao Apr 21 '16 at 9:33

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.

• I see, thanks for posting the code. This works correctly on my machine - I realized that the only difference was I defined KK=10.0 instead of KK=10. I didn't realize that the decimal point could change things so dramatically. Incidentally, do you know if there is a way to get around that? In general I will want to use non-integer values of KK, so I'll have to deal with this as well. Thanks again! – miggle Apr 21 '16 at 18:08
• @user46073 The short answer: 5.230 indicates the value $5.2$ expressed with 30 digits of precision. Alternatively you can use e.g. Rationalize[number, 0] to get an arbitrary (i.e. infinite) precision rational representation of number. Take a look at Control the Precision and Accuracy of Numerical Results and Numerical Precision Tutorial. – MarcoB Apr 21 '16 at 19:51
• @mflynn, in general, any number with a decimal point in Mathematica is inexact, and unless you specify the precision like in Marco's comment, the precision defaults to machine precision. Which will be fast, but not necessarily numerically reliable. – J. M.'s torpor Apr 22 '16 at 0:58