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

Математика на пальцах: методы наименьших квадратов

Я математик-программист. Самый большой скачок в своей карьере я совершил, когда научился говорить:«Я ничего не понимаю!» Сейчас мне не стыдно сказать светилу науки, что мне читает лекцию, что я не понимаю, о чём оно, светило, мне говорит. И это очень сложно. Да, признаться в своём неведении сложно и стыдно. Кому понравится признаваться в том, что он не знает азов чего-то-там. В силу своей профессии я должен присутствовать на большом количестве презентаций и лекций, где, признаюсь, в подавляющем большинстве случаев мне хочется спать, потому что я ничего не понимаю. А не понимаю я потому, что огромная проблема текущей ситуации в науке кроется в математике. Она предполагает, что все слушатели знакомы с абсолютно всеми областями математики (что абсурдно). Признаться в том, что вы не знаете, что такое производная (о том, что это — чуть позже) — стыдно.

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

Знаете ли вы, что такое производная? Вероятнее всего вы мне скажете про предел разностного отношения. На первом курсе матмеха СПбГУ Виктор Петрович Хавин мне определил производную как коэффициент первого члена ряда Тейлора функции в точке (это была отдельная гимнастика, чтобы определить ряд Тейлора без производных). Я долго смеялся над таким определением, покуда в итоге не понял, о чём оно. Производная не что иное, как просто мера того, насколько функция, которую мы дифференцируем, похожа на функцию y=x, y=x^2, y=x^3.

Я сейчас имею честь читать лекции студентам, которые боятся математики. Если вы боитесь математики — нам с вами по пути. Как только вы пытаетесь прочитать какой-то текст, и вам кажется, что он чрезмерно сложен, то знайте, что он хреново написан. Я утверждаю, что нет ни одной области математики, о которой нельзя говорить «на пальцах», не теряя при этом точности.

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

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

Итак, даны две точки (x0, y0), (x1, y1), например, (1,1) и (3,2), задача найти уравнение прямой, проходящей через эти две точки:

Эта прямая должна иметь уравнение типа следующего:

Здесь альфа и бета нам неизвестны, но известны две точки этой прямой:

Можно записать это уравнение в матричном виде:

Тут следует сделать лирическое отступление: что такое матрица? Матрица это не что иное, как двумерный массив. Это способ хранения данных, более никаких значений ему придавать не стоит. Это зависит от нас, как именно интерпретировать некую матрицу. Периодически я буду её интерпретировать как линейное отображение, периодически как квадратичную форму, а ещё иногда просто как набор векторов. Это всё будет уточнено в контексте.

Давайте заменим конкретные матрицы на их символьное представление:

Тогда (alpha, beta) может быть легко найдено:

Более конкретно для наших предыдущих данных:

Что ведёт к следующему уравнению прямой, проходящей через точки (1,1) и (3,2):

Окей, тут всё понятно. А давайте найдём уравнение прямой, проходящей через три точки: (x0,y0), (x1,y1) и (x2,y2):

Ой-ой-ой, а ведь у нас три уравнения на две неизвестных! Стандартный математик скажет, что решения не существует. А что скажет программист? А он для начала перепишет предыдующую систему уравнений в следующем виде:

И дальше постарается найти решение, которое меньше всего отклонится от заданных равенств. Давайте назовём вектор (x0,x1,x2) вектором i, (1,1,1) вектором j, а (y0,y1,y2) вектором b:

В нашем случае векторы i,j,b трёхмерны, следовательно, (в общем случае) решения этой системы не существует. Любой вектор (alpha\*i + beta\*j) лежит в плоскости, натянутой на векторы (i, j). Если b не принадлежит этой плоскости, то решения не существует (равенства в уравнении не достичь). Что делать? Давайте искать компромисс. Давайте обозначим через e(alpha, beta) насколько именно мы не достигли равенства:

И будем стараться минимизировать эту ошибку:

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

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

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

Совсеми иными словами (осторожно, плохо формализовано, но на пальцах должно быть ясно): мы берём все возможные прямые между всеми парами точек и ищем среднюю прямую между всеми:

Иное объяснение на пальцах: мы прикрепляем пружинку между всеми точками данных (тут у нас три) и прямой, что мы ищем, и прямая равновесного состояния есть именно то, что мы ищем.

Минимум квадратичной формы

Итак, имея данный вектор b и плоскость, натянутую на столбцы-векторы матрицы A (в данном случае (x0,x1,x2) и (1,1,1)), мы ищем вектор e с минимум квадрата длины. Очевидно, что минимум достижим только для вектора e, ортогонального плоскости, натянутой на столбцы-векторы матрицы A:

Иначе говоря, мы ищем такой вектор x=(alpha, beta), что:

Напоминаю, что этот вектор x=(alpha, beta) является минимумом квадратичной функции ||e(alpha, beta)||^2:

Тут нелишним будет вспомнить, что матрицу можно интерпретирвать в том числе как и квадратичную форму, например, единичная матрица ((1,0),(0,1)) может быть интерпретирована как функция x^2 + y^2:

Вся эта гимнастика известна под именем линейной регрессии.

Уравнение Лапласа с граничным условием Дирихле

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

Изначальный коммит доступен здесь. Для минимизации внешних зависимостей я взял код своего софтверного рендерера, уже подробно описанного на хабре. Для решения линейной системы я пользуюсь OpenNL, это отличный солвер, который, правда, очень сложно установить: нужно скопировать два файла (.h+.c) в папку с вашим проектом. Всё сглаживание делается следующим кодом:

X, Y и Z координаты отделимы, я их сглаживаю по отдельности. То есть, я решаю три системы линейных уравнений, каждое имеет количество переменных равным количеству вершин в моей модели. Первые n строк матрицы A имеют только одну единицу на строку, а первые n строк вектора b имеют оригинальные координаты модели. То есть, я привязываю по пружинке между новым положением вершины и старым положением вершины — новые не должны слишком далеко уходить от старых.

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

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

В нашей матрице A я для вершин, что находятся на краю, добавляю не строку из разряда v_i = verts[i][d], а 1000*v_i = 1000*verts[i][d]. Что это меняет? А меняет это нашу квадратичную форму ошибки. Теперь единичное отклонение от вершины на краю будет стоить не одну единицу, как раньше, а 1000*1000 единиц. То есть, мы повесили более сильную пружинку на крайние вершины, решение предпочтёт сильнее растянуть другие. Вот результат:

Давайте вдвое усилим пружинки между вершинами:

Логично, что поверхность стала более гладкой:

А теперь ещё в сто раз сильнее:

Что это? Представьте, что мы обмакнули проволочное кольцо в мыльную воду. В итоге образовавшаяся мыльная плёнка будет стараться иметь наименьшую кривизну, насколько это возможно, касаясь-таки границы — нашего проволочного кольца. Именно это мы и получили, зафиксировав границу и попросив получить гладкую поверхность внутри. Поздравляю вас, мы только что решили уравнение Лапласа с граничными условиями Дирихле. Круто звучит? А на деле всего-навсего одну систему линейных уравнений решить.

Уравнение Пуассона

Давайте ещё крутое имя вспомним.

Предположим, что у меня есть такая картинка:

Всем хороша, только стул мне не нравится.

Разрежу картинку пополам:

И выделю руками стул:

Затем всё, что белое в маске, притяну к левой части картинки, а заодно по всей картинке скажу, что разница между двумя соседними пикселями должна равняться разнице между двумя соседними пикселями правой картинки:

Код и картинки доступны здесь.

Пример из жизни

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

У меня есть некоторое количество фотографий образцов ткани типа вот такой:

Моя задача сделать бесшовные текстуры из фотографий вот такого качества. Для начала я (автоматически) ищу повторяющийся паттерн:

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

Вот фрагмент, где чётко видно шов:

Поэтому я вырезать буду не по ровной линии, вот линия разреза:

А вот повторённый четыре раза паттерн:

И его фрагмент, чтобы было виднее:

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

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

Источник

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

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

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

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

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

Матрицы и числа

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

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

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

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

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

Ну а в мире программистов она бы выглядела следующим образом:

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

Раз уж я начал про программистов, я обязан записать её код 🙂

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

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

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

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

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

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

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

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

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

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

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

Число положительно тогда и только тогда, когда для всех ненулевых вещественных выполняется условие 0$» data-tex=»inline»/>.

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

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

Квадратная матрица называется положительно определённой, если для любых ненулевых выполняется условие 0$» data-tex=»inline»/>, то есть, соответствующая квадратичная форма строго положительна везде, кроме начала координат.

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

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

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

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

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

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

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

Единственное, что нам осталось, так это убедиться, что эта эквивалентность в одномерном мире распространяется и на случай переменных. Чтобы это проверить, для начала докажем три теоремы.

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

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

Теорема 1

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

Вторая теорема позволяет дифференцировать линейные функции. В случае вещественной функции одной переменной мы знаем, что , но что происходит в случае вещественной функции переменных?

Теорема 2

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

Применив первую теорему , мы завершаем доказательство.

Теперь переходим к квадратичным формам. Мы знаем, что в случае вещественной функции одной переменной , а что будет в случае квадратичной формы?

Теорема 3

Кстати, обратите внимание, что если матрица симметрична, то теорема гласит, что .

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

А затем посмотрим, чему равна частная производная этой двойной суммы по переменной :
\begin
\frac<\partial (x^\top A x)> <\partial x_i>
&= \frac<\partial> <\partial x_i>\left(\sum\limits_\sum\limits_ a_ x_ x_\right) = \\
&= \frac<\partial> <\partial x_i>\left(
\underbrace<\sum\limits_\sum\limits_ a_x_ x_>_+\underbrace <\sum\limits_a_x_i x_>_+
\underbrace <\sum\limits_a_ x_ x_i>_+
\underbracex_i^2>_\right) = \\
& = \sum\limits_ a_x_ + \sum\limits_ a_ x_ + 2 a_ x_i = \\
& = \sum\limits_ a_x_ + \sum\limits_ a_ x_ = \\
& = \sum\limits_ (a_ + a_) x_j \\
\end
Я разбил двойную сумму на четыре случая, которые выделены фигурными скобками. Каждый из этих четырёх случаев дифференцируется тривиально. Осталось сделать последнее действие, собрать частные производные в вектор градиента:

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

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

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

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

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

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

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

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

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

У нас два уравнения с двумя неизвестными ( и ), остальное известно.

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

Получим уравнение типа , которое тривиально решается .

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

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

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

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

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

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

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

Но постойте, где же наименьшие квадраты? Сейчас будут. Функция извлечения корня монотонна, поэтому = !

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

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

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

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

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

Теорема 4

симметрична.

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

Теорема 5

положительно полуопределена:

Это следует из того факта, что 0$» data-tex=»inline»/>.

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

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

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

Переходим к практике

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

Давайте вспомним один пример из предыдущей статьи:

У нас есть обычный массив 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. Ой. А ведь эта система не имеет решения :((

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

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

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

Этот код производит список x из 32 элементов, в котором хранится решение. Мы строим матрицу , строим вектор , получаем вектор с решением , а затем в этот вектор вставляем наши фиксированные значения x0=1, x18=2 и x31=1.

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

С программистской точки зрения он гораздо приятнее, вектор 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(колво вершин). Внимание, мы солверу даём матрицу , а уж он будет строить сам, что очень удобно! Также обратите внимание на то, что я решаю не одну систему уравнений, а три последовательно для x,y и z (я соврал, что проблема трёхмерная, она по-прежнему одномерная!)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заключение

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

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

Источник