Как стать автором
Обновить

Градиентный спуск, который мы заслужили

Уровень сложности Средний
Время на прочтение 5 мин
Количество просмотров 817
Python *Машинное обучение *
Recovery mode
Из песочницы

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

Данную статью побудило написать качество современных публикаций на тему проблем современного машинного обучения (https://habr.com/ru/company/skillfactory/blog/717360/ Wouter van Heeswijk, PhD).

Что в этом не так?

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

В любой формулировке мы или исследуем все пространство малыми шагами (и не пропускаем оптимум), или приходим к оптимуму большими шагами. Вообще правильнее было бы настраивать размер шагов исходя из какого-либо анализа функции ошибки с использованием БПФ, но это в следующей статье.

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

Залог оптимальности эвристики в отсутствии занижения функции полезности. В случае попытки сэкономить миллион\миллиард шагов алгоритма можно получить субоптимальные эвристики и не получить решение и за триллион.

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

Что делать?

Давайте напишем свою версию градиентного спуска, с множеством начальных точек и очередью с приоритетом.

В общем алгоритм выглядит так:

  1. Инициализируем N начальных точек, для каждой из точек считаем ошибку, записываем все в очередь с приоритетом (можно делать и параллельно)

    Цикл:

  2. Достаем точку из очереди, считаем градиент ошибки, считаем изменение положения, считаем ошибку

  3. Если точка застряла в локальном оптимуме (последние N положений отличаются менее чем на эпсилон) - добавляем точку в список застрявших в локальном оптимуме (возможно, инициализируем новую точку)

  4. Если нет - помещаем точку обратно в очередь с приоритетом

Повторяем п. 2-4 пока не закончатся шаги или не получим требуемый уровень ошибки

Что нам это даст?

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

Сложность (количество шагов, требуемое для решения) будет меньше, чем N градиентных спусков (по M шагов), или M*N, но больше M+N(инициализация и спуск).

Код

python
from queue import PriorityQueue
import numpy as np
import random
from math import sqrt

##Parameters

#learning rate
dt = 0.01

#maximum step count
maxSteps = 1000000

#minimum_error
minError = 0.01
##Specifical parameters

#small number to calculate derivative of target function (for using script with large object functions, like neural nets or so)
delta = 0.001

#target vector len
trace_len = 20

#initilal point count
N = 100000

#size of parameters space
spaceSize=[[-1000,1000],[-1000,1000]]

#for analysis purpose
GetErrorCallCount = 0

localOptimums = []

#taked from https://github.com/behnamasadi/gradient_descent/
def objective_function(vector):
    z=-(4*np.exp(-(vector[0]-4)**2 - (vector[1]-4)**2)+2*np.exp(-(vector[0]-2)**2 - (vector[1]-2)**2))+4
    return z
    
#we will compute exact derivative
def get_objective_function_derivatives(vector,delta,error):
    res = []
    
    initial_error = error
    
    pos_vector = vector[0]

    for n in range(len(pos_vector)):
        vec = pos_vector
        vec[n] = vec[n]+delta
        
        derivative_n = -(GetError(vec)-initial_error)/delta
        res.append(derivative_n)

    
    #change to correct derivative calculation
    return res
    
#returns normalized derivative vector
def get_error_gradient(vector,delta, error):
    derivative = get_objective_function_derivatives(vector, delta, error)
    #recalculate derivative, if take 0 multiply delta by ten
    
    while all(d ==0 for d in derivative):
        # print("get zero derivative")
        # print(derivative)
        # print(delta)
        delta=delta*10
        
        #doesnt help, return ReInitialization flag
        if (delta>1000):
            return derivative
            
            
        derivative = get_objective_function_derivatives(vector, delta, error)
        
    return derivative
    
def initGradDesc(sampleVector, priorityQueue, N):
    res = []
    global spaceSize
    #Let's generate N random vectors with same length
    for num in range(N):

        #Initialize derivatives and vector randomly
        pos = []
        for i in range(len(sampleVector)):
            pos.append(random.uniform(spaceSize[i][0],spaceSize[i][1]))
        
        error = GetError(pos)

        res.append([error, [pos]])
    
    #Now we have N random vectors, please note that we must have normaziled data everywhere )
    #lets add all those vectors in priorityQueue with our GetError Function

    for x in res:
        priorityQueue.put(x)

#TODO: rewrite euler integration and calculation of derivatives

#use simpliest Euler method to find change of coordinates with some derivatives
def calculateNewPosition(stateVector, error_grad):
    delta = []

    delta = [x[0]+x[1]*dt for x in zip(stateVector[0],error_grad)]

    return delta

    
#returns error for current vector as a solution
def GetError(vector):
    #update call count
    global GetErrorCallCount
    GetErrorCallCount = GetErrorCallCount+1;
    return objective_function(vector)
    
def getLengthBetween(item, newPoint):
    l = sqrt(sum((i[0]-i[1])*(i[0]-i[1]) for i in zip(item,newPoint)))
    return l
    
def makeGradDesc(priorityQ, point):
    global delta
    
    errorInit = point[0]
    state_vector = point[1]
 
    #Note that point structure is [vec, derivatives_1,...,derivatives_N]
    error_grad = get_error_gradient(state_vector,delta,errorInit)
   
    newPoint = calculateNewPosition(state_vector, error_grad)
    
    error = GetError(newPoint)
    
    #merge data into one list
    res = [newPoint]

    curr_trace = 0;
    for positions in state_vector:
        res.append(positions)
        curr_trace = curr_trace+1
        if (curr_trace>=trace_len):
            break
        
    #calculate error 
    res = [error,res]
        
    #Check for local optimum (cost function dont change last N steps), initialize new point
    if (all(getLengthBetween(item, newPoint) < dt for item in state_vector)):
        localOptimums.append(res)
        initGradDesc(newPoint, priorityQ, 1)
    else:
        priorityQ.put(res)
    
#main gradient descent function that runs for MaxStepCount or while minimal error is not reached    
def PreformDradDesc(priorityQ, MaxStepCount,minError):
    point = priorityQ.get()
    error = point[0]
    steps = 0
    while ((abs(error)>=abs(minError)) and (steps<MaxStepCount)):
        steps = steps + 1
        makeGradDesc(priorityQ, point)
        point = priorityQ.get()
        error =  point[0]
        
    return [steps, error, point]

#Where to store statistics about steps count/result error
research_results = []

#Queue needed for our A*
priorityQ = PriorityQueue()

initGradDesc(spaceSize, priorityQ, N+1)

print("Data example: ")
point = priorityQ.get()
print(point)


#Print results
print("Cost function calls for initialization: ")
print(GetErrorCallCount)

minimum_error = PreformDradDesc(priorityQ, maxSteps, minError)
print("Minimum error: ")
print(minimum_error[1])
print("Result point: ")
print(minimum_error[2][1][1])

print("Initial points: ")
print(N)
print("Total steps: ")
print(minimum_error[0])
print("Maximum steps: ")
print(maxSteps)

print("Cost function calls: ")
print(GetErrorCallCount)

print("Local minimums queue len: ")
print(len(localOptimums))
# for item in localOptimums:
    # print(item) 
    
print("Raw results: ")
for item in research_results:
    print(item)

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

Результаты анализа

В 10 запусках из 10 был получен глобальный минимум (что уже само по себе неплохо :) ).

Детектор локальных минимумов работает как на локальный минимум, так и на плато.

Использованные материалы

  1. картинка

  2. целевая функция

Теги:
Хабы:
Всего голосов 7: ↑3 и ↓4 -1
Комментарии 2
Комментарии Комментарии 2

Публикации

Истории

Работа

Python разработчик
144 вакансии
Data Scientist
122 вакансии