Find close pairs in very high dimensional space with sparse vectors I have $N$ (~a million) feature vectors. There are $M$ (~a million) binary features, but in each vector only $K$ (~a thousand) of them would be $1$, the rest are $0$. I'm looking for the pairs of vectors that have at least $L$ (~a hundred) features in common ($1$ in both). The number of such pairs is of a similar magnitude to $N$ (~a million).
I think this could be approached as looking for close point pairs in a very high-dimensional space. The distance function could be such that it's based on how many features the two vectors have in common. But it would probably be useful with a more conventional distance metric (such as Euclidean) as well.
Which well-known algorithms would be useful for approaching this problem? Anything that is quadratic in $N$ or $M$ will not be practical.

An example real-world formulation of the problem is to consider $N$ people  moving between a number of locations. If two people were at the same location at the same time, we say they met each other. (The number of location-time combinations with at least 1 person present is $M$.) We are looking for friends: people who met at least $L$ times.
 A: 
I'm looking for the pairs of vectors that have at least $L$ features in common.

This is just an inner product of binary feature vectors. When the inner product is greater than $L-1$, the pair will have at least $L$ elements in common. This should be a relatively fast computation -- at least, faster than euclidean distance, which would be wasteful and slow for this data. Because you stipulate that you're  looking for pairs, this will inherently mean you have to do $\binom{N}{2}$ computations to compare every vector. 
Finding points that are close together is indeed a clustering problem. But the first step of the clustering algorithms that I'm familiar with is computing pairwise distances or similarities. I'm sure someone has developed more efficient alternatives. A point about terminology: having at least $L$ common neighbors is phrased as a similarity, not a distance! Inner products are, in this case, unnormalized cosine similarities.
You can make this more tractable by only performing the inner product computation when the sum of the feature vector (which is in this case the same as the norm) for an observation is greater than $L-1$, since it's impossible for that binary feature vector to have an inner product with another binary feature vector which will satisfy my criterion when this sum is less than $L$. Obviously, computing these sums is only $O(N)$ complexity, so i's a cheap way to drive down the magnitude of the inner product step.
But the classic way to reduce the scope of this problem is to do additional pre-filtering. Are you especially interested in when one, somewhat uncommon feature takes the value 1? If so, only perform the computation for those feature vectors.
Or perhaps you could benefit from re-framing your problem. For example, sampling is known to have nice properties; inferential statistics develops on this idea to quite some depth. So perhaps it's unfeasible to analyze the entire data set, but it's perfectly feasible to examine a small sample. I don't know what question you're trying to answer, but if you carefully design your experiment, you may get away with only looking at a few thousand observations, with more than enough data left over for validation purposes.
After some additional thought, I have a strong hunch that the data you're working with is some kind of graph $G$. It's very plausible that $G$ is composed of several connected components, in which case you can decompose $G$ into a set of graphs, with the happy side-effect of reducing the dimensionality of the data. Even if the graph is only two connected components of roughly the same size, that means your $O(N^2)$ pairwise comparisons has roughly $\frac{1}{4}$ the total cost! 
If the graph is symmetric, the following observations may be helpful:


*

*Define the Laplacian of your graph as $P=D-A$, where $D$ is a diagonal matrix of degree (the sum of each feature vector) and $A$ is the adjacency matrix (the stacking of feature vectors into a matrix).

*The number times $0$ appears as an eigenvalue of $P$ is the number of connected components of $G$. Decomposing the graph into its connected components  and working solely with those components will have the side-effect of reducing the dimension of your data; computing your quantity of interest will be easier. But computing the eigendecomposition will be expensive for a million vertices...

*(After a full permutation) $P$ is a block diagonal matrix of the Laplacians of the  connected components of $G$.

*$P$ is positive semidefinite. This is almost certainly useful somehow.

*The algebraic connectivity of $G$ is the value of the second-smallest eigenvalue of $P$. This tells you how well-connected $G$ is. Perhaps that will answer some of the questions you are interested in re: the vectors that have features in common. Spectral graph theory develops this idea in some more detail.



"Is this an SNA problem?" I'm not sure. In one application the features describe behavior and we're looking to connect people with similar behaviors. Does that make this an SNA problem?

If you have a bipartite graph connecting people to behaviors, you can think of this as an affiliation network $B$, with people as rows and behaviors as columns. If you want to connect people to people via the behaviors they have in common, you can compute $BB^T=A$. $A_{ij}$ is the number of behaviors the people have in common. Obviously, the set of vertices where $A_{ij}\ge L$ answers your question.
A: It looks like the approach you're looking for is a combination of minhash signatures and Locality Sensitive Hashing (LSH); the (freely available) pdf of Mining Massive Datasets describes this approach (and other similarity measures) in some detail in Chapter 3, but briefly:
A minhash signature is a condensed representation of your original matrix that is constructed by applying some number n of hash functions to features, thereby reducing the number of features per observation. This reduces the size of your data, however you'll probably notice that this still leaves you with a $O(N^2)$ problem.
To address this, MMDS advises that if all you want to find is pairs above a certain threshold of similarity (which would seem to apply in your case), then you can focus only on those pairs that are most likely to be similar - this approach is called Locality Sensitive Hashing, and in section 3.4 they walk through an example of how to combine the minhash signature approach with LSH.
In addition to the text, there are also lectures available on the Coursera course of the same name.
A: On looking for people meeting in space-time blocks:
split space into $Nspace$ blocks (city blocks, square km, whatever),
and time into $Ntime$ blocks.
There's a good chance that if people meet, they'll meet within the same block.
So run NN within each block.
Runtimes and error rates will of course depend on block sizes and shapes
(also on what you can parallelize / MapReduce),
but you have parameters to play with -- engineering, not wide-open $O( N^2 )$.
See also:
nearest-neighbors-search-for-very-high-dimensional-data on datascience.stackexchange   
pairwise.py:

uses the Python Gensim library and heapq from the standard library to make
  massively fast and scalable pairwise comparisons between an aribtrarily large number of
  documents using TF-IDF and cosine distance.

A: Inverted dictionnary! Represent a point $x$ as $feat_1:value_1, feat_{101}:value_{101}$, the keys corresponding to non zero values (i.e. the features holding true). The average size of storage of an element will be $K$. Indeed, I only need $K$ strings to store the features and $K$ floats to hold the values.
For each feature, build a dictionary holding the indexes sharing this feature. Hopefully, this number will not be too big (if you have a feature which is shared by all the indexes, this approach is ruined, you can stop reading here).
This dictionary looks like : $feat_1 : \{1,101,202\}, feat_2 : \{7,202\},feat_3 : \{202\}...feat_M:\{3,45,6\}$. If I want to gain speed and save space, I can even drop the features that are only found with one element (here:$feat_3$) as they will not produce close pairs. This dictionary is built in $O(NK)$ operations.
Now, when you want to evaluate the distance of an element $x$ to the others, generate (with the dictionary) the list of the indexes sharing at least one feature with $x$. You know that all the other elements are far from $x$ (they don't even share one feature!). If the average number of "elements per feature" is low (call it $P$), you need not to be in $O(N^2)$ any more.
Now there is another big improvement if $x$ and $y$ are represented as dictionaries as well, since $d(x,y)$ or $<x,y>$ can be evaluated iterating over the keys of $x$ and $y$, in $O(K)$ operations.
Your final complexity is $O(NPK)$ instead of the naive $O(MN^2)$ initial approach.
I applied this method to implement a KNN over large text set (train : 2 000 000 lines, test 35 000 lines, number of features : 10 000, average number of features per element : 20), which ran in about an hour...
A: I've found a reference that you might find helpful, and I believe it's asymptotically more efficient than every other solution presented so far. If I understand correctly, you can construct a $k$-nearest neighbor (KNN) graph in $O(LN\log(N))$ time.
L. Erotz, M. Steinbach, and V. Kumar. "A new shared nearest neighbor clustering algorithm and its applications." Proceedings of the 1st Workshop on Clustering High Dimensional Data and its Applications, 2002.
A: A crazy, but likely to work approach might be to go to frequency domain.  There is a crazy/sick fast fft called "sparse FFT" where you specify the number of modes you care about (your count of 100 features) and then you work in convolutions, and look for row-max greater than a threshold (look for bits in upper registers of your numbers).  It is going to be $O(k \cdot \log{n} )$ where $ k << n $. 
Given that your k is 100 and your n is 1e6, this should give you ~1e4x speed up compared with classic FFT.  
If you need another 20x in speed, and you are a risk taker, then instead of convoluting all rows against the domain and looking for the peak, you could bootstrap a subset of rows.  
You might also prefilter columns by removing columns whose sums are below 50, or some other threshold that is on the order of half the number of rows you are looking to match.  At the very least you should remove columns of all zeros and all 1's as non-informative.  Same with rows that are entirely empty or empty enough, or rows that are so full that they are irrelevant.
To-do:  I should put an example here using synthetic data, and compare some of the methods.
A: I just came across a paper that is directly relevant.

Randomized algorithms and NLP: using locality sensitive hash function for high speed noun clustering (Ravichandran et al, 2005)

It is actually implemented in https://github.com/soundcloud/cosine-lsh-join-spark which is where I found it.
It is based on locality sensitive hashing (already mentioned in other answers). After it has reduced the feature vectors to a low-dimensional space it uses a fast Hamming distance join to find the nearest neighbors.
