Algoritmul Levinson-Durbin

Algoritmul Levinson-Durbin , în algebră liniară , este utilizat pentru a calcula, cu o metodă recursivă , soluția unei ecuații care implică o matrice Toeplitz . Algoritmul rulează $\Theta (n^{2})$ ${\ displaystyle \ Theta (n ^ {2})}$ $\ Theta (n ^ {2})$ pași, ceea ce reprezintă o îmbunătățire majoră față de metoda Gaussian Elimination , care necesită $\Theta (n^{3})$ ${\ displaystyle \ Theta (n ^ {3})}$ ${\ displaystyle \ Theta (n ^ {3})}$ pași.

Algoritmul Levinson - Durbin a fost propus pentru prima dată de Norman Levinson , în 1947, și îmbunătățit de James Durbin în 1960; ulterior a fost îmbunătățit și mai mult, aducându-l din $4n^{2}$ ${\ displaystyle 4n ^ {2}}$ ${\ displaystyle 4n ^ {2}}$ până la $3n^{2}$ ${\ displaystyle 3n ^ {2}}$ ${\ displaystyle 3n ^ {2}}$ multiplicări, de la WF Trench și, respectiv, S. Zohar.

Alte metode de prelucrare a datelor includ descompunerea Schur și descompunerea Cholesky . În comparație cu acestea, recursiunea Levinson (în special recursivitatea Levinson divizată) tinde să fie mai rapidă din punct de vedere al calculului, dar mai sensibilă la inexactitățile de calcul, cum ar fi erorile de rotunjire .

Algoritmul Bareiss pentru matricile Toeplitz (nu trebuie confundat cu algoritmul general Bareiss) este la fel de rapid ca recursiunea Levinson-Durbin, dar folosește $\Theta (n^{2})$ ${\ displaystyle \ Theta (n ^ {2})}$ $\ Theta (n ^ {2})$ pași, în timp ce recurența Levinson - Durbin necesită numai $\Theta (n)$ ${\ displaystyle \ Theta (n)}$ $\ Theta (n)$ pași. Totuși, algoritmul Bareiss este stabil din punct de vedere numeric , ^[1] ^{[2] în} timp ce recursiunea Levinson-Durbin, în cel mai bun caz, este doar slab stabilă (adică arată stabilitate numerică pentru sistemele liniare bine condiționate ). ^[3]

Algoritmi mai noi, numiți asimptotic rapid sau, în unele texte, algoritmi Toeplitz super- rapizi, pot rezolva problema în $\Theta (n\log _{p}n)$ ${\ displaystyle \ Theta (n \ log _ {p} n)}$ ${\ displaystyle \ Theta (n \ log _ {p} n)}$ pentru diverse $p$ ${\ displaystyle p}$ $p$ (ex. $p=2$ ${\ displaystyle p = 2}$ $p = 2$ , ^[4] ^[5] $p=3$ ${\ displaystyle p = 3}$ $p = 3$ ^[5] ). Recursiunea Levinson-Durbin rămâne populară din mai multe motive; în primul rând, este relativ ușor de înțeles; în plus, poate fi mai rapid decât un algoritm super rapid pentru bebeluși $n$ ${\ displaystyle n}$ $n$ (obișnuit $n<256$ ${\ displaystyle n <256}$ ${\ displaystyle n <256}$ ). ^[6]

Derivare

Introducere

Ecuațiile matriciale au următoarea formă:

\mathbf {M} \ {\vec {x}}={\vec {y}}.

{\ displaystyle \ mathbf {M} \ {\ vec {x}} = {\ vec {y}}.}

{\ displaystyle \ mathbf {M} \ {\ vec {x}} = {\ vec {y}}.}

Algoritmul Levinson-Durbin poate fi utilizat pentru orice ecuație, atâta timp cât $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ este o matrice Toeplitz cunoscută, cu diagonală principală diferită de zero; unde este ${\vec {y}}$ ${\ displaystyle {\ vec {y}}}$ ${\ displaystyle {\ vec {y}}}$ este vectorul cunoscut, în timp ce ${\vec {x}}$ ${\ displaystyle {\ vec {x}}}$ ${\ vec x}$ este vectorul necunoscut al numerelor $x_{i}$ ${\ displaystyle x_ {i}}$ $x_i$ a fi determinat.

Considera ${\hat {e}}_{i}$ ${\ displaystyle {\ hat {e}} _ {i}}$ ${\ displaystyle {\ hat {e}} _ {i}}$ un vector compus în întregime din zerouri, cu excepția celui de-al i- lea termen, care conține o valoare unitară. Lungimea acestuia va fi implicit determinată de context. Termenul $Nu.$ ${\ displaystyle N}$ $Nu.$ Se referă la lățimea matricei $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ având dimensiuni $N\times N$ ${\ displaystyle N \ times N}$ $N \ ori N$ . În cele din urmă, indicii superiori se referă la un indice inductiv , în timp ce indicii indică indicii. De exemplu (și definiție) matricea $\mathbf {T} ^{n}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ este o matrice $n\times n$ ${\ displaystyle n \ times n}$ $n \ ori n$ care copiază blocul $n\times n$ ${\ displaystyle n \ times n}$ $n \ ori n$ sus în stânga $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ , adică $T_{i,j}^{n}=M_{j}^{i}$ ${\ displaystyle T_ {i, j} ^ {n} = M_ {j} ^ {i}}$ ${\ displaystyle T_ {i, j} ^ {n} = M_ {j} ^ {i}}$ .

De aceea, de asemenea $\mathbf {T} ^{n}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ este o matrice Toeplitz; în sensul că poate fi scris în următoarea formă:

\mathbf {T} ^{n}={\begin{bmatrix}t_{0}&t_{-1}&t_{-2}&\dots &t_{-n+1}\\t_{1}&t_{0}&t_{-1}&\dots &t_{-n+2}\\t_{2}&t_{1}&t_{0}&\dots &t_{-n+3}\\\vdots &\vdots &\vdots &\ddots &\vdots \\t_{n-1}&t_{n-2}&t_{n-3}&\dots &t_{0}\end{bmatrix}}.

{\ displaystyle \ mathbf {T} ^ {n} = {\ begin {bmatrix} t_ {0} & t _ {- 1} & t _ {- 2} & \ dots & t _ {- n + 1} \ \ t_ {1} & t_ {0} & t _ {- 1} & \ dots & t _ {- n + 2} \\ t_ {2} & t_ {1} & t_ {0} & \ dots & t _ {- n + 3} \\\ vdots & \ vdots & \ vdots & \ ddots & \ vdots \\ t_ {n-1} & t_ {n-2} & t_ {n-3} & \ dots & t_ {0} \ end {bmatrix}}.}

{\ displaystyle \ mathbf {T} ^ {n} = {\ begin {bmatrix} t_ {0} & t _ {- 1} & t _ {- 2} & \ dots & t _ {- n + 1} \ \ t_ {1} & t_ {0} & t _ {- 1} & \ dots & t _ {- n + 2} \\ t_ {2} & t_ {1} & t_ {0} & \ dots & t _ {- n + 3} \\\ vdots & \ vdots & \ vdots & \ ddots & \ vdots \\ t_ {n-1} & t_ {n-2} & t_ {n-3} & \ dots & t_ {0} \ end {bmatrix}}.}

Etape introductive

Algoritmul continuă urmând doi pași. În primul pas, se stabilesc două grupuri de vectori, numiți vectori înainte și invers . Vectorii înainte sunt folosiți pentru a ajuta la obținerea setului de vectori înapoi; prin urmare, ele pot fi imediat aruncate. În schimb, sunt necesari vectori înapoi pentru al doilea pas, unde sunt utilizați pentru a crea soluția dorită.

Recursiunea Levinson-Durbin definește încă un „vector înainte”, numit ${\vec {f}}^{n}$ ${\ displaystyle {\ vec {f}} ^ {n}}$ ${\ displaystyle {\ vec {f}} ^ {n}}$ , un vector de lungime n care satisface ecuația:

\mathbf {T} ^{n}{\vec {f}}^{n}={\hat {e}}_{1}.

{\ displaystyle \ mathbf {T} ^ {n} {\ vec {f}} ^ {n} = {\ hat {e}} _ {1}.}

{\ displaystyle \ mathbf {T} ^ {n} {\ vec {f}} ^ {n} = {\ hat {e}} _ {1}.}

Încă un „vector înapoi”, numit ${\vec {b}}^{n}$ ${\ displaystyle {\ vec {b}} ^ {n}}$ ${\ displaystyle {\ vec {b}} ^ {n}}$ este definit într-un mod similar; este vectorul de lungime n care satisface ecuația:

\mathbf {T} ^{n}{\vec {b}}^{n}={\hat {e}}_{n}.

{\ displaystyle \ mathbf {T} ^ {n} {\ vec {b}} ^ {n} = {\ hat {e}} _ {n}.}

{\ displaystyle \ mathbf {T} ^ {n} {\ vec {b}} ^ {n} = {\ hat {e}} _ {n}.}

O simplificare importantă poate apărea atunci când $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ este o matrice simetrică ; în acest caz cei doi vectori sunt corelați de $b_{i}^{n}=f_{n+1-i}^{n}$ ${\ displaystyle b_ {i} ^ {n} = f_ {n + 1-i} ^ {n}}$ ${\ displaystyle b_ {i} ^ {n} = f_ {n + 1-i} ^ {n}}$ , adică sunt inversări de linie una de alta. Acest lucru poate salva unele calcule în acest caz particular.

Recuperați transportatorii

Chiar dacă matricea nu este simetrică, încă un vector înainte și înapoi poate fi găsit din vectorii de lungime $n-1$ ${\ displaystyle n-1}$ $n-1$ după cum urmează. În primul rând, vectorul din spate poate fi „extins” cu un zero pentru a obține:

\mathbf {T} ^{n}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\\\end{bmatrix}}={\begin{bmatrix}\ &\ &\ &t_{-n+1}\\\ &\mathbf {T} ^{n-1}&\ &t_{-n+2}\\\ &\ &\ &\vdots \\t_{n-1}&t_{n-2}&\dots &t_{0}\\\end{bmatrix}}{\begin{bmatrix}\ \\{\vec {f}}^{n-1}\\\ \\0\\\ \\\end{bmatrix}}={\begin{bmatrix}1\\0\\\vdots \\0\\\epsilon _{f}^{n}\end{bmatrix}}.

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \\\ end {bmatrix}} = {\ begin {bmatrix} \ & \ & \ & t _ {- n + 1} \\\ & \ mathbf {T} ^ {n-1} & \ & t _ {- n + 2} \\\ & \ & \ & \ vdots \ \ t_ {n -1} & t_ {n-2} & \ dots & t_ {0} \\\ end {bmatrix}} {\ begin {bmatrix} \ \\ {\ vec {f}} ^ {n- 1} \\\ \ \ 0 \\\ \\\ end {bmatrix}} = {\ begin {bmatrix} 1 \\ 0 \\\ vdots \\ 0 \\\ epsilon _ {f} ^ {n} \ end {bmatrix}}.}

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \\\ end {bmatrix}} = {\ begin {bmatrix} \ & \ & \ & t _ {- n + 1} \\\ & \ mathbf {T} ^ {n-1} & \ & t _ {- n + 2} \\\ & \ & \ & \ vdots \ \ t_ {n -1} & t_ {n-2} & \ dots & t_ {0} \\\ end {bmatrix}} {\ begin {bmatrix} \ \\ {\ vec {f}} ^ {n- 1} \\\ \ \ 0 \\\ \\\ end {bmatrix}} = {\ begin {bmatrix} 1 \\ 0 \\\ vdots \\ 0 \\\ epsilon _ {f} ^ {n} \ end {bmatrix}}.}

Plecând de la $\mathbf {T} ^{n-1}$ ${\ displaystyle \ mathbf {T} ^ {n-1}}$ ${\ displaystyle \ mathbf {T} ^ {n-1}}$ la $\mathbf {T} ^{n}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ ${\ displaystyle \ mathbf {T} ^ {n}}$ , coloana suplimentară adăugată la matrice nu deranjează soluția atunci când se utilizează un zero pentru a extinde vectorul înainte. Cu toate acestea, linia suplimentară adăugată la matrice a perturbat soluția; și a creat un termen de eroare nedorit ε _f care apare în cele din urmă. Ecuația anterioară dă valoarea:

\epsilon _{f}^{n}\ =\ \sum _{i=1}^{n-1}\ M_{ni}\ f_{i}^{n-1}\ =\ \sum _{i=1}^{n-1}\ t_{n-i}\ f_{i}^{n-1}.

{\ displaystyle \ epsilon _ {f} ^ {n} \ = \ \ sum _ {i = 1} ^ {n-1} \ M_ {ni} \ f_ {i} ^ {n-1} \ = \ \ sum _ {i = 1} ^ {n-1} \ t_ {ni} \ f_ {i} ^ {n-1}.}

{\ displaystyle \ epsilon _ {f} ^ {n} \ = \ \ sum _ {i = 1} ^ {n-1} \ M_ {ni} \ f_ {i} ^ {n-1} \ = \ \ sum _ {i = 1} ^ {n-1} \ t_ {ni} \ f_ {i} ^ {n-1}.}

Această eroare va fi returnată în scurt timp și eliminată de noul vector forward; dar mai întâi, vectorul înapoi trebuie extins într-un mod similar (deși inversat). Pentru vectorul înapoi avem:

\mathbf {T} ^{n}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\\\end{bmatrix}}={\begin{bmatrix}t_{0}&\dots &t_{-n+2}&t_{-n+1}\\\vdots &\ &\ &\ \\t_{n-2}&\ &\mathbf {T} ^{n-1}&\ \\t_{n-1}&\ &\ &\end{bmatrix}}{\begin{bmatrix}\ \\0\\\ \\{\vec {b}}^{n-1}\\\ \\\end{bmatrix}}={\begin{bmatrix}\epsilon _{b}^{n}\\0\\\vdots \\0\\1\end{bmatrix}}.

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \\\ end {bmatrix}} = {\ begin {bmatrix} t_ {0} & \ dots & t _ {- n + 2} & t _ {- n + 1} \\\ vdots & \ & \ & \ \\ t_ {n-2} & \ & \ mathbf {T} ^ {n- 1} & \ \\ t_ {n-1} & \ & \ & \ end {bmatrix}} {\ begin {bmatrix} \ \\ 0 \\\ \\ {\ vec {b}} ^ {n-1} \\\ \\\ end {bmatrix}} = {\ begin {bmatrix} \ epsilon _ {b} ^ {n} \\ 0 \\\ vdots \\ 0 \\ 1 \ end {bmatrix }}.}

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \\\ end {bmatrix}} = {\ begin {bmatrix} t_ {0} & \ dots & t _ {- n + 2} & t _ {- n + 1} \\\ vdots & \ & \ & \ \\ t_ {n-2} & \ & \ mathbf {T} ^ {n- 1} & \ \\ t_ {n-1} & \ & \ & \ end {bmatrix}} {\ begin {bmatrix} \ \\ 0 \\\ \\ {\ vec {b}} ^ {n-1} \\\ \\\ end {bmatrix}} = {\ begin {bmatrix} \ epsilon _ {b} ^ {n} \\ 0 \\\ vdots \\ 0 \\ 1 \ end {bmatrix }}.}

La fel ca înainte, coloana suplimentară la matrice nu deranjează acest nou vector înapoi; dar linia suplimentară o face. Iată o altă eroare nedorită $\epsilon$ ${\ displaystyle \ epsilon}$ $\ epsilon$ egal cu:

\epsilon _{b}^{n}\ =\ \sum _{i=2}^{n}\ M_{1i}\ b_{i-1}^{n-1}\ =\ \sum _{i=1}^{n-1}\ t_{-i}\ b_{i}^{n-1}.\

{\ displaystyle \ epsilon _ {b} ^ {n} \ = \ \ sum _ {i = 2} ^ {n} \ M_ {1i} \ b_ {i-1} ^ {n-1} \ = \ \ sum _ {i = 1} ^ {n-1} \ t _ {- i} \ b_ {i} ^ {n-1}. \}

{\ displaystyle \ epsilon _ {b} ^ {n} \ = \ \ sum _ {i = 2} ^ {n} \ M_ {1i} \ b_ {i-1} ^ {n-1} \ = \ \ sum _ {i = 1} ^ {n-1} \ t _ {- i} \ b_ {i} ^ {n-1}. \}

Acești doi termeni de eroare pot fi folosiți pentru a forma vectori de ordine superioară înainte și înapoi descriși după cum urmează. Folosind liniaritatea matricilor, identitatea următoare este valabilă pentru toți $(\alpha ,\beta )$ ${\ displaystyle (\ alpha, \ beta)}$ $(\ alfa, \ beta)$ :

\mathbf {T} \left(\alpha {\begin{bmatrix}{\vec {f}}\\\ \\0\\\end{bmatrix}}+\beta {\begin{bmatrix}0\\\ \\{\vec {b}}\end{bmatrix}}\right)=\alpha {\begin{bmatrix}1\\0\\\vdots \\0\\\epsilon _{f}\\\end{bmatrix}}+\beta {\begin{bmatrix}\epsilon _{b}\\0\\\vdots \\0\\1\end{bmatrix}}.

{\ displaystyle \ mathbf {T} \ left (\ alpha {\ begin {bmatrix} {\ vec {f}} \\\ \\ 0 \\\ end {bmatrix}} + \ beta {\ begin {bmatrix} 0 \\\ \\ {\ vec {b}} \ end {bmatrix}} \ right) = \ alpha {\ begin {bmatrix} 1 \\ 0 \\\ vdots \\ 0 \\\ epsilon _ {f} \ \\ end {bmatrix}} + \ beta {\ begin {bmatrix} \ epsilon _ {b} \\ 0 \\\ vdots \\ 0 \\ 1 \ end {bmatrix}}.}

{\ displaystyle \ mathbf {T} \ left (\ alpha {\ begin {bmatrix} {\ vec {f}} \\\ \\ 0 \\\ end {bmatrix}} + \ beta {\ begin {bmatrix} 0 \\\ \\ {\ vec {b}} \ end {bmatrix}} \ right) = \ alpha {\ begin {bmatrix} 1 \\ 0 \\\ vdots \\ 0 \\\ epsilon _ {f} \ \\ end {bmatrix}} + \ beta {\ begin {bmatrix} \ epsilon _ {b} \\ 0 \\\ vdots \\ 0 \\ 1 \ end {bmatrix}}.}

De sine $\alpha$ ${\ displaystyle \ alpha}$ $\ alfa$ Și $\beta$ ${\ displaystyle \ beta}$ $\ beta$ sunt aleși dreptaci ${\hat {e}}_{1}$ ${\ displaystyle {\ hat {e}} _ {1}}$ ${\ displaystyle {\ hat {e}} _ {1}}$ sau ${\hat {e}}_{n,n}$ ${\ displaystyle {\ hat {e}} _ {n, n}}$ ${\ displaystyle {\ hat {e}} _ {n, n}}$ _, atunci cantitatea dintre paranteze va fi egală cu definiția celui de-al n-lea vector înainte sau înapoi, respectiv. Cu alegerea termenilor $\alpha$ ${\ displaystyle \ alpha}$ $\ alfa$ Și $\beta$ ${\ displaystyle \ beta}$ $\ beta$ , suma vectorilor din paranteze este simplă și produce rezultatul dorit.

Pentru a găsi acești coeficienți, $\alpha _{f}^{n}$ ${\ displaystyle \ alpha _ {f} ^ {n}}$ ${\ displaystyle \ alpha _ {f} ^ {n}}$ , $\beta _{f}^{n}$ ${\ displaystyle \ beta _ {f} ^ {n}}$ ${\ displaystyle \ beta _ {f} ^ {n}}$ trebuie să fie astfel încât:

{\vec {f}}_{n}=\alpha _{f}^{n}{\begin{bmatrix}{\vec {f}}_{n-1}\\0\end{bmatrix}}+\beta _{f}^{n}{\begin{bmatrix}0\\{\vec {b}}_{n-1}\end{bmatrix}}

{\ displaystyle {\ vec {f}} _ {n} = \ alpha _ {f} ^ {n} {\ begin {bmatrix} {\ vec {f}} _ {n-1} \\ 0 \ end { bmatrix}} + \ beta _ {f} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} _ {n-1} \ end {bmatrix}}}

{\ displaystyle {\ vec {f}} _ {n} = \ alpha _ {f} ^ {n} {\ begin {bmatrix} {\ vec {f}} _ {n-1} \\ 0 \ end { bmatrix}} + \ beta _ {f} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} _ {n-1} \ end {bmatrix}}}

și, respectiv, $\alpha _{b}^{n}$ ${\ displaystyle \ alpha _ {b} ^ {n}}$ ${\ displaystyle \ alpha _ {b} ^ {n}}$ , $\beta _{b}^{n}$ ${\ displaystyle \ beta _ {b} ^ {n}}$ ${\ displaystyle \ beta _ {b} ^ {n}}$ sunt astfel încât:

{\vec {b}}_{n}=\alpha _{b}^{n}{\begin{bmatrix}{\vec {f}}_{n-1}\\0\end{bmatrix}}+\beta _{b}^{n}{\begin{bmatrix}0\\{\vec {b}}_{n-1}\end{bmatrix}}.

{\ displaystyle {\ vec {b}} _ {n} = \ alpha _ {b} ^ {n} {\ begin {bmatrix} {\ vec {f}} _ {n-1} \\ 0 \ end { bmatrix}} + \ beta _ {b} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} _ {n-1} \ end {bmatrix}}.}

{\ displaystyle {\ vec {b}} _ {n} = \ alpha _ {b} ^ {n} {\ begin {bmatrix} {\ vec {f}} _ {n-1} \\ 0 \ end { bmatrix}} + \ beta _ {b} ^ {n} {\ begin {bmatrix} 0 \\ {\ vec {b}} _ {n-1} \ end {bmatrix}}.}

Înmulțind și împărțind ecuațiile de mai sus cu ${\mathbf {T} }^{n}$ ${\ displaystyle {\ mathbf {T}} ^ {n}}$ ${\ displaystyle {\ mathbf {T}} ^ {n}}$ obținem următoarea ecuație:

{\begin{bmatrix}1&\epsilon _{b}^{n}\\0&0\\\vdots &\vdots \\0&0\\\epsilon _{f}^{n}&1\end{bmatrix}}{\begin{bmatrix}\alpha _{f}^{n}&\alpha _{b}^{n}\\\beta _{f}^{n}&\beta _{b}^{n}\end{bmatrix}}={\begin{bmatrix}1&0\\0&0\\\vdots &\vdots \\0&0\\0&1\end{bmatrix}}.

{\ displaystyle {\ begin {bmatrix} 1 & \ epsilon _ {b} ^ {n} \\ 0 & 0 \\\ vdots & \ vdots \\ 0 & 0 \\\ epsilon _ {f} ^ {n} & 1 \ end {bmatrix}} {\ begin {bmatrix} \ alpha _ {f} ^ {n} & \ alpha _ {b} ^ {n} \\\ beta _ {f} ^ {n} & \ beta _ {b} ^ {n} \ end {bmatrix}} = {\ begin {bmatrix} 1 & 0 \\ 0 & 0 \\\ vdots & \ vdots \\ 0 & 0 \\ 0 & 1 \ end {bmatrix }}.}

{\ displaystyle {\ begin {bmatrix} 1 & \ epsilon _ {b} ^ {n} \\ 0 & 0 \\\ vdots & \ vdots \\ 0 & 0 \\\ epsilon _ {f} ^ {n} & 1 \ end {bmatrix}} {\ begin {bmatrix} \ alpha _ {f} ^ {n} & \ alpha _ {b} ^ {n} \\\ beta _ {f} ^ {n} & \ beta _ {b} ^ {n} \ end {bmatrix}} = {\ begin {bmatrix} 1 & 0 \\ 0 & 0 \\\ vdots & \ vdots \\ 0 & 0 \\ 0 & 1 \ end {bmatrix }}.}

Acum, toate zerourile celor doi vectori de mai sus sunt ignorate, rămâne doar următoarea ecuație:

{\begin{bmatrix}1&\epsilon _{b}^{n}\\\epsilon _{f}^{n}&1\end{bmatrix}}{\begin{bmatrix}\alpha _{f}^{n}&\alpha _{b}^{n}\\\beta _{f}^{n}&\beta _{b}^{n}\end{bmatrix}}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}.

{\ displaystyle {\ begin {bmatrix} 1 & \ epsilon _ {b} ^ {n} \\\ epsilon _ {f} ^ {n} & 1 \ end {bmatrix}} {\ begin {bmatrix} \ alpha _ {f} ^ {n} & \ alpha _ {b} ^ {n} \\\ beta _ {f} ^ {n} & \ beta _ {b} ^ {n} \ end {bmatrix}} = {\ începe {bmatrix} 1 și 0 \\ 0 & 1 \ end {bmatrix}}.}

{\ displaystyle {\ begin {bmatrix} 1 & \ epsilon _ {b} ^ {n} \\\ epsilon _ {f} ^ {n} & 1 \ end {bmatrix}} {\ begin {bmatrix} \ alpha _ {f} ^ {n} & \ alpha _ {b} ^ {n} \\\ beta _ {f} ^ {n} & \ beta _ {b} ^ {n} \ end {bmatrix}} = {\ începe {bmatrix} 1 și 0 \\ 0 & 1 \ end {bmatrix}}.}

Cu aceste soluții (folosind formula inversă a matricei Cramer $2\times 2$ ${\ displaystyle 2 \ times 2}$ $2 \ ori 2$ ), noii vectori înainte și invers sunt:

{\vec {f}}^{n}={1 \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}-{\epsilon _{f}^{n} \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}

{\ displaystyle {\ vec {f}} ^ {n} = {1 \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \ end {bmatrix}} - {\ epsilon _ {f} ^ {n} \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \ end {bmatrix}}}

{\ displaystyle {\ vec {f}} ^ {n} = {1 \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \ end {bmatrix}} - {\ epsilon _ {f} ^ {n} \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \ end {bmatrix}}}

{\vec {b}}^{n}={1 \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}0\\{\vec {b}}^{n-1}\end{bmatrix}}-{\epsilon _{b}^{n} \over {1-\epsilon _{b}^{n}\epsilon _{f}^{n}}}{\begin{bmatrix}{\vec {f}}^{n-1}\\0\end{bmatrix}}.

{\ displaystyle {\ vec {b}} ^ {n} = {1 \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \ end {bmatrix}} - {\ epsilon _ {b} ^ {n} \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \ end {bmatrix}}.}

{\ displaystyle {\ vec {b}} ^ {n} = {1 \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} 0 \\ {\ vec {b}} ^ {n-1} \ end {bmatrix}} - {\ epsilon _ {b} ^ {n} \ over {1- \ epsilon _ {b} ^ {n} \ epsilon _ {f} ^ {n}}} {\ begin {bmatrix} {\ vec {f}} ^ {n-1} \\ 0 \ end {bmatrix}}.}

Execuția acestor însumări vectoriale, prin urmare, dă al n-lea vector înainte și înapoi începând cu cele anterioare. Rămâne doar să găsim primul dintre acești vectori, apoi cu sume rapide și multiplicări rapide obținem termenii rămași. Primii vectori înainte și înapoi sunt pur și simplu:

{\vec {f}}^{1}={\vec {b}}^{1}=\left[{1 \over M_{11}}\right]=\left[{1 \over t_{0}}\right].

{\ displaystyle {\ vec {f}} ^ {1} = {\ vec {b}} ^ {1} = \ left [{1 \ over M_ {11}} \ right] = \ left [{1 \ over t_ {0}} \ dreapta].}

{\ displaystyle {\ vec {f}} ^ {1} = {\ vec {b}} ^ {1} = \ left [{1 \ over M_ {11}} \ right] = \ left [{1 \ over t_ {0}} \ dreapta].}

Folosind vectori înapoi

Pașii anteriori dau N vectori înapoi pentru $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ . De acolo, o ecuație mai arbitrară este:

{\vec {y}}=\mathbf {M} \ {\vec {x}}.

{\ displaystyle {\ vec {y}} = \ mathbf {M} \ {\ vec {x}}.}

{\ displaystyle {\ vec {y}} = \ mathbf {M} \ {\ vec {x}}.}

Soluția poate fi construită în același mod recursiv în care au fost construiți vectorii înapoi. În consecință, ${\vec {x}}$ ${\ displaystyle {\ vec {x}}}$ ${\ vec x}$ trebuie generalizată la o succesiune de ${\vec {x}}^{n}$ ${\ displaystyle {\ vec {x}} ^ {n}}$ ${\ displaystyle {\ vec {x}} ^ {n}}$ intermediar, astfel încât ${\vec {x}}^{N}={\vec {x}}$ ${\ displaystyle {\ vec {x}} ^ {N} = {\ vec {x}}}$ ${\ displaystyle {\ vec {x}} ^ {N} = {\ vec {x}}}$ .

Soluția este apoi construită recursiv observând că, dacă

\mathbf {T} ^{n-1}{\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\vdots \\x_{n-1}^{n-1}\\\end{bmatrix}}={\begin{bmatrix}y_{1}\\y_{2}\\\vdots \\y_{n-1}\end{bmatrix}},

{\ displaystyle \ mathbf {T} ^ {n-1} {\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ { n-1} ^ {n-1} \\\ end {bmatrix}} = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \ end { bmatrix}},}

{\ displaystyle \ mathbf {T} ^ {n-1} {\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ { n-1} ^ {n-1} \\\ end {bmatrix}} = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \ end { bmatrix}},}

apoi, inserând din nou un zero la sfârșitul vectorului și definind o constantă de eroare acolo unde este necesar, avem:

\mathbf {T} ^{n}{\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\vdots \\x_{n-1}^{n-1}\\0\end{bmatrix}}={\begin{bmatrix}y_{1}\\y_{2}\\\vdots \\y_{n-1}\\\epsilon _{x}^{n-1}\end{bmatrix}}.

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ {n- 1} ^ {n-1} \\ 0 \ end {bmatrix}} = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \\\ epsilon _ {x} ^ {n-1} \ end {bmatrix}}.}

{\ displaystyle \ mathbf {T} ^ {n} {\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ {n- 1} ^ {n-1} \\ 0 \ end {bmatrix}} = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \\\ epsilon _ {x} ^ {n-1} \ end {bmatrix}}.}

Putem folosi apoi al n-lea vector înapoi pentru a elimina eroarea și a o înlocui cu formula dorită după cum urmează:

\mathbf {T} ^{n}\left({\begin{bmatrix}x_{1}^{n-1}\\x_{2}^{n-1}\\\vdots \\x_{n-1}^{n-1}\\0\\\end{bmatrix}}+(y_{n}-\epsilon _{x}^{n-1})\ {\vec {b}}^{n}\right)={\begin{bmatrix}y_{1}\\y_{2}\\\vdots \\y_{n-1}\\y_{n}\end{bmatrix}}.

{\ displaystyle \ mathbf {T} ^ {n} \ left ({\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ {n-1} ^ {n-1} \\ 0 \\\ end {bmatrix}} + (y_ {n} - \ epsilon _ {x} ^ {n-1}) \ {\ vec {b}} ^ {n} \ right) = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \\ y_ {n} \ end {bmatrix}}.}

{\ displaystyle \ mathbf {T} ^ {n} \ left ({\ begin {bmatrix} x_ {1} ^ {n-1} \\ x_ {2} ^ {n-1} \\\ vdots \\ x_ {n-1} ^ {n-1} \\ 0 \\\ end {bmatrix}} + (y_ {n} - \ epsilon _ {x} ^ {n-1}) \ {\ vec {b}} ^ {n} \ right) = {\ begin {bmatrix} y_ {1} \\ y_ {2} \\\ vdots \\ y_ {n-1} \\ y_ {n} \ end {bmatrix}}.}

Extinderea acestei metode până la $n=N$ ${\ displaystyle n = N}$ ${\ displaystyle n = N}$ produce soluția dorită ${\vec {x}}$ ${\ displaystyle {\ vec {x}}}$ ${\ vec x}$ .

În practică, acești pași sunt adesea realizați concomitent cu restul procedurii, formând o unitate coerentă și merită să fie tratați individual.

Algoritmul Levinson-Durbin

De sine $\mathbf {M}$ ${\ displaystyle \ mathbf {M}}$ $\ mathbf M$ nu o matrice Toeplitz în sens strict, ci o matrice bloc Toeplitz, recursiunea Levinson poate fi derivată în același mod având în vedere submatricea Toeplitz (Musicus 1988).

Tablourile cu blocuri Toeplitz apar în mod natural în algoritmi de procesare a semnalului atunci când se tratează mai multe fluxuri de semnal (de exemplu în sistemele MIMO) sau semnale ciclostationare.

Aplicarea practică a algoritmului Levinson-Durbin

Algoritmul Levinson-Durbin este utilizat pe scară largă pentru rezoluția modelelor autoregresive $AR(p)$ ${\ displaystyle AR (p)}$ ${\ displaystyle AR (p)}$ de ordine $p$ ${\ displaystyle p}$ $p$ (utilizat în protocolul GSM ), care sunt prezentate sub următoarea ecuație a diferenței :

${\begin{aligned}t(n)+x_{1}t(n-1)+x_{2}t(n-2)+\cdots +x_{p}t(n-p)&=e(n)\\t(n)+\sum _{i=1}^{n}x_{i}t(n-i)&=e(n)\\\end{aligned}},$ ${\ displaystyle {\ begin {align} t (n) + x_ {1} t (n-1) + x_ {2} t (n-2) + \ cdots + x_ {p} t (np) & = e (n) \\ t (n) + \ sum _ {i = 1} ^ {n} x_ {i} t (ni) & = e (n) \\\ end {align}},}$ ${\ displaystyle {\ begin {align} t (n) + x_ {1} t (n-1) + x_ {2} t (n-2) + \ cdots + x_ {p} t (np) & = e (n) \\ t (n) + \ sum _ {i = 1} ^ {n} x_ {i} t (ni) & = e (n) \\\ end {align}},}$

unde este $t(n)$ ${\ displaystyle t (n)}$ ${\ displaystyle t (n)}$ este eșantionul actual al sistemului estimat de model $AR(p)$ ${\ displaystyle AR (p)}$ ${\ displaystyle AR (p)}$ , $x_{i}$ ${\ displaystyle x_ {i}}$ ${\ displaystyle x_ {i}}$ sunt $p$ ${\ displaystyle p}$ ${\ displaystyle p}$ parametrii modelului, care se aplică la $p$ ${\ displaystyle p}$ ${\ displaystyle p}$ mostre anterioare sau - în mod similar - definesc memoria modelului (în acest caz de comandă $p$ ${\ displaystyle p}$ ${\ displaystyle p}$ ). La sfarsit $e(n)$ ${\ displaystyle e (n)}$ ${\ displaystyle e (n)}$ este reziduul de predicție al probei $t(n)$ ${\ displaystyle t (n)}$ ${\ displaystyle t (n)}$ pe care doriți să o minimizați, adică: componenta imprevizibilă a sistemului.

Aceste ecuații sunt prezentate sub formă de matrice Toeplitz, prin urmare algoritmul Levinson-Durbin este utilizat pentru soluția lor.

Pseudocod pentru recursiunea Levinson-Durbin

Algoritmul se bazează pe calculul coeficienților autoregresivi ai matricilor de ordin crescător. Este împărțit în două faze: o primă inițializare pentru calcularea parametrului $x_{1}$ ${\ displaystyle x_ {1}}$ ${\ displaystyle x_ {1}}$ sau - care este același - pentru cazul elementar al matricei $1\times 1$ ${\ displaystyle 1 \ times 1}$ ${\ displaystyle 1 \ times 1}$ . Ulterior, vom continua cu calculul iterativ al parametrilor pentru matricile de ordine via, crescând treptat $2\times 2$ ${\ displaystyle 2 \ times 2}$ ${\ displaystyle 2 \ times 2}$ , $3\times 3$ ${\ displaystyle 3 \ times 3}$ ${\ displaystyle 3 \ times 3}$ , ..., $p\times p$ ${\ displaystyle p \ times p}$ ${\ displaystyle p \ times p}$ .

Folosind notația MATLAB / Octave , pseudocodul pentru calcularea recursiunii Levinson-Durbin este următorul ^[7] :

 k = M ( 2 ) / M ( 1 );                                    % Estimați primul element
X = k ;
ȘI = ( 1 - k2 ) * M ( 1 );                                % Calculează eroarea pătrată medie
pentru the = 2 : p
    k = ( M ( i + 1 ) - X * M ( 2 : i )) / Și ;                % Coeficienți de reflecție
    X = [ k , X - k * x ( i - 1 : - 1 : 1 )];                   % Estimați următoarele elemente
    ȘI = ( 1 - k2 ) * Și ;                               % Actualizează eroarea pătrată medie
Sfârșit

X = [ 1 , - x ( N : - 1 : 1 )];                               % Returnează vectorul necunoscut

Trebuie subliniat faptul că MATLAB funcționează cu vectori și matrici, prin urmare - dacă doriți să traduceți codul în limbi precum C ++ sau Java , veți obține două imbricate pentru bucle.

Notă

^ Adam W. Bojanczyk, Brent, RP, De Hoog, FR, & Sweet, DR, Despre stabilitatea Bareiss și algoritmi de factorizare Toeplitz înrudiți , în SIAM Journal on Matrix Analysis and Applications , vol. 16, n. 1, 1995, pp. 40-57.
^ Brent, Richard P., Stabilitatea algoritmilor rapide pentru sisteme liniare structurate , în Algoritmi fiabili rapid pentru matrici cu structură , Society for Industrial and Applied Mathematics, 1999, pp. 103-116.
^ Hari Krishna, Yunbiao Wang, Algoritmul Levinson divizat este slab stabil , în jurnalul SIAM privind analiza numerică , vol. 30, n. 5, 1993, pp. 1498-1508.
^ Bojanczyk, Adam W., Richard P. Brent și Frank R. De Hoog, Un algoritm slab stabil pentru sistemele generale Toeplitz , în Algoritmi numerici , vol. 1, arXiv: 1005.0503, 2010, pp. 225-244.
^ ^a ^b Stewart, Michael, Un solver Toeplitz super rapid cu stabilitate numerică îmbunătățită , în jurnalul SIAM privind analiza matricei și aplicații , vol. 25, nr. 3, 2003, pp. 669-693.
^ Ammar, Gregory S. și William B. Gragg, Soluție super rapidă a sistemelor Toeplitz definite pozitive reale , în SIAM Journal on Matrix Analysis and Applications , vol. 9, nr. 1, 1988, pp. 61-76.
^ Giacomo Alessandroni, Analiză și modele pentru monitorizarea suprafeței drumurilor ( PDF ), 2016, p. 42, DOI : 10.13140 / RG.2.1.2935.5283 . Adus la 11 octombrie 2019 .

Bibliografie

Surse pentru definiții

Levinson, N. (1947). "Criteriul de eroare Wiener RMS în proiectarea și predicția filtrelor." J. Math. Fizic. , v. 25, pp. 261-278.
Durbin, J. (1960). „Montarea modelelor de serie cronologică”. Rev. Inst. Int. Stat. , v. 28, pp. 233–243.
Trench, WF (1964). „Un algoritm pentru inversiunea matricilor Toeplitz finite.” J. Soc. Indust. Aplic. Matematica. , v. 12, pp. 515-522.
Musicus, BR (1988). „Algoritmi Levinson și Fast Choleski pentru matricele Toeplitz și Almost Toeplitz”. RLE TR nr. 538, MIT. [1]
Delsarte, P. și Genin, YV (1986). „Algoritmul Levinson divizat”. Tranzacții IEEE privind acustica, vorbirea și procesarea semnalului , v. ASSP-34 (3), pp. 470–478.

Lucrări viitoare

AW Bojanczyk, RP Brent, FR De Hoog și DR Sweet, Despre stabilitatea Bareiss și algoritmi de factorizare Toeplitz înrudiți , în SIAM Journal on Matrix Analysis and Applications , vol. 16, 1995, pp. 40–57, DOI : 10.1137 / S0895479891221563 , arXiv : 1004.5510 .
Brent RP (1999), „Stabilitatea algoritmilor rapidi pentru sisteme liniare structurate”, Algoritmi fiabili rapid pentru matrici cu structură (editori - T. Kailath, AH Sayed), cap.4 ( SIAM ).
Bunch, JR (1985). „Stabilitatea metodelor de rezolvare a sistemelor de ecuații Toeplitz”. SIAM J. Sci. Stat. Calculator. , v. 6, pp. 349-364. [2]
H. Krishna și Wang, Y., Split Levinson Algorithm este slab stabil , în SIAM Journal on Numerical Analysis , vol. 30, n. 5, 1993, pp. 1498-1508, DOI : 10.1137 / 0730078 .

Sinteză

Bäckström, T. (2004). "2.2. Levinson - Recurență Durbin." Modelarea predictivă liniară a vorbirii - constrângeri și descompunerea perechii de spectru de linie. Teză de doctorat. Raportul nr. 71 / Universitatea de Tehnologie din Helsinki, Laboratorul de acustică și procesare a semnalului audio. Espoo, Finlanda. [3]
Claerbout, Jon F. (1976). "Capitolul 7 - Aplicații ale formelor de undă ale celor mai mici pătrate." Bazele procesării datelor geofizice. Palo Alto: Blackwell Scientific Publications. [4]
WH Press, SA Teukolsky, WT Vetterling și BP Flannery, secțiunea 2.8.2. Matrici Toeplitz , în Rețete Numerice: Arta Computării Științifice , 3rd, New York, Cambridge University Press, 2007, ISBN 978-0-521-88068-8 .
Golub, GH și Loan, CF Van (1996). "Secțiunea 4.7: Toeplitz și sisteme conexe" Calcule matriciale , Johns Hopkins University Press

Elemente conexe

Portalul de matematică : accesați intrările Wikipedia care se ocupă de matematică

[1] Adam W. Bojanczyk, Brent, RP, De Hoog, FR, & Sweet, DR, Despre stabilitatea Bareiss și algoritmi de factorizare Toeplitz înrudiți , în SIAM Journal on Matrix Analysis and Applications , vol. 16, n. 1, 1995, pp. 40-57.

[2] Brent, Richard P., Stabilitatea algoritmilor rapide pentru sisteme liniare structurate , în Algoritmi fiabili rapid pentru matrici cu structură , Society for Industrial and Applied Mathematics, 1999, pp. 103-116.

[3] Hari Krishna, Yunbiao Wang, Algoritmul Levinson divizat este slab stabil , în jurnalul SIAM privind analiza numerică , vol. 30, n. 5, 1993, pp. 1498-1508.

[4] Bojanczyk, Adam W., Richard P. Brent și Frank R. De Hoog, Un algoritm slab stabil pentru sistemele generale Toeplitz , în Algoritmi numerici , vol. 1, arXiv: 1005.0503, 2010, pp. 225-244.

[stewart2003-5] Stewart, Michael, Un solver Toeplitz super rapid cu stabilitate numerică îmbunătățită , în jurnalul SIAM privind analiza matricei și aplicații , vol. 25, nr. 3, 2003, pp. 669-693.

[6] Ammar, Gregory S. și William B. Gragg, Soluție super rapidă a sistemelor Toeplitz definite pozitive reale , în SIAM Journal on Matrix Analysis and Applications , vol. 9, nr. 1, 1988, pp. 61-76.

[7] Giacomo Alessandroni, Analiză și modele pentru monitorizarea suprafeței drumurilor ( PDF ), 2016, p. 42, DOI : 10.13140 / RG.2.1.2935.5283 . Adus la 11 octombrie 2019 .

[1]

[2] în

[3]

[4]

[5]

[6]

[7]