I think the embedding idea in specific is quite good and there are several potential approaches! One is, if you let xᵢ be the song vector for user i (whose kth entry is 1 if user i has song k in their library and 0 otherwise). Then you can compute the "overlap" by something like xᵢ ⋅ xⱼ (with some normalization of course!) where ⋅ is the inner product. [0]
A simple approximation of this inner product would be to generate a random (potentially sparse!) matrix S whose nonzero entries are i.i.d. Gaussian, for example, and whose number of rows is much smaller than the number of columns [1], then you can instead store and compute
(Sxᵢ) ⋅ (Sxⱼ)
which gives you an approximate overlap, whose storage and computation requirements are much smaller for each user (since Sxᵢ is much smaller in number of entries than xᵢ).
-----
[0] Of course, there are many other similar methods! This is a particularly simple, but often fairly effective one.
[1] More specifically, it goes like O(log(n)/ε²) where ε is the error you wish to achieve. Often, a fairly large choice of ε actually will suffice. See https://en.wikipedia.org/wiki/Johnson–Lindenstrauss_lemma
A simple approximation of this inner product would be to generate a random (potentially sparse!) matrix S whose nonzero entries are i.i.d. Gaussian, for example, and whose number of rows is much smaller than the number of columns [1], then you can instead store and compute
which gives you an approximate overlap, whose storage and computation requirements are much smaller for each user (since Sxᵢ is much smaller in number of entries than xᵢ).-----
[0] Of course, there are many other similar methods! This is a particularly simple, but often fairly effective one.
[1] More specifically, it goes like O(log(n)/ε²) where ε is the error you wish to achieve. Often, a fairly large choice of ε actually will suffice. See https://en.wikipedia.org/wiki/Johnson–Lindenstrauss_lemma