Abstract We construct the general permutation invariant Gaussian two-matrix model for matrices of arbitrary size D . The parameters of the model are given in terms of variables defined using the representation theory of the symmetric group S D . A correspondence is established between the permutation invariant polynomial functions of the matrix variables (the observables of the model) and directed colored graphs, which sheds light on stability properties in the large D counting of these invariants. A refined counting of the graphs is given in terms of double cosets involving permutation groups defined by the local structure of the graphs. Linear and quadratic observables are transformed to an S D representation theoretic basis and are used to define the convergent Gaussian measure. The perturbative rules for the computation of expectation values of graph-basis observables of any degree are given in terms of the representation theoretic parameters. Explicit results for a number of observables of degree up to four are given along with a Sage programme that computes general expectation values.