Finding a large clique of an impractically large generalised Kneser graph

Question

My original problem statement is simple. Find a maximal clique of $k$-length subsets of of a set of $n$ items, clique members sharing at most $s$ items with any other.

This is the maximal clique of a generalised Kneser graph $K(n, k, s)$.

A simple and straightforward example would be:

With[{n = 6, k = 3, s = 1},
 FindClique[
  RelationGraph[Length[Intersection[#1, #2]] <= s &, 
   Subsets[Range[n], {k}]]]]

(* {{{1, 2, 3}, {1, 4, 5}, {2, 4, 6}, {3, 5, 6}}} *)

Clearly the maximal clique in this case has four items.

Problems arise when number of subsets rises to millions, and the number of edges between them would be in an order of trillions. This is clearly not practical to handle explicitly. Also, without explicit reliance on the specific properties of the graph running any sort of normal maximal clique algorithm on the graph seems futile - time complexity of the process is just too high.

My first attempt at this is a greedy sieve, taking the first item from the list of subsets, removing all entries (vertices) not fulfilling the overlap criteria with it from the rest of the list, appending the first element to the solution clique and repeating the process until the subset list is empty. This way edges are kept implicit, and even with densely connected graph of tens of millions of vertices, the storage demands can be kept at bay.

It has a problem, though. Like many greedy algorithms, the result depends on the chosen order of operations, and even with this simple example, can find two different clique sizes:

ParallelTable[
  With[{n = 6, k = 3, s = 1}, 
   NestWhile[
      Apply@Function[{list, items}, {Append[list, First@items], 
         Select[Rest@items, 
          Length[Intersection[First@items, #]] <= s &]}],
      {{}, RandomSample@Subsets[Range[n], {k}]},
      Length@#[[2]] > 0 &] // First // Length], 10000] // Tally

(* {{4, 9007}, {2, 993}} *)

So, although it does find a clique, it is not necessary a maximal one. Being able to find large cliques approaching the maximal clique size efficiently (with more intelligent sampling?) might make sense.

I also attempted looking at this problem using Boolean SatisfiabilityInstances, but the problem (time complexity) explodes even faster if the solution is sought that way.

I'm interested of any solution which would find the maximal clique when, say, $n = 40, k = 7$ and $s = 1$ or $2$. This far I've used the random sampling algorithm above, rewritten it to very optimised form (about 3500 times faster than the Mathematica version), and found out that the algorithm provides a clique of average size of about 15, 19 being the largest seen, but having a suspicion that at least a clique of size 20 would be probably found after several days of extra computing time spent.

EDIT: Clearly heuristics can help. Instead of choosing purely random vertices on the sieve step of my algorithm results are much better if ones with more overlap with earlier chosen vertices are used, providing cliques of 22 vertices with just 1000 runs, when unbiased random order got only to 19 with roughly ten million runs. Largest clique found this far with this method has 24 vertices, of which each has set intersection size of one pairwise with other vertices.