This work concerns single-trace correlations of Euclidean multi-matrix models. In the large-N limit we show that Schwinger-Dyson equations imply loop equations and non-anomalous Ward identities. Loop equations are associated to generic infinitesimal changes of matrix variables (vector fields). Ward identities correspond to vector fields preserving measure and action. The former are analogous to Makeenko-Migdal equations and the latter to Slavnov-Taylor identities. Loop equations correspond to leading large-N Schwinger-Dyson equations. Ward identities correspond to 1/N^2 suppressed Schwinger-Dyson equations. But they become leading equations since loop equations for non-anomalous vector fields are vacuous. We show that symmetries at infinite N persist at finite N, preventing mixing with multi-trace correlations. For one matrix, there are no non-anomalous infinitesimal symmetries. For two or more matrices, measure preserving vector fields form an infinite dimensional graded Lie algebra, and non-anomalous action preserving ones a subalgebra. For Gaussian, Chern-Simons and Yang-Mills models we identify up to cubic non-anomalous vector fields, though they can be arbitrarily non-linear. Ward identities are homogeneous linear equations. We use them with the loop equations to determine some correlations of these models. Ward identities alleviate the underdeterminacy of loop equations. Non-anomalous symmetries give a naturalness-type explanation for why several linear combinations of correlations in these models vanish.