I tried your code and it puzzled me too. My best guess is that perhaps the cost function is not differentiable everywhere and that the optimizer gets stuck in regions where the cost indeed isn't differentiable.
To test this hypothesis, I modified the call to the solver as follows:
options.statsfun = statsfunhelper('mindiag', @(x) min(abs(diag(x*x'))));
[X, xcost, info, options] = trustregions(problem, , options);
This makes it so that, at each iteration, the solver records the minimum absolute value of the diagonal entries of x*x', and stores that number in the info struct-array.
Now, display the results as:
and you will see that this quantity becomes dangerously close to zero. Unfortunately, that means x comes too close to a region where the cost function does not have a gradient (because your cost function divides by diag(x*x')).
You may want to reconsider the cost function to make it smooth. Alternatively, you could try nonsmooth solvers; I tried pso and neldermead, but they don't work on Stiefel (need a distance or a log function, which aren't known); I tried rlbfgs which often behaves well on nonsmooth functions, but here it goes to points with very small values on diag(x*x') even faster -- but perhaps that's not such a bad thing from the optimization point of view; this I cannot assess.