Есть ли «стандартный» способ вычисления численного градиента?

Question

Я пытаюсь вычислить численный градиент гладкой функции в C++. И значение параметра может меняться от нуля до очень большого числа (возможно, от 1 до 10 ... 1 020)

Я использовал функцию f (x, y) = 10 * x^3 + y^3 в качестве тестового стенда, но я что если x или y слишком велико, я не могу получить правильный градиент.

Вот мой код, чтобы вычислить graidient:

#include <iostream> 
#include <cmath> 
#include <cassert> 
using namespace std; 
double f(double x, double y) 
{ 
    // black box expensive function 
    return 10 * pow(x, 3) + pow(y, 3); 
} 
int main() 
{ 
    // double x = -5897182590.8347721; 
    // double y = 269857217.0017581; 
    double x = 1.13041e+19; 
    double y = -5.49756e+14; 
    const double epsi = 1e-4; 

    double f1 = f(x, y); 
    double f2 = f(x, y+epsi); 
    double f3 = f(x, y-epsi); 
    cout << f1 << endl; 
    cout << f2 << endl; 
    cout << f3 << endl; 
    cout << f1 - f2 << endl; // 0 
    cout << f2 - f3 << endl; // 0 
    return 0; 
} 

Если я использую приведенный выше код для вычисления градиента, градиент будет равен нулю!

Функция testbench, 10 * x^3 + y^3, является просто демо, реальная проблема, которую мне нужно решить, на самом деле является функцией черного ящика.

Итак, есть ли «стандартный» способ вычисления численного градиента?

Answer 1

только способ расчета градиента исчисления.

Градиент представляет собой вектор:

g(x, y) = Df/Dx i + Df/Dy j 

где (I, J) являются единичными векторами в направлениях х и у, соответственно.

Один из способов приближенного производных первых разностей порядка:

Df/Dx ~ (f(x2, y)-f(x1, y))/(x2-x1) 

и

Df/Dy ~ (f(x, y2)-f(x, y1))/(y2-y1) 

Это не похоже на то, что вы делаете.

Вы имеете замкнутую форму выражения:

g(x, y) = 30*x^2 i + 3*y^2 j 

можно подключить значения (х, у), и вычисления градиента точности в любой точке. Сравните это с вашими различиями и посмотрите, насколько хорошо ваше приближение.

Как вы его воплощаете в цифрах, это ваша ответственность. (10^19)^3 = 10^57, правильно?

Какой размер двойной на вашей машине? Это 64-битное число с плавающей запятой двойной точности IEEE?

Answer 2

Использование

dx = (1+abs(x))*eps, dfdx = (f(x+dx,y) - f(x,y))/dx 
dy = (1+abs(y))*eps, dfdy = (f(x,y+dy) - f(x,y))/dy 

, чтобы получить значимые размеры шагов для больших аргументов.

Использовать eps = 1e-8 для односторонних разностных формул, eps = 1e-5 для центральных коэффициентов разности.

Исследуйте автоматическую дифференциацию (см. Autodiff.org) для производных без разностных коэффициентов и, следовательно, намного меньше числовых ошибок.

Answer 3

Вам необходимо учитывать необходимую точность.

На первый взгляд, так как |y| = 5.49756e14 и epsi = 1e-4, что вам нужно, по крайней мере ⌈log2(5.49756e14)-log2(1e-4)⌉ = 63 биты мантиссы точности (то есть количество битов, используемых для кодирования цифр номера, также известный как мантиссы) для y и y+epsi следует рассматривать другой.

Формат с плавающей запятой с двойной точностью имеет только 53 бит значимой точности (при условии, что это 8 байтов). Таким образом, в настоящее время f1, f2 и f3 являются точно такими же, поскольку y, y+epsi и y-epsi равны.

Теперь рассмотрим предел: y = 1e20, и результат вашей функции, 10x^3 + y^3. Давайте теперь проигнорируем x, поэтому давайте возьмем f = y^3. Теперь мы можем рассчитать точность, необходимую для f(y) и f(y+epsi), чтобы быть разными: f(y) = 1e60 и f(epsi) = 1e-12. Это дает минимальную значимую точность ⌈log2(1e60)-log2(1e-12)⌉ = 240 бит.

Даже если вы должны были использовать long double типа, при условии, что 16 байт, ваши результаты не будут отличаться: f1, f2 и f3 все равно будет равно, даже если y и y+epsi не будет.

Если мы учтем x, максимальное значение f будет 11e60 (с x = y = 1e20). Таким образом, верхний предел точности равен ⌈log2(11e60)-log2(1e-12)⌉ = 243 битам или не менее 31 байта.

Один из способов решения вашей проблемы - использовать другой тип, возможно, bignum, используемый как фиксированная точка.

Другой способ - переосмыслить вашу проблему и разобраться с ней по-разному. В конечном счете, вы хотите f1 - f2. Вы можете попробовать разложить f(y+epsi). Опять же, если вы игнорируете x, f(y+epsi) = (y+epsi)^3 = y^3 + 3*y^2*epsi + 3*y*epsi^2 + epsi^3. Итак, f(y+epsi) - f(y) = 3*y^2*epsi + 3*y*epsi^2 + epsi^3.

Answer 4

Во-первых, вы должны использовать центральную разностную схему, которая является более точной (путем отмены еще одного срока развития Тейлора).

(f(x + h) - f(x - h))/2h 

вместо

(f(x + h) - f(x))/h 

Тогда выбор h имеет решающее значение и использование фиксированной константой является худшее, что вы можете сделать. Поскольку для маленьких x, h будет слишком большим, чтобы формула приближения больше не работала, а для больших x, h будет слишком маленькой, что приведет к серьезной ошибке усечения.

Гораздо лучший выбор - принять относительное значение, h = x√ε, где ε - это машина epsilon (1 ulp), которая дает хороший компромисс.

(f(x(1 + √ε)) - f(x(1 - √ε)))/2x√ε 

берегись, когда x = 0, относительное значение не может работать, и вы должны падать обратно на константу. Но тогда ничто не говорит вам, что использовать!

Answer 5

Мы можем исследовать поведение ошибки в производной с использованием следующей программы - она ​​вычисляет одностороннюю производную и центральную производную, основанную на разности, с использованием разного размера шага. Здесь я использую x и y ~ 10^10, что меньше, чем вы использовали, но должны проиллюстрировать одну и ту же точку.

#include <iostream> 
#include <cmath> 
#include <cassert> 
using namespace std; 
double f(double x, double y) { 
    return 10 * pow(x, 3) + pow(y, 3); 
} 

double f_x(double x, double y) { 
    return 3 * 10 * pow(x,2); 
} 

double f_y(double x, double y) { 
    return 3 * pow(y,2); 
} 

int main() 
{ 
    // double x = -5897182590.8347721; 
    // double y = 269857217.0017581; 
    double x = 1.13041e+10; 
    double y = -5.49756e+10; 
    //double x = 10.1; 
    //double y = -5.2; 

    double epsi = 1e8; 
    for(int i=0; i<60; ++i) { 
    double dfx_n = (f(x+epsi,y) - f(x,y))/epsi; 
    double dfx_cd = (f(x+epsi,y) - f(x-epsi,y))/(2*epsi); 
    double dfx = f_x(x,y); 
    cout<<epsi<<" "<<fabs(dfx-dfx_n)<<" "<<fabs(dfx - dfx_cd)<<std::endl; 
    epsi/=1.5; 
    } 
    return 0; 
} 

Выходные данные показывают, что 1-стороннее разница получает нас оптимальную ошибку около 1.37034e+13 на этапе длиной около 100,0. Обратите внимание, что в то время как эта ошибка выглядит большой, как относительная погрешность это 3.5746632302764072e-09 (так как точное значение 3.833e+21)

В сравнении 2-сторонняя разница получает оптимальную ошибку около 1.89493e+10 с размером шага около 45109.3. Это на три порядка лучше (с гораздо большим размером шага).

Как мы можем определить размер шага? Ссылка на комментарии Ив Даостса дает нам приблизительную оценку:

h=x_c sqrt(eps) для односторонней связи и h=x_c cbrt(eps) для 2-сторонней.

Но в любом случае, если требуемый размер шага для достойной точности при x ~ 10^10 равен 100,0, размер шага с x ~ 10^20 также будет на 10 10 больше. Таким образом, проблема заключается только в том, что ваш размер шага способ слишком мал.

Это можно проверить, увеличив начальный шаг в приведенном выше коде и сбросив значения x/y до исходных значений.

Тогда как ожидаются производная O(1e39), лучше 1-стороннее ошибка около O(1e31) происходит вблизи шаг длиной 5.9e10, лучше 2-сторонняя ошибка о O(1e29) происходит вблизи шаг длиной 6.1e13.

Answer 6

Поскольку числовое дифференцирование плохо обусловлено (что означает, что небольшая ошибка может существенно повлиять на ваш результат), вы должны рассмотреть возможность использования Cauchy's integral formula. Таким образом, вы можете вычислить n-ю производную с интегралом. Это приведет к меньшим проблемам с учетом точности и стабильности.