# Finding probability vectors from an implicit equation

I have $q$ $n$-dimensional vectors $\vec y_i$ and a matrix $\hat B$ of shape $n\times m$. I'm looking for $q$ $m$-dimensional vectors $\vec x_i$ such that:

1. $\vec y_i=\hat B \vec x_i$
2. each vector $\vec x_i$ represents probability distribution (each entry is non-negative and the $L_1$ norm is 1)

Associated issues are that:

• vectors $\vec y_i$ are very noisy
• matrix $\hat B$ can be estimated with bigger accuracy, but it may be rather sparse

How does one approach such problem?

• Welcome to SciComp.SE! You can formulate this as a simple (set of) quadratic optimization problems with linear constraints, but you can get better results if you have additional knowledge of the desired $x_i$ -- are they independent? Do you have a prior distribution over them? Do they have some specific features (sparse, integral entries)? – Christian Clason Jul 27 '18 at 11:29
• I hope they are independent. But unfortunately I have as little information as stated in the problem. – Pawel Jul 27 '18 at 15:39

I have some ideas, you could try. I will fix one $i$, to simplify the notation. If the $(\vec{x}_i)_{i=1}^q$ are independent, this should be no issue.
1. Least squares: Solve the minimisation problem $$\min \|B\vec{x} - \vec{y}\|_2^2,$$ $$\text{s.t.: } \vec{x} \geq 0, \|\vec{x}\|_1 = 1.$$ This could be done with probably any standard non-linear optimisation library.
2. Bayesian approach: In the Bayesian approach, you would try to compute a posterior measure on the vector $\vec{x}$. You assume you have some prior information about the vector and then incorporate the information given by $\vec{y}$ into this prior information. The prior information you have is that $\vec{x}$ is a probability vector with $m$ entries. Luckily, there are probability measures representing such probability vectors, such as the Dirichlet distribution. Say we have a prior measure $p_0$ on $\mathbb{R}^m$ representing the prior knowledge. For the Bayesian procedure you need to construct a likelihood. If we assume that the noisy data is generated by $$\vec{y} \leftarrow B(\vec{x}_{\rm true}) +\eta,$$ where $\eta \sim \mathrm{N}(0, \Gamma)$ is Gaussian noise. The likelihood would be $$L(\vec{y}|\vec{x}) = \exp(-(1/2)\|\Gamma^{-1/2}(\vec{y} \leftarrow B\vec{x})\|_2^2).$$ You can now use Bayes' rule and for instance importance sampling or Markov chain Monte Carlo to compute the posterior measure. But for that, you should probably look at the Wikipedia pages: Bayesian inference and Markov chain Monte Carlo. The posterior mean would give you the $L^2$-optimal estimator for $\vec{x}_{\rm true}$, the measure itself would quantify the uncertainty around the vectors $\vec{x}$.