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.