Метод градиентного спуска для n переменных

[vsrmis]

Столкнулся я тут с проблемой. Есть набор значений функции x[n] и y[n];
Функция одной переменной и n коэффициентов. Предположим:
$a+b\cdot (1 - e^{c\cdot x})$

И мне неизвестны коэффициенты a,b,c;
Я заглянул в гугл и подумал, что эту задачу лучше всего решать методом градиентного спуска.

Я понимаю так этот алгоритм:
Есть x[n],y[n] и F(x,args[m]), которая, предположительно, описывает зависимость y[n] от x[n].
Шаг 1. - устанавливаем начальные значения args[m]. Вычисляем новый Y, назовём его Ych
Шаг 2. - вычисляем сумму квадратов невязок (Y[i]-Ych[i])^2. если она не равна 0 то вычисляем градиент и вектор градиента.
Шаг 3.Если модуль градиента не равен нулю. Изменяем каждый аргумент на значение вектора градиента.
Шаг 4. переходим к шагу 2

Наваял я такое изобилие

Используется синтаксис C++

#include <math.h>
//#include "Solover.h"
#include <iostream>

using namespace std;

struct GARDIENT {
        double* vector;//массив вектора гардиента
        int len;// длинна вектора
        double Norm;//значение гардиента
};

double deltaX = 1e-6; //шаг для дифиринцирования 
int len = 10; // длинна массива исходных точек 
int lenargs = 3;// длинна массива аргументов
double *x = new double[len]; //ИМХО массив иксов
double *y = new double[len]; //ИМХО массив Y
double *args = new double[lenargs]; // массив аргументов
double *Ych;// массив апроксимированных значений


double Funk(double x,double* args){
        // тестовая функция
        double a = args[0];
        double b = args[1];
        double c = args[2];     
        return a+b*(1-exp(c*x));
}
                
void countYch(){        
        // функция считает Y с текущими коеффициентами
                Ych = new double[len];
                for(int i =0;i<len;i++)
                        Ych[i]=Funk(x[i],args);
}
        
double SummAbsValsResiduals(){
// вычисление суммы квадратов невязок
                double Summ = 0;
                countYch();
                for (int i =0;i<len;i++){       
                        Summ+=pow((y[i]-Ych[i]),2);
                }
                return Summ;
}
        
double Derivative(int DerivativeVar){
        //производная
                args[DerivativeVar]-=deltaX;
                double Res1 = SummAbsValsResiduals();
                args[DerivativeVar]+=2*deltaX;
                double Res2 = SummAbsValsResiduals();
                double deltaY = Res2-Res1;
                args[DerivativeVar]-=deltaX;
                return deltaY/(2*deltaX);
}
        
        
GARDIENT Gardient(){
        //гардиент
                GARDIENT result;
                result.Norm = 0;
                result.len = lenargs;
                result.vector = new double[lenargs];
                for(int i =0;i<result.len ;i++){
                        double derivative = Derivative(i);
                        result.Norm = result.Norm+derivative*derivative;
                        result.vector[i]= derivative*deltaX;
                }
                result.Norm = sqrt(result.Norm);
                return result;
}


int main(){
        double a = 1;
        double b = 2;
        double c = 2;
        // исходные коэффициенты
        args[0]=a;
        args[1]=b;
        args[2]=c;
        for(int i = 0;i<len;i++){
                x[i]=i;
                y[i] = Funk(x[i],args);
        }
        
         //коефициенты поменялись
        args[0]=1;
        args[0]=2;
        args[0]=1;
        //их надо восстановить
        countYch();
        GARDIENT gar = Gardient();
        int n = 0;
        while((gar.Norm>0.1)&&(n<1e6)){
                // пока мы не пришли к локальному минимуму функции
                n++;
                countYch();
                gar = Gardient();
                for(int i =0;i<gar.len;i++){
                        args[i]=args[i]-gar.vector[i];
                }
                cout<<"----"<<gar.Norm<<endl;
                for(int i = 0;i<gar.len;i++){
                                cout<<"---->"<<gar.vector[i]<<endl;//значение вектор гардиента для каждой переменной 
                }
                cout<<"SUMM : "<<SummAbsValsResiduals()<<endl;// сумма на текущем этапе
        }
        return 1;
}

 

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

Производную я вычисляю вычисляя значения функции на отрезке + dx и -dx , и делю эти (dy1+dy2)/2dx

И порой, даже когда я в локальном минимуме, этот метод всё-равно мне даёт очень не хорошие значения , которые порой, становтся критическими для работы программы

Я подозреваю примерно в таких случаях:

как с этим бороться? Как правильно минимизировать сумму квадратом невязок?

[ewert]

vsrmis в сообщении #642761 писал(а):

Изменяем каждый аргумент на значение вектора градиента.

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

Формально на каждом шагу надо минимизировать по $t$ функцию $F(\vec a-t\,\nabla F)$, где $\vec a$ -- вектор текущих значений параметров, после чего брать $\vec a-t\,\nabla F$ в качестве нового приближения для этих параметров. Но точная минимизация даже в одномерном случае -- задача довольно нетривиальная, и обойти её можно, скажем, так.

Если бы целевая функция была чисто квадратичной, то было бы $F(\vec a-t\,\nabla F)=F(\vec a)-t\,|\nabla F|^2+\gamma\,t^2$, где $\gamma$ -- половина второй производной по направлению градиента. Искать вторую производную -- некоторая морока, да и всё равно в неквадратичном случае к точной минимизации это не приведёт. Однако можно заметить, что для чисто квадратичного случая (к которому мы будем приближаться по мере приближения к точке минимума) для оптимального $t$ будет выполняться равенство $F(\vec a-t\,\nabla F)-\left(F(\vec a)-\frac12t\,|\nabla F|^2\right)=0$ (таковы уж свойства параболы), причём слева от точки минимума эта разность будет отрицательной, справа -- положительной. Так вот, надо ввести некий фиксированный настроечный параметр $\theta$ (скажем, $\theta=2$) и на каждом шаге делать следующее: сначала увеличивать текущее значение $t$ в $\theta$ раз до тех пор, пока та разность не станет положительной, затем -- уменьшать во столько же раз до тех пор, пока разность не станет отрицательной. Последнее значение $\vec a-t\,\nabla F$ и брать в качестве нового приближения.

[vsrmis]

Цитата:

Во-первых, не градиента, а антиградиента.

Прошу прощения, оговорился. Теперь ближе к теме.
До этого момента, я брал за t deltaX. То есть от вектора аргументов я отнимал вектор градиента помноженный на шаг по иксам... и это не работало.

Цитата:

Формально на каждом шагу надо минимизировать по функцию по t

Совершенно не понял, что вы имеете в виду. На каждом шагу я минимизирую функцию сразу по всем значениям её аргументов. Главная идея, которую я подхватил из вашего поста - изменяемое t, в зависимости от приближения к точке минимума. Вы говорите про квадратичную функцию и очень любопытные её свойства. Но у меня функция, я даже затрудняюсь сказать какая она. Это сумма квадратов невязок идеального значения и значения вычесленого с помощью произвольной функции. Я поробую ваш метод вариирования t относительно нуля с помощью дополнительного параметра, но я не понимаю каким брать начальное t?

[AleksandrPavlovich]

Советы
1. Шаги для дифференцирования берите по-меньше, но так, чтобы погрешность округления была не больше (меньше) погрешности апроксимации.

2. Можно настраивать размер t шага по антиградиенту адаптивно:
Если шаг привёл к уменьшению функции ($F(x_i+t\cdot g_i)<F(x_i)$),
ТО $x_{i+1}:=x_i+t\cdot g_i; t:=1.1\cdot t;$
Иначе $t:=t/2$ пока $F(x_i+t\cdot g_i)>F(x_i)$

Такая рекомендация есть в книге
Васильев Ф.П. Методы оптимизации (насколько я помню, именно в этой)

[Lyoha]

vsrmis в сообщении #642761 писал(а):

Шаг 3.Если модуль градиента не равен нулю. Изменяем каждый аргумент на значение вектора градиента.

Как уже заметили, минус градиента, но проблема в другом. Для использования градиентных методов градиент должен удовлетворять условию Липшица $||\nabla f(x)-\nabla f(y)||<L||x-y||$, а шаг выбираться с условием
$x_{n+1}=x_n - \gamma_n\nabla f(x_n)$, где $\gamma_n<2/L$.

[ewert]

Lyoha в сообщении #643535 писал(а):

Для использования градиентных методов градиент должен удовлетворять условию Липшица

Да функция попросту должна быть дважды гладкой, иначе метод градиентного спуска превращается в филькину грамоту ("в случае общего положения"). От липшицевой же константы нет никакого практического проку: глобально она заведомо завышена, к тому же в боевых условиях заведомо неизвестна, локально же (т.е. после достаточно близкого приближения к искомой точке) определяется собственными числами матрицы Гессе для целевой функции, которые тоже в здравом рассудке никому даже и в голову не придёт находить.

-- Пн ноя 12, 2012 19:57:09 --

Апдейт. Нет, матрица Гессе полезна, конечно. Но только сама по себе, а вовсе не её спектр. Полезна просто потому, что на методе градиентного спуска свет клином не сошёлся, и полезно перемежать его время от времени методом Ньютона.

[vsrmis]

Цитата:

Такая рекомендация есть в книге
Васильев Ф.П. Методы оптимизации (насколько я помню, именно в этой)

Огромное спасибо за эту информацию градиент действительно стал сходиться.
Но теперь стоит вторая проблема - автоматический приблизительный подбор коэфициентов.

Потому, что если начальный подбор не верен, то метод градиентного спуска ведёт меня в локальный минимум, который вовсе не является глобальным.
Я пробовал простым перебором параметров с неким шагом ( 0.7)
Но уже от -50 до 50 это работает достаточно долго.
И к тому-же, при функции
$a+b\cdot e^{x\cdot c}$
с начальными параметрами 0.4 ; 1.1 ; 2.6; оно совсем не находит приближение а как самый близкий мне выдаёт 49.4 ; 2.5 ; 2.5 С суммой квадратов не вязок = 1.47696e+18;

Как иначе найти некое приближение к глобальному минимуму?

[vsrmis]

Цитата:

Такая рекомендация есть в книге
Васильев Ф.П. Методы оптимизации (насколько я помню, именно в этой)

Огромное спасибо за эту информацию градиент действительно стал сходиться.
Но теперь стоит вторая проблема - автоматический приблизительный подбор коэфициентов.

Потому, что если начальный подбор не верен, то метод градиентного спуска ведёт меня в локальный минимум, который вовсе не является глобальным.
Я пробовал простым перебором параметров с неким шагом ( 0.7)
Но уже от -50 до 50 это работает достаточно долго.
И к тому-же, при функции
$a+b\cdot e^{x\cdot c}$
с начальными параметрами 0.4 ; 1.1 ; 2.6; оно совсем не находит приближение а как самый близкий мне выдаёт 49.4 ; 2.5 ; 2.5 С суммой квадратов не вязок = 1.47696e+18;

Как иначе найти некое приближение к глобальному минимуму?

[vsrmis]

Пробовал только что ещё метод Монте карло.
Но он совсем плохо работает. Может у меня в нём ошибка

Используется синтаксис C++

void MonteCarlo(){
        double idealsumm = 1e+4;// Приближусь хотя-бы к такой сумме квадратов невязок
        MinSumm = SummAbsValsResiduals(); // считаю начальную сумму квадратов невязок
        
        while(MinSumm>idealsumm){ // Пока текущая сумма больше к той к которой я стремлюсь буду крутится в цикле 
                for(int i = 0 ;i<lenargs;i++){ 
                        args[i]=downLimit+rand()%((int)(uperLimit-downLimit));// Все аргументы принимают случайное значение
                }
                double summ = SummAbsValsResiduals(); 
                        if(MinSumm>summ){ // если уменьшизи сумму - запомним.
                                
                                for(int i = 0 ;i<lenargs;i++){
                                        MinArgs[i]=args[i];
                                }
                                MinSumm = summ;
                        }
        }
}
int main(){
 

[ewert]

vsrmis в сообщении #644264 писал(а):

метод градиентного спуска ведёт меня в локальный минимум, который вовсе не является глобальным.

Никакой регулярный метод не способен выделить глобальный минимум из множества (если оно есть) локальных.

[vsrmis]

Цитата:

Никакой регулярный метод не способен выделить глобальный минимум из множества (если оно есть) локальных.

Вы хотите сказать, что задача хотя бы приближения к глобальному минимуму не решаема?
Ведь если я приближусь к нему так, что от остальных локальных минимумов менч будут отделать точки максимума, метод Градинетов приведёт меня к этому минимуму.

А в данной задаче глобальный минимум абсолютно точно существует. Это же сумма квадратов не вязок.

[Lyoha]

ewert в сообщении #643668 писал(а):

Lyoha в сообщении #643535 писал(а):

Для использования градиентных методов градиент должен удовлетворять условию Липшица

Да функция попросту должна быть дважды гладкой, иначе метод градиентного спуска превращается в филькину грамоту ("в случае общего положения"). От липшицевой же константы нет никакого практического проку: глобально она заведомо завышена, к тому же в боевых условиях заведомо неизвестна, локально же (т.е. после достаточно близкого приближения к искомой точке) определяется собственными числами матрицы Гессе для целевой функции, которые тоже в здравом рассудке никому даже и в голову не придёт находить.

Вторая производная нас может интересовать только как способ оценки той самой константы Липшица. Попросту дважды гладкая нам ничем не поможет, потому как надо выбрать шаг градиентного метода.

-- 14 ноя 2012, 11:51 --

vsrmis в сообщении #644265 писал(а):

Как иначе найти некое приближение к глобальному минимуму?

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

[Yuri Gendelman]

vsrmis в сообщении #644359 писал(а):

Цитата:

Никакой регулярный метод не способен выделить глобальный минимум из множества (если оно есть) локальных.

Вы хотите сказать, что задача хотя бы приближения к глобальному минимуму не решаема?
Ведь если я приближусь к нему так, что от остальных локальных минимумов менч будут отделать точки максимума, метод Градинетов приведёт меня к этому минимуму.

Вот именно для такого приближения и нет регулярных методов.
Обычно используемые методы на основе вычисляемых или приближенно оцениваемых (как в методе Нелдера-Мида) производных - принципиально локальные. Но даже для поиска локального минимума возможны ситуации (типа "банан Розенброка"), когда алгоритм ошибается.
Для поиска глобального минимума используются эвристики:
- случайный поиск
- повторение поиска из существенно разных начальных точек.

Еще надо учесть, что задача поиска глобального минимума неустойчива (некорректна?) - при малом изменении параметров целевой функции положение глобального минимума может скачком измениться (перескочить к другому лок.минимуму).

[AleksandrPavlovich]

vsrmis в сообщении #644265 писал(а):

Как иначе найти некое приближение к глобальному минимуму?

Как вариант - запустить несколько раз из разных начальных точек.
Варианты выбора начальных точек:
1. Из сетки
2. Случайным образом
3=1+2, т.е. Комбинировать 1 и 3

[vjung]

исходная функция
$f(x) = a + b(1-e^{-cx})$

эквивалентна такой:
$f(x) = \hat{a} - e^{\hat{b}x + \hat{c}}$,

где $\hat{a} = a+b, \hat{b} = -c, \hat{c} = \log(b)$. Отсюда следует, что:
$\log(y_i - \hat{a}) = \hat{b}x_i + \hat{c} $.

Таким образом, на надо найти такую $\hat{a}$, что набор точек $(\log(y_i - \hat{a}), x_i)$ даёт наименьшую ошибку при аппроксимации прямой. Что-то мне подсказывает, что в этой формулировке минимум один и задача должна решаться бисекцией очень эффективно.