# Recent posts

## Run-time limits of comparison-based sorting algorithms

Posted by Diego Assencio on 2017.10.01 under Computer science (Algorithms)

Let $T(n)$ be the run-time of some comparison-based sorting algorithm for a sequence of length $n$. By "comparison-based", we mean a sorting algorithm which accesses input array elements only via comparisons as is the case for general-purpose sorting algorithms such as merge sort, quicksort, and heapsort. Such an algorithm does not know in advance which values the input array can store (as is for instance the case with bucket sort), but simply orders them using a predefined "less than" comparison function. In this post, we will show that regardless of how well designed an algorithm of this type is, its worst-case running time cannot be faster than $c n\log_2(n)$ as $n$ becomes large (for some constant $c \gt 0$). In more technical words, there exist constants $n_0 \geq 0$ and $c \gt 0$ such that $T(n) \geq cn\log_2(n)$ for $n \geq n_0$, and therefore $T(n) = \Omega(n\log_2(n))$.

This post will only deal with deterministic algorithms, but the results extend to randomized algorithms as well. The main idea of the proof is to consider only input arrays which are permutations of $S_n = \{1, 2, \ldots, n\}$. The total number of distinct arrays which represent permutations of $S_n$ is $n!$.

For a given deterministic, comparison-based sorting algorithm, let $k$ be the largest number of comparisons within which it can sort any of the $n!$ arrays mentioned above. Because the execution of the algorithm is based on the results of the value comparisons it performs, there are at most $2^k$ possible distinct execution paths for it. To clarify, the algorithm sorts an input array by repeatedly comparing pairs of values and deciding on how to proceed depending on which value is greater according to the comparison function. Each comparison generates two possible branches of execution, and since we know the total number of comparisons performed cannot exceed $k$, there can be at most $2^k$ distinct execution paths for the arrays we are considering.

Suppose now that $2^k \lt n!$. If that is the case, then the number of possible execution paths of the algorithm is smaller than the total number of input arrays we are considering. Therefore, there must be at least two distinct input arrays for which the algorithm executes following the exact same path (this is referred to as the pigeonhole principle in mathematics). This implies that when the algorithm is applied to these arrays, the same transformations (e.g. elements swaps) are applied to each of them, with the final result being the same: the values of $S_n$ ordered according to the predefined comparison function. This is not possible since the algorithm is deterministic: by inverting its steps from the final (sorted) result, we cannot arrive at two distinct initial inputs. Therefore $2^k \geq n!$.

The final step in the proof needs a simple fact about $n!$: $$\displaystyle n! = n(n-1)(n-2)\ldots\frac{n}{2}\left(\frac{n}{2} - 1\right) \ldots1 \geq \left(\frac{n}{2}\right)^{\frac{n}{2}}$$ Hence, since $2^k \geq n!$: $$\log_2(2^k) \geq \log_2(n!) \Longrightarrow k \geq \frac{n}{2}\log\left(\frac{n}{2}\right) \Longrightarrow k = \Omega(n\log_2(n))$$ Therefore, the best lower bound for the worst-case running time of the sorting algorithm is $T(n) = \Omega(n\log_2(n))$.

## The DFT of an image: phase vs. magnitude

Posted by Diego Assencio on 2017.09.01 under Computer science (Digital signal processing)

Let $X$ be an $M \times N$ matrix (a two-dimensional array) of complex numbers $x_{mn}$ for $m = 0, 1, \ldots, M - 1$ and $n = 0, 1, \ldots, N - 1$. The discrete Fourier transform of $X$ is another $M \times N$ matrix $\tilde{X}$ of complex numbers $\tilde{x}_{pq}$ such that: $$\tilde{x}_{pq} = \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} x_{mn} e^{-i2\pi pm/M} e^{-i2\pi qn/N} \label{post_81059bc883bd6156c2e83f83e2e38c2b_dft_2d}$$ for $p = 0, 1, \ldots, M - 1$ and $q = 0, 1, \ldots, N - 1$. The original matrix $X$ can be recovered from $\tilde{X}$ through its inverse discrete Fourier transform: $$x_{mn} = \displaystyle\frac{1}{MN} \sum_{p = 0}^{M - 1} \sum_{q = 0}^{N - 1} \tilde{x}_{pq} e^{i2\pi pm/M} e^{i2\pi qn/N}$$ Since a grayscale image can be represented as a matrix whose elements are 8-bit integers representing pixel values ranging from $0$ (black) to $255$ (white), we can compute its DFT using equation \eqref{post_81059bc883bd6156c2e83f83e2e38c2b_dft_2d}. Interestingly, the magnitude of the DFT of such an image does not contain much information, but its phase does. To clarify what this means, let us define two $M \times N$ matrices $\tilde{X}_{\textrm{mag}}$ and $\tilde{X}_{\textrm{phase}}$ whose elements are given by: $$\begin{eqnarray} (\tilde{X}_{\textrm{mag}})_{pq} &=& \left|\tilde{x}_{pq}\right| \label{post_81059bc883bd6156c2e83f83e2e38c2b_M}\\[5pt] (\tilde{X}_{\textrm{phase}})_{pq} &=& \tilde{x}_{pq} / \left|\tilde{x}_{pq}\right| \label{post_81059bc883bd6156c2e83f83e2e38c2b_P} \end{eqnarray}$$ Being a complex number, $\tilde{x}_{pq} = \rho_{pq} e^{i\phi_{pq}}$, where $\rho_{pq}$ is its magnitude and $\phi_{pq}$ is its phase. From equations \eqref{post_81059bc883bd6156c2e83f83e2e38c2b_M} and \eqref{post_81059bc883bd6156c2e83f83e2e38c2b_P}, we have that $(\tilde{X}_{\textrm{mag}})_{pq} = \rho_{pq}$ and $(\tilde{X}_{\textrm{phase}})_{pq} = e^{i\phi_{pq}}$, so the matrices $\tilde{X}_{\textrm{mag}}$ and $\tilde{X}_{\textrm{phase}}$ contain only magnitude and phase information of the elements of $\tilde{X}$ respectively.

If we compute the inverse DFTs of $\tilde{X}_{\textrm{mag}}$ and $\tilde{X}_{\textrm{phase}}$, what type of "images" will we recover? Our original image had only real values, but these inverse DFTs will in general not, so we will consider only their real components. Additionally, the pixel values of the resulting images will not necessarily fall within the closed interval $[0, 255]$, but we can clip the negative values (set them to zero) and rescale the nonnegative values to have them fall within this interval. An example of what the resulting images look like is shown on figure 1 (all figures in this post were generated using the script available here).

 Fig. 1: An original image (left) and the images produced by computing the inverse DFTs of the magnitude (center) and phase (right) components of the original image's DFT.

From figure 1, we can see that the inverse DFT of the magnitude matrix $\tilde{X}_{\textrm{mag}}$ produces a nearly black image, but the inverse DFT of the phase matrix $\tilde{X}_{\textrm{phase}}$ shows well-defined contours from the original image (if you cannot see them, try increasing the brightness of your screen or click on the figure to see a larger version of it). This indicates that the phase component of the image's DFT carries more information than the magnitude component. To better see this, we can apply histogram equalization to the images produced from the inverse DFTs. The results, shown in figure 2, serve as further confirmation of our hypothesis.

 Fig. 2: An original image (left) and the images produced (after histogram equalization) by computing the inverse DFTs of the magnitude (center) and phase (right) components of the original image's DFT.

## Determining types deduced by the compiler in C++

Posted by Diego Assencio on 2017.08.23 under Programming (C/C++)

Suppose you are debugging a C++ program and need to know what is the type which the compiler deduces for a certain expression. There are many complicated ways of doing this, but in this post, I will show you a very simple trick which allows you to determine types directly at compile time.

The trick is this: declare a class template but do not define it, then attempt to instantiate this class template with the expression whose type you are trying to determine. Here is an example:

/* class template declaration (no definition available) */
template<typename T>
class ShowType;

int main()
{
signed int x = 1;
unsigned int y = 2;

/* (signed int) + (unsigned int): what is the resulting type? */
ShowType<decltype(x + y)> dummy;

return 0;
}


In the code above, we are trying to determine the type deduced by the compiler when we add a signed int (x) and an unsigned int (y). This type is decltype(x + y). When the compiler attempts to create an instance of ShowType<decltype(x + y)>, it realizes this is not possible and indicates the problem with a very helpful error message:

error: aggregate ‘ShowType<unsigned int> dummy’ has incomplete type and cannot
be defined


In this message, the compiler (in my case, gcc) is telling us that it tried to create an instance of ShowType<unsigned int> but failed at it. Therefore, the type of the expression x + y is decltype(x + y) = unsigned int. This ŕesulting type comes directly from the integer addition rules specified in the C++ language.

Let's try a more interesting example. In C++, the type deduction rules for template parameters are complex. When in doubt, you can use the trick above to determine which type is deduced by the compiler for a certain template parameter:

/* class template declaration (no definition available) */
template<typename T>
class ShowType;

template<typename T>
void my_function(T x)
{
ShowType<T> dummy;

/* do something with x */
}

int main()
{
const int x = 3;
my_function(x);

return 0;
}


One common doubt which developers often have regarding the type T on my_function is: will it be deduced to be int or const int? As it turns out, since we are passing x by value, the compiler will deduce T to be int:

error: ‘ShowType<int> dummy’ has incomplete type


As a final example, let's consider auto. The rules for auto type deduction are usually the same as the ones for template types, but auto type deduction assumes that initializing expressions such as {1,2,3} represent initializer lists. Let's show that in practice:

#include <initializer_list>

template<typename T>
class ShowType;

int main()
{
auto x = {1,2,3};

ShowType<decltype(x)> dummy;

return 0;
}


The error message tells us what we expect:

error: aggregate ‘ShowType<std::initializer_list<int>> dummy’ has incomplete
type and cannot be defined


Notice that we need to pass an actual type (not an expression) to ShowType<T>, so in the examples above in which we wanted to determine the type of a certain expression (e.g. x + y), we needed to enclose the expression with decltype. On the second example, we already had the desired type T available on the definition of my_function, so we could use it directly.