¿Por qué hay coeficientes grandes para polinomios de orden superior?

En el libro de Bishop sobre aprendizaje automático, se analiza el problema de ajustar una función polinómica a un conjunto de puntos de datos.

Sea M el orden del polinomio ajustado. Dice como eso

Vemos que, a medida que M aumenta, la magnitud de los coeficientes generalmente aumenta. En particular para el polinomio M = 9, los coeficientes se han ajustado finamente a los datos desarrollando grandes valores positivos y negativos para que la función polinómica correspondiente coincida exactamente con cada uno de los puntos de datos, pero entre puntos de datos (particularmente cerca de los extremos de rango) la función exhibe las grandes oscilaciones.

No entiendo por qué los valores grandes implican un ajuste más preciso de los puntos de datos. Creo que los valores serían más precisos después del decimal para un mejor ajuste.

Este es un problema bien conocido con polinomios de alto orden, conocido como fenómeno de Runge . Numéricamente está asociado con un mal acondicionamiento de la matriz de Vandermonde , lo que hace que los coeficientes sean muy sensibles a pequeñas variaciones en los datos y / o redondeo en los cálculos (es decir, el modelo no es establemente identificable ). Vea también esta respuesta en el SciComp SE.

Hay muchas soluciones a este problema, por ejemplo , la aproximación de Chebyshev , el suavizado de splines y la regularización de Tikhonov . La regularización de Tikhonov es una generalización de la regresión de crestas , penalizando una norma del vector de coeficientes θ , donde para alisar la matriz de peso Λ es un operador derivado. Para penalizar las oscilaciones, puede usar Λ θ = p ′ ′ [ x ] , donde p [ x $||\Lambda \theta]||$$\theta$$\Lambda$$\Lambda \theta=p^{\prime\prime}[x]$$p[x]$ es el polinomio evaluado en los datos.

EDITAR: La respuesta del usuario hxd1011 señala que algunos de los problemas numéricos de mal acondicionamiento pueden abordarse utilizando polinomios ortogonales, lo cual es un buen punto. Sin embargo, quisiera señalar que aún persisten los problemas de identificación con los polinomios de alto orden. Es decir, el mal acondicionamiento numérico se asocia con la sensibilidad a las perturbaciones "infinitesimales" (por ejemplo, redondeo), mientras que el mal acondicionamiento "estadístico" se refiere a la sensibilidad a las perturbaciones "finitas" (por ejemplo, valores atípicos; el problema inverso está mal planteado ).

Los métodos mencionados en mi segundo párrafo están relacionados con esta sensibilidad atípica . Puede pensar en esta sensibilidad como una violación del modelo de regresión lineal estándar, que al usar un desajuste asume implícitamente que los datos son gaussianos. La regularización de Splines y Tikhonov trata esta sensibilidad atípica al imponer una suavidad previa al ajuste. Ofertas de aproximación de Chebyshev con este mediante el uso de una L ∞ inadaptado aplican sobre el dominio continuo , es decir, no sólo en los puntos de datos. Aunque los polinomios de Chebyshev son ortogonales (wrt un cierto producto interno ponderado), creo que si se usan con un desajuste L 2 sobre los datos,$L_2$$L_{\infty}$$L_2$ todavía tener sensibilidad atípica.

Lo primero que debe verificar es si el autor está hablando de polinomios sin procesar frente a polinomios ortogonales .

Para polinomios ortogonales. el coeficiente no se está volviendo "más grande".

Aquí hay dos ejemplos de expansión polinómica de segundo y décimo quinto orden. Primero mostramos el coeficiente para la expansión de segundo orden.

summary(lm(mpg~poly(wt,2),mtcars))

Call:
lm(formula = mpg ~ poly(wt, 2), data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-3.483 -1.998 -0.773  1.462  6.238 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   20.0906     0.4686  42.877  < 2e-16 ***
poly(wt, 2)1 -29.1157     2.6506 -10.985 7.52e-12 ***
poly(wt, 2)2   8.6358     2.6506   3.258  0.00286 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.651 on 29 degrees of freedom
Multiple R-squared:  0.8191,    Adjusted R-squared:  0.8066 
F-statistic: 65.64 on 2 and 29 DF,  p-value: 1.715e-11

Luego mostramos el decimoquinto orden.

summary(lm(mpg~poly(wt,15),mtcars))

Call:
lm(formula = mpg ~ poly(wt, 15), data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3233 -0.4641  0.0072  0.6401  4.0394 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     20.0906     0.4551  44.147  < 2e-16 ***
poly(wt, 15)1  -29.1157     2.5743 -11.310 4.83e-09 ***
poly(wt, 15)2    8.6358     2.5743   3.355  0.00403 ** 
poly(wt, 15)3    0.2749     2.5743   0.107  0.91629    
poly(wt, 15)4   -1.7891     2.5743  -0.695  0.49705    
poly(wt, 15)5    1.8797     2.5743   0.730  0.47584    
poly(wt, 15)6   -2.8354     2.5743  -1.101  0.28702    
poly(wt, 15)7    2.5613     2.5743   0.995  0.33459    
poly(wt, 15)8    1.5772     2.5743   0.613  0.54872    
poly(wt, 15)9   -5.2412     2.5743  -2.036  0.05866 .  
poly(wt, 15)10  -2.4959     2.5743  -0.970  0.34672    
poly(wt, 15)11   2.5007     2.5743   0.971  0.34580    
poly(wt, 15)12   2.4263     2.5743   0.942  0.35996    
poly(wt, 15)13  -2.0134     2.5743  -0.782  0.44559    
poly(wt, 15)14   3.3994     2.5743   1.320  0.20525    
poly(wt, 15)15  -3.5161     2.5743  -1.366  0.19089    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.574 on 16 degrees of freedom
Multiple R-squared:  0.9058,    Adjusted R-squared:  0.8176 
F-statistic: 10.26 on 15 and 16 DF,  p-value: 1.558e-05

Tenga en cuenta que estamos utilizando polinomios ortogonales , por lo que el coeficiente del orden inferior es exactamente el mismo que los términos correspondientes en los resultados del orden superior. Por ejemplo, la intersección y el coeficiente para primer orden es 20.09 y -29.11 para ambos modelos.

$10^6$

> summary(lm(mpg~poly(wt,15, raw=T),mtcars))

Call:
lm(formula = mpg ~ poly(wt, 15, raw = T), data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6217 -0.7544  0.0306  1.1678  5.4308 

Coefficients: (3 not defined because of singularities)
                          Estimate Std. Error t value Pr(>|t|)
(Intercept)              6.287e+05  9.991e+05   0.629    0.537
poly(wt, 15, raw = T)1  -2.713e+06  4.195e+06  -0.647    0.526
poly(wt, 15, raw = T)2   5.246e+06  7.893e+06   0.665    0.514
poly(wt, 15, raw = T)3  -6.001e+06  8.784e+06  -0.683    0.503
poly(wt, 15, raw = T)4   4.512e+06  6.427e+06   0.702    0.491
poly(wt, 15, raw = T)5  -2.340e+06  3.246e+06  -0.721    0.480
poly(wt, 15, raw = T)6   8.537e+05  1.154e+06   0.740    0.468
poly(wt, 15, raw = T)7  -2.184e+05  2.880e+05  -0.758    0.458
poly(wt, 15, raw = T)8   3.809e+04  4.910e+04   0.776    0.447
poly(wt, 15, raw = T)9  -4.212e+03  5.314e+03  -0.793    0.438
poly(wt, 15, raw = T)10  2.382e+02  2.947e+02   0.809    0.429
poly(wt, 15, raw = T)11         NA         NA      NA       NA
poly(wt, 15, raw = T)12 -5.642e-01  6.742e-01  -0.837    0.413
poly(wt, 15, raw = T)13         NA         NA      NA       NA
poly(wt, 15, raw = T)14         NA         NA      NA       NA
poly(wt, 15, raw = T)15  1.259e-04  1.447e-04   0.870    0.395

Residual standard error: 2.659 on 19 degrees of freedom
Multiple R-squared:  0.8807,    Adjusted R-squared:  0.8053 
F-statistic: 11.68 on 12 and 19 DF,  p-value: 2.362e-06

Abhishek, tiene razón en que mejorar la precisión de los coeficientes mejorará la precisión.

Vemos que, a medida que M aumenta, la magnitud de los coeficientes generalmente aumenta. En particular, para el polinomio M = 9, los coeficientes se han ajustado con precisión a los datos mediante el desarrollo de grandes valores positivos y negativos para que la función polinómica correspondiente coincida exactamente con cada uno de los puntos de datos, pero entre puntos de datos (particularmente cerca de los extremos de los datos). rango) la función exhibe grandes oscilaciones.

Creo que el problema de la magnitud es bastante irrelevante para el punto general de Bishop: que el uso de un modelo complicado con datos limitados conduce a un 'sobreajuste'. En su ejemplo, se usan 10 puntos de datos para estimar un polinomio de 9 dimensiones (es decir, 10 variables y 10 incógnitas).

Si ajustamos una onda sinusoidal (sin ruido), entonces el ajuste funciona perfectamente, ya que las ondas sinusoidales [en un intervalo fijo] pueden aproximarse con precisión arbitraria utilizando polinomios. Sin embargo, en el ejemplo de Bishop tenemos una cierta cantidad de "ruido" que no deberíamos encajar. La forma en que hacemos esto es manteniendo el número de puntos de datos al número de variables del modelo (coeficientes polinómicos) grandes o mediante la regularización.

La regularización impone restricciones 'blandas' al modelo (por ejemplo, en la regresión de cresta), la función de costo que intenta minimizar es una combinación de 'error de ajuste' y complejidad del modelo: por ejemplo, en la regresión de cresta, la complejidad se mide por la suma de coeficientes cuadrados En efecto, esto impone un costo en la reducción del error: el aumento de los coeficientes solo se permitirá si tiene una reducción lo suficientemente grande en el error de ajuste [qué tan grande es lo suficientemente grande se especifica mediante un multiplicador en el término de complejidad del modelo]. Por lo tanto, la esperanza es que al elegir el multiplicador apropiado no se ajustará a un término de ruido pequeño adicional, ya que la mejora en el ajuste no justifica el aumento de los coeficientes.

Preguntó por qué los coeficientes grandes mejoran la calidad del ajuste. Esencialmente, la razón es que la función estimada (sin + ruido) no es un polinomio, y los grandes cambios en la curvatura requeridos para aproximar el efecto de ruido con polinomios requieren coeficientes grandes.

Tenga en cuenta que el uso de polinomios ortogonales no tiene ningún efecto (he agregado un desplazamiento de 0.1 solo para que los polinomios ortogonales y sin procesar no estén uno encima del otro)

require (penalized)
poly_order<-9
x_long<-seq(0,1, length.out = 100)
nx<-10
x<-seq(0,1, length.out = nx)
noise<- rnorm(nx, 0, 1)
noise_scale<-0.2
y<-sin(2*pi*x)+noise_scale*noise

training_data<-data.frame(x=x,y=y)
y_long<-sin(2*pi*x_long)

plot(x,y, col ='blue',ylim=c(-1.5,1.5))
lines(x_long,y_long,col='green')

polyfit_raw<-lm(y~poly(x,poly_order,raw=TRUE),data=training_data)
summary(polyfit_raw)

polyfit_raw_ridge1<-penalized(y,~poly(x,poly_order,raw=TRUE), model='linear', data=training_data, lambda2=0.0001, maxiter=10000, standardize=TRUE)

polyfit_orthog<-lm(y~poly(x,poly_order),data=training_data)
summary(polyfit_orthog)

pred_raw<-predict(polyfit_raw,data.frame(x=x_long))
pred_ortho<-predict(polyfit_orthog,data.frame(x=x_long))
pred_raw_ridge<-predict(polyfit_raw_ridge1,data=data.frame(x=x_long))[,'mu']
lines(x_long,pred_raw,col='red')
# add 0.1 offset to make visible
lines(x_long,pred_ortho+0.1,col='black')
lines(x_long,pred_raw_ridge,col='purple')
legend("bottomleft",legend=c('data sin(2 pi x) + noise','sin(2 pi x)', 
                             'raw poly','orthog poly +0.1 offset','raw poly + ridge regression'),
       fill=c('blue','green','red','black','purple'))