Algoritmo di Levinson-Durbin

L'algoritmo di Levinson-Durbin, in algebra lineare, serve per calcolare, con metodo ricorsivo, la soluzione di un'equazione che coinvolge una matrice di Toeplitz. L'algoritmo viene eseguito in $\Theta (n^{2})$ passi, il che rappresenta un forte miglioramento rispetto al Metodo di eliminazione di Gauss, il quale richiede $\Theta (n^{3})$ passaggi.

L'algoritmo Levinson–Durbin fu proposto prima da Norman Levinson, nel 1947, e migliorato da James Durbin nel 1960; successivamente fu ulteriormente migliorato, portandolo da $4n^{2}$ fino a $3n^{2}$ moltiplicazioni, da W. F. Trench e S. Zohar, rispettivamente.

Altri metodi per elaborare i dati, includono la decomposizione di Schur e la decomposizione di Cholesky. In confronto a questi, la ricorsione di Levinson (in particolare la ricorsione di Levinson suddivisa) tende a essere più veloce dal punto di vista computazionale, ma più sensibile a inesattezze computazionali come gli errori di arrotondamento.

L'algoritmo di Bareiss per le matrici di Toeplitz (da non confondere con l'algoritmo di Bareiss generale) è veloce quanto la ricorsione di Levinson-Durbin, ma utilizza $\Theta (n^{2})$ passaggi, mentre la ricorsione di Levinson-Durbin richiede solo $\Theta (n)$ passaggi. L'algoritmo di Bareiss, tuttavia, è numericamente stabile, mentre la ricorsione di Levinson-Durbin, nella migliore delle ipotesi, è solo debolmente stabile (cioè mostra stabilità numerica per sistemi lineari ben condizionati).

Gli algoritmi più recenti, chiamati algoritmi di Toeplitz asintoticamente veloci o, in alcuni testi, superveloci, possono risolvere il problema in $\Theta (n\log _{p}n)$ per vari $p$ (es. $p=2$, $p=3$). La ricorsione di Levinson-Durbin rimane popolare per diversi motivi; primo, è relativamente facile da comprendere; inoltre, può essere più veloce di un algoritmo superveloce per piccoli $n$ (di solito $n<256$).

Derivazione

Introduzione

Le equazioni matriciali sono della seguente forma:

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

L'algoritmo di Levinson-Durbin può essere usato per qualsiasi equazione, purché $\mathbf {M}$ sia una matrice di Toeplitz nota, con diagonale principale diversa da zero; dove ${\vec {y}}$ è il vettore noto, mentre ${\vec {x}}$ è il vettore incognito di numeri $x_{i}$ da determinare.

Si consideri ${\hat {e}}_{i}$un vettore composto interamente da zeri, a eccezione del suo termine i-esimo, il quale contiene un valore unitario. La sua lunghezza sarà implicitamente determinata dal contesto. Il termine $N$ si riferisce alla larghezza della matrice $\mathbf {M}$ avente dimensioni $N\times N$. Infine, gli apici fanno riferimento a un indice induttivo, mentre i pedici indicano gli indici. Per esempio (e definizione) la matrice $\mathbf {T} ^{n}$ è una matrice $n\times n$ che copia il blocco $n\times n$ in alto a sinistra da $\mathbf {M}$, ovvero $T_{i,j}^{n}=M_{j}^{i}$.

Pertanto, anche $\mathbf {T} ^{n}$ è una matrice di Toeplitz; nel senso che può essere scritta nella seguente forma:

$\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}}.$

Passaggi introduttivi

L'algoritmo procede seguendo due passaggi. Nel primo passaggio vengono stabiliti due gruppi di vettori, chiamati vettori avanti e indietro. I vettori avanti sono usati per aiutare a ottenere l'insieme dei vettori indietro; pertanto, possono essere immediatamente scartati. Viceversa, i vettori all'indietro sono necessari per il secondo passaggio, dove vengono utilizzati per creare la soluzione desiderata.

La ricorsione di Levinson-Durbin definisce l'ennesimo "vettore avanti", chiamato ${\vec {f}}^{n}$, un vettore di lunghezza n che soddisfa l'equazione:

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

L'ennesimo "vettore indietro", chiamato ${\vec {b}}^{n}$è definito in modo analogo; è il vettore di lunghezza n che soddisfa l'equazione:

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

Una semplificazione importante può verificarsi quando $\mathbf {M}$ è una matrice simmetrica; in questo caso i due vettori sono correlati da $b_{i}^{n}=f_{n+1-i}^{n}$, ovvero sono le inversioni di riga l'uno dell'altra. Questo può risparmiare alcuni calcoli in questo caso particolare.

Ottenere i vettori indietro

Anche se la matrice non è simmetrica, l'ennesimo vettore avanti e indietro può essere trovato dai vettori di lunghezza $n-1$ come segue. Innanzitutto, il vettore indietro può essere "esteso" con uno zero al fine di ottenere:

$\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}}.$

Andando da $\mathbf {T} ^{n-1}$ a $\mathbf {T} ^{n}$, la colonna aggiuntiva aggiunta alla matrice non disturba la soluzione quando si usa uno zero per estendere il vettore in avanti. Tuttavia, la riga aggiuntiva aggiunta alla matrice ha perturbato la soluzione; e ha creato un termine di errore indesiderato ε_f che si verifica in ultimo luogo. L'equazione precedente fornisce il valore di:

$\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}.$

Questo errore verrà restituito a breve ed eliminato dal nuovo vettore avanti; ma prima, il vettore all'indietro deve essere esteso in modo simile (anche se invertito). Per il vettore all'indietro si ha:

$\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}}.$

Come prima, la colonna aggiuntiva alla matrice non disturba questo nuovo vettore all'indietro; ma la riga aggiuntiva lo fa. Qui si ha un altro errore indesiderato $\epsilon$ pari a:

$\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}.\$

Questi due termini di errore possono essere utilizzati per formare vettori avanti e indietro di ordine superiore descritti come segue. Usando la linearità delle matrici, la seguente identità vale per tutti gli $(\alpha ,\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}}.$

Se $\alpha$ e $\beta$ sono scelti in modo destrorso ${\hat {e}}_{1}$ o ${\hat {e}}_{n,n}$_, allora la quantità in parentesi sarà pari alla definizione dell'ennesimo vettori avanti o indietro, rispettivamente. Con la scelta dei termini $\alpha$ e $\beta$, la somma dei vettori tra parentesi è semplice e produce il risultato desiderato.

Per trovare questi coefficienti, $\alpha _{f}^{n}$, $\beta _{f}^{n}$ devono essere tali che:

${\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}}$

e, rispettivamente, $\alpha _{b}^{n}$, $\beta _{b}^{n}$ sono tali che:

${\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}}.$

Moltiplicando e dividendo le precedenti equazioni per ${\mathbf {T} }^{n}$si ottiene la seguente equazione:

${\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}}.$

Ora, tutti gli zeri dei due vettori sopra vengono ignorati, rimane solo la seguente equazione:

${\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}}.$

Con questi soluzione (utilizzando la formula inversa della matrice di Cramer $2\times 2$), i nuovi vettori avanti e indietro sono:

${\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}}.$

L'esecuzione di queste sommatorie vettoriali, quindi, dà gli ennesimi vettori avanti e indietro partendo da quelli precedenti. Non resta che trovare il primo di questi vettori, quindi con rapide somme e moltiplicazioni rapide si ottengono i termini restanti. I primi vettori avanti e indietro sono semplicemente:

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

Utilizzo dei vettori all'indietro

I passaggi precedenti danno agli N vettori all'indietro per $\mathbf {M}$. Da lì, un'equazione più arbitraria è:

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

La soluzione può essere costruita nello stesso modo ricorsivo in cui sono stati costruiti i vettori all'indietro. Di conseguenza, ${\vec {x}}$ deve essere generalizzato a una sequenza di ${\vec {x}}^{n}$intermedi, tali che ${\vec {x}}^{N}={\vec {x}}$.

La soluzione viene quindi costruita ricorsivamente notando che, se

$\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}},$

allora, inserendo nuovamente uno zero al termine del vettore, e definendo una costante di errore dove necessario, si ha:

$\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}}.$

Possiamo quindi utilizzare l'ennesimo vettore all'indietro per eliminare l'errore e sostituirlo con la formula desiderata come segue:

$\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}}.$

L'estensione di questo metodo fino a quando $n=N$ produce la soluzione desiderata ${\vec {x}}$ .

In pratica, questi passaggi vengono spesso eseguiti in concomitanza con il resto della procedura, formando un'unità coerente e meritano di essere trattati singolarmente.

Algoritmo di Levinson-Durbin

Se $\mathbf {M}$ non è una matrice di Toeplitz in senso stretto, ma una matrice a blocchi di Toeplitz, la ricorsione di Levinson può essere derivata più o meno allo stesso modo considerando la sottomatrice di Toeplitz (Musicus 1988).

Le matrici con blocchi di Toeplitz sorgono naturalmente negli algoritmi di elaborazione del segnale quando si trattano flussi di segnali multipli (per esempio nei sistemi MIMO) o segnali ciclostazionari.

Applicazione pratica dell'algoritmo di Levinson-Durbin

L'algoritmo di Levinson-Durbin è molto utilizzato per la risoluzione dei modelli autoregressivi $AR(p)$ di ordine $p$ (utilizzati nel protocollo GSM), i quali si presentano sotto la seguente equazione alle differenze:

${\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}},$

dove $t(n)$ è il campione attuale del sistema stimato tramite il modello $AR(p)$, $x_{i}$ sono i $p$ parametri del modello, i quali vengono applicati ai suoi $p$ campioni precedenti, o - analogamente - definiscono la memoria del modello (in questo caso di ordine $p$). Infine $e(n)$ è il residuo di predizione del campione $t(n)$ che si desidera minimizzare, ovvero: la componente non predicibile del sistema.

Queste equazioni si presentano sotto forma di matrici di Toeplitz, pertanto per la loro soluzione si ricorre all'algoritmo di Levinson-Durbin.

Pseudocodice per la ricorsione Levinson-Durbin

L'algoritmo si basa sul calcolo dei coefficienti autoregressivi di matrici di ordine crescente. Si divide in due fasi: una prima di inizializzazione per il calcolo del parametro $x_{1}$ o - che è lo stesso - per il caso elementare della matrice $1\times 1$. Successivamente, si procede con il calcolo iterativo dei parametri per matrici di ordine via, via crescente $2\times 2$, $3\times 3$,..., $p\times p$.

Utilizzando la notazione MATLAB/Octave, il pseudocodice per il calcolo della ricorsione di Levinson-Durbin è il seguente:

k = M(2) / M(1);                                    % Stima il primo elemento 
x = k;
E = (1 - k2) * M(1);                                % Calcola l'errore quadratico medio
for i = 2:p
    k = (M(i + 1) - x * M(2:i)) / E;                % Coefficienti di riflessione
    x = [k, x - k * x(i-1:-1:1)];                   % Stima gli elementi successivi
    E = (1 - k2) * E;                               % Aggiorna l'errore quadratico medio
end

x = [1, - x(N:-1:1)];                               % Restituisce vettore incognito

Da sottolineare che MATLAB opera con vettori e matrici, pertanto - volendo tradurre il codice in linguaggi come C++ o Java, si otterranno due cicli for annidati.

Bibliografia

Fonti per le definizioni

• Levinson, N. (1947). "The Wiener RMS error criterion in filter design and prediction." J. Math. Phys., v. 25, pp. 261–278.
• Durbin, J. (1960). "The fitting of time series models." Rev. Inst. Int. Stat., v. 28, pp. 233–243.
• Trench, W. F. (1964). "An algorithm for the inversion of finite Toeplitz matrices." J. Soc. Indust. Appl. Math., v. 12, pp. 515–522.
• Musicus, B. R. (1988). "Levinson and Fast Choleski Algorithms for Toeplitz and Almost Toeplitz Matrices." RLE TR No. 538, MIT. [1]
• Delsarte, P. and Genin, Y. V. (1986). "The split Levinson algorithm." IEEE Transactions on Acoustics, Speech, and Signal Processing, v. ASSP-34(3), pp. 470–478.

Lavori futuri

• A.W. Bojanczyk, R.P. Brent, F.R. De Hoog e D.R. Sweet, On the stability of the Bareiss and related Toeplitz factorization algorithms, in SIAM Journal on Matrix Analysis and Applications, vol. 16, 1995, pp. 40–57, DOI:10.1137/S0895479891221563, arXiv:1004.5510.
• Brent R.P. (1999), "Stability of fast algorithms for structured linear systems", Fast Reliable Algorithms for Matrices with Structure (editors—T. Kailath, A.H. Sayed), ch.4 (SIAM).
• Bunch, J. R. (1985). "Stability of methods for solving Toeplitz systems of equations." SIAM J. Sci. Stat. Comput., v. 6, pp. 349–364. [2]
• H. Krishna e Wang, Y., The Split Levinson Algorithm is weakly stable, in SIAM Journal on Numerical Analysis, vol. 30, n. 5, 1993, pp. 1498–1508, DOI:10.1137/0730078.

Sintesi

• Bäckström, T. (2004). "2.2. Levinson–Durbin Recursion." Linear Predictive Modelling of Speech – Constraints and Line Spectrum Pair Decomposition. Doctoral thesis. Report no. 71 / Helsinki University of Technology, Laboratory of Acoustics and Audio Signal Processing. Espoo, Finland. [3]
• Claerbout, Jon F. (1976). "Chapter 7 – Waveform Applications of Least-Squares." Fundamentals of Geophysical Data Processing. Palo Alto: Blackwell Scientific Publications. [4]
• WH Press, SA Teukolsky, WT Vetterling e BP Flannery, Section 2.8.2. Toeplitz Matrices, in Numerical Recipes: The Art of Scientific Computing, 3rd, New York, Cambridge University Press, 2007, ISBN 978-0-521-88068-8.
• Golub, G.H., and Loan, C.F. Van (1996). "Section 4.7 : Toeplitz and related Systems" Matrix Computations, Johns Hopkins University Press

Voci correlate

• Equazioni di Yule-Walker
• Modello autoregressivo