Heapsort 堆排序算法詳解（Java實現)

Heapsort (堆排序)是最經典的排序算法之一，在google或者百度中搜一下可以搜到很多非常詳細的解析。同樣好的排序算法還有quicksort（快速排序）和merge sort（歸並排序），選擇對這個算法進行分析主要是因為它用到了一個非常有意思的算法技巧：數據結構 - 堆。而且堆排其實是一個看起來復雜其實並不復雜的排序算法，個人認為heapsort在機器學習中也有重要作用。這里重新詳解下關於Heapsort的方方面面，也是為了自己鞏固一下這方面知識，有可能和其他的文章有不同的入手點，如有錯誤，還請指出。文中引用的referecne會再結尾標注。

p.s. 個人認為所謂詳解是你在看相關wiki或者算法書看不懂的時候看通俗易懂的解釋，不過最佳方案還是去看教授們的講解，推薦reference[1]中的heapsort章節。

 

以上是廢話，可以不看

 

---

 Section 1 - 簡介

Heapsort是一個comparison-based的排序算法（快排，歸並，插入都是；counting sort不是），也是一種選擇排序算法（selection sort），一個選擇算法（selection algorithm）的定義是找到一個序列的k-th order statistic（統計學中的術語），直白的說就是找到一個list中第k-th小的元素。以上都可以大不用懂，heapsort都理解了回來看一下是這回事就是了。同樣，插值排序也是一種選擇排序算法。

 

Heapsort的時間復雜度在worst-case是\(O(nlgn)\)，average-case是\(O(nlgn)\)；空間復雜度在worst-case是\(O(1)\)，也就是說heapsort可以in-place實現；heapsort不穩定。

 

以下順便附上幾種排序算法的時間復雜度比較（\(\Theta-notation\)比\(O-notation\)更准確的定義了漸進分析（asymptotic analysis）的上下界限，詳細了解可以自行google）：

Table 1 - 四種排序算法的running time比較

Algorithm	Worst-case	Average-case/expected
Insertion sort（插值排序）	\(\Theta (n^2)\)	\(\Theta (n^2)\)
Merge sort（歸並排序）	\(\Theta (nlgn)\)	\(\Theta (nlgn)\)
Heapsort（堆排序）	\(O(nlgn)\)	\(O(nlgn)\)
Quicksort（快速排序）	\(\Theta (n^2)\)	\(\Theta (n^2)\) (expected)

*Additional Part - KNN

heapsort在實踐中的表現經常不如quicksort（盡管quicksort最差表現為 \(\Theta (n^2)\)，但quicksort 99%情況下的runtime complexity為 \(\Theta (nlgn)\)），但heapsort的\(O(nlgn)\)的上限以及固定的空間使用經常被運作在嵌入式系統。在搜索或機器學習中經常也有重要的作用，它可以只返回k個排序需要的值而不管其他元素的值。例如KNN（K-nearest-neighbour）中只需返回K個最小值即可滿足需求而並不用對全局進行排序。當然，也可以使用divide-and-conquer的思想找最大/小的K個值，這是一個題外話，以后有機會做一個專題比較下。

以下程序為一個簡單的在python中調用heapq進行heapsort取得k個最小值，可以大概體現上面所述的特性：

'''
Created On 15-09-2014

@author: Jetpie

'''


import heapq, time
import scipy.spatial.distance as spd
import numpy as np

pool_size = 100000

#generate an 3-d random array of size 10,000
# data = np.array([[2,3,2],[3,2,1],[2,1,3],[2,3,2]])
data = np.random.random_sample((pool_size,3))
#generate a random input
input = np.random.random_sample()
#calculate the distance list
dist_list =  [spd.euclidean(input,datum) for datum in data]

#find k nearest neighbours
k = 10

#use heapsort
start = time.time()
heap_sorted = heapq.nsmallest(k, range(len(dist_list)), key = lambda x: dist_list[x])
print('Elasped time for heapsort to return %s smallest: %s'%(k,(time.time() - start))) 

#find k nearest neighbours
k = 10000

#use heapsort
start = time.time()
heap_sorted = heapq.nsmallest(k, range(len(dist_list)), key = lambda x: dist_list[x])
print('Elasped time for heapsort to return %s smallest: %s'%(k,(time.time() - start)))

 運行結果為:

Elasped time for heapsort to return 10 smallest: 0.0350000858307
Elasped time for heapsort to return 10000 smallest: 0.0899999141693

 

---

 Section 2 - 算法過程理解

2.1 二叉堆

　　　

在“堆排序”中的“堆”通常指“二叉堆（binary heap)”，許多不正規的說法說“二叉堆”其實就是一個完全二叉樹（complete binary tree），這個說法正確但不准確。但在這基礎上理解“二叉堆”就非常的容易了，二叉堆主要滿足以下兩項屬性（properties）：

#1 - Shape Property: 它是一個完全二叉樹。

#2 - Heap Property: 主要分為max-heap property和min-heap property（這就是我以前說過的術語，很重要）

　　|--max-heap property ：對於所有除了根節點（root）的節點 i，\(A[Parent] \geq A[i]\)

　　|--min-heap property ：對於所有除了根節點（root）的節點 i，\(A[Parent] \leq A[i]\)

上圖中的兩個二叉樹結構均是完全二叉樹，但右邊的才是滿足max-heap property的二叉堆。

在以下的描述中，為了方便，我們還是用堆來說heapsort中用到的二叉堆。

 

2.2 一個初步的構想

有了這樣一個看似簡單的結構，我們可以產生以下初步構想來對數組A做排序：

1.將A構建成一個最大堆（符合max-heap property,也就是根節點最大）；

2.取出根節點（how?）；

3.將剩下的數組元素在建成一個最大二叉堆，返回第2步，直到所有元素都被取光。

 

如果已經想到了以上這些，那么就差不多把heapsort完成了，剩下的就是怎么術語以及有邏輯、程序化的表達這個算法了。

 

2.3 有邏輯、程序化的表達

通常，heapsort使用的是最大堆（max-heap）。給一個數組A（我們使用 Java序列[0...n]），我們按順序將它初始化成一個堆：

Input:　

Initialization: 

*堆的根節點（root）為A[0];

對這個堆中index為\(i\)的節點，我們可以得到它的parent, left child and right child，有以下操作：

Parent(i): \(parent(i)\gets A[floor((i-1)/2)]\)

Left(i): \(left(i)\gets A[2*i + 1]\)

Right(i): \(right(i)\gets A[2*i + 2]\)

通過以上操作，我們可以在任意index-\(i\)得到與其相關的其他節點（parent/child)。

在heapsort中，還有三個非常重要的基礎操作（basic procedures):

Max-Heapify(A , i): 維持堆的#2 - Heap Property，別忘了在heapsort中我們指的是max-heap property（min-heap property通常是用來實現priority heap的，我們稍后提及）。

Build-Max-Heap(A): 顧名思義，構建一個最大堆（max-heap）。

Heapsort(A): 在Build-Max-Heap(A)的基礎上實現我們2.2構想中得第2-3步。

其實這三個操作每一個都是后面操作的一部分。

下面我們對這三個非常關鍵的步驟進行詳細的解釋。

 

Max-Heapify(A , i)

+Max-Heapify的輸入是當前的堆A和index-\(i\)，在實際的in-place實現中，往往需要一個heapsize也就是當前在堆中的元素個數。

+Max-Heapify有一個重要的假設：以Left(\(i\))和Right(\(i\))為根節點的subtree都是最大堆（如果樹的知識很好這里就很好理解了，但為什么這么假設呢？在Build-Max-Heap的部分會解釋）。 

+有了以上的輸入以及假設，那么只要對A[i], A[Left(i)]和A[Right(i)]進行比較，那么會產生兩種情況：

　　-第一種，最大值（\(largest\)）是A[i]，那么基於之前的重要假設，以\(i\)為根節點的樹就已經符合#2 - Heap Property了。

　　-第二種，最大值（\(largest\)）是A[Left(i)]或A[Right(i)]，那么交換A[i]與A[\(largest\)]，這樣的結果是以\(largest\)為根節點的subtree有可能打破了#2 - Heap Property，那么對以\(largest\)為根節點的樹進行Max-Heapify(A, largest)的操作。

+以上所述的操作有一個形象的描述叫做A[i] “float down", 使以\(i\)為根節點的樹是符合#2 - Heap Property的，以下的圖例為A[0] ”float down"的過程（注意，以A[1]和A[2]為根節點的樹均是最大堆）。

 　　　　　　

以下附上reference[1]中的Psudocode：

MAX-HEAPIFY（A, i） 
    l = LEFT(i)
    r = RIGHT(i)
    if <= heapsize and A[l] > A[i]
        largest = l
    else largest = i
    if r <= heapsize and A[r] > A[largest]
        largest = r
    if not largest = i
        exchange A[i] with a[largest]
        MAX-HEAPIFY（A, largest） 

 

Build-Max-Heap(A)

先附上reference[1]中的Psudocode（做了部分修改，這樣更明白），因為非常簡單：

BUILD-MAX-HEAP(A)
    heapsize = A.length
    for i = PARENT(A.length-1) downto 0
        MAX-HEAPIFY(A , i)

+Build-Max-Heap首先找到最后一個有子節點的節點 \(i = PARENT(A.length -1)\) 作為初始化（Initialization），因為比 i 大的其他節點都沒有子節點了所以都是最大堆。

+對 i 進行降序loop並對每個 i 都進行Max-Heapify的操作。由於比 i 大的節點都進行過Max-Heapify操作而且 i 的子節點一定比 i 大， 因此符合了Max-Heapify的假設（以Left(\(i\))和Right(\(i\))為根節點的subtree都是最大堆）。

下圖為對我們的輸入進行Build-Max-Heap的過程：

　　　　　　

　　　　

 

Heapsort(A)

到現在為止我們已經完成了2.2中構想的第一步，A[0]也就是root節點是數組中的最大值。如果直接將root節點取出，會破壞堆的結構，heapsort算法使用了一種非常聰明的方法。

+將root節點A[0]和堆中最后一個葉節點（leaf）進行交換，然后取出葉節點。這樣，堆中除了以A[0]為root的樹破壞了#2 - Heap Property，其他subtree仍然是最大堆。只需對A[0]進行Max-Heapify的操作。

+這個過程中將root節點取出的方法也很簡單，只需將\(heapsize\gets heapsize -1\)。

下面是reference[1]中的Psudocode：

HEAPSORT(A):
    BUILD-MAX-HEAP(A)
    for i = A.length downto 1
        exchange A[0] with A[i]
        heapsize =  heapsize -1
        MAX-HEAPIFY(A , 0)

 

到此為止就是整個heapsort算法的流程了。注意，如果你是要閉眼睛也能寫出一個堆排，最好的方法就是理解以上六個重要的操作。 

 

---

 Section 3 - runtime復雜度分析

這一個section，我們對heapsort算法過程中的操作進行復雜度分析。

首先一個總結：

• Max-Heapify ~ \(O(lgn)\)
• Build-Max-Heap ~ \(O(n)\)
• Heapsort ~ \(O(nlgn)\)

然后我們分析一下為什么是這樣的。在以下的分析中，我們所指的所有節點\(i\)都是從1開始的。

Max-Heapify 

這個不難推導，堆中任意節點 i 到葉節點的高度（height）是\(lgn\)。要專業的推導，可以參考使用master theorem。

Build-Max-Heap

在分析heapsort復雜度的時候，最有趣的就是這一步了。

如果堆的大小為\(n\)，那么堆的高度為\(\lfloor lgn\rfloor\)；

對於任意節點\(i\)，\(i\)到葉節點的高度是\(h\)，那么高度為\(h\)的的節點最多有\(\lceil n /2^{h+1}\rceil\)個，下面是一個大概的直觀證明：

　　-首先，一個大小為\(n\)的堆的葉節點（leaf）個數為\(\lceil n/2\rceil\):

　　　　--還記不記得最后一個有子節點的節點parent(length - 1)是第\(\lfloor n/2\rfloor\)（注意這里不是java序號，是第幾個），由此可證葉節點的個數為n - \(\lfloor n/2\rfloor\);

　　-那么如果去掉葉節點，剩下的堆的節點個數為\(n - \lceil n/2\rceil = \lfloor n/2\rfloor\)，這個新樹去掉葉節點后節點個數為\(\lfloor \lfloor n/2\rfloor /2\rfloor\) ;

　　-（這需要好好想一想）以此類推，最后一個樹的葉節點個數即為高度為\(h\)的節點的個數，一定小於\(\lceil (n/2)/2^h\rceil\)，也就是\(\lceil n/2^{h+1}\rceil\)。

對於任意節點\(i\)，\(i\)到葉節點的高度是\(h\)，運行Max-Heapify所需要的時間為\(O(h)\)，上面證明過。

那么Build-Max-Heap的上限時間為（參考reference[1]）：

$\sum_{h=0}^{\lfloor lgn\rfloor } \lceil \frac{n}{2^{h+1}}\rceil O(h) = O\left(n\sum_{h=0}^{\lfloor lgn\rfloor }\frac{h}{2^h}\right)$

根據以下定理：

$\sum_{k=0}^{\infty } kx^k = \frac{x}{(1-x)^2} for \quad |x| < 1$

我們用$x = \frac{1}{2}$替換求和的部分得到：

$\sum_{h=0}^{\infty } \frac{h}{2^h} = \frac{1/2}{(1-1/2)^2} = 2$

綜上所述，我們可以求得：

$O\left(n\sum_{h=0}^{\lfloor lgn\rfloor }\frac{h}{2^h}\right) = O\left(n\sum_{h=0}^{\infty}\frac{h}{2^h}\right) = O(2n) = O(n)$

Heapsort

由於Build-Max-Heap復雜度為$O(n)$，有n-1次調用Max-Heapify（復雜度為$O(lgn)$），所有總的復雜度為$O(nlgn)$

 

到此為止，所有functions的運行復雜度都分析完了，下面的章節就是使用Java的實現了。

 

---

 Section 4 - Java Implementation

這個Section一共有兩個內容，一個簡單的Java實現（只有對key排序功能）和一個Priority Queue。

Parameters & Constructors:

protected double A[];
protected int heapsize;

//constructors
public MaxHeap(){}
public MaxHeap(double A[]){
    buildMaxHeap(A);
}

求parent/left child/right child:

protected int parent(int i) {return (i - 1) / 2;}
protected int left(int i) {return 2 * i + 1;}
protected int right(int i) {return 2 * i + 2;}

保持最大堆特性：

protected void maxHeapify(int i){
    int l = left(i);
    int r = right(i);
    int largest = i;
    if (l <= heapsize - 1 && A[l] > A[i])
        largest = l;
    if (r <= heapsize - 1 && A[r] > A[largest])
        largest = r;
    if (largest != i) {
        double temp = A[i];
        // swap
        A[i] = A[largest];
        A[largest] = temp;
        this.maxHeapify(largest);
    }
}

 構造一個“最大堆”：

public void buildMaxHeap(double [] A){
    this.A = A;
    this.heapsize = A.length;
    
    for (int i = parent(heapsize - 1); i >= 0; i--)
        maxHeapify(i);        
}

對一個array使用heapsort:

public void heapsort(double [] A){
    buildMaxHeap(A);
    
    int step = 1;
    for (int i = A.length - 1; i > 0; i--) {
        double temp = A[i];
        A[i] = A[0];
        A[0] = temp;
        heapsize--;
        System.out.println("Step: " + (step++) + Arrays.toString(A));
        maxHeapify(0);
    }        
}

main函數:

public static void main(String[] args) {
    //a sample input
    double [] A = {3, 7, 2, 11, 3, 4, 9, 2, 18, 0};
    System.out.println("Input: " + Arrays.toString(A));
    MaxHeap maxhp = new MaxHeap();
    maxhp.heapsort(A);
    System.out.println("Output: " + Arrays.toString(A));

}

運行結果:

Input: [3.0, 7.0, 2.0, 11.0, 3.0, 4.0, 9.0, 2.0, 18.0, 0.0]
Step: 1[0.0, 11.0, 9.0, 7.0, 3.0, 4.0, 2.0, 2.0, 3.0, 18.0]
Step: 2[0.0, 7.0, 9.0, 3.0, 3.0, 4.0, 2.0, 2.0, 11.0, 18.0]
Step: 3[2.0, 7.0, 4.0, 3.0, 3.0, 0.0, 2.0, 9.0, 11.0, 18.0]
Step: 4[2.0, 3.0, 4.0, 2.0, 3.0, 0.0, 7.0, 9.0, 11.0, 18.0]
Step: 5[0.0, 3.0, 2.0, 2.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Step: 6[0.0, 3.0, 2.0, 2.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Step: 7[0.0, 2.0, 2.0, 3.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Step: 8[2.0, 0.0, 2.0, 3.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Step: 9[0.0, 2.0, 2.0, 3.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Step: 10[0.0, 2.0, 2.0, 3.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]
Output: [0.0, 2.0, 2.0, 3.0, 3.0, 4.0, 7.0, 9.0, 11.0, 18.0]

 

heapsort在實踐中經常被一個實現的很好的快排打敗，但heap有另外一個重要的應用，就是Priority Queue。這篇文章只做拓展內容提及，簡單得說，一個priority queue就是一組帶key的element，通過key來構造堆結構。通常，priority queue使用的是min-heap，例如按時間順序處理某些應用中的objects。

為了方便，我用Inheritance實現一個priority queue：

package heapsort;

import java.util.Arrays;

public class PriorityQueue extends MaxHeap{

    public PriorityQueue(){super();}
    public PriorityQueue(double [] A){super(A);}
    
    public double maximum(){
        return A[0];
    }
    
    public double extractMax(){
        if(heapsize<1)
            System.err.println("no element in the heap");
        double max = A[0];
        A[0] = A[heapsize-1];
        heapsize--;
        this.maxHeapify(0);
        return max;
    }
    
    public void increaseKey(int i,double key){
        if(key < A[i])
            System.err.println("new key should be greater than old one");
        
        A[i] = key;
        while(i>0 && A[parent(i)] <A[i]){
            double temp = A[i];
            A[i] = A[parent(i)];
            A[parent(i)] = temp;
            i = parent(i);
        }
    }
    
    public void insert(double key){
        heapsize++;
        A[heapsize - 1] = Double.MIN_VALUE;
        increaseKey(heapsize - 1, key);
    }
    
    public static void main(String[] args) {
        //a sample input
        double [] A = {3, 7, 2, 11, 3, 4, 9, 2, 18, 0};
        System.out.println("Input: " + Arrays.toString(A));
        PriorityQueue pq = new PriorityQueue();
        pq.buildMaxHeap(A);
        System.out.println("Output: " + Arrays.toString(A));
        pq.increaseKey(2, 100);
        System.out.println("Output: " + Arrays.toString(A));
        System.out.println("maximum extracted: " + pq.extractMax());
        pq.insert(33);
        System.out.println("Output: " + Arrays.toString(A));    

    }
}

 運行結果：

Input: [3.0, 7.0, 2.0, 11.0, 3.0, 4.0, 9.0, 2.0, 18.0, 0.0]
Output: [18.0, 11.0, 9.0, 7.0, 3.0, 4.0, 2.0, 2.0, 3.0, 0.0]
Output: [100.0, 11.0, 18.0, 7.0, 3.0, 4.0, 2.0, 2.0, 3.0, 0.0]
maximum extracted: 100.0
Output: [33.0, 18.0, 4.0, 7.0, 11.0, 0.0, 2.0, 2.0, 3.0, 3.0]

 

 

---

Section 5 - 小結

首先要說本文全部是原創，如需要使用，只需要引用一下並不需要通知我。

寫到最后發現有很多寫得很冗余，也有啰嗦的地方。感覺表達出來對自己的知識鞏固很有幫助。Heapsort真是一個非常有意思的排序方法，是一個通用而不算復雜的算法，這是決定開始寫blogs后的第一篇文章，一定有很多不足，歡迎討論！之后打算寫一些機器學習和計算機視覺方面的來拋磚引玉。希望通過博客園這個平台可以交到更多有鑽研精神的朋友。

 

Bibliography

[1] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein.
Introduction to Algorithms. MIT Press, 3rd edition, 2009.

[2] Wikipedia. Heapsort — wikipedia, the free encyclopedia, 2014. [Online; accessed
15-September-2014].