An Algorithm for Finding Positive
Solutions to Polynomial Equations

It suffices to prove this when $a$ and $b$ are scalars. We let $t=b/a$, and,

$$D\left(a\|b\right)=b\log\left(\frac{b}{a}\right)-b+a=a(t\log t-t+1)=a\int_{1}^{t}\log s\,ds.$$

It is clear that $\int_{1}^{t}\log s\,ds$ is non-negative and equal to zero if and only if $t=1$. ∎

Again, we assume that $a$ and $b$ are scalars and then it becomes a straightforward computation. ∎

We let $A=\sum_{i=1}^{m}a_{i}$ and $B=\sum_{i=1}^{m}b_{i}$, and apply Lemma 4 to the last term:

	$\displaystyle D\left(\frac{B}{A}a\|b\right)$	$\displaystyle=D\left(a\|b\right)+\left(\frac{B}{A}-1\right)A-B\log\frac{B}{A}$
		$\displaystyle=D\left(a\|b\right)-B\log\frac{B}{A}+B-A=D\left(a\|b\right)-D\left(A\|B\right).$

After rearranging, we get the desired expression. ∎

First, since the statement only depends on the $j$th row of the matrix $g$, we can assume that $g$ is a row vector and we drop $j$ from future subscripts. Second, we can assume that $d=1$ by replacing $g_{i}$ with $g_{i}/d$.

Third, we reduce to the case when $g_{i}=1$ for all $i$. We define a new set of exponents $\tilde{\alpha}$ and coefficients $\tilde{a}_{\tilde{\alpha}}$ by $\tilde{\alpha}_{i}=g_{i}\alpha_{i}$ and $\tilde{a}_{\alpha}=a_{\alpha}\prod x_{i}$, where the product is taken over all indices $i$ such that $g_{i}=0$. We take $\tilde{x}$ to be a vector indexed by those $i$ such that $g_{i}\neq 0$. Then, under the change of coordinates $\tilde{x}_{i}=x_{i}^{1/g_{i}}$, we have $a_{\alpha}x^{\alpha}=\tilde{a}_{\alpha}\tilde{x}^{\tilde{\alpha}}$ and the update rule in (4) is the same for the new system with coefficients $\tilde{a}_{\tilde{\alpha}}$ and exponents $\tilde{\alpha}$. Furthermore, if all entries of $\tilde{\alpha}$ are zero, then $\tilde{x}^{\tilde{\alpha}}=1$ for all vectors $x$ and so we can drop $\tilde{\alpha}$ from our exponent set. Therefore, for the rest of the proof, we drop the tildes, and assume that $\sum_{i}\alpha_{i}=1$ for all $\alpha\in S$ and $g_{i}=1$ for all $i$, in which case $g$ drops out of the equations.

To prove the desired inequality, we substitute the updated assignment $x^{\prime}$ into the definition of Kullback-Leibler divergence:

	$\displaystyle D\left(a_{\alpha}(x^{\prime})^{\alpha}\|w_{\alpha}\right)$	$\displaystyle=w_{\alpha}\log\left(\frac{w_{\alpha}}{a_{\alpha}x^{\alpha}\prod_{i=1}^{n}\left(\frac{\sum_{\beta}w_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}\right)^{\alpha_{i}}}\right)-w_{\alpha}+a_{\alpha}(x^{\prime})^{\alpha}$
		$\displaystyle=w_{\alpha}\log\frac{w_{\alpha}}{a_{\alpha}x^{\alpha}}-\sum_{i=1}^{n}\alpha_{i}w_{\alpha}\log\frac{\sum_{\beta}\beta_{i}w_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}-w_{\alpha}+a_{\alpha}(x^{\prime})^{\alpha}$
(7)			$\displaystyle=D\left(a_{\alpha}x^{\alpha}\|w_{\alpha}\right)-\sum_{i=1}^{n}\alpha_{i}w_{\alpha}\log\frac{\sum_{\beta}\beta_{i}w_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}-a_{\alpha}x^{\alpha}+a_{\alpha}(x^{\prime})^{\alpha}.$

On the other hand, let $C$ denote the last term of (6), which we can expand as,

	$\displaystyle C$	$\displaystyle=\sum_{i=1}^{n}D\left(\sum_{\alpha}\alpha_{i}a_{\alpha}x^{\alpha}\|\sum_{\alpha}\alpha_{i}w_{\alpha}\right)$
		$\displaystyle=\sum_{i=1}^{n}\left(\left(\sum_{\alpha}\alpha_{i}w_{\alpha}\right)\log\frac{\sum_{\alpha}\alpha_{i}w_{\alpha}}{\sum_{\alpha}\alpha_{i}a_{\alpha}x^{\alpha}}-\sum_{\alpha}\alpha_{i}w_{\alpha}+\sum_{\alpha}\alpha_{i}a_{\alpha}x^{\alpha}\right)$
		$\displaystyle=\sum_{i=1}^{n}\sum_{\alpha}\left(\alpha_{i}w_{\alpha}\log\frac{\sum_{\beta}\beta_{i}w_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}-\alpha_{i}w_{\alpha}+\alpha_{i}a_{\alpha}x^{\alpha}\right)$
(8)			$\displaystyle=\sum_{\alpha}\left(\sum_{i=1}^{n}\alpha_{i}w_{\alpha}\log\frac{\sum_{\beta}\beta_{i}w_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}\right)-w_{\alpha}+a_{\alpha}x^{\alpha},$

where the last step follows from the assumption that that $\sum_{i}\alpha_{i}=1$ for all $\alpha\in S$. We take the sum of (7) over all $\alpha\in S$ and add it to (8) to get,

(9)

$$\sum_{\alpha}D\left(a_{\alpha}(x^{\prime})^{\alpha}\|w_{\alpha}\right)+C=\sum_{\alpha}D\left(a_{\alpha}x^{\alpha}\|w_{\alpha}\right)-\sum_{\alpha}b_{\alpha}+\sum_{\alpha}a_{\alpha}(x^{\prime})^{\alpha}.$$

Finally, we expand the last term of (8) using the definition of $x^{\prime}$ and apply the arithmetic-geometric mean inequality,

	$\displaystyle\sum_{\alpha}a_{\alpha}(x^{\prime})^{\alpha}$	$\displaystyle=\sum_{\alpha}a_{\alpha}x^{\alpha}\prod_{i=1}^{n}\left(\frac{\sum_{\beta}\beta_{i}b_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}\right)^{\alpha_{i}}$
		$\displaystyle\leq\sum_{\alpha}a_{\alpha}x^{\alpha}\sum_{i=1}^{n}\alpha_{i}\frac{\sum_{\beta}\beta_{i}b_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}$
		$\displaystyle=\sum_{i=1}^{n}\left(\sum_{\alpha}\alpha_{i}a_{\alpha}x^{\alpha}\right)\frac{\sum_{\beta}\beta_{i}b_{\beta}}{\sum_{\beta}\beta_{i}a_{\beta}x^{\beta}}$
		$\displaystyle=\sum_{i=1}^{n}\sum_{\beta}\beta_{i}b_{\beta}=\sum_{\beta}b_{\beta}.$

Together with (9), this gives the desired inequality. ∎

For the update rule to be constant means that the numerator and denominator in (4) are equal, i.e.

(10)

$$\sum_{\alpha}\alpha_{i}g_{ji}a_{\alpha}x^{\alpha}=\sum_{\alpha}\alpha_{i}g_{ji}w_{\alpha}\quad\mbox{for all $i$ and $j$.}$$

By our assumption on $g$, for each $i$, some $g_{ji}$ is non-zero, so (10) is equivalent to

(11)

$$\sum_{\alpha}\alpha_{i}a_{\alpha}x^{\alpha}=\sum_{\alpha}\alpha_{i}w_{\alpha}\quad\mbox{for all }i.$$

On the other hand, we compute the partial derivative

$$\frac{\partial}{\partial x_{i}}\sum_{\alpha}D\left(a_{\alpha}x^{\alpha}\|w_{\alpha}\right)=\sum_{\alpha}-w_{\alpha}\frac{\alpha_{i}}{x_{i}}+\alpha_{i}a_{\alpha}\frac{x^{\alpha}}{x_{i}}.$$

Since each $x_{i}$ is assumed to be non-zero, it is clear that all partial derivatives being zero is equivalent to (11). ∎

We consider

(12)

$$\sum_{i,\alpha}D\left(a_{i\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right)\qquad\mbox{where }w_{i\alpha}=\frac{b_{i}a_{i\alpha}x^{\alpha}}{\sum_{\beta}a_{i\beta}x^{\beta}},$$

and apply Lemma 5 in two different ways. First, by applying Lemma 5 to each group of (12) with fixed $\alpha$, we get

$$\sum_{i,\alpha}D\left(a_{i\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right)=\sum_{\alpha}D\left(a_{\alpha}(x^{\prime})^{\alpha}\|w_{\alpha}\right)+\sum_{i,\alpha}D\left(\frac{\sum_{j}w_{j\alpha}}{\sum_{j}a_{\alpha}(x^{\prime})^{\alpha}}a_{i\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right).$$

In the last term, the monomials $(x^{\prime})^{\alpha}$ cancel and so it is a constant independent of $x^{\prime}$ which we denote $E$. On the other hand, we can apply Lemma 5 to each group in (12) with fixed $i$,

$$\sum_{i,\alpha}D\left(a_{i\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right)=\sum_{i}D\left(\sum_{\alpha}a_{i\alpha}(x^{\prime})^{\alpha}\|b_{i}\right)+\sum_{i,\alpha}D\left(\frac{b_{i}a_{i\alpha}(x^{\prime})^{\alpha}}{\sum_{\beta}a_{i\beta}(x^{\prime})^{\beta}}\|w_{i\alpha}\right).$$

We can combine these equations to get

(13)

$$\sum_{i}D\left(\sum_{\alpha}a_{i\alpha}(x^{\prime})^{\alpha}\|b_{i}\right)=\sum_{\alpha}D\left(a_{\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right)+E-\sum_{i,\alpha}D\left(\frac{b_{i}a_{i\alpha}(x^{\prime})^{\alpha}}{\sum_{\beta}a_{i\beta}(x^{\prime})^{\beta}}\|w_{i\alpha}\right).$$

By Proposition 3, the last term of 13 is non-negative, and by the definition of $w_{i\alpha}$, it is zero for $x^{\prime}=x$. Therefore, any value of $x^{\prime}$ which decreases the first term compared to $x$ will also decrease the left hand side by at least as much, which i sthe desired inequality.

In order to prove the statement about the derivative, we consider the derivative of (13) at $x^{\prime}=x$. Because the last term is mimimized at $x^{\prime}=x$, its derivative is zero, so

$$\left.\frac{\partial}{\partial x_{j}^{\prime}}\right|_{x^{\prime}=x}\sum_{i}D\left(\sum_{\alpha}a_{i\alpha}(x^{\prime})^{\alpha}\|b_{i}\right)=\left.\frac{\partial}{\partial x_{j}^{\prime}}\right|_{x^{\prime}=x}\sum_{i,\alpha}D\left(a_{i\alpha}(x^{\prime})^{\alpha}\|w_{i\alpha}\right).$$

By Proposition 7, a positive vector $x$ is a fixed point of the inner loop if and only if these partial derivatives on the right are zero for all indices $j$, which is the definition of a critical point. ∎

The Kullback-Leibler divergence $\sum_{\alpha}D\left(a_{\alpha}x^{\alpha}\|w_{\alpha}\right)$ decreases at each step of the inner loop by Lemma 6. Thus, by Lemma 8, the divergence

(14)

$$\sum_{i=1}^{n}D\left(\sum_{\alpha}a_{i\alpha}x^{\alpha}\|b_{i}\right)$$

decreases at least as much. However, the divergence (14) is non-negative according to Proposition 3. Therefore, the magnitude of the decreases in divergence must approach zero over the course of the algorithm. By Lemma 6, this means that the quantity $C$ in that theorem must approach zero. By Proposition 3, this means that the quantities in that divergence approach each other. However, up to a power, these are the numerator and denominator of the factor in the update rule (4), so the difference between consecutive vectors $x$ approaches zero.

Thus, we just need to show that $x$ remains bounded. However, since some power of each variable $x_{i}$ occurs in some equation, as $x_{i}$ gets large, the divergence for that equation also gets arbitrarily large. Therefore, each $x_{i}$ must remain bounded, so the vector $x$ must have a limit as the algorithm is iterated. If this limit is in the interior of the positive orthant, then it must be a fixed point. By Lemma 8 and Proposition 7, this fixed point must be a critical point of the divergence (5). ∎

We write our system of equations as $\sum_{\alpha\in S}a_{i\alpha}x^{\alpha}=0$ for $1\leq i\leq\ell$, where $S\subset\mathbb{N}^{n}$ is an arbitrary finite set of exponents and $a_{i\alpha}$ are any real numbers. We let $d$ be the maximum degree of any monomial $x^{\alpha}$ for $\alpha\in S$. We homogenize the equations with a new variable $x_{n+1}$. Explicitly, define $S^{\prime}\subset\mathbb{N}^{n+1}$ to consist of $\alpha^{\prime}=(\alpha,d-\sum_{i}\alpha_{i})$ for all $\alpha$ in $S$ and we write $a_{i\alpha^{\prime}}=a_{i\alpha}$. We add a new equation with coefficients $a_{\ell+1,\alpha}=1$ for all $\alpha\in S^{\prime}$ and $b_{\ell+1}=1$. For this system, we can clearly take $g_{1i}=1$ and $d_{1}=d$ to satisfy the condition on exponents. Furthermore, for any positive solution $(x_{1},\ldots,x_{n})$ to the original system of equations, $(x_{1}^{\prime},\ldots,x_{n+1}^{\prime})$ with $x_{i}^{\prime}=x_{i}/\big{(}\sum_{\alpha}x^{\alpha}\big{)}^{1/d}$ and $x_{n+1}^{\prime}=1/\big{(}\sum_{\alpha}x^{\alpha}\big{)}^{1/d}$ is a solution to the homogenized system of equations.

Next, we add a multiple of the last equation to each of the others in order to make all the coefficients positive. For each $1\leq i\leq\ell$, choose a positive $b_{i}>-\min_{\alpha}\{a_{i\alpha}\mid\alpha\in S^{\prime}\}$, and define $a_{i\alpha}^{\prime}=a_{i\alpha}+b_{i}$. By construction, the resulting system has all positive coefficients, and since the equations are formed from the previous equations by elementary linear transformations, the set of solutions are the same. ∎

We use a variation on the technique used to prove Proposition 9.

First, we pick $2m-2$ generic homogeneous linear functions $b_{1},\ldots,b_{2m-2}$ on $m$ variables. By generic, we mean for any $m$ of the $b_{k}$, the only simultaneous solution of all $m$ linear equations is the trivial one. This genericity implies that any $m-1$ of the $b_{k}$ define a point in $\mathbb{P}^{m-1}$ By taking a linear changes of coordinates in each set of variables, we can assume that all of these points are positive, i.e. have a representative consisting of all positive real numbers.

Then we consider the system of equations

(15)		$\displaystyle b_{k}(x_{1},\ldots,x_{m})\cdot b_{k}(x_{m+1},\ldots,x_{n})$	$\displaystyle=0,\mbox{ for }1\leq k\leq 2m-2$
(16)		$\displaystyle(x_{1}+\ldots+x_{m})(x_{m+1}+\ldots+x_{2m})$	$\displaystyle=1$
(17)		$\displaystyle x_{1}+\ldots+x_{m}$	$\displaystyle=1.$

The equations (15) are bihomogeneous and so we can think of their solutions in $\mathbb{P}^{m-1}\times\mathbb{P}^{m-1}$. There are exactly ${2m-2\choose m-1}$ positive real solutions, corresponding to the subsets $A\subset[2m-2]$ of size $m-1$. For any such $A$, there is a unique, distinct solution satisfying $b_{k}(x_{1},\ldots,x_{m})=0$ for all $k$ in $A$ and $b_{k}(x_{m+1},\ldots,x_{2m})=0$ for all $k$ not in $A$. By assumption, for each solution, all the coordinates can be chosen to be positive. The last two equations (16) and (17) dehomogenize the system in a way such that there are ${2m-2\choose m-1}$ positive real solutions. Finally, as in the last paragraph of the proof of Proposition 9, we can add multiples of (16) to the equations (15) in order to make all the coefficients positive. ∎