In this paper, based on the paradigm of supplementary variable method (SVM), we reformulate the Klein–Gordon–Zakharov system into equivalent optimization problem subject to PDE constraints, and then present a novel class of high-order energy-preserving numerical algorithms to solve it numerically. The optimization model is discretized by applying the Gauss collocation method as well as the prediction–correction technique in time and the sine pseudo-spectral method in space. Experimental results and some comparisons between this method and other reported ones are provided to favor the suggested method in the overall performance.