Минимизация функционала, регрессия

[vlad_light]

Подскажите, пожалуйста, правильно ли я делаю. Хочу научиться предсказывать распределение целевой переменной по вектору признаков. Пусть настоящее распределение $p(y)$ хорошо приближается смесью равномерных распределений $\sum _{k=1..K} \varphi _k p(y | a_k, b_k)$, где $Y_k \sim \mathcal U (a_i, b_i)$, $[a_k, b_k)\cap [a_j, b_j) = \empty, i \neq j$. Тогда по выборке $(\mathbf x_n, y_n)_{n=1..N}$ я хочу минимизировать функционал

$$\mathbb E [(Y - y)^2] \to \min,$$

где $Y$ из нашего приближенного распределения. Пусть наша модель $f(\mathbf x | \theta)$ будет предсказывать вектор $\mathbf z \in \mathbb R^K$. Тогда наш функционал примет вид:

$$\frac 1n\sum _{n=1..N} \sum _{k=1..K} s_k(\mathbf z_n) \int (x - y_n)^2p(x| a_k, b_k)dx=$$

$$\frac 1n\sum _{n=1..N} \sum _{k=1..K} s_k(\mathbf z_n) \left[\frac {b_k-a_k}{12} + \left(\frac {a_k + b_k}{2} - y_n\right)^2\right],$$

где $s_k(\mathbf z) = \frac {\exp(z_k)}{\sum _i \exp (z_i)}$. Если мы потребуем, чтоб дисперсия каждой компоненты смеси была одинаковой, и обозначим $\mu _k = \frac {a_k + b_k}{2}$, то получим более простой вид функционала:

$$\frac 1n\sum _{n=1..N} \sum _{k=1..K} s_k(\mathbf z_n) \left(\mu _k - y_n\right)^2$$

Вопрос: я правильно понимаю, что последний функционал не будет выпуклым и находить его минимум градиентным спуском не получится? Что тогда делать, какой-то ЕМ алгоритм?

[ipgmvq]

Немного запутался. Вектора $\mathbf a$ и $\mathbf b$ это данность, и единственные параметры для оптимизации тут это вектор $\mathbf \theta$, верно?

На всякий случай хотел уточнить: $p(y)$ и $p(x)$ это что? Плотность вероятности? Вероятность на интервале от $a_k$ до $b_k$? Функция распределения? Или в разных местах разные сущности?

[vlad_light]

Да, все, кроме $\theta$ у нас известно. $p(x)$ это плотность распределения, в данном случае

$$p(x|a,b)=\frac{1}{b-a}\mathbf 1_{[a, b)}$$

Также, явно не сказал, что $\mathbf z =\mathbf z(\theta)$.

[ipgmvq]

А исходя из названия темы ("регрессия") следует предполагать, что $f(x)$ это чисто линейная история?
$\mathbf x \times \theta = \mathbf z$
И если $\mathbf x$ это не матрица, но вектор, то $\theta$ это, получается, вектор с одним элементом? И у оптимизации только один параметр?

[vlad_light]

$f(\mathbf x | \theta )$ может быть линейной, например: $f(\mathbf x | \theta ) = (\theta _{0,k} + \sum _{m=1..M}x_m\theta _{m,k})_{k=1..K}$. Но хотелось бы использовать градиентный бустинг или нейронную сеть.
В примере выше, $\theta $ -- это матрица $(M+1)\times K$. Но в целом, это может быть произвольный набор параметров модели.

[ipgmvq]

vlad_light в сообщении #1586504 писал(а):

В примере выше, $\theta $ -- это матрица $(M+1)\times K$.

В плане выпуклости (и возможности оптимизировать путем аналитического нахождения производной), мне кажется, важно принять во внимание конкретную реализацию Вашей функции f(x).
Но предлагаемый функционал мне кажется странным.
Возмем к примеру упрощенную функцию f:
$\mathbf \Theta \times \mathbf x = \mathbf z $
Сразу выразим выборку $\mathbf y$ не в виде вектора абсолютных значений случайной величины, а в виде "индикаторной" матрицы, которая показывает, в какой интервал k попало индивидуальное значение в выборке:
$\mathbf Y = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \end{bmatrix}$
Также возмем вспомогательный вектор
$\mathbf i = \begin{bmatrix} -1 \\ -1 \\ -1 \\ -1 \end{bmatrix}$
Мне кажется, функционал должен быть:
$f(\mathbf \Theta) = \mathbf i^\intercal \times \mathbf Y \times \ln(s(\mathbf \Theta \times \mathbf x ))$

или (что эквивалентно)

$f(\mathbf \Theta) = \mathbf i^\intercal \times \mathbf Y \times \mathbf \Theta \times \mathbf x + n \times \ln{\sum\limits_{k}^{} e^{\mathbf \Theta \times \mathbf x} }$

Как Вам кажется?

Производная такого функционала вроде как аналитически нулевому вектору (тут матрице) равняться не может. А значит, у функции нет минимума и не помешали бы какие-то подходящие ограничения-неравенства для $\mathbf \Theta$

[ipgmvq]

vlad_light в сообщении #1586457 писал(а):

Тогда по выборке $(\mathbf x_n, y_n)_{n=1..N}$

vlad_light в сообщении #1586504 писал(а):

В примере выше, $\theta $ -- это матрица $(M+1)\times K$.

Неверно выразил идею.
Упрощенная функция f будет:
$ \mathbf{X} \times \mathbf{\Theta} = \mathbf{Z} $
И функционал соответственно:
$f(\mathbf \Theta) = -(\mathbf Y \cdot \ln ( s( \mathbf{X} \times \mathbf{\Theta} ))) $
где $\cdot$ — скалярное произведение матриц.

Ваш функционал (при отсутствии предсказательной ценности $ \mathbf{X} $) пытается минимизировать 2-й момент разницы двух независимых случайных величин, который как обычно равен сумме дисперсии этой разницы и квадрату мат ожидания. При этом дисперсия этой разницы равна сумме дисперсий выборки $ \mathbf{y} $ и дисперсии моделируемого распределения. А мат ожидание этой разницы по модулю равно модулю разницы выборочного среднего $ \mathbf{y} $ и мат ожидания моделируемого распределения. Соответственно оптимизация будет подтягивать матожидание моделируемого распределения к выборочной средней $ \mathbf{y} $ (что ок), но будет пытаться уменьшить дисперсию моделируемого распределения до нуля (а это совсем не то, что мы хотели).

[ipgmvq]

Иллюстрация применения на примере равномерного распределения и нормально распределенного экзогенного регрессора $\mathbf{X}$ без какой-либо предсказательной ценности. При $m$ равному или превосходящему $k$ на низких $n$ даст мощнейший overfitting.

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

import numpy as np
from scipy.optimize import minimize

n = 3000
m = 8
k = 10

s = lambda x: np.exp(x) / np.exp(x).sum(axis=1)

def f(x):
    A = np.matrix(x.reshape(m, k))
    return -np.tensordot(Y, np.log(s(X * A)))

X = np.matrix(np.random.normal(size=(n, m)))
y = np.random.uniform(low=0.0, high=1.0, size=n)

Y = np.matrix(np.zeros(shape=(n, k)))
for loc in zip(range(n),(y * 10.0).astype(int)):
    Y[loc[0], loc[1]] = 1.0

результат = minimize(fun=f, x0=np.ones(k * m))
    
прогноз = np.asarray(np.exp(X * np.matrix(результат.x.reshape(m, k))))
прогноз /= прогноз.sum(axis=1).reshape(n, 1)
print(прогноз.round(1))

print(f"Показатель отсутствия overfitting {np.abs(прогноз.round(0) - Y).sum()}")
 

[vlad_light]

Прошу прощения, что долго не отвечал.

Вы правы, ошибка в самом начале. Нужно минимизировать не $\mathbb E[(Y - y)^2]$, а правдоподобие модели $-\ln \mathcal L(Y, y)$.

Но в таком случае, если носители распределений $p_k(y) = p(y|a_k, b_k)$ не будут пересекаться, то задача превратиться в классификацию, что не совсем то, что я хотел.

У меня удалось добиться похожего эффекта, если взять в качестве $p_k = \mathcal N_k$, но это уже другая история

Спасибо большое за ответы, проблему нашли, исправили, тему можно закрывать.

(Оффтоп)

P.S. если вдруг интересно, я делаю custom loss для xgboost, который позволяет в качестве $f$ брать деревья, что значительно сложнее, чем линейная функция. Похожую штуку можно сделать для нейронных сетей, вот пример.