Are extremely small singular values (less than ~10^-15) trustworthy?

Question

I'm using iTensor to simulate a certain 1D process with MPS. The Schmidt values (i.e. singular values) across any cut for this particular process decay really fast. There are typically on the order of only 5 values that are smaller than 10^-5, and then they quickly decay to values smaller than 10^-25. I'm interested in studying the functional form of the spectrum of these singular values if possible.

I noticed that iTensor seems to have a hardcoded "truncation error" cutoff of 10^-13 if the truncation error parameter is not specified, for example in the .orthogonalize function. However, if the truncation error cutoff parameter is set to 0, then no singular values are automatically truncated. For my particular simulation, singular values as small as ~10^-55 can be recovered if the truncation error is set to zero.

My question is: if an MPS simulation is performed with zero truncation error as above, and the list of Schmidt/singular values across some cut is obtained, can extremely tiny values (smaller than 10^-15 or so) be considered reliable? Or should these extremely small values be considered garbage due to numerical precision issues? I'm wondering if 10^-13 is the default truncation error because anything less should not be trusted.

Thank you!

Answer 1

Good question, though I don't feel comfortable giving a totally absolute answer. It's one of those "it depends" things as with any very technical matter.

But as a rule of thumb, when working with double-precision floating point numbers most operations are only trustworthy to about a precision of 10^-13 or so. This is because, among other reasons, that while double can store numbers much more precise than that, precision losses occur when subtracting two similar numbers for example. It's hard to avoid these kinds of operations in most algorithms of course.

Here is a sample calculation I just did in the Julia REPL to illustrate this:

Start by generating some random numbers:

julia> r1 = rand()
0.034238108102722764

julia> r2 = rand()
0.7802465824048777

julia> r3 = rand()
0.09507249159735309

Now do some arithmetic and then undo it

julia> x = (r1/r2)*r3-100
-99.9958281108584

julia> y = (100+x)/r3*r2 # y should equal r1 if using exact arithmetic
0.03423810810274365

julia> r1
0.034238108102722764 # you can see the last few digits differ from y

julia> y-r1
2.0886070650760757e-14