This is a second post on sorting algorithms. In this post we talk about four common algorithms that have average time complexity \(O(n*log(n))\), which is a lot better than \(O(n^2)\).
Most of the sorting methods in this post are by comparion, except bucket sort.
Volume 3 of Don Knuth’s book series “The Art of Computer Programming” compares sorting speed for an “artificial typical” computer invented by the author. Some of the average case results are in the following plot:
Quicksort
The quicksort algorithm reminds me of binary search.
In binary search on a sorted array, we half the search range by comparing target with the median of the search range.
Quicksort works by sorting as if we only care about sorting into two parts: larger than (something) and smaller than (something). The element that is used to split the original array into two parts is commonly called the “pivot”.
Then we apply the same method to each of the two parts, and continue until each of the parts has only 1 element (cannot be splited anymore).
Because we are divide-sorting in half every iteration, time complexity is \(n*log(n)\) on average.
Worst case time complexity
In the unfortunate case that the pivot is the largest or the smallest, then our intent to half the array fails:
• Then one subarray is always empty.
• The second subarray contains n − 1 elements, i.e. all the elements other than the pivot.
• Quicksort is recursively called only on this second group.
To do the first partitioning, there are \(n − 1\)$$ comparisons.
• At each subsequent step, the number of comparisons is one less, so that
\(T(n) = T(n − 1) + (n − 1)\);
\(T(1) = 0\).
• “Telescoping” \(T(n) − T(n − 1) = n − 1\):
\[T(n)+T(n − 1)+T(n − 2)+. . .+T(3)+T(2) −T(n − 1)−T(n − 2)−. . .−T(3)−T(2)− T(1) = (n − 1) + (n − 2) +. . .+ 2 + 1 − 0\] \[T(n)= (n − 1) + (n − 2) +. . .+ 2 + 1 =(n−1)n/2\]This yields that \(T(n) ∈ Ω(n^2)\).
Quicksort is slow on sorted or nearly sorted data.
It is fast on the “randomly scattered” pivots, excpet with no guarantees.
Quicksort implementations
Below quicksort implementation is simple but not optimized. We use recursion. The “initial” condition is that if input has only 1 element, then stop. Otherwise, we split input in half according to whether elements are bigger or smaller than pivot.
1. First shot
An error I made before was to use A[-1] as the pivot instead of A.pop().
A[-1] does not remove the last element (used as pivot) from the comparison.
A.pop() removes the last element.
To make sure we fully understand what’s going on, we print out each step.
In the second iteration of the sorting, the right array has only 1 element, 100. So recursion stops for the right array.
The new left array L is now [3, 9, 1, 0, -10]. Because nothing is smaller than the pivot -10, L is empty. The new right array has to carry on and do all the work.
Each time we see A is an empty array, we know it was recursed from an empty left or right array.
def quickSort(A):
print("\nA:", A)
if len(A) < 2: # base case
return A
else:
L, R = [], []
print("pivot is ", A[-1])
pivot = A.pop()
for num in A: # pivot is not in A anymore
if num < pivot:
L.append(num)
else:
R.append(num)
print("L:", L)
print("R:", R)
return quickSort(L) + [pivot] + quickSort(R)
nums = [100, 3, 9, 1, 0, -10, 87]
print(quickSort(nums))
# A: [100, 3, 9, 1, 0, -10, 87]
# pivot is 87
# L: [3, 9, 1, 0, -10]
# R: [100]
# A: [3, 9, 1, 0, -10]
# pivot is -10
# L: []
# R: [3, 9, 1, 0]
# A: []
# A: [3, 9, 1, 0]
# pivot is 0
# L: []
# R: [3, 9, 1]
# A: []
# A: [3, 9, 1]
# pivot is 1
# L: []
# R: [3, 9]
# A: []
# A: [3, 9]
# pivot is 9
# L: [3]
# R: []
# A: [3]
# A: []
# A: [100]
# [-10, 0, 1, 3, 9, 87, 100]
2. quicksort using list comprehension
One problem with the above code is that at the end, our input array is altered: it is missing the last element because of pop().
It is easy to fix that, we use slicing instead of pop().
Below snippet is adapted from the book “Python Cookbook”.
It uses only 2 lines of code. Simple!
# use the rightmost as pivot
def qsort(A):
if len(A) <= 1: return A
return qsort([lt for lt in A[:-1] if lt < A[-1]]) + A[-1:] + \
qsort([ge for ge in A[:-1] if ge >= A[-1]])
print(qsort(nums))
# use the leftmost as pivot
def qsort(A):
if len(A) <= 1: return A
return qsort([lt for lt in A[1:] if lt < A[0]]) + A[0:1] + \
qsort([ge for ge in A[1:] if ge >= A[0]])
print(qsort(nums))
Below quicksort implementation has improvement in space complexity. It uses the rightmost of range as pivot, which is not the best method because the rightmost can be the largest or smallest, and has higher likelihood to be the extreme values than if we use some methods to prevent it.
def partition(A, l, r):
i = l - 1
pivot = A[r] # use the rightmost of range as pivot (not the best method)
for j in range(l, r):
if A[j] <= pivot:
i += 1 # if a number is smaller than pivot, then pointer increments by 1
A[i], A[j] = A[j], A[i] # swap with that smaller number with the one pointer is at
A[i+1], A[r] = A[r], A[i+1]
return i+1
def quicksort(A, l=0, r = None):
if r==None:
r = len(A) - 1
if l < r:
pivot = partition(A, l,r)
quicksort(A, l, pivot-1)
quicksort(A, pivot+1, r)
return A
nums = [100, 3, 9, 1, 0]
quicksort(nums,l=0)
Preventing worst case
Given the importance of pivot, below code takes a random element as pivot instead of what happens to be at the leftmost or rightmost position. This code is adapted from this stackoverflow post.
This example implements quicksort inplace (\(O(n)\) space complexity). The pivot is picked randomly pivotIdx = random.randint(l, r). Another popular method is to take the median of 3 numbers.
It has two parts:
1. quicksort body. It takes input array, left, right position, randomly pick an element as the pivot.
-
It gives the function partitionSort input array, left, right and pivot indices, to do its partition-sort. Although I almost never made any errors on the index of l begins at 0, I had often made the mistakes of the right index initial value, which is the last index of array, len(A) - 1, not the number of elements.
-
The partitionSort function sorts input array into two subarrays: smaller than pivot, and greater than pivot, and returns the position of the pivot where it ended up.
- The condition if r == None: r = len(A) - 1 takes care of both worlds: initial parameters and subsequent recursion parameters. When we run quickSort(nums) , we only provide the input array. The inital run uses the default values l = 0, r = None. herefore, r = the last index of input array.
- It is critical to have if r <= l: return to control where to stop.
It then recurse on the two subarrays.
Note: when recurse, do think beyond the first run and what parameters work for the general case in recursive runs.
quickSort(A, l, p - 1) works for the genearl case. Whereas quickSort(A, 0, p - 1) won’t. Similary, quickSort(A, p + 1, r) works for the genearal case. Whereas quickSort(A, p + 1) won’t.
2. partitionSort: It takes input array, left, right and pivot indices, to do its partition-sort.
- Station the pivot at the leftmost place by swapping A[l] with pivot.
- Then it does its partition-sort by keeping track of two things: i and j both begin as l + 1 because they need to leave the pivot alone.
- i keeps track of the less-than-pivot elements, and moves by 1 to the right each time a less-than-pivot element is moved to the left subarray.
- j does the looping from l + 1 all the way to r (inclusive), and checks one by one if less than pivot. If less than pivot, A[j] swaps with A[i], and i increments by \(1\).
Note: Do not omit the \(=\) in while j <= r:. Because we may have duplicates in the array. And because everyone in the subarray needs to compare with the pivot, including the one in the rightest.
On the other hand, result is correct whether we put \(=\) in if A[j] < pivot:.
import random
def partitionSort(A, l, r, pivotIdx):
'''
sort into two parts, smaller than pivot,
'''
if not (l <= pivotIdx <= r):
raise ValueError("pivot index must be within l and r")
# move pivot out of the way
pivot = A[pivotIdx]
# now in A, pivot is moved to the leftmost place
A[pivotIdx], A[l] = A[l], A[pivotIdx]
# adding 1 because we don't touch the pivot.
i = l + 1
j = l + 1
while j <= r:
if A[j] < pivot:
A[i], A[j] = A[j], A[i]
i += 1 # i represent small item count, but we would need to decrement it by 1 for pivot in the end of the function
j += 1
# put pivot where it belongs because we have i-1 items small or equal to pivot
A[l], A[i -1] = A[i -1], A[l]
return i - 1
def quickSort(A, l = 0, r = None):
if r == None: # takes care of both worlds: initial parameters and subsequent recursion parameters
r = len(A) - 1
if r <= l: # controls where to stop: r cannot be less or equal to l
return
pivotIdx = random.randint(l, r) # starting pivot index
p = partitionSort(A, l, r, pivotIdx) # ending pivot index after sorting
quickSort(A, l, p - 1)
quickSort(A, p + 1, r)
return A
nums = [100, 3, 9, 1, 0, -10, 87]
print(quickSort(nums))
print(nums)
The error catching code if not (l <= pivotIdx <= r): was never used, but it is nevertheless good to put into place in case of any misuse.
Heapsort 堆排序
Like Selection sort, heapsort partitions its input into a sorted and an unsorted region, and it iteratively shrinks the unsorted region by selecting the largest element from it and inserting it into the sorted region in order.
Unlike selection sort, heapsort does not waste time with a linear-time scan (no pointers) of the unsorted region; rather, heap sort maintains the unsorted region in a max heap data structure to more quickly find the largest element in each step.
In practice on most machines it is slower than a well-implemented quicksort (\(O(n^2)\)), it has better worst-case \(O(n*log(n))\) runtime than quicksort.
Heapsort is an in-place algorithm, but it is not a stable sort.
Max heap data structrue
The heap data structure is a fusion of array and tree.
The max heap data structure is a complete binary tree with the max heap property: At each node, the parent is larger than its chidren.
-
Complete means that at all levels of the tree, except possibly the last one (deepest) are fully filled, and, if the last level of the tree is not complete, the nodes of that level are filled from left to right (left-aligned).
-
$$ \forall \text{ complete binary tree, }\exists\text{ an unique array, and vice versa.}
- \[Left(i)=2*i+1\]
- \[Right(i)=2*i+2\]
- \[Parent(j)=[(i-2)/2]\]
The “parent-children” relationship are defined by the array indices in the above three functions. There is no need of any pointers. There is no need for attributes such as left or right children either.
Implement max heap in Python
Below implementation of the maxHeap object has 3 public functions: push, peek, and pop.
When a maxHeap object is instantiated, the input is transformed into a max heap using the __heapifyUp function.
Both __heapifyUp and __heapifyDown functions are recursive functions. The former compares with “parent” whereas the latter compares with “children”.
The __heapifyUp compares the value at the current index with the value at its parent index, if the value at current index is bigger, then it swaps with “its parent” so that the value at the parent index is bigger.
The __heapifyDown compares the value at the current index with the values at its children indices. If the current value is smaller, then it swaps with “one of its children” so that the value at the child index is smaller, or equivalently the value at the current index is bigger.
.
class maxHeap:
"""
3 private functions:
__swap
__heapifyUp
__heapifyDown
"""
def __init__(self, items= []):
# def __init__(self, items= object):
self.heap = []
for i in items:
self.heap.append(i)
self.__heapifyUp(len(self.heap) - 1)
def push(self, data):
self.heap.append(data)
self.__heapifyUp(len(self.heap) - 1)
def peek(self):
if self.heap[0]:
return self.heap[0]
else:
return False
def pop(self):
if len(self.heap) > 1:
self.__swap(0, len(self.heap) - 1) # swap the max and last element
max = self.heap.pop()
self.__heapifyDown(0) # heapify down the last element that got swapped to the max position
elif len(self.heap) ==1:
max = self.heap.pop()
else:
max = False
return max
def __swap(self, i, j):
self.heap[i], self.heap[j] = self.heap[j], self.heap[i]
def __heapifyUp(self, idx):
parent = (idx - 1) //2 # parent index
if idx <= 0:
return # do nothing if already at the top
elif self.heap[idx] > self.heap[parent]: # if larger than parent, then swap with parent
self.__swap(idx, parent)
self.__heapifyUp(parent) # recurse
def __heapifyDown(self, idx): # heapify down until it is not smaller than its children
left = idx *2 +1 # left child index
right = idx * 2 + 2 # right child index
largest = idx # assume current idx holds the largest
if len(self.heap) > left and self.heap[largest] < self.heap[left]:
largest = left
if len(self.heap) > right and self.heap[largest] < self.heap[right]:
largest = right
if largest != idx:
self.__swap(idx, largest)
self.__heapifyDown(largest) # recurse
Below code instantiates a maxHeap object we defined above, and repeatedly pops the maximum.
A = [5, 7,3,9, 2,8, 1, 4,6,]
maxHeap(A).heap
# [9, 7, 8, 6, 2, 3, 1, 4, 5]
M = maxHeap(A)
print(M.heap)
# [9, 7, 8, 6, 2, 3, 1, 4, 5]
M.push(10)
print(M.heap)
# [10, 9, 8, 6, 7, 3, 1, 4, 5, 2]
from collections import deque
A = deque()
for i in range(len(M.heap)):
A.appendleft(M.pop())
print(A) # sorted
# [10, 9, 8, 6, 7, 3, 1, 4, 5, 2]
# deque([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
The heapq module
We can also use the Python heapq module. The heapq module is part of Python standard library, with well-documented Github source code.
Note that it implements min heap, which means that the root node is the smallest.
The _heapify_max function converts an input to a max heap in-place, in O(n) time.
import heapq
A = [5, 7 ,3, 9, 2, 8, 1, 4, 6]
heapq.heapify(A)
print(A)
# [1, 2, 3, 4, 7, 8, 5, 9, 6]
A = [5, 7 ,3, 9, 2, 8, 1, 4, 6]
heapq._heapify_max(A)
print("max heap", A)
# max heap [9, 7, 8, 6, 2, 3, 1, 4, 5]
print(heapq.heappop(A))
print("after poping max")
print(A)
print(heapq.heappop(A))
print("after poping the second number, we see that it is no longer max heap")
print(A)
print("need to heapify again")
heapq._heapify_max(A)
print("max heap", A)
Bucket sort (Radix sort)
Radix sort is a very fast sorting algorithm for integers. Unlike other sorting methods, it does no comparisons. Digits of integers are slotted into their respective buckets (0, 1, 2, 3, …, 9).
There is no comparions and no if-else branching
Bucket sort is very fast for large quantities of small integers.
Even though the algorithm is for integers only, other objects can be mapped to integers to use this sorting method.
Note the order of operations in i//10**(digit)%10:
The exponential is computed first even if you put a space between like i//10** (digit)%10! To avoid any confusion, it may be better to write i//(10**digit)%10
Also, pay attention to that \(A\) is modified with each loop.
In the first outer loop, \(A\) is the input. In all subsequent outer loops, \(A\) is the flattened buckets \(B\).
def bucketSort(A):
num_digits = len(str(max(A)))
for digit in range(0, num_digits):
B = [[] for i in range(10)] # list of list
for i in A:
num = i//10**(digit)%10
B[num].append(i)
# flatten
A = []
for i in B:
A.extend(i)
return A
A = [4, 9, 100, 1, 598, 0, 8]
print(bucketSort(A))
# [0, 1, 4, 8, 9, 100, 598]