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⁻¹

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

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

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

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

Filter

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ᵢ⁻¹ᵀ