# Sparse least-squares fitting of discrete probability distributions

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.

You just want quadratic programming and an implementation which supports sparse matrices. Most established optimisation packages support this, including those of scipy and MATLAB.

I wouldn't bother implementing it yourself unless you really need to do it fast and really don't trust others implementations. I would look at the documentation for the various packages, as the specific algorithm used can usually be specified. You'll almost certainly be able to find a good one, especially if it has been around for a while.

• Please note that the product $A_0A_1$ is a huge dense matrix, but $A_0$ and $A_1$ themselves are sparse. Jun 21, 2013 at 10:10
• Yeah, well I guess all I was saying is give it a go with whatever package you have to hand. The old do first optimise later maxim. That's been my experience with these things, though not a great answer admittedly. Perhaps someone can do better, maybe on cs.SE or the scientific computing one? Have you estimated how long it will take with standard algorithms? Jun 21, 2013 at 12:47
• Thanks for the link to the wikipedia page. Actually, there is a list of quite a few implementations and if I am lucky, I might find what I need. Jun 21, 2013 at 13:30
• @Piostruk do let me know how you get on. I would be interested the approach you end up taking. Jun 25, 2013 at 21:42