Esta es una respuesta larga . Así que, vamos a dar una versión corta de la historia aquí.
- No hay una buena solución algebraica para este problema de encontrar la raíz, así que necesitamos un algoritmo numérico.
- La función $ \mathrm {df}( \lambda )$ tiene muchas propiedades bonitas. Podemos aprovecharlas para crear una versión especializada del método de Newton para este problema con garantía monótona convergencia a cada raíz.
- Incluso con muerte cerebral
R sin ningún intento de optimización puede calcular una cuadrícula de tamaño 100 con $p = 100\,000$ en unos pocos segundos. Cuidadosamente escrito C el código reduciría esto en al menos 2-3 órdenes de magnitud.
A continuación se presentan dos esquemas para garantizar la convergencia monótona. Uno utiliza los límites que se muestran a continuación, que parecen ayudar a salvar un paso o dos de Newton en ocasiones.
Ejemplo : $p = 100\,000$ y una cuadrícula uniforme para los grados de libertad del tamaño 100. Los valores propios están distribuidos por Pareto, por lo que están muy sesgados. Abajo hay tablas con el número de pasos de Newton para encontrar cada raíz.
# Table of Newton iterations per root.
# Without using lower-bound check.
1 3 4 5 6
1 28 65 5 1
# Table with lower-bound check.
1 2 3
1 14 85
No habrá una solución de forma cerrada para esto en general, pero hay es mucha estructura presente que puede ser utilizada para producir soluciones muy efectivas y seguras usando métodos estándar de búsqueda de raíces.
Antes de profundizar en las cosas, vamos a recoger algunas propiedades y consecuencias de la función $$ \newcommand { \df }{ \mathrm {df}} \df ( \lambda ) = \sum_ {i=1}^p \frac {d_i^2}{d_i^2 + \lambda } \>. $$
Propiedad 0 : $ \df $ es una función racional de $ \lambda $ . (Esto se desprende de la definición.)
Consecuencia 0 : No existe una solución algebraica general para encontrar la raíz $ \df ( \lambda ) - y = 0$ . Esto es porque hay un equivalente polinomio el problema de la raíz del grado $p$ y por lo tanto si $p$ no es extremadamente pequeño (es decir, menos de cinco), no existirá una solución general. Por lo tanto, necesitaremos un método numérico.
Propiedad 1 : La función $ \df $ es convexo y está disminuyendo en $ \lambda \geq 0$ . (Toma los derivados.)
Consecuencia 1(a) : El algoritmo de búsqueda de raíces de Newton se comportará muy muy bien en esta situación. Deje que $y$ los grados de libertad deseados y $ \lambda_0 $ la raíz correspondiente, es decir, $y = \df ( \lambda_0 )$ . En particular, si empezamos con cualquier valor inicial $ \lambda_1 < \lambda_0 $ (así que, $ \df ( \lambda_1 ) > y$ ), luego la secuencia de iteraciones de pasos de Newton $ \lambda_1 , \lambda_2 , \ldots $ convergerán monótonamente a la solución única $ \lambda_0 $ .
Consecuencia 1(b) : Además, si empezamos con $ \lambda_1 > \lambda_0 $ entonces el primero paso daría lugar a $ \lambda_2 \leq \lambda_0 $ de donde aumentará monótonamente a la solución por la consecuencia anterior (véase la advertencia más abajo). Intuitivamente, este último hecho sigue porque si comenzamos a la derecha de la raíz, el derivado es "demasiado" superficial debido a la convexidad de $ \df $ y así el primer paso de Newton nos llevará a algún lugar a la izquierda de la raíz. NB Desde $ \df $ es no en general convexo por negativo $ \lambda $ esto proporciona una fuerte razón para preferir comenzar a la izquierda de la raíz deseada. De lo contrario, tenemos que comprobar que el paso de Newton no ha dado lugar a un valor negativo para la raíz estimada, lo que puede situarnos en algún lugar de una porción no convexa de $ \df $ .
Consecuencia 1(c) : Una vez que hayamos encontrado la raíz de algunos $y_1$ y luego están buscando la raíz de algunos $y_2 < y_1$ usando $ \lambda_1 $ de tal manera que $ \df ( \lambda_1 ) = y_1$ como nuestra suposición inicial garantiza que empezamos a la izquierda de la segunda raíz. Por lo tanto, nuestra convergencia está garantizada para ser monótona desde allí.
Propiedad 2 : Existen límites razonables para dar puntos de partida "seguros". Usando argumentos de convexidad y la desigualdad de Jensen, tenemos los siguientes límites $$ \frac {p}{1+ \frac { \lambda }{p} \sum d_i^{-2}} \leq \df ( \lambda ) \leq \frac {p \sum_i d_i^2}{ \sum_i d_i^2 + p \lambda } \>. $$ Consecuencia 2 : Esto nos dice que la raíz $ \lambda_0 $ satisfactoria $ \df ( \lambda_0 ) = y$ obedece $$ \frac {1}{ \frac {1}{p} \sum_i d_i^{-2}} \left ( \frac {p - y}{y} \right ) \leq \lambda_0 \leq \left ( \frac {1}{p} \sum_i d_i^2 \right ) \left ( \frac {p - y}{y} \right ) \>. \tag {$ \star $} $$ Así que, hasta una constante común, hemos metido la raíz entre los medios armónicos y aritméticos de la $d_i^2$ .
Esto supone que $d_i > 0$ para todos $i$ . Si este no es el caso, entonces el mismo límite se mantiene considerando sólo lo positivo $d_i$ y reemplazando $p$ por el número de positivos $d_i$ . NB : Desde $ \df (0) = p$ asumiendo que todos $d_i > 0$ Entonces $y \in (0,p]$ de donde los límites son siempre no triviales (por ejemplo, el límite inferior es siempre no negativo).
Aquí está la trama de un ejemplo "típico" de $ \df ( \lambda )$ con $p = 400$ . Hemos superpuesto una cuadrícula de tamaño 10 para los grados de libertad. Estas son las líneas horizontales en la trama. Las líneas verdes verticales corresponden al límite inferior en $( \star )$ .
![Example dof plot with grid and bounds]()
Un algoritmo y algún código R de ejemplo
Un algoritmo muy eficiente dado una cuadrícula de grados de libertad deseados $y_1, \ldots y_n$ en $(0,p]$ es clasificarlos en orden decreciente y luego secuencialmente encontrar la raíz de cada uno, usando la raíz anterior como punto de partida para la siguiente. Podemos afinar más esto comprobando si cada raíz es mayor que el límite inferior de la siguiente y, si no, podemos empezar la siguiente iteración en el límite inferior.
Aquí hay un ejemplo de código en R sin intentar optimizarlo. Como se ve a continuación, sigue siendo bastante rápido aunque R es, por decirlo educadamente, horriblemente, terriblemente lento en los bucles.
# Newton's step for finding solutions to regularization dof.
dof <- function(lambda, d) { sum(1/(1+lambda / (d[d>0])^2)) }
dof.prime <- function(lambda, d) { -sum(1/(d[d>0]+lambda / d[d>0])^2) }
newton.step <- function(lambda, y, d)
{ lambda - (dof(lambda,d)-y)/dof.prime(lambda,d) }
# Full Newton step; Finds the root of y = dof(lambda, d).
newton <- function(y, d, lambda = NA, tol=1e-10, smart.start=T)
{
if( is.na(lambda) || smart.start )
lambda <- max(ifelse(is.na(lambda),0,lambda), (sum(d>0)/y-1)/mean(1/(d[d>0])^2))
iter <- 0
yn <- Inf
while( abs(y-yn) > tol )
{
lambda <- max(0, newton.step(lambda, y, d)) # max = pedantically safe
yn <- dof(lambda,d)
iter = iter + 1
}
return(list(lambda=lambda, dof=y, iter=iter, err=abs(y-yn)))
}
A continuación se muestra el algoritmo completo final que toma una cuadrícula de puntos, y un vector de la $d_i$ ( no $d_i^2$ !).
newton.grid <- function(ygrid, d, lambda=NA, tol=1e-10, smart.start=TRUE)
{
p <- sum(d>0)
if( any(d < 0) || all(d==0) || any(ygrid > p)
|| any(ygrid <= 0) || (!is.na(lambda) && lambda < 0) )
stop("Don't try to fool me. That's not nice. Give me valid inputs, please.")
ygrid <- sort(ygrid, decreasing=TRUE)
out <- data.frame()
lambda <- NA
for(y in ygrid)
{
out <- rbind(out, newton(y,d,lambda, smart.start=smart.start))
lambda <- out$lambda[nrow(out)]
}
out
}
Muestra de llamada de función
set.seed(17)
p <- 100000
d <- sqrt(sort(exp(rexp(p, 10)),decr=T))
ygrid <- p*(1:100)/100
# Should take ten seconds or so.
out <- newton.grid(ygrid,d)
