Algebra « Rod Carvalho

Archive for the ‘Algebra’ Category

Linear dynamical systems over finite rings II

July 28, 2011

Last week, we studied some properties of the following finite linear dynamical system

$\left[\begin{array}{cc} x_1^+\\ x_2^+\end{array}\right] = \left[\begin{array}{cc} 2 & 3\\ 1 & 0\\\end{array}\right] \left[\begin{array}{c} x_1 \\ x_2\\\end{array}\right]$

where $x := (x_1, x_2)$ lives in $\mathbb{Z}_4^2$ , where $\mathbb{Z}_4 := \{0, 1, 2, 3\}$ . Let $A \in \mathbb{Z}_4^{2 \times 2}$ be defined by

$A := \left[\begin{array}{cc} 2 & 3\\ 1 & 0\\\end{array}\right]$

so that we can rewrite the state update equation as $x^+ = A x$ .

In the 1st installment, we obtained qualitative information about the dynamical system from visual inspection of its state transition diagram, which is clearly impossible for finite dynamical systems with state sets of sufficiently large cardinality. Thus, in this post, we will study the same dynamical system using only algebra. The goal is to see how far we can go without knowing what the state transition diagram looks like.

Before we embark on this journey, please take a good look at the addition and multiplication tables of the finite ring $(\mathbb{Z}_4, +, \times)$

__________

Time reversibility

Suppose there exists a matrix $B \in \mathbb{Z}_4^{2 \times 2}$ such that $B A = I$ , where $I \in \mathbb{Z}_4^{2 \times 2}$ is the identity matrix. Then, left-multiplying both sides of the state update equation by $B$ , we obtain $B x^+ = B A x = I x$ , which is equivalent to $x = B x^+$ . Hence, if matrix $B$ exists, we can propagate the state backwards in time, and the dynamical system is time reversible.

Matrix equation $B A = I$ can be written as

$\left[\begin{array}{cc} b_{11} & b_{12}\\ b_{21} & b_{22}\\\end{array}\right] \left[\begin{array}{cc} 2 & 3\\ 1 & 0\\\end{array}\right] = \left[\begin{array}{cc} 1 & 0\\ 0 & 1\\\end{array}\right]$

and, vectorizing both sides, we obtain a linear system of equations $(A^T \otimes I)\, \mbox{vec}(B) = \mbox{vec}(I)$ , i.e., we have

$\left[\begin{array}{cccc} 2 & 0 & 1 & 0\\ 0 & 2 & 0 & 1\\ 3 & 0 & 0 & 0\\ 0 & 3 & 0 & 0\\\end{array}\right] \left[\begin{array}{c} b_{11}\\ b_{21}\\ b_{12}\\ b_{22}\\\end{array}\right] = \left[\begin{array}{c} 1\\ 0\\ 0\\ 1\\\end{array}\right]$

From the 4th and last equation, $3 b_{21} = 1$ , we get $b_{21} = 3$ , as $3 \times 3 = 1$ in modulo 4 arithmetic. From the 3rd equation, $3 b_{11} = 0$ , we get $b_{11} = 0$ . The 2nd equation, $2 b_{21} + b_{22} = 0$ , becomes then $2 + b_{22} = 0$ , which yields $b_{22} = 2$ . Finally, the 1st equation, $2 b_{11} + b_{12} = 1$ , becomes $b_{12} = 1$ . Let us define $A^{-1} := B$ . Hence, we have

$A^{-1} = \left[\begin{array}{cc} 0 & 1\\ 3 & 2\\\end{array}\right]$

which is the (left and right) inverse of matrix $A$ . We then have the (backwards in time) state update equation $x = A^{-1} x^+$ . If we draw the state transition diagram for this linear dynamical system, we obtain the diagram we obtained last week, but with the direction of the arrows reversed. If you do not believe me, then define $g (x) := A^{-1} x$ and run the following Python script:

# define function g
def g (x):
    g1 = x[1] % 4
    g2 = (3 * x[0] + 2 * x[1]) % 4
    return (g1, g2)

print "State transitions:\n"
for x1 in [0,1,2,3]:
    for x2 in [0,1,2,3]:
        print "g(%d,%d) = (%d,%d)" % (x1, x2, g((x1,x2))[0], g((x1,x2))[1])

__________

Fixed points

As mentioned last week, we can find fixed points of the dynamical system $x^+ = A x$ by solving the equation $A x = x$ . Since $-1 = 3$ in modulo 4 arithmetic, let us add $3 x$ to both sides of the equation, which yields $(A + 3 I) x = 0_2$ , i.e.,

$\left[\begin{array}{cc} 1 & 3\\ 1 & 3\end{array}\right] \left[\begin{array}{c} x_1\\ x_2 \end{array}\right] = \left[\begin{array}{c} 0\\ 0 \end{array}\right]$

and, since one of the equations is redundant, we are left with $x_1 + 3 x_2 = 0$ , which is equivalent to $x_1 = x_2$ . The set of fixed points is, thus,

$\{ (0,0), (1,1), (2,2), (3,3)\}$ .

Fixed points are periodic points of period equal to 1. Let us now find periodic points of period greater than 1.

__________

Periodic points of period equal to 2

To find periodic points of period equal to 2, we must solve the equation $A^2 x = x$ , where

$A^2 = \left[\begin{array}{cc} 3 & 2\\ 2 & 3\\\end{array}\right]$ .

Once again, let us add $3 x$ to both sides of the equation, which yields $(A^2 + 3 I) x = 0_2$ , i.e.,

$\left[\begin{array}{cc} 2 & 2\\ 2 & 2\end{array}\right] \left[\begin{array}{c} x_1\\ x_2 \end{array}\right] = \left[\begin{array}{c} 0\\ 0 \end{array}\right]$

and, eliminating the redundant equation, we are left with $2 x_1 + 2 x_2 = 0$ . Since $-2 = 2$ in modulo 4 arithmetic, we finally obtain the equation $2 x_1 = 2 x_2$ . If we were solving equations over $\mathbb{R}$ , we could just multiply both sides of $2 x_1 = 2 x_2$ by $1 / 2$ and obtain $x_1 = x_2$ . However, we are working with the finite ring $(\mathbb{Z}_4, +, \times)$ , for which the element 2 does not have a multiplicative inverse! Let us write the equation for various values of $x_2$ :

if $x_2 = 0$ , then we obtain the equation $2 x_1 = 0$ , which yields $x_1 \in \{0, 2\}$ . We thus have points $(0, 0)$ and $(2,0)$ .

if $x_2 = 1$ , then we obtain the equation $2 x_1 = 2$ , which yields $x_1 \in \{1, 3\}$ . We thus have points $(1, 1)$ and $(3,1)$ .

if $x_2 = 2$ , then we obtain the equation $2 x_1 = 0$ , which yields $x_1 \in \{0, 2\}$ . We thus have points $(0, 2)$ and $(2,2)$ .

if $x_2 = 3$ , then we obtain the equation $2 x_1 = 2$ , which yields $x_1 \in \{1, 3\}$ . We thus have points $(1, 3)$ and $(3,3)$ .

To summarize, we have eight periodic points of period equal to 2

$\{ (0,0), (2,0), (1,1), (3,1), (2,2), (0,2), (3,3), (1,3)\}$

of which, 4 points are fixed points. Removing these fixed points, we are left with 4 periodic points of fundamental period equal to 2

$\{ (2,0), (3,1), (0,2), (1,3)\}$ .

We finally conclude that the finite dynamical system under study has two cycles of length equal to 2, namely

$(2,0) \mapsto (0,2) \mapsto (2,0)$

$(3,1) \mapsto (1,3) \mapsto (3,1)$

__________

Periodic points of period equal to 3

To find periodic points of period equal to 3, we must solve the equation $A^3 x = x$ , where

$A^3 = \left[\begin{array}{cc} 0 & 1\\ 3 & 2\\\end{array}\right]$ .

Note that $A^3 = A^{-1}$ . Hence, $A^4 = A^{-1} A = I$ , and $A^5 = A$ . In order to solve equation $A^3 x = x$ , let us (yet once again) add $3 x$ to both sides of the equation, which yields $(A^3 + 3 I) x = 0_2$ , i.e.,

$\left[\begin{array}{cc} 3 & 1\\ 3 & 1\end{array}\right] \left[\begin{array}{c} x_1\\ x_2 \end{array}\right] = \left[\begin{array}{c} 0\\ 0 \end{array}\right]$

and, eliminating the redundant equation, we are left with $3 x_1 + x_2 = 0$ . Since $-1 = 3$ in modulo 4 arithmetic, we finally obtain the equation $3 x_1 = 3 x_2$ . Fortunately, element 3 has a multiplicative inverse, which is itself (as $3 \times 3 = 1$ ). Multiplying both sides by 3, we obtain $x_1 = x_2$ , an equation we solved already when looking for the fixed points. Since the only solutions of this equation are the four fixed points, we conclude that there are no periodic points of fundamental period equal to 3, i.e., there are no cycles of length 3.

__________

Periodic points of period equal to 4

To find periodic points of period equal to 4, we must solve the equation $A^4 x = x$ . Nevertheless, we know that $A^4 = I$ and, hence, $A^4 x = x$ reduces to $x = x$ , which is satisfied by every $x \in \mathbb{Z}_4^2$ , of course. In other words, every single point in the state set is a periodic point of period equal to 4. Some of them have fundamental period equal to 1, some others have fundamental period equal to 2. Removing from the state set these two sets, we are left with the periodic points of fundamental period equal to 4

$\{ (1,0), (2,1), (3,2), (0,3), (3,0), (2,3), (1,2), (0,1)\}$ .

We finally conclude that the finite dynamical system under study has two cycles of length equal to 4, namely

$(1,0) \mapsto (2,1) \mapsto (3,2) \mapsto (0,3) \mapsto (1,0)$

$(3,0) \mapsto (2,3) \mapsto (1,2) \mapsto (0,1) \mapsto (3,0)$

The hunt for periodic points is over. Since $A^5 = A$ , equation $A^5 x = x$ reduces to equation $A x = x$ , which we already solved. Likewise, $A^6 x = x$ boils down to $A^2 x = x$ , which we already solved too. We have found all the cycles.

__________

Concluding remarks

All the qualitative information about the finite dynamical system that we were able to obtain by taking a look at the state transition diagram can instead be obtained algebraically.

Last week, we found the solutions of the equations $A^k x = x$ using sheer brute force, i.e., via exhaustive search over the entire state set. By contrast, in this post we actually did solve the equations using algebra. Intuitively, the latter should be computationally cheaper. If this is correct, then the question is: how much cheaper?

Tags:Discrete Dynamical Systems, Finite Dynamical Systems, Finite Rings, Fixed Points, Modular Arithmetic, Periodic Points
Posted in Algebra, Dynamical Systems | 2 Comments »

Linear dynamical systems over finite rings

July 18, 2011

Consider the following linear dynamical system (example 3.4 in [1])

$\left[\begin{array}{cc} x_1^+\\ x_2^+\end{array}\right] = \left[\begin{array}{cc} 2 & 3\\ 1 & 0\\\end{array}\right] \left[\begin{array}{c} x_1 \\ x_2\\\end{array}\right]$

where the state vector $x := (x_1, x_2)$ travels in the finite state space $\mathbb{Z}_4^2$ , where $\mathbb{Z}_4 := \{0, 1, 2, 3\}$ . Since the state space is finite, let us henceforth call it state set. Let $A \in \mathbb{Z}_4^{2 \times 2}$ be defined by

$A := \left[\begin{array}{cc} 2 & 3\\ 1 & 0\\\end{array}\right]$

so that we can write the state update equation in the more compressed form $x^+ = A x$ . For simplicity, let us introduce function $f : \mathbb{Z}_4^2 \to \mathbb{Z}_4^2$ defined by $f (x) := A x$ . Hence, we have $x^+ = f (x)$ .

Do note that we are working with the finite ring $(\mathbb{Z}_4, +, \times)$ , whose addition and multiplication tables are as follows

and, therefore, when we update each state via $x_i^+ = a_{i1} x_1 + a_{i2} x_2$ , we use addition and multiplication modulo 4.

When we study discrete-time dynamical systems over $\mathbb{R}$ , we are usually interested in finding equilibrium points, fixed points and limit cycles. For dynamical systems over a finite ring, what are the properties of interest? And how can one study them?

__________

State transition graph

The state set, $\mathbb{Z}_4^2$ , which has $4^2 = 16$ elements, is a subset of the 2-dimensional integer lattice $\mathbb{Z}^2$ , as depicted below

Computing $x^+ = f(x)$ for each $x$ in the state set and forming ordered pairs $(x, f(x))$ , we can build the state transition graph, whose vertex set is the state set $\mathbb{Z}_4^2$ , and whose edge set is the set $\{(x, f(x))\}_{x \in \mathbb{Z}_4^2}$ . A pictorial representation of this graph is shown below

Note that the out-degree of each vertex is 1. This is due to the fact that function $f$ maps states to states, not states to sets of states. In other words, we do not allow non-determinism.

Visual inspection of the state transition diagram allows us to conclude that, interestingly, the in-degree of each vertex is also 1, i.e., function $f$ is injective and surjective. Thus, $f$ is invertible, which means that the dynamical system we are studying is time reversible! Switch the direction of each arrow in the diagram, and we obtain a new state transition diagram. Hence, if are given the state at a certain time step, we can determine the entire history of the system, both past and future!

A state $x$ is a fixed point if $f(x) = x$ , i.e., if $x^+ = x$ . Taking a look at the state transition diagram, we see that there are four vertices with self-loops, which means that the system has four fixed points

$\{ (0,0), (1,1), (2,2), (3,3)\}$ .

If the system starts at a fixed point $x$ , it will remain there forever

$x \mapsto f(x) = x \mapsto f(x) = x \mapsto \dots$

and, thus, the set of all fixed points is an invariant set. Note that a fixed point is a periodic point of period equal to 1.

Observing the state transition diagram again, we notice the existence of loops of various lengths. We have two loops of length equal to 2 and two loops of length equal to 4, depicted below in blue and magenta, respectively

Do note, however, that we also have four loops of length equal to 1, which are the self-loops corresponding to the fixed points, but we already mentioned those.

For the dynamical system under study, we could obtain qualitative information about its behavior by visual inspection of its state transition diagram. Suppose, for example, that we want to study a dynamical system whose state set is $\mathbb{Z}_4^{16}$ , which has over four billion states. It is evident that drawing diagrams and counting loops is an approach that does not scale! Can one use Algebra to obtain qualitative information about the dynamical system being studied?

__________

Finding periodic points

If state $x$ is a fixed point, then $f (x) = x$ or, equivalently, $A x = x$ , as we defined $f (x) := A x$ . Hence, we can find fixed points by solving the equation $A x = x$ .

A state $x$ is a periodic point of period equal to $k$ if $f^{(k)} (x) = x$ , where $f^{(k)} = f \circ f^{(k-1)}$ and $f^{(1)} = f$ . Since $f (x) := A x$ , we have that $f^{(k)} (x) = A^k x$ . Thus, we can find periodic points of period equal to $k$ by solving the equation $A^k x = x$ .

However, do note that fixed points, i.e., periodic points of period equal to 1, are also periodic points of period equal to 2, equal to 3, equal to 4, etc. Therefore, let the fundamental period be the smallest natural number $k$ for which $A^k x = x$ . For example, the set of all periodic points of fundamental period equal to 2 is

$\{x \in \mathbb{Z}_4^2 \mid A^2 x = x \land A x \neq x\}$

and the set of all periodic points of fundamental period equal to 3 is

$\{x \in \mathbb{Z}_4^2 \mid A^3 x = x \land A x \neq x\}$ .

However, the periodic points of period equal to 2 are also periodic points of period equal to 4. Hence, the set of periodic points of fundamental period equal to 4 is

$\{x \in \mathbb{Z}_4^2 \mid A^4 x = x \land A^2 x \neq x \land A x \neq x\}$ .

But, how do we solve equations of the form $A^k x = x$ over the ring $(\mathbb{Z}_4, +, \times)$ ? A possibility is to search the whole state set and find the states for which the equations hold. That is precisely what the following Python script does:

# define function f
def f (x):
    f1 = (2 * x[0] + 3 * x[1]) % 4
    f2 = x[0] % 4
    return (f1, f2)

# create state set
X = []
for x1 in [0,1,2,3]:
    for x2 in [0,1,2,3]:
        X.append((x1,x2))

print "State transitions:\n"
for x in X:
    x1, x2, f1, f2 = x[0], x[1], f(x)[0], f(x)[1]
    print "f(%d,%d) = (%d,%d)" % (x1, x2, f1, f2)

# find periodic points of fundamental periods 1, 2, 3, 4
PP1, PP2, PP3, PP4 = [], [], [], []
for x in X:
    if f(x) == x:
        PP1.append(x)
    elif f(f(x)) == x:
        PP2.append(x)
    elif f(f(f(x))) == x:
        PP3.append(x)
    elif f(f(f(f(x)))) == x:
        PP4.append(x)

print "\nLists of periodic points:\n"
print "\nFundamental period equal to 1:\n %s" % PP1
print "\nFundamental period equal to 2:\n %s" % PP2
print "\nFundamental period equal to 3:\n %s" % PP3
print "\nFundamental period equal to 4:\n %s" % PP4

The output of this script is the following:

State transitions:

f(0,0) = (0,0)
f(0,1) = (3,0)
f(0,2) = (2,0)
f(0,3) = (1,0)
f(1,0) = (2,1)
f(1,1) = (1,1)
f(1,2) = (0,1)
f(1,3) = (3,1)
f(2,0) = (0,2)
f(2,1) = (3,2)
f(2,2) = (2,2)
f(2,3) = (1,2)
f(3,0) = (2,3)
f(3,1) = (1,3)
f(3,2) = (0,3)
f(3,3) = (3,3)

Lists of periodic points:

Fundamental period equal to 1:
[(0, 0), (1, 1), (2, 2), (3, 3)]

Fundamental period equal to 2:
[(0, 2), (1, 3), (2, 0), (3, 1)]

Fundamental period equal to 3:
[]

Fundamental period equal to 4:
[(0, 1), (0, 3), (1, 0), (1, 2), (2, 1), (2, 3), (3, 0), (3, 2)]

Note that there are no periodic points of fundamental period equal to 3, as a quick glance at the state transition diagram will tell.

__________

References

[1] O. Colón-Reyes, A. Jarrah, R. Laubenbacher, B. Sturmfels, Monomial Dynamical Systems over Finite Fields, May 2006.

Tags:Discrete Dynamical Systems, Finite Dynamical Systems, Finite Rings, Fixed Points, Periodic Points
Posted in Algebra, Dynamical Systems | 2 Comments »

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 »

Rod Carvalho

Archive for the ‘Algebra’ Category

Linear dynamical systems over finite rings II

Linear dynamical systems over finite rings

Sum of squares decomposition II

Sum of squares decomposition

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 ‘Algebra’ Category

Linear dynamical systems over finite rings II

Linear dynamical systems over finite rings

Sum of squares decomposition II

Sum of squares decomposition

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”