  We present iDARR, a scalable iterative Data-Adaptive RKHS Regularization
method, for solving ill-posed linear inverse problems. The method searches for
solutions in subspaces where the true solution can be identified, with the
data-adaptive RKHS penalizing the spaces of small singular values. At the core
of the method is a new generalized Golub-Kahan bidiagonalization procedure that
recursively constructs orthonormal bases for a sequence of RKHS-restricted
Krylov subspaces. The method is scalable with a complexity of $O(kmn)$ for
$m$-by-$n$ matrices with $k$ denoting the iteration numbers. Numerical tests on
the Fredholm integral equation and 2D image deblurring show that it outperforms
the widely used $L^2$ and $l^2$ norms, producing stable accurate solutions
consistently converging when the noise level decays.
