We consider a class of matrix integrals over the unitary group $U(N)$ with an infinite set of couplings characterized by a series $f(q)=\sum_{n \ge 1} a_n q^n$, with $a_n \in \mathbb{Z}$. Such integrals arise in physics as the partition functions of free four-dimensional gauge theories on $S^3$ and, in particular, as the superconformal index of super Yang-Mills theory. We show that any such model can be expressed in terms of a system of free fermions in an ensemble parameterized by the infinite set of couplings. Integrating out the fermions in a given quantum state leads to a convergent expansion as a series of determinants, as shown by Borodin-Okounkov many years ago. By further averaging over the ensemble, we obtain a formula for the matrix integral as a $q$-series with successive terms suppressed by $q^{αN + β}$ where $α$, $β$ do not depend on $N$. This provides a matrix-model explanation of the giant graviton expansion that has been observed recently in the literature.