Archive for August, 2011

Feasibility via quantifier elimination II

August 29, 2011

Suppose we want to decide whether the following semialgebraic set

S := \{(x,y) \in \mathbb{R}^2 \mid x - y^2 + 1 \geq 0 \land -x - y^2 + 1 \geq 0\}

is nonempty. This is the same as determining whether the following system of polynomial inequalities

\begin{array}{rrl} g_1 (x, y) :=& x - y^2 + 1 &\geq 0\\ g_2 (x, y) :=& -x - y^2 +1 &\geq 0\end{array}

is feasible. We now define S_1 := \{(x,y) \in \mathbb{R}^2 \mid g_1 (x, y) \geq 0\}, and S_2 := \{(x,y) \in \mathbb{R}^2 \mid g_2 (x, y) \geq 0\}. Set S_1 is depicted below

and set S_2 is depicted below

Note that S = S_1 \cap S_2. We depict this intersection below

Visual inspection of the plot of the intersection S = S_1 \cap S_2 allows us to conclude that S in nonempty. Asserting that S is nonempty is the same as stating that the proposition

\exists{x}\exists{y} \left(g_1 (x,y) \geq 0 \land g_2 (x,y) \geq 0\right)

is true. We can determine the truth value of this proposition via quantifier elimination using the following REDUCE + REDLOG script:

% feasibility via quantifier elimination

load_package redlog;
rlset ofsf;

% define polynomial functions
g1 :=  x-y^2+1;
g2 := -x-y^2+1;

% define existential formula
phi := ex({x,y}, g1>=0 and g2>=0);

% perform quantifier elimination
rlqe phi;

end;

which produces the truth value true. This means that there exists a (x,y) that satisfies the polynomial inequalities, i.e., S is nonempty.

__________

Remark: the plots in this post were generated by Wolfram Alpha.

Feasibility via quantifier elimination

August 27, 2011

Consider the system of polynomial equations and inequalities [1]

\begin{array}{rrl} f_1 (x_1, x_2) :=& x_1^2 + x_2^2 - 1 &= 0\\ g_1 (x_1, x_2) :=& 3 x_2 - x_1^3 - 2 &\geq 0\\ g_2 (x_1, x_2) :=& x_1 - 8 x_2^3 &\geq 0\end{array}

whose solution set is the following semialgebraic set

S := \{ (x_1,x_2) \in \mathbb{R}^2 \mid f_1 (x_1,x_2) = 0 \land g_1 (x_1,x_2) \geq 0 \land g_2 (x_1,x_2) \geq 0\}

Parrilo showed in [1] that this system is infeasible, i.e., set S is empty, using the Positivstellensatz [2]. Alternatively, one could use quantifier elimination to conclude infeasibility.

Consider the following proposition

\exists{x_1}\exists{x_2} \left(f_1 (x_1,x_2) = 0 \land g_1 (x_1,x_2) \geq 0 \land g_2 (x_1,x_2) \geq 0\right).

If this proposition is true, then the system is feasible and S is nonempty. If the proposition is false, then the system is infeasible and S is empty. We can determine the truth value of the proposition using REDUCE + REDLOG to perform quantifier elimination:

% feasibility via quantifier elimination

load_package redlog;
rlset ofsf;

% define polynomial functions
f1 := x1^2+x2^2-1;
g1 := 3*x2-x1^3-2;
g2 := x1-8*x2^3;

% define existential formula
phi := ex({x1,x2}, f1=0 and g1>=0 and g2>=0);

% perform quantifier elimination
rlqe phi;

end;

The output of this script is the following:

The proposition is false and, thus, the system is infeasible.

__________

References

[1] Pablo Parrilo, Sum of Squares Programs and Polynomial Inequalities.

[2] Jacek Bochnak, Michel Coste, Marie-Françoise Roy, Real Algebraic Geometry, Springer, 1998.

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?

Like that Roman soldier in Pompeii

August 14, 2011

A passage from Oswald Spengler‘s Der Mensch und die Technik:

Wir sind in diese Zeit geboren und müssen tapfer den Weg zu Ende gehen, der uns bestimmt ist. Es gibt keinen andern. Auf dem verlorenen Posten ausharren ohne Hoffnung, ohne Rettung, ist Pflicht. Ausharren wie jener römische Soldat, dessen Gebeine man vor einem Tor in Pompeji gefunden hat, der starb, weil man beim Ausbruch des Vesuv vergessen hatte, ihn abzulösen. Das ist Größe, das heißt Rasse haben. Dieses ehrliche Ende ist das einzige, das man dem Menschen nicht nehmen kann.

A possible translation would be:

We are born into this time and must bravely follow the path to our destined end. There is no other way. Our duty is to hold on to the lost position, without hope, without rescue, like that Roman soldier whose bones were found in front of a door in Pompeii, who, during the eruption of Vesuvius, died at his post because they forgot to relieve him. That is greatness. That is what it means to be a thoroughbred. The honorable end is the one thing that can not be taken from a man.

__________

Source:

Oswald Spengler, Der Mensch und die Technik. Beitrag zu einer Philosophie des Lebens, München 1931.

LMI representation of convex polytopes

August 7, 2011

Suppose we are given a convex polytope P \subset \mathbb{R}^n defined by the intersection of m closed half-spaces

P := \{ x \in \mathbb{R}^n \mid A x \leq b \}

where A \in \mathbb{R}^{m \times n} and b \in \mathbb{R}^m. Note that set P is the solution set of a system of m non-strict linear inequalities in n variables. Let

p_i (x) := b_i - a_{i1} x_1 - a_{i2} x_2 - \dots - a_{in} x_n

so that set P can be described as follows

P = \{ x \in \mathbb{R}^n \mid p_1 (x) \geq 0 \land p_2 (x) \geq 0 \land \dots \land p_m (x) \geq 0\}

which allows us to conclude that P is a semialgebraic set [1]. The conjunction of the m non-strict linear inequalities p_i (x) \geq 0 is equivalent to the following linear matrix inequality (LMI) [2]

F (x) := \left[\begin{array}{cccc} p_1 (x) & 0 & \ldots & 0\\ 0 & p_2 (x) & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & p_m (x)\end{array}\right] \succeq 0

which, since the p_i polynomials are of degree 1, can be written as

F (x) = F_0 + x_1 F_1 + x_2 F_2 + \dots + x_n F_n \succeq 0

where the real diagonal m \times m matrices F_0, F_1, \dots, F_n are given by

F_0 := \left[\begin{array}{cccc} b_1 & 0 & \ldots & 0\\ 0 & b_2 & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & b_m\end{array}\right],    F_j = \left[\begin{array}{cccc} -a_{1j} & 0 & \ldots & 0\\ 0 & -a_{2j} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & -a_{m j}\end{array}\right]

for j = 1, 2,\dots, n. Thus, convex polytope P can be defined using the LMI F (x) \succeq 0, as follows

P = \{ x \in \mathbb{R}^n \mid F (x) \succeq 0 \}

where the F_0, F_1, \dots, F_n matrices are diagonal.

__________

References

[1] Jacek Bochnak, Michel Coste, Marie-Françoise Roy, Real Algebraic Geometry, Springer, 1998.

[2] Stephen Boyd, Laurent El Ghaoui, Eric Feron, and V. Balakrishnan, Linear Matrix Inequalities in System and Control Theory, SIAM, 1994.


Follow

Get every new post delivered to your Inbox.

Join 77 other followers