Why this choice?
In my tests for #11271 I found that the computation time slowly increase with the memory size. It's the first few plots of this comment: #11271 (comment).
I'm mostly curious about this choice. I don't think this should prevent this PR from being merged.

A fixed 10Mib maxmem can lead to very small chunks when the number of features is very large (extreme case I agree). In that case the drop of performance can be quite severe.

Another possibility could be to have a fixed maxmem with a min chunk_size. I think it's kind of equivalent, but in this case you can have a memory increase bigger than expected whereas in the previous case the memory increase is bounded. That's why I went for the first one.

Large chunks would also decrease cache efficiency. Although I'm not a big fan of having an unbounded amount of additional memory, 10% of the memory already used is probably fair.

When I look at my old benchmarks, I think that there might be something going on under the hood. Maybe the optimal block size is actually constant. Moreover, if the number of features is large, it can also be chunked. That might deserve to be explored as well.
I think we can add to the todo-list to benchmark the influence of the memory size and chunking to get a proper understanding of the performance.

If the input data is memmapped, calculating this as a function of n_features may not be fair.

I think my comment about memmapping was incorrect in any case.

This logic seems to be more complicated than necessary, but okay.

I am fine to merge with this heuristic.

In the future, I am wondering if we could expose a chunk_size parameter which by default do such heuristic and could be fine-tuned by the user. @jeremiedbb mentioned to me that it might not be that easy since it will impact pairwise_distance where we would need to give also this parameter, etc. Just a thought.

Can we not use utils.gen_even_slices here?

gen_batches since we want a specific batch size, but it's cleaner indeed.

Nit: X.nnz maybe

"with float32 features"

maybe more relevant: "and therefore several estimators with metric='euclidean'"

I think correct syntax is just :user:`Celelibi`

"_fast" is a bit obscure

it's a residue of my previous PR implementing fast and exact method. I removed it.

If the input data is memmapped, calculating this as a function of n_features may not be fair.

add "in 8-byte blocks"... in fact where do you account for this?

maxmem is the max allowed number of 8-byte blocks for extra memory usage.
10Mib is 10 * 2**20 bits = 10 * 2**17 8-byte blocks.

I think my comment about memmapping was incorrect in any case.

should we raise a warning that these cannot be used??

Why not but I'm afraid it could generate a lot of warnings. For instance when pairwise_distances is called in a loop.

One thing we should certainly do is adding an option to row_norms to use a float64 accumulator. That's on my todo list :)

why is this out of the if? It's only applicable if the if holds.

(Curiously we currently appear to duplicate this logic in pairwise_distances_chunked but not in pairwise_distances)

you mean zeroing the diagonal ? It applies using both methods.

It's necessary to also do it in pairwise_distances_chunked, because we lose the diagonal info when passing chunks to pairwise_distances.

Yes, but I don't think we zero the diagonal in pairwise_distances at all for other metrics

I was looking at euclidean_distances, not pairwise_distances... That's strange indeed. Even for metric=euclidean it's not enforced when n_jobs > 1. Let me put that also on my todo list :)

This logic seems to be more complicated than necessary, but okay.

When a grade school student ask, "Why are we learning the quadratic equation?", we have an answer.

At first glance, it is not obvious that None is used by _euclidean_distances_upcast.

What do you mean ? XX and YY are used by _euclidean_distances_upcast no matter their value

Setting XX or YY to None tells _euclidean_distances_upcast to call row_norms(*_chunk, squared=True)[:, np.newaxis].

Without reading _euclidean_distances_upcast, it is difficult to tell why XX or YY are set to None.

I think that we would need a small comment to clarify this part.

fast method?

I made the comment clearer

I put the comment here but I think that it should be written above.

I would either add a note (using sphinx syntax) or within the parameter description that X_norm_squared and Y_norm_squared will be ignored if they are passed as float32 since they will lead to inaccurate distances.

I think that we would need a small comment to clarify this part.

I am fine to merge with this heuristic.

In the future, I am wondering if we could expose a chunk_size parameter which by default do such heuristic and could be fine-tuned by the user. @jeremiedbb mentioned to me that it might not be that easy since it will impact pairwise_distance where we would need to give also this parameter, etc. Just a thought.

a small comments regarding rtol=1e-6

Maybe we could have a better name for s

Add copy=False, in the case when d = distances[y_slice, x_slice].T it is already float32, and astype make a copy by default.

this function only computes euclidean distances on float64 and downcast it to float32 at the end so actually there will always be a copy. Actually I find it more clear to not add the copy=False parameter to emphasize that

It must be correct since the corresponding tests pass, but after thinking about it for 20s I'm still not confident that when we compute distances[x_slice, y_slice], distances[y_slice, x_slice] is already computed, is not null and can be copied from.

How about just adding a continue here, then add a second loop below,

if X is Y:
   # the result is symmetric, copy distances from the lower triangular matrix
	for i, x_slice in enumerate(x_batches):
    	  for j, y_slice in enumerate(y_batches):
              if j < i:
              	  distances[x_slice, y_slice] = distances[x_slice, y_slice].T

it shouldn't matter for run time too much time I think, but it might make things a bit more explicit?

still not confident that when we compute distances[x_slice, y_slice], distances[y_slice, x_slice] is already computed

It is since we compute the chunks with i (rows) as the outer loop. It implies that all chunks of the upper right triangular part of the distance matrix are computed before their symmetric chunk in the lower left triangular part.

Can we use MiB (or MB) as all other cache sizes in scikit-learn? Mib unit can lead to confusion when comparing to the CPU L3 cache size or to working_memory.

This makes the assumption that sizeof(dtype) == 8 without making it appear in the computation, I think?

Maybe we could add X.dtype.itemsize / 8 even if that is equal to 1, to make it easier to follow what is going on.

Also not that it matters too much, but

1. for dense isn't that its 5% of (X, Y, distance) as those are float32. i.e. X.nbytes == x_density * n_samples_X * 4 / 8?
2. for sparse it's a roughly 10% if it is CSR and we account for X.data, X.indptr in 32 bit, but again it's far from obvious for a casual reader.

Actually, maybe,

def get_array_nbytes(X):
   if issparse(X):
      if hasattr(X, 'indices'):
         # CSR or CSC
     	 return X.data.nbytes  + X.data.nbytes + X.data.indptr
      else:
         # some other sparse format, assume 8 bytes to index
         # one non zero element (e.g. COO)
         return X.data.nbytes + X.nnz*8
   else:
       return X.nbytes

maxmem = (get_array_size(X) + get_array_size(Y) + get_array_size(distances)) / 10

This makes the assumption that sizeof(dtype) == 8 without making it appear in the computation, I think?

This function is only used for float32 X and Y. I've explained it in its docstring.

for sparse it's a roughly 10% if it is CSR and we account for X.data, X.indptr in 32 bit, but again it's far from obvious for a casual reader.

When you change the type of the csr matrix, only a copy of data is made. indices and indptr are still int, so I think this formula is correct

When you change the type of the csr matrix, only a copy of data is made.

Are you sure?

>>> import scipy.sparse
>>> import numpy as nop
>>> import numpy as np
>>> X = scipy.sparse.csr_matrix(np.ones((10, 10)))
>>> X.dtype
dtype('float64')
>>> Y = X.astype(np.float32)
>>> np.shares_memory(Y.data, Y.data)
True
>>> np.shares_memory(Y.data, X.data)
False
>>> np.shares_memory(Y.indices, X.indices)
False

That's not critical, my point here is that it might be beter (not necessarily in this PR) to use ndarray.nbytes than trying to re-compute that from scratch. If we get this calculation wrong we won't likely know, since no test will fail and we will only reason in wrong sizes of memory cache.