The authors establish a large-scale regularity theory for the steady Stokes problem in presence of a random suspension of small and rigid particles \( I_{n}\): \(-\Delta u_{g}+\nabla P_{g}=\operatorname{div}(g)\), \(\operatorname{div}(u_{g})=0\), in \(\mathbb{R} ^{d}\setminus \mathcal{I}\), \(D(u_{g})=0\) in \(\mathcal{I}\), \(\int_{\partial I_{n}}(g+\sigma (u_{g},P_{g}))\nu =0\), \(\int_{\partial I_{n}}\Theta (x-x_{n})\cdot \sigma (g+\sigma (u_{g},P_{g}))\nu =0\) \(\forall n\), \(\forall \Theta \in \mathbb{M}^{\mathrm{skew}}\), the subset of skew-symmetric matrices. Here \( \sigma (u,P)=2D(u)-PId=(\partial _{j}u_{i}+\partial _{i}u_{j})-PId\). The authors also prove moment bounds for the associated correctors and optimal estimates on the homogenization error. They introduce a probability space \( (\Omega ;\mathbb{P})\), a random point process \(\mathcal{P}=\{x_{n}\}_{n}\) on \(\mathbb{R}^{d}\), \(d\geq 2\), and a collection of random shapes \( \{I_{n}^{o}\}_{n}\), where \(I_{n}^{o}\) is a connected random Borel subset of the unit ball \(B\) which is centered at 0, that is which satisfies \( \int_{I_{n}^{o}}xdx=0\). They define the corresponding inclusions \( I_{n}=x_{n}+I_{n}^{o}\) centered at the points of \(\mathcal{P}\). The authors assume that the random set \(\mathcal{I}=\cup _{n}I_{n}\) is stationary and ergodic, that the random shapes \(\{I_{n}^{o}\}_{n}\) satisfy interior and exterior ball conditions with radius \(\delta >0\), a fixed deterministic constant, almost surely, and that \((I_{n}^{+}+\delta B)\cap (I_{m}^{+}+\delta B)=\varnothing \), almost surely for all \(n\neq m\), where \( I_{n}^{+}\) denotes the convex hull of \(I_{n}\). They also assume that there exists a non-increasing weight function \(\pi :\mathbb{R}^{+}\rightarrow \mathbb{R}^{+}\) with superalgebraic decay: \(\pi (l)\leq C_{p}\left\langle l\right\rangle ^{p}\) for all \(p<1\), such that the random set \(\mathcal{I}\) is such that \(\mathrm{Var}[Y(\mathcal{I})]\) is bounded from above for all \(\sigma ( \mathcal{I})\)-measurable random variables \(Y(\mathcal{I})\). They quote from their previous paper [Arch. Ration. Mech. Anal. 239, No. 2, 1025–1060 (2021; Zbl 1456.76134)] the existence and properties of correctors \((\psi _{E},\Sigma _{E})\in H_{\mathrm{loc}}^{1}(\mathbb{R}^{d})^{d}\times L_{\mathrm{loc}}^{2}( \mathbb{R}^{d}\setminus \mathcal{I})\) solution to the infinite-volume corrector problem: \(-\Delta \psi _{E}+\nabla \psi _{E}=0\), \(\operatorname{div}(\psi _{E})=0\) in \(\mathbb{R}^{d}\setminus \mathcal{I}\), \(D(\psi _{E}+Ex)=0\) in \(\mathcal{I} \), \(\int_{\partial I_{n}}\sigma (\psi _{E}+E(x-x_{n}),\Sigma _{E})\nu =0\), \( \int_{\partial I_{n}}\Theta (x-x_{n})\cdot \sigma (\psi _{E}+E(x-x_{n}),\Sigma _{E})\nu =0\) \(\forall n\), \(\forall \Theta \in \mathbb{ M}^{\mathrm{skew}}\). The first main result proves uniform estimates for the \(L^{q}\)-norm of \([(\nabla \psi _{E},\Sigma _{E}1_{\mathbb{R}^{d}\setminus \mathcal{I }})]_{2}\) and \([\psi _{E}(x)]_{2}\), where \([f]_{2}(x)=(\frac{1}{\left\vert B(x)\right\vert }\int_{B(x)}\left\vert f\right\vert ^{2}))^{1/2}\). The authors then prove a quenched large-scale Schauder theory, a quenched large-scale \(L^{p}\)-regularity and an annealed \(L^{p}\)-regularity for the solution \((u_{g},P_{g})\) to the above steady Stokes problem. Finally, they prove a quantitative homogenization result for this steady Stokes problem, giving the expressions of the effective viscosity tensor and of the corrector in the present context. For the proof, the authors establish a Meyers-type perturbative result, they transform the above Stokes and corrector problems in single equations in the whole space \(\mathbb{R}^{d}\) and they prove localized pressure estimates. They also use properties of the Bogovskii operator, a dual Calderon-Zygmund Lemma, Gehring’s Lemma and a reverse Jensen’s inequality.