Brain function networks (BFN) are widely used in the diagnosis of electroencephalography (EEG)-based major depressive disorder (MDD). Typically, a BFN is constructed by calculating the functional connectivity (FC) between each pair of channels. However, it ignores high-order relationships (e.g., relationships among multiple channels), making it a low-order network. To address this issue, a novel classification framework, based on matrix variate normal distribution (MVND), is proposed in this study. The framework can simultaneously generate high-and low-order BFN and has a distinct mathematical interpretation. Specifically, the entire time series is first divided into multiple epochs. For each epoch, a BFN is constructed by calculating the phase lag index (PLI) between different EEG channels. The BFNs are then used as samples, maximizing the likelihood of MVND to simultaneously estimate its low-order BFN (Lo-BFN) and high-order BFN (Ho-BFN). In addition, to solve the problem of the excessively high dimensionality of Ho-BFN, Kronecker product decomposition is used for dimensionality reduction while retaining the original high-order information. The experimental results verified the effectiveness of Ho-BFN for MDD diagnosis in 24 patients and 24 normal controls. We further investigated the selected discriminative Lo-BFN and Ho-BFN features and revealed that those extracted from different networks can provide complementary information, which is beneficial for MDD diagnosis.