Python « Rod Carvalho

Archive for the ‘Python’ Category

Fair play

March 31, 2012

Last February’s Ponder This challenge was the following:

To decide who would play first in a human vs. computer game, the human suggested to flip a coin it produced from its pocket, stating “heads you win, tails I win.” The computer replied, “Humans are very devious. I know that your coin is not fair, and that the probability of heads on that coin is slightly less than one half. To make this perfectly fair, flip the coin repeatedly, and I win if the coin first shows heads on toss number 1, 14, 15, 19, 20, or 23. You win if it first shows heads on any other toss.”

If the computer is correct and $p$ is the probability that the coin shows heads, give the best rational approximation to $p$ with a denominator having fewer than 10 digits.

Here is the solution.

__________

My solution #1:

Let $\mathbb{P} (H) = p$ be the probability that the coin will land with the heads side up. If the coin is biased in favor of the human player, then $p < 1/2$ . The probability that the coin toss will first produce “heads” on the $n$ -th trial (where $n \geq 1$ ) is $(1-p)^{n-1} p$ . Let us define $\mathcal{I} := \{1, 14, 15, 19, 20, 23\}$ . The probability that the computer will win is thus

$\displaystyle \sum_{n \in \mathcal{I}} (1-p)^{n-1} p$ .

We now assume that the computer knows the exact value of $p$ and, thus, is able to compute its probability of winning. Since the computer is interested in making this game “perfectly fair”, then his probability of winning should equal the human’s probability of winning, i.e., it should equal $1/2$ . Hence, we obtain the equation

$p + (1-p)^{13} p + (1-p)^{14} p + (1-p)^{18} p + (1-p)^{19} p +(1-p)^{22} p = \frac{1}{2}$

which we can solve using the following Python + SymPy script:

from sympy import *

# create symbolic variable
p = Symbol('p')

# solve polyomial equation
f = p + p*(1-p)**13+p*(1-p)**14+p*(1-p)**18+p*(1-p)**19+p*(1-p)**22
Solutions = solve(f-0.5,p)

# find (purely) real solutions
RealSols = filter(lambda z: im(z)==0, Solutions)

# print solutions
print "\nAll solutions:"
for sol in Solutions:
    print sol

# print real solutions
print "\nReal solutions:"
for sol in RealSols:
	print sol

This script outputs the following list of solutions:

All solutions:
0.334938061741437 - 1.03291350457235*I
0.334938061741437 + 1.03291350457235*I
1.97217279050061 + 0.143849563303172*I
1.50662666735084 - 0.797180194029411*I
1.92419089868505 + 0.454390584787364*I
0.0780115260218319 + 0.464383156442821*I
0.00197767511752256 - 0.101887310197153*I
1.27721993255011 + 0.951650156864179*I
1.27721993255184 - 0.951650156863701*I
0.247036637335883 - 0.624592516289942*I
0.625456509765409 - 0.845921918004056*I
0.0780115260218319 - 0.464383156442821*I
0.00197767511752256 + 0.101887310197153*I
0.991590970537258 - 0.923668626671437*I
1.79082570892033 - 0.599036471538201*I
1.92419089870934 - 0.454390584808249*I
1.79082570891289 + 0.599036471546724*I
0.499905242928390
0.6254565097654 + 0.845921918004058*I
0.991590970537348 + 0.923668626671824*I
0.247036637335883 + 0.624592516289942*I
1.50662666734992 + 0.797180194030292*I
1.97217279048996 - 0.14384956326771*I

Real solutions:
0.499905242928390

Hence, we have that the probability of heads is

$p = 0.499905242928390$ .

We can convert this number to a ratio of integers in Haskell using function approxRational in the library Data.Ratio:

import Data.Ratio

-- probability of heads
p = 0.499905242928390

-- create (infinite) list of epsilons
epsilons = map (10**) [-1,-2..]

-- create (infinite) list of ratios
ratios = map (approxRational p) epsilons

-- list of ratios whose denominators have less than 10 digits
ratios' = takeWhile (\x -> denominator x < 10^10) ratios

Running the script above in GHCi, we obtain a list of ratios:

*Main Data.Ratio> ratios'
[1 % 2,1 % 2,1 % 2,1 % 2,2386 % 4773,2611 % 5223,2636 % 5273,
2638 % 5277,13189 % 26383,44843 % 89703,166183 % 332429,
393036 % 786221,1405961 % 2812455,6243733 % 12489833,
14906352 % 29818355,69693988 % 139414397]

and we finally obtain the following rational approximation

$p \approx \displaystyle\frac{69693988}{139414397}$

whose denominator has nine digits. According to the webmaster, the rational approximation I presented above is not the best one. Hence, let us now try to obtain a better one.

__________

My solution #2:

Most likely, the reason why I did not obtain the best rational approximation is that the real solution I obtained was not precise enough. Hence, let us use the following MATLAB command to solve the polynomial equation:

solve('p + p*(1-p)^13+p*(1-p)^14+p*(1-p)^18+p*(1-p)^19+p*(1-p)^22 = 0.5')

This command produces the following list of solutions:

ans =

 .19776751175225606651411342151624e-2-.10188731019715301414082803578298*i
 .19776751175225606651411342151624e-2+.10188731019715301414082803578298*i
 .78011526021831885333432975572559e-1-.46438315644282071005250209738747*i
 .78011526021831885333432975572559e-1+.46438315644282071005250209738747*i
    .24703663733588276238861533742100-.62459251628994218281562354280746*i
    .24703663733588276238861533742100+.62459251628994218281562354280746*i
    .33493806174143698737516495548579-1.0329135045723528674007841220342*i
    .33493806174143698737516495548579+1.0329135045723528674007841220342*i
                                        .49990524292838950452962194522516
    .62545650976540488339104798827398-.84592191800405765193777444233924*i
    .62545650976540488339104798827398+.84592191800405765193777444233924*i
    .99159097053730697830960408610308-.92366862667162608535284662817351*i
    .99159097053730697830960408610308+.92366862667162608535284662817351*i
    1.2772199325509563050241902871670-.95165015686389448554291677056666*i
    1.2772199325509563050241902871670+.95165015686389448554291677056666*i
    1.5066266673486985232104204448828-.79718019402868789409470957216463*i
    1.5066266673486985232104204448828+.79718019402868789409470957216463*i
    1.7908257089204728277307542114304-.59903647154154201043833983995659*i
    1.7908257089204728277307542114304+.59903647154154201043833983995659*i
    1.9241908987061350990877558129099-.45439058479360516453156947604463*i
    1.9241908987061350990877558129099+.45439058479360516453156947604463*i
    1.9721727904901564352190617939256-.14384956327517564484908725760533*i
    1.9721727904901564352190617939256+.14384956327517564484908725760533*i

We now obtain a much more precise real solution

$p =0.49990524292838950452962194522516$ .

We modify the Haskell script, using the new value of $p$ and a much finer list of epsilons (the precisions):

import Data.Ratio

-- probability of heads
p = 0.49990524292838950452962194522516

-- create (infinite) list of epsilons
epsilons = map ((1.005)**) [-1,-2..]

-- create (infinite) list of ratios
ratios = map (approxRational p) epsilons

-- list of ratios whose denominators have less than 10 digits
ratios' = takeWhile (\x -> denominator x < 10^10) ratios

Which yields the following rational approximation

$p \approx \displaystyle\frac{61031369}{122085875}$

whose denominator has nine digits. Although this rational approximation does not match the one in the official solution, I still got my name on the “people who answered correctly” list. Surprisingly, the value of $p$ which I obtained using MATLAB is quite different from the value of $p$ used by the webmaster

$p = 0.499905242928506776611191007553$ .

Which value of $p$ is “more correct”? I do not know. Do you?

Tags:Data.Ratio, Haskell, Ponder This, Probability Theory, Puzzles, Python, SymPy
Posted in Haskell, Puzzles, Python | Leave a Comment »

Wythoff sequences

August 21, 2011

On November 20, 2007, the Putnam problem of the day was:

The set of all positive integers is the union of two disjoint subsets $\{f(1), f(2), f(3), \dots\}$ , $\{g(1), g(2), g(3), \dots\}$ , where $f(1) < f(2) < f(3) < \dots$ , and $g(1) < g(2) < g(3) < \dots$ , and $g(n) = f(f(n))+1$ for $n=1, 2, 3,\dots$ . Determine $f(240)$ .

Back in November 2007 I could not solve the problem. Now, almost four years later, I will give it another try.

__________

My “solution”:

Let $\mathbb{N} := \{1, 2, 3, \dots\}$ be the set of natural numbers (positive integers). We define $F := \{f(n)\}_{n \in \mathbb{N}}$ and $G := \{g(n)\}_{n \in \mathbb{N}}$ . Note that we must have $F \cup G = \mathbb{N}$ and $F \cap G = \emptyset$ .

Let $f(1) = 1$ . Then, we obtain

$g(1) = f(f(1)) + 1 = f(1) + 1 = 2$ .

Next, let $f(2) = 3$ and $f(3) = 4$ . Then,

$g(2) = f(f(2)) + 1 = f(3) + 1 = 5$ .

Let $f(4) = 6$ . Then, we obtain

$g(3) = f(f(3)) + 1 = f(4) + 1 = 7$ .

Next, let $f(5) = 8$ and $f(6) = 9$ . Then,

$g(4) = f(f(4)) + 1 = f(6) + 1 = 10$ .

Next, let $f(7) = 11$ and $f(8) = 12$ . Then,

$g(5) = f(f(5)) + 1 = f(8) + 1 = 13$ .

Let $f(9) = 14$ . Then, we obtain

$g(6) = f(f(6)) + 1 = f(9) + 1 = 15$ .

Next, let $f(10) = 16$ and $f(11) = 17$ . Then,

$g(7) = f(f(7)) + 1 = f(11) + 1 = 18$ .

Let $f(12) = 19$ . Then, we obtain

$g(8) = f(f(8)) + 1 = f(12) + 1 = 20$ .

Next, let $f(13) = 21$ and $f(14) = 22$ . Then,

$g(9) = f(f(9)) + 1 = f(14) + 1 = 23$ .

Let us stop iterating here. So far, we have

$\left(f(n)\right)_{n \in \mathbb{N}} = \left(1, 3, 4, 6, 8, 9, 11, 12, 14, 16, 17, 19, 21, 22, \dots\right)$

$\left(g(n)\right)_{n \in \mathbb{N}} = \left(2, 5, 7, 10, 13, 15, 18, 20, 23, \dots\right)$ .

Consulting the On-Line Encyclopedia of Integer Sequences (OEIS), we find out that $\left(f(n)\right)_{n \in \mathbb{N}}$ and $\left(g(n)\right)_{n \in \mathbb{N}}$ are the lower Wythoff (A000201) and upper Wythoff (A001950) sequences, respectively. These two are Beatty sequences generated by the golden ratio $\varphi := \frac{1 + \sqrt{5}}{2}$ , as follows

$f (n) = \lfloor n \varphi \rfloor$ , $g (n) = \lfloor n \varphi^2 \rfloor$

where $\lfloor \cdot \rfloor$ is the floor function. Recall that $\varphi^2 = \varphi + 1$ , which allows us to write $g (n) = \lfloor n (\varphi + 1) \rfloor$ .

The following Python script produces the first twenty terms of the lower and upper Wythoff sequences:

from math import floor
from math import sqrt

# golden ratio
phi = (1 + sqrt(5)) / 2.0

# define function f
def f (n):
    return int(floor(phi * n))

# define function g
def g (n):
    return int(floor((phi+1) * n))

print "Lower and upper Wythoff sequences:\n"
for n in range(1,21):
    print "f(%d) = %d \t g(%d) = %d" % (n, f(n), n, g(n))

This script produces the output:

Lower and upper Wythoff sequences:

f(1) = 1         g(1) = 2
f(2) = 3         g(2) = 5
f(3) = 4         g(3) = 7
f(4) = 6         g(4) = 10
f(5) = 8         g(5) = 13
f(6) = 9         g(6) = 15
f(7) = 11        g(7) = 18
f(8) = 12        g(8) = 20
f(9) = 14        g(9) = 23
f(10) = 16       g(10) = 26
f(11) = 17       g(11) = 28
f(12) = 19       g(12) = 31
f(13) = 21       g(13) = 34
f(14) = 22       g(14) = 36
f(15) = 24       g(15) = 39
f(16) = 25       g(16) = 41
f(17) = 27       g(17) = 44
f(18) = 29       g(18) = 47
f(19) = 30       g(19) = 49
f(20) = 32       g(20) = 52

Finally, we have $f (240) = \lfloor 240 \varphi \rfloor = 388$ .

__________

I am not happy with this “solution”. I obviously cheated when I consulted the OEIS. Supposing one is not acquainted with Beatty sequences, is there any hope of solving the problem? If one happens to be acquainted with Beatty sequences and postulates that

$f (n) = \lfloor n r \rfloor$ , $g (n) = \lfloor n s \rfloor$

where $\frac{1}{r} + \frac{1}{s} = 1$ , how does one arrive at the conclusion that $r = \varphi$ and $s = \varphi^2$ ? In other words, what would a proper solution look like?

Tags:Beatty Sequences, Integer Sequences, Putnam problem of the day, Python, Wythoff Sequences
Posted in Puzzles, Python | 1 Comment »

Sum of squares decomposition II

July 5, 2011

Last week, we concluded that the quartic form

$F(x, y) = 2 x^4 + 2 x^3 y - x^2 y^2 + 5 y^4$

has the following sum of squares decomposition

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 -3 y^2 + x y\right)^2 + \frac{1}{2} \left(y^2 + 3 x y\right)^2$

and, thus, $F$ is globally nonnegative. More generally, we saw that $F$ can be written as follows

$F(x,y) = \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]^T \left[\begin{array}{ccc} 2 & -\lambda & 1\\ -\lambda & 5 & 0\\ 1 & 0 & 2 \lambda - 1\end{array}\right] \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]$

where $\lambda \in \mathbb{R}$ is a parameter. Let us define $z := (x^2, y^2, x y)$ and

$Q := \left[\begin{array}{ccc} 2 & -\lambda & 1\\ -\lambda & 5 & 0\\ 1 & 0 & 2 \lambda - 1\end{array}\right]$

so that we can write polynomial $F$ in the more compressed form $F(x, y) = z^T Q z$ . A sufficient condition for the existence of a SOS decomposition is that $Q$ is positive semidefinite, i.e., $Q \succeq 0$ . Why? If $Q$ is positive semidefinite, then it has a Cholesky factorization $Q = L^T L$ and, thus, $F(x, y) = z^T Q z = z^T L^T L z = \| L z \|_2^2$ , which finally yields the SOS decomposition

$F(x, y) = \displaystyle \sum_{i =1}^3 \left(L_{i1} x^2 + L_{i2} y^2 + L_{i3} xy \right)^2$

where $L_{ij}$ is the $(i,j)$ entry of matrix $L$ . One can also think of the Cholesky factorization as $Q = B B^T$ , where $B = L^T$ , so that we obtain $F(x, y) = \| B^T z \|_2^2$ . However, let us not lose sight of where we are going. The question to be answered is: how do we find a $\lambda \in \mathbb{R}$ such that $Q \succeq 0$ ?

__________

Admissible values of $\lambda$

Saying that $Q$ is positive semidefinite is equivalent to saying that the $2^3 - 1 = 7$ principal minors of $Q$ are nonnegative. If $Q \succeq 0$ , then the three 1st degree principal minors (which are the diagonal entries of $Q$ ) lead to the three inequalities

$2 \geq 0$ , $5 \geq 0$ , $2 \lambda - 1 \geq 0$ .

The first two inequalities are redundant, and the third yields $\lambda \geq \frac{1}{2}$ . The three 2nd degree principal minors (determinants of $2 \times 2$ submatrices of $Q$ ) and the 3rd degree principal minor (determinant of $Q$ ) can be computed with the following Python / SymPy script:


from sympy import *

# symbolic variable
Lambda = Symbol('Lambda')

# build matrix Q
Q = Matrix([[2,-Lambda,1], [-Lambda,5,0], [1,0,2*Lambda-1]])

print "\n2nd degree principal minors:"
print Q.extract([0,1],[0,1]).det()
print Q.extract([0,2],[0,2]).det()
print Q.extract([1,2],[1,2]).det()

print "\n3rd degree principal minor:"
print Q.det().factor()

which produces the following output:

2nd degree principal minors:
-Lambda**2 + 10
4*Lambda - 3
10*Lambda - 5

3rd degree principal minor:
-(Lambda - 3)*(2*Lambda**2 + 5*Lambda - 5)

From the three 2nd degree principal minors, we then obtain the three inequalities

$- \lambda^2 + 10 \geq 0$ , $4 \lambda - 3 \geq 0$ , $10 \lambda - 5 \geq 0$

which are equivalent to

$\lambda^2 \leq 10$ , $\lambda \geq \frac{3}{4}$ , $\lambda \geq \frac{1}{2}$

which yields $\lambda \in [\frac{3}{4}, \sqrt{10}]$ . From the 3rd degree principal minor, we obtain the inequality

$- 2 \lambda^{3} + \lambda^{2} + 20 \lambda - 15 = - \left(\lambda -3\right) \left(2 \lambda^{2} + 5 \lambda -5\right) \geq 0$ .

Hence, we have that admissible values of $\lambda$ must satisfy

$\lambda \in [\frac{3}{4}, \sqrt{10}]$ , $\left(\lambda -3\right) \left(2 \lambda^{2} + 5 \lambda -5\right) \leq 0$ .

The roots of polynomial $2 \lambda^{2} + 5 \lambda -5$ are $(-5 \pm \sqrt{65}) / 4$ . Unfortunately $(-5 + \sqrt{65}) / 4$ is in $[\frac{3}{4}, \sqrt{10}]$ . Hence, let us partition the interval $[\frac{3}{4}, \sqrt{10}]$ into two intervals:

for $\lambda \in [\frac{3}{4}, \frac{-5 + \sqrt{65}}{4})$ , we have that $2 \lambda^{2} + 5 \lambda -5$ is negative, and thus $\left(\lambda -3\right) \left(2 \lambda^{2} + 5 \lambda -5\right) \leq 0$ implies that $\lambda \geq 3$ , but

$\{ \lambda \in \mathbb{R} : \lambda \in [\frac{3}{4}, \frac{-5 + \sqrt{65}}{4}) \land \lambda \geq 3\} = \emptyset$ .

for $\lambda \in [\frac{-5 + \sqrt{65}}{4}, \sqrt{10}]$ , we have that polynomial $2 \lambda^{2} + 5 \lambda -5$ is nonnegative, and thus $\left(\lambda -3\right) \left(2 \lambda^{2} + 5 \lambda -5\right) \leq 0$ implies that $\lambda \leq 3$ . Hence,

$\{ \lambda \in \mathbb{R} : \lambda \in [\frac{-5 + \sqrt{65}}{4}, \sqrt{10}] \land \lambda \leq 3\} = [\frac{-5 + \sqrt{65}}{4}, 3]$ .

Finally, we conclude that the set of admissible values of $\lambda$ is $\Lambda := [\frac{-5 + \sqrt{65}}{4}, 3]$ . Note that $\frac{-5 + \sqrt{65}}{4} \approx 0.765$ .

__________

An ensemble of SOS decompositions

Now that we know for which values of $\lambda$ matrix $Q$ is positive semidefinite, we can compute various SOS decompositions. The following Python / SymPy script computes an SOS decomposition for $\lambda = 3$ :


from sympy import *

# symbolic variables
Lambda = Symbol('Lambda')
x = Symbol('x')
y = Symbol('y')

# define polynomial
F = 2*(x**4) + 2*(x**3)*y - (x**2)*(y**2) + 5*(y**4)

# build matrix Q
Q = Matrix([[2,-Lambda,1], [-Lambda,5,0], [1,0,2*Lambda-1]])

# build vector of monomials
z = Matrix(3, 1, [x**2, y**2, x*y])

# compute Cholesky factorization
# for a particular value of Lambda
B = Q.subs(Lambda, Rational(3,1)).cholesky()

# obtain vector of SOS polynomials
f = B.T * z

print "\nSOS polynomials:"
pprint(f[0])
pprint(f[1])
pprint(f[2])

# compute residual polynomial
R = F - sum([f[i]**2 for i in [0,1,2]])

print "\nResidual polynomial: %s" % R.expand()

and it produces the SOS decomposition

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 -3 y^2 + x y\right)^2 + \frac{1}{2} \left(y^2 + 3 x y\right)^2$

which is the one we obtained last week. For $\lambda = 3/2$ , we get

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 - \frac{3}{2} y^2 + x y\right)^2 + 62 \left(\frac{1}{4} y^2 + \frac{3}{62} x y\right)^2 +\frac{1302}{961} x^2 y^2$

and to verify that this is correct, we can run the script:

h0 = Rational(1,2) * (2*x**2 - Rational(3,2)*y**2 + x*y)**2
h1 = 62 * (Rational(3,62)*x*y + Rational(1,4)*y**2)**2
h2 = Rational(1302,961)*x**2*y**2
pprint((F - (h0 + h1 + h2)).expand())

For $\lambda = 5/2$ , we obtain the SOS decomposition

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 - \frac{5}{2} y^2 + x y\right)^2 + \frac{30}{16} \left(y^2 + \frac{2}{3} x y\right)^2 +\frac{8}{3} x^2 y^2$ .

But, Доверяй, но проверяй!

h0 = Rational(1,2) * (2*x**2 - Rational(5,2)*y**2 + x*y)**2
h1 = Rational(30,16) * (Rational(2,3)*x*y + y**2)**2
h2 = Rational(8,3)*x**2*y**2
pprint((F - (h0 + h1 + h2)).expand())

__________

What happens for non-admissible values of $\lambda$ ?

We know that if $\lambda \in \Lambda$ , then $Q \succeq 0$ . What if $\lambda \notin \Lambda$ ? Then matrix $Q$ will be negative definite and will, or will not, have a Cholesky factorization, depending on what convention one uses.

One can think of matrices $L$ and $B$ in the Cholesky factorizations $Q = L^T L$ and $Q = B B^T$ , respectively, as square roots of matrix $Q$ . Do negative reals have square roots? They do if you are willing to work with complex numbers. Likewise, a negative definite $Q$ will have a Cholesky factorization if one is willing to work with complex matrices. For example, for $\lambda = 0$ we obtain the SOS decomposition

$F(x, y) = \displaystyle \sum_{i=0}^2 f_i^2 (x, y)$

where the $f_i$ polynomials are:

$f_0 (x,y) = \frac{\sqrt{2}}{2} (2 x^2 + x y)$ , $f_1 (x,y) = \sqrt{5} y^2$ , $f_2 (x,y) = i \frac{\sqrt{6}}{2} x y$

where $i := \sqrt{-1}$ . Note that $f_0, f_1 \in \mathbb{R}[x, y]$ , but $f_2 \in \mathbb{C}[x, y]$ . If we square $f_2$ , we get $f_2^2 = -\frac{3}{2} x^2 y^2$ . This term is being subtracted from the sum of squares, which prevents us from concluding global nonnegativity. That is why we impose the condition $Q \succeq 0$ .

Tags:Algebra, Global Nonnegativity, Polynomials, Real Algebraic Geometry, Sum of Squares, Sum of Squares Decomposition, SymPy
Posted in Algebra, Algebraic Geometry, Python | Leave a Comment »

Sum of squares decomposition

June 30, 2011

Consider the following quartic form in $\mathbb{R}[x, y]$

$F(x, y) = 2 x^4 + 2 x^3 y - x^2 y^2 + 5 y^4$

and suppose we would like to decide if $F$ is globally nonnegative. In other words, we want to know whether the proposition

$F(x, y) \geq 0, \qquad{} \forall x, y \in \mathbb{R}$

is true. A sufficient condition for $F(x,y)$ to be globally nonnegative is the existence of a sum of squares decomposition [1, 2]

$F(x, y) = \displaystyle \sum_{i} f_i^2 (x, y)$

where $f_i \in \mathbb{R}[x, y]$ . If such a decomposition does indeed exist, how do we compute it? How do we find the $f_i$ polynomials?

Note that $F(x, y)$ is a quartic form, i.e., a polynomial whose monomials all have total degree 4. Using Parrilo’s method [1, 2], let us try to write $F(x, y)$ as a quadratic form in all the monomials of total degree equal to 2, as follows

$F(x,y) = \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]^T \left[\begin{array}{ccc} q_{11} & q_{12} & q_{13}\\ q_{12} & q_{22} & q_{23}\\ q_{13} & q_{23} & q_{33}\end{array}\right] \left[\begin{array}{c} x^2\\ y^2\\ x y\end{array}\right]$

where the real matrix

$Q := \left[\begin{array}{ccc} q_{11} & q_{12} & q_{13}\\ q_{12} & q_{22} & q_{23}\\ q_{13} & q_{23} & q_{33}\end{array}\right]$

is symmetric. Hence, we obtain

$F(x,y) = q_{11} x^4 + 2 q_{13} x^3 y + (2 q_{12} + q_{33}) x^2 y^2 + 2 q_{23} x y^3 + q_{22} y^4$

and, matching coefficients, we get the equality constraints

$q_{11} = 2$ , $q_{13} = 1$ , $2 q_{12} + q_{33} = -1$ , $q_{23} = 0$ , $q_{22} = 5$ .

Let us make $q_{12} = - \lambda$ , where $\lambda \in \mathbb{R}$ , so that $q_{33} = 2 \lambda - 1$ . Hence, we obtain a family of symmetric $Q$ matrices parameterized by $\lambda$

$Q := \left[\begin{array}{ccc} 2 & -\lambda & 1\\ -\lambda & 5 & 0\\ 1 & 0 & 2 \lambda - 1\end{array}\right]$

For each choice of $\lambda$ , we get a matrix $Q$ . We then conclude that the representation $F(x, y) = z^T Q z$ , where $z := (x^2, y^2, x y)$ , is not unique. Since there are infinitely many admissible $Q$ matrices, if we are able to find a positive semidefinite one, then $F(x, y) = z^T Q z \geq 0$ for all $z$ , and global nonnegativity of the given polynomial $F$ is guaranteed!

For instance, let us try $\lambda = 3$ . We then get

$Q = \left[\begin{array}{ccc} 2 & -3 & 1\\ -3 & 5 & 0\\ 1 & 0 & 5\end{array}\right]$

which is positive semidefinite. Computing a Cholesky factorization of $Q$ , i.e., finding a matrix $L$ such that $Q = L^T L$ , we get

$L = \displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{ccc} 2 & -3 & 1\\ 0 & 1 & 3\end{array}\right]$

and, thus, $F(x, y) = z^T Q z = z^T L^T L z = \| L z \|_2^2$ , which finally yields a sum of squares (SOS) decomposition

$F(x, y) = \displaystyle \frac{1}{2} \left(2 x^2 -3 y^2 + x y\right)^2 + \frac{1}{2} \left(y^2 + 3 x y\right)^2$ .

To make sure that this SOS decomposition is correct, we can use a little Python / SymPy script:


from sympy import *

# symbolic variables
x = Symbol('x')
y = Symbol('y')

# define polynomial
F = 2*(x**4) + 2*(x**3)*y - (x**2)*(y**2) + 5*(y**4)

# define SOS polynomials
f1 = (sqrt(2)/2) * (2*(x**2) - 3*(y**2) + x*y)
f2 = (sqrt(2)/2) * (y**2 + 3*x*y)

# construct residual polynomial
R = F - (f1**2 + f2**2)

print "Residual polynomial: %s" % R.expand()

The residual polynomial is $R (x, y) = 0$ . The SOS decomposition is correct! We can then conclude that the given polynomial $F$ is globally nonnegative.

__________

Remark: this post is based on example 4.1 of Parrilo’s doctoral dissertation [1] and also on example 3.5 of Parrilo’s paper [2].

__________

References

[1] Pablo Parrilo, Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization, Dissertation (Ph.D.), California Institute of Technology, 2000.

[2] Pablo Parrilo, Semidefinite programming relaxations for semialgebraic problems, Mathematical Programming, 2003.

Tags:Algebra, Global Nonnegativity, Polynomials, Real Algebraic Geometry, Sum of Squares, Sum of Squares Decomposition, SymPy
Posted in Algebra, Algebraic Geometry, Python | Leave a Comment »

Fixed-Odds Betting Arbitrage II

December 28, 2010

My recent post on fixed-odds betting arbitrage left some rather intriguing questions unanswered. Although I used Linear Programming to find an arbitrage opportunity, the focus was on constraint satisfaction, rather than optimality. But, then, if betting can be viewed as playing a sequential game with fate, what would the word “optimal” even mean in this context?

In this post, I will focus on optimality. I will assume the reader is already acquainted with my previous post and, therefore, for convenience I will abstain from “beginning at the beginning”. In other words, this post is not self-contained.

__________

Arbitrage

Given the odds matrix $\Omega$ , the set of admissible $\theta$ vectors $\Theta \subseteq \mathbb{F}_2^m$ , and a budget $c > 0$ , we would like to find a money allocation matrix $X$ such that $1_{m n}^T \text{vec}(X) \leq c$ and

$P (\theta, X; \Omega) \geq 0, \quad{} \forall \theta \in \Theta$

so that, regardless of what $\theta$ vector fate happens to choose, we make a profit. Let us enumerate the $\theta$ vectors, $\Theta = \{\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(|\Theta|)}\}$ , where $m \leq |\Theta|\leq 2^m$ . Since $P (\theta, X; \Omega) = \text{vec}^T \left( R (\theta; \Omega) \right) \text{vec}(X)$ , the $|\Theta|$ arbitrage inequalities can be written as

$\left[\begin{array}{c} \text{vec}^T \left( R (\theta^{(1)}; \Omega) \right) \\ \text{vec}^T \left( R (\theta^{(2)}; \Omega) \right)\\ \vdots \\ \text{vec}^T \left( R (\theta^{|\Theta|}; \Omega) \right)\end{array}\right] \text{vec}(X) \geq \left[\begin{array}{c} 0 \\ 0 \\ \vdots \\ 0\end{array}\right]$ .

Let $\tilde{x} := \text{vec}(X)$ . Hence, the inequality $1_{m n}^T \text{vec}(X) \leq c$ can be rewritten as $1_{m n}^T \tilde{x} \leq c$ . Let us define $q_k := \text{vec} \left( R (\theta^{(k)}; \Omega) \right)$ , and let

$Q := \left[\begin{array}{cccc} q_1 & q_2 & \ldots & q_{|\Theta|}\end{array}\right]$

be the $m n \times |\Theta|$ matrix whose $k$ -th column is $q_k$ . Thus, the $|\Theta|$ arbitrage inequalities can be written more compactly as $Q^T \tilde{x} \geq 0_{|\Theta|}$ or, equivalently, as $-Q^T \tilde{x} \leq 0_{|\Theta|}$ . Hence, so far, we have $1+|\Theta|$ inequality constraints

$\left[\begin{array}{c} 1_{mn}^T\\ -Q^T\end{array}\right] \tilde{x} \leq \left[\begin{array}{c} c \\ 0_{|\Theta|}\end{array}\right]$ .

Since the $m n$ entries of matrix $X$ must be non-negative, we have $m n$ additional inequality constraints $\tilde{x} \geq 0_{m n}$ or, equivalently, $- I_{m n} \tilde{x} \leq 0_{m n}$ , where $I_{m n}$ is the $m n \times m n$ identity matrix.

To make the model more realistic, let us suppose that the amount of money being betted on outcome $i \in [m]$ at bookmaker $j \in [n]$ is bounded above by $b_{ij} \geq 0$ , i.e., $0 \leq x_{ij} \leq b_{ij}$ . This is useful not only because bookmakers impose upper bounds on how much money one can bet, but also because it may happen that not all bookmakers are posting odds on the very same outcomes. By making $b_{ij} = 0$ one can make sure that no money will be betted on a certain outcome at a certain bookmaker. If bookmaker $j$ is not imposing any upper bound on $x_{ij}$ then we can make $b_{ij} = c$ , i.e., the most we can bet on outcome $i$ at bookmaker $j$ is our budget $c > 0$ . Let

$B = \left[\begin{array}{cccc}b_{11} & b_{12} & \ldots & b_{1n}\\b_{21} & b_{22} & \ldots & b_{2n}\\\vdots & \vdots & \ddots & \vdots\\b_{m1} & b_{m2} & \ldots & b_{mn}\\\end{array}\right]$

be the matrix of upper bounds, and let $\tilde{b} := \text{vec}\left(B\right)$ . Then, the $m n$ inequalities $x_{ij} \leq b_{ij}$ can be written more compactly as $\tilde{x} \leq \tilde{b}$ or, equivalently, as $I_{m n} \tilde{x} \leq \tilde{b}$ .

Finally, we have a total of $1+|\Theta| + 2 m n$ inequality constraints

$\left[\begin{array}{c} 1_{mn}^T\\ -Q^T\\ - I_{m n}\\ I_{m n}\end{array}\right] \tilde{x} \leq \left[\begin{array}{c} c \\ 0_{|\Theta|} \\ 0_{m n} \\ \tilde{b}\end{array}\right]$ .

If we can find an $m n$ -dimensional decision vector $\tilde{x}$ that satisfies the $1+|\Theta| + 2 m n$ inequality constraints above, then we conclude that an arbitrage opportunity does exist. This is a constraint satisfaction problem.

As we discussed before, the $1+|\Theta| + 2 m n$ inequality constraints define a polytope in $\mathbb{R}^{m n}$ . If this polytope is non-empty, then there exists at least one arbitrage opportunity. We can determine if the polytope is non-empty by solving a linear program with an arbitrary objective function

$\begin{array}{ll} \displaystyle\min_{\tilde{x} \in \mathbb{R}^{mn}} & 0_{m n}^T \tilde{x} \\ \text{subject to} & \left[\begin{array}{c} 1_{mn}^T\\ -Q^T\\ - I_{m n}\\ I_{m n}\end{array}\right] \tilde{x} \leq \left[\begin{array}{c} c \\ 0_{|\Theta|} \\ 0_{m n} \\ \tilde{b}\end{array}\right]\end{array}$

where we chose $0_{m n}^T \tilde{x}$ , the most nihilistic of all linear objective functions. If this linear program is feasible, then an arbitrage opportunity does exist. In other words, we solve the constraint satisfaction problem via linear programming.

__________

Minimum Guaranteed Profit

We can go beyond mere arbitrage. Suppose we want to guarantee a minimum return of $\alpha \geq 0$ regardless of what $\theta$ is chosen by fate, i.e.,

$P (\theta, X; \Omega) \geq \alpha c, \quad{} \forall \theta \in \Theta$ .

Hence, instead of the linear inequality constraints $Q^T \tilde{x} \geq 0_{|\Theta|}$ , we now have the linear inequality constraints $Q^T \tilde{x} \geq \alpha c 1_{|\Theta|}$ or, equivalently, $-Q^T \tilde{x} \leq -\alpha c 1_{|\Theta|}$ . Suppose also that the bookmakers impose no upper bounds on the $x_{ij}$ . Hence, we have $b_{ij} = c$ , and $\tilde{b} = c 1_{m n}$ . We then have the following $1+|\Theta| + 2 m n$ inequality constraints

$\left[\begin{array}{c} 1_{mn}^T\\ -Q^T\\ - I_{m n}\\ I_{m n}\end{array}\right] \tilde{x} \leq \left[\begin{array}{c} c \\ -\alpha c 1_{|\Theta|} \\ 0_{m n} \\ c 1_{m n}\end{array}\right]$ .

Yet once again, we have a constraint satisfaction problem that can be solved via linear programming with an arbitrary objective function

$\begin{array}{ll} \displaystyle\min_{\tilde{x} \in \mathbb{R}^{mn}} & 0_{m n}^T \tilde{x} \\ \text{subject to} & \left[\begin{array}{c} 1_{mn}^T\\ -Q^T\\ - I_{m n}\\ I_{m n}\end{array}\right] \tilde{x} \leq \left[\begin{array}{c} c \\ -\alpha c 1_{|\Theta|} \\ 0_{m n} \\ c 1_{m n}\end{array}\right]\end{array}$ .

Starting with a small $\alpha$ , say $\alpha = 0.01$ , and increasing it in a successive fashion until the linear program is no longer feasible, one can obtain money allocation matrices that yield higher and higher minimum guaranteed profits.

__________

Beyond Constraint Satisfaction

Instead of trying various values of $\alpha$ until the linear program that solves the constraint satisfaction problem is no longer feasible, we can find the maximum $\alpha$ via optimization. Note that $Q^T \tilde{x} \geq \alpha c 1_{|\Theta|}$ can be rewritten as

$\alpha c 1_{|\Theta|} - Q^T \tilde{x} \leq 0_{|\Theta|}$ .

Hence, we can find the maximum $\alpha$ that yields a feasible linear program by solving the following maximization problem in $(\alpha, \tilde{x})$ with the additional inequality constraint $\alpha \geq 0$ or, equivalently, $-\alpha \leq 0$

$\begin{array}{ll} \displaystyle\max_{(\alpha, \tilde{x}) \in \mathbb{R}^{1 + mn}} & \left[\begin{array}{cc} 1 & 0_{m n}^T \end{array}\right]\left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \\\\ \text{subject to} & \left[\begin{array}{cc} -1 & 0_{m n}^T\\ 0 & 1_{mn}^T\\ c 1_{|\Theta|} & -Q^T\\ 0_{m n} & - I_{m n}\\ 0_{m n} & I_{m n}\end{array}\right] \left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \leq \left[\begin{array}{c} 0\\ c \\ 0_{|\Theta|} \\ 0_{m n} \\ c 1_{m n}\end{array}\right]\end{array}$

Note that now we do have an actual optimization problem, i.e., the linear program above is a maximization problem, not merely a constraint satisfaction problem. Since maximizing $\alpha$ is the same as minimizing $-\alpha$ , we can rewrite the maximization problem as a minimization problem, as follows

$\begin{array}{ll} \displaystyle\min_{(\alpha, \tilde{x}) \in \mathbb{R}^{1 + mn}} & \left[\begin{array}{cc} -1 & 0_{m n}^T \end{array}\right]\left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \\\\ \text{subject to} & \left[\begin{array}{cc} -1 & 0_{m n}^T\\ 0 & 1_{mn}^T\\ c 1_{|\Theta|} & -Q^T\\ 0_{m n} & - I_{m n}\\ 0_{m n} & I_{m n}\end{array}\right] \left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \leq \left[\begin{array}{c} 0\\ c \\ 0_{|\Theta|} \\ 0_{m n} \\ c 1_{m n}\end{array}\right]\end{array}$

The solutions of the maximization / minimization problems are maximin solutions, i.e., we are maximizing the minimum profit. In other words, we are pessimists and assume that fate will always choose the most damaging $\theta \in \Theta$ . Hence, we choose a money allocation matrix $X$ such that the minimum profit will be maximized. We are looking for the best worst-case scenario.

__________

Example: Sharapova versus Kirilenko

Let us revisit the tennis match we considered previously. Suppose that $n = 2$ bookmakers are posting odds for a Sharapova versus Kirilenko match, and that the bookies are offering odds for $m = 2$ mutually exclusive outcomes: Sharapova wins, or Kirilenko wins. Imagine that the bookies posted the following decimal odds

$\Omega = \left[\begin{array}{cc} 1.25 & 1.43\\ 3.90 & 2.85\end{array}\right]$

which means that bookie #1 is offering odds $1.25$ for a Sharapova victory, whereas bookie #2 is offering odds $1.43$ for the same outcome. Note that $|\Theta| = 2$ and that the $Q$ matrix is

$Q = \left[\begin{array}{cc} 0.25 & -1\\ -1 & 2.90\\ 0.43 & -1\\ -1 & 1.85\end{array}\right]$ .

Assuming that the bookies impose no upper bounds on the $x_{ij}$ , then the maximin return and the maximin money allocation matrix can be found by solving the linear program

$\begin{array}{ll} \displaystyle\min_{(\alpha, \tilde{x}) \in \mathbb{R}^{5}} & \left[\begin{array}{cc} -1 & 0_{4}^T \end{array}\right]\left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \\\\ \text{subject to} & \left[\begin{array}{cc} -1 & 0_{4}^T\\ 0 & 1_{4}^T\\ c 1_{2} & -Q^T\\ 0_{4} & - I_{4}\\ 0_{4} & I_{4}\end{array}\right] \left[\begin{array}{c} \alpha \\ \tilde{x}\end{array}\right] \leq \left[\begin{array}{c} 0\\ c \\ 0_{2} \\ 0_{4} \\ c 1_{4}\end{array}\right]\end{array}$

where $c > 0$ is our budget. The following Python 2.5 / CVXOPT script solves a linear program to determine whether there are arbitrage opportunities and, if so, what the maximin solution is:

 from cvxopt import matrix
 from cvxopt import spmatrix
 from cvxopt import solvers

 # available budget
 c = 100

 # build (decimal) odds matrix
 Omega = matrix([[1.25, 3.90],
                 [1.43, 2.85]])

 # build q vectors and Q matrix
 q1 = matrix([0.25, -1.0, 0.43, -1.0])
 q2 = matrix([-1.0, 2.90, -1.0, 1.85])
 Q = matrix([[q1], [q2]])

 # build 4x4 identity matrix
 Id = spmatrix(1.0, range(4), range(4))

 # ---------------
 # linear program
 # ---------------

 # build objective vector
 f = matrix([-1.0, matrix(0.0, (4,1))])

 # build inequality constraint matrices
 A = matrix([[-1.0, 0.0, c*matrix(1.0, (2,1)), matrix(0.0, (8,1))],
 [matrix(0.0, (1,4)), matrix(1.0, (1,4)), -Q.T, -Id, Id]])
 b = matrix([0.0, c, matrix(0.0, (6,1)), c*matrix(1.0, (4,1))])

 # solve linear program
 solvers.options['show_progress'] = True
 solution = solvers.lp(f, A, b)

 # print solution and profits
 sol = solution['x']
 alpha = sol[0]
 x = sol[1:5]

 # build money allocation matrix
 X = matrix([[sol[1], sol[2]],[sol[3], sol[4]]])

 print "\nMaximin return (%):"
 print 100 * alpha
 print "\nMaximin money allocation matrix:"
 print X
 print "Profit if Sharapova wins:"
 print q1.T * x
 print "Profit if Kirilenko wins:"
 print q2.T * x
 print "Total amount at stake:"
 print sum(x)

The output of this script is the following:

     pcost       dcost         gap    pres   dres   k/t
 0:  1.7194e-002 -5.0054e+002  7e+002 4e-001 2e+002 1e+000
 1:  1.9534e-001 -1.7513e+001  2e+001 2e-002 7e+000 1e+000
 2:  1.3994e-001 -2.3748e+000  2e+000 2e-003 1e+000 2e-001
 3:  1.4821e-002 -8.9421e-001  9e-001 8e-004 4e-001 7e-002
 4: -1.7326e-002 -9.9421e-002  7e-002 7e-005 3e-002 9e-004
 5: -3.8136e-002 -5.1946e-002  1e-002 1e-005 6e-003 3e-004
 6: -4.6244e-002 -4.6469e-002  2e-004 2e-007 9e-005 8e-006
 7: -4.6340e-002 -4.6343e-002  2e-006 2e-009 9e-007 8e-008
 8: -4.6341e-002 -4.6341e-002  2e-008 2e-011 9e-009 8e-010
Optimal solution found.

Maximin return (%):
4.63414537166

Maximin money allocation matrix:
[ 2.31e-006  7.32e+001]
[ 2.68e+001  1.53e-006]

Profit if Sharapova wins:
[ 4.63e+000]

Profit if Kirilenko wins:
[ 4.63e+000]

Total amount at stake:
99.999998456

Since a solution has been found, we conclude that at least one arbitrage opportunity exists. Moreover, the maximin return is $4.63\%$ and, hence, no matter what $\theta$ fate does happen to choose, we have a guaranteed profit of at least $4.63\%$ . Indeed, note that regardless of which Maria wins the tennis match, the profit will be $4.63$ for a budget of $c = 100$ . The money allocation matrix is approximately

$X = \left[\begin{array}{cc} 0.00 & 73.2\\ 26.8 & 0.00\end{array}\right]$ .

As we have rounded off the money allocation matrix produced by the CVXOPT script, the profits will deviate slightly from the maximin ones:

Profit if Sharapova wins:
[ 4.68e+000]

Profit if Kirilenko wins:
[ 4.52e+000]

Total amount at stake:
100.0

Note that the profit in case Sharapova wins is slightly increased, whereas the profit in case Kirilenko wins is slightly decreased.

How “good” is the maximin solution? Let us compare the maximin solution we have just obtained with the mere constraint-satisfying solution we found before:

Profit if Sharapova wins:
[ 2.71e+000]

Profit if Kirilenko wins:
[ 9.14e-001]

Total amount at stake:
80.4618403492

where we have a return of $3.37\%$ in case Sharapova wins, and a return of $1.14\%$ in case Kirilenko wins. It is clear that the maximin solution is superior, not only because the returns are higher, but also because the variance of the returns is considerably lower.

Tags:Arbitrage, Betting, Betting Arbitrage, CVXOPT, Fixed-Odds Betting, Linear Programming, Operations Research
Posted in Operations Research, Python | Leave a Comment »

Rod Carvalho

Archive for the ‘Python’ Category

Fair play

Wythoff sequences

Sum of squares decomposition II

Sum of squares decomposition

Fixed-Odds Betting Arbitrage II

Pages

Categories

Top Posts

Recent Posts

Recent Comments

Tag Cloud

Archives

Blog Stats

Artists

Biology

Blogs

Companies

Haskell

Math Blogs

Mathematics

MATLAB

Optimization

Physics

Python

Robotics

Science

Software

Symbolic Computation

Technology

Powered by FeedBurner

Meta

Follow “Rod Carvalho”

May 2012
S	M	T	W	T	F	S
« Apr
		1	2	3	4	5
6	7	8	9	10	11	12
13	14	15	16	17	18	19
20	21	22	23	24	25	26
27	28	29	30	31

Rod Carvalho

Archive for the ‘Python’ Category

Fair play

Wythoff sequences

Sum of squares decomposition II

Sum of squares decomposition

Fixed-Odds Betting Arbitrage II

Pages

Categories

Top Posts

Recent Posts

Recent Comments

Tag Cloud

Archives

Blog Stats

Artists

Biology

Blogs

Companies

Haskell

Math Blogs

Mathematics

MATLAB

Optimization

Physics

Python

Robotics

Science

Software

Symbolic Computation

Technology

Category Cloud

Powered by FeedBurner

Meta

Follow “Rod Carvalho”