In this paper, we propose a numerical method to solve stochastic elliptic interface
problems with random interfaces. Shape calculus is first employed to derive the
shape-Taylor expansion in the framework of the asymptotic perturbation approach. Given the
mean field and the two-point correlation function of the random interface, we can thus
quantify the mean field and the variance of the random solution in terms of certain orders
of the perturbation amplitude by solving a deterministic elliptic interface problem and
its tensorized counterpart with respect to the reference interface. Error estimates are
derived for the interface-resolved finite element approximation in both, the physical and
the stochastic dimension. In particular, a fast finite difference scheme is proposed to
compute the variance of random solutions by using a low-rank approximation based on the
pivoted Cholesky decomposition. Numerical experiments are presented to validate and
quantify the method.