Аппроксимация функции с помощью функционала scipy.optimize.minimize

Question

Пытаюсь аппроксимировать функцию с помощью полинома.

Допустим реальные данные задаются следующей функцией:

# Функция генерирует реальные значения 
def data_gen(x, alpha=1, beta=10, gama=5, per=1.75):
    data = alpha*x + beta*np.sin(per*x) + gama  + x**2 - x**3
    return data

Тогда возьмем несколько точек и составим датасет

import numpy as np

# Ось Х
dots = np.linspace(1,8,8)
# Число точек
num_points = len(dots)
# Реальные данные
data = data_gen(dots)

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


# Список степеней полинома
set_of_powers = list(range(0, num_points-3))
# Функция аппроксимирующего полинома 
def poly_func_2(w, dots=dots, set_of_powers=set_of_powers):
    features = np.vstack([dots ** i for i in set_of_powers]).transpose()
    power = features.shape[1]
    poly = features @ w.reshape([power,1])
    return poly

А также напишем функционал ошибки, который будем минимизировать

# Функция минимизации для поиска коэф-ов полинома
def minimize(w, data=data):
    return np.sum((data - poly_func_2(w))**2)/len(data)

На выходе получаем веса полинома

# Процесс минимизации
w0 = np.ones(shape=len(set_of_powers))
result = scipy.optimize.minimize(fun=minimize, x0=w0, method='BFGS')

# Итог
w0 = result.x

К сожалению, веса подобраны неверно. Ниже приведен график реальных данных, линии полинома, которая получается при подборе параметров с помощью scipy.optimize.minimize, а также линии полинома, которые получается при подоборе параметров с помощью LinearRegression.

Что мне необходимо изменить в коде или в функционале ошибки, чтобы веса подбирались корректно?

Реализация в коллабе: https://colab.research.google.com/drive/1UvTor_ewTR877A7d6I20rkx49PifLaCx?usp=sharing

Answer 1

Результат вашей функции poly_func_2 - столбец. Поэтому разность data - poly_func_2(w) в minimize - квадратная матрица. Сумма всех элементов этой матрицы совсем не похожа на сумму квадратов разностей.

Достаточно, чтобы ваша функция вычисления полинома возвращала строку, а не столбец:

def poly_func_2(w, dots=dots, set_of_powers=set_of_powers):
    features = np.vstack([dots ** i for i in set_of_powers]).transpose()
    power = features.shape[1]
    poly = features @ w.reshape([power,1])
    return poly.T

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

def poly_func_2(w, dots=dots, set_of_powers=set_of_powers):
    return np.polyval(np.flip(w), dots)

Это даёт более чем 10% прирост скорости.

Функция numpy.polyval предполагает, что коэффициенты идут в порядке от старшей степени к младшей, а у вас наоборот, поэтому перед вычислением полинома нужно поменять порядок на обратный вызовом numpy.flip.

А второе, что хочется сделать - заменить вызов minimize на curve_fit - специализированную функцию аппроксимации.

def model(x, *w):
    return np.polyval(w, x)
w, wcov = spo.curve_fit(model, dots, data, p0=np.ones(shape=num_points-3))
w_error = np.sqrt(np.diag(wcov))

Функция возвращает оптимальный набор параметров и соответствующую ковариационную матрицу. Квадратный корень из диагональных элементов даёт оценку ошибки аппроксимации (w_error в примере).

Вариант с curve_fit работает в 40 раз быстрее, чем ваш вариант с minimize.

БЛОГ НА HUSL

Аппроксимация функции с помощью функционала scipy.optimize.minimize

Ответы (1 шт):