We present a probabilistic Las Vegas algorithm for computing the local zeta function of a hyperelliptic curve of genus $g$ defined over $\mathbb{F}_q$. It is based on the approaches by Schoof and Pila combined with a modeling of the $\ell$-torsion by structured polynomial systems. Our main result improves on previously known complexity bounds by showing that there exists a constant $c > 0$ such that, for any fixed $g$, this algorithm has expected time and space complexity $O((\log q)^{cg})$ as $q$ grows and the characteristic is large enough.