You can think of the $x_{i}^{2}$ as nonnegative weights that sum to one. This opens up the whole world of the generalized mean inequality.

If the matrix isn't sparse, and the cost of getting individual matrix entries is large compared to the cost of accessing an element of a matrix stored in conventional dense matrix form, then iterative methods are going to be horribly slow in practice.

Let me clarify what I meant here- "being able to get an arbitrary element M(i,j) at little cost" isn't very useful. If you don't know where the nonzero elements are in the matrix, then you have to check every single one to find the nonzeros. If you do happen to know where the nonzero elements are, and you can compute them quickly, then you could use this as a way to do matrix vector multiplications in an iterative method.

See ams.org/samplings/feature-column/fcarc-stern-brocot for a discussion of the Stern-Brocot tree which has been used to design optimal gear trains. Tables of optimal gear trains based on this method have been around for a long time.

On the other hand, if this really an unconstrained problem, I'd simply use BFGS. The computations inside the BFGS method won't be time consuming for a small problem like this, so the big concern will be the cost of function and derivative computations. You'll probably want to experiment with finite difference derivatives vs. derivatives by automatic differentiation vs. analytical formulas for the derivatives to see what works best (both in terms of accuracy and speed.) You can also play with the convergence criteria- perhaps a relatively imprecise solution will be adequate for your work.

Now I'm even more confused about your problem. You say that x lives "in a non-euclidean manifold", so it appears that you don't have an unconstrained optimization problem. In that case, you'd have a constrained nonlinear optimization problem and you'd need to consider algorithms for constrained problems rather than methods for unconstrained optimization. A good (but somewhat old) textbook that deals with mathematics and many practical aspects of optimization is "Practical Optimization" by Gill, Murray, and Wright. You'd find lots of good advice in that book on how to approach this.

@Alex; The $J^{T}J$ approximation to the hessian of $c(x)$ works well when $f_{i}(x)-y_{i}$ is small, but if the residuals are large then the approximation can degrade. This is because the second order terms that go into the Hessian are dropped when we approximate it with $J^{T}J$, and those terms have $f_{i}(x)-y_{i}$ factors.

This is confusing- are you minimizing a sum of squares or are you using some other objective (such as the Huber measure)? If you're not minimizing a sum of squares than LM isn't appropriate in the first place. Assuming that you are doing nonlinear least squares, and assuming that you have a reasonable assessment of the uncertainty in the measurements (e.g. $y_{i}$ is known with uncertainty $\sigma_{i}$, then you should normalize by $c(x)=\sum ((f_{i}(x)-y_{i})/\sigma_{i})^{2}$. You should also scale the parameters $x$ so that they're all of similar magnitude.

By the way, it's easy to compute the gradient of $c(x)$ given $J$ and the values of $f_{i}(x)-y_{i}$. You should always check that the gradient is reasonably small as part of your termination criteria.

"It just means that the problem doesn't depend on the corresponding variable" isn't a correct statement. Consider the example of minimizing $(x_{1}^{2}-0)^{2}+(x_{2}^{2}-0)^{2}$. Clearly, $x=0$ is the unique optimal solution, but $J^{T}J$ is 0 at $x=0$.

You haven't mentioned this, but it would seem obvious that you want the additional constraint that $\sum_{b} z_{j}(b)=1$. That is, document j has exactly one rank. Once you've added that constraint, the second constraint $z_{ij}^{1}(a)+z_{ij}^{2}(a)=z_{i}(a)$ is a trivial consequence of the definitions of $z_{ij}(a)^{1}$ and $z_{ij}(a)^{2}$.