ENH dwidenoise: Optimal shrinkage

[Lestropie]

Not familiar with what may be considered to be best practise here, but I at least have a basic sense of the concept.

Depending on where the upper threshold of the MP distribution is determined to be relative to the component eigenvalues, inclusion of components in the output DWI series could be fractional.

[Lestropie]

Animation below shows raw data, then denoising using eigenspectrum truncation, then optimal shrinkage.

@dchristiaens @jdtournier Would appreciate some input on the optimal shrinkage implementation:

1. To make sure that my implementation is correct (some of the variable handling was already different to that in published manuscript)
2. Given that the noise level threshold is determined based on choosing a specific component index, we want to make sure that the truncation / alignment of the optimal shrinkage function relative to that selection is correct; ie. check for off-by-one errors.

[Lestropie]

The equation currently in use here is based on:
"Optimal Shrinkage of Singular Values"
https://ieeexplore.ieee.org/document/7820186/?arnumber=7820186

I have only just realised that this was not the reference in Lucilio's paper:
"Optimal Shrinkage of Eigenvalues in the Spiked Covariance Model"
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6152949/

[Lestropie]

I believe that the implementation as at 471be56 in #3029 is consistent with the tensor MPPCA implementation:

• Raw eigenvalues are divided by the size of the max matrix dimension and sigma^2 (as opposed to lambda_plus)
• The same non-linear transform function is applied

The subsequent rescaling differs due to other implementation details:

• They immediately produce rescaled versions of the eigenvalues by multiplying by the size of the max matrix dimension and sigma*2.
• Our 3.0.4 implementation was based on producing a weights vector, originally with 1.0 for preserved components and 0.0 for truncated components. In my updated code I produce weights con a continuous scale, as the ratio of the input to and output from the non-linear transform function.

So I think that makes me confident enough that the current implementation is correct and adequate to close off this Issue. Any sanity checking from someone with experience in the mathematics in the course of reviewing #3029 once it's out of draft form would nevertheless be appreciated.