How fast can we compute the set inclusion poset of a set family?
If randomness is in bounds, one rough idea would be to generate a bunch of "random monotonic signature" functions and use them to approximate the subset relation (a la Bloom filters).  Unfortunately, I don't know how to make this into a practical algorithm, but here are some estimates that don't immediately prove the idea impossible.  This is very far from a useful solution, but I'll write it out in case it helps.

Assume for simplicity that the sets are all nearly the same size, say $|S| = s \pm O(1)$, and that $s = o(u)$.  We can assume $1 \ll s$, otherwise we're done.  Define
$$\begin{aligned}
q &amp;= [s/2] \\
p &amp;= \left[\frac{u \choose q}{s \choose q}\right]
\end{aligned}$$
Note that $p \gg 1$.

Here is the wildly impractical part.  Randomly choose $p$ subsets $A_1, \ldots, A_p \subset U$ with replacement, each of size $q$, and define a function $f : 2^U \to \{0,1\}$ by $f(S) = 1$ iff $A_i \subset S$ for some $i$.  With $S$ fixed and $A_i,f$ varying randomly, we have
$$\begin{aligned}
\Pr(f(S) = 0) &amp;= \Pr(\forall i. A_i \not\subset S) \\
              &amp;= \Pr(A_1 \not\subset S)^p \\
              &amp;= \left(1 - {s \choose q}/{u \choose q}\right)^p \\
              &amp;= e^{-\Theta(1)}
\end{aligned} $$
Since $f(S)$ is monotonic, $S \subset T$ implies $f(S) \le f(T)$.  If $T \not\subset S$, fix some $t \in T-S$.  The probability that $f$ detects $T \not\subset S$ is
$$\begin{aligned}
\Pr(f(S) = 0 &lt; 1 = f(T)) &amp;= \Pr(f(S) = 0) \Pr(f(T) = 1 | f(S) = 0) \\
&amp;= e^{-\Theta(1)} \Pr(\exists i. A_i \subset T, A_i \cap T-S \ne 0 | f(S) = 0) \\
&amp;= e^{-\Theta(1)} \Pr(\exists i. t \in A_i \subset T | f(S) = 0) \\
&amp;\le e^{-\Theta(1)} \Pr(\exists i. t \in A_i \subset T) \\
&amp;\approx e^{-\Theta(1)} p \Pr(t \in A_1 \subset T) \\
&amp;\le e^{-\Theta(1)} p {s \choose q-1} / {u \choose q} \\
&amp;\approx e^{-\Theta(1)} p \frac{q}{s-q} {s \choose q} / {u \choose q} \\
&amp;= e^{-\Theta(1)}
\end{aligned}$$
Some of those steps are pretty tenuous, but I don't have time to improve them tonight.  In any case, if they all hold, then at least it's not clearly impossible to randomly generate signature functions that have reasonable likelihood of distinguishing subsets from nonsubsets.  A logarithmic number of such functions would then distinguish all pairs correctly.  If generating a signature function $f$ and computing $f(S)$ could be reduced to $\tilde{O}(n+u)$ time, the result would be an overall $\tilde{O}(n^2+u^2)$ algorithm.

Even if the above calculations are correct, I have no idea how to generate monotonic signature functions with the desired features quickly.  It's also likely that this technique doesn't extend to significantly different set sizes.