Raise the default of lambda_cutoff from 1e-12 to 1e-7

[aoymt]

Summary

The simple update gave wrong results, without any warning, where a bond had fewer nonzero Schmidt values than its dimension. The default of [parameter.simple_update] lambda_cutoff was too small to do what it is there for.

Case	Exact	Before	After
Decoupled dimers at beta = 1, D = 4 on every bond	-0.124700	+0.391957	-0.124700
Ising antiferromagnet in a staggered field, honeycomb lattice, one bond of each site along the diagonal	-0.625	reached by 6 of 16 seeds	16 of 16
The same model with nearest-neighbour bonds only	-0.625	12 of 16	15 of 16

The remaining seed of the last row ends in the inverted Neel state, which is a local minimum of a classical model.

Cause

The simple update multiplies a tensor by the mean field lambda of its other bonds, decomposes it, and divides by lambda again. lambda is the square root of a Schmidt value, and both ends of a bond carry one factor of it.

A component whose Schmidt value is below the machine epsilon, relative to the largest, is lost in the round-off of the decomposition once it is multiplied by lambda. The division then returns an error of eps / lambda in place of a component of size lambda. The error is larger than the component for lambda below sqrt(eps) = 1.5e-8. It enters the next decomposition of that bond as a Schmidt value, and grows.

lambda_cutoff takes such a lambda for zero. At 1e-12 it corresponds to a Schmidt value of 1e-24, which an SVD in double precision does not return for anything but an exact zero.

The two cases

• Finite temperature, a bond that no term of the Hamiltonian uses. The run starts from the identity, with one component on every bond. This case appeared with 7a4fb38, which puts an identity gate on such bonds: an executable built before that commit gives the same +0.391957 when the six identity gates are written into the input by hand.
• Ground state of a Hamiltonian diagonal in the site basis. Its ground state is a product state, so every bond ends at Schmidt rank one. It looked like a defect of the gate chain that carries a non-nearest-neighbour bond, since the chain is where most of the round-off comes from; the gates of a nearest-neighbour bond of such a model are diagonal and keep exact zeros exact.

Choice of the value

lambda_cutoff	Honeycomb, seeds reaching -0.625	Dimers
1e-12 (former default)	6 of 16	wrong
1e-11		wrong
1e-10	15 of 16	right
1e-9	16 of 16	wrong
1e-8	16 of 16	right
1e-7 (new default)	16 of 16	right
1e-6	16 of 16	right

Below the threshold the result is right or wrong by chance. 1e-7 leaves a factor of seven to sqrt(eps).

With the new default the honeycomb model also reaches -0.625 for 6 of 6 seeds at D = 3 and at D = 4, where 3 and 2 of 6 did before.

The tables were measured with a build of a4ad153, by setting lambda_cutoff in the input. The new test below shows the same failures with the current develop.

Behavior changes

• The default of lambda_cutoff is 1e-7. Components with a Schmidt value below about 1e-14 of the largest are dropped in the simple update.
• The reference results of the test suite do not change.
• lambda_cutoff = 1e-12 in the input gives the former behavior.

Note for reviewers: two existing tests are changed

test/fermion/longrange_gate.cpp and test/fermion/unit_d1.cpp require every lambda after an exact gate to be zero, more than a thousand times the cutoff, or below 1e-8 of the largest.

• The zeros of the SVD come out at 9.8e-9, 9.4e-9 and the like, just above that 1e-8, and six orders below the smallest kept value (0.03).
• With the former default they passed as "more than a thousand times 1e-12", that is, as kept values.
• With the new default they are what the cutoff drops. The third alternative is now "at or below the cutoff".

The comparison of the chain with the exact operator, which follows this check, holds as before.

Documentation

• parameter_section.rst: the new default, and a description of lambda_cutoff, which had none.
• output_format.rst: the example of parameters.dat.
• NEWS.md.

English and Japanese.

Tests

• New: rank_deficient_bond, an end-to-end test with the dimers and with the smallest input that showed the second case: one Ising bond over two lattice spacings, six seeds. With a build of develop the dimers and five of the six seeds fail.
• test/input.cpp: the default.
• ctest: 78/78 (Apple clang, no MPI, Release).

Not tested: an MPI build.

🤖 Generated with Claude Code

[yomichi]

Thanks. The analysis is convincing, and I agree with 1e-7. rank_deficient_bond fails with lambda_cutoff = 1e-12 injected (the dimers, and 3 of the 6 seeds of the chain), and passes for 40 of 40 seeds with the new default. The tighter check in longrange_gate.cpp and unit_d1.cpp is an improvement: with the former cutoff, the forbidden band [1e-8 * lmax, 1e-9] was empty, since lmax >= D^(-1/4).

1. Normalization of lambda in the docs.
parameter_section.rst says lambda is "the square root of a Schmidt value normalized to the largest scale of the bond", and the description and NEWS speak of "the largest Schmidt value". The code (simple_update.cpp, around line 225) divides by the 2-norm: lambda_i = sqrt(s_i / ||s||_2). The conclusion does not change, but could the English docs say "divided by the 2-norm of the Schmidt values of the bond"? The Japanese text ("規格化された Schmidt 係数の平方根") is fine as it is.

2. (FYI) The table of cutoffs depends on the build.
With my build (Apple clang, Release, no MPI), 1e-10 gives a wrong value for the dimers, while the table says right. 1e-8, 3e-8, 1e-6 and 1e-5 were all right. This only confirms that the result below the threshold is a matter of chance; no change is needed.

Not tested: an MPI build.

[aoymt]

Thank you for checking this, and for the 40 seeds.

1. Fixed in b8f8868: parameter_section.rst and NEWS.md now say that lambda is the square root of a Schmidt value divided by the 2-norm of the Schmidt values of the bond.
2. Noted. I leave the table as it is; your 1e-10 and my 1e-10 disagreeing is the point it makes.

On #123: I reproduced it, and the fix is in fix/gauge-fix-rank-deficient (PR to follow). As you note there, the dimers with gauge_fix = true run under 1e-12 and crash under 1e-7, so that fix should go in before this PR or together with it.

[yomichi]

Thanks, the normalization now matches simple_update.cpp, and leaving the table as it is is fine. Approving; to be merged after #124, as you note, so that gauge_fix = true on a rank-deficient bond keeps running.