How to become an author
Search
Write a publication
Pull to refresh

Приводим уравнение линейной регрессии в матричный вид

Reading time6 min
Views24K
Mathematics*Machine learning*Studying in IT


Цель статьи — оказание поддержки начинающим датасайнтистам. В предыдущей статье мы на пальцах разобрали три способа решения уравнения линейной регрессии: аналитическое решение, градиентный спуск, стохастический градиентный спуск. Тогда для аналитического решения мы применили формулу $X^T X \vec{w} = X^T \vec{y}$. В этой статье, как следует из заголовка, мы обоснуем применение данной формулы или другими словами, самостоятельно ее выведем.

Почему имеет смысл уделить повышенное внимание к формуле $X^T X \vec{w} = X^T \vec{y}$?

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

Например, на курсах по машинному обучению от Яндекса, когда слушателей знакомят с регуляризацией, то предлагают воспользоваться функциями из библиотеки sklearn, при этом ни слова не упоминается о матричном представлении алгоритма. Именно в этот момент у некоторых слушателей может появится желание разобраться в этом вопросе подробнее — написать код без использования готовых функций. А для этого, надо сначала представить уравнение с регуляризатором в матричном виде. Данная статья, как раз, позволит желающим овладеть такими умениями. Приступим.

Исходные условия


Целевые показатели


У нас имеется ряд значений целевого показателя. Например, целевым показателем может быть цена на какой-либо актив: нефть, золото, пшеница, доллар и т.д. При этом, под рядом значений целевого показателя мы понимаем количество наблюдений. Такими наблюдениями могут быть, например, ежемесячные цены на нефть за год, то есть у нас будет 12 значений целевого показателя. Начнем вводить обозначения. Обозначим каждое значение целевого показателя как $y_i$. Всего мы имеем $n$ наблюдений, а значит можно представить наши наблюдения как $y_1, y_2, y_3 ... y_n$.

Регрессоры


Будем считать, что существуют факторы, которые в определенной степени объясняют значения целевого показателя. Например, на курс пары доллар/рубль сильное влияние оказывает цена на нефть, ставка ФРС и др. Такие факторы называются регрессорами. При этом, каждому значению целевого показателя должно соответствовать значение регрессора, то есть, если у нас имеется 12 целевых показателей за каждый месяц в 2018 году, то и значений регрессоров у нас тоже должно быть 12 за тот же период. Обозначим значения каждого регрессора через $x_i: x_1, x_2, x_3 ... x_n$. Пусть в нашем случае имеется $k$ регрессоров (т.е. $k$ факторов, которые оказывают влияние на значения целевого показателя). Значит наши регрессоры можно представить следующим образом: для 1-го регрессора (например, цена на нефть): $x_{11}, x_{12}, x_{13} ... x_{1n}$, для 2-го регрессора (например, ставка ФРС): $x_{21}, x_{22}, x_{23} ... x_{2n}$, для "$k$-го" регрессора: $x_{k1}, x_{k2}, x_{k3} ... x_{kn}$

Зависимость целевых показателей от регрессоров


Предположим, что зависимость целевого показателя $y_i$ от регрессоров "$i$-го" наблюдения может быть выражена через уравнение линейной регрессии вида:

$ f(w,x_i) = w_0 + w_1 x_{1i} + ... + w_k x_{ki} $



, где $x_i$ — "$i$-ое" значение регрессора от 1 до $n$,

$k$ — количество регрессоров от 1 до $k$

$w$ — угловые коэффициенты, которые представляют величину, на которую изменится расчетный целевой показатель в среднем при изменении регрессора.

Другими словами, мы для каждого (за исключением $w_0$) регрессора определяем «свой» коэффициент $w$, затем перемножаем коэффициенты на значения регрессоров "$i$-го" наблюдения, в результате получаем некое приближение "$i$-го" целевого показателя.

Следовательно, нам нужно подобрать такие коэффициенты $w$, при которых значения нашей апроксимирующей функции $f(w,x_i)$ будут расположены максимально близко к значениям целевых показателей.

Оценка качества апроксиммирующей функции


Будем определять оценку качества апроксимирующей функции методом наименьших квадратов. Функция оценки качества в таком случае примет следующий вид:

$Err= \sum\limits_{i=1}^n(y_i-f(x_i))^2 \rightarrow min$



Нам требуется подобрать такие значения коэффициентов $w$, при которых значение $Err$ будет наименьшим.

Переводим уравнение в матричный вид


Векторное представление


Для начала, чтобы облегчить себе жизнь, следует обратить внимание на уравнение линейной регрессии и заметить, что первый коэффициент $w_0$ не умножается ни на один регрессор. При этом, когда мы переведем данные в матричный вид, вышеобозначенное обстоятельство будет серьезно осложнять расчеты. В этой связи предлагается ввести еще один регрессор для первого коэффициента $w_0$ и приравнять его единице. Вернее, каждое "$i$-ое" значение этого регрессора приравнять единице — ведь при умножении на единицу у нас с точки зрения результата вычислений ничего не изменится, а с точки зрения правил произведения матриц, существенно сократятся наши мучения.

Теперь, на некоторое время, с целью упрощения материала, предположим, что у нас только одно "$i$-ое" наблюдение. Тогда, представим значения регрессоров "$i$-ого" наблюдения в качестве вектора $\vec{x_i}$. Вектор $\vec{x_i}$ имеет размерность $(k \times 1)$, то есть $k$ строк и 1 столбец:

$\vec{x_i} = \begin{pmatrix} x_{0i} \\ x_{1i} \\ ... \\ x_{ki} \end{pmatrix} \qquad$



Искомые коэффициенты представим в виде вектора $\vec{w}$, имеющего размерность $(k \times 1)$:

$\vec{w}=\begin{pmatrix} w_0 \\ w_1 \\ ... \\ w_k \end{pmatrix} \qquad$



Уравнение линейной регрессии для "$i$-го" наблюдения примет вид:

$f(w,x_i) = \vec{x_i}^T \vec{w}$



Функция оценки качества линейной модели примет вид:

$Err= \sum\limits_{i=1}^n(y_i-\vec{x_i}^T \vec{w})^2 \rightarrow min$



Обратим внимание, что в соответствии с правилами умножения матриц, нам потребовалось транспонировать вектор $\vec{x_i}$.

Матричное представление


В результате умножения векторов, мы получим число: $(1 \times k) \centerdot (k \times 1) = 1 \times 1$, что и следовало ожидать. Это число и есть приближение "$i$-го" целевого показателя. Но нам-то нужно приближение не одного значения целевого показателя, а всех. Для этого запишем все "$i$-ые" регрессоры в формате матрицы $X$. Полученная матрица имеет размерность $(n \times k)$:

$X=\begin{pmatrix} x_{00} & x_{01} & ... & x_{0k} \\ x_{10} & x_{11} & ... & x_{1k} \\ ... & ... & ... & ... \\ x_{n0} & x_{n1} & ... & x_{nk} \end{pmatrix} \qquad$



Теперь уравнение линейной регрессии примет вид:

$f(w,X) = X \vec{w}$



Обозначим значения целевых показателей (все $y_i$) за вектор $\vec{y}$ размерностью $(n \times 1)$:

$\vec{y} = \begin{pmatrix} y_{0} \\ y_{1} \\ ... \\ y_{n} \end{pmatrix} \qquad$



Теперь мы можем записать в матричном формате уравнение оценки качества линейной модели:

$Err = (X \vec{w} - \vec{y})^2\rightarrow min$



Собственно, из этой формулы далее получают известную нам формулу $X^T X w = X^T y$

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

Матричные преобразования


Раскроем скобки


$(X \vec{w} - \vec{y})^2 = (X \vec{w} - \vec{y})^T(X \vec{w} - \vec{y})$

$=(X\vec{w})^TX\vec{w} - \vec{y}^TX\vec{w} - (X\vec{w})^T\vec{y} + \vec{y}^T\vec{y}$

Подготовим уравнение для дифференцирования


Для этого проведем некоторые преобразования. В последующих расчетах нам будет удобнее, если вектор $\vec{w}^T$ будет представлен в начале каждого произведения в уравнении.

Преобразование 1


$\vec{y}^TX\vec{w} = (X\vec{w})^T\vec{y} = \vec{w}^TX^T\vec{y}$

Как это получилось? Для ответа на этот вопрос достаточно посмотреть на размеры умножаемых матриц и увидеть, что на выходе мы получаем число или иначе $const$.

Запишем размеры матричных выражений.

$\vec{y}^TX\vec{w} : (1 \times n) \centerdot (n \times k) \centerdot (k \times 1) = (1 \times 1) = const$

$(X\vec{w})^T\vec{y} : ((n \times k) \centerdot (k \times 1))^T \centerdot (n \times 1) = (1 \times n) \centerdot (n \times 1) = (1 \times 1) = const$

$\vec{w}^TX^T\vec{y} : (1 \times k) \centerdot (k \times n) \centerdot (n \times 1) = (1 \times 1) = const$

Преобразование 2


$(X\vec{w})^TX\vec{w} = \vec{w}^TX^TX\vec{w}$

Распишем аналогично преобразованию 1

$(X\vec{w})^TX\vec{w} : ((n \times k) \centerdot (k \times 1))^T \centerdot (n \times k) \centerdot (k \times 1) = (1 \times 1) = const$

$\vec{w}^TX^TX\vec{w} : (1 \times k) \centerdot (k \times n) \centerdot (n \times k) \centerdot (k \times 1) = (1 \times 1) = const$

На выходе получаем уравнение, которое нам предстоит продифференцировать:
$Err = \vec{w}^TX^TX\vec{w} - 2\vec{w}^TX^T\vec{y} + \vec{y}^T\vec{y}$

Дифференцируем функцию оценки качества модели


Продифференцируем по вектору $\vec{w}$:

$\frac{d(\vec{w}^TX^TX\vec{w} - 2\vec{w}^TX^T\vec{y} + \vec{y}^T\vec{y})}{d\vec{w}}$


$(\vec{w}^TX^TX\vec{w})' - (2\vec{w}^TX^T\vec{y})' + (\vec{y}^T\vec{y})' = 0$

$2X^TX\vec{w} - 2X^T\vec{y} + 0 = 0$

$X^TX\vec{w}=X^T\vec{y}$

Вопросов почему $(\vec{y}^T\vec{y})' = 0$ быть не должно, а вот операции по определению производных в двух других выражениях мы разберем подробнее.

Дифференцирование 1


Раскроем дифференцирование: $\frac{d(\vec{w}^TX^TX\vec{w})}{d\vec{w}} = 2X^TX\vec{w}$

Для того, чтобы определить производную от матрицы или вектора требуется посмотреть, что у них там внутри. Смотрим:

$\vec{w}^T=\begin{pmatrix} w_0 & w_1 & ... & w_k \end{pmatrix} \qquad$

$\vec{w}=\begin{pmatrix} w_0 \\ w_1 \\ ... \\ w_k \end{pmatrix} \qquad$

$X^T = \begin{pmatrix} x_{00} & x_{10} & ... & x_{n0} \\ x_{01} & x_{11} & ... & x_{n1} \\ ... & ... & ... & ... \\ x_{0k} & x_{1k} & ... & x_{nk} \end{pmatrix} \qquad$ $X = \begin{pmatrix} x_{00} & x_{01} & ... & x_{0k} \\ x_{10} & x_{11} & ... & x_{1k} \\ ... & ... & ... & ... \\ x_{n0} & x_{n1} & ... & x_{nk} \end{pmatrix} \qquad$

Обозначим произведение матриц $X^TX$ через матрицу $A$. Матрица $A$ квадратная и более того, она симметричная. Эти свойства нам пригодятся далее, запомним их. Матрица $A$ имеет размерность $(k \times k)$:

$A = \begin{pmatrix} a_{00} & a_{01} & ... & a_{0k} \\ a_{10} & a_{11} & ... & a_{1k} \\ ... & ... & ... & ... \\ a_{k0} & a_{k1} & ... & a_{kk} \end{pmatrix} \qquad$

Теперь наша задача правильно перемножить вектора на матрицу и не получить «дважды два пять», поэтому сосредоточимся и будем предельно внимательны.

$\vec{w}^TA\vec{w}=\begin{pmatrix} w_0 & w_1 & ... & w_k \end{pmatrix} \qquad \times \begin{pmatrix} a_{00} & a_{01} & ... & a_{0k} \\ a_{10} & a_{11} & ... & a_{1k} \\ ... & ... & ... & ... \\ a_{k0} & a_{k1} & ... & a_{kk} \end{pmatrix} \qquad \times \begin{pmatrix} w_0 \\ w_1 \\ ... \\ w_k \end{pmatrix} \qquad = $

$= \begin{pmatrix} w_0a_{00}+w_1a_{10}+...+w_ka_{k0} & ... & w_0a_{0k}+w_1a_{1k}+...+w_ka_{kk} \end{pmatrix} \times \begin{pmatrix} w_0 \\ w_1 \\ ... \\ w_k \end{pmatrix} \qquad =$

$= \begin{pmatrix} (w_0a_{00}+w_1a_{10} +...+w_ka_{k0})w_0 \mkern 10mu + \mkern 10mu ... \mkern 10mu + \mkern 10mu (w_0a_{0k}+w_1a_{1k}+...+w_ka_{kk})w_k \end{pmatrix} =$

$= w_0^2a_{00}+w_1a_{10}w_0+w_ka_{k0}w_0 \mkern 10mu + \mkern 10mu ... \mkern 10mu + \mkern 10mu w_0a_{0k}w_k+w_1a_{1k}w_k+...+w_k^2a_{kk}$

Однако, замысловатое выражение у нас получилось! На самом деле мы получили число — скаляр. И теперь, уже по-настоящему, переходим к дифференцированию. Необходимо найти производную полученного выражения по каждому коэффициенту $w_0 w_1 ... w_k$ и получить на выходе вектор размерности $(k \times 1)$. На всякий случай распишу процедуры по действиям:

1) продифференцируем по $w_o$, получим: $2w_0a_{00}+w_1a_{10}+w_2a_{20} + ... + w_ka_{k0}+a_{01}w_1+a_{02}w_2+...+a_{0k}w_{k}$

2) продифференцируем по $w_1$, получим: $w_0a_{01}+2w_1a_{11}+w_2a_{21} + ... + w_ka_{k1}+a_{10}w_0+a_{12}w_2+...+a_{1k}w_{k}$

3) продифференцируем по $w_k$, получим: $w_0a_{0k}+w_1a_{1k}+w_2a_{2k} + ... + w_{(k-1)}a_{(k-1)k}+a_{k0}w_0+a_{k1}w_1+a_{k2}w_2+...+2w_ka_{kk}$

На выходе — обещанный вектор размером $(k \times 1)$:

$\begin{pmatrix} 2w_0a_{00}+w_1a_{10}+w_2a_{20} + ... + w_ka_{k0}+a_{01}w_1+a_{02}w_2+...+a_{0k}w_{k} \\ w_0a_{01}+2w_1a_{11}+w_2a_{21} + ... + w_ka_{k1}+a_{10}w_0+a_{12}w_2+...+a_{1k}w_{k} \\ ... \\ ... \\ ... \\ w_0a_{0k}+w_1a_{1k}+w_2a_{2k} + ... + w_{(k-1)}a_{(k-1)k}+a_{k0}w_0+a_{k1}w_1+a_{k2}w_2+...+2w_ka_{kk} \end{pmatrix}$



Если присмотреться к вектору повнимательнее, то можно заметить, что левые и соответствующие правые элементы вектора можно сгруппировать таким образом, что в итоге из представленного вектора можно выделить вектор $\vec{w}$ размера $(k \times 1)$. Например, $w_1a_{10}$ (левый элемент верхней строчки вектора) $+ a_{01}w_1$ (правый элемент верхней строчки вектора) можно представить как $w_1(a_{10}+a_{01})$, а $w_2a_{20}+a_{02}w_2$ — как $w_2(a_{20}+a_{02})$ и т.д. по каждой строчке. Сгруппируем:

$\begin{pmatrix} 2w_0a_{00}+w_1(a_{10}+a_{01})+w_2(a_{20}+a_{02})+...+w_k(a_{k0}+a_{0k}) \\ w_0(a_{01}+a_{10})+2w_1a_{11}+w_2(a_{21}+a_{12})+...+w_k(a_{k1}+a_{1k}) \\ ... \\ ... \\ ... \\ w_0(a_{0k}+a_{k0})+w_1(a_{1k}+a_{k1})+w_2(a_{2k}+a_{k2})+...+2w_ka_{kk} \end{pmatrix}$



Вынесем вектор $\vec{w}$ и на выходе получим:

$\begin{pmatrix} 2a_{00} & a_{10}+a_{01} & a_{20}+a_{02} &...& a_{k0}+a_{0k} \\ a_{01}+a_{10} & 2a_{11} & a_{21}+a_{12} &...& a_{k1}+a_{1k} \\ ...&...&...&...&... \\ ...&...&...&...&... \\ ...&...&...&...&... \\ a_{0k}+a_{k0} & a_{1k}+a_{k1} & a_{2k}+a_{k2} &...& 2a_{kk} \end{pmatrix} \times \begin{pmatrix} w_0 \\ w_1 \\ ... \\ ... \\ ... \\ w_k \end{pmatrix} \qquad$



Теперь, присмотримся к получившейся матрице. Матрица представляет собой сумму двух матриц $A+A^T$:

$\begin{pmatrix} a_{00} & a_{01} & a_{02} &...& a_{0k} \\ a_{10} & a_{11} & a_{12} &...& a_{1k} \\ ...&...&...&...&... \\ a_{k0} & a_{k1} & a_{k2} &...& a_{kk} \end{pmatrix} + \begin{pmatrix} a_{00} & a_{10} & a_{20} &...& a_{k0} \\ a_{01} & a_{11} & a_{21} &...& a_{k1} \\ ...&...&...&...&... \\ a_{0k} & a_{1k} & a_{2k} &...& a_{kk} \end{pmatrix} \qquad$



Вспомним, что несколько ранее, мы отметили одно важное свойство матрицы $A$ — она симметричная. Исходя из этого свойства, мы можем с уверенностью заявить, что выражение $A + A^T$ равняется $2A$. Это легко проверить, раскрыв поэлементно произведение матриц $X^TX$. Мы не будем делать этого здесь, желающие могут провести проверку самостоятельно.

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

$(A+A^T) \times \begin{pmatrix} w_0 \\ w_1 \\ ... \\ w_k \end{pmatrix} \qquad = 2A \vec{w} = 2X^TX\vec{w}$



Итак, с первым дифференцированием мы справились. Переходим ко второму выражению.

Дифференцирование 2


$\frac{d(2\vec{w}^TX^T\vec{y})}{d\vec{w}} = 2X^T\vec{y}$

Пойдем по протоптанной дорожке. Она будет намного короче предыдущей, так что не уходите далеко от экрана.

Раскроем поэлементно вектора и матрицу:

$\vec{w}^T=\begin{pmatrix} w_0 & w_1 & ... & w_k \end{pmatrix} \qquad$

$X^T = \begin{pmatrix} x_{00} & x_{10} & ... & x_{n0} \\ x_{01} & x_{11} & ... & x_{n1} \\ ... & ... & ... & ... \\ x_{0k} & x_{1k} & ... & x_{nk} \end{pmatrix} \qquad$

$\vec{y}=\begin{pmatrix} y_0 \\ y_1 \\ ... \\ y_n \end{pmatrix} \qquad$

На время уберем из расчетов двойку — она большой роли не играет, потом вернем ее на место. Перемножим вектора на матрицу. В первую очередь умножим матрицу $X^T$ на вектор $\vec{y}$, здесь у нас нет никаких ограничений. Получим вектор размера $(k \times 1)$:

$\begin{pmatrix} x_{00}y_0+x_{10}y_1+...+x_{n0}y_n \\ x_{01}y_0+x_{11}y_1+...+x_{n1}y_n \\ ... \\ x_{0k}y_0+x_{1k}y_1+...+x_{nk}y_n \end{pmatrix} \qquad$



Выполним следующее действие — умножим вектор $\vec{w}$ на полученный вектор. На выходе нас будет ждать число:

$\begin{pmatrix} w_0(x_{00}y_0+x_{10}y_1+...+x_{n0}y_n)+ w_1(x_{01}y_0+x_{11}y_1+...+x_{n1}y_n) \mkern 10mu + \mkern 10mu ... \mkern 10mu + \mkern 10mu w_k(x_{0k}y_0+x_{1k}y_1+...+x_{nk}y_n) \end{pmatrix} \qquad$



Его то мы и продифференцируем. На выходе получим вектор размерности $(k \times 1)$:

$\begin{pmatrix} x_{00}y_0+x_{10}y_1+...+x_{n0}y_n \\ x_{01}y_0+x_{11}y_1+...+x_{n1}y_n \\ ... \\ x_{0k}y_0+x_{1k}y_1+...+x_{nk}y_n \end{pmatrix} \qquad$



Что-то напоминает? Все верно! Это произведение матрицы $X^T$ на вектор $\vec{y}$.

Таким образом, второе дифференцирование успешно завершено.

Вместо заключения


Теперь мы знаем, как получилось равенство $X^T X \vec{w} = X^T \vec{y}$.

Напоследок опишем быстрый путь преобразований основных формул.

Оценим качество модели в соответствии с методом наименьших квадратов:
$\sum\limits_{i=1}^n(y_i-f(x_i))^2 \mkern 20mu = \mkern 20mu \sum\limits_{i=1}^n(y_i-\vec{x_i}^T \vec{w})^2=$

$= (X \vec{w} - \vec{y})^2 \mkern 20mu = \mkern 20mu (X \vec{w} - \vec{y})^T(X \vec{w} - \vec{y}) \mkern 20mu = \mkern 20mu \vec{w}^TX^TX\vec{w} - 2\vec{w}^TX^T\vec{y} + \vec{y}^T\vec{y}$

Дифференцируем полученное выражение:
$\frac{d(\vec{w}^TX^TX\vec{w} - 2\vec{w}^TX^T\vec{y} + \vec{y}^T\vec{y})}{d\vec{w}}=$ $2X^TX\vec{w} - 2X^T\vec{y} = 0$

$X^TX\vec{w}=X^T\vec{y}$

$\leftarrow $ Предыдущая работа автора — «Решаем уравнение простой линейной регрессии»
$\rightarrow$ Следующая работа автора — «Пережевывая логистическую регрессию»

Литература


Интернет источники:

1) habr.com/ru/post/278513
2) habr.com/ru/company/ods/blog/322076
3) habr.com/ru/post/307004
4) nabatchikov.com/blog/view/matrix_der

Учебники, сборники задач:

1) Конспект лекций по высшей математике: полный курс / Д.Т. Письменный – 4-е изд. – М.: Айрис-пресс, 2006
2) Прикладной регрессионный анализ / Н. Дрейпер, Г. Смит – 2-е изд. – М.: Финансы и статистика, 1986 (перевод с английского)
3) Задачи на решение матричных уравнений:
function-x.ru/matrix_equations.html
mathprofi.ru/deistviya_s_matricami.html
Tags:
Hubs:
Total votes 11: ↑11 and ↓0+11
Comments6
18
Karma
0
Rating
Александр Петренко @AlexanderPetrenko

Машинное обучение

Send message
Instagram

Articles

Your account

Sections

Information

Services

Support
© 2006–2024, Habr