SICP:符號求導、集合表示和Huffman樹(Python實現)

orion發表於2023-01-04

緒論

到目前為止,我們已經使用過的所有複合資料,最終都是從數值出發構造起來的(比如我們在上一篇部落格《SICP 2.2: 層次性資料和閉包性質(Python實現)》所介紹的連結串列和樹就基於數來進行層次化構造)。在這一節裡,我們要擴充所用語言的表達能力,引進將任意符號作為資料的功能。

2.3.1 Scheme語言中的引號

在《SICP》原書採用的Scheme語言(Lisp的一種方言)中,要想表示諸如(a b c d)這種包含著符號的表非常簡單,可以直接對資料物件加引號。例如在Scheme語言中,若我們定義變數

(define a 1)
(define b 2)

則我們對(list a b)(list 'a 'b)(list 'a b)分別進行求值的結果如下

(list a b)
(1 2)

(list 'a 'b)
(a b)

(list 'a b)
(a 2)

可見直接採用(list a b)構造出的是ab的值的表,而不是這兩個符號本身的表。若想要將符號作為資料物件看待,而不是作為應該求值的表示式,可以在被引物件的前面放一個引號,如'a'b

這裡的引號非常有趣,因為在自然語言環境中這種情況也很常見,在那裡的單詞和句子可能看作語義實體,也可以看作是字元的序列(語法實體)。在自然語言中,我們常常用引號表明一個詞或者句子作為文字看待,將它們直接作為字元的序列。例如我們對某人說“大聲說你的名字”,此時希望聽到的是那個人的名字。如果說“大聲說‘你的名字’”,此時希望聽到的就是片語“你的名字”。請注意,我們在這裡不得不用巢狀的引號去描述別人應該說的東西。

然而,在Python語言中,卻不能如此優雅地像資料一樣運算子號/程式碼本身了,這也是Scheme/Lisp本身的魅力之所在,即所謂萬物皆資料(參看知乎回答《Lisp 值得去學習嗎?》)。

由於我們採用Python語言進行實現,故接下我們統一採用字串的形式來表示符號

2.3.2 例項:符號求導

為了闡釋符號操作的情況,並進一步闡釋資料抽象的思想,現在考慮設計一個執行代數表示式的符號求導的過程。我們希望該過程以一個代數表示式和一個變數作為引數,返回這個表示式相對於該變數的導數。例如,如果送給這個過程的引數是\(ax^2+bx+c\)\(x\),它應該返回\(2ax+b\)。符號求導對於Lisp有著特殊的歷史意義,它正是推動人們去為符號操作開發計算機語言的重要例項之一。進一步說,它是人們為符號數學開發強有力系統的開端,今天已經有越來越多的應用數學家和物理學家們正在使用這類系統。

為了開發出一個符號計算程式,我們將按照2.1.1節開發有理數系統那樣,採用同樣的資料抽象策略。也就是說,首先定義一個求導演算法,令它在一些抽象物件上操作,例如“和”、“乘積”和“變數”,並不考慮這些物件實際上如何表示,以後才去關心具體表示的問題。

對抽象資料的求導程式
為了使有關的討論簡單化,我們在這裡考慮一個非常簡單的符號求導程式,它處理的表示式都是由對於兩個引數的加和乘運算構造起來的。對於這種表示式求導的工作可以透過下面幾條規約規則完成:

\[\begin{aligned} & \frac{\mathrm{d} c}{\mathrm{~d} x}=0\quad \text{當}c\text{是一個常量, 或者一個與}x{不同的變數}\\ & \frac{\mathrm{d} x}{\mathrm{~d} x}=1 \\ & \frac{\mathrm{d}(u+v)}{\mathrm{d} x}=\frac{\mathrm{d} u}{\mathrm{~d} x}+\frac{\mathrm{d} v}{\mathrm{~d} x} \\ & \frac{\mathrm{d}(u v)}{\mathrm{d} x}=u\left(\frac{\mathrm{d} v}{\mathrm{~d} x}\right)+v\left(\frac{\mathrm{d} u}{\mathrm{~d} x}\right) \end{aligned} \]

可以看到,這裡的最後兩條規則具有遞迴的性質,也就是說,想要得到一個和式的導數,我們首先要找到其中各個項的導數,而後將它們相加。這裡的每個項又可能是需要進一步分解的表示式。透過這種分解,最終將產生出常量和變數,它們的導數就是0或者1.

為了能在一個過程中體現這些規則,我麼就像在前面設計有理數的實現時所做的那樣,採用同樣的資料抽象策略。給定一個代數表示式,我們可以判斷出某個表示式是否為一個和式、乘式、常量或變數,也能夠提取出表示式裡的各個部分。對於一個和式,我們可能希望取得其被加項(第一個項)和加項(第二個項)。我們假定現在已經有了一些過程,它們實現了下述的建構函式、選擇函式和謂詞:

is_number(e)  # e是數值嗎?
is_variable(e)  # e是變數嗎?
is_same_variable(v1, v2) # v1和v2是同一個變數嗎
is_sum(e) # e是和式嗎?
addend(e) # 取e中的被加數
augend(e) # 取e中的加數
make_sum(a1, a2) # 構造起a1與a2的和式

is_product(e) # e是乘式嗎?
multiplier(e) # 取e中的被乘數
multiplicand(e) # 取e中的乘數
make_product(m1, m2) # 構造起m1與m2的乘式

利用這些過程,我們就可以將各種求導規則用下面的過程表示出來了:

# 引數exp為表示式,var為欲求導的自變數
def deriv(exp, var):
    if is_number(exp):
        return "0"
    elif is_variable(exp):
        if is_same_variable(exp, var):
            return "1"
        else:
            return "0"
    elif is_sum(exp):
        return make_sum(deriv(addend(exp), var), \
                        deriv(augend(exp), var))
    elif is_product(exp):
        return make_sum(\
            make_product(multiplier(exp), \
                        deriv(multiplicand(exp), var)), \
            make_product(deriv(multiplier(exp), var), \
                        multiplicand(exp))
            )
    else:
        raise ValueError("unknown expression type") 

過程deriv裡包含了一個完整的求導演算法。因為它是基於抽象資料表述的,因此,無論我們如何選擇代數表示式exp的具體表示,只要正確地設計了我們前面列出的選擇函式、謂詞函式和建構函式,這個過程都可以工作。下面我們來看如何對這些函式進行設計。

函式設計

在對函式進行設計之前,我們需要先確定代數表示式的表示。我們固然可以用字串來直接表達代數表示式,如將表示式\(ax+b\)表示為字串"a*x+b"。然而,這種表示方式體現不出表示式樹的層次,不利於我們後面對錶達式的操作。回想一下,我們在部落格《Python:實現簡單的遞迴下降Parser》中介紹瞭如何將表示式解析為表示式樹,也即將\(ax+b\)表示為("+", ("*", "a", "x"), "b"),這樣我們就能夠對上述的函式進行進一步設計了。注意,在Lisp中組合式本身就採用帶括號的字首形式來表示,因此若用Lisp語言實現則不需要事先寫Parser(其實推廣一步說,Lisp的整個程式都是一個由巢狀括號組成的表示式樹,也就是說Lisp語言的的直譯器也是不需要寫Parser的)。

  • 判斷是否為數值:
def is_number(x):
    return isinstance(x, str) and x.isnumeric()
  • 判斷是否為變數:
def is_variable(x):
    return isinstance(x, str) and x.isalpha()
  • 判斷兩個變數是否相同:
def is_same_variable(v1, v2):
    return is_variable(v1) and is_variable(v2) and v1 == v2
  • 和式與乘式都構造為表
# 構造起a1和a2的和式
def make_sum(a1, a2):
    return ("+", a1, a2)

# 構造起m1和m2的乘式
def make_product(m1, m2):
    return ("*", m1, m2)
  • 和式就是第0個元素(下標從0開始)為符號為+的元組
def is_sum(x):
    return isinstance(x, tuple) and x[0] == "+"
  • 被加數是表示和式的表裡的第1個元素(下標從0開始):
def addend(s):
    return s[1]
  • 加數是表達和式元組裡的第2個元素(下標從0開始):
def augend(s):
    return s[2]
  • 乘式就是第0個元素(下標從0開始)為符號為*的元組
def is_product(x):
    return isinstance(x, tuple) and x[0] == "*"
  • 被乘數是表示乘式的元組裡的第1個元素(下標從0開始)
def multiplier(p):
    return p[1]
  • 乘數是表示乘式的元組裡的第2個元素(下標從0開始)
def multiplicand(p):
    return p[2]

這樣,為了得到一個能夠工作的符號求導程式,我們只需要將這些過程與deriv裝在一起,現在讓我們看幾個表現這一程式的行為的例項:

# 輸入為元組,預設已經完成了表示式解析
print(deriv(("+", "x", "3"), "x"))
# ('+', '1', '0')

print(deriv(("*", "x", "y"), "x")) 
# ('+', ('*', 'x', '0'), ('*', '1', 'y'))

print(deriv(("*", ("*", "x", "y"), ("+", "x", "3")), "x"))
# ('+', ('*', ('*', 'x', 'y'), ('+', '1', '0')), ('*', ('+', ('*', 'x', '0'), ('*', '1', 'y')), ('+', 'x', '3')))

程式產生出的這些結果是對的,但是它們沒有經過化簡。我們確實有:

\[\frac{d(x y)}{d x}=x \cdot 0+1 \cdot y \]

當然,我們也可能希望這一程式能夠知道\(x\cdot 0 = 0, 1\cdot y = y\)以及\(0+y=y\)。因此,第二個例子的結果就應該是簡單的\(y\)。正如上面的第三個例子所顯示的,當表示式變得更加複雜時,這一情況也可能變成嚴重的問題。

現在所面臨的困難很像我們在做有理數時所遇到的問題:希望將結果化簡到最簡單的形式。為了完成有理數的化簡,我們只需要修改建構函式和選擇函式的實現。這裡也可以採取同樣的策略。我們在這裡也完全不必修改deriv,只需要修改make_sum,使得當兩個求和物件都是數時,make_sum求出它們的和返回。此外,如果其中的一個求和物件是0,那麼make_sum就直接返回另一個物件。

def make_sum(a1, a2):
    if eq_number(a1, "0"):
        return a2
    elif eq_number(a2, "0"):
        return a1
    elif is_number(a1) and is_number(a2):
        return str(int(a1) + int(a2))
    else:
        return ("+", a1, a2)

在這個實現裡用到了過程eq_number,它檢查某個表示式是否等於一個給點的數。

def eq_number(exp, num):
    return is_number(exp) and exp == num

與此類似,我們也需要修改make-product,設法引進下面的規則:0與任何東西的乘積都是0,1與任何東西的乘積總是那個東西:

def make_product(m1, m2):
    if eq_number(m1, "0") or eq_number(m2, "0"):
        return "0"
    elif eq_number(m1, "1"):
        return m2
    elif eq_number(m2, "1"):
        return m1
    elif is_number(m1) and is_number(m2):
        return str(int(m1) * int(m2))
    else:
        return ("*", m1, m2)

下面是這一新過程版本對前面三個例子的結果:

print(deriv(("+", "x", "3"), "x")) # 1

print(deriv(("*", "x", "y"), "x")) # y

print(deriv(("*", ("*", "x", "y"), ("+", "x", "3")), "x"))
# ('+', ('*', 'x', 'y'), ('*', 'y', ('+', 'x', '3')))
# 即對(x * y) * (x + 3)求導後得到  (x * y) + y * ( x + 3)

顯然情況已經大大改觀。但是,第三個例子還是說明,想要做出一個程式,使它能將表示式做成我們都能同意的“最簡單”形式,前面還有很長的路要走。代數化簡是一個非常複雜的問題,除了其他各種因素之外,還有另外一個根本性問題:對於某種用途的最簡形式,對於另一用途可能就不是最簡形式。

2.3.3 例項:集合的表示

用元組來表示代數表示式是直截了當的,現在我們要轉到集合的表示問題,此時表示方式的選擇就不那麼顯然了。實際上,在這裡存在幾種選擇,而且它們之間在幾個方面存在明顯的不同。

非形式化地說,一個集合就是一些不同物件的彙集。要給出精確定義,我們可以利用資料抽象的方法,也就是說,用一組可以作用於“集合”的操作來定義它們。這些操作是union_setintersection_setis_element_of_setadjoin_setis_element_of_set是一個謂詞,用於確定某個給定元素是不是某個給定集合的成員。adjoin-set以一個物件和一個集合為引數,返回一個集合,其中包含了原集合的所有元素與新加入的元素。union-set為計算兩個集合的並集,intersection-set為計算兩個集合的交集。

集合作為未排序的表
集合的一種表示形式就是用其元素的表,其中任何元素的出現都不超過一次。這樣,空集就用空表來表示。對於這種表示形式,is_element_of_set定義如下:

def is_element_of_set(x, set):
    if not set:
        return False
    elif x == set[0]:
        return True
    else:
        return is_element_of_set(x, set[1: ]) 

利用它就能寫出adjoin-set。如果要加入的物件已經在相應集合裡,那麼就返回那個集合;否則就將這一物件加入表示集合的表裡:

def adjoin_set(x, set):
    if is_element_of_set(x, set):
        return set
    else:
        return [x] + set

實現intersection_set時可以採用遞迴策略:如果我們已知如何做出set2set1[1: ]的交集,那麼就只需要確定是否將set1[0]包含到結果之中了,而這依賴於set1[0]是否也在set2中。下面是這樣寫出的過程:

def intersection_set(set1, set2):
    if not set1 or not set2:
        return []
    elif is_element_of_set(set1[0], set2):
        return [set1[0]] + intersection_set(set1[1: ], set2)
    else:
        return intersection_set(set1[1: ], set2)

在設計一種表示形式時,有一件必須關注的事情是效率問題。為考慮這一問題,就需要考慮上面定義的各集合操作所需要的工作步數。因為它們都使用了is_element_of_set,這一操作的速度對整個集合的實現效率將有重大影響。在上面這個實現裡,如果集合有\(n\)個元素,is_element_of_set就可能需要\(n\)步才能完成。這樣,這一操作所需的步數將以\(\Theta(n)\)的速度增長。adjoin_set使用了這個操作,因此它所需的步數也已\(\Theta(n)\)的速度增長。而對於intersection_set,它需要對set_1的每個元素做一次is_element_of_set檢查,因此所需步數在兩個集合大小都為\(n\)時將以\(\Theta(n^2)\)增長。union_set的情況也是如此。

集合作為排序的表

加速集合操作的一種方式是改變表示方式,使集合元素在表中按照上升序排列。為此,我們就需要有某種方式來比較兩個元素,以便確定哪個元素更大。如按照字典序做符號的比較,或者同意採用某種方式為每個物件關聯一個唯一的數,在比較元素的時候比較該數即可。為了簡化討論,我們僅考慮元素是數值的情況。下面將數的集合表示為元素按照升序排列的表。在前面第一種表示方式下,集合\(\{1, 3, 6, 10\}\)的元素在相應的表裡可以任意排列,而在現在新的表示方式中,我們就只允許用表[1, 3, 6, 10]

從操作is_element_of_set可以看到採用有序表的一個優勢:為了檢查一個項的存在性,現在就不必掃描整個表了。如果檢查中遇到的某個元素大於當時要找的東西,那麼就可以斷定這個東西根本不在表裡:

def is_element_of_set(x, set):
    if not set:
        return False
    elif x == set[0]:
        return True
    elif x < set[0]:
        return False
    else:
        return is_element_of_set(x, set[1: ]) 

這樣能節約多少步數呢?在最壞情況下,我們要找的專案可能是集合中的最大元素,此時所需步數與採用未排序的表示時一樣。不過在平均情況下,我們可以期望需要檢查表中的一半元素,所需步數大約是\(n/2\),這仍然是\(\Theta(n)\)的增長速度,但與前一實現相比,平均來說節約了大約一半的步數(注意,原書此處有誤,前面說未排序表需要檢查整個表,考慮的只是一種特殊情況:查詢沒有出現在表裡的元素。如果查詢的是表裡存在的元素,即使採用未排序的表,平均查詢長度也是表元素的一半)。

操作intersection_set的加速情況更使人印象深刻。在未排序的表示方式裡,這一操作需要\(\Theta(n)^2\)的步數,因為對set1中的每個元素,我們都需要對set2做一次完全的掃描。對於排序表示則可以有一種更聰明的方法。我們在開始時比較兩個集合的起始元素,例如x1x2。如果x1等於x2。那麼這樣就得到了交集的一個元素,而交集的其他元素就是這兩個集合[1:]的交集。如果x1小於x2,由於x2是集合set2的最小元素,我們立即可以斷定x1不會出現在集合set2的任何地方,因此它不應該在交集裡。這樣,兩集合的交集就等於集合set2set1[1:]的交集。與此類似,如果x2小於x1,兩集合的交集就等於集合set1set2[1:]的交集。下面是按這種方式寫出的過程:

def intersection_set(set1, set2):
    if not set1 or not set2:
        return []
    else:
        x1, x2 = set1[0], set2[0]
        if x1 == x2:
            return [x1] + intersection_set(set1[1:], set2[1: ])
        elif x1 < x2:
            return intersection_set(set1[1: ], set2)
        elif x2 < x1:
            return intersection_set(set1, set2[1: ])

現估計這一過程所需的步數。注意,在每個步驟中,我們都將求交集問題歸結到更小集合的交集計算問題——去掉了set1set2之一或者是兩者的第一個元素。這樣,所需步數至多等於set1set2的大小之和,而不像在未排序表示中它們的乘積。這也就是\(\Theta(n)\)的增長速度,而不是\(\Theta(n^2)\)——這一加速非常明顯,即使對中等大小的集合也是如此。

集合作為二叉樹

如果將集合元素安排成一棵樹的形式,我們還可以得到比排序表表示更好的結果。樹中每個節點儲存集合中的一個元素,稱為該節點的“資料項”,它還連結到另外的兩個節點(可能為空)。其中“左邊”的連結所指向的所有元素均小於本節點的元素,而“右邊”連結到的元素都大於本節點裡的元素。如下圖所示是一棵表示集合的樹(同一個集合表示為樹可以有多種不同的方式)。

SICP:符號求導、集合表示和Huffman樹(Python實現)

樹方法的優點在於,對於檢查某個數\(x\)是否在集合裡這個任務,我們可以用x與樹頂點的資料項進行比較。如果x小於它,我麼就知道只需要搜尋左子樹;如果x比較大,那麼就只需要搜尋右子樹。在這樣做時,如果該樹是“平衡的”,也即每棵子樹大約是整個樹的一半大,那麼經過這樣一步我們就將搜尋規模為\(n\)的樹的問題,歸約為搜尋規模為\(n/2\)的樹的問題。由於每次經過這個步驟能夠使樹的大小減小一半,我們可以期望搜尋規模為\(n\)的樹的計算步數以\(\Theta (\log n)\)速度增長。在集合很大時,相對於原來的表示,現在的操作速度將明顯快很多。

我們可以用表來表示樹,講節點表示為三個元素的表:本節點中的資料項,其左子樹和右子樹。以空表作為左子樹或者右子樹,就表示沒有子樹連線在那裡。我們可以用下面過程描述這種表示:

def entry(tree):
    return tree[0]

def left_branch(tree):
    return tree[1]

def right_branch(tree):
    return tree[2]

def make_tree(entry, left, right):
    return [entry, left, right]

現在,我們就可以採用上面描述的方式實現過程is_element_of_set了:

def is_element_of_set(x, set):
    if not set:
        return False
    elif entry(set) == x:
        return True
    elif x < entry(set):
        return is_element_of_set(x, left_branch(set))
    elif x > entry(set):
        return is_element_of_set(x, right_branch(set))

我們對is_element_of_set操作的測試結果如下:

set = [7, [3, [1, [], []], [5, [], []]], [9, [], [11, [], []]]]
print(is_element_of_set(5, set)) # True
print(is_element_of_set(99, set)) # False

像集合裡插入一個項的實現方式與此類似,也需要\(\Theta(\log n)\)的步數。為了加入元素x,我們需要將x與節點資料項比較,以便確定x應該加入右子樹還是左子樹中。下面是這個過程:

def adjoin_set(x, set):
    if not set:
        return make_tree(x, [], [])
    elif x == entry(set):
        return set
    elif x < entry(set):
        return make_tree(entry(set), \
                        adjoin_set(x, left_branch(set)), \
                        right_branch(set))
    elif x > entry(set):
        return make_tree(entry(set), \
                        left_branch(set), \
                        adjoin_set(x, right_branch(set)))

我們對adjoin_set操作的測試結果如下:

set = [7, [3, [1, [], []], [5, [], []]], [9, [], [11, [], []]]]
print(adjoin_set(6, set))
# [7, [3, [1, [], []], [5, [], [6, [], []]]], [9, [], [11, [], []]]]

注意,我們在上面斷言,搜尋樹的操作在對數步數中完成,這實際上依賴於樹“平衡”的假設。也就是說,每個樹的左右子樹中的節點大致上一樣多,因此每棵子樹中包含的節點大約就是其父的一半。但是我們怎樣才能確保構造出的樹是平衡的呢?即使是從一顆平衡的樹開始工作,採用adjoin_set加入元素也可能產生出不平衡的的結果。因為新加入的元素的位置依賴於它與當時已經在樹中的那些項比較的情況。我們可以期望,如果“隨機地”將元素加入樹中,平均而言將會使樹趨於平衡。但在這裡並沒有任何保證。例如,如果我們從空集出發,順序將數值17加入其中,我們就會得到如下圖所示的高度不平衡的樹。

SICP:符號求導、集合表示和Huffman樹(Python實現)

在這個樹裡,所有的左子樹都為空,所以它與簡單排序表相比一點優勢也沒有。解決這個問題的一種方式是定義一個操作,它可以將任意的樹變換為一棵具有同樣元素的平衡樹。在每執行過幾次adjoin_set操作之後,我們就可以透過執行它來保持樹的平衡。當然,解決這一問題的方法還有很多,大部分這類方法都涉及到設計一種新的資料結構,設法使這種資料結構上的搜尋和插入操作都在\(\Theta(\log n)\)步數內完成(這種結構的例子如B樹和紅黑樹,參見《演算法導論》)。

集合與資訊檢索
我們考察了用表表示集合的各種選擇,並看到了資料物件表示的選擇可能如何深刻地影響到使用資料的程式的效能。關注集合的另一個原因是,這裡所討論的技術在設計資訊檢索的應用中將會一次又一次地出現。

現在考慮一個包含大量獨立記錄的資料庫,我們可以將這個資料庫表示為一個記錄的集合,將每個記錄中的一部分當做標識key(鍵值)。為了根據給定鍵值確定相關記錄的位置,我們用一個過程lookup,它以一個鍵值和資料庫為引數,返回具有這個鍵值的記錄,或者在找不到相應記錄時報告失敗。lookup的實現方式幾乎與is_element_of_set一模一樣,如果記錄的集合被表示為未排序的表,我們就可以用:

def lookup(given_key, set_of_records):
    if not set_of_records:
        return False
    elif given_key == key(set_of_records[0]):  
        return set_of_records[0]
    else:
        return lookup(given_key, set_of_records[1: ])

不言而喻,還有比未排序的表更好的表示大集合的方法。常常需要“隨機訪問”其中記錄的資訊檢索系統通常用某種基於樹的方式實現,例如前面討論過的二叉樹。

2.3.4 例項:Huffman編碼樹

本節將給出一個實際使用表結構和資料抽象去操作集合與樹的例子。這一應用是想確定一些用0和1(二進位制位)的序列表示資料的方法。舉例說,用於在計算機裡表示文字的ASCII標準編碼將每個字元表示為一個包含\(7\)個二進位制位的序列,採用\(7\)個二進位制位能夠區分\(2^7\)種不同的情況,即\(128\)個可能不同的字元。一般而言,如果我們需要區分\(n\)個不同字元,那麼就需要為每個字元使用\(\log_2n\)個二進位制位。假設我們的所有資訊都是用A、B、C、D、E、F、G和H這樣8個字元構成的,那麼就可以選擇為每個字元用3個二進位制位,例如:

\[\begin{aligned} \text{A 000} \quad \text{C 010} \quad \text{E 100} \quad \text{G 110} \\ \text{B 001} \quad \text{D 011} \quad \text{F 101} \quad \text{H 111} \end{aligned} \]

採用這種編碼方式,訊息

\[\text{BACADAEAFABBAAAGAH} \]

將編碼為54個二進位制位

\[\text{001000010000011000100000101000001001000000000110000111} \]

像ASCII碼和上面A到H這樣的方式稱為定長編碼,因為它們採用同樣數目的二進位制位表示訊息中的每一個字元。變長編碼方式就是用不同數目的二進位制位表示不同的字元,這種方式有時也可能有些優勢。舉例說,莫爾斯電報碼對於字母表中各個字母就沒有采用同樣數目的點和劃,特別是最常見的字母E只用一個點表示。一般而言,如果在我們的訊息裡,某些符號出現得頻繁,而另一些卻很少見,那麼如果為這些頻繁出現的字元指定較短的碼字,我們就可能更有效地完成資料的編碼(對於同樣訊息使用更少的二進位制位)。請考慮下面對於字母A到H的另一種編碼:

\[\begin{aligned} \text{A 0} \quad \text{C 1010} \quad \text{E 1100} \quad \text{G 1110} \\ \text{B 100} \quad \text{D 1011} \quad \text{F 1101} \quad \text{H 1111} \end{aligned} \]

採用這種編碼方式,上面的同樣資訊將編碼為如下的串:

\[\text{100010100101101100011010100100000111001111} \]

這個串種只包含42個二進位制位,也就是說,與上面的定長編碼相比,現在的這種方式節約了超過\(20\%\)的空間。

採用變長編碼有一個困難,那就是在讀0/1序列的過程中確定何時到達了一個字元的結束。莫爾斯碼解決這一問題的方式是在每個字母的電劃序列之後用一個特殊的分隔符(它用的是一個間歇)。另一種解決方式是以某種方式設計編碼,使得其中每個字元的完整編碼都不是另一個字元編碼的開始一段(或稱字首)。這樣的編碼稱為字首碼。在上面的例子中,A編碼為0而B編碼為100,沒有其他字元的編碼由0或者100開始。

一般而言,如果能夠透過變長字首碼去利用被編碼訊息中符號出現的相對頻度,那麼就能夠明顯地節約空間。完成這件事情的一種特定方式稱為Huffman編碼。一個Huffman編碼可以表示為一顆二叉樹,如下就是上面給出的A到H編碼所對應的Huffman樹:

SICP:符號求導、集合表示和Huffman樹(Python實現)

可以看到,位於樹葉的每一個符號被賦予一個相對權重(也就是它的相對頻度),非葉節點所包含的權重是位於它之下的所有葉節點的權重之和。這種權重只在構造樹的過程中使用,而在編碼和解碼的過程中並不使用。

給定了一顆Huffman樹,要找出任一符號的編碼,我們只需從樹根開始往下,知道找到存有這一符號的樹葉為止,每次向左給編碼加一個0,向右加一個1。比如從上圖中的樹根開始,到達D的葉節點的方式是走右-左-右-右,因此其程式碼為1011。

對應地,解碼時,我們也從樹根開始,透過位序列中的0或1確定是移向左分支還是右分支。當我們到達一個葉節點時,就生成出了訊息中的一個符號。比如如果給我們的是上面的樹和序列10001010,則操作序列為右-左-左找到B,然後再此從根開始左找到A,然後再此從根開始右-左-右-左找到C。這樣整個訊息也就是ABC。

生成Huffman樹
我們上面只說了編碼和解碼的過程,那麼給定了符號的“字母表”和它們的相對頻度,我們怎麼才能構造出“最好的編碼”,也即使訊息編碼的位數達到最少呢?Huffman給出了完成這件事的第一個演算法, 並且證明了若根據符號所出現的相對頻度來建樹, 這樣產出的編碼確實是最好的變長編碼。我們這裡略過Huffman編碼最優性質的證明,將直接展示如何構造Huffman樹。

生成Huffman樹的演算法實際上十分簡單,其想法就是設法安排這棵樹,使得那些帶有最低頻度的符號出現在離樹根最遠的地方。這一過程從葉節點的集合開始,找出兩個具有最低權重的葉進行歸併,產生出一個以這兩個節點為左右分支的節點,新節點的權重就是那兩個點的圈子之和。然後從原集合裡刪除前面的兩個葉節點,利用新節點替代它們。之後以此類推繼續這一過程。當集合中只剩下一個節點時,這一過程終止,而這個節點就是樹根。下面顯示的是上圖中的Huffman樹的生成過程:

SICP:符號求導、集合表示和Huffman樹(Python實現)

這一演算法並不總能描述一棵唯一的樹,因為選出的最小權重節點有可能不唯一。還有,在做歸併時,兩個節點的順序也是任一的,也就是說隨便哪個都可以作為左分支或者右分支。

Huffman樹的表示
在下面的練習中,我們將要做出一個使用Huffman樹完成訊息編碼和編碼、並能根據上面給出的梗概生成Huffman樹的系統。開始還是討論這種樹的表示。

將一棵樹的樹葉表示為包含符號leaf、葉中符號和權重的表:

def make_leaf(symbol, weight):
    return ["leaf", symbol, weight]

def is_leaf(object):
    return object[0] == "leaf"

def symbol_leaf(x):
    return x[1]

def weight_leaf(x):
    return x[2]

在歸併兩個節點做出一棵樹時,樹的權重也就是,樹的權重也就是這兩個節點的權重之和,其符號集就是兩個節點的符號集的並集。這裡符號集我們直接用Python的字串來表示,直接相加即可:

def make_code_tree(left, right):
    return [left, \
            right, \
            symbols(left) + symbols(right), \
            weight(left) + weight(right)]

如果以這種方式構造,我們就需要採用下面的選擇函式:

def left_branch(tree):
    return tree[0]

def right_branch(tree):
    return tree[1]

def symbols(tree):
    if is_leaf(tree):
        return symbol_leaf(tree)
    else:
        return tree[2]
    
def weight(tree):
    if is_leaf(tree):
        return weight_leaf(tree)
    else:
        return tree[3]

在對樹葉或者一般樹呼叫過程symbolsweight時,它們需要做的事情有一點不同。這些不過是通用型過程(可以處理多於一種資料的過程)的簡單示例,有關這方面的情況我們在2.4節和2.4節將有更多討論。

解碼過程
下面的過程實現解碼演算法,它以一個0/1的Python字串和一棵Huffman樹為引數:

def decode(bits, tree):
    def decode_1(bits, current_branch):
        if not bits:
            return ""
        else:
            next_branch = choose_branch(bits[0], current_branch)
            if is_leaf(next_branch):
                return symbol_leaf(next_branch) + decode_1(bits[1: ], tree)
            else:
                return decode_1(bits[1: ], next_branch)                
    return decode_1(bits, tree)

def choose_branch(bit, branch):
    if bit == "0":
        return left_branch(branch)
    elif bit == "1":
        return right_branch(branch)
    else:
        raise ValueError("error: bad bit -- CHOOSE-BRANCH %s" % bit)

下面我們進行一下測試。定義一棵如下所示的編碼樹:

sample_tree = make_code_tree(make_leaf("A", 4),  
                             make_code_tree(
                                 make_leaf("B", 2),
                                 make_code_tree(make_leaf("D", 1), make_leaf("C", 1))
                             ))
print(sample_tree)
# [['leaf', 'A', 4], [['leaf', 'B', 2], [['leaf', 'D', 1], ['leaf', 'C', 1], 'DC', 2], 'BDC', 4], 'ABDC', 8]

這棵樹實際上如下圖所示:

SICP:符號求導、集合表示和Huffman樹(Python實現)

再定義如下的樣例訊息:

sample_message = "0110010101110"

我們用decode完成該訊息的解碼,得到解碼結果:

decoded_message = decode(sample_message, sample_tree)
print(decoded_message)

程式列印輸出的編碼結果為:

ADABBCA

編碼過程
下面的過程實現編碼演算法,它以一個字母組成的Python字串和一棵Huffman樹為引數:

def encode(message, tree):
    if not message:
        return ""
    else:
       return encode_symbol(message[0], tree) + encode(message[1:], tree)
    
def encode_symbol(symbol, current_branch):
    if is_leaf(current_branch):
        return ""
    else:
        left, right = left_branch(current_branch), right_branch(current_branch)
        if symbol in symbols(left):
            return "0" + encode_symbol(symbol, left)
        elif symbol in symbols(right):
            return "1" + encode_symbol(symbol, right)
        else:
            raise ValueError("error: bad symbol -- CHOOSE-BRANCH %s" % symbol)

編碼樹仍採用我們上面解碼部分所定義的sample_tree,然後定義如下的樣例訊息:

sample_message = "ADABBCA"

我們用encode完成該訊息的解碼,得到解碼結果:

encoded_message = encode(sample_message, sample_tree)
print(encoded_message)

程式列印輸出的編碼結果為:

0110010101110

建樹過程

下面的過程建樹編碼演算法,它以一個符號-頻度對偶表位引數(其中沒有任何符號出現在多於一個對偶中),並根據Huffman演算法生產出Huffman編碼樹:

def generate_huffman_tree(paris):
    return successive_merge(make_leaf_set(paris))

def make_leaf_set(pairs):
    if not pairs:
        return []
    else:
        pair = pairs[0]
        return adjoin_set(make_leaf(pair[0], pair[1]), \
                          make_leaf_set(pairs[1: ]))
        
def successive_merge(leaf_set):
    if len(leaf_set) == 1:
        return leaf_set
    left, right = leaf_set[0], leaf_set[1]
    return successive_merge(adjoin_set(make_code_tree(left, right), leaf_set[2: ]))
    
def adjoin_set(x, set):
    if not set:
        return [x]
    if weight(x) < weight(set[0]):
        return [x] + set 
    else:
        return [set[0]] + adjoin_set(x, set[1: ])

注意這裡的make_leaf_set過程將對偶表變換為葉的有序集,successive_merge使用make_code_tree反覆歸併集合中具有最小權重的元素,直至集合中只剩下一個元素為止。這個元素就是我們所需要的Huffman樹(這一過程稍微有點技巧性,但並不複雜,因為我們利用了有序集合表示這一事實)。

然後定義如下的對偶表:

pairs = ([("A", 4), ("B", 2), ("D", 1), ("C", 1)])

我們用generate_huffman_tree完成由該對偶表到Huffman樹的構建,得到Huffman樹的構建結果:

tree = generate_huffman_tree(pairs)
print(tree)

程式列印輸出的Huffman樹構建結果為:

# [['leaf', 'A', 4], [['leaf', 'B', 2], [['leaf', 'C', 1], ['leaf', 'D', 1], 'CD', 2], 'BCD', 4], 'ABCD', 8]

這實質上和我們前面所畫出來的Huffman樹在功能上是等效的,但具體的某些左右子樹可能有差異,因為Huffman樹本身就不唯一。

參考

  • [1] Abelson H, Sussman G J. Structure and interpretation of computer programs[M]. The MIT Press, 1996.
  • [2] Cormen T H, Leiserson C E, Rivest R L, et al. Introduction to algorithms[M]. MIT press, 2022.

相關文章