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

в 13:45, , рубрики: c++, Алгоритмы, геометрия для пятого класса, математика, математика для программистов, математика на пальцах, наименьшие квадраты, Программирование, Работа с анимацией и 3D-графикой, уравнение лапласа, уравнение Пуассона, метки:

Введение

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

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

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

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

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

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

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

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

иллюстрация

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Почему квадрат?

Мы ищем не просто минимум нормы, а минимум квадрата нормы. Почему? Сама точка минимума совпадает, а квадрат даёт гладкую функцию (квадратичную функцию от агрументов (alpha,beta)), в то время как просто длина даёт функцию в виде конуса, недифференцируемую в точке минимума. Брр. Квадрат удобнее.

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

Иллюстрация

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

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

Иллюстрация

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

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

Иллюстрация

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

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

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

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

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

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

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

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

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

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

квадратичная форма

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

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

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

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

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

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

    for (int d=0; d<3; d++) {
        nlNewContext();
        nlSolverParameteri(NL_NB_VARIABLES, verts.size());
        nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
        nlBegin(NL_SYSTEM);
        nlBegin(NL_MATRIX);

        for (int i=0; i<(int)verts.size(); i++) {
            nlBegin(NL_ROW);
            nlCoefficient(i, 1);
            nlRightHandSide(verts[i][d]);
            nlEnd(NL_ROW);
        }

        for (unsigned int i=0; i<faces.size(); i++) {
            std::vector<int> &face = faces[i];
            for (int j=0; j<3; j++) {
                nlBegin(NL_ROW);
                nlCoefficient(face[ j     ],  1);
                nlCoefficient(face[(j+1)%3], -1);
                nlEnd(NL_ROW);
            }
        }

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

        for (int i=0; i<(int)verts.size(); i++) {
            verts[i][d] = nlGetVariable(i);
        }
    }

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

Все последующие строки матрицы A (faces.size()*3 = количеству рёбер всех треугольников в сетке) имеют одно вхождение 1 и одно вхождение -1, причём вектор b имеет нулевые компоненты напротив. Это значит, я вешаю пружинку на каждое ребро нашей треугольной сетки: все рёбра стараются получить одну и ту же вершину в качестве отправной и финальной точки.

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

Вот результат:

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

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

        for (int i=0; i<(int)verts.size(); i++) {
            float scale = border[i] ? 1000 : 1;
            nlBegin(NL_ROW);
            nlCoefficient(i, scale);
            nlRightHandSide(scale*verts[i][d]);
            nlEnd(NL_ROW);
        }

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

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

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

                nlCoefficient(face[ j     ],  2);
                nlCoefficient(face[(j+1)%3], -2);

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

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

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

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

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

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

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

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

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

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

Разрежу картинку пополам:
Математика на пальцах: методы наименьших квадратов - 30
Математика на пальцах: методы наименьших квадратов - 31

И выделю руками стул:
Математика на пальцах: методы наименьших квадратов - 32

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

    for (int i=0; i<w*h; i++) {
        if (m.get(i%w, i/w)[0]<128) continue;
        nlBegin(NL_ROW);
        nlCoefficient(i, 1);
        nlRightHandSide(a.get(i%w, i/w)[0]);
        nlScaleRow(100);
        nlEnd(NL_ROW);
    }

    for (int i=0; i<w-1; i++) {
        for (int j=0; j<h-1; j++) {
            for (int d=0; d<2; d++) {
                nlBegin(NL_ROW);
                int v1 = b.get(i,   j    )[0];
                int v2 = b.get(i+d, j+1-d)[0];
                nlCoefficient( i   + j     *w,  1);
                nlCoefficient((i+d)+(j+1-d)*w, -1);
                nlRightHandSide(v1-v2);
                nlScaleRow(10);
                nlEnd(NL_ROW);
            }
        }
    }

Вот результат:

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

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

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

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

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

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

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

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

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

Скрытый текст

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

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

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

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

Скрытый текст

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

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

Скрытый текст

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

И его фрагмент, чтобы было виднее:
Математика на пальцах: методы наименьших квадратов - 40

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

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

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

Автор: haqreu

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js