A systematic method is presented to obtain a large class of constraint equations for matrix models at finite N. These constraints are associated with the higher order differential operators with respect to the eigenvalues of the matrices, {{ λ n ∂ λ m }} with n, m⩾0. In the case of the one-matrix model, we find that the constraints for lower m are reducible to the Virasoro constraints. We derive a class of new constraint equations for the matrix chains.