## Archive for September, 2011

### Circular convolution in Haskell

September 5, 2011

Given two finite sequences, $x$ and $y$, of length $n$, their $n$-point circular convolution can be written in matrix form [1] as follows

$\left[\begin{array}{c} z_0\\ z_1\\ z_2\\ \vdots\\ z_{n-2}\\ z_{n-1}\end{array}\right] = \left[\begin{array}{cccccc} y_0 & y_{n-1} & y_{n-2} & \ldots & y_2 & y_1\\ y_1 & y_0 & y_{n-1} & \ldots & y_3 & y_2\\ y_2 & y_1 & y_0 & \ldots & y_4 & y_3\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ y_{n-2} & y_{n-3} & y_{n-4} & \ldots & y_0 & y_{n-1}\\ y_{n-1} & y_{n-2} & y_{n-3} & \ldots & y_1 & y_0\end{array}\right] \left[\begin{array}{c} x_0\\ x_1\\ x_2\\ \vdots\\ x_{n-2}\\ x_{n-1}\end{array}\right]$

where $z$ is also a sequence of length $n$. Note that the $n \times n$ matrix above is a circulant matrix, as each row is a circularly right-shifted version of the previous one. Each $z_k$ is obtained by taking the inner product of the $k$-th row of the matrix with the $x$ vector.

In Digital Signal Processing (DSP), we often think of finite sequences as vectors, which allows us to compute convolutions using matrix operations. Taking a more Computer Science-ish approach, one can think of finite sequences as finite lists, and compute convolutions using list operations. The matrix above can be thought of as a list of lists (where each row is a list). We know how to multiply a matrix by a vector, but how do we “multiply” a list of lists by a list?

I started teaching myself Haskell this Summer, and I am (obviously) looking for interesting problems at which to throw some functional firepower. Solutions looking for a problem. But, first, let us “begin at the beginning”…

__________

Reverse

Take another look at the matrix above and note that its last row, when viewed as a list, is the reversal of list $y$, which means that we need a function that reverses lists. Fortunately, we do not need to implement such a function, as there is the function reverse in Prelude that uses foldl and flip:

reverse :: [a] -> [a]
reverse = foldl (flip (:)) []

__________

Circular shifts

Yet once again, let us look at the matrix at the top of this post. Note that the 2nd row is a circularly right-shifted version of the 1st one, and that the 3rd row is a circularly right-shifted version of the 2nd one, et cetera. Also, note that the $(n-1)$-th row is a circularly left-shifted version of the $n$-th row, and so and so on. Thus, we need to implement functions that perform circular shifts.

The following function circularly right-shifts a list:

circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

This function merely takes the last element of the list and moves it to the beginning. We can circularly left-shift a list using the function:

circShiftL :: [a] -> [a]
circShiftL [] = []
circShiftL xs = (tail xs) ++ [head xs]

which moves the first element of a list to the end. Let us now test these functions on some simple lists using the GHCi interpreter:

*Main> circShiftR [0..9]
[9,0,1,2,3,4,5,6,7,8]
*Main> circShiftR (circShiftR [0..9])
[8,9,0,1,2,3,4,5,6,7]
*Main> (circShiftR . circShiftR) [0..9]
[8,9,0,1,2,3,4,5,6,7]
*Main> circShiftL [0..9]
[1,2,3,4,5,6,7,8,9,0]
*Main> (circShiftL . circShiftL) [0..9]
[2,3,4,5,6,7,8,9,0,1]

Note that in the 3rd and 5th inputs we composed the circular shifting functions to perform double right- and left-shifts.

__________

Iterate

Since each row of the matrix is a right-shifted version of the previous one or a left-shifted version of the next one, we can generate the matrix from its first or last rows. Putting it in terms of lists, we can generate a list of lists that represents the matrix from a single list using the iterate function in Prelude:

iterate :: (a -> a) -> a -> [a]
iterate f x =  x : iterate f (f x)

and the circular shifting functions we implemented before. To make things concrete, we can generate all four circular right-shifts of a list of length 4, as follows:

*Main> take 4 (iterate circShiftR [0..3])
[[0,1,2,3],[3,0,1,2],[2,3,0,1],[1,2,3,0]]

A word of caution is in order. Note that we take the first four lists, as

iterate f x

does generate an infinite list containing $x$, $f (x)$, $(f \circ f ) (x)$, $(f \circ f \circ f ) (x)$, et cetera. We do not need an infinite list of lists.

We now know how to generate a list of lists representing the matrix.

__________

Inner product

To obtain the circular convolution of $x$ and $y$, we multiply a matrix by a (column) vector, which consists of computing $n$ inner products. However, in this post we are thinking in terms of lists. We thus need to implement a function that computes the inner product of two (finite) lists, as follows:

innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

where we first use zipWith to obtain a list of the $x_i y_i$ products, and then use sum to compute the sum $\sum_i x_i y_i$. Let us test the inner product function:

*Main> let xs = [1..4]
*Main> let ys = [0..3]
*Main> zipWith (*) xs ys
[0,2,6,12]
*Main> sum(zipWith (*) xs ys)
20

__________

Circular convolution

We finally have all we need to compute the circular convolution.

For example, suppose we want to compute the circular convolution of two 4-point sequences, say, $x = (1,1,1,1)$ and $y = (0,1,2,3)$. Using matrix-vector multiplication, we obtain:

$\left[\begin{array}{cccc} 0 & 3 & 2 & 1\\ 1 & 0 & 3 & 2\\ 2 & 1 & 0 & 3\\ 3 & 2 & 1 & 0\end{array}\right] \left[\begin{array}{c} 1\\ 1\\ 1\\ 1\end{array}\right] = \left[\begin{array}{c} 6\\ 6\\ 6\\ 6\end{array}\right]$

Alternatively, we can compute the circular convolution by representing sequences by lists and using the functions we have implemented before:

*Main> -- define lists
*Main> let xs = [1,1,1,1]
*Main> let ys = [0,1,2,3]
*Main> -- reverse and right-shift ys
*Main> let ys' = (circShiftR . reverse) ys
*Main> ys'
[0,3,2,1]
*Main> -- compute list of lists representing the matrix
*Main> let yss = take 4 (iterate circShiftR ys')
*Main> yss
[[0,3,2,1],[1,0,3,2],[2,1,0,3],[3,2,1,0]]
*Main> -- compute inner product of xs with each list in yss
*Main> map (innerProd xs) yss
[6,6,6,6]

Note that I wrote comments on the GHCi command line in order to improve readability. We are now ready to implement the function that computes the circular convolution:

circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
where
n   = length xs
ys' = (circShiftR . reverse) ys
yss = take n (iterate circShiftR ys')

How does it work? We first take sequence $y$, reverse and right-shift it (to obtain the 1st row of the matrix), generate the remaining rows by right-shifting the previous ones, and finally compute the inner product of each row with sequence $x$ using map. Note that

innerProd xs

is a partial application of the inner product function, as we provide only the first argument. We could define a function $g$ as follows:

*Main> let g = innerProd [1,1,1,1]
*Main> :type g
g :: [Integer] -> Integer
*Main> g [0..3]
6

Note that $g$ takes a list of integers and returns its inner product with the list of four ones. Putting it all together in a .hs file, we have:

-- circular right-shift
circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

-- inner product of two lists
innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

-- circular convolution of two lists
circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
where
n   = length xs
ys' = (circShiftR . reverse) ys
yss = take n (iterate circShiftR ys')

We use several higher-order functions in this implementation: map, iterate, zipWith. We also use function composition and partial function application. A lot of beautiful ideas in only a dozen lines of code!

I am a Haskell neophyte, and I make no claims that the code above cannot be improved. If you have any constructive suggestions, I would be happy to hear them. Please use the HTML tags <pre> and </pre> to post code in the comments.

__________

Testing

Last but not least, let us carry out a couple of tests using examples from Mitra’s book [1]:

*Main> -- problem 5.2 b)
*Main> let xs = [-1,5,3,0,3]
*Main> let hs = [-2,0,5,3,-2]
*Main> circConv xs hs
[1,-1,-2,16,26]
*Main> -- problem 5.45 b)
*Main> let gs = [2,-1,3,0]
*Main> let hs = [-2,4,2,-1]
*Main> circConv gs hs
[3,7,-6,8]

It appears to be working! Please let me know if you find any bugs.

__________

References

[1] Sanjit K. Mitra, Digital Signal Processing: a computer-based approach, 3rd edition, McGraw-Hill, 2005.

### Deciding 2-colorability via quantifier elimination

September 1, 2011

A vertex coloring of a graph $G = (V,E)$ is a map $\kappa : V \to C$ that assigns a color to each vertex of $G$. The coloring is proper if adjacent vertices are assigned distinct colors. The graph is 2-colorable if there exists a proper coloring whose color set $C$ has two elements.

For instance, consider a graph $G = (V,E)$ whose vertex set is $V := \{1, 2, 3\}$ and whose edge set is $E := \{\{1,2\}, \{1,3\}, \{2,3\}\}$. A pictorial representation of this graph is depicted below

Note that this graph is complete. We can assign two distinct colors to any two vertices, but we cannot assign a color to the third vertex that will yield a proper coloring. Therefore, the graph is not 2-colorable. To illustrate, we use the color set $C := \{\text{pink}, \text{blue}\}$ and generate all $2^3 = 8$ possible vertex 2-colorings, as depicted below

Notice that none of these 2-colorings is proper, as expected.

As discussed by De Loera et alia [1], the graph vertex-coloring problem can be modeled using polynomial equations. We introduce variables $x_1, x_2, x_3$, where $x_i = -1$ if vertex $i$ is assigned the color $\text{pink}$, and $x_i = 1$ if vertex $i$ is assigned the color $\text{blue}$. Given an edge $\{i,j\} \in E$, we have that:

• if both vertices are pink, we get $x_i + x_j = -2$.
• if both vertices are blue, we get $x_i + x_j = 2$.
• if the vertices have distinct colors, we get $x_i + x_j = 0$.

Thus, stating that adjacent vertices must be assigned distinct colors is equivalent to imposing the constraints $x_i + x_j = 0$ for all $\{i,j\} \in E$. We thus obtain a system of three equations

$\begin{array}{rl} x_1 + x_2 &= 0\\ x_1 + x_3 &= 0\\ x_2 + x_3 &= 0\end{array}$

where $x_1, x_2, x_3 \in \{-1,1\}$. This is still a combinatorial problem, though. We can transform it into an algebraic-geometric problem. Here is a truly beautiful trick [1]: the constraint $x_i \in \{-1,1\}$ can be encoded as $x_i^2 = 1$ with $x_i \in \mathbb{R}$. How quaint!! We then obtain a system of six equations with polynomials in $\mathbb{R}[x_1, x_2, x_3]$

$\begin{array}{rl} x_1 + x_2 &= 0\\ x_1 + x_3 &= 0\\ x_2 + x_3 &= 0\\ x_1^2 - 1 &= 0\\ x_2^2 - 1 &= 0\\ x_3^2 -1 &= 0\end{array}$

whose solution set, which we denote by $S$, is algebraic [2]. The given graph is 2-colorable if and only if the system of polynomial equations above is feasible, i.e., if set $S$ is nonempty [1]. The following REDUCE + REDLOG script uses quantifier elimination to decide if $S$ is nonempty:

% feasibility via quantifier elimination

load_package redlog;
rlset ofsf;

% define polynomial functions
f1 := x1+x2;
f2 := x1+x3;
f3 := x2+x3;
f4 := x1^2-1;
f5 := x2^2-1;
f6 := x3^2-1;

% define existential formula
phi := ex({x1,x2,x3}, f1=0 and f2=0 and f3=0 and
f4=0 and f5=0 and f6=0);

% perform quantifier elimination
rlqe phi;

end;

This script produces the truth value false. Thus, the solution set $S$ is empty, which means that the system of polynomial equations is infeasible. We conclude that the given graph is not 2-colorable.

__________

References

[1] Jesús De Loera, Jon Lee, Susan Margulies, Shmuel Onn, Expressing Combinatorial Optimization Problems by Systems of Polynomial Equations and the Nullstellensatz, 2007.

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