We study a two-locus model of a quantitative trait with a continuum-of alleles and multilinear epistasis that evolves under mutation, selection, and genetic drift. We derive analytical results based on the so-called House of Gauss approximation for the genetic variance, the mean phenotype, and the mutational variance in the balance of the evolutionary forces. The analytical work is complemented by extensive individual-based computer simulations. We find that (1) analytical results are accurate in a large parameter space; (2) epistasis always reduces the equilibrium genetic variance, as predicted in earlier studies that exclude drift; (3) large-scale stochastic fluctuations and non-equilibrium phenomena like adaptive inertia can strongly influence the evolution of the genetic architecture of the trait.