Let the set of feasible solution be the set of all row-stochastic $n \times k$ matrices $P = [p_{ij}]$, that is $\mathcal{P} := \{P \in \mathbb{R}^{n \times k} \ | \ P \mathbf{1} = \mathbf{1}, P \geq 0\}$, where $\mathbf{1}$ denotes the all-ones-vector and $\geq$ is applied element-wise.
Further let $C = [c_{ij}] \in \mathbb{R}^{n \times n}$ be given.
How to solve the following, in general non-convex, maximization problem efficiently?
$$\arg\max_{P \in \mathcal{P}} \sum_{i,j} c_{ij} \cdot \sum_{t} p_{it} \cdot p_{jt} \tag{$\star$}$$
This can be re-stated with vectors $p_i := (p_{i1} \ldots p_{ik})^T$, and $\left<\cdot, \cdot\right>$ the standard inner product, as
$$\arg\max_{P \in \mathcal{P}} \sum_{i,j} c_{ij} \cdot \left<p_i, p_j\right>$$
This reminds me to semi-definite-programming, and also to the downhill-simplex-algorithm, but none of these approaches matches exactly. For example, one can also state it as
$$\arg\max_{Q} \sum_{i,j} c_{ij} \cdot q_{ij}$$
where $Q$ is restricted to be positive semi-definite (i.e., a Gram matrix, that is $Q=X^T X$ for some matrix $X$), plus being further restricted to satisfy $Q = P P^T$ for row-stochastic $P \in \mathcal{P}$. This additional restriction on $Q$'s implicit embedding (rows of $P$ to stay within the $(k-1)$-Simplex) seems not to transform into linear constraints on $Q$ (?), but probably it even allows to simplify the problem analytically? If it helps, one might alternatively consider the substochastic case $P\mathbf{1} \leq \mathbf{1}$, or some alternating maximization approach...
So, can this quadratic problem $(\star)$ be transformed into a simpler form? How would you try to solve it numerically for, say, $n\leq400, k\leq50$? Is there some non-convex global optimization strategy that applies here, or do I have to apply a heuristic?