Блог

Я работаю над реализацией числового трехдиагонального гессиана в Python. Для этого я использую метод конечных разностей, это то, что я пробовал до сих пор:

 import numpy as np

def principal_diagonal(x, fun):
    h = 10e-6
    result = list()
    for i in range(len(x)):
        
        temp_forward = np.array(x)
        temp_backward = np.array(x)
        temp_forward[i] = temp_forward[i]   h
        temp_backward[i] = temp_backward[i] - h
        value = (fun(temp_forward) - 2 * fun(x)   fun(temp_backward)) / (h ** 2)
        result.append(value)
    return(np.array(result))

def one_above_diagonal(x, fun):
    h =  10e-6 
    result = list()
    for i in range(len(x)-1):
        clean = np.array(x)
        forward_first = np.array(x)
        forward_first[i] = forward_first[i]   h
        forward_second = np.array(forward_first)
        forward_second[i 1] = forward_second[i 1]   h
        clean[i 1] = clean[i 1]   h
        value =  ((fun(forward_second) - fun(clean)) - (fun(forward_first) - fun(x))) / (h*h)
        result.append(value)
    return(np.array(result))

def tridiag_hessian(x, fun):
    new_array = np.ndarray((len(x), len(x)))
    di = np.diag_indices(len(x))
    new_array[di] = principal_diagonal(x, fun)
    new_array[(np.arange(len(x)-1), np.arange(1,len(x)))] = one_above_diagonal(x, fun)
    new_array[(np.arange(1,len(x)), np.arange(len(x)-1))] = one_above_diagonal(x, fun)
    return(new_array)
 

Я протестировал свои функции с помощью функции Розенброка, и мои результаты были следующими:

 from scipy.optimize import rosen
tridiag_hessian(np.array([1, 2], "float64"), rosen)

array([[ 401.99964246, -400.00216472],
       [-400.00216472,  200.00001655]])
 

Результаты довольно близки к фактическому решению, что вполне прилично, но для более сложных задач я хотел бы иметь большую точность. Есть ли в моем коде что-нибудь, что можно было бы скорректировать? Или, может быть, альтернативный метод? Я попытался использовать numdifftools , а затем приступил к «тридиагонализации» результатов этой функции, но это было невероятно медленно, когда я попробовал этот подход в задаче, в которой использовалась функция с более чем 40 переменными. Я также рассматривал возможность использования np.sqrt(np.finfo(float).eps) вместо 0,000001, но я не получил никакой точности.