  Integral equations are commonly encountered when solving complex physical
problems. Their discretization leads to a dense kernel matrix that is block or
hierarchically low-rank. This paper proposes a new way to build a low-rank
factorization of those low-rank blocks at a nearly optimal cost of
$\mathcal{O}(nr)$ for a $n \times n$ block submatrix of rank r. This is done by
first sampling the kernel function at new interpolation points, then selecting
a subset of those using a CUR decomposition and finally using this reduced set
of points as pivots for a RRLU-type factorization. We also explain how this
implicitly builds an optimal interpolation basis for the Kernel under
consideration. We show the asymptotic convergence of the algorithm, explain his
stability and demonstrate on numerical examples that it performs very well in
practice, allowing to obtain rank nearly equal to the optimal rank at a
fraction of the cost of the naive algorithm.
