I imagine your notation means theta_j is a complex vector in R^M, so that theta_j(i) is its i-th component. And you have N such vectors. if that's correct, than your search space seems to be the set of complex matricex of size M x N whose individual entries have modulus 1.
That is indeed a perfect fit for complexcirclefactory(M, N).
Your problem structure should then look like this:
problem = struct();
problem.M = complexcirclefactory(M, N); % this is the manifold, with your values of M and N
problem.cost = @(X) ... your cost function here, where X is an MxN matrix with unit-modulus entries ... ;
For quick prototyping, you can then try to call
problem = manoptAD(problem);
to see if automatic differentation can figure out your gradient.
This sometimes fails (and sometimes can be fixed with tinkering, but it takes a bit of trial and error).
In any case, eventually you'll want to implement the gradient of your cost function, by defining the field problem.egrad or problem.grad, and then running checkgradient(problem) to make sure it's correct. You can find more info in the tutorial and in my book (pdf on my webpage).