Sequence

Sequence

「數列」。一串數字。

(4 -1 6 0 9)

Sequence運算

 	result (noun)
---	---------------------------------
+	direct sum            直和(加法)
×	direct product        直積(乘法)
·	dot product           點積
∗	convolution           摺積        ⇒ Toeplitz matrix  常對角矩陣
⊛	circular convolution  循環摺積     ≡ circulant matrix 循環矩陣
註:芒星asterisk、角星star是不同東西。
  芒星*(U+002A)、角星★(U+2605)、芒星運算子∗(U+2217)、角星運算子⋆(U+22C6)
  摺積運算符號是六芒星,最理想的字元是芒星運算子∗(U+2217)。
  根據字型選擇,芒星運算子可能顯示五芒星或六芒星。

加減乘除模

對應項加減乘除模。時間複雜度O(N)。

(1 2 3) + (4 5 6) = (1+4 2+5 3+6) = (5 7 9)

點積

點積是對應項相乘,然後加總。時間複雜度O(N)。

(1 2 3) · (4 5 6) = 1×4 + 2×5 + 3×6 = 32

摺積

摺積是很多次點積。

第二串數列頭尾顛倒(為了配合多項式乘法),窮舉各種對齊方式,各做一次點積。沒有數字的地方視作零。時間複雜度O(AB)。

運用傅立葉轉換或數論轉換,時間複雜度O(NlogN)。

(1 2 3 4 5) ∗ (4 5 6) = (4 13 28 43 58 49 30)

(- - 1 2 3 4 5 - -) · (6 5 4 - - - - - -) = 4
(- - 1 2 3 4 5 - -) · (- 6 5 4 - - - - -) = 13
(- - 1 2 3 4 5 - -) · (- - 6 5 4 - - - -) = 28
(- - 1 2 3 4 5 - -) · (- - - 6 5 4 - - -) = 43
(- - 1 2 3 4 5 - -) · (- - - - 6 5 4 - -) = 58
(- - 1 2 3 4 5 - -) · (- - - - - 6 5 4 -) = 49
(- - 1 2 3 4 5 - -) · (- - - - - - 6 5 4) = 30

循環摺積

數列頭尾循環。時間複雜度O(N²)。

運用傅立葉轉換或數論轉換,時間複雜度O(NlogN)。

(1 2 3 4 5) ⍟ (5 6 7 8 9) = (105 110 110 105 95)

(1 2 3 4 5) · (5 9 8 7 6) = 105
(1 2 3 4 5) · (6 5 9 8 7) = 110
(1 2 3 4 5) · (7 6 5 9 8) = 110
(1 2 3 4 5) · (8 7 6 5 9) = 105
(1 2 3 4 5) · (9 8 7 6 5) = 95

循環矩陣

每行每列皆循環的方陣。

[ 5 9 8 7 6 ]
[ 6 5 9 8 7 ]
[ 7 6 5 9 8 ]
[ 8 7 6 5 9 ]
[ 9 8 7 6 5 ]

循環摺積=循環矩陣求值。

                        [ 5 9 8 7 6 ] [ 1 ]
    (1 2 3 4 5)         [ 6 5 9 8 7 ] [ 2 ]
         ⊛       <--->  [ 7 6 5 9 8 ] [ 3 ]
    (5 6 7 8 9)         [ 8 7 6 5 9 ] [ 4 ]
                        [ 9 8 7 6 5 ] [ 5 ]
circular convolution  circulant matrix evaluation

求值、求解、乘法、反矩陣,運用傅立葉轉換或數論轉換,時間複雜度為O(NlogN)。

常對角矩陣

每一條左上右下斜線是相同元素的矩陣。

[ 1 5 3 2 ]
[ 2 1 5 3 ]
[ 7 2 1 5 ]
[ 0 7 2 1 ]
[ 8 0 7 2 ]
[ 9 8 0 7 ]

摺積可以化作常對角矩陣求值。反之則不然。

                                         [ 0 ]
                   [ 8 7 6 0 0 0 0 0 0 ] [ 0 ]
                   [ 0 8 7 6 0 0 0 0 0 ] [ 1 ]
(1 2 3 4 5)        [ 0 0 8 7 6 0 0 0 0 ] [ 2 ]
     ∗       --->  [ 0 0 0 8 7 6 0 0 0 ] [ 3 ]
  (6 7 8)          [ 0 0 0 0 8 7 6 0 0 ] [ 4 ]
                   [ 0 0 0 0 0 8 7 6 0 ] [ 5 ]
                   [ 0 0 0 0 0 0 8 7 6 ] [ 0 ]
                                         [ 0 ]
convolution        Toeplitz matrix evaluation
                   [ 6 0 0 0 0 ]
                   [ 7 6 0 0 0 ] [ 1 ]
(1 2 3 4 5)        [ 8 7 6 0 0 ] [ 2 ]
     ∗       --->  [ 0 8 7 6 0 ] [ 3 ]
  (6 7 8)          [ 0 0 8 7 6 ] [ 4 ]
                   [ 0 0 0 8 7 ] [ 5 ]
                   [ 0 0 0 0 8 ]
convolution    Toeplitz matrix evaluation

加法:兩個常對角矩陣相加。O(N)。

2N-1個數字就能記錄一個常對角矩陣。

乘法:兩個常對角矩陣相乘。O(N²)。

http://stackoverflow.com/questions/15889521/

反矩陣(乘法反元素):O(N²)。高斯喬登消去法是O(N³)。

http://www.sciencedirect.com/science/article/pii/S0893965907000535

求值:O(NlogN)。普通矩陣求值是O(N²)。

填充元素成為循環矩陣,化為循環摺積,套用傅立葉轉換、數論轉換。最後從中抽取正確答案。

                     [ 1 5 3 7 2 ] [ 1 ]
[ 1 5 3 ] [ 1 ]      [ 2 1 5 3 7 ] [ 2 ]
[ 2 1 5 ] [ 2 ] ---> [ 7 2 1 5 3 ] [ 3 ]
[ 7 2 1 ] [ 3 ]      [ 3 7 2 1 5 ] [ 0 ]
                     [ 5 3 7 2 1 ] [ 0 ]
Toeplitz matrix      circulant matrix   

求解(求值反運算):O(N²)。高斯消去法是O(N³)。

https://en.wikipedia.org/wiki/Levinson_recursion

求解也有時間複雜度O(N(logN)²)的演算法。我沒有涉獵。

Recurrence

Recurrence

「遞迴數列」。以遞迴函數得到的數列。

recursive function:
{ f(0) = 1
{ f(n) = 2 f(n-1)² - 4

recursive sequence:
1 -2 4 28 1564 ......

數學當中,遞迴數列與遞迴函數一體兩面,同稱Recurrence。計算學當中,則是各飾一角。本文為了方便講解,暫將Recurrence視作遞迴數列。

Linear Recurrence

「線性遞迴數列」。以線性遞迴函數得到的數列。

{ f(0) = c₀
{ f(1) = c₁
{   :     :
{ f(n) = a₁ f(n-1) + a₂ f(n-2) + ... + aₖ f(n-k)

知名範例是Fibonacci Sequence。

recursive function:          f(0) = 0
{ f(0) = 0                   f(1) = 1
{ f(1) = 1                   f(2) = f(1) + f(0) = 1
{ f(n) = f(n-1) + f(n-2)     f(3) = f(2) + f(1) = 2
                             f(4) = f(3) + f(2) = 3
recursive sequence:          f(5) = f(4) + f(3) = 5
0 1 1 2 3 5 ...                :      :      :    :

從頭開始一項一項算。時間複雜度O(K + (N-K)K),可以簡單寫成O(NK)。

Companion Matrix

線性函數,即是矩陣。線性遞迴函數,即是「相伴矩陣」。

線性遞迴函數。
{ f(0) = 0
{ f(1) = 1
{ f(2) = 2
{ f(n+3) = 7 f(n) + 8 f(n+1) + 9 f(n+2)

構造向量F₀、F₁、...、Fn來表示函數輸入。
     [ f(0) ]   [ 0 ]        [ f(1) ]              [ f(n)   ]
F₀ = [ f(1) ] = [ 1 ]   F₁ = [ f(2) ]   ...   Fn = [ f(n+1) ]
     [ f(2) ]   [ 2 ] ,      [ f(3) ] ,     ,      [ f(n+2) ]

構造相伴矩陣A來表示線性遞迴函數。
    [ 0 1 0 ]
A = [ 0 0 1 ]
    [ 7 8 9 ]

F₀通過A就變成F₁,F₁通過A就變成F₂,以此類推。
重點在於A最底端的橫條,其他只是簡單位移。
       [ 0 1 0 ] [ f(0) ]   [   0   + 1f(1) +   0   ]   [ f(1) ]
A F₀ = [ 0 0 1 ] [ f(1) ] = [   0   +   0   + 1f(2) ] = [ f(2) ] = F₁
       [ 7 8 9 ] [ f(2) ]   [ 7f(0) + 8f(1) + 9f(2) ]   [ f(3) ]

以此類推的過程,簡化成A的次方,最後得到Fn。
A F₀   = A¹ F₀ = F₁
A F₁   = A² F₀ = F₂
A F₂   = A³ F₀ = F₃
 :         :      :
A Fn-1 = Aⁿ F₀ = Fn

我們想要的第N項,位於Fn。
                      [ f(n)   ]   <------- here is f(n)
A Fn-1 = Aⁿ F₀ = Fn = [ f(n+1) ]
                      [ f(n+2) ]

也可以提早一點,位於Fn-2。
                          [ f(n-2) ]
A Fn-3 = Aⁿ⁻² F₀ = Fn-2 = [ f(n-1) ]
                          [ f(n)   ]   <--- here is f(n)

遞迴函數,即是函數複合許多次。線性遞迴函數,即是線性函數複合許多次,即是相伴矩陣矩陣次方。

recursive function:        linear recursive function:
{ f(0) = 1                 { f(0) = 0
{ f(n) = 2 f(n-1)² - 4     { f(1) = 1
                           { f(2) = 2
function compositions:     { f(n+3) = 7 f(n) + 8 f(n+1) + 9 f(n+2)
let g(x) = 2x² - 4
f(0) = 1                   companion matrix exponentiation:
f(1) = g(1)                [ f(n)   ]   [ 0 1 0 ]ⁿ [ f(0) ]
f(2) = g(g(1))             [ f(n+1) ] = [ 0 0 1 ]  [ f(1) ]
f(3) = g(g(g(1)))          [ f(n+2) ]   [ 7 8 9 ]  [ f(2) ]

矩陣次方有快速演算法,原理是Divide and Conquer。

線性遞迴數列,直接求得第N項,僅需O(log(N-K))次矩陣乘法。時間複雜度O(K² + K³log(N-K)),可以簡單寫成O(K³logN)。此處假設矩陣相乘是O(K³)。

Polynomial Multiplication and Polynomial Modulo

2015年由競賽選手台灣大學李瑞珉《High Speed Linear Recursion》提出。我不清楚是否已有學術論文。

[shift]   Fn = w1 Fn-1 + ... + wk Fn-k
             = (w1 w1 Fn-2 + ... + w1 wk Fn-k-1) + ... + wk Fn-k
             = v2 Fn-2 + ... + vk-1 Fn-k-1

 RHS第一項再展開,然後同類項的係數相加。
 特色是,我們可以持續位移到我們喜歡的那k項(連續的)!
 位移的時間複雜度O(k),移到我們喜歡的那k項是 O(nk)。
 (但是以下不會用到。)

[expand]   Fn = w1 Fn-1 + ... + wk Fn-k
              = (w1 w1 Fn-2   + ... + w1 wk Fk-2  ) + ...... +
                (wk w1 Fn-k-1 + ... + wk wk Fn-k-k)
              = v2 Fn-2 + ............ + v2k Fn-2k

 RHS每個項各自再展開,然後同類項的係數相加。
 直接窮舉,時間複雜度 O(k²)。

 其實這就是多項式直式乘法(摺積),摺積有O(klogk)算法。
 (如果用快速數論轉換去算摺積,如何保證w1...wk不會爆表呢?)
 (所以只好請出中國餘數定理)

[expand]   Fn = w1 Fn-x-1   + ... + wk  Fn-x-k
              = v2 Fn-x-y-2 + ... + v2k Fn-x-y-2k   第二次展開,也可以用別的
              = uk Fn-y-k-1 + ... + u2k Fn-x-y-2k   縮成剛好k項

 因為有了[shift]的概念,所以我們取的那k項
 沒必要是 Fn-1  ... Fn-k
 也可以是 Fn-x-1 ... Fn-x-k
 最後再用 k-1次 [shrink] 縮成剛好k項

[shrink]   Fn = v1 Fn-1 + v2 Fn-2 + ... + vk Fn-k + ... + vK Fn-K
              =           u2 Fn-2 + ... + uk Fn-k + ... + uK Fn-K

 項數大於k的時候,RHS第一項再展開,可以縮減一項。
 縮一項的複雜度 O(k)。
 縮到剩下k項的時間複雜度是 O(k(K-k))。
 
 其實這就是多項式長除法(減法變加法。除數乘上一個負號)
 最後要取餘數。
 既然已經採取快速數論轉換,那就是要取模了。
 取模的多項式餘數,有O(klogk)算法。

令那k項是 Fn/2-1 ... Fn/2-k,總是取一半的位置。
如此一來,仿照整數次方(Exp)或餘數次方(ModExp)的演算法。
利用 [expand] [shrink] ,分割 logn 次。
     O(klogk) O(klogk)
求第n項只需要 O(klogk logn) 時間。

總結

線性遞迴函數K項,求數列第N項:

一、動態規劃,第0項到第N項。O(NK)。

二、相伴矩陣N次方。O(K³logN)。假設矩陣相乘O(K³)。

三、多項式乘法與餘數。O(K²logN)甚至O(KlogKlogN)。

四、生成函數。手工推導公式解。詳情請查閱離散數學書籍。

當K是常數,一變成O(N),二三變成O(logN)。

UVa 10229 10518 10655 10743 10754 10870 12653

Filter

Filter

「濾波器」。數列到數列的函數。輸入數列、輸出數列。

每個輸出變數分開來看,濾波器其實是由許多函數組成。

簡單的例子是每項加1的濾波器、每項延遲1時刻的濾波器。

當全部函數都相同,僅索引值(時刻)不同,可以簡化成一個函數。

Linear Time-Invariant Filter(LTI Filter)

「線性非時變濾波器」。濾波器是相同的線性函數。

數學家和計算學家擅長線性函數,容易計算。

LTI Filter可以寫成線性組合形式

數列x,通過LTI Filter f,得到數列y:
y(n) = a(0)x(n) + a(1)x(n-1) + a(2)x(n-2) + ... + a(k)x(n-k)

教科書最喜歡使用的形式。

LTI Filter可以寫成矩陣形式

[ a(0) 0    0    0 ... ]
[ a(1) a(0) 0    0 ... ]
[ a(2) a(1) a(0) 0 ... ]
[        :             ] [ x(0) ]   [ y(0) ]
[ a(k)~a(0) 0 0 0  ... ] [ x(1) ]   [ y(1) ]
[ 0 a(k)~a(0) 0 0  ... ] [   :  ] = [   :  ]
[ 0 0 a(k)~a(0) 0  ... ] [   :  ]   [   :  ]
[         :            ] [ x(n) ]   [ y(n) ]
[ ...  0 a(k)~a(0) 0 0 ]
[ ...  0 0 a(k)~a(0) 0 ]
[ ...  0 0 0 a(k)~a(0) ]
           A                 x          y

數列看做向量,就可以寫成矩陣形式。

可援引線性代數。例如一連串的LTI Filter,複合成一個LTI Filter!又例如計算eigenvector,找出不受LTI Filter影響的數列!

LTI Filter可以寫成摺積形式

{ x(0) ... x(n)     } dot { a(0)                           } = y(0)
{ x(0) ... x(n)     } dot { a(1) a(0)                      } = y(1)
        :                                  :                      :
{ x(0) ... x(n)     } dot {      a(k) .......... a(0)      } = y(n-1)
{ x(0) ... x(n)     } dot {           a(k) .......... a(0) } = y(n)
—————————————————————————————————————————————————————————————————————
{ x(0) ... x(n)     } dot {                a(k) ..... a(1) } = 0
        :                                              :      :
{ x(0) ... x(n)     } dot {                     a(k) a(k-1)} = 0
{ x(0) ... x(n)     } dot {                           a(k) } = 0

    (x append 0s)   convolution           a           = (y append 0s)

數列末端補零K個,就可以寫成摺積形式。

可援引傅立葉轉換。時域循環摺積=頻域乘法。

LTI Filter運算

evaluation f(x) = ?              find y     (moving average)
resolution f(?) = y              find x     (deconvolution)
regression ?(x) = y              find f     (Wiener filter)
inversion  f(x) = y, x = ?(y)    find f⁻¹

四種運算:求值(通過濾波器)、求解(恢復原數列)、迴歸(求濾波器)、反函數(求反濾波器)。

四種運算都有時域和頻域兩種解法。儘管頻域解法的時間複雜度較低,但是沒人用!一、濾波器長度通常很短(函數項數很少),傅立葉轉換反而浪費時間。二、濾波器無法完美轉換到頻域。

求值(通過濾波器)

evaluation f(x) = ?              find y     (moving average)

寫成矩陣形式,化作常對角矩陣求值。時間複雜度O(N²)。

[ a(0)   0    0  ...   0    0    0  ] [ x(0) ]   [ y(0) ]
[ a(1) a(0)   0  ...   0    0    0  ] [ x(1) ]   [ y(1) ]
[ a(2) a(1) a(0) ...   0    0    0  ] [   :  ]   [   :  ]
[   :    :    :        :    :    :  ] [   :  ] = [   :  ]
[   0    0    0  ... a(0)   0    0  ] [   :  ]   [   :  ]
[   0    0    0  ... a(1) a(0)   0  ] [   :  ]   [   :  ]
[   0    0    0  ... a(2) a(1) a(0) ] [ x(n) ]   [ y(n) ]
                  A                       x          y

滑動視窗,取加權平均值。時間複雜度O(NK)。

求解(恢復原數列)

resolution f(?) = y              find x     (deconvolution)

寫成矩陣形式,化作常對角矩陣求解。時間複雜度O(N²)。

滑動視窗,解加權平均值方程式。時間複雜度O(NK)。

迴歸(求濾波器)

regression ?(x) = y              find f

寫成另一種矩陣形式,化作常對角矩陣求解。

等式太多,通常無解。折衷令誤差最小,誤差是Least Squares,Xa和y對應項的差的平方的總和。套用「Normal Equation」,化作一般矩陣求解,時間複雜度O(N³)。

[ x(0)     0    ...   0        0      ]            [ y(0) ]
[ x(1)   x(0)   ...   0        0      ] [ a(0) ]   [ y(1) ]
[   :      :          :        :      ] [ a(1) ]   [   :  ]
[   :      :          :        :      ] [   :  ] = [   :  ]
[   :      :          :        :      ] [   :  ]   [   :  ]
[ x(n-1) x(n-2) ... x(n-k)   x(n-k-1) ] [ a(k) ]   [   :  ]
[ x(n)   x(n-1) ... x(n-k+1) x(n-k)   ]            [ y(n) ]
                 X                          a          y

對應項平方誤差總和最小,即是滿足 Xᵀ X a = Xᵀ y。
a = (Xᵀ X)⁻¹ Xᵀ y

迴歸(求濾波器):旁門左道解法

regression ?(x) = y              find f     (Wiener filter)

x y末端補零K個。缺點:不合邏輯,得到歪解。優點:數學式子非常美觀、非常好算,稱做Wiener-Hopf Equation。

一、建方程式:摺積。O(NK)。

二、解方程式:常對角矩陣求解。O(K²)。

總時間複雜度O(NK + K²),或者簡單寫成O(N²)。

[ x(0)                      ]            [ y(0) ]
[ :    x(0)                 ]            [   :  ]
[ :    :                    ] [ a(0) ]   [   :  ]
[ :    :          x(0)      ] [   :  ]   [   :  ]
[ :    :    ..... :    x(0) ] [   :  ] = [   :  ]
[ x(n) :          :    :    ] [   :  ]   [ y(n) ]
[      x(n)       :    :    ] [ a(k) ]   [   0  ]
[                 :    :    ]            [   :  ]
[                 x(n) :    ]            [   :  ]
[                      x(n) ]            [   0  ]
              X                   a          y

[ xx(0)   xx(1)   ... xx(k)   ] [ a(0) ]   [ xy(0) ]
[ xx(1)   xx(0)   ... xx(k-1) ] [ a(1) ]   [ xy(1) ]
[    :       :           :    ] [   :  ] = [    :  ]
[ xx(k-1) xx(k-2) ... xx(1)   ] [   :  ]   [    :  ]
[ xx(k)   xx(k-1) ... xx(0)   ] [ a(k) ]   [ xy(k) ]
             Xᵀ X                   a         Xᵀ y

xx(t) = ∑ x(i-t)⋅x(i)   autocorrelation of x(i-t) and x(i).
xy(t) = ∑ x(i-t)⋅y(i)   cross-correlation of x(i-t) and y(i).

對應項平方誤差總和最小,即是滿足 Xᵀ X a = Xᵀ y。
Xᵀ X 和 Xᵀ y 乘起來,稱做Wiener-Hopf Equation。
xx(t):數列x延遲t時刻、數列x,兩者的自相關數(相差t時刻,兩串相同數列的點積)。
xy(t):數列x延遲t時刻、數列y,兩者的互相關數(相差t時刻,兩串不同數列的點積)。

另一種演算法是使用循環摺積、循環矩陣。

一、建方程式:摺積再度補零=循環摺積=頻域的乘法。O(NlogN)。

二、解方程式:循環矩陣求解=循環反摺積=頻域的除法。O(KlogK)。

總時間複雜度O(NlogN + KlogK),或者簡單寫成O(NlogN)。

甚至打從一開始就用傅立葉轉換,將數列和濾波器轉換到頻域,合併兩步驟,得到公式解â = x̂ŷ / x̂²,稱做Wiener Filter。

反函數(求反濾波器)

inversion  f(x) = y, x = ?(y)    find f⁻¹

濾波器通常不是一對一函數,通常不存在反函數。什麼時候存在反函數呢?我查不到相關資料!

我也沒找到反函數演算法!目前有兩種做法:

一、隨便找一組輸出輸入,對調一下,然後迴歸。

二、濾波器,轉換到頻域,取倒數,轉換回時域。

Recurrence暨LTI Filter運算

autoevaluation f(?) = ? >> 1     find x
autoregression ?(x) = x >> 1     find f
resolution     { f(?) = ? >> 1   find x     (Kalman filter)
               { h(?) = y

遞迴數列可以視作數列暨濾波器:輸入與輸出是同一數列,但是輸出延遲1時刻。

遞迴數列可以做為濾波器的輸入。

引入遞迴數列,多了三種運算:自求值(求遞迴數列)、自迴歸(求遞迴函數)、求解(恢復原數列)。

三種運算應該是首次歸類在一起,名稱是我暫訂的。

自求值(求遞迴數列)

autoevaluation f(?) = ? >> 1     find x

前面章節介紹過了。線性遞迴函數K項,求數列第N項:

一、動態規劃,第0項到第N項。O(NK)。

二、相伴矩陣,矩陣的N次方。O(KωlogN)。

三、多項式乘法、多項式餘數。O(K²logN)或O(KlogKlogN)。

四、生成函數,手工推導公式解。

當K是常數,一變成O(N),二三變成O(logN)。

自迴歸(求遞迴函數)

autoregression ?(x) = x >> 1     find f

每一個數值,皆是先前緊鄰的K個數值的加權總和,權重是定值。反覆套用線性遞迴函數、代入數列最後N個數值,就能預測下一個即將出現的數值。此行為稱作「線性預測Linear Prediction」。

一串長長的數列,壓縮成一個線性遞迴函數,只需儲存K個函數係數、K個初始數值。反覆套用函數、代入數列最後K個數值,解壓縮得到一串長長的數列。此行為稱作「線性預測編碼Linear Predictive Coding」。

旁門左道的解法,得到Yule-Walker Equation。

視作普通摺積、常對角矩陣,時間複雜度簡單寫成O(N²)。

視作循環摺積、循環矩陣,時間複雜度簡單寫成O(NlogN)。

x(n) = a(1)x(n-1) + a(2)x(n-2) + ... + a(k)x(n-k)

[ xx(0)   xx(1)   ... xx(k-1) ] [ a(1)   ]   [ xx(1)   ]
[ xx(1)   xx(0)   ... xx(k-2) ] [ a(2)   ]   [ xx(2)   ]
[    :       :           :    ] [   :    ] = [    :    ]
[ xx(k-2) xx(k-1) ... xx(1)   ] [ a(k-1) ]   [ xx(k-1) ]
[ xx(k-1) xx(k-2) ... xx(0)   ] [ a(k)   ]   [ xx(k)   ]

求解(恢復原數列)

resolution     { f(?) = ? >> 1   find x     (Kalman filter)
               { h(?) = y

使用測量儀器h,時時收集數據y,時時推敲現況x。推定現況是自迴歸濾波器f,合乎自然規律。自創測量儀器是濾波器h,合乎電路原理。根據結果,推敲原因,此行為稱作「線性二次估計Linear Quadratic Estimation」。

輸出誤差,通過反濾波器,得到輸入誤差,用以修正輸入。

為什麼要大費周章,將觀察值的誤差通過反濾波器,而不是直截了當,將觀察值通過反濾波器?我也不清楚。

h(x) = y
狀態 x₀ x₁ x₂ ...   未知(一串數列,線性遞迴關係)
函數 h₀ h₁ h₂ ...   已知(數列每一項皆有一個函數)
觀察 y₀ y₁ y₂ ...   已知(數列每一項皆有一個輸出)

函數可以一樣               h₀ = h₁ = h₂ = ...
                          
狀態是自迴歸的LTI filter   xₙ = f(xₙ₋₁, xₙ₋₂, ..., xₙ₋ₖ)
                           xₙ = a₁xₙ₋₁ + a₂xₙ₋₂ + ... + aₚxₙ₋ₖ
                
函數可以推廣成LTI filter   yₙ = hₙ(xₙ, xₙ₋₁, xₙ₋₂, ..., xₙ₋ₚ)
                           yₙ = bₙ₀xₙ + bₙ₁xₙ₋₁ + ... + bₙₚxₙ₋ₚ

從 x₁ 開始一個一個算到 xᵢ,現在要算 xᵢ 是多少:

1. x̂ᵢ = f(xᵢ₋₁)              x 數列估計下一項。 (x₀ = 0)
                             此處以一階為例。

2. ŷᵢ = hᵢ(x̂ᵢ)               y 數列估計下一項。
                            盜版 ŷᵢ 通常不等於正版 yᵢ。
                             這表示 x̂ᵢ 不夠準,需要修正。

3. kᵢ⁻¹                      求得誤差的反濾波器 kᵢ⁻¹。稱做Kalman gain。
                             kᵢ(x̂ᵢ - xᵢ) = (ŷᵢ - yᵢ)

 (1) pᵢ = (x̂ᵢ - xᵢ)²         為了計算反濾波器,需要 xᵢ 的誤差的平方。
        = cov(x̂ᵢ - xᵢ)       即是誤差的covariance。
     sᵢ = (ŷᵢ - yᵢ)²         不過由於我們還沒求出 xᵢ,所以不能這樣算。
        = cov(ŷᵢ - yᵢ)

 (2) p̂ᵢ = f²(pᵢ₋₁)           p 數列估計下一項。 (p₀ = 0)
        = f pᵢ₋₁ fᵀ          原理是 Var[k⋅X] = k² ⋅ Var[X]。

 (3) ŝᵢ = hᵢ²(p̂ᵢ)            s 數列估計下一項。
        = hᵢ p̂ᵢ hᵢᵀ          如果 f h 各有誤差項:常態分布 q r,平均值0。
                             那麼 p̂ᵢ ŝᵢ 各加上 var(q) var(r)。

 (4) kᵢ⁻¹ = p̂ᵢ hᵢᵀ ŝᵢ⁻¹      反濾波器是變異數的比值。這是什麼鬼啦?

4. xᵢ = x̂ᵢ - kᵢ⁻¹(ŷᵢ - yᵢ)   修正 x。
                             yᵢ 的誤差,通過反濾波器,還原成 xᵢ 的誤差。
                             x̂ᵢ 補上誤差,答案就對了。

5. pᵢ = p̂ᵢ - kᵢ⁻²(ŝᵢ)        修正 p。
      = p̂ᵢ - kᵢ⁻¹ ŝᵢ kᵢ⁻¹ᵀ

Circular Convolution

Circular Convolution

數列、多項式、多項式函數!

數列循環摺積的高速演算法:O(NlogN)!

Duality: Circular Convolution ⟺ Multiplication

(2 -5 1 0 4)  <———>  2x⁰ - 5x¹ + 1x² + 0x³ + 4x⁴

數列與多項式互相轉換!

               2x⁰ - 5x¹ + 1x² + 0x³ + 4x⁴
(2 -5 1 0 4)  <———————————————————————————> (2 2 60 320 5192)
                   x = {0, 1, 2, 3, 6}

係數與函數值互相轉換!

正向轉換:多項式求值N回合O(N²)。

反向轉換:多項式內插O(N³)。

[ 0⁰ 0¹ 0² 0³ 0⁴ ] [  2 ]   [    2 ]
[ 1⁰ 1¹ 1² 1³ 1⁴ ] [ -5 ]   [    2 ]
[ 2⁰ 2¹ 2² 2³ 2⁴ ] [  1 ] = [   60 ]
[ 3⁰ 3¹ 3² 3³ 3⁴ ] [  0 ]   [  320 ]
[ 6⁰ 6¹ 6² 6³ 6⁴ ] [  4 ]   [ 5192 ]

改寫成矩陣形式。稱做Vandermonde Matrix。

正向轉換:矩陣求值O(N²)。

反向轉換:先求反矩陣O(N³),再求值O(N²)。

    [ x₀⁰ x₀¹ x₀² x₀³ ] [  2 ]   [ y₀ ]
 1  [ x₁⁰ x₁¹ x₁² x₁³ ] [ -5 ]   [ y₁ ]
——— [ x₂⁰ x₂¹ x₂² x₂³ ] [  1 ] = [ y₂ ]
 k  [ x₃⁰ x₃¹ x₃² x₃³ ] [  0 ]   [ y₃ ]

Fourier transform:
x = {e-i⋅(2π/N)⋅0, e-i⋅(2π/N)⋅1, ..., e-i⋅(2π/N)⋅(N-1)}
e = Euler number
k = 1 / √n

number theoretic transform:
x = {r⁰, r¹, ……, rN-1}
r = anyone primitive root 
k = 1

x設定成特殊數值,矩陣大小恰是2的次方,有高速演算法!

已知兩種設定方式:複數系統傅立葉轉換、餘數系統數論轉換。

正向轉換:O(NlogN)。

反向轉換:O(NlogN)。

               2x⁰ - 5x¹ + 1x² + 0x³ + 4x⁴
(2 -5 1 0 4)  <———————————————————————————> (2 2 60 320 5192)
     ∗         1x⁰ - 1x¹ + 0x² + 1x³ + 2x⁴          ×
(1 -1 0 1 -2) <———————————————————————————> (1 -1 -25 -137 -2381)
               2x⁰ - 7x¹ + 6x² + 1x³ - 5x⁴
     ‖         + 7x⁵ - 2x⁶ + 4x⁷ - 8x⁸              ‖
(2 -7 6 1 -5  <———————————————————————————> (2 -2 -1500 -43840 -12362152)
 7 -2 4 -8)        x = {0, 1, 2, 3, 6}

係數與函數值互相轉換之特色:

一、係數加減=多項式加減=函數值加減。(這不是重點)

二、係數摺積與反摺積=多項式乘除=函數值乘除。

               (2x⁰-5x¹+1x²+0x³+4x⁴) ÷ √5
(2 -5 1 0 4)  <——————————————————————————> (0.8, 3.5-0.7i, 0.5+3.0i,
                                            0.5-3.0i, 3.5+0.7i)
      ×        (1x⁰-1x¹+0x²+1x³+2x⁴) ÷ √5             ⊛
(1 -1 0 1 -2) <——————————————————————————> (-0.4, -0.2-0.2i, -1.7+0.4i,
                                            -1.7-0.4i, -0.2+0.2i)
      ‖        (2x⁰+5x¹+0x²+0x³-8x⁴) ÷ √5             ‖
(2 5 0 0 -8)  <——————————————————————————> (-0.4, -5.1+2.1i, -3.6-1.6i,
               x = {e-i⋅(2π/5)⋅0, e-i⋅(2π/5)⋅1   -3.6+1.6i, -5.1-2.1i)
                    e-i⋅(2π/5)⋅2, e-i⋅(2π/5)⋅3
                    e-i⋅(2π/5)⋅4}

係數摺積=函數值乘法。

如果是對偶轉換(正向轉換=反向轉換),那麼性質就對稱了:係數乘法=函數值摺積。

改成對偶轉換,目前沒有辦法。改成幾乎是對偶變換,則有一種辦法:令多項式次方循環。

已知兩種設定方式:複數系統傅立葉轉換、餘數系統數論轉換。

多項式次方循環,導致性質稍有變化:

一、係數頭尾循環=多項式次方循環=函數值頭尾循環。

二、係數循環摺積=函數值乘法。

三、係數乘法=函數值循環摺積。

            fft
(2 -5 1 0) —————> (? ? ? ?)
     ⊛      fft       ×
(1 -1 0 1) —————> (? ? ? ?)
     ‖      ifft      ↓
(? ? ? ?)  <————— (? ? ? ?)

當數列長度恰是2的次方,循環摺積有高速演算法!

循環摺積:正向傅立葉轉換、相乘、逆向傅立葉轉換。

總時間複雜度O(NlogN)。

以上只是綱要。以下則是細節。

點積、摺積、循環摺積

點積是對應項相乘,然後加總。

(a₀ a₁ a₂) dot (b₀ b₁ b₂) = a₀b₀ + a₁b₁ + a₂b₂

摺積是很多次點積。

第二串數列頭尾顛倒,窮舉各種對齊方式,各做一次點積。

第二串數列頭尾顛倒,是為了配合多項式乘法。

數列摺積
  (a₀ a₁ a₂) convolution (b₀ b₁ b₂) = (c₀ c₁ c₂)

其中
c₀ = a₀b₀
c₁ = a₀b₁ + a₁b₀
c₂ = a₀b₂ + a₁b₁ + a₂b₀
c₃ = a₁b₂ + a₂b₁
c₄ = a₂b₂

直式
c₀:                     c₁:                     c₂:
(-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )
(b₂ b₁ b₀ -  -  -  - )  (-  b₂ b₁ b₀ -  -  - )  (-  -  b₂ b₁ b₀ -  - )

c₃:                     c₄:
(-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )
(-  -  -  b₂ b₁ b₀ - )  (-  -  -  -  b₂ b₁ b₀)

橫式
c₀ = (-  -  a₀ a₁ a₂ -  - ) dot (b₂ b₁ b₀ -  -  -  - )
c₁ = (-  -  a₀ a₁ a₂ -  - ) dot (-  b₂ b₁ b₀ -  -  - )
c₂ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  b₂ b₁ b₀ -  - )
c₃ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  -  b₂ b₁ b₀ - )
c₄ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  -  -  b₂ b₁ b₀)

循環摺積,超出數列的部分,改成頭尾循環。

                  circular
(a₀ a₁ a₂ a₃ a₄) convolution (b₀ b₁ b₂ b₃ b₄) = (c₀ c₁ c₂ c₃ c₄)

直式
c₀:               c₁:               c₂:
(a₀ a₁ a₂ a₃ a₄)  (a₀ a₁ a₂ a₃ a₄)  (a₀ a₁ a₂ a₃ a₄)
(b₀ b₄ b₃ b₂ b₁)  (b₁ b₀ b₄ b₃ b₂)  (b₂ b₁ b₀ b₄ b₃)

c₃:               c₄:
(a₀ a₁ a₂ a₃ a₄)  (a₀ a₁ a₂ a₃ a₄)
(b₃ b₂ b₁ b₀ b₄)  (b₄ b₃ b₂ b₁ b₀)

橫式
c₀ = (a₀ a₁ a₂ a₃ a₄) dot (b₀ b₄ b₃ b₂ b₁)
c₁ = (a₀ a₁ a₂ a₃ a₄) dot (b₁ b₀ b₄ b₃ b₂)
c₂ = (a₀ a₁ a₂ a₃ a₄) dot (b₂ b₁ b₀ b₄ b₃)
c₃ = (a₀ a₁ a₂ a₃ a₄) dot (b₃ b₂ b₁ b₀ b₄)
c₄ = (a₀ a₁ a₂ a₃ a₄) dot (b₄ b₃ b₂ b₁ b₀)

數列=多項式

(a₀ a₁ a₂)   <--->   a₀x⁰ + a₁x¹ + a₂x²

數列加法與減法=多項式加法與減法

數列加法
(a₀ a₁ a₂) + (b₀ b₁ b₂) = (a₀+b₀ a₁+b₁ a₂+b₂)
多項式加法
        a₀ x⁰ +       a₁ x¹ +       a₂ x²
+)      b₀ x⁰ +       b₁ x¹ +       b₂ x²
—————————————————————————————————————————
   (a₀+b₀) x⁰ +  (a₁+b₁) x¹ +  (a₂+b₂) x²

數列摺積與反摺積=多項式乘法與除法

多項式乘法
                         a₀ x⁰ +   a₁ x¹ +   a₂ x²
×)                       b₀ x⁰ +   b₁ x¹ +   b₂ x²
——————————————————————————————————————————————————
                       a₀b₂ x² + a₁b₂ x³ + a₂b₂ x⁴
             a₀b₁ x¹ + a₁b₁ x² + a₂b₁ x³
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
——————————————————————————————————————————————————
     c₀ x⁰ +   c₁ x¹ +   c₂ x² +   c₃ x³ +   c₄ x⁴
數列摺積
(a₀ a₁ a₂) convolution (b₀ b₁ b₂) = (c₀ c₁ c₂ c₃ c₄)
其中
c₀ = a₀b₀
c₁ = a₀b₁ + a₁b₀
c₂ = a₀b₂ + a₁b₁ + a₂b₀
c₃ = a₁b₂ + a₂b₁
c₄ = a₂b₂

直式
c₀:                     c₁:                     c₂:
(-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )
(b₂ b₁ b₀ -  -  -  - )  (-  b₂ b₁ b₀ -  -  - )  (-  -  b₂ b₁ b₀ -  - )

c₃:                     c₄:
(-  -  a₀ a₁ a₂ -  - )  (-  -  a₀ a₁ a₂ -  - )
(-  -  -  b₂ b₁ b₀ - )  (-  -  -  -  b₂ b₁ b₀)

橫式
c₀ = (-  -  a₀ a₁ a₂ -  - ) dot (b₂ b₁ b₀ -  -  -  - )
c₁ = (-  -  a₀ a₁ a₂ -  - ) dot (-  b₂ b₁ b₀ -  -  - )
c₂ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  b₂ b₁ b₀ -  - )
c₃ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  -  b₂ b₁ b₀ - )
c₄ = (-  -  a₀ a₁ a₂ -  - ) dot (-  -  -  -  b₂ b₁ b₀)

數列循環摺積=多項式循環乘法

多項式循環乘法
                         a₀ x⁰ +   a₁ x¹ +   a₂ x²
×)                       b₀ x⁰ +   b₁ x¹ +   b₂ x²
—————————————————————————————————————————————————————
                       a₀b₂ x² + a₁b₂ x³ + a₂b₂ x⁴
             a₀b₁ x¹ + a₁b₁ x² + a₂b₁ x³
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
—————————————————————————————————————————————————————
   a₁b₂ x⁰ + a₂b₂ x¹ + a₀b₂ x²
   a₂b₁ x⁰ + a₀b₁ x¹ + a₁b₁ x²
   a₀b₀ x⁰ + a₁b₀ x¹ + a₂b₀ x²
—————————————————————————————————————————————————————
     c₀ x⁰ +   c₁ x¹ +   c₂ x²
數列循環摺積
(a₀ a₁ a₂) circular convolution (b₀ b₁ b₂) = (c₀ c₁ c₂)
其中
c₀ = a₁b₂ + a₂b₁ + a₀b₀
c₁ = a₂b₂ + a₀b₁ + a₁b₀
c₂ = a₀b₂ + a₁b₁ + a₂b₀

直式
c₀:         c₁:         c₂:
(a₀ a₁ a₂)  (a₀ a₁ a₂)  (a₀ a₁ a₂)
(b₀ b₂ b₁)  (b₁ b₀ b₂)  (b₂ b₁ b₀)

橫式
c₀ = (a₀ a₁ a₂) dot (b₀ b₂ b₁)
c₁ = (a₀ a₁ a₂) dot (b₁ b₀ b₂)
c₂ = (a₀ a₁ a₂) dot (b₂ b₁ b₀)

Convolution Theorem【尚無正式名稱】

數列轉換成多項式。

(a₀ a₁ a₂)  ---->  a₀x⁰ + a₁x¹ + a₂x²

數列轉換成多項式,可以看成點積,可以看成線性函數。

線性函數矩陣
A = [ x⁰ x¹ x² ]

數列
    [ a₀ ]
a = [ a₁ ]
    [ a₂ ]

數列轉換成多項式
                  [ a₀ ]
Aa = [ x⁰ x¹ x² ] [ a₁ ] = a₀x⁰ + a₁x¹ + a₂x²
                  [ a₂ ]

這種線性函數有個特性:變換前的摺積=變換後的乘法。

A = [ x⁰ x¹ x² ]

    [ a₀ ]       [ b₀ ]
a = [ a₁ ]   b = [ b₁ ]
    [ a₂ ]       [ b₂ ]

p = Aa = (a₀x⁰ + a₁x¹ + a₂x²)
q = Ab = (b₀x⁰ + b₁x¹ + b₂x²)

pq = Aa Ab = (a₀x⁰ + a₁x¹ + a₂x²) (b₀x⁰ + b₁x¹ + b₂x²)
           = (a convolution b) 然後添上 [ x⁰ x¹ x² x³ x⁴ ]
           = A(a convolution b)    a與b末端補零,A末端補項。

次方值乘上任意倍率,此特性一樣成立。

A = [ x⁰ x¹ x² ]

A = [ x⁰ x⁵ x¹⁰ ]

A = [ x⁰ x⁻⁷ x⁻²¹ ]

A = [ x⁰ x⁰ x⁰ ]

所有元素一齊乘上任意倍率,此特性一樣成立。

A = [ 7x⁰ 7x¹ 7x² ]

A = [ -x⁰ -x⁵ -x¹⁰ ]

A = [ 0 0 0 ]

從一個橫列推廣到一個矩陣,此特性一樣成立。

    [ 7x⁰ 7x¹  7x²   ]
A = [ -x⁰ -x⁵  -x¹⁰  ]
    [  x⁰  x⁻⁷  x⁻²¹ ]

Aa multiply Ab = A(a convolution b)

結論就是:

Circular Convolution Theorem【尚無正式名稱】

接下來繼續補強矩陣,除了滿足「變換前的摺積=變換後的乘法」,也要滿足「變換前的乘法=變換後的摺積」!

從數學來看,補強性質,達成了對稱之美;從計算學來看,補強限制,極可能產生特殊演算法。

那麼就得讓A擁有反矩陣A⁻¹,而且A和A⁻¹都具備上個段落提到的特性。一種嘗試是正規正交矩陣:A⁻¹ = Aᵀ,前述特性變成了行與列同時成立,容易觀察。

實數系統下,x次方漸增,數值越來越大,導致基底不可能垂直(內積不可能為零)。很不幸的,這種矩陣不存在。

             [ x⁰ x⁰ x⁰ x⁰ .. ]
     -1   T  [ x⁰ x¹ x² x³ .. ]
A = A  = A = [ x⁰ x² x⁴ x⁶ .. ]
             [ x⁰ x³ x⁶ x⁹ .. ]
             [ :  :  :  :     ]

Aa multiply Ab = A(a convolution b)
Aa convolution Ab = A(a multiply b)

折衷的方式是令x的次方產生循環,讓數值能夠變小,使得基底互相垂直、內積是零。(用數學術語來說:從open set改成closed set。)

複數系統下,把x設定成ei⋅2π/N(或其倒數e-i⋅2π/N),令x的次方產生循環。此即「傅立葉轉換」。

Fourier Transform: x = e-i⋅2π/N; k = 1/sqrt(N); N is size of matrix

      [ x⁰ x⁰ x⁰ x⁰ .. ]          [ x⁻⁰ x⁻⁰ x⁻⁰ x⁻⁰ .. ]
      [ x⁰ x¹ x² x³ .. ]    -1    [ x⁻⁰ x⁻¹ x⁻² x⁻³ .. ]
A = k [ x⁰ x² x⁴ x⁶ .. ]   A  = k [ x⁻⁰ x⁻² x⁻⁴ x⁻⁶ .. ]
      [ x⁰ x³ x⁶ x⁹ .. ]          [ x⁻⁰ x⁻³ x⁻⁶ x⁻⁹ .. ]
      [ :  :  :  :     ]          [ :   :   :   :      ]

Aa multiply Ab = A(a circular convolution b)
Aa circular convolution Ab = A(a multiply b)

餘數系統下,則是把x設定成任意一個原根(或其倒數),令x的次方產生循環。此即「數論轉換」。

Number Theoretic Transform: x = primitive root (mod n)

    [ x⁰ x⁰ x⁰ x⁰ .. ]        [ x⁻⁰ x⁻⁰ x⁻⁰ x⁻⁰ .. ]
    [ x⁰ x¹ x² x³ .. ]    -1  [ x⁻⁰ x⁻¹ x⁻² x⁻³ .. ]
A = [ x⁰ x² x⁴ x⁶ .. ]   A  = [ x⁻⁰ x⁻² x⁻⁴ x⁻⁶ .. ]
    [ x⁰ x³ x⁶ x⁹ .. ]        [ x⁻⁰ x⁻³ x⁻⁶ x⁻⁹ .. ]
    [ :  :  :  :     ]        [ :   :   :   :      ]

Aa multiply Ab = A(a circular convolution b)
Aa circular convolution Ab = A(a multiply b)

原本的特性,也變得循環:「變換前的循環摺積=變換後的乘法」、「變換前的乘法=變換後的循環摺積」!

引入時域與頻域的觀念:「頻域的乘法=時域的循環摺積」、「頻域的循環摺積=時域的乘法」!

Fourier Transform

請參考本站文件「Fourier Transform」。

Circular Convolution快速演算法
Polynomial Circular Multiplication快速演算法

正向傅立葉轉換、數論轉換需時O(NlogN),對應項相乘需時O(N),逆向傅立葉轉換、數論轉換需時O(NlogN)。總時間複雜度為O(NlogN)。

傅立葉轉換的弱點是浮點數誤差。數論轉換的弱點是數值大小不得超過模數大小。

注意到,快速傅立葉轉換、快速數論轉換,數列長度必須是2的次方。當數列長度不是2的次方,千萬不能直接補零到2的次方。

正確的循環摺積計算結果:
            circular
(a₀ a₁ a₂) convolution (b₀ b₁ b₂) = (c₀ c₁ c₂)

c₀ = a₀b₀ + a₁b₂ + a₂b₁
c₁ = a₀b₁ + a₁b₀ + a₂b₂
c₂ = a₀b₂ + a₁b₁ + a₂b₀

長度補到2的次方,計算結果完全不同:
              circular
(a₀ a₁ a₂ 0) convolution (b₀ b₁ b₂ 0) = (d₀ d₁ d₂ d₃)

d₀ = a₀b₀ + a₂b₂
d₁ = a₀b₁ + a₁b₀
d₂ = a₀b₂ + a₁b₁ + a₂b₀
d₃ = a₁b₂ + a₂b₁

正確的方式:先補零直到不受循環影響,再補零直到長度是2的次方,最後讓輸出數列循環。

想要計算      circular
  (a₀ a₁ a₂) convolution (b₀ b₁ b₂) = (c₀ c₁ c₂)

改為計算          circular
 (a₀ a₁ a₂ 0 0) convolution (b₀ b₁ b₂ 0 0) = (d₀ d₁ d₂ d₃ d₄)

長度補到2的次方         circular
 (a₀ a₁ a₂ 0 0 0 0 0) convolution (b₀ b₁ b₂ 0 0 0 0 0)
= (d₀ d₁ d₂ d₃ d₄ 0 0 0)

最後讓輸出數列循環
  c₀ = d₀ + d₃
  c₁ = d₁ + d₄
  c₂ = d₂

Convolution快速演算法
Polynomial Multiplication快速演算法

運用循環摺積,計算摺積。

想要計算
  (a₀ a₁ a₂ a₃) convolution (b₀ b₁ b₂)
= (c₀ c₁ c₂ c₃ c₄ c₅)

改為計算             circular
  (a₀ a₁ a₂ a₃ 0 0) convolution (b₀ b₁ b₂ 0 0 0)
= (c₀ c₁ c₂ c₃ c₄ c₅ 0)

長度補到2的次方          circular
  (a₀ a₁ a₂ a₃ 0 0 0 0) convolution (b₀ b₁ b₂ 0 0 0 0 0)
= (c₀ c₁ c₂ c₃ c₄ c₅ 0 0)

截斷輸出數列至正確長度
  (c₀ c₁ c₂ c₃ c₄ c₅ 0 0) -> (c₀ c₁ c₂ c₃ c₄ c₅)

範例:大數乘法

大數乘法即是多項式乘法!

傅立葉轉換、數論轉換得以計算大數乘法,時間複雜度從O(N²)降為O(NlogN)。

範例:兩兩和(Pairwise Sum)(X+Y Problem)

甲集合、乙集合,只有整數。甲取一個數,乙取一個數,相加,會是哪些數?

集合資料結構是循序儲存:窮舉。O(N²)。

集合資料結構是索引儲存:摺積。O(RlogR)。R為數字範圍。

多項式觀點:整數是指數,該整數出現與否(未出現0,出現>0)是係數,兩兩和是多項式乘法。

摺積觀點:整數是索引值,相加是索引值位移,兩兩和是摺積。

範例:樣式匹配(Pattern Matching)

經典問題「一條陣列,尋找一個值」,解法是循序搜尋、先排序再二分搜尋。

進階問題「一條陣列,尋找一串連續數列」,解法是字串匹配、傅立葉轉換。此處討論傅立葉轉換。

兩串數列 a b
a = (1 2 3 4)
b = (1 3 5 7)

定義兩串數列 a b 的差距:對應項的差的平方的總和。
(1 - 1)² + (2 - 3)² + (3 - 5)² + (4 - 7)²

定義兩串數列 a b 的差距:對應項的差的平方的總和。
sum (a[i] - b[i])²

當差距等於零,則兩串數列相同。
sum (a[i] - b[i])² = 0  ---->  a = b

展開之後,重新整理,以數列為主角:數列每項平方和、數列點積。
  sum (a[i] - b[i])²
= sum (a[i]² + b[i]² - 2 a[i] b[i])
= sum a[i]² + sum b[i]² - 2 sum a[i] b[i]

縮寫。
(a - b)² = a² + b² - 2 (a · b)
兩串數列
(1 2 3 4 5 6)
(3 4 5)

匹配:窮舉各種對齊方式,判斷是否相符。
(1 2 3 4 5 6)  (1 2 3 4 5 6)  (1 2 3 4 5 6)  (1 2 3 4 5 6)
 | | |            | | |            | | |            | | |
(3 4 5)          (3 4 5)          (3 4 5)          (3 4 5)

換句話說,判斷差距是否等於零。
        (a - b)²              a²         b²     2(a·b)
(1-3)² + (2-4)² + (3-5)² = 1²+2²+3² + 3²+4²+5² - 2×26 = 12
(2-3)² + (3-4)² + (4-5)² = 2²+3²+4² + 3²+4²+5² - 2×38 = 3
(3-3)² + (4-4)² + (5-5)² = 3²+4²+5² + 3²+4²+5² - 2×50 = 0  <--- here!
(4-3)² + (5-4)² + (6-5)² = 4²+5²+6² + 3²+4²+5² - 2×62 = 5
                           ^^^^^^^^                ^^
                           moving window           convolution

預先計算每個數字的平方、預先計算各種對齊方式的點積(就是摺積),就可以快速求得a² + b² - 2ab。時間複雜度O(NlogN)。

匹配(找到零)、最相似(找最小值)、最相異(找最大值)。

UVa 12298 ICPC 4671 5705 7159