derivata numpy seconda di un array ndimensionale

Question

Ho un set di dati di simulazione dove vorrei trovare la pendenza più bassa in n dimensioni. La spaziatura dei dati è costante lungo ogni dimensione, ma non tutti uguali (potrei cambiarla per semplicità).

Posso vivere con qualche imprecisione numerica, specialmente verso i bordi. Preferirei fortemente non generare una spline e usare quella derivata; solo sui valori grezzi sarebbe sufficiente.

È possibile calcolare la prima derivata con numpy utilizzando la funzione numpy.gradient().

import numpy as np 

data = np.random.rand(30,50,40,20) 
first_derivative = np.gradient(data) 
# second_derivative = ??? <--- there be kudos (: 

---

Questo è un commento per quanto riguarda Laplace rispetto alla matrice di iuta; questa non è più una domanda, ma ha lo scopo di aiutare la comprensione dei futuri lettori.

Uso come testcase una funzione 2D per determinare l'area "piatta" al di sotto di una soglia. Le seguenti immagini mostrano la differenza nei risultati tra utilizzando il minimo second_derivative_abs = np.abs(laplace(data)) e il minimo dei seguenti:

second_derivative_abs = np.zeros(data.shape) 
hess = hessian(data) 
# based on the function description; would [-1] be more appropriate? 
for i in hess[0]: # calculate a norm 
    for j in i[0]: 
     second_derivative_abs += j*j 

La scala di colore rappresenta i valori funzioni, le frecce rappresentano la derivata prima (pendenza), il punto rosso la punto più vicino a zero e la linea rossa la soglia.

La funzione di generatore per i dati era (1-np.exp(-10*xi**2 - yi**2))/100.0 con xi, yi generato con np.meshgrid.

Laplace:

dell'Assia:

Answer 1

derivate seconde sono dati dal Hessian matrix. Ecco un'implementazione Python per array ND, che consiste nell'applicare il np.gradient due volte e memorizzare l'uscita appropriato,

import numpy as np 

def hessian(x): 
    """ 
    Calculate the hessian matrix with finite differences 
    Parameters: 
     - x : ndarray 
    Returns: 
     an array of shape (x.dim, x.ndim) + x.shape 
     where the array[i, j, ...] corresponds to the second derivative x_ij 
    """ 
    x_grad = np.gradient(x) 
    hessian = np.empty((x.ndim, x.ndim) + x.shape, dtype=x.dtype) 
    for k, grad_k in enumerate(x_grad): 
     # iterate over dimensions 
     # apply gradient again to every component of the first derivative. 
     tmp_grad = np.gradient(grad_k) 
     for l, grad_kl in enumerate(tmp_grad): 
      hessian[k, l, :, :] = grad_kl 
    return hessian 

x = np.random.randn(100, 100, 100) 
hessian(x) 

nota che se si è interessati solo nella grandezza delle derivate seconde, è possibile utilizzare il Laplace operator attuato di scipy.ndimage.filters.laplace, che è la traccia (somma di elementi diagonali) della matrice hessiana.

Prendere l'elemento più piccolo della matrice hessiana potrebbe essere utilizzato per stimare la pendenza più bassa in qualsiasi direzione spaziale.

Answer 2

È possibile visualizzare la matrice di Hessian come un gradiente di gradiente, in cui si applica il gradiente una seconda volta per ciascun componente del primo gradiente calcolato qui è una wikipedia link definig matrice Hessian e si può vedere chiaramente che è un gradiente di gradiente , ecco un'implementazione pitone definente gradiente poi iuta:

import numpy as np 
#Gradient Function 
def gradient_f(x, f): 
    assert (x.shape[0] >= x.shape[1]), "the vector should be a column vector" 
    x = x.astype(float) 
    N = x.shape[0] 
    gradient = [] 
    for i in range(N): 
    eps = abs(x[i]) * np.finfo(np.float32).eps 
    xx0 = 1. * x[i] 
    f0 = f(x) 
    x[i] = x[i] + eps 
    f1 = f(x) 
    gradient.append(np.asscalar(np.array([f1 - f0]))/eps) 
    x[i] = xx0 
    return np.array(gradient).reshape(x.shape) 

#Hessian Matrix 
def hessian (x, the_func): 
    N = x.shape[0] 
    hessian = np.zeros((N,N)) 
    gd_0 = gradient_f(x, the_func) 
    eps = np.linalg.norm(gd_0) * np.finfo(np.float32).eps 
    for i in range(N): 
    xx0 = 1.*x[i] 
    x[i] = xx0 + eps 
    gd_1 = gradient_f(x, the_func) 
    hessian[:,i] = ((gd_1 - gd_0)/eps).reshape(x.shape[0]) 
    x[i] =xx0 
    return hessian 

Come test, la matrice Hessiana di (x^2 + y^2) è 2 * I_2 dove I_2 è la matrice identità di dimensione 2