  We present a many-body $GW$ formalism for quantum subsystems embedded in
discrete polarizable environments containing up to several hundred thousand
atoms described at a fully ab initio random phase approximation level. Our
approach is based on a fragment approximation in the construction of the
Green's function and independent-electron susceptibilities. Further, the
environing fragments susceptibility matrices are reduced to a minimal but
accurate representation preserving low order polarizability tensors through a
constrained minimization scheme. This approach dramatically reduces the cost
associated with inverting the Dyson equation for the screened Coulomb potential
$W$, while preserving the description of short to long-range screening effects.
The efficiency and accuracy of the present scheme is exemplified in the
paradigmatic cases of fullerene bulk, surface, subsurface, and slabs with
varying number of layers.
