numpy druga pochodna wielowymiarowej tablicy

Question

Mam zestaw danych symulacyjnych, w których chciałbym znaleźć najniższe nachylenie w n wymiarach. Rozmieszczenie danych jest stałe wzdłuż każdego wymiaru, ale nie wszystkie takie same (mogłem to zmienić ze względu na prostotę).

Mogę żyć z pewną niedokładnością numeryczną, szczególnie w kierunku krawędzi. Wolałbym nie generować splajnu i używać tej pochodnej; tylko na surowych wartościach byłoby wystarczające.

Możliwe jest obliczenie pierwszej pochodnej za pomocą numpy przy użyciu funkcji numpy.gradient().

import numpy as np 

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

---

To uwagi odnośnie Laplace stosunku do matrycy juty; to już nie jest pytanie, ale ma pomóc zrozumieć przyszłych czytelników.

Jako funkcję testową stosuję funkcję 2D w celu określenia obszaru "najbardziej płaskiego" poniżej progu. Poniższe wyniki pokazują różnicę między wynikami przy użyciu minimum second_derivative_abs = np.abs(laplace(data)) i co najmniej następujące elementy:

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 

Skala kolorów przedstawia wartości funkcji, strzałki pokazują pierwszą pochodną (gradientu), Czerwony kropka punkt najbliższy zeru, a linia czerwona próg.

Funkcja generatora dla danych to (1-np.exp(-10*xi**2 - yi**2))/100.0 z XI, yi generowane z np.meshgrid.

Laplace'a:

Heskie:

Answer 1

Drugi pochodne są przez Hessian matrix. Oto implementacja Pythona do tablic ND, który polega na stosowaniu np.gradient dwukrotnie i przechowywania danych wyjściowych odpowiednio,

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) 

Należy pamiętać, że jeśli jesteś zainteresowany tylko w wielkości drugich pochodnych, można użyć Laplace operator realizowane przez scipy.ndimage.filters.laplace, który jest śladem (suma elementów diagonalnych) heskiej macierzy.

Pobranie najmniejszego elementu macierzy Hesji może posłużyć do oszacowania najniższego nachylenia w dowolnym kierunku przestrzennym.

Answer 2

Możesz zobaczyć hessian Matrix jako gradient gradientu, gdzie zastosujesz gradient po raz drugi dla każdego komponentu pierwszego gradientu obliczonego tutaj jest wikipedią link definiującą heską matrycę i widać wyraźnie, że jest to gradient gradientu tutaj jest realizacja pyton określenie gradient, następnie juty:

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 

w teście Heskie matrycy (x^2 + y^2) 2 * I_2 gdzie I_2 jest macierzą tożsamości wymiaru 2