programacion - método del gradiente matlab
El descenso de gradiente y el método de ecuación normal para resolver la regresión lineal ofrecen diferentes soluciones (4)
Estoy trabajando en un problema de aprendizaje automático y quiero usar la regresión lineal como algoritmo de aprendizaje. He implementado 2 métodos diferentes para encontrar los parámetros theta del modelo de regresión lineal: gradiente (más pronunciado) descenso y ecuación normal. En los mismos datos, ambos deberían dar un vector theta aproximadamente igual. Sin embargo no lo hacen.
Ambos vectores theta son muy similares en todos los elementos, pero el primero. Ese es el que se usa para multiplicar el vector de todos los 1 agregados a los datos.
Así es como se ven los theta (la primera columna es la salida del descenso del gradiente, la segunda salida de la ecuación normal):
Grad desc Norm eq
-237.7752 -4.6736
-5.8471 -5.8467
9.9174 9.9178
2.1135 2.1134
-1.5001 -1.5003
-37.8558 -37.8505
-1.1024 -1.1116
-19.2969 -19.2956
66.6423 66.6447
297.3666 296.7604
-741.9281 -744.1541
296.4649 296.3494
146.0304 144.4158
-2.9978 -2.9976
-0.8190 -0.8189
¿Qué puede causar la diferencia en theta(1, 1) devuelta por el descenso del gradiente en comparación con theta(1, 1) devuelta por la ecuación normal? ¿Tengo errores en mi código?
Aquí está mi implementación de la ecuación normal en Matlab:
function theta = normalEque(X, y)
[m, n] = size(X);
X = [ones(m, 1), X];
theta = pinv(X''*X)*X''*y;
end
Aquí hay un código para el descenso del gradiente:
function theta = gradientDesc(X, y)
options = optimset(''GradObj'', ''on'', ''MaxIter'', 9999);
[theta, ~, ~] = fminunc(@(t)(cost(t, X, y)),...
zeros(size(X, 2), 1), options);
end
function [J, grad] = cost(theta, X, y)
m = size(X, 1);
X = [ones(m, 1), X];
J = sum((X * theta - y) .^ 2) ./ (2*m);
for i = 1:size(theta, 1)
grad(i, 1) = sum((X * theta - y) .* X(:, i)) ./ m;
end
end
Paso exactamente los mismos datos X e y a ambas funciones (no normalizo X ).
Editar 1:
Basado en respuestas y comentarios, revisé mi código y realicé algunas pruebas.
Primero quiero verificar si el problema puede ser causado porque X está cerca de singular como lo sugiere la respuesta de @ user1489497 . Así que reemplacé pinv by inv - y cuando lo ejecuté realmente advertí que Matrix is close to singular or badly scaled. . Para estar seguro de que ese no es el problema, obtuve un conjunto de datos mucho más grande y realizo pruebas con este nuevo conjunto de datos. Esta vez, inv(X) no pinv la advertencia y el uso de pinv e inv dieron los mismos resultados. Así que espero que X ya no esté cerca del singular .
Luego cambié el código normalEque como lo sugieren las astillas de madera, así que ahora parece que:
function theta = normalEque(X, y)
X = [ones(size(X, 1), 1), X];
theta = pinv(X)*y;
end
Sin embargo, el problema sigue ahí . La nueva función normalEque en los datos nuevos que no están cerca de singular da diferente theta como gradientDesc .
Para averiguar qué algoritmo tiene errores, he ejecutado el algoritmo de regresión lineal del software Weka de minería de datos con los mismos datos. Weka calculó theta muy similar a la salida de normalEque pero diferente a la salida de gradientDesc . Así que supongo que normalEque es correcto y hay un error en gradientDesc .
Aquí hay una comparación de normalEque calculada por Weka, normalEque y GradientDesc :
Weka(correct) normalEque gradientDesc
779.8229 779.8163 302.7994
1.6571 1.6571 1.7064
1.8430 1.8431 2.3809
-1.5945 -1.5945 -1.5964
3.8190 3.8195 5.7486
-4.8265 -4.8284 -11.1071
-6.9000 -6.9006 -11.8924
-15.6956 -15.6958 -13.5411
43.5561 43.5571 31.5036
-44.5380 -44.5386 -26.5137
0.9935 0.9926 1.2153
-3.1556 -3.1576 -1.8517
-0.1927 -0.1919 -0.6583
2.9207 2.9227 1.5632
1.1713 1.1710 1.1622
0.1091 0.1093 0.0084
1.5768 1.5762 1.6318
-1.3968 -1.3958 -2.1131
0.6966 0.6963 0.5630
0.1990 0.1990 -0.2521
0.4624 0.4624 0.2921
-12.6013 -12.6014 -12.2014
-0.1328 -0.1328 -0.1359
También calculé los errores según lo sugerido por la respuesta de Justin Peel . La salida de normalEque da un error al cuadrado ligeramente menor, pero la diferencia es marginal. Lo que es más cuando calculo el gradiente de costo de theta usando el cost función (el mismo que el usado por gradientDesc ) obtuve un gradiente cercano a cero . Lo mismo hecho en la salida de gradientDesc no da gradiente cerca de cero. Esto es lo que quiero decir:
>> [J_gd, grad_gd] = cost(theta_gd, X, y, size(X, 1));
>> [J_ne, grad_ne] = cost(theta_ne, X, y, size(X, 1));
>> disp([J_gd, J_ne])
120.9932 119.1469
>> disp([grad_gd, grad_ne])
-0.005172856743846 -0.000000000908598
-0.026126463200876 -0.000000135414602
-0.008365136595272 -0.000000140327001
-0.094516503056041 -0.000000169627717
-0.028805977931093 -0.000000045136985
-0.004761477661464 -0.000000005065103
-0.007389474786628 -0.000000005010731
0.065544198835505 -0.000000046847073
0.044205371015018 -0.000000046169012
0.089237705611538 -0.000000046081288
-0.042549228192766 -0.000000051458654
0.016339232547159 -0.000000037654965
-0.043200042729041 -0.000000051748545
0.013669010209370 -0.000000037399261
-0.036586854750176 -0.000000027931617
-0.004761447097231 -0.000000027168798
0.017311225027280 -0.000000039099380
0.005650124339593 -0.000000037005759
0.016225097484138 -0.000000039060168
-0.009176443862037 -0.000000012831350
0.055653840638386 -0.000000020855391
-0.002834810081935 -0.000000006540702
0.002794661393905 -0.000000032878097
Esto sugeriría que el descenso de gradiente simplemente no converge al mínimo global ... Pero ese no es el caso, ya que lo ejecuto durante miles de iteraciones. Entonces, ¿dónde está el error?
Como otros mencionaron, una matriz de arpillera mal condicionada es probablemente la causa de su problema.
El número de pasos que toman los algoritmos de gradiente de gradiente estándar para alcanzar un óptimo local es una función del valor propio más grande de la arpillera dividido por el más pequeño (esto se conoce como el número de condición del hessiano). Por lo tanto, si su matriz está mal acondicionada, podría tomar un número extremadamente grande de iteraciones para que el gradiente de descenso converja a un nivel óptimo. (Para el caso singular, podría converger a muchos puntos, por supuesto).
Sugeriría probar tres cosas diferentes para verificar que un algoritmo de optimización no restringido funcione para su problema (lo cual debería): 1) Generar algunos datos sintéticos calculando el resultado de una función lineal conocida para entradas aleatorias y agregando una pequeña cantidad de ruido gaussiano . Asegúrese de tener muchos más puntos de datos que dimensiones. Esto debería producir una arpillera no singular. 2) Agregue un término de regularización a su función de error para aumentar el número de condición de la arpillera. 3) Use un método de segundo orden como degradado conjugado o L-BFGS en lugar de descenso de gradiente para reducir el número de pasos necesarios para que el algoritmo converja. (Probablemente necesite hacer esto junto con el n. ° 2).
¿Podría publicar un poco más sobre cómo luce X? Estás usando pinv () que es pseudo inverso de Moore-Penrose. Si la matriz está mal acondicionada, esto podría causar problemas para obtener la inversa. Apuesto a que el método de gradiente de descenso está más cerca de la marca.
Debería ver qué método le está dando el error más pequeño. Eso indicará qué método está luchando. Sospecho que el método de ecuación normal es la solución problemática porque si X está mal acondicionado, entonces puede haber algunos problemas allí.
Probablemente deberías reemplazar tu solución de ecuación normal con theta = X/y que usará un método de descomposición QR para resolverlo.
Finalmente tuve tiempo de volver a esto. No hay "error".
Si la matriz es singular, entonces hay infinitas soluciones. Puede elegir cualquier solución de ese conjunto y obtener igualmente una buena respuesta. La solución pinv (X) * y es una buena que muchos prefieren porque es la solución de norma mínima.
NUNCA hay una buena razón para usar inv (X) * y. Peor aún, es usar inversa en las ecuaciones normales, por lo tanto inv (X ''* X) * X'' * y es simplemente basura numérica. No me importa quién te haya dicho que lo uses, te están guiando al lugar equivocado. (Sí, funcionará aceptablemente para problemas que están bien acondicionados, pero la mayoría de las veces no sabes cuándo te va a dar una porquería. Entonces, ¿para qué usarlo?)
Las ecuaciones normales son en general algo malo, INCLUSO si estás resolviendo un problema regularizado. Hay formas de hacerlo que evitan cuadrar el número de condición del sistema, aunque no las explicaré a menos que se me pida, ya que esta respuesta ha sido lo suficientemente larga.
X / y también producirá un resultado que es razonable.
ABSOLUTAMENTE no hay una buena razón para lanzar un optimizador sin restricciones al problema, ya que esto arrojará resultados que son inestables, completamente dependientes de sus valores iniciales.
Como ejemplo, comenzaré con un problema singular.
X = repmat([1 2],5,1);
y = rand(5,1);
>> X/y
Warning: Rank deficient, rank = 1, tol = 2.220446e-15.
ans =
0
0.258777984694222
>> pinv(X)*y
ans =
0.103511193877689
0.207022387755377
pinv y backslash devuelven soluciones ligeramente diferentes. Como resultado, hay una solución básica, a la que podemos agregar CUALQUIER cantidad del vector de espacio nulo para el espacio de fila de X.
null(X)
ans =
0.894427190999916
-0.447213595499958
pinv genera la solución de norma mínima. De todas las soluciones que podrían haber resultado, esta tiene una norma mínima de 2.
Por el contrario, la barra diagonal inversa genera una solución que tendrá una o más variables establecidas en cero.
Pero si usa un optimizador sin restricciones, generará una solución que depende completamente de sus valores iniciales. Una vez más, ANY puede agregar a su solución la cantidad de ese vector nulo, y todavía tiene una solución completamente válida, con el mismo valor de la suma de cuadrados de errores.
Tenga en cuenta que aunque no se devuelve ningún waring de singularidad, esto no necesariamente significa que su matriz no está cerca del singular. Has cambiado poco sobre el problema, por lo que TODAVÍA está cerca, pero no lo suficiente como para activar la advertencia.