Compositionally complex or high entropy alloys have been extensively investigated for a plethora of applications due to their outstanding thermal, 2 mechanical, energy storage, and catalytic properties, but rational design of these materials is challenged by the combinatorial explosion of the compositional space from which they can be synthesized. Accurate and tractable modeling strategies are therefore critical for discovering and/or optimizing high entropy alloys to direct the synthesis and testing of promising materials and conserve expensive experimental resources. Crucially, any modeling approach must be accurate and data-efficient in treating solid solutions with high compositional complexity. Another critical consideration is that the model's input needs to be trivial-to-generate such that hundreds of thousands or millions of energy calculations can then be performed. This enables convergence of materials' property predictions through first principles simulations like Monte Carlo or direct statistical sampling. The introduction of the cluster expansion (CE) method often satisfies the aforementioned requirements and its initial applications in the prediction of phase stability of complex alloys marked a significant milestone in the ability to rationally design complex inorganic materials and impact a host of technological applications. Its widespread popularity mainly stems from its ability to accurately predict complex energy landscapes in materials exhibiting significant configurational degrees of freedom while requiring only a relatively small (or at least tractable) number of expensive first-principles electronic structure calculations. The cluster expansion formalism has enabled first-principles accuracy calculations of temperature-composition phase diagrams, activation barriers, and even electronic excitations 18 across diverse materials applications, ranging from hydrogen storage to battery 20 to electronic device materials. Furthermore, along with the recent explosion of machine learning (ML) to accelerate various tasks in materials science, ML improvements on the original cluster expansion formalism have been proposed to ameliorate some of its shortcomings. Parallely, in recent years, graph neural networks (GNNs) have undergone widespread adoption in ma-terials science thanks to their automated feature extraction capabilities. They generally remove the need to manually derive input features (e.g., the process of defining the cluster functions in CE) that suitably describe a material of high compositional complexity. They can serve as direct surrogate models for highly demanding computational tasks (e.g., bulk modulus and equation-of-state fitting, supercell defect relaxations, or phonon density of state calculations ) and more recently have been developed towards "universal" force fields to enable firstprinciples accuracy simulations across the periodic table. Here we demonstrate how minimal-complexity GNNs (i.e., ∼2,000 parameter models) with identical architectures and training hyperparameters are applicable to relaxed formation energy predictions from un-relaxed crystallographic representations to aid the design of compositionally complex materials with limited training data. For example, this approach permits the high-throughput screening needed to predict and therefore compositionally optimize phase stability of low density high-entropy alloys and hydrogen equilibria in metal hydrides. Critically, these models' predictions extrapolate to configurations with formation energies outside the range of target values in the training data, and small training datasets can be rationally and efficiently constructed to maximize performance. In the future, we expect such an approach to be useful in any number of applications that CE has typically been used to study, while alleviating challenges of constructing cluster functions in high compositional complexity systems.