# Optimization on the manifold of stochastic matrices

So I have an optimization problem of the form $$\text{maximize}\hspace{3mm}f(A):{\bf R}^{K\times K}\rightarrow{\bf R}$$ $$\text{subject to}\hspace{19mm}A^T{\bf 1}=\bf{1}$$ $$\hspace{33mm}A\geq 0$$ where the constraint set is thus the manifold of left stochastic matrices (each column sums to 1). The objective function $f$ is not concave and thus my goal is to very quickly find a local maxima (speed is important).

The function $f$ is a nasty piece of work, nevertheless the gradient $\nabla f$ is analytically computable. However even without the constraint, setting the gradient equal to zero and solving for the extrema is intractable, and thus I can't imagine that back solving for Lagrange multipliers is doable.

The method I've tried is manifold gradient ascent using some code found here from the Manopt package. They refer to the manifold of left stochastic matrices as the multinomial manifold, and the code provides methods for basically projecting the gradient onto the tangent space, and then retracting this new gradient onto the manifold itself.

Unfortunately as far as retractions go it's fairly primitive I think, and it involves entry-wise dividing the projected gradient by the current point and then exponentiating it, which causes overflow unless you make your step-sizes quite small (however I don't fully understand the intricacies of manifold optimization so I could be wrong about this).

I could compute the Hessian numerically and try to run a 2nd order optimization on the multinomial manifold, but I'd have to do some reading to figure out exactly how to do this. Basically I need something that is fast, speed is more important than accuracy, as I have to do this thousands of times and it's only 1 step in a larger coordinate ascent algorithm.

The gradient $\nabla f$ is complicated and thus each evaluation of it is quite costly, however the one saving grace is that the dimension of $A$ is relatively small, almost certainly $<100$.

What would be the best approach here?

Update

As requested, here is the function and its gradient: \begin{align} f(A )&= \Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K} \cdot \big(\sum_{n=2}^NA^{n-2}\big)\cdot {\bf x}\\ &+ \sum_{n=1}^N[\log B_{in}]_{i=1:K}^T\cdot A^{n-1}\cdot{\bf x} \end{align}

\begin{align} \nabla f(A) &= [\psi(C_{ij}) - \log A_{ij} -1]_{i=1:K,\hspace{1mm}j=1:K}\cdot \text{diagm}\Big(\big(\sum_{n=2}^NA^{n-2}\big)\cdot{\bf x}\Big)\\ &+ \sum_{n=3}^{N}\sum_{r=0}^{n-3}\Big((A^r)^T\cdot\Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K}^T \cdot {\bf x}^T\cdot (A^{n-3-r})^T\Big)\\ &+ \sum_{n=2}^N\sum_{r=0}^{n-2}(A^r)^T\cdot[\log B_{in}]_{i=1:K}\cdot{\bf x}^T\cdot(A^{n-2-r})^T \end{align} Where $\psi$ is the digamma function, $A$ and $C$ are $K\times K$ matrices, $B$ a $K\times N$ matrix and ${\bf x}$ is a length $K$ vector, all real-valued. The function $\text{diagm}$ converts a vector to a diagonal matrix. Also I checked my derivation of the gradient numerically, so you can be sure its correct.

I'm writing this in Julia, and I tried using the NLopt package both using a gradient-based method and a derivative free method, but they were even slower than my original manifold approach.

• All you've told us about your objective function f is that is it not concave and is a "nasty piece of work" whose gradient is analytically computable. Perhaps more information on f would facilitate more meaningful responses. Regardless of the terminology "manifold of left stochastic matrices", you have a linearly constrained problem. There are plenty of local solvers for linearly constrained nonlinear objective problems. As to whether this is a good problem to be solving within a larger algorithm and what attributes its solution should have, that is another story. – Mark L. Stone May 9 '16 at 12:57
• @MarkL.Stone, as per your request I added $f$ and $\nabla f$. – Thoth May 9 '16 at 23:59
• Are the problem inputs for the thousands of times you have to solve the problem closely related? For example, can you use the optimal A from a previous problem as a good starting value for the current problem? If so, you might get massive speedups on the subsequent problems. – Mark L. Stone May 10 '16 at 1:36
• It looks like f(A) is a matrix polynomial. Is that right? – Mark L. Stone May 10 '16 at 1:37
• I would put the semidefinite constraint into the objective using a log-det barrier, and then solve the resulting equality-constrained nonlinear optimization problem using any standard technique (e.g. augmented Lagrangian, sequential quadratic programming). Since the dimensionality of the problem is so small, you can afford to use a second-order (i.e. Newton-based) method, even if you have to compute the Hessians numerically. – Richard Zhang May 12 '16 at 15:11