29.10.2025 21:57
lapkin25 4 880 Источник

Поставим задачу автоматического подбора весовых коэффициентов вместе с настройкой функции связи обобщенной линейной модели.

Регрессионная модель:

y=f(w_1x_1+\ldots+w_dx_d)

Неизвестными считаются не только весовые коэффициенты w_1,\ldots,w_d , но и функция f .

Метод

Как выбрать веса одновременно с функцией связи? Если перейти к одной координате z=w_1x_1+\ldots+w_dx_d , потребуется провести через точки (z_i,y_i) кривую, принимающую вблизи каждого z очень близкие значения.

Для нахождения функции f используем кубические сплайны класса C^2. Известно, что (при соответствующем выборе краевых условий) кубические сплайны минимизируют функционал

J(f)=\int_a^b (f''(x))^2\,dx

Значение функционала J на кубическом сплайне мы возьмем в качестве меры качества подбора весовых коэффициентов. Однако, поскольку при растяжении точек \mathbf x_i в C раз функционал уменьшается в C^3 раз, наша целевая функция будет такой:

J_0(\mathbf w)=(z_{\max}-z_{\min})^3\int_{z_{\min}}^{z_{\max}} (f''(z))^2\,dz

В этой формуле функция f - это кубический сплайн, построенный по точкам (z_i, y_i).

Эксперименты показали, что для оптимизации целевой функции J_0 не подходят стандартные алгоритмы BFGS и Nelder-Mead. В связи с этим будем использовать генетический алгоритм, реализованный в библиотеке PyGAD.

Пример

Сгенерируем в 3-мерном пространстве набор случайных точек \mathbf x_i. Зададим вектор w = (1, 2, 3) и вычислим y_i = f(\mathbf w \cdot \mathbf x_i)+\varepsilon, где f(z) = \sqrt{z}, \varepsilon - случайный шум. Посмотрим, сможет ли метод восстановить веса w вместе с функцией f.

Код
import numpy as np
from scipy.optimize import minimize
from scipy.interpolate import CubicSpline
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import pygad
import math
import random


def J(z, y):
    ind = np.argsort(z)
    cs = CubicSpline(z[ind], y[ind])
    d2 = cs.derivative(nu=2)
    integrand = lambda x: d2(x) ** 2
    I = sum([integrate.quad(integrand, d2.x[i], d2.x[i + 1])[0] for i in range(len(d2.x) - 1)])
    return I * (d2.x[-1] - d2.x[0]) ** 3


def train_spline(x, y):
    N, d = x.shape
    def f(w):
        res = J(np.dot(x, w), y)
        #print(res)
        return res

    num_generations = 100
    num_parents_mating = 10

    sol_per_pop = 20
    num_genes = d

    # веса от 0 до 5 (для простоты)
    gene_space = [{'low': 0, 'high': 5} for _ in range(d)]

    def on_generation(ga_instance):
        print(f"Generation = {ga_instance.generations_completed}")
        print(f"Fitness    = {ga_instance.best_solution(pop_fitness=ga_instance.last_generation_fitness)[1]}")

    def fitness_func(ga_instance, solution, solution_idx):
        return -f(solution)

    ga_instance = pygad.GA(num_generations=num_generations,
                           num_parents_mating=num_parents_mating,
                           sol_per_pop=sol_per_pop,
                           num_genes=num_genes,
                           fitness_func=fitness_func,
                           on_generation=on_generation,
                           gene_space=gene_space)

    ga_instance.run()

    solution, solution_fitness, solution_idx = ga_instance.best_solution(ga_instance.last_generation_fitness)

    return solution


def show_spline(z, y):
    ind = np.argsort(z)
    cs = CubicSpline(z[ind], y[ind])
    z_grid = np.linspace(np.min(z), np.max(z), 1000)
    y_grid = cs(z_grid)
    vfunc = np.vectorize(test_f)
    y_test = vfunc(z_grid)
    plt.figure(figsize=(8, 5))
    plt.plot(z, y, 'o', label='Data Points')
    plt.plot(z_grid, y_grid, label='Cubic Spline Interpolation')
    plt.plot(z_grid, y_test, label='Test')
    plt.xlabel('z')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()



test_f = lambda x: math.sqrt(x)
test_w = [1, 2, 3]
test_w /= np.linalg.norm(test_w)
N = 100
d = 3
x = np.zeros((N, d))
y = np.zeros(N)
np.random.seed(123)
random.seed(123)
for i in range(N):
    for j in range(d):
        x[i, j] = np.random.random()
    y[i] = test_f(np.dot(x[i, :], test_w)) + 0.01 * (np.random.random() - 0.5)

w = train_spline(x, y)
print(w / w[0])

show_spline(np.dot(x, w) / np.linalg.norm(w), y)
Восстановленная зависимость
Восстановленная зависимость

На графике можно видеть зависимость y от z при найденных весах w = (1.0, 2.015, 2.987). Зеленым цветом выделена идеальная кривая, оранжевым - кубический сплайн.

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

Таким образом, эксперименты показывают работоспособность данного метода на тесте.

Комментарии (4)


  1. Daddy_Cool
    29.10.2025 23:00
    #29036824

    Ну здесь всё же не dxdy. Расписать бы это всё попопулярней и привести жизненный пример...
    ---
    "используем кубические сплайны класса C^2"

    " при растяжении точек \mathbf x_i в C раз функционал уменьшается в C^3 раз "

    Вам букв не хватило? ;)

    "В этой формуле функция f - это кубический сплайн"

    так кроме сплайнов-то f может еще чем-то быть или нет?


    1. lapkin25 Автор
      29.10.2025 23:00
      #29037026

      Данная задача возникла на таком примере. С помощью численных методов рассчитывалась зависимость температуры двух тепловых источников от их мощности. В результате получилась зависимость, похожая на линейную. Точнее, линии уровня на координатной плоскости были почти прямыми линиями, но не равноотстоящими. В связи с этим возникла идея: как оценить аналитически такую зависимость. Если наклон изолиний заранее неизвестен, а также неизвестна функция связи.

      так кроме сплайнов-то f может еще чем-то быть или нет?

      В примере f = sqrt(z), то есть это может быть любая функция.


  1. interesting-cs-math
    29.10.2025 23:00
    #29036866

    Интереснее было бы даже не сплайнами аппроксимировать. Сплайны больше для интерполяции.


    1. lapkin25 Автор
      29.10.2025 23:00
      #29037016

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