  $f(R)$ gravity, capable of driving the late-time acceleration of the
universe, is emerging as a promising alternative to dark energy. Various $f(R)$
gravity models have been intensively tested against probes of the expansion
history, including type Ia supernovae (SNIa), the cosmic microwave background
(CMB) and baryon acoustic oscillations (BAO). In this paper we propose to use
the statistical lens sample from Sloan Digital Sky Survey Quasar Lens Search
Data Release 3 (SQLS DR3) to constrain $f(R)$ gravity models. This sample can
probe the expansion history up to $z\sim2.2$, higher than what probed by
current SNIa and BAO data. We adopt a typical parameterization of the form
$f(R)=R-\alpha H^2_0(-\frac{R}{H^2_0})^\beta$ with $\alpha$ and $\beta$
constants. For $\beta=0$ ($\Lambda$CDM), we obtain the best-fit value of the
parameter $\alpha=-4.193$, for which the 95% confidence interval that is
[-4.633, -3.754]. This best-fit value of $\alpha$ corresponds to the matter
density parameter $\Omega_{m0}=0.301$, consistent with constraints from other
probes. Allowing $\beta$ to be free, the best-fit parameters are $(\alpha,
\beta)=(-3.777, 0.06195)$. Consequently, we give $\Omega_{m0}=0.285$ and the
deceleration parameter $q_0=-0.544$. At the 95% confidence level, $\alpha$ and
$\beta$ are constrained to [-4.67, -2.89] and [-0.078, 0.202] respectively.
Clearly, given the currently limited sample size, we can only constrain $\beta$
within the accuracy of $\Delta\beta\sim 0.1$ and thus can not distinguish
between $\Lambda$CDM and $f(R)$ gravity with high significance, and actually,
the former lies in the 68% confidence contour. We expect that the extension of
the SQLS DR3 lens sample to the SDSS DR5 and SDSS-II will make constraints on
the model more stringent.
