Consider a function of $N_a$ discrete variables (on the order of 1000), each of which can take one of $N_b$ different states (on the order of 10). The function result is a single, real number. For example, with $N_a=5$ and $N_b=3$:
$f(1,2,3,1,2)=5.234$
$f(1,2,3,3,3)=-1.68$
I am trying to optimize this function, so I am looking into optimization methods. Unfortunately my background is not in mathematics and I don't know exactly what to search for. I am familiar with non-discrete optimization methods like gradient descent and the Metropolis algorithm, etc. I have looked into graph coloring, barely scratching the surface of combinatorial optimization, but I haven't found an example where my problem fits. I'd be obliged if someone could point me in the right direction by telling me what type of problem this is.
Edit:
The underlying problem is kind of a graph-coloring problem; I have a graph of $N_a$ vertices, each one can take a "color" from one of $N_b$ values. The function result is a kind of "interaction energy" depending on the color of each neighbor of each vertex. It's kind of like an Ising model on a semi-random graph. There are some rules on the distribution of "colors", like each color must be represented by at least one vertex. I started out by using Metropolis sampling/simulated annealing, but I wasn't sure if there is a more efficient route to optimization already traced for me, based on purer mathematics. Because the set of "colors" and vertices are discrete, it's a little out of my wheelhouse.