Matriz de Hesse inversa de retorno al final del entrenamiento de DNN y derivados parciales wrt las entradas
Utilizando Keras y Tensorflow como back-end, he creado un DNN que toma espectros estelares como entrada (7213 puntos de datos) y genera tres parámetros estelares (temperatura, gravedad y metalicidad). La red se entrena bien y predice bien en mis conjuntos de prueba, pero para que los resultados sean científicamente útiles, necesito poder estimar mis errores. El primer paso para hacer esto es obtener la matriz de Hesse inversa, que no parece posible usando solo Keras. Por lo tanto, intento crear una solución alternativa con scipy, utilizando scipy.optimize.minimize con BFGS, L-BFGS-B o Netwon-CG como método. Cualquiera de estos devolverá la matriz de Hesse inversa.
La idea es entrenar el modelo usando el optimizador Adam para 100 épocas (o hasta que el modelo converja) y luego ejecutar una única iteración o función de BFGS (o una de las otras) para devolver la matriz de Hesse de mi modelo.
Aquí está mi código:
from scipy.optimize import minimize
import numpy as np
from keras.models import Sequential
from keras.layers import Dense, Activation
from keras.optimizers import Adam
# Define vars
activation = 'relu'
init = 'he_normal'
beta_1 = 0.9
beta_2 = 0.999
epsilon = 1e-08
input_shape = (None,n)
n_hidden = [2048,1024,512,256,128,32]
output_dim = 3
epochs = 100
lr = 0.0008
batch_size = 64
decay = 0.00
# Design DNN Layers
model = Sequential([
Dense(n_hidden[0], batch_input_shape=input_shape, init=init, activation=activation),
Dense(n_hidden[1], init=init, activation=activation),
Dense(n_hidden[2], init=init, activation=activation),
Dense(n_hidden[3], init=init, activation=activation),
Dense(n_hidden[4], init=init, activation=activation),
Dense(n_hidden[5], init=init, activation=activation),
Dense(output_dim, init=init, activation='linear'),
])
# Optimization function
optimizer = Adam(lr=lr, beta_1=beta_1, beta_2=beta_2, epsilon=epsilon, decay=decay)
# Compile and train network
model.compile(optimizer=optimizer, loss='mean_squared_error')
#train_X.shape = (50000,7213)
#train_Y.shape = (50000,3)
#cv_X.shape = (10000,7213)
#cv_Y.shape = (10000,3)
history = model.fit(train_X, train_Y, validation_data=(cv_X, cv_Y),
nb_epoch=epochs, batch_size=batch_size, verbose=2)
weights = []
for layer in model.layers:
weights.append(layer.get_weights())
def loss(W):
weightsList = W
weightsList = np.array(W)
new_weights = []
for i, layer in enumerate((weightsList)):
new_weights.append(np.array(weightsList[i]))
model.set_weights(np.array(new_weights))
preds = model.predict(train_X)
mse = np.sum(np.square(np.subtract(preds,train_Y)))/len(train_X[:,0])
print(mse)
return mse
x0=weights
res = minimize(loss, x0, args=(), method = 'BFGS', options={'maxiter':1,'eps':1e-6,'disp':True})
#res = minimize(loss, x0, method='L-BFGS-B', options={'disp': True, 'maxls': 1, 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 1, 'ftol': 0.5, 'maxcor': 1, 'maxfun': 1})
#res = minimize(loss, x0, args=(), method='Newton-CG', jac=None, hess=None, hessp=None, tol=None, callback=None, options={'disp': False, 'xtol': 1e-05, 'eps': 1.4901161193847656e-08, 'return_all': False, 'maxiter': 1})
inv_hess = res['hess_inv']
1) El modelo entrena extremadamente bien, pero cuando intento ejecutar el minimizador scipy para una sola iteración con los pesos previamente entrenados, me encuentro con problemas.
Salida al intentar método = BFGS:
0.458706819754
0.457811632697
0.458706716791
...
0.350124572422
0.350186770445
0.350125320636
ValueErrorTraceback (most recent call last)
---> 19 res = minimize(loss, x0, args=(), method = 'BFGS', tol=1, options={'maxiter':1,'eps':1e-6,'disp':True})#,'gtol':0.1}, tol=5)
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
442 return _minimize_cg(fun, x0, args, jac, callback, **options)
443 elif meth == 'bfgs':
--> 444 return _minimize_bfgs(fun, x0, args, jac, callback, **options)
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
963 try: # this was handled in numeric, let it remaines for more safety
--> 964 rhok = 1.0 / (numpy.dot(yk, sk))
965 except ZeroDivisionError:
966 rhok = 1000.0
ValueError: operands could not be broadcast together with shapes (7213,2048) (2048,1024)
Salida al intentar método = L-BFGS-B:
ValueErrorTraceback (most recent call last)
---> 20 res = minimize(loss, x0, method='L-BFGS-B', options={'disp': True, 'maxls': 1, 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 1, 'ftol': 0.5, 'maxcor': 1, 'maxfun': 1})
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
448 elif meth == 'l-bfgs-b':
449 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
--> 450 callback=callback, **options)
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/lbfgsb.pyc in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)
300 raise ValueError('maxls must be positive.')
301
--> 302 x = array(x0, float64)
303 f = array(0.0, float64)
304 g = zeros((n,), float64)
ValueError: setting an array element with a sequence.
Salida al intentar método = Newton-CG
ValueErrorTraceback (most recent call last)
---> 21 res = minimize(loss, x0, args=(), method='Newton-CG', jac=None, hess=None, hessp=None, tol=None, callback=None, options={'disp': False, 'xtol': 1e-05, 'eps': 1.4901161193847656e-08, 'return_all': False, 'maxiter': 1})
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
445 elif meth == 'newton-cg':
446 return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
--> 447 **options)
448 elif meth == 'l-bfgs-b':
449 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
/opt/anaconda3/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback, xtol, eps, maxiter, disp, return_all, **unknown_options)
1438 _check_unknown_options(unknown_options)
1439 if jac is None:
-> 1440 raise ValueError('Jacobian is required for Newton-CG method')
ValueError: Jacobian is required for Newton-CG method
2) La siguiente tarea es obtener la derivada de las salidas del modelo con respecto a las entradas del modelo. Por ejemplo, para un parámetro estelar (una de las salidas), digamos Temperatura, necesito encontrar las derivadas parciales con respecto a cada una de las 7213 entradas. Y luego haga lo mismo para cada una de las 3 salidas.
Básicamente, mi primera tarea (1) es encontrar una forma de devolver la matriz de Hesse inversa de mi modelo y luego (2) necesito encontrar una forma de devolver las derivadas parciales de primer orden de mis salidas con respecto a mis entradas .
¿Alguien tiene alguna idea sobre alguna de estas dos tareas? Gracias.
EDITAR
Estoy tratando de usar theano.gradient.jacobian () para devolver la matriz jacobiana de mi salida w.r.t. mis entradas He convertido mi modelo en una función de los pesos del modelo y usé esa función como primer parámetro en theano.gradient.jacobian (). Mi problema surge cuando intento ejecutar el gradiente con matrices multidimensionales en las que los pesos de mi modelo y los datos de entrada están en forma.
import theano.tensor as T
weights_in_model = T.dvector('model_weights')
x = T.dvector('x')
def pred(x,weights_in_model):
weights = T.stack((weights_in_model[0],weights_in_model[1]), axis=0)
x = T.shape_padright(x, n_ones=1)
prediction=T.dot(x, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[2],weights_in_model[3]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[4],weights_in_model[5]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[6],weights_in_model[7]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[8],weights_in_model[9]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[10],weights_in_model[11]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
prediction = T.clip(prediction, 0, 9999.)
weights = T.stack((weights_in_model[12],weights_in_model[13]), axis=0)
prediction = T.shape_padright(prediction, n_ones=1)
prediction = T.dot(prediction, weights)
T.flatten(prediction)
return prediction
f=theano.gradient.jacobian(pred(x,weights_in_model),wrt=x)
h=theano.function([x,weights_in_model],f,allow_input_downcast=True)
x = train_X
weights_in_model = model.get_weights()
h(x,weights_in_model)
Esta última línea da el error:
TypeError: ('Bad input argument to theano function with name "<ipython-input-365-a1ab256aa220>:1" at index 0(0-based)', 'Wrong number of dimensions: expected 1, got 2 with shape (2000, 7213).')
Pero cuando cambio las entradas a:
weights_in_model = T.matrix('model_weights')
x = T.matrix('x')
Me sale un error de la línea:
f=theano.gradient.jacobian(pred(x,weights_in_model),wrt=x)
leyendo:
AssertionError: tensor.jacobian expects a 1 dimensional variable as `expression`. If not use flatten to make it a vector
¿Alguna idea sobre cómo solucionar esto?