El algoritmo MSEG

Las ecuaciones de Wiener-Hopf dan el mejor ecualizador de un saque, $c_{\mathrm{opt}}={\mathrm{\Phi }}^{-1}\alpha$, pero piden invertir $\mathrm{\Phi }$. El algoritmo MSEG (minimum square error gradient) llega al mismo punto de otra manera: arranca en cualquier lado y baja por el tazón del MSE dando pasos contra el gradiente. Si llega y cuánto tarda lo decide la forma del tazón, o sea los autovalores de $\mathrm{\Phi }$. Cada uno aporta un modo que se achica como ${(1-\beta {\lambda }_i)}^j$: el autovalor más grande pone el techo del paso, y el más chico, la velocidad. En la cadena canal → AGC → ecualizador → slicer, es la regla con la que el ecualizador movería sus taps si conociera la estadística de su entrada. El LMS es la misma regla sin conocerla.

01 · Dónde se ubica

La filmina de introducción a la optimización ordena los métodos por dos preguntas: cuántas variables tiene la función y si se dispone de sus derivadas. Con una sola variable están la sección áurea, la interpolación cuadrática y el método de Newton. Con varias hay métodos directos, que sólo evalúan la función —búsqueda por malla, búsqueda aleatoria, univariadas, Powell—, y métodos con gradiente, que usan las derivadas para elegir hacia dónde moverse: el steepest ascent para máximos y el steepest descent para mínimos.

El gradiente apunta hacia donde la función crece más rápido y es perpendicular a las curvas de nivel. El steepest descent da un paso en la dirección contraria, vuelve a evaluar el gradiente y repite. La filmina distingue dos maneras de elegir el largo del paso: uno arbitrario y fijo, o el $h^{\ast }$ que minimiza la función a lo largo de esa dirección, que es una búsqueda unidimensional en cada paso. El MSEG es el primero: descenso por gradiente con paso fijo sobre el MSE de un ecualizador de $N$ taps. Un método directo no es opción: una malla de 10 valores por tap sobre 11 taps son 10¹¹ evaluaciones.

02 · El gradiente del MSE

El MSE del ecualizador es una forma cuadrática en los taps. De dónde sale, y por qué tiene un solo mínimo, está en el criterio MSE:

$$E[{|E_k|}^2]=E[{|A_k|}^2]-2\mathrm{\Re }(c^H\alpha )+c^H\mathrm{\Phi }c$$

con $\mathrm{\Phi }=E[r_k^{\ast }{r}_k^T]$ y $\alpha =E[A_k{r}_k^{\ast }]$. Con señales reales —PAM-2, o la entrada de MSEG.m— los asteriscos y la parte real no hacen nada, y queda $E[A_k^2]-2c^T\alpha +c^T\mathrm{\Phi }c$. El gradiente es el vector de las derivadas parciales respecto de cada tap $c_n$, y sale término por término.

Paso 1: el término constante. $E[A_k^2]$ no depende de los taps: su gradiente es cero.

Paso 2: el término lineal. $c^T\alpha ={\sum }_n{c}_n{\alpha }_n$. Cada tap aparece una sola vez, multiplicado por su ${\alpha }_n$, así que la derivada respecto de $c_n$ es ${\alpha }_n$, y el gradiente es $\alpha$ entero.

Paso 3: el término cuadrático. Escrito como doble suma, $c$ aparece dos veces, una como índice de fila y otra como índice de columna:

$$c^T\mathrm{\Phi }c=\sum _m\sum _n{c}_m{[\mathrm{\Phi }]}_{mn}{c}_n\Rightarrow \frac{\partial }{\partial c_n}{c}^T\mathrm{\Phi }c=\sum _m{[\mathrm{\Phi }]}_{nm}{c}_m+\sum _m{c}_m{[\mathrm{\Phi }]}_{mn}=2{(\mathrm{\Phi }c)}_n$$

Las dos sumas son iguales porque $\mathrm{\Phi }$ es simétrica: la correlación entre los taps $m$ y $n$ es la misma que entre $n$ y $m$. Juntando los tres pasos:

$${\nabla }_{c}E[{|E_k|}^2]=2\mathrm{\Phi }c-2\alpha$$

Es el gradiente que escribe la filmina 4, que vale igual con señales complejas. Se puede leer de dos maneras. Igualado a cero da las ecuaciones de Wiener-Hopf, $\mathrm{\Phi }c=\alpha$: es el camino B de la filmina hacia $c_{\mathrm{opt}}$. Y como $\alpha -\mathrm{\Phi }c=E[E_k{r}_k^{\ast }]$, el gradiente es $-2$ veces la correlación entre el error y las muestras de la línea. Esa correlación se puede medir sobre las señales, y de ahí sale el LMS.

Hay una tercera forma del MSE que conviene tener a mano. Con $q=c-c_{\mathrm{opt}}$, completar el cuadrado da $E[{|E_k|}^2]={\xi }_{\text{min}}+q^H\mathrm{\Phi }q$. Las curvas de nivel son elipses centradas en $c_{\mathrm{opt}}$, con los ejes sobre los autovectores de $\mathrm{\Phi }$: empinadas en la dirección del autovalor más grande, casi planas en la del más chico.

03 · La actualización

Partiendo de un $c_0$ cualquiera, el MSEG da un paso contra el gradiente evaluado donde está:

$$c_{j+1}=c_j-\frac{\beta }{2}{\nabla }_{c_j}E[{|E_k|}^2]=c_j-\frac{\beta }{2}(2\mathrm{\Phi }{c}_j-2\alpha )=(I-\beta \mathrm{\Phi })c_j+\beta \alpha$$

$\beta$ es el paso de adaptación, y el $\beta /2$ está para que se vaya el 2 del gradiente. El índice $j$ cuenta iteraciones del algoritmo, no símbolos: el tiempo del canal es $k$, y el MSEG no lo mira, porque trabaja con la esperanza. Como $\mathrm{\Phi }$ y $\alpha$ se suponen conocidas, el gradiente es exacto y el algoritmo es determinístico: la misma $c_0$ da siempre la misma trayectoria. La versión que estima el gradiente con cada símbolo es el LMS.

Cada iteración cuesta un producto $\mathrm{\Phi }{c}_j$ y ninguna inversión. El método de Newton de la filmina de optimización hace lo contrario: multiplica el gradiente por la inversa del Hessiano, que acá es $2\mathrm{\Phi }$, y cae en el fondo en un paso, $c-{(2\mathrm{\Phi })}^{-1}(2\mathrm{\Phi }c-2\alpha )={\mathrm{\Phi }}^{-1}\alpha$. Newton sobre una cuadrática es Wiener-Hopf con otro nombre. El MSEG cambia esa inversa por un número, $\beta$, y lo paga en iteraciones. Cuántas, depende de cuánto se parezca $\mathrm{\Phi }$ a un múltiplo de la identidad.

04 · El error de coeficientes y sus modos

Paso 1: medir desde el fondo. Se define el error de coeficientes $q_j=c_j-c_{\mathrm{opt}}$. Como $\alpha =\mathrm{\Phi }{c}_{\mathrm{opt}}$, el término independiente de la actualización es $\beta \mathrm{\Phi }{c}_{\mathrm{opt}}$, y al restar $c_{\mathrm{opt}}$ de los dos lados desaparece:

$$c_{j+1}-c_{\mathrm{opt}}=(I-\beta \mathrm{\Phi })c_j+\beta \mathrm{\Phi }{c}_{\mathrm{opt}}-c_{\mathrm{opt}}\Rightarrow q_{j+1}=(I-\beta \mathrm{\Phi })q_j={(I-\beta \mathrm{\Phi })}^{j+1}{q}_0$$

Una recursión lineal homogénea: todo depende de las potencias de $I-\beta \mathrm{\Phi }$.

Paso 2: diagonalizar. $\mathrm{\Phi }$ es hermítica, así que sus autovalores son reales y sus autovectores $v_i$ ortonormales; es definida positiva, así que los ${\lambda }_i$ son positivos. Es la descomposición espectral de la filmina, $\mathrm{\Phi }=V\mathrm{\Lambda }{V}^H={\sum }_i{\lambda }_i{v}_i{v}_i^H$. $I-\beta \mathrm{\Phi }$ tiene los mismos autovectores, con autovalores $1-\beta {\lambda }_i$, y elevarla a la $j$ es elevar cada uno:

$${(I-\beta \mathrm{\Phi })}^j=\sum _{i=1}^N{(1-\beta {\lambda }_i)}^j{v}_i{v}_i^H$$

Paso 3: mirar el error en esa base. La coordenada $u_{i,j}=v_i^H{q}_j$ dice cuánto error queda a lo largo del autovector $v_i$, y con la suma de arriba evoluciona sola:

$$u_{i,j}={(1-\beta {\lambda }_i)}^j{u}_{i,0}$$

Son los modos de convergencia de la filmina. Un sistema de $N$ ecuaciones acopladas se vuelve $N$ geométricas sueltas, cada una con su razón $1-\beta {\lambda }_i$, y la trayectoria de $q_j$ es la suma de las $N$. El MSE también se separa, porque $q^H\mathrm{\Phi }q={\sum }_i{\lambda }_i{|u_i|}^2$:

$$E[{|E_k|}^2]{|}_{c_j}={\xi }_{\text{min}}+\sum _{i=1}^N{\lambda }_i{(1-\beta {\lambda }_i)}^{2j}{|u_{i,0}|}^2$$

En dB, el exceso de MSE de cada modo baja $20{\log }_{10}|1-\beta {\lambda }_i|$ por iteración: una recta, que es lo que dibuja MSEG.m en su figura «Avance del MSE» (sección 08).

Con $\chi =1$, $\mathrm{\Phi }=I$ y el gradiente, $2(c-c_{\mathrm{opt}})$, apunta derecho al fondo desde cualquier lado. Los dos modos tienen la misma razón, las dos rectas de abajo son paralelas y la trayectoria es una recta; con $\beta =1$ llega en un paso, que es el paso de Newton. Al subir $\chi$ las elipses se estiran sobre la dirección de ${\lambda }_{\text{min}}$ y el gradiente deja de apuntar al fondo: apunta a la pared más cercana. El modo de ${\lambda }_{\text{max}}$ se va en pocos pasos y lo que queda es recorrer el valle con la razón de ${\lambda }_{\text{min}}$, por eso la trayectoria dobla: primero cae a lo ancho del valle y después camina a lo largo. Con $\beta$ entre $1/{\lambda }_{\text{max}}$ y $2/{\lambda }_{\text{max}}$ el modo rápido cambia de signo en cada paso y la trayectoria zigzaguea de pared a pared, como en la figura de la filmina. Arriba de $2/{\lambda }_{\text{max}}$ cada rebote es más grande que el anterior. En el widget ${\phi }_0=1$, así que ${\lambda }_{\text{min}}+{\lambda }_{\text{max}}=2$ y el paso óptimo de la sección 06 vale 1 para cualquier $\chi$.

05 · Cuándo converge

$q_j$ tiende a cero para cualquier punto de partida si y sólo si todos los modos se achican: $|1-\beta {\lambda }_i|<1$ para todo $i$. Como los ${\lambda }_i$ son positivos, eso es $0<\beta {\lambda }_i<2$, y el más exigente es el del autovalor más grande:

$$0<\beta <\frac{2}{{\lambda }_{\text{max}}}$$

La dirección más empinada del tazón es la primera que se pasa de largo. Si el paso es tan grande que en esa dirección el algoritmo salta más allá de la pared de enfrente, cada rebote es mayor que el anterior, aunque en todas las otras direcciones esté convergiendo: alcanza con un modo que crezca. Entre $1/{\lambda }_{\text{max}}$ y $2/{\lambda }_{\text{max}}$ converge oscilando; justo en $2/{\lambda }_{\text{max}}$ rebota sin achicarse.

La filmina escribe primero la condición como $|1-\beta {\phi }_0|<1$. Ese es el caso de entrada blanca, donde $\mathrm{\Phi }={\phi }_{0}I$ y todos los autovalores valen ${\phi }_0$; con una entrada cualquiera la condición es sobre ${\lambda }_{\text{max}}$, que es la que la filmina da tres filminas después.

Hay una cota que no pide autovalores. Los autovalores suman la traza, y como son todos positivos, ${\lambda }_{\text{max}}\le \mathrm{tr}\mathrm{\Phi }=N{\phi }_0$. Así, $\beta <2/(N{\phi }_0)$ es seguro y sólo pide el largo del ecualizador y la potencia de entrada ${\phi }_0$. Si la señal llega el doble de fuerte, $\mathrm{\Phi }$ se duplica y el paso máximo se parte a la mitad: por eso el control automático de ganancia va antes del ecualizador.

En MSEG.m la entrada es randn: blanca y de potencia 1, así que $\mathrm{\Phi }=I$ y hay un solo modo, con $\lambda =1$. El script actualiza c = c - beta*gradient, sin el $1/2$ de la filmina: su beta es la mitad del $\beta$ de la teoría. El beta = 0.001 de la clase es $\beta =\mathrm{0,002}$, con razón $\mathrm{0,998}$ por iteración, y el límite $\beta <2$ es beta < 1 en el script. La filmina pide probar varios beta. Con el lazo de MSEG.m, sin tocar, lo que se achica el error en la primera iteración:

El único que se aparta es $\beta =1$: la teoría dice que llega en un paso, y el script se queda a un 3,9 % del error inicial. Es que el script no tiene el gradiente exacto. Lo calcula sobre un lote de 1000 muestras, y la $\mathrm{\Phi }$ de un lote no es exactamente $I$: en 200 lotes, su autovalor más chico promedia 0,916 y el más grande 1,075. Con $\beta =1$ lo único que queda es esa diferencia.

06 · Qué tan rápido: el modo lento

Dentro del rango estable, la velocidad la pone el modo que se achica más despacio:

$$\rho (\beta )=\underset{i}{\max }|1-\beta {\lambda }_i|,\parallel q_j\parallel \le {\rho }^j\parallel q_0\parallel$$

Cada $|1-\beta {\lambda }_i|$, como función de $\beta$, es una V con el vértice en $1/{\lambda }_i$. Las V de los autovalores intermedios quedan entre las de los extremos, así que el máximo lo ponen ${\lambda }_{\text{min}}$ y ${\lambda }_{\text{max}}$, que es el gráfico de la filmina. Con $\beta$ chico manda ${\lambda }_{\text{min}}$, la V que baja más despacio; con $\beta$ grande manda ${\lambda }_{\text{max}}$, la que primero vuelve a subir. El mínimo del máximo está donde se cruzan la rama que baja de una y la que sube de la otra:

$$1-\beta {\lambda }_{\text{min}}=\beta {\lambda }_{\text{max}}-1\Rightarrow {\beta }_{\mathrm{opt}}=\frac{2}{{\lambda }_{\text{min}}+{\lambda }_{\text{max}}}$$

Reemplazando, la mejor razón posible depende de un solo número, la dispersión de autovalores $\chi ={\lambda }_{\text{max}}/{\lambda }_{\text{min}}$:

$$\rho ({\beta }_{\mathrm{opt}})=\frac{{\lambda }_{\text{max}}-{\lambda }_{\text{min}}}{{\lambda }_{\text{max}}+{\lambda }_{\text{min}}}=\frac{\chi -1}{\chi +1}$$

En ${\beta }_{\mathrm{opt}}$ los dos modos extremos empatan: el lento se acerca de un solo lado y el rápido oscilando, con la misma razón en módulo. Achicar el paso ayuda al rápido y frena al lento; agrandarlo, al revés. Ninguna elección de $\beta$ hace rápido un tazón alargado. Ni siquiera en el borde de estabilidad: si $\beta <2/{\lambda }_{\text{max}}$, entonces $\beta {\lambda }_{\text{min}}<2/\chi$, y el modo lento nunca baja más de $2/\chi$ por iteración. Cuánto tarda la cota ${\rho }^j$ en bajar a 10⁻³ con el mejor paso:

Es $\ln {10}^{-3}/\ln \rho$, en forma cerrada. Para $\chi$ grande, $\ln \rho \approx -2/\chi$ y hacen falta unas $\mathrm{3,45}\chi$ iteraciones: el costo crece lineal con la dispersión. Las filas de 13,41, 49,75 y 126,8 son el canal de la sección siguiente con 5, 11 y 21 taps.

07 · La dispersión la pone el espectro

$\mathrm{\Phi }$ es la autocorrelación de lo que entra al ecualizador, o sea de lo que dejó el canal, y sus autovalores quedan encerrados por el espectro de esa señal. Para un autovector $v$ de norma 1, $\lambda =v^H\mathrm{\Phi }v={\sum }_m{\sum }_n{v}_m^{\ast }{\phi }_{m-n}{v}_n$. Escribiendo cada ${\phi }_{m-n}$ como la antitransformada de la PSD $S_R$ de $R_k$, la doble suma se vuelve el módulo al cuadrado de la DTFT de los componentes de $v$:

$$\lambda =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{S}_R(e^{j\omega }){|V(e^{j\omega })|}^{2}d\omega ,\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{|V(e^{j\omega })|}^{2}d\omega =1$$

La segunda igualdad es Parseval. Cada autovalor es un promedio de $S_R$ con pesos que suman 1, así que cae entre el mínimo y el máximo del espectro, y $\chi \le \max S_R/\min S_R$. Es lo que dice la filmina: un espectro casi plano converge rápido, uno con picos y pozos converge lento.

El caso más simple: símbolos blancos de potencia 1 por un canal [1 a], $R_k=A_k+a{A}_{k-1}$. Entonces ${\phi }_0=1+a^2$, ${\phi }_1=a$, y $S_R=1+a^2+2a\cos \omega$ va de ${(1-a)}^2$ a ${(1+a)}^2$. $\mathrm{\Phi }$ es tridiagonal y Toeplitz, y sus autovalores tienen forma cerrada:

$${\lambda }_i=1+a^2+2a\cos \frac{i\pi }{N+1},i=1,\dots ,N$$

que es muestrear $S_R$ en $N$ frecuencias. Contra eig difieren en 2·10⁻¹⁵ como máximo. La dispersión, según el eco y el largo del ecualizador:

Con $a=\mathrm{0,9}$ el canal tiene un pozo de 25,6 dB en $\omega =\pi$, y la cota es 361. Con más taps el ecualizador ve el espectro con más resolución y $\chi$ se acerca a la cota: 13,4 con 5 taps, 126,8 con 21. Con el mejor paso, el de 21 taps tarda 438 iteraciones en bajar la cota del error mil veces; con el canal plano, una. El canal que más necesita ecualización es el que más tarda en aprenderse, y no hay paso que lo arregle, porque es una propiedad de la entrada y no del algoritmo. Con dos taps, este canal no pasa de $\chi =3$: las dispersiones más grandes del widget piden muestras vecinas más parecidas entre sí, como las que deja un canal más largo.

08 · El análisis en MATLAB

El script de la clase, MSEG.m, identifica un filtro de referencia h_ref = fir1(4, 0.3)' de 5 taps: la misma entrada blanca pasa por h_ref y por el filtro c que se adapta, y el error es la diferencia de las salidas. No usa $\mathrm{\Phi }$ ni $\alpha$: en cada iteración mide el MSE sobre un lote de 1000 muestras, suma Delta_tap = 1e-4 a un tap por vez y estima cada derivada parcial como un cociente incremental. Corre sin cambios en R2026a y tarda unos 5 s. Abajo, el script corrido con semilla fija, y después su mismo lazo con otros pasos y con la entrada coloreada. MSEG.m tiene que estar en la carpeta. Lo que figura como comentario es la salida real.

%% 1. El script de la clase, tal cual (tarda unos 5 s)
rng(1)
run('MSEG.m')        % deja h_ref, c, beta, Niters, Ncoeff, Delta_tap, mse_logeo
close all

%% 2. La predicción. Entrada blanca de potencia 1: Phi = I, un solo modo con lambda = 1.
% El script hace c = c - beta*gradient, sin el 1/2 de la filmina:
% su beta es la mitad del beta de la teoría, así que el factor es 1 - 2*beta.
factor = 1 - 2*beta;
j = (1:Niters)';
mse_dB = 10*log10(mse_logeo);
pred_dB = 10*log10(factor.^(2*(j-1)) * sum(h_ref.^2));   % MSE = ||q_j||^2
p = polyfit(j(200:3000), mse_dB(200:3000), 1);
fprintf('pendiente: medida %.2f, predicha %.2f dB cada 1000 iteraciones\n', ...
        1000*p(1), 1000*20*log10(factor))
fprintf('MSE en j = 1000: medido %.2f dB, predicho %.2f dB\n', mse_dB(1000), pred_dB(1000))
% pendiente: medida -17.28, predicha -17.39 dB cada 1000 iteraciones
% MSE en j = 1000: medido -21.63 dB, predicho -21.86 dB

%% 3. El piso no es el mínimo: la derivada hacia adelante corre el gradiente
fprintf('c - h_ref = %s\n', mat2str((c - h_ref)', 3))
fprintf('piso: medido %.2f dB, N*(Delta/2)^2 = %.2f dB\n', ...
        mean(mse_dB(8001:end)), 10*log10(Ncoeff*(Delta_tap/2)^2))
% c - h_ref = [-5e-05 -5.01e-05 -5.01e-05 -5.01e-05 -5e-05]
% piso: medido -79.05 dB, N*(Delta/2)^2 = -79.03 dB

%% 4. Probar varios betas (lo que pide la filmina), con el lazo de MSEG.m
h = fir1(4, 0.3)';
for b = [0.01 0.25 0.5 0.75 0.95 1.05]          % beta del script
    rng(2)
    cl = mseg_clase(h, b, 40, 1);               % entrada blanca
    nq = vecnorm(cl - h', 2, 2);                % ||q_j||, j = 0..40
    fprintf('beta script %.2f: ||q1||/||q0|| = %.3f, predicho |1-2beta| = %.2f\n', ...
            b, nq(2)/nq(1), abs(1 - 2*b))
end
% beta script 0.01: ||q1||/||q0|| = 0.980, predicho |1-2beta| = 0.98
% beta script 0.25: ||q1||/||q0|| = 0.498, predicho |1-2beta| = 0.50
% beta script 0.50: ||q1||/||q0|| = 0.039, predicho |1-2beta| = 0.00
% beta script 0.75: ||q1||/||q0|| = 0.511, predicho |1-2beta| = 0.50
% beta script 0.95: ||q1||/||q0|| = 0.913, predicho |1-2beta| = 0.90
% beta script 1.05: ||q1||/||q0|| = 1.114, predicho |1-2beta| = 1.10

%% 5. Entrada coloreada: x pasa por [1 0.9] antes de entrar
a = 0.9;  N = 5;
Phi = toeplitz([1+a^2, a, zeros(1, N-2)]);      % autocorrelación exacta de x
[V, L] = eig(Phi);  lam = diag(L);
fprintf('lambda = %s, dispersion %.2f\n', mat2str(lam', 4), max(lam)/min(lam))
% lambda = [0.2512 0.91 1.81 2.71 3.369], dispersion 13.41
beta_opt = 2/(min(lam) + max(lam));
for B = [0.2, beta_opt]                         % beta de la filmina
    rng(3)
    cl = mseg_clase(h, B/2, 120, [1 a]);        % el script usa beta/2
    nq = vecnorm(cl - h', 2, 2);
    u = (cl - h') * V;                          % coordenadas modales de q_j
    pl = polyfit((0:30)', log(abs(u(1:31, 1))), 1);      % modo de lambda_min
    fprintf('beta %.4f: modo lento %.4f (pred %.4f), ||q|| baja 100 veces en %d iteraciones\n', ...
            B, exp(pl(1)), abs(1 - B*min(lam)), find(nq < nq(1)/100, 1) - 1)
end
% beta 0.2000: modo lento 0.9496 (pred 0.9498), ||q|| baja 100 veces en 48 iteraciones
% beta 0.5525: modo lento 0.8677 (pred 0.8612), ||q|| baja 100 veces en 30 iteraciones
B = 1.1*2/max(lam);                             % 10 % arriba del limite
rng(3)
cl = mseg_clase(h, B/2, 60, [1 a]);
u = (cl - h') * V;
pr = polyfit((20:60)', log(abs(u(21:61, end))), 1);      % modo de lambda_max
fprintf('beta %.4f > 2/lambda_max = %.4f: el modo rapido crece %.3f por iteracion (pred %.3f)\n', ...
        B, 2/max(lam), exp(pr(1)), abs(1 - B*max(lam)))
% beta 0.6530 > 2/lambda_max = 0.5937: el modo rapido crece 1.162 por iteracion (pred 1.200)
rng(3)
cl = mseg_clase(h, 0.5, 20, 1);                 % entrada blanca, beta = 1 = su beta_opt
nq = vecnorm(cl - h', 2, 2);
fprintf('entrada blanca, beta = 1: ||q|| baja 100 veces en %d iteraciones\n', find(nq < nq(1)/100, 1) - 1)
% entrada blanca, beta = 1: ||q|| baja 100 veces en 2 iteraciones

function c_log = mseg_clase(h_ref, beta, Niters, g)
% El lazo de MSEG.m, sin tocar, con la entrada pasada por g
Ncoeff = length(h_ref);  Lbatch = 1e3;  Delta_tap = 1e-4;
c = zeros(Ncoeff, 1);  gradient = zeros(Ncoeff, 1);
c_log = zeros(Niters+1, Ncoeff);
for iter = 1:Niters
    x = filter(g, 1, randn(Lbatch, 1));
    y_ref = filter(h_ref, 1, x);
    y = filter(c, 1, x);
    mse0 = mean(abs(y - y_ref).^2);
    for nc = 1:Ncoeff
        caux = c;  caux(nc) = caux(nc) + Delta_tap;
        y = filter(caux, 1, x);
        mse1 = mean(abs(y - y_ref).^2);
        gradient(nc) = (mse1 - mse0)/Delta_tap;
    end
    c = c - beta*gradient;
    c_log(iter+1, :) = c;
end
end

El bloque 2 cruza la simulación contra la forma cerrada. Con $\mathrm{\Phi }=I$ y ${\xi }_{\text{min}}=0$ —el filtro de 5 taps puede copiar exactamente a h_ref—, el MSE es ${\parallel q_j\parallel }^2$ y cae como ${\mathrm{0,998}}^{2j}{\parallel h_{\text{ref}}\parallel }^2$, desde ${\parallel h_{\text{ref}}\parallel }^2=\mathrm{0,3555}$ (−4,49 dB): una recta de $20{\log }_{10}\mathrm{0,998}$ por iteración, −17,39 dB cada mil. Medida entre las iteraciones 200 y 3000 da −17,28. El bloque 5 hace lo mismo con la entrada coloreada del canal [1 0.9], $\chi =\mathrm{13,41}$: el modo lento se achica 0,9496 por iteración con $\beta =\mathrm{0,2}$ (predicho $1-\mathrm{0,2}·\mathrm{0,2512}=\mathrm{0,9498}$) y 0,8677 con ${\beta }_{\mathrm{opt}}$ (predicho 0,8612). Un 10 % por encima de $2/{\lambda }_{\text{max}}=\mathrm{0,5937}$ el modo rápido crece 1,162 por iteración, contra 1,200 predicho. Para bajar el error cien veces hacen falta 30 iteraciones con el mejor paso; con entrada blanca, 2.

Cuatro detalles que muerden:

• El beta del script es la mitad del $\beta$ de la filmina. La filmina resta $(\beta /2)\nabla$ y el script beta*gradient. Con entrada blanca la razón es 1 - 2*beta y el límite, beta < 1/lambda_max. Probar beta = 1.5 «porque el límite es 2» diverge.
• El piso de −79 dB no es ruido: es sesgo. El MSE de un lote es cuadrático en los taps, así que la diferencia hacia adelante da el gradiente más un término fijo, $\Delta ·{[\hat{\mathrm{\Phi }}]}_{nn}$. El algoritmo converge a donde ese gradiente corrido se anula, $\Delta /2$ por debajo de cada tap con entrada blanca, y el MSE se estanca en $N{(\Delta /2)}^2$: −79,03 dB, y se mide −79,05. Con la diferencia centrada, (mse(c + Delta) - mse(c - Delta))/(2*Delta), el sesgo desaparece, porque para una cuadrática es exacta.
• Las $N+1$ evaluaciones del MSE de una iteración usan el mismo lote. Con lotes distintos, la diferencia mse1 - mse0 la dominaría lo que varía el MSE de un lote a otro, muchísimo más grande que el efecto de mover un tap en 10⁻⁴.
• El gradiente de un lote no es el exacto. Por eso $\beta =1$ no llega en un paso con entrada blanca, y por eso el modo rápido divergente crece 1,162 y no 1,200.

Para tener a mano

1. El gradiente del MSE es $2\mathrm{\Phi }c-2\alpha =-2E[E_k{r}_k^{\ast }]$: sale de derivar el término lineal y el cuadrático, con $\mathrm{\Phi }$ simétrica.
2. El MSEG: $c_{j+1}=c_j-\frac{\beta }{2}\nabla =(I-\beta \mathrm{\Phi })c_j+\beta \alpha$. $j$ cuenta iteraciones, no símbolos.
3. El error de coeficientes: $q_j={(I-\beta \mathrm{\Phi })}^j{q}_0$. En la base de autovectores, cada modo es ${(1-\beta {\lambda }_i)}^j{u}_{i,0}$, y el exceso de MSE es ${\sum }_i{\lambda }_i{(1-\beta {\lambda }_i)}^{2j}{|u_{i,0}|}^2$.
4. Converge si y sólo si $0<\beta <2/{\lambda }_{\text{max}}$. Cota segura sin autovalores: $\beta <2/(N{\phi }_0)$.
5. El paso más rápido es ${\beta }_{\mathrm{opt}}=2/({\lambda }_{\text{min}}+{\lambda }_{\text{max}})$, con razón $(\chi -1)/(\chi +1)$. Para $\chi$ grande, unas $\mathrm{3,45}\chi$ iteraciones por cada factor mil.
6. $\chi$ está acotada por el cociente entre el máximo y el mínimo del espectro de la entrada: un canal con pozos profundos se aprende lento, y ningún $\beta$ lo arregla.
7. En MSEG.m, beta es $\beta /2$, el gradiente sale de un lote de 1000 muestras, y el piso de −79 dB es el sesgo de la diferencia hacia adelante.