...

воскресенье, 24 ноября 2019 г.

Методы наименьших квадратов: текст, написанный программистом для программистов

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

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

Текущий план лекций:

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

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


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

Различные интерпретации матриц


Ответ очень простой. Матрица это просто шкафчик, в котором хранятся хреновины. Каждая хреновина лежит в своей ячейке, ячейки группируются рядами в строки и столбцы. В нашем конкретном случае хреновинами будут обычные вещественные числа; для программиста проще всего представлять матрицу $A$ как нечто навроде 
float A[m][n];

Зачем же такие хранилища нужны? Что они описывают? Может быть, я вас расстрою, но сама по себе матрица не описывает ничего, она хранит. Например, в ней можно хранить коэффициенты всяких функций. Давайте на секунду забудем про матрицы, и представим, что у нас есть число $a$. Что оно означает? Да чёрт его знает… Например, это может быть коэффициент внутри функции, которая на вход берёт одно число, и на выход даёт другое число:

$ f(x) : \mathbb R \rightarrow \mathbb R $

Одну версию такой функции математик мог бы записать как

$ f(x) = ax $


Ну а в мире программистов она бы выглядела следующим образом: 
float f(float x) {
    return a*x;
}

С другой стороны, а почему такая функция, а не совсем другая? Давайте возьмём другую!

$ f(x) = ax^2 $


Раз уж я начал про программистов, я обязан записать её код :) 
float f(float x) {
    return x*a*x;
}

Одна из этих функций линейная, а вторая квадратичная. Какая из них правильная? Да никакая, само по себе число $a$ не определяет этого, оно просто хранит какое-то значение! Какую вам надо функцию, такую и стройте.

То же самое происходит и с матрицами, это хранилища, нужные на случай, когда одиночных чисел (скаляров) не хватает, своего рода расширение чисел. Над матрицами, ровно как и над числами, определены операции сложения и умножения.

Давайте представим, что у нас есть матрица $A$, например, размера 2x2:

$ A=\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22}\end{bmatrix} $

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

$ f(x) : \mathbb R^2 \rightarrow \mathbb R^2, \quad f(x) = Ax $

vector<float> f(vector<float> x) {
    return vector<float>{a11*x[0] + a12*x[1],  a21*x[0] + a22*x[1]};
}

Эта функция преобразует двумерный вектор в двумерный вектор. Графически это удобно представлять как преобразование изображения: на вход даём изображение, а на выходе получаем его растянутую и/или повёрнутую (возможно даже зеркально отражённую!) версию. Очень скоро я приведу картинку с различными примерами такой интерпретации матриц.

А можно матрицу $A$ представлять себе как функцию, которая двумерный вектор преобразует в скаляр:

$ f(x) : \mathbb R^2 \rightarrow \mathbb R, \quad f(x) = x^\top A x = \sum\limits_i\sum\limits_j a_{ij}x_i x_j $

Обратите внимание, что с векторами возведение в степень не очень-то определено, поэтому я не могу написать $x^2$, как писал в случае с обычными числами. Очень рекомендую тем, кто не привык с лёгкостью жонглировать матричными умножениями, ещё раз вспомнить правило умножения матриц, и проверить, что выражение $x^\top A x$ вообще разрешено и действительно даёт скаляр на выходе. Для этого можно, например, явно поставить скобки $x^\top A x = (x^\top A) x$. Напоминаю, что у нас $x$ — вектор размерности 2 (сохранённый в матрице размерности 2x1), выпишем явно все размерности:

$ \underbrace{\underbrace{\left(\underbrace{x^\top}_{1\times 2} \times \underbrace{A}_{2\times 2}\right)}_{1\times 2} \times \underbrace{x}_{2\times 1}}_{1 \times 1} $

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

float f(vector<float> x) {
    return x[0]*a11*x[0] + x[0]*a12*x[1] + x[1]*a21*x[0] + x[1]*a22*x[1];
}

Что такое положительное число?


Теперь я задам очень глупый вопрос: что такое положительное число? У нас есть отличный инструмент, называется предикат больше $>$. Но не торопитесь отвечать, что число $a$ положительно тогда и только тогда, когда $a>0$, это было бы слишком просто. Давайте определим положительность следущим образом:

Определение 1


Число $a$ положительно тогда и только тогда, когда для всех ненулевых вещественных $x\in\mathbb R,\ x\neq 0$ выполняется условие $ax^2>0$.

Выглядит довольно мудрёно, но зато отлично применяется к матрицам:

Определение 2


Квадратная матрица $A$ называется положительно определённой, если для любых ненулевых $x$ выполняется условие $x^\top A x > 0$, то есть, соответствующая квадратичная форма строго положительна везде, кроме начала координат.

Для чего мне нужна положительность? Как я уже упоминал в начале статьи, моим основным инструментом будет поиск минимумов квадратичных функций. А ведь для того, чтобы минимум искать, было бы недурно, если бы он существовал! Например, у функции $f(x) = - x^2$ минимума очевидно не существует, посколькуо число -1 не является положительным, и $f(x)$ бесконечно убывает с ростом $x$: ветви параболы $f(x)$ смотрят вниз. Положительно определённые матрицы хороши тем, что соответствующие квадратичные формы образуют параболоид, имеющие (единственный) минимум. Следующая иллюстрация показывает семь различных примеров матриц $2\times 2$, как положительно определённых и симметричных, так и не очень. Верхний ряд: интерпретация матриц как функций $f(x):\mathbb R^2 \rightarrow \mathbb R^2$. Средний ряд: интерпретация матриц как функций $f(x):\mathbb R^2 \rightarrow \mathbb R$.

Таким образом, я буду работать с положительно определёнными матрицами, которые являются обобщением положительных чисел. Более того, конкретно в моём случае матрицы будут ещё и симметричными! Любопытно, что довольно часто, когда люди говорят про положительную определённость, они подразумевают ещё и симметричность, что может быть косвенно объяснено следующим (необязательным для понимания последующего текста) наблюдением.

Лирическое отступление о симметричности матриц квадратичных форм


Давайте рассмотрим квадратичную форму $x^\top M x$ для произвольной матрицы $M$; затем к $M$ добавим и сразу же отнимем половину её транспонированной версии:

$ M = \underbrace{\frac{1}{2} (M+M^\top)}_{:=M_s} + \underbrace{\frac{1}{2} (M-M^\top)}_{:=M_a} = M_s + M_a $

Матрица $M_s$ симметрична: $M_s^\top = M_s$; матрица $M_a$ антисимметрична: $M_a^\top=-M_a$. Примечательным фактом является то, что для любой антисимметричной матрицы соответствующая квадратичная форма тождественно равна нулю. Это вытекает из следующего наблюдения:

$ q = x^\top M_a x = (x^\top M_a^\top x)^\top = - (x^\top M_a x)^\top = -q $


То есть, квадратичная форма $x^\top M_a x$ одновременно равна $q$ и $-q$, что возможно только в том случае, когда $q\equiv 0$. Из этого факта вытекает то, что для произвольной матрицы $M$ соответствующая квадратичная форма $x^\top M x$ может быть выражена и через симметричную матрицу $M_s$:

$ x^\top M x = x^\top (M_s + M_a) x = x^\top M_s x + x^\top M_a x = x^\top M_s x. $

Ищем минимум квадратичной функции


Вернёмся ненадолго в одномерный мир; я хочу найти минимум функции $f(x) = ax^2 - 2bx$. Число $a$ положительно, поэтому минимум существует; чтобы его найти, приравняем нулю соответствующую производную: $\frac{d}{dx}f(x) = 0$. Продифференцировать одномерную квадратичную функцию труда не составляет: $\frac{d}{dx}f(x) = 2ax - 2b = 0$; поэтому наша проблема сводится к решению уравнения $ax-b=0$, откуда путём страшных усилий получам решение $x^* = b/a$. Следующий рисунок иллюстрирует эквивалентность двух проблем: решение $x^*$ уравнения $ax-b=0$ совпадает с решением уравнения $\mathrm{argmin}_x(ax^2 - 2bx)$.

Я клоню к тому, что нашей глобальной целью является минимизация квадратичных функций, мы же про наименьшие квадраты говорим. Но при этом что мы действительно умеем делать хорошо, так это решать линейные уравнения (и системы линейных уравнений). И очень здорово, что одно эквивалентно другому!
Единственное, что нам осталось, так это убедиться, что эта эквивалентность в одномерном мире распространяется и на случай $n$ переменных. Чтобы это проверить, для начала докажем три теоремы.

Три теоремы, или дифференцируем матричные выражения


Первая теорема гласит о том, что матрицы $1\times 1$ инварианты относительно транспонирования:

Теорема 1


$x\in \mathbb R \Rightarrow x^\top = x$

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

Вторая теорема позволяет дифференцировать линейные функции. В случае вещественной функции одной переменной мы знаем, что $\frac{d}{dx}(bx) = b$, но что происходит в случае вещественной функции $n$ переменных?

Теорема 2


$\nabla b^\top x = \nabla x^\top b = b$

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

$\nabla(b^\top x) = \begin{bmatrix}\frac{\partial (b^\top x)}{\partial x_1} \\ \vdots \\ \frac{\partial (b^\top x)}{\partial x_n} \end{bmatrix} = \begin{bmatrix}\frac{\partial (b_1 x_1 + \dots + b_n x_n)}{\partial x_1} \\ \vdots \\ \frac{\partial (b_1 x_1 + \dots + b_n x_n)}{\partial x_n} \end{bmatrix} = \begin{bmatrix}b_1 \\ \vdots \\ b_n \end{bmatrix} = b$


Применив первую теорему $b^\top x = x^\top b$, мы завершаем доказательство.

Теперь переходим к квадратичным формам. Мы знаем, что в случае вещественной функции одной переменной $\frac{d}{dx}(ax^2) = 2ax$, а что будет в случае квадратичной формы?

Теорема 3


$\nabla (x^\top A x) = (A+A^\top)x$

Кстати, обратите внимание, что если матрица $A$ симметрична, то теорема гласит, что $\nabla (x^\top A x) = 2Ax$.

Это доказательство тоже прямолинейно, просто запишем квадратичную форму как двойную сумму:

$x^\top A x = \sum\limits_i\sum\limits_j a_{ij} x_i x_j$


А затем посмотрим, чему равна частная производная этой двойной суммы по переменной $x_i$:
\begin{align*}
\frac{\partial (x^\top A x)}{\partial x_i}
&= \frac{\partial}{\partial x_i} \left(\sum\limits_{k_1}\sum\limits_{k_2} a_{k_1 k_2} x_{k_1} x_{k_2}\right) = \\
&= \frac{\partial}{\partial x_i} \left(
\underbrace{\sum\limits_{k_1\neq i}\sum\limits_{k_2\neq i} a_{ik_2}x_{k_1} x_{k_2}}_{k_1 \neq i, k_2 \neq i}+\underbrace{\sum\limits_{k_2\neq i} a_{ik_2}x_i x_{k_2}}_{k_1 = i, k_2\neq i}+
\underbrace{\sum\limits_{k_1\neq i} a_{k_1 i} x_{k_1} x_i}_{k_1 \neq i, k_2 = i}+
\underbrace{a_{ii}x_i^2}_{k_1 = i, k_2 = i}\right) = \\
& = \sum\limits_{k_2\neq i} a_{ik_2}x_{k_2} + \sum\limits_{k_1\neq i} a_{k_1 i} x_{k_1} + 2 a_{ii} x_i = \\
& = \sum\limits_{k_2} a_{ik_2}x_{k_2} + \sum\limits_{k_1} a_{k_1 i} x_{k_1} = \\
& = \sum\limits_{j} (a_{ij} + a_{ji}) x_j \\
\end{align*}
Я разбил двойную сумму на четыре случая, которые выделены фигурными скобками. Каждый из этих четырёх случаев дифференцируется тривиально. Осталось сделать последнее действие, собрать частные производные в вектор градиента:

$\nabla(x^\top A x) = \begin{bmatrix}\frac{\partial (x^\top Ax)}{\partial x_1} \\ \vdots \\ \frac{\partial (x^\top A x)}{\partial x_n} \end{bmatrix} = \begin{bmatrix}\sum\limits_{j} (a_{1j} + a_{j1}) x_j \\ \vdots \\ \sum\limits_{j} (a_{nj} + a_{jn}) x_j \end{bmatrix} = (A+A^\top)x $

Минимум квадратичной функции и решение линейной системы


Давайте не забывать направление движения. Мы видели, что при положительном числе $a$ в случае вещественных функций одной перменной решить уравнение $ax=b$ или минимизировать квадратичную функцию $\mathrm{argmin}_x(ax^2 - 2bx)$ — это одно и то же.

Мы хотим показать соответствующую связь в случае симметричной положительно определённой матрицы $A$.
Итак, мы хотим найти минимум квадратичной функции

$\mathrm{argmin}_{x\in\mathbb R^n} (x^\top A x - 2b^\top x).$


Ровно как и в школе, приравняем нулю производную:

$\nabla (x^\top A x - 2b^\top x) = \begin{bmatrix}0 \\ \vdots \\ 0 \end{bmatrix}.$


Оператор Гамильтона линеен, поэтому мы можем переписать наше уравнение в следующем виде:

$\nabla (x^\top A x) - 2\nabla(b^\top x) = \begin{bmatrix}0 \\ \vdots \\ 0 \end{bmatrix}$


Теперь применим вторую и третью теоремы о дифференцировании:

$(A+A^\top)x - 2b = \begin{bmatrix}0 \\ \vdots \\ 0 \end{bmatrix}$


Вспомним, что $A$ симметрична, и поделим уравнение на два, получаем нужную нам систему линейных уравнений:

$Ax = b$


Ура-ура, перейдя от одной переменной к многим, мы не потеряли ровным счётом ничего, и можем эффективно минимизировать квадратичные функции!
Наконец мы можем перейти к основному содержанию этой лекции. Представьте, что у нас есть две точки на плоскости $(x_1, y_1)$ и $(x_2, y_2)$, и мы хотим найти уравнение прямой, проходящей через эти две точки. Уравнение прямой можно записать в виде $y = \alpha x + \beta$, то есть, наша задача это найти коэффициенты $\alpha$ и $\beta$. Это упражнение для седьмого-восьмого класса школы, запишем систему уравнений:

$ \left\{ \begin{split} \alpha x_1 + \beta &= y_1\\ \alpha x_2 + \beta &= y_2\\ \end{split} \right. $

У нас два уравнения с двумя неизвестными ($\alpha$ и $\beta$), остальное известно.
В общем случае решение существует и единственно. Для удобства перепишем ту же самую систему в матричном виде:

$ \underbrace{\begin{bmatrix}x_1 & 1 \\ x_2 & 1 \end{bmatrix}}_{:=A} \underbrace{\begin{bmatrix} \alpha \\ \beta \end{bmatrix}}_{:=x} = \underbrace{\begin{bmatrix} y_1 \\ y_2 \end{bmatrix}}_{:=b} $

Получим уравнение типа $Ax = b$, которое тривиально решается $x^* = A^{-1}b$.

А теперь представим, что у нас есть три точки, через которые нужно провести прямую:

$ \left\{ \begin{split} \alpha x_1 + \beta &= y_1\\ \alpha x_2 + \beta &= y_2\\ \alpha x_3 + \beta &= y_3\\ \end{split} \right. $

В матричном виде эта система запишется следующим образом:

$ \underbrace{\begin{bmatrix}x_1 & 1 \\ x_2 & 1 \\x_3 & 1 \end{bmatrix} }_{:= A\,(3\times 2)} \underbrace{\begin{bmatrix} \alpha \\ \beta \end{bmatrix}}_{:=x\,(2\times 1)} = \underbrace{\begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}}_{:=b\, (3\times 1)} $


Теперь у нас матрица $A$ прямоугольная, и у неё просто не существует обратной! Это совершенно нормально, так как у нас всего две переменных и три уравнения, и в общем случае эта система не имеет решения. Это соверешенно нормальная ситуация в реальном мире, когда точки — это данные зашумлённых измерений, и нам нужно найти параметры $\alpha$ и $\beta$, которые наилучшим образом приближают данные измерений. Мы этот пример уже рассматривали в первой лекции, когда калибровали безмен. Но тогда у нас решение было чисто алгебраическим и очень громоздким. Давайте попробуем более интуитивный способ.

Нашу систему можно записать в следующем виде:

$ \alpha \underbrace{\begin{bmatrix}x_1 \\ x_2 \\x_3 \end{bmatrix} }_{:=\vec{i}} +\beta \underbrace{\begin{bmatrix}1 \\ 1 \\1 \end{bmatrix} }_{:=\vec{j}} = \begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix} $


Или, более кратко,

$ \alpha \vec{i} + \beta\vec{j} = \vec{b}. $


Векторы $\vec i$, $\vec j$ и $\vec b$ известны, нужно найти скаляры $\alpha$ и $\beta$.
Очевидно, что линейная комбинация $\alpha \vec{i} + \beta\vec{j}$ порождает плоскость, и если вектор $\vec b$ в этой плоскости не лежит, то точного решения не существует; однако же мы ищем приближённое решение.

Давайте определим ошибку решения как $\vec{e}(\alpha, \beta) := \alpha \vec{i} + \beta\vec{j} - \vec b$. Нашей зачей является найти такие значения параметров $\alpha$ и $\beta$, которые минимизируют длину вектора $\vec{e}(\alpha, \beta)$. Иначе говоря, проблема записывается как $\mathrm{argmin}_{\alpha, \beta} \|\vec{e}(\alpha, \beta)\|$.
Вот иллюстрация:.

Имея заданные векторы $\vec i$, $\vec j$ и $\vec b$, мы стараемся минимизировать длину вектора ошибки $\vec e$. Очевидно, что его длина минимизируется, если он перпендикулярен плоскости, натянутой на векторы $\vec i$ и $\vec j$.

Но постойте, где же наименьшие квадраты? Сейчас будут. Функция извлечения корня $\sqrt{\cdot}$ монотонна, поэтому $\mathrm{argmin}_{\alpha, \beta} \|\vec{e}(\alpha, \beta)\|$ = $\mathrm{argmin}_{\alpha, \beta} \|\vec{e}(\alpha, \beta)\|^2$!

Вполне очевидно, что длина вектора ошибки минимизируется, если он перпендикулярен плоскости, натянутой на векторы $\vec i$ и $\vec j$, что можно записать, приравняв нулю соответствующие скалярные произведения:

$ \left\{ \begin{split}\vec{i}^\top \vec{e}(\alpha, \beta) &= 0\\ \vec{j}^\top \vec{e}(\alpha, \beta) &= 0 \end{split} \right. $

В матричном виде эту же самую систему можно записать как

$ \begin{bmatrix}x_1 & x_2 & x_3 \\ 1 & 1 & 1 \end{bmatrix} \left(\alpha \begin{bmatrix}x_1 \\ x_2 \\x_3 \end{bmatrix} +\beta \begin{bmatrix}1 \\ 1 \\1 \end{bmatrix} - \begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix}\right) = \begin{bmatrix}0\\0\end{bmatrix} $


или же

$ \begin{bmatrix}x_1 & x_2 & x_3 \\ 1 & 1 & 1 \end{bmatrix} \left( \begin{bmatrix}x_1 & 1 \\ x_2 & 1 \\x_3 & 1 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}- \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} \right) = \begin{bmatrix}0\\0\end{bmatrix} $


Но мы на этом не остановимся, так как запись можно ещё больше сократить:

$ A^\top (Ax - b)= \begin{bmatrix}0\\0\end{bmatrix} $


И самая последняя трансформация:

$ A^\top Ax = A^\top b. $


Давайте посчитаем размерности. Матрица $A$ имеет размер $3\times 2$, поэтому $A^\top A$ имеет размер $2\times 2$. Матрица $b$ имеет размер $3\times 1$, но вектор $A^\top b$ имеет размер $2\times 1$. То есть, в общем случае матрица $A^\top A$ обратима! Более того, $A^\top A$ имеет ещё пару приятных свойств.

Теорема 4


$A^\top A$ симметрична.

Это совсем нетрудно показать:

$ (A^\top A)^\top = A^\top (A^\top)^\top = A^\top A. $

Теорема 5


$A^\top A$ положительно полуопределена: $\forall x\in \mathbb R^n\quad x^\top A^\top A x \geq 0.$

Это следует из того факта, что $x^\top A^\top A x = (A x)^\top A x > 0$.

Кроме того, $A^\top A$ положительно определена в том случае, если $A$ имеет линейно независимые столбцы (ранг $A$ равен количеству переменных в системе).

Итого, для системы с двумя неизвестными мы доказали, что минимизация квадратичной функции $\mathrm{argmin}_{\alpha, \beta} \|\vec{e}(\alpha, \beta)\|^2$ эквивалентна решению системы линейных уравнений $A^\top A x = A^\top b$. Разумеется, всё это рассуждение применимо и к любому другому количеству переменных, но давайте ещё раз компактно запишем всё вместе простым алгебраическим подсчётом. Мы начнём с проблемы наименьших квадратов, придём к минимизации квадратичной функции знакомого нам вида, и из этого сделаем вывод об эквивалентности решению системы линейных уравнений. Итак, поехали:

$ \begin{split} \mathrm{argmin} \| Ax - b \|^2 &= \mathrm{argmin} (Ax-b)^\top (Ax-b) =\\ & = \mathrm{argmin}(x^\top A^\top - b^\top)(Ax-b) = \\ & = \mathrm{argmin}(x^\top A^\top A x - b^\top Ax - x^\top A^\top b + \underbrace{b^\top b}_{\rm const})=\\ & = \mathrm{argmin}(x^\top A^\top A x - 2b^\top Ax) = \\ & = \mathrm{argmin}(x^\top \underbrace{(A^\top A)}_{:=A'} x - 2\underbrace{(A^\top b)}_{:=b'}\phantom{}^\top x) \end{split} $


Таким образом, проблема наименьших квадратов $\mathrm{argmin} \| Ax - b \|^2$ эквивалентна минимизации квадратичной функции $\mathrm{argmin} (x^\top A' x - 2b'^\top x)$ с (в общем случае) симметричной положительно определённой матрицей $A'$, что, в свою очередь, эквивалентно решению системы линейных уравнений $A'x = b'$. Уфф. Теория закончилась.

Пример первый, одномерный


Давайте вспомним один пример из предыдущей статьи: 
# initialize the data array
x = [0.] * 32
x[0] = x[31] = 1.
x[18] = 2.

# smooth the array values other than at indices 0,18,31
for iter in range(1000):
    for i in range(len(x)):
        if i in [0,18,31]: continue
        x[i] = (x[i-1] + x[i+1])/2.

У нас есть обычный массив x из 32 элементов, инициализированный следующим образом:

А затем мы тысячу раз выполним следующую процедуру: для каждой ячейки x[i] мы запишем в неё среднее значение соседних ячеек: x[i] = ( x[i-1] + x[i+1] )/2. Единственное исключение: мы не будем переписывать элементы массива под номерами 0, 18 и 31. Вот так выглядит эволюция массива на первых ста пятидесяти итерациях:

Как мы выяснили в предыдущей статье, оказывается, наш питоновский мегакод на четыре строчки решает следующиую систему линейных уравнений методом Гаусса-Зейделя:

Выяснить-то выяснили, но откуда вообще взялась эта система? Что за магия? Давайте отвлечёмся на секунду, и попробуем решить немного другую систему. У меня по-прежнему есть 32 переменных, но теперь я хочу, чтобы они все-все-все были равны одна другой. Это нетрудно записать, достаточно записать систему уравнений, которая гласит, что все соседние элементы попарно равны:

Вышеприведённый питоновский код в принципе не трогает значения трёх переменных: x0, x18, x31, поэтому давайте их перенесём в правую часть системы уравнений, т.е., исключим из множества неизвестных:

Я фиксирую начальное значение x0=1, первое уравнение говорит, что x1 = x0 = 1, второе уравнение говорит, что x2=x1=x0=1 и так далее, покуда мы не дойдём до уравнения 1 = x0 = x17 = x18 = 2. Ой. А ведь эта система не имеет решения :((

А и пофиг, мы же программисты! Давайте звать на помощь наименьшие квадраты! Наша нерешаемая система может быть записана в матричном виде $Ax=b$; чтобы решить нашу систему приближённо, достаточно решить систему $A^\top Ax=A^\top b$. А тут оказывается, что это ровно, один-в-один та самая система уравнений из предыдущей статьи!

Интуитивно можно представлять себе процесс создания матрицы системы следующим образом: мы соединяем пружинками значения, равенства которых хотим добиться. Например, в нашем случае мы хотим, чтобы $x_i$ был равен $x_{i+1}$, поэтому мы добавляем уравнение $x_i - x_{i+1}=0$, вешая между этими значениями пружинку, которая будет их стягивать. Но поскольку значения x0, x18 и x31 жёстко зафиксированы, свободным значениям не остаётся ничего другого, как равномерно растянуть пружинки, произведя (в данном случае) кусочно-линейное решение.

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

import numpy as np

n = 32 # size of the vector to produce

# fill the matrix with 2nd order finite differences
A = np.matrix(np.zeros((n-1, n)))
for i in range(0,n-1):
    A[i, i]   = -1.
    A[i, i+1] =  1.

# eliminate the constrained variables from the matrix
A = A[:,[*range(1,18)] + [*range(19,31)]]

b = np.matrix(np.zeros((n-1, 1)))
b[0,0]   =  1.
b[17,0] = -2.
b[18,0]   =  2.
b[n-2,0] = -1.

# solve the system
x = np.linalg.inv(A.transpose()*A)*A.transpose()*b
x = x.transpose().tolist()[0]

# re-introduce the constrained variables
x.insert(0, 1.)
x.insert(18, 2.)
x.append(1.)

Этот код производит список x из 32 элементов, в котором хранится решение. Мы строим матрицу $A$, строим вектор $b$, получаем вектор с решением $x = (A^\top A)^{-1}A^\top b$, а затем в этот вектор вставляем наши фиксированные значения x0=1, x18=2 и x31=1.

Этот код выполняет необходимую работу, но довольно неприятно переносить в правую часть значения фиксированных переменных. Нельзя ли схитрить? Можно! Давайте посмотрим на следующий код:

import numpy as np

n = 32 # size of the vector to produce

# fill the matrix with 2nd order finite differences
A = np.matrix(np.zeros((n-1+3, n)))
for i in range(0,n-1):
    A[i, i]   = -1.
    A[i, i+1] =  1.

# enforce the constraints through the quadratic penalty
A[n-1+0, 0] = 100.
A[n-1+1,18] = 100.
A[n-1+2,31] = 100.

b = np.matrix(np.zeros((n-1+3, 1)))
b[n-1+0,0] = 100.*1.
b[n-1+1,0] = 100.*2.
b[n-1+2,0] = 100.*1.

# solve the system
x = np.linalg.inv(A.transpose()*A)*A.transpose()*b
x = x.transpose().tolist()[0]

С программистской точки зрения он гораздо приятнее, вектор x получается сразу нужной размерности, да и матрицы A и b заполняются куда как проще. Но в чём же хитрость? Давайте запишем систему уравнений, которую мы решаем:

Посмотрите внимательно на последние три уравнения:
100 x0 = 100 * 1
100 x18 = 100 * 2
100 x31 = 100 * 1

Не проще ли было их записать следующим образом?
x0 = 1
x18 = 2
x31 = 1

Нет, не проще. Как я уже говорил, каждое уравнение навешивает пружинки, в данном случае пружинку между x0 и значением 1, пружинку между x18 и значением 2 и, наконец, переменная x31 тоже будет притягиваться к единице.

Но коэффициент 100 там стоит не просто так. Мы же помним, что эта система просто так не решается, мы решаем её в смысле наименьших квадратов, минимизируя соответствующую квадратичную функцию. Умножив на 100 каждый из коэффициентов последних трёх уравнений, мы вводим штраф за отклонение x0, x18 и x32 от нужных значений, в десять тысяч раз (100^2) превышающий штраф за отклонение от равенства $x_i = x_{i+1}$. Грубо говоря, пружинки на последних трёх уравнениях в десять тысяч раз более жёсткие, нежели на всех остальных. Этот метод квадратичного штрафа — очень удобный инструмент введения ограничений, гораздо проще (хотя и не без недостатков), нежели редуцирование набора переменных.

Пример второй, трёхмерный


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

Код доступен тут, но на всякий случай приведу полный листинг главного файла:

#include <vector>
#include <iostream>

#include "OpenNL_psm.h"
#include "geometry.h"
#include "model.h"

int main(int argc, char** argv) {
    if (argc<2) {
        std::cerr << "Usage: " << argv[0] << " obj/model.obj" << std::endl;
        return 1;
    }

    Model m(argv[1]);

    for (int d=0; d<3; d++) { // solve for x, y and z separately
        nlNewContext();
        nlSolverParameteri(NL_NB_VARIABLES, m.nverts());
        nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
        nlBegin(NL_SYSTEM);
        nlBegin(NL_MATRIX);

        for (int i=0; i<m.nhalfedges(); i++) { // fix the boundary vertices
            if (m.opp(i)!=-1) continue;
            int v = m.from(i);
            nlRowScaling(100.);
            nlBegin(NL_ROW);
            nlCoefficient(v,  1);
            nlRightHandSide(m.point(v)[d]);
            nlEnd(NL_ROW);
        }

        for (int i=0; i<m.nhalfedges(); i++) { // smooth the interior
            if (m.opp(i)==-1) continue;
            int v1 = m.from(i);
            int v2 = m.to(i);

            nlRowScaling(1.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlEnd(NL_ROW);
        }

        nlEnd(NL_MATRIX);
        nlEnd(NL_SYSTEM);
        nlSolve();

        for (int i=0; i<m.nverts(); i++) {
            m.point(i)[d] = nlGetVariable(i);
        }
    }

    std::cout << m;
    return 0;
}

Я написал простейший парсер wavefront .obj файлов, поэтому моделька подгружается в одну строчку. В качестве решателя я буду использовать OpenNL, это очень мощный решатель для проблем наименьших квадратов с разреженными матрицами. Обратите внимание, что размеры матриц у нас могут быть гигантскими: например, если мы хотим получить ч/б картинку размером 1000 x 1000 пикселей, то матрица $A^\top A$ у нас будет размером миллион на миллион! Однако же на практике чаще всего подавляющее большинство элементов этой матрицы нулевые, поэтому всё не так страшно, но всё же нужно использовать специальные солверы.

Итак, давайте смотреть на код. Для начала я объявляю солверу, что я хочу делать:

        nlNewContext();
        nlSolverParameteri(NL_NB_VARIABLES, m.nverts());
        nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
        nlBegin(NL_SYSTEM);
        nlBegin(NL_MATRIX);

Я хочу иметь матрицу $A^\top A$ размером (колво вершин)x(колво вершин). Внимание, мы солверу даём матрицу $A$, а уж $A^\top A$ он будет строить сам, что очень удобно! Также обратите внимание на то, что я решаю не одну систему уравнений, а три последовательно для x,y и z (я соврал, что проблема трёхмерная, она по-прежнему одномерная!)

А затем я объявляю строку за строкой нашей матрицы. Для начала я фиксирую вершины, которые находятся на краю поверхности:

        for (int i=0; i<m.nhalfedges(); i++) { // fix the boundary vertices
            if (m.opp(i)!=-1) continue;
            int v = m.from(i);
            nlRowScaling(100.);
            nlBegin(NL_ROW);
            nlCoefficient(v,  1);
            nlRightHandSide(m.point(v)[d]);
            nlEnd(NL_ROW);
        }

Можно видеть, что я использую квадратичный штраф в 100^2 для фиксирования краевых вершин.

Ну а затем для каждой пары смежных вершин я вешаю пружинку жёсткости 1 между ними:

        for (int i=0; i<m.nhalfedges(); i++) { // smooth the interior
            if (m.opp(i)==-1) continue;
            int v1 = m.from(i);
            int v2 = m.to(i);

            nlRowScaling(1.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlEnd(NL_ROW);
        }

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

Усиляем характеристические черты


Давайте превратим наше лицо в гротескную маску:

В предыдущей статье я уже показывал, как это сделать на коленке, здесь же я привожу настоящий код, использующий солвер наименьших квадратов:

    for (int d=0; d<3; d++) { // solve for x, y and z separately
        nlNewContext();
        nlSolverParameteri(NL_NB_VARIABLES, m.nverts());
        nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
        nlBegin(NL_SYSTEM);
        nlBegin(NL_MATRIX);

        for (int i=0; i<m.nhalfedges(); i++) { // fix the boundary vertices
            if (m.opp(i)!=-1) continue;
            int v = m.from(i);
            nlRowScaling(100.);
            nlBegin(NL_ROW);
            nlCoefficient(v,  1);
            nlRightHandSide(m.point(v)[d]);
            nlEnd(NL_ROW);
        }
        for (int i=0; i<m.nhalfedges(); i++) { // light attachment to the original geometry
            if (m.opp(i)==-1) continue;
            int v1 = m.from(i);
            int v2 = m.to(i);

            nlRowScaling(.2);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlRightHandSide(m.point(v1)[d] - m.point(v2)[d]);
            nlEnd(NL_ROW);
        }

        for (int v=0; v<m.nverts(); v++) { // amplify the curvature
            nlRowScaling(1.);
            nlBegin(NL_ROW);
            int nneigh = m.incident_halfedges(v).size();
            nlCoefficient(v,  nneigh);
            Vec3f curvature = m.point(v)*nneigh;
            for (int hid : m.incident_halfedges(v)) {
                nlCoefficient(m.to(hid),  -1);
                curvature = curvature - m.point(m.to(hid));
            }
            nlRightHandSide(2.5*curvature[d]);
            nlEnd(NL_ROW);
        }

        nlEnd(NL_MATRIX);
        nlEnd(NL_SYSTEM);
        nlSolve();

        for (int i=0; i<m.nverts(); i++) {
            m.point(i)[d] = nlGetVariable(i);
        }
    }

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

        for (int i=0; i<m.nhalfedges(); i++) { // light attachment to the original geometry
            if (m.opp(i)==-1) continue;
            int v1 = m.from(i);
            int v2 = m.to(i);

            nlRowScaling(.2);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlRightHandSide(m.point(v1)[d] - m.point(v2)[d]);
            nlEnd(NL_ROW);
        }

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

Но мы на этом не останавливемся!

        for (int v=0; v<m.nverts(); v++) { // amplify the curvature
            nlRowScaling(1.);
            nlBegin(NL_ROW);
            int nneigh = m.incident_halfedges(v).size();
            nlCoefficient(v,  nneigh);
            Vec3f curvature = m.point(v)*nneigh;
            for (int hid : m.incident_halfedges(v)) {
                nlCoefficient(m.to(hid),  -1);
                curvature = curvature - m.point(m.to(hid));
            }
            nlRightHandSide(2.5*curvature[d]);
            nlEnd(NL_ROW);
        }

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

Затем вызываем наш решатель и получаем неплохую карикатуру.

Последний пример на сегодня


До этого момента мы не видели ни одного нового примера, всё уже было показано в предыдущей статье. Давайте попробуем сделать что-то менее тривиальное, нежели простейшее сглаживание, которое легко можно получить, не используя никаких солверов. Код доступен тут, но давайте я приведу листинг: 
#include <vector>
#include <iostream>

#include "OpenNL_psm.h"
#include "geometry.h"
#include "model.h"

const Vec3f axes[] = {Vec3f(1,0,0), Vec3f(-1,0,0), Vec3f(0,1,0), Vec3f(0,-1,0), Vec3f(0,0,1), Vec3f(0,0,-1)};
int snap(Vec3f n) { // returns the coordinate axis closest to the given normal
    double nmin = -2.0;
    int    imin = -1;
    for (int i=0; i<6; i++) {
        double t = n*axes[i];
        if (t>nmin) {
            nmin = t;
            imin = i;
        }
    }
    return imin;
}

int main(int argc, char** argv) {
    if (argc<2) {
        std::cerr << "Usage: " << argv[0] << " obj/model.obj" << std::endl;
        return 1;
    }

    Model m(argv[1]);

    std::vector<int> nearest_axis(m.nfaces());
    for (int i=0; i<m.nfaces(); i++) {
        Vec3f v[3];
        for (int j=0; j<3; j++) v[j] = m.point(m.vert(i, j));
        Vec3f n = cross(v[1]-v[0], v[2]-v[0]).normalize();
        nearest_axis[i] = snap(n);
    }

    for (int d=0; d<3; d++) { // solve for x, y and z separately
        nlNewContext();
        nlSolverParameteri(NL_NB_VARIABLES, m.nverts());
        nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
        nlBegin(NL_SYSTEM);
        nlBegin(NL_MATRIX);

        for (int i=0; i<m.nhalfedges(); i++) {
            int v1 = m.from(i);
            int v2 = m.to(i);

            nlRowScaling(1.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlRightHandSide(m.point(v1)[d] - m.point(v2)[d]);
            nlEnd(NL_ROW);

            int axis = nearest_axis[i/3]/2;
            if (d!=axis) continue;

            nlRowScaling(2.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlEnd(NL_ROW);
        }

        nlEnd(NL_MATRIX);
        nlEnd(NL_SYSTEM);
        nlSolve();

        for (int i=0; i<m.nverts(); i++) {
            m.point(i)[d] = nlGetVariable(i);
        }
    }

    std::cout << m;
    return 0;
}

Мы по-прежнему работаем с той же поверхностью, что и раньше; мы вызываем функцию snap() для каждого треугольника. Эта функция говорит, к какой оси системы координат ближе всего нормаль этого треугольника. То есть, мы хотим знать, этот треугольник скорее вертикален или горизонтален. Ну а затем мы хотим изменить геометрию нашей поверхности, чтобы то, что близко к горизонтали, стало ещё ближе! Левая часть вот этой картинки показывает результат работы функции snap():

Итак, мы раскрасили нашу поверхность в три цвета, соответствующие трём осям системы координат. Ну а затем для каждого ребра нашей системы мы добавляем следующие уравнения:

            nlRowScaling(1.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlRightHandSide(m.point(v1)[d] - m.point(v2)[d]);
            nlEnd(NL_ROW);

То есть, мы предписываем каждому ребру выходной поверхности походить на соответствующее ребро входной поверхности, жёсткость пружинки 1. А затем, если мы разрешаем, например, размерность x, а ребро должно быть вертикально, то просто говорим, что одна вершина равна другой с жёсткостью пружины 2:

            nlRowScaling(2.);
            nlBegin(NL_ROW);
            nlCoefficient(v1,  1);
            nlCoefficient(v2, -1);
            nlEnd(NL_ROW);

Ну а вот и результат работы нашей программы:

Вполне резонный вопрос: истуканы с острова Пасхи, это, конечно, круто, но при чём тут геология? Есть ли примеры реальных задач?

Да, конечно. Давайте возьмём один срез земной коры:

Для геологов очень важны слои разных пород, границы (горизонты) между которыми показаны зелёным, и геологические разломы, которые показаны красным. У нас есть на вход сетка треугольников (тетраэдров в 3д) той области, которая нас интересует, но для моделирования определённых процессов необходимы квады (гексаэдры в 3д). Мы деформируем нашу модель так, чтобы разломы стали вертикальными, а горизонты (сюрприз) горизонтальными:

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

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

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

Let's block ads! (Why?)

на
Ярлыки: ,

Комментариев нет:

Отправить комментарий

Подписаться на: Комментарии к сообщению (Atom)