To compactly represent a graph signal in the frequency domain, critically sampled biorthogonal wavelet filterbanks have been proposed to decompose signals on bipartite graphs. However, in practice graph signals often reside on general graph structures that are not bipartite. Thus, an original nonbipartite graph must be decomposed into a sequence of bipartite graph approximations, so that the filterbanks can be applied successively for signal decomposition. In this paper, unlike previous proposals that are heuristic in nature, we design new bipartite approximation strategies for model-based and empirically derived probability distributions. In the first case, a signal prior assumes a Gaussian Markov Random Field (GMRF) model parametrized by the original graph. In the second case, beyond the GMRF signal prior, empirical signal observations are available to compute a posterior probability distribution. In both cases, we optimize for energy compaction in the bipartite subgraphs with two criteria: the Kullback–Leibler divergence metric, which encourages preserving the spectral characteristics of the original graph; and multiplicity of eigenvalue at frequency 1 for graph Laplacian, which is the frequency with minimal energy discrimination. The relative importance between the two criteria is determined by a proposed measure that evaluates the degree of mismatch between the signal prior and posterior. Given these criteria, we first design a global numerical optimization, and then propose a local heuristic approach for fast implementation with simplifications leading to local metric computation. Experimental results show that our proposed bipartite subgraph decomposition outperforms competing proposals in terms of energy compaction.