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


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


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.

Pictures of relevant matrices: the first is the exact input and the second is the reconstructed matrix from the singular value decomposition.

  • 1
    $\begingroup$ Change the precision in N to an appropriate value, e.g. 20 gets 4.85165*10^8 in LL corner... $\endgroup$
    – ciao
    Apr 21, 2016 at 8:47
  • $\begingroup$ 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. $\endgroup$
    – miggle
    Apr 21, 2016 at 8:52
  • $\begingroup$ Beats me then - try restarting your kernel - looking at the result with appropriate N setting matches your top (input matrix) exactly when chopped... $\endgroup$
    – ciao
    Apr 21, 2016 at 9:14
  • $\begingroup$ 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... $\endgroup$
    – miggle
    Apr 21, 2016 at 9:25
  • $\begingroup$ I'll post my run as answer (too big for comment). $\endgroup$
    – ciao
    Apr 21, 2016 at 9:33

1 Answer 1

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}]], 
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.

  • 1
    $\begingroup$ 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! $\endgroup$
    – miggle
    Apr 21, 2016 at 18:08
  • $\begingroup$ @user46073 The short answer: 5.2`30 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. $\endgroup$
    – MarcoB
    Apr 21, 2016 at 19:51
  • $\begingroup$ @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. $\endgroup$ Apr 22, 2016 at 0:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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