How fast can we compute the set inclusion poset of a set family? Given a set family $\mathcal{F}$ of subsets of a universe $U$.
Let $S_1,S_2 \in \mathcal F$ and we want to answer is $S_1 \subseteq S_2$.

I am looking for a data-structure that will allow me to quickly answer this.
My application is from graph theory where I want to see if deleting a vertex and its neighbourhood leaves any isolated vertices, and for each vertex list all isolated vertices it leaves.

I want to create the complete poset or eventually a $|\mathcal{F}|^2$ table storing true false telling exactly which sets are subset of eachother.

Let $m = \sum_{S\in \mathcal{F}} |S|$, $u = |U|$ and $n = |\mathcal{F}|$, assume $u,n \leq m$

We can generate the $n \times u$ containment matrix (the bipartite graph) in $O(un)$ time and then can create the table of all $n^2$ comparisons in $O(nm)$ time by for each set $S \in \mathcal{F}$, loop through all elements of all other sets and mark the set as not a subset of $S$ if they the element is not in $S$.
In total $O(nm)$ time.

Can we do anything faster?
In particular, is $O((n+u)^2)$ time possible or not?

I found some related articles:

A Simple Sub-Quadratic Algorithm for Computing the Subset Partial Order (1995)
which give an $O(m^2 / log(m))$ algorithm.

The Subset Partial Order: Computing and Combinatorics slightly improves the above but also claim that the above paper solves the problem in $O(md)$ time where $d$ is the max number of sets sharing a common element, but I could not understand this result.

In the article Between $O(nm)$ and $O(n^{\alpha})$ the authors show how to in a graph find the connected components after deleting the closed neighbourhood of a vertex by using matrix multiplication. This can be used to compute the set inclusion poset by finding all the components which are singletons with a runtime of $O((n+u)^{2.79})$.

Also this forum discussion is related: 
What is the fastest way to check for set inclusion?
which implies a lower bound of $O(n^{2-\epsilon})$.
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.