We propose an algebraic method for extracting excited states from lattice gauge theory correlation functions. Instead of fitting to a sum of decaying exponentials, we adopt a variant of Prony's method to obtain $M$ energies (exponential decay rates) by finding the roots of an $M^{\rm th}$ order polynomial, and then solving for the amplitudes linearly. The resulting states tend to have overlapping error ellipses, making identification of states ambiguous. This is especially problematic at large Euclidean times where the signal to noise may be low, as well as when many states are extracted. We propose a variation of K-means clustering to identify each extracted state.