Many combinatorial optimization problems can be reformulated as finding the ground state of the Ising model. Existing Ising solvers are mostly inspired by simulated annealing. Although annealing techniques offer scalability, they lack convergence guarantees and are sensitive to the cooling schedule. We propose solving the Ising problem by relaxing the binary spins to continuous variables and introducing an attraction potential that steers the solution toward binary spin configurations. A key property of this potential is that its combination with the Ising energy produces a Hamiltonian that can be written as a difference of convex polynomials. This enables us to design efficient iterative algorithms that require a single matrix-vector multiplication per iteration and provide convergence guarantees. We implement our Ising solver on a wide range of GPU platforms, from edge devices to high-performance computing clusters, and demonstrate that it consistently outperforms existing solvers across problem sizes ranging from small ($10^3$ spins) to ultra-large ($10^8$ spins).