I am trying to use TFOCS for a problem where I tried interior point methods before, as well as with crude alternating projection based methods. The first class was extremely slow for my problem sizes (in the range of 1e6 unknowns, dense problem) and sometimes suffered numerical issues where very similar formulations seemed to lead the solver to solutions of very different quality. Hand-made alternating projections with some improvised gradient descent worked better, but convergence could be ridiculously slow there as well.

I am fairly sure that my problem is indeed convex, but my objective function (f in my tfocs call) is only defined in R+. As far as I understand, this should still be supported? Several of the included examples for smooth functions seem to have similar restrictions, and based on those I return Inf for function value and NaN for gradient outside of this region.

In practice, I end up with the method terminating with the following message:

“Finished: Step size tolerance reached (||dx||=0)”

Sometimes the solution seems fine, but in purely simulated cases, where I can plug in the correct solution, the objective function value of that solution is clearly superior, with a relative error > 1 %. I’ve verified specifically that the trajectory of linear combinations between my x0 and the correct solution are all within the defined region and that it is numerically convex.

I realized that this could be due to the optimum for some variables being exactly at the boundary of the domain. Any non-zero step along the gradient will bring the function outside the domain. Things can be improved by adding a convex barrier term, at the cost of some optimality. However, I still get early termination with weird results, maybe due to my projection in g interacting in weird ways.

Is the answer simply that such domain restrictions are not allowed? There might be numerical inaccuracy issues at play as well. Initially, I thought that this was all about f(y)-f(y + dy) = 0 and that things might be improved by allowing the smooth function to specifically compute this difference (this would require an API change). However, I am not all that sure anymore that this would help.

I would be happy to provide more details, but some dos and donts for the R+ restriction seems to be the place to start.