[MRG] MNT Fixes for PCA with n_components='mle'

[NicolasHug]

Fixes #16730
Closes #16546

• infer_dimension() now returns a value in [1, n_features - 1]. rank = 0 isn't possible anymore (it's a bug fix, as this would transform into an empty array).
• rank as passed to _assess_dimension() is actually the rank, not an index as previously. This fixes the infamous off-by-one error. Note also that the method isn't defined if we pass it rank == n_features (division by zero)
• a few fixes in the formula

[NicolasHug]

ping @adrinjalali @jnothman @agramfort @lschwetlick

[lschwetlick]

This looks good to me, thanks for fixing!

[NicolasHug]

Thanks for checking @lschwetlick . If you have time / feel like it, it'd be great if you could check with your notebook if the results make sense now? (but no worries otherwise)

[lschwetlick]

Well... not sure if this is maybe quite edge case-y but

b = np.ones((9, 6))
print("rank=", np.linalg.matrix_rank(b))
u2, s2, vh2 = np.linalg.svd(b, full_matrices=True)
ll = np.empty_like(s2)
ll[0] = -np.inf
for r in range(1,s2.shape[0]):
    ll[r] = (_assess_dimension(np.asarray(s2), r, 9))
print(ll)

-> [-inf nan -inf -inf -inf -inf]

but then funnily enough it gives the right number of dimensions (1) because apparently argmax treats nan>-inf. Its giving nan because if all subsequent values in the spectrum are 0 then in line 76 v is also 0...

[lschwetlick]

I'm also still slightly confused to be honest: we dont test full rank because the method isn't defined for testing full rank. But what if the number of dimensions is full rank?

[NicolasHug]

Thanks @lschwetlick , I fixed that and added a test (though indeed it's an unlikely case I think)

But what if the number of dimensions is full rank?

I guess users use mle when they want some dimensionality reduction, but just don't know by how much. So that's fine not to consider rank==n_features (we can't anyway). There's always the solution to directly provide n_components=n_features.

[lschwetlick]

Alright cool, other than that my notebook and I are happy :)

[agramfort]

thanks a lot @NicolasHug and @lschwetlick for the team work !