Optimization on the multinomial manifolds of stochastic non-square matrices

Thanks for note! So I have an optimization problem with simple form but the decision variable is a large-scale matrix. My problem is similar to a existing problem here about multinomial manifolds and my objective function is simpler than the counterpart in that problem. The decision variable $$\mathbf T \in {\bf R}^{1300\times 3100}$$ is a stochastic matrix (each row sums to 1) with non-negative elements. $$\mathbf S \in {\bf R}^{600\times 1300}$$ is a known constant matrix with non-negative elements. $$\mathbf L \in {\bf R}^{3100\times 9000}$$ is a known constant sparse matrix. There are exactly one $$1$$ and one $$-1$$ in each column of $$\mathbf L$$, so the column sum of $$\mathbf L$$ is equal to $$0$$ for every column. Let $$\mathbf X$$ denote the product of above three matrices, i.e. $$\mathbf X=\mathbf S\mathbf T\mathbf L$$. The optimization problem is below. $$\text{minimize}\hspace{3mm}f(\mathbf T)=\sum_{i=1}^{600}\sum_{j=1}^{9000}|\mathrm X_{ij}|$$ $$\text{subject to}\hspace{19mm}\mathbf T{\bf 1}=\bf{1}$$ $$\hspace{39mm}\mathrm T_{ij}>0$$ The optimization function means I expect to obtain a matrix $$\mathbf X$$ as similar to zero matrix as possible, so it can be replaced by minimizing Frobenius norm of $$\mathbf X$$ or other similar forms. As far as I know, my problem resembles optimal transport in some ways. In addition, there may be some trick to reduce the amount of computation and cost of storage given the special nature of matrix $$\mathbf L$$. I'm doubting whether the manopt toolbox in Matlab here is capable of dealing with this large-scale optimization problem (speed is not so important). I have thought that a python package named pymanopt here might work but the multinomial manifolds is not currently supported by pymanopt. I'm planning to employ a feasible algorithm or a computation environment to solve this optimization task. If there are ideas to make my mission possible, what would be the best approach here? Thanks very much for any advice or comments.

New contributor
Hu jiaran is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• You'll want to replace the last constraint by $T_{ij}\ge 0$. Without the inequality, you are optimizing over an open set and you can't guarantee that there is a solution. – Wolfgang Bangerth Nov 22 at 4:58
• Furthermore, this is just a linear program, though a fairly large one. If you introduce slack variables, you'll realize that you can rewrite the problem in terms of the sum of slack variables, and you end up with a whole bunch of inequalities that are all linear in the $X_{ij}$ and because $X$ is linear in your primary variable $T$, these are all linear constraints in $T$. – Wolfgang Bangerth Nov 22 at 5:00
• What you said is very helpful, I need to add slack variables $u_{ij}$, $v_{ij}$, and constraints $u_{ij}\geq 0$, $v_{ij}\geq 0$, $u_{ij}-v_{ij}=X_{ij}$. Furthermore, I should rewrite $f(\mathbf T)$ as $f(\mathbf T)=\sum_{i=1}^{600}\sum_{j=1}^{9000}(u_{ij}+v_{ij})$. – Hu jiaran Nov 22 at 7:53
• I will replace the last constraint by $T_{ij}\geq 0$. Thanks. It now appears that I should look for a software to implement the computation. I will try to do, and any advice for efficient computation is adorable. – Hu jiaran Nov 22 at 8:28
• It's just a linear program after reformulation. There should be plenty of software around for that. – Wolfgang Bangerth 2 days ago