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

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

В простейшем случае метод наименьших квадратов применяется для проведения прямой линии через набор экспериментальных точек и состоит в минимизации суммы квадратов отклонений $$\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)$$ от точек (x_i, y_i) до проводимой прямой. Он дает нелинейную систему уравнений, которую можно решать численно. Однако этот подход тяжело обобщается на интересующий меня многомерный случай. Поэтому мы с самого начала будем рассматривать задачу во многомерном пространстве.

Задача

Пусть задан набор точек $$\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 количество степеней свободы в неизвестных величинах n_i, 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.

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

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

$${\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$$ быстро сходятся к искомому вектору нормали.