See http://winvector.github.io/freq/minimax.pdf for background and details
From: http://mathoverflow.net/questions/177574/existence-of-solutions-of-a-polynomial-system
Fix \(k \in \mathbb{N}\), \(k \geq 1\). Let \(p \in [0,1]\) and \(x = (x_0, \ldots, x_k)\) be a \((k+1)\)-dimensional real vector, and define \[S_k(p,x) = -x_0^2 + \sum_{i=0}^k {k \choose i} p^i (1 - p)^{k - i} \cdot (x_i - p)^2.\] Experiments show that for small values of \(k\) \[\exists x \in \text{ interior of } [0,1]^{k+1} \,.\, \forall p \in [0,1] \,.\, S_k(p,x) = 0.\] In other words, there are \(x_i\)'s such that \(S_k(x,p)\) is identically zero as a polynomial in \(p\).
For a given \(k\) we can expand \(S_k(x,p)\) as a polynomial in \(p\) and equate the coefficients to \(0\). For \(k = 2\) we getand this has two real solutions: \[x = (\frac{1}{2} (-1-\sqrt{2}),\frac{1}{2},\frac{1}{2} (3+\sqrt{2}))\] and \[x = (\frac{1}{2} (-1+\sqrt{2}),\frac{1}{2},\frac{1}{2} (3-\sqrt{2})).\]
One of which satisifies our conditions.
The problem arises in statistics, see John Mount's blog post for background.
Question: Is there a solution for every \(k\)?
Addendum: John says he wants soltions in the interor of \([0,1]^{k+1}\)...
Solution submitted 8-4-2013 by me (John Mount), having a lot of trouble with links and formatting. Enough of that, submitting it here.
This is some background to the question and the solution (minus one check mentioned at the end).
Define \[S(k,p,x) = \sum_{i=0}^k {k \choose i} p^i (1-p)^{k-i} (x_i-p)^2.\] Define \[f_k(k) = \mathrm{argmin}_x \max_p S(k,p,x).\] Then \(f_k(k)\) is the minimax square-loss solution to trying to estimate the win rate of a random process by observing \(k\) results (Wald wrote on this). The neat thing is: we can show if there is a real solution \(x\) in \([0,1]^{k+1}\) to \(S(k,p,x) = x_0^2\) then \(x=f_k(k)\). Meaning we avoided two nasty quantifiers. See this file for some experimental examples.
We know there is only one connected component of solutions in the interior of the unit cube because these solutions represent extreme points of the minimax estimation problem. We show that there is a diversity of gradients by reflecting coordinates of \(x\) around \(p\) (and thus we have an extreme point).
From the original problem we expect a lot of symmetries. Also, a change of variables \(z = p/(1-p)\) makes collecting terms easier. In fact I now have a conjectured exact solution, I now only need a proof that it always works (cancels the p's, is real and in the interior of \([0,1]^{k+1}\); I already have a proof that such a solution when it exists solves the original estimation problem). The conjectured solution for \(k>1\) (for \(k=1\) the solution is \([1/4,3/4]\)) is:
\[ \begin{align} f_k(0) &= (\sqrt{k}-1)/(2 (k-1)) \\ f_k(1) &= \sqrt{f_k(0)^2+2 f_k(0)/k} \\ & \text{ for } h>1: \\ f_k(h)^2 &= (k+2) (k+1) (f_k(0)^2)/((k+2-h) (k+1-h)) \\ & + 2 h f_k(h-1) (1-f_k(h-1))/(k+1-h) \\ & - h (h-1) ((f_k(h-2)-1)^2)/((k+2-h) (k+1-h)) \end{align} \]
A Python implementation, demonstration and check of this solution up through \(k=8\) is given here. So really all that is left to prove is the right hand side of \(f_k(h)^2\) is always positive and in the interior of \([0,1]\) for all \(k,h\).
Note 8-4-2014: Vladimir Dotsenko finished the solution by adding the important insight that the \(f_k(h)\) are evenly spaced when \(k\) is held constant. This lets him get a closed form solution for each \(f_k(h)\) (without having to refer to ealier \(h\)).
import sympy
# expecting a dictionary solution
def isGoodSoln(si):
def isGoodVal(x):
xn = complex(x)
xr = xn.real
xi = xn.imag
return (abs(xi)<1.0e-6) and (xr>0.0) and (xr<1.0)
return all([ isGoodVal(xi) for xi in si.values() ])
# only good for k>=1
def solveKz(k):
vars = sympy.symbols(['phi' + str(i) for i in range((k+1)/2)])
if k%2!=0:
phis = vars + [1-varsi for varsi in reversed(vars) ]
else:
phis = vars + [sympy.Rational(1,2)] + [1-varsi for varsi in reversed(vars) ]
z = sympy.symbols('z')
poly = sum([ sympy.binomial(k,h) * z**h * ((1+z)*phis[h] -z)**2 for h in range(k+1)]) - phis[0]**2 * (1+z)**(k+2)
polyTerms = poly.expand().collect(z,evaluate=False)
eqns = [ polyTerms[ki] for ki in polyTerms.keys() if (not ki==1) ]
solns = sympy.solve(eqns,vars,dict=True)
soln1 = [ si for si in solns if isGoodSoln(si)][0]
solnv = [ soln1[vi] for vi in vars ]
if k%2!=0:
xs = solnv + [1-solni for solni in reversed(solnv) ]
else:
xs = solnv + [sympy.Rational(1,2)] + [1-solni for solni in reversed(solnv) ]
return xs
# original substitution from inspecting tri-diagonal recurrance
# only good for k>=1
def conjectureKa(k,numeric=False):
if k<=1:
return [sympy.Rational(1,4),sympy.Rational(3,4)]
phi = [ 0 for i in range(k+1) ]
phi[0] = (sympy.sqrt(k)-1)/(2*(k-1))
phi[1] = (sympy.sqrt((phi[0]**2+2*phi[0]/k).expand())).simplify()
if numeric:
for h in range(2):
phi[h] = float(phi[h])
for h in range(2,(k+1)):
phi[h] = sympy.sqrt(( (k+2)*(k+1)*(phi[0]**2)/((k+2-h)*(k+1-h)) + 2*h*phi[h-1]*(1-phi[h-1])/(k+1-h) - h*(h-1)*((phi[h-2]-1)**2)/((k+2-h)*(k+1-h)) ))
return phi
# simplified in pseudo-observation form
def conjectureK(k,numeric=False):
sqrtk = sympy.sqrt(k)
if numeric:
sqrtk = float(sqrtk)
return [(sqrtk/2 + h)/(sqrtk+k) for h in range(k+1) ]
p = sympy.symbols('p')
for k in range(1,9):
print
print 'k',k
solnk = solveKz(k)
print 'soln ',solnk
poly = sum([ p**h * (1-p)**(k-h) * sympy.binomial(k,h) * (solnk[h]-p)**2 for h in range(k+1) ]).expand()
print 'check poly',poly
conjk = conjectureK(k)
print 'conjecture:',conjk
print 'max difference:',max([ abs(complex(solnk[i]-conjk[i])) for i in range(len(solnk)) ])
print '1/k for scale:',1/float(k)
print
p = sympy.symbols('p')
for k in range(1,21):
print
print 'k',k
conjk = conjectureK(k,numeric=True)
print 'conjecture:',conjk
polyc = sum([ p**h * (1-p)**(k-h) * sympy.binomial(k,h) * (conjk[h]-p)**2 for h in range(k+1) ]).expand()
print 'conjecture check poly',polyc
print '1/k for scale:',1/float(k)
print