The digit similarity matrix

Flatten the 1,797 real handwritten digits to (1797, 64) and compute the similarity between every pair with a single einsum:

np.einsum('id,jd->ij', D, D)    # (1797, 1797)

d is the pixel index and it is contracted; i and j both survive, one per digit, which is why the result is square. That is 3.2 million dot products from one line — and it is exactly the operation a recommender runs to rank items for a user.