  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\'e approximation technique introduced in [6]
following [11]. 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 [6], 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\'e 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 [6] for differential operators with normal and compact
resolvent.
