Заметки о теоретической физике → 2013 → 11 → 17
Роман Парпалак

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

17 ноября 2013 года, 17:35

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

В простейшем случае метод наименьших квадратов применяется для проведения прямой линии через набор экспериментальных точек и состоит в минимизации суммы квадратов отклонений $$\inline\sum(y_i-ax_i-b)^2$$, которые списываются на погрешность измерений. В результате минимизации для коэффициентов a и b получается простая система линейных уравнений. Здесь важно предположение о том, что ошибки по оси x пренебрежимо малы по сравнению с ошибками по оси y. Если это не так, то минимизировать нужно более сложное выражение.

Иногда возникает задача другого рода — провести геометрическую прямую через набор геометрических точек «наилучшим образом». Для этой задачи метод наименьших квадратов нужно адаптировать, так как поспешное применение формул для коэффициентов a и b будет давать разные прямые в разных системах координат. Теперь отклонения по осям должны быть одинаковы. Правильный подход заключается в минимизации суммы квадратов расстояний $$\inline\sum(y_i-ax_i-b)^2/(1+a^2)$$ от точек (xi, yi) до проводимой прямой. Он дает нелинейную систему уравнений, которую можно решать численно. Однако этот подход тяжело обобщается на интересующий меня многомерный случай. Поэтому мы с самого начала будем рассматривать задачу во многомерном пространстве.

Задача

Пусть задан набор точек $$\vec{x}^k$$. Мы хотим провести гиперплоскость $$(\vec{n}\cdot\vec{x}) = d$$ такую, что сумма квадратов расстояний от точек $$\vec{x}^k$$ до нее будет минимальна. Расстояние до гиперплоскости находится с помощью проекции на единичный вектор нормали $$\vec{n}$$, и выражение для минимизации принимает вид

$$\sum_k\left((\vec{n}\cdot\vec{x}^k)-d\right)^2\to\text{min}.$$

При этом нужно учитывать уравнение связи $$(\vec{n}\cdot\vec{n}) = 1$$, которое уменьшает на 1 количество степеней свободы в неизвестных величинах ni, d. Учет связи выполняется с помощью метода множителей Лагранжа. Однако мы пойдем другим путем, который сократит выкладки и напрямую приведет к выражениям, подходящим для численного счета. Мы разрешим вектору $$\vec{n}$$ иметь произвольную длину, и введем явную нормировку:

$$\sum_k\left({(\vec{n}\cdot\vec{x}^k)\over|\vec{n}|}-d\right)^2\to\text{min}.$$

Параллельный перенос

Продифференцируем по d:

$$\sum_k\left({(\vec{n}\cdot\vec{x}^k)\over|\vec{n}|}-d\right)={(\vec{n}\cdot\sum_k\vec{x}^k)\over|\vec{n}|}-\sum_kd=0.$$

Как видим, «центр масс» набора точек $$\inline\sum\vec{x}^k/\sum 1$$ находится на искомой плоскости. Выполним параллельный перенос системы координат таким образом, чтобы ее начало совпало с центром набора точек $$\inline\sum\vec{x}^k=0$$. В этой системе координат d=0.

Условие на вектор нормали

Перейдем к индексным обозначениям и продифференцируем по na:

$${\partial\over\partial n_a}\left({n_ix_i^k\,n_jx_j^k\over n_ln_l}\right)={2x_a^k\,n_jx_j^k\over n_ln_l}-{2n_a\,n_ix_i^k\,n_jx_j^k\over n_ln_l\,n_pn_p}=0,$$

$$n_jx_j^k\left[x_a^k(n_pn_p)-n_a(n_ix_i^k)\right]=0,$$

$$\sum_k\vec{n}\cdot\vec{x}^k\left[\vec{x}^k(\vec{n}\cdot\vec{n})-\vec{n}(\vec{n}\cdot\vec{x}^k)\right]=0.$$

Вычислительный аспект

Нелинейное уравнение относительно вектора $$\vec{n}$$ можно решать методом итераций:

$$\sum_k\vec{n}_{i+1}\cdot\vec{x}^k\left[\vec{x}^k(\vec{n}_i\cdot\vec{n}_i)-\vec{n}_i(\vec{n}_i\cdot\vec{x}^k)\right]=0.$$

С помощью матрицы $$A_{aj}(\vec{n})=\left[x_a^k(n_pn_p)-n_a(n_ix_i^k)\right]x_j^k$$ оно представляется в виде

$$A(\vec{n}_i)\,\vec{n}_{i+1}=0$$

и сводится к поиску ядра линейного оператора. Нетривиальность ядра связана с «лишней» степенью свободы, появившейся из-за отбрасывания условия нормировки вектора нормали. Читатели могут самостоятельно проверить с помощью формулы Бине — Коши, что определитель матрицы $$A(\vec{n}_i)$$ равен нулю.

По теореме Фредгольма ядро оператора ортогонально образу сопряженного оператора, то есть линейной оболочке, натянутой на строки $$\vec{a}_a$$ матрицы $$A_{aj}(\vec{n}_i)$$. Алгоритм поиска ортогонального дополнения состоит в выборе произвольного вектора $$\vec{r}$$ и ортогонализации набора векторов $$\vec{r}, \vec{a}_a$$:

$$\vec{r}^{\,\prime}=\vec{r}-\vec{a}_1{(\vec{r}\cdot\vec{a}_1)\over(\vec{a}_1\cdot\vec{a}_1)},$$

$$\vec{a}_2^{\,\prime}=\vec{a}_2-\vec{a}_1{(\vec{a}_2\cdot\vec{a}_1)\over(\vec{a}_1\cdot\vec{a}_1)},\quad\vec{r}^{\,\prime\prime}=\vec{r}^{\,\prime}-\vec{a}_2^{\,\prime}{(\vec{r}^{\,\prime}\cdot\vec{a}_2^{\,\prime})\over(\vec{a}_2^{\,\prime}\cdot\vec{a}_2^{\,\prime})}\ldots$$

Так как строки матрицы $$A(\vec{n}_i)$$ линейно зависимы, один из векторов $$\vec{a}_a$$ при ортогонализации из набора исключается. Для большей определенности алгоритма в качестве начального приближения перебираем базисные векторы, пока в результате ортогонализации не получится ненулевой вектор следующего приближения $$\vec{n}_{i+1}$$. В двумерном и трехмерном случае процесс ортогонализации значительно упрощается. Например, в трехмерном случае нетривиальный элемент ядра найдется среди тройки векторов $$\vec{a}_1\times\vec{a}_2, \vec{a}_1\times\vec{a}_3, \vec{a}_2\times\vec{a}_3$$.

Как показывают практические вычисления, последовательные приближения $$\vec{n}_i$$ быстро сходятся к искомому вектору нормали.

Ключевые слова: геометрия

Комментарии

#1. 27 декабря 2013 года, 22:29. ud1 пишет:
Думаю, что можно решать намного проще из общефизических соображений. Т.к. центр масс нашли, и он замечательно оказался лежащем на искомой плоскости, то дальше логично было бы найти тензор инерции относительно этого центра масс. И найти систему координат, в которой этот тензор бы имел бы диагональный вид. Тогда, если я не ошибаюсь, ось координат для которой момент инерции принимает наибольшее значение и будет нормалью искомой плоскости.
#2. 29 декабря 2013 года, 14:02. пишет:
Интересная догадка. Момент инерции системы есть $$J_{ij}=\delta_{ij}x_p^kx_p^k-x_i^kx_j^k$$. Наибольшее собственное число соответствует максимуму квадратичной формы $$n_iJ_{ij}n_j=n_in_ix_p^kx_p^k-n_ix_i^kn_jx_j^k$$ на единичном векторе $$\vec{n}$$. Первая часть не зависит от $$\vec{n}$$, максимум второй части соответствует искомому в задаче минимуму.

Тут всплывает другой момент, который я не заметил. Минимизация выражения $$n_ix_i^kn_jx_j^k$$ есть поиск минимального значения квадратичной формы с матрицей $$x_i^kx_j^k$$. Оно достигается на собственном векторе этой матрицы. На этом математическую задачу можно считать решенной. Что касается вычислительной стороны, то мои итерации всё равно проще (и быстрее) диагонализации симметричных матриц, обычно выполняемой методом Якоби: http://en.wikipedia.org/wiki/Jacobi_eigenv … _algorithm
#3. 7 августа 2014 года, 23:03. Михаил пишет:
Хороший материал, прочел с удовольствием. Мелкие добавления:
1) В 3-м абзаце — надо убрать лишний корень (из 1+a^2), ведь всё в квадрате!
2) В случае плоскости итерации не нужны — есть простое конечное решение (см., напр., Chernov N. Circular and linear regression: fitting circles and lines by least squares // Monographs on statistics and applied probability 117; CRC Press, Taylor & Francis Group, 2011, 254 pp.).
3) Интересно было бы также провести через эти точки «наилучшую» прямую (в n-мерном пространстве). Я нигде не встречал даже постановки такой задачи (даже для n=3), не говоря о решении. А оно как-то не укладывается в стандартные регрессионные модели.
#4. 7 августа 2014 года, 23:34. пишет:
И ещё: с общих позиций ортогональная регрессия хорошо изложена в:
Суслов В.И., Ибрагимов Н.М., Талышева Л.П., Цыплаков А.А.
Эконометрия: Учебное пособие. — Новосибирск: Издательство СО РАН, 2005. — 744 с.
http://institutiones.com/general/1502-ekon … uslov.html
#5. 8 августа 2014 года, 00:33. пишет:
Михаил, спасибо за исправление. Корень действительно был лишний.

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

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

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

Ваше имя:

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

Формулы на латехе: $$f(x) = x^2-\sqrt{x}$$ превратится в $$f(x) = x^2-\sqrt{x}$$.
Для выделения используйте следующий код: [i]курсив[/i], [b]жирный[/b].
Цитату оформляйте так: [q = имя автора]цитата[/q] или [q]еще цитата[/q].
Ссылку начните с http://. Других команд или HTML-тегов здесь нет.

Сколько будет 47+2?