Основы линейной регрессии

    Здравствуй, Хабр!

    Цель этой статьи — рассказать о линейной регрессии, а именно собрать и показать формулировки и интерпретации задачи регрессии с точки зрения математического анализа, статистики, линейной алгебры и теории вероятностей. Хотя в учебниках эта тема изложена строго и исчерпывающе, ещё одна научно-популярная статья не помешает.

    ! Осторожно, трафик! В статье присутствует заметное число изображений для иллюстраций, часть в формате gif.

    Содержание




    Введение


    Есть три сходных между собой понятия, три сестры: интерполяция, аппроксимация и регрессия.
    У них общая цель: из семейства функций выбрать ту, которая обладает определенным свойством.


    Интерполяция — способ выбрать из семейства функций ту, которая проходит через заданные точки. Затем функцию обычно используют для вычисления в промежуточных точках. Например, мы вручную задаем цвет нескольким точкам и хотим чтобы цвета остальных точек образовали плавные переходы между заданными. Или задаем ключевые кадры анимации и хотим плавные переходы между ними. Классические примеры: интерполяция полиномами Лагранжа, сплайн-интерполяция, многомерная интерполяция (билинейная, трилинейная, методом ближайшего соседа и т.д). Есть также родственное понятие экстраполяции — предсказание поведения функции вне интервала. Например, предсказание курса доллара на основании предыдущих колебаний — экстраполяция.

    Аппроксимация — способ выбрать из семейства «простых» функций приближение для «сложной» функции на отрезке, при этом ошибка не должна превышать определенного предела. Аппроксимацию используют, когда нужно получить функцию, похожую на данную, но более удобную для вычислений и манипуляций (дифференцирования, интегрирования и т.п). При оптимизации критических участков кода часто используют аппроксимацию: если значение функции вычисляется много раз в секунду и не нужна абсолютная точность, то можно обойтись более простым аппроксимантом с меньшей «ценой» вычисления. Классические примеры включают ряд Тейлора на отрезке, аппроксимацию ортогональными многочленами, аппроксимацию Паде, аппроксимацию синуса Бхаскара и т.п.

    Регрессия — способ выбрать из семейства функций ту, которая минимизирует функцию потерь. Последняя характеризует насколько сильно пробная функция отклоняется от значений в заданных точках. Если точки получены в эксперименте, они неизбежно содержат ошибку измерений, шум, поэтому разумнее требовать, чтобы функция передавала общую тенденцию, а не точно проходила через все точки. В каком-то смысле регрессия — это «интерполирующая аппроксимация»: мы хотим провести кривую как можно ближе к точкам и при этом сохранить ее максимально простой чтобы уловить общую тенденцию. За баланс между этими противоречивыми желаниями как-раз отвечает функция потерь (в английской литературе «loss function» или «cost function»).

    В этой статье мы рассмотрим линейную регрессию. Это означает, что семейство функций, из которых мы выбираем, представляет собой линейную комбинацию наперед заданных базисных функций ${f_i}$

    $ f = \sum_i w_i f_i. $

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

    Регрессия с нами уже давно: впервые метод опубликовал Лежандр в 1805 году, хотя Гаусс пришел к нему раньше и успешно использовал для предсказания орбиты «кометы» (на самом деле карликовой планеты) Цереры. Существует множество вариантов и обобщений линейной регрессии: LAD, метод наименьших квадратов, Ridge регрессия, Lasso регрессия, ElasticNet и многие другие.

    гифка
    Точки генерируются случайно по распределению Гаусса с заданным средним и вариациями. Синяя линия — регрессионная прямая.

    Можно поиграться с демонстрацией в GoogleColab.
    Много других материалов по классическому машинному обучению на соответствующей страничке на GitHub


    Метод наименьших квадратов


    Начнём с простейшего двумерного случая. Пусть нам даны точки на плоскости $\{(x_1,y_1),\cdots,(x_N,y_N)\}$ и мы ищем такую аффинную функцию

    $ f(x) = a + b \cdot x, $


    чтобы ее график ближе всего находился к точкам. Таким образом, наш базис состоит из константной функции и линейной $(1, x)$.

    Как видно из иллюстрации, расстояние от точки до прямой можно понимать по-разному, например геометрически — это длина перпендикуляра. Однако в контексте нашей задачи нам нужно функциональное расстояние, а не геометрическое. Нас интересует разница между экспериментальным значением и предсказанием модели для каждого $x_i,$ поэтому измерять нужно вдоль оси $y$.

    Первое, что приходит в голову, в качестве функции потерь попробовать выражение, зависящее от абсолютных значений разниц $|f(x_i) - y_i|$. Простейший вариант — сумма модулей отклонений $\sum_i |f(x_i) - y_i|$ приводит к Least Absolute Distance (LAD) регрессии.

    Впрочем, более популярная функция потерь — сумма квадратов отклонений регрессанта от модели. В англоязычной литературе она носит название Sum of Squared Errors (SSE)

    $ \text{SSE}(a,b)=\text{SS}_{res[iduals]}=\sum_{i=1}^N{\text{отклонение}_i}^2=\sum_{i=1}^N(y_i-f(x_i))^2=\sum_{i=1}^N(y_i-a-b\cdot x_i)^2, $

    Метод наименьших квадратов (по англ. OLS) — линейная регрессия c $\text{SSE}(a,b)$ в качестве функции потерь.

    Такой выбор прежде всего удобен: производная квадратичной функции — линейная функция, а линейные уравнения легко решаются. Впрочем, далее я укажу и другие соображения в пользу $\text{SSE}(a,b)$.
    гифка
    Регрессионная прямая (синяя) и пробная прямая (зеленая). Справа показана функция потерь и точки соответствующие параметра пробной и регрессионной прямых.

    Можно поиграться с демонстрацией в GoogleColab.
    Много других материалов по классическому машинному обучению на соответствующей страничке на GitHub


    Математический анализ


    Простейший способ найти $\text{argmin}_{a,b} \, \text{SSE}(a,b)$ — вычислить частные производные по $ a $ и $ b $, приравнять их нулю и решить систему линейных уравнений

    $ \begin{aligned} \frac{\partial}{\partial a}\text{SSE}(a,b)&=-2\sum_{i=1}^N(y_i-a-bx_i), \\ \frac{\partial}{\partial b}\text{SSE}(a,b)&=-2\sum_{i=1}^N(y_i-a-bx_i)x_i. \end{aligned} $

    Значения параметров, минимизирующие функцию потерь, удовлетворяют уравнениям

    $ \begin{aligned} 0 &= -2\sum_{i=1}^N(y_i-\hat{a}-\hat{b}x_i), \\ 0 &= -2\sum_{i=1}^N(y_i-\hat{a}-\hat{b}x_i)x_i, \end{aligned} $

    которые легко решить

    $ \begin{aligned} \hat{a}&=\frac{\sum_i y_i}{N}-\hat{b}\frac{\sum_i x_i}{N},\\ \hat{b}&=\frac{\frac{\sum_i x_i y_i}{N}-\frac{\sum_i x_i\sum_i y_i}{N^2}}{\frac{\sum_i x_i^2}{N}-\left(\frac{\sum_i x_i^2}{N}\right)^2}. \end{aligned} $

    Мы получили громоздкие и неструктурированные выражения. Сейчас мы их облагородим и вдохнем в них смысл.

    Статистика


    Полученные формулы можно компактно записать с помощью статистических эстиматоров: среднего $\langle{\cdot}\rangle$, вариации $\sigma_{\cdot}$ (стандартного отклонения), ковариации $\sigma({\cdot},{\cdot})$ и корреляции $\rho({\cdot},{\cdot})$

    $ \begin{aligned} \hat{a}&=\langle{y}\rangle-\hat{b}\langle{x}\rangle, \\ \hat{b}&=\frac{\langle{xy}\rangle-\langle{x}\rangle\langle{y}\rangle}{\langle{x^2}\rangle-\langle{x}\rangle^2}. \end{aligned} $

    Перепишем $\hat{b}$ как

    $ \hat{b} = \frac{\sigma(x,y)}{\sigma_x^2}, $

    где $\sigma_x$ это нескорректированное (смещенное) стандартное выборочное отклонение, а $\sigma(x,y)$ — ковариация. Теперь вспомним, что коэффициент корреляции (коэффициент корреляции Пирсона)

    $ \rho(x,y)=\frac{\sigma(x,y)}{\sigma_x \sigma_y} $

    и запишем

    $ \hat{b}=\rho(x,y)\frac{\sigma_y}{\sigma_x}. $



    Теперь мы можем оценить все изящество дескриптивной статистики, записав уравнение регрессионной прямой так

    $ \boxed{y-\langle {y} \rangle = \rho(x,y)\frac{\sigma_y}{\sigma_x}(x-\langle {x} \rangle)}. $

    Во-первых, это уравнение сразу указывает на два свойства регрессионной прямой:
    • прямая проходит через центр масс $(\langle{x}\rangle, \langle{y}\rangle)$;
    • если по оси $x$ за единицу длины выбрать $\sigma_x$, а по оси $y$$\sigma_y$, то угол наклона прямой будет от $-45^\circ$ до $45^\circ$. Это связано с тем, что $-1 \leq\rho(x,y)\leq 1$.

    Во-вторых, теперь становится понятно, почему метод регрессии называется именно так. В единицах стандартного отклонения $y$ отклоняется от своего среднего значения меньше чем $x$, потому что $|\rho(x,y)|\leq1$. Это называется регрессией(от лат. regressus — «возвращение») по отношению к среднему. Это явление было описано сэром Фрэнсисом Гальтоном в конце XIX века в его статье «Регрессия к посредственности при наследовании роста». В статье показано, что черты (такие как рост), сильно отклоняющиеся от средних, редко передаются по наследству. Характеристики потомства как-бы стремятся к среднему — на детях гениев природа отдыхает.

    Возведя коэффициент корреляции в квадрат, получим коэффициент детерминации $R = \rho^2$. Квадрат этой статистической меры показывает насколько хорошо регрессионная модель описывает данные. $R^2$, равный $1$, означает что функция идеально ложится на все точки — данные идеально скоррелированны. Можно доказать, что $R^2$ показывает какая доля вариативности в данных объясняется лучшей из линейных моделей. Чтобы понять, что это значит, введем определения

    $ \begin{aligned} \text{Var}_{data} &= \frac{1}{N}\sum_i (y_i-\langle y \rangle)^2, \\ \text{Var}_{res} &= \frac{1}{N} \sum_i (y_i-\text{модель}(x_i))^2, \\ \text{Var}_{reg} &= \frac{1}{N} \sum_i (\text{модель}(x_i)-\langle y \rangle)^2. \end{aligned} $

    $\text{Var}_{data}$ — вариация исходных данных (вариация точек $y_i$).

    $\text{Var}_{res}$ — вариация остатков, то есть вариация отклонений от регрессионной модели — от $y_i$ нужно отнять предсказание модели и найти вариацию.

    $\text{Var}_{reg}$ — вариация регрессии, то есть вариация предсказаний регрессионной модели в точках $x_i$ (обратите внимание, что среднее предсказаний модели совпадает с $\langle y \rangle$).

    Дело в том, что вариация исходных данных разлагается в сумму двух других вариаций: вариации, которая объясняется моделью, и вариации случайного шума (остатков)

    $ \boxed{{\color{red}{\text{Var}_{data}}} ={\color{green}{\text{Var}_{res}}}+ {\color{blue}{\text{Var}_{reg}}}.} $

    или

    $ \sigma^2_{data} =\sigma^2_{res}+ \sigma^2_{reg}. $

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


    Мы стремимся избавиться от вариативности связанной с шумом и оставить лишь вариативность, которая объясняется моделью, — хотим отделить зерна от плевел. О том, насколько это удалось лучшей из линейных моделей, свидетельствует $R^2$, равный единице минус доля вариации ошибок в суммарной вариации

    $ R^2=\frac{\text{Var}_{data}-\text{Var}_{res}}{\text{Var}_{data}}=1-\frac{\color{green}{\text{Var}_{res}}}{\color{red}{\text{Var}_{data}}} $

    или доле объясненной вариации (доля вариации регрессии в полной вариации)

    $ R^2=\frac{\color{blue}{\text{Var}_{reg}}}{\color{red}{\text{Var}_{data}}}. $

    $R$ равен косинусу угла в прямоугольном треугольнике $(\sigma_{data}, \sigma_{reg}, \sigma_{res})$. Кстати, иногда вводят долю необъясненной вариации $FUV=1-R^2$ и она равна квадрату синуса в этом треугольнике. Если коэффициент детерминации мал, возможно мы выбрали неудачные базисные функции, линейная регрессия неприменима вовсе и т.п.

    Теория вероятностей


    Ранее мы пришли к функции потерь $\text{SSE}(a,b)$ из соображений удобства, но к ней же можно прийти с помощью теории вероятностей и метода максимального правдоподобия (ММП). Напомню вкратце его суть. Предположим, у нас есть $N$ независимых одинаково распределенных случайных величин (в нашем случае — результатов измерений). Мы знаем вид функции распределения (напр. нормальное распределение), но хотим определить параметры, которые в нее входят (например $\mu$ и $\sigma$). Для этого нужно вычислить вероятность получить $N$ датапоинтов в предположении постоянных, но пока неизвестных параметров. Благодаря независимости измерений, мы получим произведение вероятностей реализации каждого измерения. Если мыслить полученную величину как функцию параметров (функция правдоподобия) и найти её максимум, мы получим оценку параметров. Зачастую вместо функции правдоподобия используют ее логарифм — дифференцировать его проще, а результат — тот же.

    Вернемся к задаче простой регрессии. Допустим, что значения $x$ нам известны точно, а в измерении $y$ присутствует случайный шум (свойство слабой экзогенности). Более того, положим, что все отклонения от прямой (свойство линейности) вызваны шумом с постоянным распределением (постоянство распределения). Тогда

    $ y = a + bx + \epsilon, $

    где $\epsilon$ — нормально распределенная случайная величина

    $ \epsilon \sim \mathcal{N}(0,\,\sigma^{2}), \qquad p(\epsilon) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{\epsilon^2}{2\sigma^2}}. $




    Исходя из предположений выше, запишем функцию правдоподобия

    $ \begin{aligned} L(a,b|\mathbf{y})&=P(\mathbf{y}|a,b)=\prod_i P(y_i|a,b)=\prod_i p(y_i-a-bx|a,b)=\\ &= \prod_i \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(y_i-a-bx)^2}{2\sigma^2}}= \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{\sum_i (y_i-a-bx)^2}{2 \sigma^2}}=\\ &= \frac{1}{\sqrt{2 \pi \sigma^2}}e^{-\frac{\text{SSE}(a,b)}{2 \sigma^2}} \end{aligned} $

    и ее логарифм

    $ l(a,b|\mathbf{y})=\log{L(a,b|\mathbf{y})}=-\text{SSE}(a,b)+const. $

    Таким образом, максимум правдоподобия достигается при минимуме $\text{SSE}$

    $ (\hat{a},\hat{b})=\text{argmax}_{a,b} \, l(a,b|\mathbf{y}) = \text{argmin}_{a,b} \, \text{SSE}(a,b), $

    что дает основание принять ее в качестве функции потерь. Кстати, если

    $ \begin{aligned} \epsilon \sim \text{Laplace}(0, \alpha), \qquad p_{L}(\epsilon; \mu, \alpha) =\frac{\alpha}{2}e^{-\alpha |\epsilon-\mu|} \end{aligned} $

    мы получим функцию потерь LAD регрессии

    $ E_{LAD}(a,b)=\sum_i |y_i-a-bx_i|, $

    которую мы упоминали ранее.

    Подход, который мы использовали в этом разделе — один из возможных. Можно прийти к такому же результату, используя более общие свойства. В частности, свойство постоянства распределения можно ослабить, заменив на свойства независимости, постоянства вариации (гомоскедастичность) и отсутствия мультиколлинеарности. Также вместо ММП эстимации можно воспользоваться другими методами, например линейной MMSE эстимацией.

    Мультилинейная регрессия


    До сих пор мы рассматривали задачу регрессии для одного скалярного признака $x$, однако обычно регрессор — это $n$-мерный вектор $\mathbf{x}$. Другими словами, для каждого измерения мы регистрируем $n$ фич, объединяя их в вектор. В этом случае логично принять модель с $n+1$ независимыми базисными функциями векторного аргумента — $n$ степеней свободы соответствуют $n$ фичам и еще одна — регрессанту $y$. Простейший выбор — линейные базисные функции $(1, x_1, \cdots, x_n)$. При $n = 1$ получим уже знакомый нам базис $(1, x)$.

    Итак, мы хотим найти такой вектор (набор коэффициентов) $\mathbf{w}$, что

    $ \sum_{j=0}^n w_j x_j^{(i)}= \mathbf{w}^{\top}\mathbf{x}^{(i)} \simeq y_i, \qquad \qquad \qquad \qquad i=1\dots N. $

    Знак "$\simeq$" означает, что мы ищем решение, которое минимизирует сумму квадратов ошибок

    $ \hat{\mathbf{w}}=\text{argmin}_\mathbf{w} \, \sum_{i=1}^N \left({y_i - \mathbf{w}^{\top}\mathbf{x}^{(i)}}\right)^2 $

    Последнее уравнение можно переписать более удобным образом. Для этого расположим $\mathbf{x}^{(i)}$ в строках матрицы (матрицы информации)

    $ X= \begin{pmatrix} - & \mathbf{x}^{(1)\top} & - \\ \cdots & \cdots & \cdots\\ - & \mathbf{x}^{(N)\top} & - \end{pmatrix} = \begin{pmatrix} | & | & & | \\ \mathbf{x}_0 & \mathbf{x}_1 & \cdots & \mathbf{x}_n \\ | & | & & | \end{pmatrix} = \begin{pmatrix} 1 & x^{(1)}_{1} & \cdots & x^{(1)}_{n} \\ \cdots & \cdots & \cdots & \cdots\\ 1 & x^{(N)}_{1} & \cdots & x^{(N)}_{n} \end{pmatrix}. $

    Тогда столбцы матрицы $\mathbf{x}_{i}$ отвечают измерениям $i$-ой фичи. Здесь важно не запутаться: $N$ — количество измерений, $n$ — количество признаков (фич), которые мы регистрируем. Систему можно записать как

    $ X \, \mathbf{w} \simeq \mathbf{y}. $

    Квадрат нормы разности векторов в правой и левой частях уравнения образует функцию потерь

    $ \text{SSE}(\mathbf{w}) = {\|\mathbf{y}-X \mathbf{w}\|}^2, \qquad \qquad \mathbf{w} \in \mathbb{R}^{n+1}; \, \mathbf{y} \in \mathbb{R}^{N}, $

    которую мы намерены минимизировать

    $ \begin{aligned} \hat{\mathbf{w}}&=\text{argmin}_\mathbf{w} \, \text{SSE}(\mathbf{w}) = \text{argmin}_\mathbf{w} \, (\mathbf{y}-X \mathbf{w})^{\top}(\mathbf{y}-X \mathbf{w})=\\ &= \text{argmin}_\mathbf{w} \,(\mathbf{y}^{\top}\mathbf{y}-2\mathbf{w}^{\top}X^{\top}\mathbf{y}+\mathbf{w}^{\top}X^{\top}X\mathbf{w}). \end{aligned} $

    Продифференцируем финальное выражение по $\mathbf{w}$ (если забыли как это делается — загляните в Matrix cookbook)

    $ \frac{\partial \, \text{SSE}(\mathbf{w})}{\partial \mathbf{w}}=-2 X^{\top}\mathbf{y}+2 X^{\top}X\mathbf{w}, $

    приравняем производную к $\mathbf{0}$ и получим т.н. нормальные уравнения

    $ X^{\top}X \, \hat{\mathbf{w}}=X^{\top}\mathbf{y}. $

    Если столбцы матрицы информации $X$ линейно независимы (нет идеально скоррелированных фич), то матрица $X^{\top}X$ имеет обратную (доказательство можно посмотреть, например, в видео академии Хана). Тогда можно записать

    $ \boxed{\hat{\mathbf{w}} = (X^{\top}X)^{-1}X^{\top}\mathbf{y}=X^{+}\mathbf{y}}, $

    где

    $ X^{+}=(X^{\top}X)^{-1}X^{\top} $


    псевдообратная к $X$. Понятие псевдообратной матрицы введено в 1903 году Фредгольмом, она сыграла важную роль в работах Мура и Пенроуза.

    Напомню, что обратить $X^{\top}X$ и найти $X^{+}$ можно только если столбцы $X$ линейно независимы. Впрочем, если столбцы $X$ близки к линейной зависимости, вычисление $(X^{\top}X)^{-1}$ уже становится численно нестабильным. Степень линейной зависимости признаков в $X$ или, как говорят, мультиколлинеарности матрицы $X^{\top}X$, можно измерить числом обусловленности — отношением максимального собственного значения к минимальному. Чем оно больше, тем ближе $X^{\top}X$ к вырожденной и неустойчивее вычисление псевдообратной.


    Линейная алгебра


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

    $ X \, \mathbf{w} \simeq \mathbf{y}. $

    Если количество переменных равно количеству неизвестных и уравнения линейно независимы, то система имеет единственное решение. Однако, если число измерений превосходит число признаков, то есть уравнений больше чем неизвестных — система становится несовместной, переопределенной. В этом случае лучшее, что мы можем сделать — выбрать вектор $\mathbf{w}$, образ которого $X\mathbf{w}$ ближе остальных к $\mathbf{y}$. Напомню, что множество образов или колоночное пространство $\mathcal{C}(X)$ — это линейная комбинация вектор-столбцов матрицы $X$

    $ \begin{pmatrix} | & | & & | \\ \mathbf{x}_0 & \mathbf{x}_1 & \cdots & \mathbf{x}_n \\ | & | & & | \end{pmatrix} \mathbf{w} = w_0 \mathbf{x}_0 + w_1 \mathbf{x}_1 + \cdots w_n \mathbf{x}_n . $

    $\mathcal{C}(X)$$n+1$-мерное линейное подпространство (мы считаем фичи линейно независимыми), линейная оболочка вектор-столбцов $X$. Итак, если $\mathbf{y}$ принадлежит $\mathcal{C}(X)$, то мы можем найти решение, если нет — будем искать, так сказать, лучшее из нерешений.

    Если в дополнение к векторам $\mathcal{C}(X)$ мы рассмотрим все вектора им перпендикулярные, то получим еще одно подпространство и сможем любой вектор из $\mathbb{R}^{N}$ разложить на две компоненты, каждая из которых живет в своем подпространстве. Второе, перпендикулярное пространство, можно характеризовать следующим образом (нам это понадобится в дальнейшем). Пускай $\mathbf{v} \in \mathbb{R}^{N}$, тогда

    $ X^\top \mathbf{v} = \begin{pmatrix} - & \mathbf{x}_0^{\top} & - \\ \cdots & \cdots & \cdots\\ - & \mathbf{x}_n^{\top} & - \end{pmatrix} \mathbf{v} = \begin{pmatrix} \mathbf{x}_0^{\top} \cdot \mathbf{v} \\ \cdots \\ \mathbf{x}_n^{\top} \cdot \mathbf{v} \\ \end{pmatrix} $

    равен нулю в том и только в том случае, если $\mathbf{v}$ перпендикулярен всем $\mathbf{x}_i$, а значит и целому $\mathcal{C}(X)$. Таким образом, мы нашли два перпендикулярных линейных подпространства, линейные комбинации векторов из которых полностью, без дыр, «покрывают» все $\mathbb{R}^N$. Иногда это обозначают c помощью символа ортогональной прямой суммы


    где $\text{ker}(X^{\top})=\{\mathbf{v}|X^{\top}\mathbf{v}=\mathbf{0}\}$. В каждое из подпространств можно попасть с помощью соответствующего оператора проекции, но об этом ниже.

    Теперь представим $\mathbf{y}$ в виде разложения

    $ \mathbf{y} = \mathbf{y}_{\text{proj}} + \mathbf{y}_{\perp}, \qquad \mathbf{y}_{\text{proj}} \in \mathcal{C}(X), \qquad \mathbf{y}_{\perp} \in \text{ker}(X^{\top}). $



    Если мы ищем решение $\hat{\mathbf{w}}$, то естественно потребовать, чтобы $|| \mathbf{y} - X\mathbf{w} ||$ была минимальна, ведь это длина вектора-остатка. Учитывая перпендикулярность подпространств и теорему Пифагора

    $ \text{argmin}_\mathbf{w} || \mathbf{y} - X\mathbf{w} || = \text{argmin}_\mathbf{w} || \mathbf{y}_{\perp} + \mathbf{y}_{\text{proj}} - X\mathbf{w} || = \text{argmin}_\mathbf{w} \sqrt{|| \mathbf{y}_{\perp} ||^2 + || \mathbf{y}_{\text{proj}} - X\mathbf{w} ||^2}, $

    но поскольку, выбрав подходящий $\mathbf{w}$, я могу получить любой вектор колоночного пространства, то задача сводится к

    $ X\hat{\mathbf{w}} = \mathbf{y}_{\text{proj}}, $

    а $\mathbf{y}_{\perp}$ останется в качестве неустранимой ошибки. Любой другой выбор $\hat{\mathbf{w}}$ сделает ошибку только больше.

    Если теперь вспомнить, что $X^{\top} \mathbf{y}_{\perp} = \mathbf{0}$, то легко видеть

    $ X^\top X \mathbf{w} = X^{\top} \mathbf{y}_{\text{proj}} = X^{\top} \mathbf{y}_{\text{proj}} + X^{\top} \mathbf{y}_{\perp} = X^{\top} \mathbf{y}, $

    что очень удобно, так как $\mathbf{y}_{\text{proj}}$ у нас нет, а вот $\mathbf{y}$ — есть. Вспомним из предыдущего параграфа, что $X^{\top} X$ имеет обратную при условии линейной независимости признаков и запишем решение

    $ \mathbf{w} = (X^\top X)^{-1} X^\top \mathbf{y} = X^{+} \mathbf{y}, $

    где $X^{+}$ уже знакомая нам псевдообратная матрица. Если нам интересна проекция $\mathbf{y}_{\text{proj}}$, то можно записать

    $ \mathbf{y}_{\text{proj}} = X \mathbf{w} = X X^{+} \mathbf{y} = \text{Proj}_X \mathbf{y}, $

    где $\text{Proj}_X$ — оператор проекции на колоночное пространство.

    Выясним геометрический смысл коэффициента детерминации.

    Заметьте, что фиолетовый вектор $\bar{y} \cdot \boldsymbol{1}=\bar{y} \cdot (1,1,\dots,1)^{\top}$ пропорционален первому столбцу матрицы информации $X$, который состоит из одних единиц согласно нашему выбору базисных функций. В RGB треугольнике

    $ {\color{red}{\mathbf{y}-\hat{y} \cdot \boldsymbol{1}}}={\color{green}{\mathbf{y}-\bar{\mathbf{y}}}}+{\color{blue}{\hat{\mathbf{y}}-\bar{y} \cdot \boldsymbol{1}}}. $

    Так как этот треугольник прямоугольный, то по теореме Пифагора

    $ {\color{red}{\|\mathbf{y}-\hat{y} \cdot \boldsymbol{1}\|^2}}={\color{green}{\|\mathbf{y}-\bar{\mathbf{y}}\|^2}}+{\color{blue}{\|\hat{\mathbf{y}}-\bar{y} \cdot \boldsymbol{1}\|^2}}. $

    Это геометрическая интерпретация уже известного нам факта, что

    $ {\color{red}{\text{Var}_{data}}} = {\color{green}{\text{Var}_{res}}}+{\color{blue}{\text{Var}_{reg}}}. $

    Мы знаем, что

    $ R^2=\frac{\color{blue}{\text{Var}_{reg}}}{\color{red}{\text{Var}_{data}}}, $

    а значит

    $ R=\cos{\theta}. $

    Красиво, не правда ли?


    Произвольный базис


    Как мы знаем, регрессия выполняется на базисных функциях $f_i$ и её результатом есть модель

    $ f = \sum_i w_i f_i, $

    но до сих пор мы использовали простейшие $f_i$, которые просто ретранслировали изначальные признаки без изменений, ну разве что дополняли их постоянной фичей $f_0(\mathbf{x}) = 1$. Как можно было заметить, на самом деле ни вид $f_i$, ни их количество ничем не ограничены — главное чтобы функции в базисе были линейно независимы. Обычно, выбор делается исходя из предположений о природе процесса, который мы моделируем. Если у нас есть основания полагать, что точки $\{(x_1,y_1),\cdots,(x_N,y_N)\}$ ложатся на параболу, а не на прямую, то стоит выбрать базис $(1, x, x^2)$. Количество базисных функций может быть как меньшим, так и большим, чем количество изначальных фич.
    гифка
    Регрессия в полиномиальном базисе. Выделенная часть кода демонстрирует использование стандартных функций scikit-learn для выполнения регрессии полиномами разной степени, снизу — визуализация результата работы.

    Можно поиграться с демонстрацией в GoogleColab.
    Много других материалов по классическому машинному обучению на соответствующей страничке на GitHub

    Если мы определились с базисом, то дальше действуем следующим образом. Мы формируем матрицу информации

    $ \Phi = \begin{pmatrix} - & \boldsymbol{f}^{(1)\top} & - \\ \cdots & \cdots & \cdots\\ - & \boldsymbol{f}^{(N)\top} & - \end{pmatrix} = \begin{pmatrix} {f}_{0}\left(\mathbf{x}^{(1)}\right) & {f}_{1}\left(\mathbf{x}^{(1)}\right) & \cdots & {f}_{n}\left(\mathbf{x}^{(1)}\right) \\ \cdots & \cdots & \cdots & \cdots\\ {f}_{0}\left(\mathbf{x}^{(N)}\right) & {f}_{1}\left(\mathbf{x}^{(N)}\right) & \cdots & {f}_{n}\left(\mathbf{x}^{(N)}\right) \end{pmatrix}, $

    записываем функцию потерь

    $ E(\mathbf{w})={\|{\boldsymbol{\epsilon}}(\mathbf{w})\|}^2={\|\mathbf{y}-\Phi \, \mathbf{w}\|}^2 $

    и находим её минимум, например с помощью псевдообратной матрицы

    $ \hat{\mathbf{w}} = \text{argmin}_\mathbf{w} \,E(\mathbf{w}) = (\Phi^{\top}\Phi)^{-1}\Phi^{\top}\mathbf{y}=\Phi^{+}\mathbf{y} $

    или другим методом.


    Заключительные замечания


    Проблема выбора размерности


    На практике часто приходится самостоятельно построить модель явления, то есть определиться сколько и каких нужно взять базисных функций. Первый порыв «набрать побольше» может сыграть злую шутку: модель окажется слишком чувствительной к шумам в данных (переобучение). С другой стороны, если излишне ограничить модель, она будет слишком грубой (недообучение).

    Есть два способа выйти из ситуации. Первый: последовательно наращивать количество базисных функций, проверять качество регрессии и вовремя остановиться. Или же второй: выбрать функцию потерь, которая определит число степеней свободы автоматически. В качестве критерия успешности регрессии можно использовать коэффициент детерминации, о котором уже упоминалось выше, однако, проблема в том, что $R^2$ монотонно растет с ростом размерности базиса. Поэтому вводят скорректированный коэффициент

    $ \bar{R}^2=1-(1-R^2)\left[\frac{N-1}{N-(n+1)}\right], $

    где $N$ — размер выборки, $n$ — количество независимых переменных. Следя за $\bar{R}^2$, мы можем вовремя остановиться и перестать добавлять дополнительные степени свободы.

    Вторая группа подходов — регуляризации, самые известные из которых Ridge($L_2$/гребневая/Тихоновская регуляризация), Lasso($L_1$ регуляризация) и Elastic Net(Ridge+Lasso). Главная идея этих методов: модифицировать функцию потерь дополнительными слагаемыми, которые не позволят вектору коэффициентов $\mathbf{w}$ неограниченно расти и тем самым воспрепятствуют переобучению

    $ \begin{aligned} E_{\text{Ridge}}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\alpha \sum_i |w_i|^2 = \text{SSE}(\mathbf{w})+\alpha \| \mathbf{w}\|_{L_2}^2,\\ E_{\text{Lasso}}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\beta \sum_i |w_i| =\text{SSE}(\mathbf{w})+\beta \| \mathbf{w}\|_{L_1},\\ E_{\text{EN}}(\mathbf{w})&=\text{SSE}(\mathbf{w})+\alpha \| \mathbf{w}\|_{L_2}^2+\beta \| \mathbf{w}\|_{L_1}, \\ \end{aligned} $

    где $\alpha$ и $\beta$ — параметры, которые регулируют «силу» регуляризации. Это обширная тема с красивой геометрией, которая заслуживает отдельного рассмотрения. Упомяну кстати, что для случая двух переменных при помощи вероятностной интерпретации можно получить Ridge и Lasso регрессии, удачно выбрав априорное распределения для коэффициента $b$

    $ y = a + bx + \epsilon,\qquad \epsilon \sim \mathcal{N}(0,\,\sigma^{2}),\qquad \left\{\begin{aligned} &b \sim \mathcal{N}(0,\,\tau^{2})&\leftarrow\text{Ridge},\\ &b \sim \text{Laplace} (0,\,\alpha)&\leftarrow\text{Lasso}. \end{aligned}\right. $





    Численные методы


    Скажу пару слов, как минимизировать функцию потерь на практике. SSE — это обычная квадратичная функция, которая параметризируется входными данными, так что принципиально ее можно минимизировать методом скорейшего спуска или другими методами оптимизации. Разумеется, лучшие результаты показывают алгоритмы, которые учитывают вид функции SSE, например метод стохастического градиентного спуска. Реализация Lasso регрессии в scikit-learn использует метод координатного спуска.

    Также можно решить нормальные уравнения с помощью численных методов линейной алгебры. Эффективный метод, который используется в scikit-learn для МНК — нахождение псевдообратной матрицы с помощью сингулярного разложения. Поля этой статьи слишком узки, чтобы касаться этой темы, за подробностями советую обратиться к курсу лекций К.В.Воронцова.

    Реклама и заключение


    Эта статья — сокращенный пересказ одной из глав курса по классическому машинному обучению в Киевском академическом университете (преемник Киевского отделения Московского физико-технического института, КО МФТИ). Автор статьи помогал в создании этого курса. Технически курс выполнен на платформе Google Colab, что позволяет совмещать формулы, форматированные LaTeX, исполняемый код Python и интерактивные демонстрации на Python+JavaScript, так что студенты могут работать с материалами курса и запускать код с любого компьютера, на котором есть браузер. На главной странице собраны ссылки на конспекты, «рабочие тетради» для практик и дополнительные ресурсы. В основу курса положены следующие принципы:

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

    Если хотите посмотреть на результат — загляните на страничку курса на GitHub.

    Надеюсь вам было интересно, спасибо за внимание.
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама

    Комментарии 34

      +1
      Всегда думал, что интерполяция — это выбор из семейства функций уже проходящих через заданные точки, поскольку их также бесконечное множество.
        0
        Нам когда-то давали определение, что «интер» значит приближение внутри интервала, в котором рассыпаны точки. «экстра» — выход за пределы этого интервала.
          0

          Да, Вы правы, я внес правку для прояснения этого момента

            0
            Я не совсем понял ваш ответ, потому что экстраполяция — это вопрос интерпретации конкретной мат.модели и неважно, проходит она через какие-то явно заданные точки или нет. Никто не мешает интерполяцию использовать для экстраполяции:


            А в ЦОС подобное разделение так и вообще не имеет смысла — там что интерполяция, что экстраполяция равнозначны.
              0

              Кстати, есть еще Безье-интерполяция, которая не укладывается в данное определение интерполяции, ведь часть точек — управляющие.

                0
                Управляющие точки в Безье — это чисто визуальная абстракция, задающие значения производной, их можно вполне считать автоматически используя какую-нибудь взвешенную функцию. А по (математической) сути Безье-интерполяция — это всё та же кусочно-линейная полиномиальная интерполяция.
                  0

                  Я понимаю, что кривую Безье в каждой точке можно вычислить с помощью последовательностей кусочно-линейных интерполяций и что в конце получится функция линейная по точкам и полиномиальная по параметру, но исходная задача не подпадает ни под определение интерполяции ни под определение кусочной интерполяции.

                0
                Никто не мешает интерполяцию использовать для экстраполяции:


                Главное, чтобы был результат и понимание процесса. Поэтому автор поступил грамотно — дал определение терминов в начале.
            +3
            И если разрешите позанудствовать, то мне искренне непонятна тенденция излагать предельно простые вещи предельно тяжёлым матаном, и ваша статья среди прочих статей о линейной регрессии в этом соревновании безусловно лидирует. Всегда думал, что научно-популярное — это когда сложное излагается простым языком, а не наоборот.
            Как весь этот матан в статье позволит перейти к аппроксимации комплексными синусоидами, например? Как он он поможет победить численную нестабильность в расчётах на большом количестве данных и полиномах степени выше 1?
              +1
              До этого в науке часто доходят путём «написал просто, а потом Вася прочёл и указал на одну небольшую неточность, так что я уточнил вот такой формулой» (повторить 45 раз). В итоге все профессиональные Васи счастливы точности формулировок, но по ходу дела мы потеряли 45 000 потенциальных читателей.
                0
                Редко, но бывает, что иногда очень не хватает именно матана. В сложных случаях. Когда недостаточно основ и примеров из википедии.
                Мне статья понравилась.
                  +1
                  Статья действительно предполагает несколько подготовленного читателя, но всё же она далека от жёсткого матана, когда изложение строится на формальных определениях и доказательствах, сильно зацепляющихся друг за друга, и где ослабив внимание или недопоняв хоть что-то одно — всё, до свидания, ну или возвращайся назад и разбирайся заново. При написании я ориентировался на сложность/строгость научно-популярных статей журнала «Квант» и брошюр серии «Математическое просвещение» — там отнюдь на сюсюкаются с читателем.

                  Если уж задаться вопросом, кому адресована эта статья, то это скорее читатель, который уже видел линейную регрессию, но знает её «только с одной стороны» (например, один раз делал на лабе по готовым формулам). Здесь же удалось собрать целый ряд разных подходов: с точки зрения матана, теорвера, линейки и т.д. Чем с большего числа сторон Вы взглянете на явление, тем глубже его поймете. Иногда задача решается удачным выбором точки зрения.

                  По поводу регрессии синусоидами, то процитирую себя же из статьи: «Отмечу, что линейную регрессию называют линейной именно из-за линейной комбинации базисных функций — это не связано с самыми базисными функциями (они могут быть линейными или нет).» Берите синусы в качестве базисных функций и будут Вам синусы. Более того, гифка в разделе «произвольный базис» показывает регрессию в полиномиальном базисе.

                  По поводу нестабильности: это отдельная тема и здесь я лишь упомянул, что она существует. Охватить всё в одной статье невозможно, да и чтоб начать разбираться с нестабильностью, неплохо бы сначала разобраться в основах.
                    +1
                    Берите синусы в качестве базисных функций и будут Вам синусы
                    То есть готовых, конкретных и пригодных к кодированию формул опять не будет. И не «синусы», а «комплексные синусоиды» — это не одно и то же. Если бы я на просьбу математика о программе для численных вычислений ответил «в учебнике [Страуструпа] всё написано, берите компилятор с++ и будет вам счастье» — вряд ли ему такой ответ понравился бы. И вряд ли бы он смог аналог матлаба наваять.
                      0

                      Я не уверен что такое "комплексная синусоида". Это exp(i x), где x — действительное число, а i — мнимая единица, или sin(Z), где Z комплексное, или Z sin(x), где x — действительное число, а Z — комплексное, или это что-то ещё? Сначала сформулируйте адекватно задачу: укажите что и на каком семействе функций Вы хотите "регрессировать" и тогда я смогу ответить, знаю ли я как такое решается и если да, то можно ли получить этот ответ из статьи. Кстати, если я расскажу про функции комплексного переменного, то будут претензии, что я ничего не написал про кватернионные синусы? Так нужно сразу написать про функции, действующие на любом линейном пространстве? Но ведь были претензии, что статья недостаточно элементарная. Теперь она уже недостаточно продвинутая?


                      Раз уже затронули комплексные числа, то на всякий случай прокомментирую их. Задача линейной регрессии совершенно "прямолинейно" обобщается и на комплексные числа. Если переписать всю статью в терминах комплексных чисел, то единственная заметная разница — вместо транспонированных матриц будут стоять эрмитово-сопряжённые, ну может ещё где-то появятся знаки сопряжения. Но приводить линейную регрессию для комплексных чисел я считаю излишеством: во-первых, это для понимания сложнее обычной регрессии (хотя формулы почти те же) и начинать с такого явно не стоит, во-вторых, это всё-таки довольно редкая задача — даже sklearn.LinearRegression выдаёт "Complex data not supported\n"


                      По поводу готовых формул для кодирования. Видите ли, уметь что-то запрограммировать по готовым формулам и понимать математический аппарат за ними — это два разных навыка. Эта статья ориентирована на развитие второго и для задачи "просто надо чтоб работало" будет малополезной. Но и так есть довольно много других хороших статей и лекций на ютубе, где рассматривается именно практическая часть, будь то самостоятельная имплементация регрессии или использование готовых библиотечных функций. А вот статей, которые показывают "математическое закулисье" мало, вот на эту нишу я и ориентировался.

                        0
                        Я не уверен что такое «комплексная синусоида». Это exp(i x), где x — действительное число, а i — мнимая единица, или sin(Z), где Z комплексное, или Z sin(x), где x — действительное число, а Z — комплексное, или это что-то ещё?
                        Гугл даёт вполне однозначный ответ — хотя я предпочитаю более компактную запись в виде ix (наверняка её можно записать ещё более компактно с помощью статистических эстиматоров — только у меня не получилось). Конечно, можно рассмотреть действительную и мнимые компоненты по отдельности — но это лишает возможности контролировать амплитуду и фазу при аппроксимации.

                        Сначала сформулируйте адекватно задачу
                        Задачи типичные — шумоподавление квадратурного сигнала, аппроксимация спектра импульсной характеристики инвариантной к временной задержке, определение частоты.
                          0

                          Я не занимаюсь теорией обработки сигналов и могу не совсем понимать условие. Правильно ли я понимаю, что задача состоит, например, в следующем: для некоторой функции f, которая являет собой линейную комбинацию A_1 exp(i w_1 t) + A_2 exp(i w_2 t) +...+ A_2 exp(i w_N t) (сигнал) плюс случайная величина e (шум), необходимо найти A_1,...,A_N при известных w_1,...,w_n? Если да, то нужно воспользоваться разделом "Произвольный базис" и установить там базисные функции exp(i w_1 t),...,exp(i w_N t), а потом воспользоваться приведённой там формулой с учётом моего предыдущего комментария (транспонирование -> эрмитово сопряжение). Long story short, вот демонстрационный кусок кода для Mathematica


                          Remove["Global`*"]
                          basisF[t_] := {Exp[0.1*t*I], Exp[0.2*t*I], Exp[0.15*t*I]};
                          signalA = {5, 2, 3};
                          noiseA = 0.1/Sqrt[2];
                          sampleN = 100;
                          
                          signal[t_] := signalA.basisF[t];
                          noise[t_] := 
                            noiseA*(RandomVariate[NormalDistribution[0, 1]] + 
                               I*RandomVariate[NormalDistribution[0, 1]]);
                          Print["Signal to noise ratio = ", Norm[signalA]/noiseA]
                          (*Generate noisy signal*)
                          y = Table[N[signal[t] + noise[t]], {t, 0, sampleN}];
                          (*Construct features*)
                          \[CapitalPhi] = Table[N[basisF[t]], {t, 0, sampleN}];
                          (*Apply formula*)
                          w[y_, \[CapitalPhi]_] := 
                            Inverse[ConjugateTranspose[\[CapitalPhi]].\[CapitalPhi]].\
                          ConjugateTranspose[\[CapitalPhi]].y;
                          (*Result*)
                          recovered = Abs[w[y, \[CapitalPhi]]];
                          Print["Initial amplitudes:", signalA]
                          Print["Recovered amplitudes:", recovered]

                          Signal to noise ratio = 87.178
                          Initial amplitudes:{5,2,3}
                          Recovered amplitudes:{5.00467,2.01432,2.99681}

                          Ясно что при определённых условиях здесь могут быть проблемы, например кратные частоты при слишком малой частоте дискретизации могут привести к двум линейно зависимым столбцам матрицы. В статье эта проблема упоминается. Если это не та задача, что Вам интересна, то сформулируйте условие более понятно для неспециалиста в обработке сигналов.

                            0
                            Спасибо за код. Там отсутствует один принципиальный момент — инициальная фаза, то есть базис имеет вид
                            a*Exp[k*t*I+phase]
                            при аппроксимации которым ищется не только a, но и phase.

                            скрытый текст
                            Всё потому, что ЦОС рассматривает бесконечные во времени сигналы и инициальная фаза передатчика нам неизвестна.
                              0

                              На самом деле, она там есть. w[y, [CapitalPhi]] — комплексное число и как у каждого уважающего себя комплексного числа у него есть и модуль Abs[w[y, [CapitalPhi]]], который я вывожу, и фаза, которую печатать не стал (всё равно в начальных условиях она везде ноль). Вот очень минималистичный пример с фазой


                              sampleN = 100;
                              signal[t_] = 
                                5*Exp[0.1*t*I + 0.5*I] + 2*Exp[0.2*t*I + 0.8*I] + 
                                 3*Exp[0.15*t*I + 1.5*I];
                              noise[t_] := 
                                0.1*(RandomVariate[NormalDistribution[0, 1]] + 
                                   I*RandomVariate[NormalDistribution[0, 1]]);
                              y = Table[N[signal[t] + noise[t]], {t, 0, sampleN}];
                              \[CapitalPhi] = Table[N[basisF[t]], {t, 0, sampleN}];
                              w[y_, \[CapitalPhi]_] := 
                                Inverse[ConjugateTranspose[\[CapitalPhi]].\[CapitalPhi]].\
                              ConjugateTranspose[\[CapitalPhi]].y;
                              AbsArg[w[y, \[CapitalPhi]]]

                              {{4.99847, 0.501885}, {1.9858, 0.802154}, {3.01153, 1.50344}}
                                0

                                UPD базисные функции надо взять как в предыдущем примере


                                basisF[t_] := {Exp[0.1*t*I], Exp[0.2*t*I], Exp[0.15*t*I]};

                                (потерял строчку при копировании кода)

                          –1
                          Раз уже затронули комплексные числа, то на всякий случай прокомментирую их
                          Давайте и я прокомментирую — всё, что может быть выражено через комплексные (и прочие) числа, должно через них и выражаться. Gimbal lock — наглядный пример ситуации, когда матричная алгебра заходит в тупик. Ещё наглядный пример — вот в этой статье вы можете сравнить решение автора традиционным подходом и решения через комплексные числа в комментариях.
                            0

                            1) Некоторые задачи естественно решаются с помощью одного аппарата, другие — с помощью другого. В частности, комплексные числа особенно хороши для двумерных задач, красиво решают задачи электротехники и т.п. При этом комплексные числа вовсе не исключают векторы, так например, при решении задачи движения нерелятивистской частицы со спином удобно объединить два комплексных числа в вектор.
                            Вообще, матрицы сами по себе, это просто таблицы с числами. Обычно эти числа — координаты чего-то (например тензора) в определенном базисе. Смысл в эти таблички вдыхает линейная алгебра (тензорный анализ).


                            2) Мне понравились Ваши комплексные гильоши ^^


                            3) Gimbal lock возникает не столько из-за матриц как таковых, сколько от неудачной параметризации (углы Эйлера). Можно выбрать другую параметризацию в матрицах или использовать единичные кватернионы. Кватернионы и правда удобны для представления вращений, но тоже не лишены проблем параметризации (двойное покрытие). Вообще, споры о том, что лучше — кватернионы или 3D векторы — восходит к Гиббсу и Гамильтону и мне тоже немного жалко, что сейчас кватернионы используются намного реже векторов.

                              0
                              Можно выбрать другую параметризацию в матрицах или использовать единичные кватернионы
                              Так в этом-то и дело — чтобы описать через матрицы кватернионы, нужно заранее знать про кватернионы. Имея на руках матрицы поворота, полученные из комплексных чисел, нельзя получить из них матрицу для умножения кватернионов, просто увеличив размерность.

                              Также вычисления матрицами носит дискретный характер и посчитать производную или промежуточные значения будет непросто.

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

                              что лучше — кватернионы или 3D векторы
                              Если поставить задачу о построении спирали, закрученной в спираль, закрученной в спираль (повторить n раз) и посчитать её длину, то ответ будет однозначным.
                                0
                                Также вычисления матрицами носит дискретный характер и посчитать производную или промежуточные значения будет непросто.

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


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

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


                                Если поставить задачу о построении спирали, закрученной в спираль, закрученной в спираль (повторить n раз) и посчитать её длину, то ответ будет однозначным.

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

                                  0
                                  Я не утверждаю, что матрицы хуже — я убеждён, что матрицы вторичны, и вполне допускаю, что моё убеждение ошибочно. В частности, я не знал, что матрицы можно дифференцировать.

                                  Окей, я понимаю что раньше матрицы сильно облегчали вычисления, но сейчас же у нас есть системы компьютерной алгебры! Когда мне нужно решить систему уравнений, я и вызываю метод с названием «решить систему уравнений», а СКА уже сама определяет, является ли она линейной, недо/переопределённой и т.д. Нет никакой необходимости заполнять матрицу вручную.

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

                                  Complex a;
                                  Dual b;
                                  Dual<Complex> c;


                                  а не

                                  Matrix<2,1> complex_a;
                                  Matrix<2,1> dual_b;
                                  Matrix<2,2> dual_of_complex_c;


                                  чтобы при их умножении соответствующие матрицы автоматически выбирались в зависимости от типа.
                              0
                              Мне понравились Ваши комплексные гильоши ^^
                              Спасибо)

                              Если вместо комплексных чисел взять кватернионы, то можно строить узоры на поверхности сферы. В статике они, правда, выглядят не особо интересно:






                              А делалось это не ради узоров, а ради синтеза звука.
                        –1
                        гифка в разделе «произвольный базис» показывает регрессию в полиномиальном базисе
                        Эта гифка сделана в стороннем инструменте используя метод под названием «полиномиальная регрессия», а вовсе не комбинацию статистических эстиматоров. Это не выглядит похожим на численную проверку ваших теоретических изысканий.
                          0

                          К сожалению, это не моё теоретическое изыскание, а вполне стандартный подход, преподаваемый во всех курсах ML. Вот, например, этот подход в известном курсе Воронцова http://www.machinelearning.ru/wiki/images/a/a2/Voron-ML-regression-slides.pdf (4-й слайд). На самом деле, я был бы искренне горд собой, будь это моим теоретическим изысканием.


                          По поводу стороннего инструмента, претензия мне совсем непонятна. С чего вдруг стало зазорно при иллюстрации подхода приводить примеры известных библиотек где он используется? И в выделеном куске кода как раз прекрасно видно, что "Полиномиальная регрессия" — это обычная линейная регрессия с базисом из полиномов (команда: make_pipeline(PolynomialFeatures(.), LinearRegression())). И даже можно догадаться, что PolynomialFeatures(.) строит нам строчки матрицы Phi (см. раздел перед гифкой). Или описание теории можно сопровождать исключительно своими велосипедами с квадратными колёсами? Или при упоминании третьесторонних библиотек ни в коем случае нельзя описывать на чём они базируются — преподавать как магию, а код как заклинания?

                      0

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

                        0
                        Могу предложить вам ещё одну интересную задачу: посчитать линейную регрессию над произвольным (в том числе и превышающим доступную оперативную память компьютера) количеством данных за O(1). Решение, естественно, должно выводится из матриц.
                          0

                          Я полагаю, под О(1) имеется в виду независимость от количества датапоинтов N, но не от размерности задачи n. В противном случае, я сомневаюсь в существовании подобных алгоритмов — результат работы регрессии составляют n чисел и независимый от размерности алгоритм не смог бы их даже распечатать "в один присест".


                          Если же достаточно независимости от N, то в некоторых разумных предположениях можно придумать подобный алгоритм, как раз основываясь на столь любимых Вами матрицах. Для этого обратимся к формуле раздела "Произвольный базис" и обратим внимание, что часть Phi^T y имеет размерность 1 х (n+1)



                          и всё остальное под минус первой степенью имеет размерность (n+1) х (n+1).



                          Основываясь на матричном умножении, легко доказать следующие два утверждения, проиллюстрированные графически




                          Строгих доказательств приводить не буду — любой видевший хоть обложку книги по матанализу должен уметь доказывать такое самостоятельно. Вместо этого приведу демонстрационный код


                          import numpy as np
                          
                          # does regression
                          class Regressor:
                           def __init__(self, dimensions):
                             self.XTy = np.zeros(dimensions)
                             self.XTX = np.zeros((dimensions, dimensions))
                          
                           def push_features(self, f, y):
                             self.XTy += f * y
                             self.XTX += f * f[np.newaxis].T
                          
                           def get_regression_params(self):
                             return np.dot(np.linalg.inv(self.XTX), self.XTy)
                          
                          # generates features from raw data
                          features = lambda x: np.concatenate(([1.0], x))
                          
                          f = 10 # number of dimensions
                          r = Regressor(f + 1)
                          
                          # create a dataset
                          N = 100 # numer of datapoints
                          np.random.seed(42)
                          X = np.random.rand(N, f)
                          Y = np.random.rand(N, 1) # should always be 1
                          
                          # test regressor
                          for x,y in zip(X, Y):
                           r.push_features(features(x), y)
                          
                          w = r.get_regression_params()
                          print('my:', w)
                          
                          # do the same with scikit-learn
                          from sklearn.linear_model import LinearRegression
                          model = LinearRegression(fit_intercept=True) # Regressor fits with intersept
                          model.fit(X, Y)
                          print('sklearn:', model.intercept_, model.coef_)

                          my: [ 0.50998109  0.04065423 -0.22417411  0.12527675 -0.01315162  0.06585738
                            0.05293676  0.06372162 -0.03954708  0.07962786 -0.05955656]
                          sklearn: [0.50998109] [[ 0.04065423 -0.22417411  0.12527675 -0.01315162  0.06585738  0.05293676
                             0.06372162 -0.03954708  0.07962786 -0.05955656]]

                          Вся процедура происходит в классе Regressor. Он специально так спроектирован, чтоб подчёркивать независимость от количества датапоинтов N — их даже подавать туда надо по одной через push_features, потому что класс без понятия, сколько их будет. Вся память класса резервируется в конструкторе и больше никогда не увеличивается. С количеством датапоинтов я определяюсь уже после создания экземпляра класса. Не хотел усложнять код, но всё же написал сразу с генерацией фич features. Можно заказать r = Regressor(f) и подавать сырые иксы r.push_features(x, y), тогда результатом будет регрессия без свободного члена, эквивалент LinearRegression(fit_intercept=False).


                          По поводу сделанных предположений. Я полагаю, что выполняемые мной арифметические действия не приводят к переполнениям или большой потере точности. Ясно, что на практике это потребует дополниительного исследования и хотя бы "обвязки ифами", чтоб предупредить пользователя, что произошло переполнение или точность работы скомпрометирована.

                            0
                            Задача имеет решение при некоторых дополнительных условиях, когда её можно свести к рекурсивному фильтру. В самом простом случае получится скользящее среднее, которое на каждом этапе вычисления пересчивается как
                            avg += input[value]
                            avg -= delayed[value]

                            вне зависимости от количества усредняемых данных.

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

                            Для обеспечения устойчивости вычислений нужно использовать либо целочисленную арифметику, либо вводить обратную связь для компенсации накопленных погрешностей.
                          0
                          Спасибо, хоть наконец-то нормальное определение регрессии встретил. Хотя до сих пор не понимаю, почему такое слово используется. Регресс — это противоположность прогресса. То есть что-то типа примитивизации. А тут больше похоже не редукцию, то есть сокращение количества параметров и приближение ортогональными функциами. При этом набор ортогональных функций может быть оочень далек от примитива.
                            0

                            В разделе "Статистика" я пояснил откуда взялось такое название ("регрессия к среднему"). Вот кстати неплохое науч-поп видео про регрессию к среднему. Слово "редукция" тоже подходит, как по мне, но про dimension reduction обычно говорят в контексте метода главных компонент (PCA).

                          Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                          Самое читаемое