Abstract: It is shown how large N properties of multi-matrix systems can be obtained by minimization of a loop truncated effective Hamiltonian expressed directly in terms of gauge invariant operators. The large N loop space constraints are handled by the use of master variables. For two and three massless Yang-Mills coupled matrices, highly accurate results for large N planar correlators, as well as spectrum, are presented. Generalization to larger number of matrices, relevant for systems such as BFSS, are discussed.