In this work, we consider the approximation of Hilbert space-valued meromorphic functions that arise as solution maps of parametric PDEs whose operator is the shift of an operator with normal and compact resolvent, e.g., the Helmholtz equation. In this restrictive setting, we propose a simplified version of the Least-Squares Padé approximation technique studied in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284] following [J. Approx. Theory 95 (1998), pp. 203–2124]. In particular, the estimation of the poles of the target function reduces to a low-dimensional eigenproblem for a Gramian matrix, allowing for a robust and efficient numerical implementation (hence the “fast” in the name). Moreover, we prove several theoretical results that improve and extend those in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284], including the exponential decay of the error in the approximation of the poles, and the convergence in measure of the approximant to the target function. The latter result extends the classical one for scalar Padé approximation to our functional framework. We provide numerical results that confirm the improved accuracy of the proposed method with respect to the one introduced in [ESAIM Math. Model. Numer. Anal. 52 (2018), pp. 1261–1284] for differential operators with normal and compact resolvent.