The representation of ground states of fermionic quantum impurity problems as superpositions of Gaussian states has recently been given a rigorous mathematical foundation. [S. Bravyi and D. Gosset, Comm. Math. Phys. 356, 451 (2017)]. It is natural to ask how many parameters are required for an efficient variational scheme based on this representation. An upper bound is O(N 2 ), where N is the system size, which corresponds to the number parameters needed to specify an arbitrary Gaussian state. We provide an alternative representation, with more favorable scaling, only requiring O(N ) parameters, that we illustrate for the interacting resonant level model. We achieve the reduction by associating mean-field-like parent Hamiltonians with the individual terms in the superposition, using physical insight to retain only the most relevant channels in each parent Hamiltonian. We benchmark our variational ansatz against the Numerical Renormalization Group, and compare our results to existing variational schemes of a similar nature to ours. Apart from the ground state energy, we also study the spectrum of the correlation matrix -a very stringent measure of accuracy. Our approach outperforms some existing schemes and remains quantitatively accurate in the numerically challenging near-critical regime.