Potential INF in LNCC loss

[mathpluscode]

Describe the bug

The current implementation of LNCC loss https://github.com/Project-MONAI/MONAI/blob/master/monai/losses/image_dissimilarity.py comes from DeepReg.

However, recently in DeepReg, we encountered INF loss, due to a negative variance at certain voxels. The variance should be zero, but somehow it became around -1e-5 and sometimes it eventually caused INF.

The fix in DeepReg is

• First normalize the kernel so that the sum of the kernel weights equals 1. For some reason, this makes the variance be around machine error (-1e-11) and by adding an EPS, we will not encounter INF anymore.
• We've increased the EPS from 1e-7 to 1e-5, as under mixed-precision, we think 1e-5 is safer.
• Further, we clip the variance to zero, so that the variance must be >=0 in following calculation.

Unfortunately, I could not provide some simple tests to reproduce the inf error (due to the difficult debugging in TensorFlow).

More details are in

• Nan/inf loss encountered DeepRegNet/DeepReg#690
• 690 nan inf loss DeepRegNet/DeepReg#719

[YipengHu]

FYI @tvercaut

[tvercaut]

I haven't looked at the code but a scenario like the one you describe arises frequently if you are using a naive algorithm to compute the variance, e.g. something like:
∑x² - (∑x)²

While correct from a mathematical standpoint, a naive computation of the variance is numerically unstable. You can either use a 2-pass algortithm or or slightly more involved online algorith:
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

I would assume that pytorch uses a stable algoruthm for this by default
https://pytorch.org/docs/stable/generated/torch.var_mean.html
but if not, this would need to be raised upstream.

[mathpluscode]

I haven't looked at the code but a scenario like the one you describe arises frequently if you are using a naive algorithm to compute the variance, e.g. something like:
∑x² - (∑x)²

While correct from a mathematical standpoint, a naive computation of the variance is numerically unstable. You can either use a 2-pass algortithm or or slightly more involved online algorith:
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

I would assume that pytorch uses a stable algoruthm for this by default
https://pytorch.org/docs/stable/generated/torch.var_mean.html
but if not, this would need to be raised upstream.

Thanks, @tvercaut

• In LNCC case, we are calculating weighted variance for each voxel via Convolution. Thus, the online algo mentioned in wiki might not be an option :( I may need to dig further if there are more detailed explanations for the instability.

• Regarding the bias/unbiased, this will not solve the stability issue as we are adjusting the variance by a constant without changing the sign? But still good to correct the variance.

[tvercaut]

The online algorithm (Welford) might be an overkill anyway here, see e.g.:
numpy/numpy#6231

A basic two-step approach is typically enough:
μ=(∑x)/N
σ²=(∑x-μ)²/(N-1)

[mathpluscode]

The online algorithm (Welford) might be an overkill anyway here, see e.g.:
numpy/numpy#6231

A basic two-step approach is typically enough:
μ=(∑x)/N
σ²=(∑x-μ)²/(N-1)

Cool, then we will do this. Gonna try to find some reference to support the claim that

∑x² - (∑x)² is not safe

[tvercaut]

The wikipedia link I referred to points to this article:

Chan, Tony F., Gene H. Golub, and Randall J. LeVeque. "Algorithms for computing the sample variance: Analysis and recommendations." The American Statistician 37.3 (1983): 242-247.

Can't think of a more definitive support for this :)

[wyli]

came across a similar issue a while ago in the batch norm/tf.nn.moments, may need to make sure it's using a high precision for the intermediate values: https://github.com/tensorflow/tensorflow/blob/85c8b2a817f95a3e979ecd1ed95bff1dc1335cff/tensorflow/python/ops/nn_impl.py#L1330-L1344

[mathpluscode]

came across a similar issue a while ago in the batch norm/tf.nn.moments, may need to make sure it's using a high precision for the intermediate values: https://github.com/tensorflow/tensorflow/blob/85c8b2a817f95a3e979ecd1ed95bff1dc1335cff/tensorflow/python/ops/nn_impl.py#L1330-L1344

Currently, DeepReg does not support mixed precision yet. So the problems occurred under float32. But indeed, we should enforce float32 for this variance/mean.

[wyli]

fixed by #1915, please reopen if needed