Closed-form solution to minimize error for one pdf in a weighted sum of pdfs? Let $\theta_3$ be a categorical probability distribution that is the weighted sum of two categorical probability distributions $\theta_1$ and $\theta_2$ plus error.
$$\theta_3=\alpha_1\theta_1+\alpha_2\theta_2+\epsilon$$
$$\sum\alpha=1$$
Assume that $\theta_3, \theta_3, \alpha_1, \alpha_2$ are all known. I want to find the probability distribution $\theta_1$ which minimizes the error.
Sometimes it is possible that the error can be zero; in this case, there is a trivial solution, e.g:
$$\alpha_1=0.5,\alpha_2=0.5$$
$$\theta_2=[0.25,0.25,0.5]$$
$$\theta_3=[0.5,0.25,0.25]$$
Here, $\theta_1=\frac{\theta_3-\alpha_2\theta_2}{\alpha_1}=[0.75,0.25,0]$
In many cases, it is not possible for the error to be zero, e.g.:
$$\alpha_1=0.5,\alpha_2=0.5$$
$$\theta_2=[0.25,0.25,0.5]$$
$$\theta_3=[0.5,0.4,0.1]$$
I can define a function $\sum (\theta_3-\alpha_2\theta_2-\alpha_1\theta_1)^2$ and find the values of $\theta_1$ that minimize this using a nonlinear programming solver, e.g. in R:
sqerr = function(theta1){
  sum((c(0.5,0.4,0.1)-0.5*c(0.25,0.25,0.5)-0.5*theta1)^2)
}

Rsolnp::solnp(pars=extraDistr::rdirichlet(1,c(1,1,1)),
      fun=sqerr,
      eqfun=function(x){sum(x)},
      eqB=1,
      LB=c(0,0,0),
      UB=c(1,1,1))

This is well and good for small problems like this, but if $\theta$ has thousands of categories, solnp takes way too long to solve this for my purposes (I want to nest it will be nested in another optimization). Is there a closed-form solution to this problem?
 A: Solution 1: Quadratic programming
Your optimization problem is a quadratic program (QP). Using a dedicated QP solver should be faster than your general nonlinear programming solver, because it's explicitly designed for this type of problem. If your current method doesn't take advantage of the analytical gradient/Hessian, this would also slow things down dramatically.
The problem can be written in standard QP form as follows. The function to minimize is:
$$f(\theta_1) = \|\alpha_1 \theta_1 + \alpha_2 \theta_2 - \theta_3\|_2^2$$
where $\|\cdot\|_2^2$ is the squared Euclidean norm (so we're minimizing the sum of squares). Note that the solution won't change if we multiply by a constant, so we can equivalently minimize:
$$\frac{1}{2 \alpha_1^2} f(\theta_1)
= \frac{1}{2} \| \theta_1 - y \|_2^2
= \frac{1}{2} \theta_1^T \theta_1 - y^T \theta_1 + \frac{1}{2} y^T y$$
where $y = \frac{1}{\alpha_1} (\theta_3 - \alpha_2 \theta_2)$ is a known constant. Finally, note that $\frac{1}{2} y^T y$ is constant, so we can drop it from the problem. The QP is then:
$$\min_{\theta_1} \ \frac{1}{2} \theta_1^T \theta_1 - y^T \theta_1 \quad
\quad \text{s.t.} \quad
\begin{matrix} \vec{1}^T \theta_1 = 1 \\ \theta_1 \ge \vec{0} \end{matrix}$$
The constraints ensure that the entries of $\theta_1$ are nonnegative and sum to one (which also implies that no entry is greater than one). This problem can be fed into a standard QP solver.
Solution 2: Euclidean projection
It may be possible to solve the problem even more efficiently, because there's additional structure beyond a generic QP. In particular, the solution can be described geometrically as the Euclidean projection of $y$ onto the probability simplex (the set of vectors whose entries are nonnegative and sum to one). That is, the point $\theta_1$ on the probability simplex that's closest to $y$ in terms of Euclidean distance. The following paper describes an efficient algorithm for solving this problem:
Wang, W., & Carreira-Perpinán, M. A. (2013). Projection onto the probability simplex: An efficient algorithm with a simple proof, and an application. arXiv preprint arXiv:1309.1541.
