I have an optimization problem involving $n$ discrete probability distributions and I am looking for a suitable solution for this problem that I can implement as computer program.
Let the vector $\mathbf{x} \in \mathbb{R}^N$ hold the discrete probabilities of all probability distributions and let $S_k$ with $k \in \{1, \ldots, n\}$ hold the indices refering to elements of $\mathbf{x}$ for the $k^{\textrm{th}}$ distribution.
The problem is the following:
$\underset{\mathbf{x} \in \mathbb{R}^N}{\text{Minimize}} \quad \left\| A_0 A_1 \mathbf{x} - B \right\|_2^2 $
$\text{Subject to} \quad \forall i \in \{1, .., N\} : x_i \geqslant 0 \quad \text{and} \quad \forall k \in \{1, \ldots, n\} : \sum_{j \in S_k} x_j = 1$
This is a linear least-squares problem subject to linear inequality and equality constraints. These constraints simply enforce the probabilities to be valid, that is, every probability should be non-negative and the total probability for every probability distribution is 1.
The matrix $A_0$ is sparse and has $M$ rows. The matrix $A_1$ is also sparse and have $N$ columns. Their product $A_0A_1$, however, is dense, because $A_0$ has a couple of dense columns for which the corresponding rows of $A_1$ are also dense. Since both $M$ and $N$ are in the order of 100000, computing $A_0A_1$ explicitly and storing it in memory is inconvenient. Therefore, computing $A_0\cdot(A_1 \cdot \mathbf{x})$ is efficient, whereas computing $(A_0\cdot A_1 )\cdot \mathbf{x}$ is highly inefficient. The matrix $B$ has dimensions $M \times 1$.
Is there a suitable algorithm for solving the problem above? I do not want to reinvent the wheel trying to solve this optimization problem. For instance, there may be similar optimization problems in finance, structural mechanics, machine learning, fluid dynamics etc, for which there exists efficient algorithms that I am not aware of.
What I have looked at:
In the unconstrained case, the LSQR algoritm would be suitable, since it does not require $A_0A_1$ to be known and for many objective functions where the non-negativity constraints never become active, this algorithm would work, since the remaining problem (least-squares and linear equality constraints) can be reformulated as an unconstrained one. In my case, however, the inequality constraints will become active.
In the constrained case, the best option I can think of would be the SMO algorithm that is used to train support vector machines. It solves a problem similar to the problem that I would like to solve. I am still worried that it will be very time-consuming to find an acceptable solution.
