  Low frequency gravitational waves (GWs) are keys to understanding
cosmological inflation and super massive blackhole (SMBH) formation via
blackhole mergers, while it is difficult to identify the low frequency GWs with
ground-based GW experiments such as the advanced LIGO (aLIGO) and VIRGO due to
the seismic noise. Although quasi-stellar object (QSO) proper motions produced
by the low frequency GWs are measured by pioneering studies of very long
baseline interferometry (VLBI) observations with good positional accuracy, the
low frequency GWs are not strongly constrained by the small statistics with 711
QSOs (Darling et al. 2018). Here we present the proper motion field map of
400,894 QSOs of the Sloan Digital Sky Survey (SDSS) with optical {\it Gaia}
EDR3 proper motion measurements whose positional accuracy is $< 0.4$
milli-arcsec comparable with the one of the radio VLBI observations. We obtain
the best-fit spherical harmonics with the typical field strength of
$\mathcal{O}(0.1)\, \mu$arcsec, and place a tight constraint on the energy
density of GWs, $\Omega_{\rm gw}=(0.964 \pm 3.804) \times 10^{-4}$ (95 \%
confidence level), that is significantly stronger than the one of the previous
VLBI study by two orders of magnitude at the low frequency regime of $f
<10^{-9}\,{\rm [Hz]}\simeq (30\,{\rm yr})^{-1}$ unexplored by the pulsar timing
technique. Our upper limit rules out the existence of SMBH binary systems at
the distance $r < 400$ kpc from the Earth where the Milky Way center and local
group galaxies are included. Demonstrating the limit given by our optical QSO
study, we claim that astrometric satellite data including the forthcoming {\it
Gaia} DR5 data with small systematic errors are powerful to constrain low
frequency GWs.
