<jupyter_start><jupyter_text>Estimadores adaptivos parte IIEn esta lección veremos algunos estimadores adaptivos que extienden el filtro LMS que revisamos en la lección anterior  Algoritmo de Mínimos Cuadrados Recursivos El algoritmo LMS minimiza el error instantaneo y es simple y eficiente. Pero en algunos casos su convergencia es demasiado lenta:::{tip}Podemos obtener un filtro adaptivo que converge más rápido si utilizamos el error histórico en lugar del error instantaneo:::Sigamos considerando un filtro tipo FIR con $L+1$ pesos que se actualizan de en cada época $n$$$y_n = \sum_{k=0}^L w_{n, k} u_{n-k}$$El algoritmo RLS (*Recursive Least Squares*) es un método online que minimiza el error histórico, es decir la suma de errores desde la muestra inicial hasta la actual$$\begin{align}J^H_n(\textbf{w}) &= \sum_{i=L}^n   \beta^{n-i} |e_i|^2 \nonumber \\&= \sum_{i=L}^n \beta^{n-i} (d_i - \sum_{k=0}^{L} w_{i, k} u_{i-k} )^2, \nonumber\end{align}$$donde $n$ es el índice del instante actual y $\beta \in [0, 1]$ es el "factor de olvido", que usualmente es un valor cercano pero ligeramente menor que $1$.Adicionalmente se agrega un regularizador a los pesos$$J^w_n = \lambda  \| \textbf{w}_{n} \|^2 = \lambda \sum_{k=1} w_{n, k}^2$$Para evitar divergencias en el proceso de entrenamiento:::{important}La función de costo total del filtro RLS es la suma de error histórico y el regularizador:::  Solución exacta del filtro RLSSi derivamos la función de costo e igualamos a cero obtenemos la siguiente regla$$\begin{align}\textbf{w}_n &= (U_n^T \pmb{\beta} U_n + \lambda I)^{-1}  U_n^T \pmb{\beta} \textbf{d}_n \nonumber \\&= \Phi_n^{-1} \theta_n \nonumber\end{align}$$donde reconocemos los siguientes términos- Matriz de correlación ponderada y regularizada: $\Phi_n = U_n^T \pmb{\beta} U_n + \lambda I$- Vector de correalación cruzada ponderada:  $\theta_n = U_n^T \pmb{\beta} \textbf{d}_n$que se definen en función$$\textbf{d}_n = \begin{pmatrix}  d_n \\ d_{n-1} \\ \vdots \\ d_{L+1} \end{pmatrix} \quad\textbf{u}_n = \begin{pmatrix}  u_n \\ u_{n-1} \\ \vdots \\ u_{n-(L+1)} \end{pmatrix} \quad\pmb{\beta} = I \begin{pmatrix} \beta \\ \beta^{1} \\ \beta^{2}  \vdots \\ \beta^{n-L-1} \end{pmatrix}\quad U_n = \begin{pmatrix}\textbf{u}_n^T \\ \textbf{u}_{n-1}^T \\ \vdots \\ \textbf{u}_{L+1}^T \\\end{pmatrix} \in \mathbb{R}^{n - (L+1) \times L+1}$$donde $I$ es la matriz identidad :::{note}Esta solución es similar a la del filtro de Wiener. Es difícil actualizarla a medida que llegan nuevas observaciones y además es muy costosa debido al cálculo del inverso de la matriz de correlación:::  Solución recursiva del filtro RLSEn lugar de la solución cerrada, es más conveniente actualizar los pesos de forma recursiva. Las condiciones iniciales son - $\Phi_0 = \lambda^{-1} I$- $\theta_0 = 0$y luego la actualización viene dada por - $\Phi_{n} = \beta \Phi_{n-1} + \textbf{u}_n \textbf{u}_n^T$ - $\theta_{n} = \beta \theta_{n-1} + \textbf{u}_n d_n $ - $\textbf{w}_n = \Phi_n^{-1} \theta_n$Podemos evitar invertir la matriz de correlación si usamos el lema de inversión de matrices $$(A + UCV)^{-1} = A^{-1} - A^{-1} U (C^{-1} + VA^{-1} U)^{-1} V A^{-1}$$con $A = \Phi_{n-1}^{-1}$, $C=1$, $U= \textbf{u}_n$ y $V = \textbf{u}_n^T$. De esta forma podemos actualizar $\Phi_{n}^{-1}$ directamente sin tener que invertirla, como se muestra a continuación $$\begin{align}\Phi_{n}^{-1} &= \left(\beta \Phi_{n-1} + \textbf{u}_n \textbf{u}_n^T \right)^{-1} \nonumber \\&= \beta^{-1} \Phi_{n-1}^{-1} - \beta^{-2} \frac{\Phi_{n-1}^{-1} \textbf{u}_n \textbf{u}_n^T \Phi_{n-1}^{-1} }{1 + \beta^{-1} \textbf{u}_n^T \Phi_{n-1}^{-1} \textbf{u}_n} \nonumber \\&= \beta^{-1} \Phi_{n-1}^{-1} - \beta^{-1} \textbf{k}_n \textbf{u}_n^T \Phi_{n-1}^{-1}, \nonumber \end{align}$$donde llamamos **ganancia** al factor$$\textbf{k}_n =  \frac{\beta^{-1} \Phi_{n-1}^{-1} \textbf{u}_n }{1 + \beta^{-1} \textbf{u}_n^T \Phi_{n-1}^{-1} \textbf{u}_n}$$ El último paso es obtener al regla de actualización de pesos$$\begin{align}\textbf{w}_n &= \Phi_n^{-1} \theta_n \nonumber \\&=  \Phi_n^{-1} \left ( \beta \theta_{n-1} + \textbf{u}_n d_n \right) \nonumber \\&=  \left ( \beta^{-1} \Phi_{n-1}^{-1} - \beta^{-1} \textbf{k}_n \textbf{u}_n^T \Phi_{n-1}^{-1} \right ) \beta \theta_{n-1} + \Phi_n^{-1} \textbf{u}_n d_n \nonumber \\&=  \textbf{w}_{n-1} - \textbf{k}_n \textbf{u}_n^T  \textbf{w}_{n-1} + \Phi_n^{-1} \textbf{u}_n d_n \nonumber \\&=  \textbf{w}_{n-1} + \textbf{k}_n ( d_n - \textbf{u}_n^T  \textbf{w}_{n-1} ) \nonumber \\&=  \textbf{w}_{n-1} + \textbf{k}_n e_n \nonumber \end{align}$$donde usamos que $\textbf{w}_{n-1} = \Phi_{n-1}^{-1} \theta_{n-1}$ y $\textbf{k}_n = \Phi_n^{-1} \textbf{u}_n$ :::{note} Con esto tenemos un algoritmo de orden cuadrático en lugar de orden cúbico. Esto sigue siendo mayor que LMS que era de orden lineal pero tiene la ventaja de converger más rapidamente.:::  Resumen del algoritmo RLS e implementación en Python```{prf:algorithm} Algoritmo RLS:nonumber:**Hyper-parámetros:** $L$, $\lambda$, $\beta$1. Inicializar $\Phi_0^{-1} = \lambda I$ y $\textbf{w}_0 = 0$2. Para $n \in [1, \infty]$    - Calcular la ganancia        $$    \textbf{k}_n =  \frac{\Phi_{n-1}^{-1} \textbf{u}_n }{\beta + \textbf{u}_n^T \Phi_{n-1}^{-1} \textbf{u}_n}    $$            - Calcular el error            $$    e_n = d_n - \textbf{u}_n^T  \textbf{w}_{n-1}     $$        - Actualizar el error de pesos            $$    \textbf{w}_n = \textbf{w}_{n-1} + \textbf{k}_n e_n     $$            - Actualizar el inverso de la matriz de correlación            $$    \Phi_{n}^{-1} = \beta^{-1} \Phi_{n-1}^{-1} - \beta^{-1} \textbf{k}_n \textbf{u}_n^T \Phi_{n-1}^{-1}    $$``` - Hiperparámetro $\beta$: Define la memoría efectiva del sistema y repercute en la convergencia y estabilidad del filtro. Como punto de partida se sugiere un valor de $\beta \approx 0.99$. En general $\beta \in [0.9, 1.0)$- Hiperparámetro $\lambda$: Mientras más pequeño sea su valor mayor será la regularización. Se recomienda $\lambda < 0.01/\sigma_u^2$ donde $\sigma_u$ es la desviación estándar de la señal de entrada. En la práctica se pueden calibrar con validación cruzada al igual que $L$<jupyter_code>class Filtro_RLS:
    
    def __init__(self, L, beta=0.99, lamb=1e-2):
        self.L = L
        self.w = np.zeros(shape=(L+1, ))
        self.beta = beta
        self.lamb = lamb
        self.Phi_inv = lamb*np.eye(L+1)
        
    def update(self, un, dn):
        # Cálculo de la ganancia
        pi = np.dot(un.T, self.Phi_inv)
        kn = pi.T/(self.beta + np.inner(pi, un))
        # Actualizar el vector de pesos
        error = dn - np.dot(self.w, un)
        self.w += kn*error
        # Actualizar el inverso de Phi
        self.Phi_inv = (self.Phi_inv - np.outer(kn, pi))*self.beta**-1
        return np.dot(self.w, un)<jupyter_output><empty_output><jupyter_text>Ejemplo: ALE con filtro RLSVeamos como reacciona el filtro RLS ante cambios bruscos usando el ejemplo de la lección pasada. Comparemos con el filtro NLMS<jupyter_code>np.random.seed(12345)
Fs, f0 =  100, 5
t = np.arange(0, 3, 1/Fs)
s = np.sin(2.0*np.pi*t*f0)
s[t>1] += 5
u = s + 0.5*np.random.randn(len(t))

class Filtro_NLMS:
    
    def __init__(self, L, mu, delta=1e-6):
        self.L = L
        self.w = np.zeros(shape=(L+1, ))
        self.mu = mu
        self.delta = delta
        
    def update(self, un, dn):
        unorm = np.dot(un, un) + self.delta
        error = dn - np.dot(self.w, un)
        self.w += 2*self.mu*error*(un/unorm)
        return np.dot(self.w, un)
L = 20
lms = Filtro_NLMS(L, 0.02)
rls = Filtro_RLS(L, 0.99, 1e-2)

u_pred = np.zeros(shape=(len(u), 2))
for k in range(L+1, len(u)):
    u_pred[k, 0] = lms.update(u[k-L-1:k][::-1], u[k])
    u_pred[k, 1] = rls.update(u[k-L-1:k][::-1], u[k])
s1 = hv.Curve((t, s), 'Tiempo', 'Señal', label='Limpia')
s2 = hv.Scatter((t, u), 'Tiempo', 'Señal', label='Contaminada')
s3 = hv.Curve((t, u_pred[:, 0]), 'Tiempo', 'Señal', label='Filtrada (LMS)')
s4 = hv.Curve((t, u_pred[:, 1]), 'Tiempo', 'Señal', label='Filtrada (RLS)')
hv.Overlay([s1, s2, s3, s4]).opts(hv.opts.Overlay(legend_position='top'), 
                                  hv.opts.Curve(ylim=(-5, 10), height=350))<jupyter_output><empty_output><jupyter_text>:::{note} RLS es capaz de seguir los cambios de la señal en menos tiempo que el filtro LMS:::  Perceptrón El perceptrón es un filtro adaptivo desarrollado por [<NAME> en 1962](https://en.wikipedia.org/wiki/Frank_Rosenblatt) con el objetivo de hacer **clasificación supervisada de patrones**Asumiremos que- La respuesta deseada tiene dos categorías: $d_n \in \{-1, +1\}$. El perceptrón resuelve un problema de **clasificación binario**- La entrada es continua y de $L$ dimensiones: $u_n \in \mathbb{R}^L$- Se tienen $N$ tuplas $(u_n, d_n)$ para entrenar el filtroEl filtro tiene arquitectura FIR con $L+1$ coeficientes pero se agrega una función no lineal $\phi(\cdot)$ en la salida del filtro$$\begin{align}y_n &=  \phi \left(b + \sum_{k=1}^{M} w_k u_{nk} \right) \nonumber \\&= \phi \left(b + \langle \textbf{w}, \textbf{u}_n \rangle \right), \nonumber \end{align}$$Los coeficientes del filtro son el escalar $b$ y el vector $\textbf{w}$. :::{note}Este filtro corresponde al modelo matemático de una neurona de [McCulloch y Pitts](https://link.springer.com/article/10.1007/BF02478259), el antecesor de las actuales redes neuronales profundas:::En la implementación original se utilizó la siguiente función no lineal o función de activación$$\phi(z) = \text{sign}(z) = \begin{cases} +1 & z > 0 \\0 & z=0\\-1 & z<0 \end{cases}$$La siguiente figura esquematiza el modelo y su inspiración biológica- Las coeficientes del filtro simulan la importancia o peso de las dendritas- La función no lineal simula el axón que dispara un estímulo eléctrico cuando el voltaje supera un umbral   Ajuste de la neurona artificial: Algoritmo perceptronLa neurona ajusta sus parámetros con cada ejemplo que recibe. Asumiendo que tenemos un dataset con $N$ tuplas $(d_i, u_i)$ donde $d_i$ es la etiqueta y $u_i$ es el vector de entrada```{prf:algorithm} Algoritmo Perceptron:nonumber:**Hyper-parámetros:** $\mu$1. Inicialización de los parámetros: $b=0$, $w_i=0, \forall i = 1, \ldots, L$1. Inicialización del contador: $c=0$2. Para $\text{nepoch}=1,2,\ldots, \infty$    1. Hacer una permutación del dataset    2. $c = 0$    3. Para $n=1,2,\ldots,N$        - Si $\text{sign} \left(b + \langle \textbf{w}, \textbf{u}_n \rangle \right) \neq d_i$ entonces                    $$            \textbf{w} =  \textbf{w} + \mu d_n \textbf{u}_n            $$            y            $$            b = b + \mu d_n            $$                            de lo contrario                    $$            c = c + 1             $$                        - Si $c=N$ entonces detener el entrenamiento        ```Notas:- Se completa una época de entrenamiento cuando se han presentado los $N$ ejemplos del conjunto de entrenamiento- El hiperparámetro $\mu$ es la tasa de aprendizaje de la neurona- Detenemos el entrenamiento cuando todos los ejemplos están bien clasificados. También se puede detener el entrenamiento si se cumple un cierto número fijo de épocas o al cumplir un cierto número de épocas sin cambios en $b$ y $w$- La permutación de los ejemplos en cada época puede evitar sesgos y acelerar la convergencia- El algoritmo perceptrón está garantizado a converger en tiempo finito si el problema es **linealmente separable**. Si el problema no es **linealmente separable** la convergencia se puede forzar disminuyendo gradualmente $\mu$  Interpretación como una aplicación de gradiente descendente estocástico (SGD)El algoritmo de ajuste de la neurona puede considerarse como una minimización de la siguiente función de costo$$\mathcal{L}(b, \textbf{w} ) = \text{max} \Big(0 ~, - d_n ( b + \langle \textbf{w}, \textbf{u}_n \rangle) \Big)$$cuya derivada es $$\frac{d \mathcal{L}}{d \textbf{w}}  = \begin{cases} 0 & d_n ( b + \langle \textbf{w}, \textbf{u}_n \rangle)  \geq 0 \\ - d_n \textbf{u}_n  & d_n ( b + \langle \textbf{w}, \textbf{u}_n \rangle)  < 0\end{cases}$$$$\frac{d \mathcal{L}}{d b}  = \begin{cases} 0 & d_n ( b + \langle \textbf{w}, \textbf{u}_n \rangle)  \geq 0 \\ - d_n   & d_n ( b + \langle \textbf{w}, \textbf{u}_n \rangle)  < 0\end{cases}$$es decir que la derivada es cero si el ejemplo está bien clasificadoNotemos que si aplicamos SGD sobre esta función de costo$$\textbf{w} = \textbf{w} - \mu \frac{d \mathcal{L}}{d \textbf{w}}$$$$b = b - \mu \frac{d \mathcal{L}}{db}$$se recuperan las reglas de ajuste vistas anteriormente  Ejemplo: Clasificación binaria con perceptrónA continuación se muestra como ajustar un perceptrón a medida que se presentan los ejemplos en un problema de clasificación linealmente separable<jupyter_code>N = 5 # Ejemplos por clase
L = 2 # Dimensión de los datos
np.random.seed(1234)
u = np.concatenate((np.random.randn(N, L), 
                    4 + np.random.randn(N, L)))
d = np.ones(shape=(2*N,)); 
d[:N] = -1.
# Parámetros e hiperparámetros
b, w = 0, np.zeros(shape=(L, ))
mu = 1e-6
max_epoch = 2

w_history = np.zeros(shape=(max_epoch*len(d), L))
b_history = np.zeros(shape=(max_epoch*len(d),))
u_history = np.zeros(shape=(max_epoch*len(d), 2))
# Entrenamiento
for nepoch in range(max_epoch):
    idx = np.random.permutation(len(d))
    for n, (un, dn) in enumerate(zip(u[idx], d[idx])):
        if dn*(b+np.inner(w, un)) <= 0.:
            w += mu*dn*un
            b += mu*dn 
            
        u_history[nepoch*len(d)+n, :] = un
        w_history[nepoch*len(d)+n, :] = w
        b_history[nepoch*len(d)+n] = b<jupyter_output><empty_output><jupyter_text>La neurona define un hiperplano que separa el espacio en dos clases. En la siguiente animación se muestra el ajuste de la neurona y en consecuencia la modificación del hiperplano.<jupyter_code>x_plot = np.linspace(-2, 7, num=10)
hiperplano = lambda x, w, b : -b/(w[1]) - x*w[0]/(w[1])

c1 = hv.Points((u[:N, 0], u[:N, 1]), label='Clase 1').opts(size=10, height=350, xlim=(-2, 7), ylim=(-3, 7))
c2 = hv.Points((u[N:, 0], u[N:, 1]), label='Clase 2').opts(size=10)
p = hv.HoloMap(kdims='Iteración')
for i in range(len(b_history)):
    plane = hv.Curve((x_plot, hiperplano(x_plot, w_history[i, :], b_history[i])), label='Hiperplano').opts(color='k')
    selected = hv.Points((u_history[i, 0], u_history[i, 1])).opts(color='k', size=10)
    p[i] = c1 * c2 * plane * selected

hv.output(p.opts(legend_position='top'), holomap='gif', fps=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16749516708091688, 'fraction_numerical': 0.01988400994200497, 'mean_word_length': 3.5110557458735596, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 20, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '268899', 'license': 'MIT', 'n_tokens_mistral': 5705, 'n_tokens_neox': 5549, 'n_words': 1600, 'path': 'lectures/unit3/lecture3.ipynb', 'repo_name': 'phuijse/UACH-INFO183', 'size': '14494'}