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

В прикладной статистике метод наименьших полных квадратов (МНПК, TLS — англ.  Total Least Squares) — это вид регрессии с ошибками в переменных^[en], техника моделирования данных с помощью метода наименьших квадратов, в которой принимаются во внимание ошибки как в зависимых, так и в независимых переменных. Метод является обобщением регрессии Деминга и ортогональной регрессии и может быть применён как к линейным, так и нелинейным моделям.

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

Аппроксимация данных методом наименьших полных квадратов в общем виде эквивалентна лучшей по норме Фробениуса малоранговой аппрокимации^[en] матрицы данных.

Линейная модельПравить

ОсновыПравить

В методе наименьших квадратов моделирования данных минимизируется функция потерь S,

$S={\mathbf{r}}^{\mathbf{T}}\mathbf{W}\mathbf{r},$$S=\mathbf {r^{T}Wr} ,$ 

где r — вектор отклонений^[en], а W — весовая матрица. В линейном методе наименьших квадратов модель содержит уравнения, которые линейны от параметров в векторе $\mathit{\beta <!--\; \beta \; -->}$${\boldsymbol {\beta }}$ , так что отклонения вычисляются по формуле

$\mathbf{r}\mathbf{=}\mathbf{y}\mathbf{-<!--\; -\; -->}\mathbf{X}\mathit{\beta <!--\; \beta \; -->}.$$\mathbf {r=y-X{\boldsymbol {\beta }}} .$ 

Имеется m наблюдений в векторе y и n параметров в β при m>n. X является m×n матрицей, элементы которой либо константы, либо функции от независимых переменных x. Весовая матрица W, в идеале, является обратной к дисперсионно-ковариационной матрице ${\mathbf{M}}_y$$\mathbf {M} _{y}$  наблюдений y. Предполагается, что независимые переменные ошибок не имеют. Параметры оценки находятся путём приравнивания градиента нулю, что приводит к уравнению

${\mathbf{X}}^{\mathbf{T}}\mathbf{W}\mathbf{X}\mathit{\beta <!--\; \beta \; -->}\mathbf{=}{\mathbf{X}}^{\mathbf{T}}\mathbf{W}\mathbf{y}$$\mathbf {X^{T}WX{\boldsymbol {\beta }}=X^{T}Wy}$ 

Возможность ошибок наблюдений для всех переменныхПравить

Предположим теперь, что как x, так и y наблюдаются с ошибками с дисперсионно-ковариационными матрицами ${\mathbf{M}}_x$$\mathbf {M} _{x}$  и ${\mathbf{M}}_y$$\mathbf {M} _{y}$  соответственно. В этом случае функция потерь записывается как

$S={\mathbf{r}}_{\mathbf{x}}^{\mathbf{T}}{\mathbf{M}}_{\mathbf{x}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}{\mathbf{r}}_{\mathbf{x}}\mathbf{+}{\mathbf{r}}_{\mathbf{y}}^{\mathbf{T}}{\mathbf{M}}_{\mathbf{y}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}{\mathbf{r}}_{\mathbf{y}}$$S=\mathbf {r_{x}^{T}M_{x}^{-1}r_{x}+r_{y}^{T}M_{y}^{-1}r_{y}}$ ,

где ${\mathbf{r}}_x$$\mathbf {r} _{x}$  и ${\mathbf{r}}_y$$\mathbf {r} _{y}$  являются отклонениями для x и y соответственно. Ясно, что эти отклонения не могут быть независимыми и между ними должна быть какая-то связь. Если записать функцию как $\mathbf{f}\mathbf{(}{\mathbf{r}}_{\mathbf{x}}\mathbf{,}{\mathbf{r}}_{\mathbf{y}}\mathbf{,}\mathit{\beta <!--\; \beta \; -->}\mathbf{)}$$\mathbf {f(r_{x},r_{y},{\boldsymbol {\beta }})}$ , ограничения выражаются m условиями.

$\mathbf{F}\mathbf{=}\mathbf{\Delta <!--\; \Delta \; -->}\mathbf{y}\mathbf{-<!--\; -\; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{f}}{\mathrm{\partial <!--\; \partial \; -->}{\mathbf{r}}_{\mathbf{x}}}{\mathbf{r}}_{\mathbf{x}}\mathbf{-<!--\; -\; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{f}}{\mathrm{\partial <!--\; \partial \; -->}{\mathbf{r}}_{\mathbf{y}}}{\mathbf{r}}_{\mathbf{y}}\mathbf{-<!--\; -\; -->}\mathbf{X}\mathbf{\Delta <!--\; \Delta \; -->}\mathit{\beta <!--\; \beta \; -->}\mathbf{=}\mathbf{0}$$\mathbf {F=\Delta y-{\frac {\partial f}{\partial r_{x}}}r_{x}-{\frac {\partial f}{\partial r_{y}}}r_{y}-X\Delta {\boldsymbol {\beta }}=0}$ 

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

${\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{X}\mathbf{\Delta <!--\; \Delta \; -->}\mathit{\beta <!--\; \beta \; -->}\mathbf{=}{\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{\Delta <!--\; \Delta \; -->}\mathbf{y},$$\mathbf {X^{T}M^{-1}X\Delta {\boldsymbol {\beta }}=X^{T}M^{-1}\Delta y} ,$ 

или, альтернативно, ${\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{X}\mathit{\beta <!--\; \beta \; -->}\mathbf{=}{\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{y}$$\mathbf {X^{T}M^{-1}X{\boldsymbol {\beta }}=X^{T}M^{-1}y}$ 

Здесь M — дисперсионно-ковариационная матрица, относящаяся как к независимым, так и зависимым переменным.

$\mathbf{M}\mathbf{=}{\mathbf{K}}_{\mathbf{x}}{\mathbf{M}}_{\mathbf{x}}{\mathbf{K}}_{\mathbf{x}}^{\mathbf{T}}\mathbf{+}{\mathbf{K}}_{\mathbf{y}}{\mathbf{M}}_{\mathbf{y}}{\mathbf{K}}_{\mathbf{y}}^{\mathbf{T}}\mathbf{;}{\mathbf{K}}_{\mathbf{x}}\mathbf{=}\mathbf{-<!--\; -\; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{f}}{\mathrm{\partial <!--\; \partial \; -->}{\mathbf{r}}_{\mathbf{x}}}\mathbf{,}{\mathbf{K}}_{\mathbf{y}}\mathbf{=}\mathbf{-<!--\; -\; -->}\frac{\mathrm{\partial <!--\; \partial \; -->}\mathbf{f}}{\mathrm{\partial <!--\; \partial \; -->}{\mathbf{r}}_{\mathbf{y}}}$$\mathbf {M=K_{x}M_{x}K_{x}^{T}+K_{y}M_{y}K_{y}^{T};\ K_{x}=-{\frac {\partial f}{\partial r_{x}}},\ K_{y}=-{\frac {\partial f}{\partial r_{y}}}}$ 

ПримерПравить

В случае, когда ошибки данных не коррелируют, все матрицы M и W диагональны. Тогда используем построение прямой по точкам.

$f(x_i,\mathrm{\beta <!--\; \beta \; -->})=\mathrm{\alpha <!--\; \alpha \; -->}+\mathrm{\beta <!--\; \beta \; -->}{x}_i$$f(x_{i},\beta )=\alpha +\beta x_{i}\!$ 

И в этом случае

$M_{ii}={\mathrm{\sigma <!--\; \sigma \; -->}}_{y,i}^2+{\mathrm{\beta <!--\; \beta \; -->}}^2{\mathrm{\sigma <!--\; \sigma \; -->}}_{x,i}^2$$M_{ii}=\sigma _{y,i}^{2}+\beta ^{2}\sigma _{x,i}^{2}$ 

что показывает, как дисперсия в i-ой точке определяется дисперсией независимых и зависимых переменных, а также моделью, используемой для согласования данных. Выражение можно обобщить, если заметить, что параметр $\mathrm{\beta <!--\; \beta \; -->}$$\beta$  является наклоном прямой.

$M_{ii}={\mathrm{\sigma <!--\; \sigma \; -->}}_{y,i}^2+{\left(\frac{dy}{dx}\right)}_i^2{\mathrm{\sigma <!--\; \sigma \; -->}}_{x,i}^2$$M_{ii}=\sigma _{y,i}^{2}+\left({\frac {dy}{dx}}\right)_{i}^{2}\sigma _{x,i}^{2}$ 

Выражение такого вида используется для аппроксимации данные титрования pH^[en], когда малые ошибки в x дают большие ошибки y в случае большого наклона.

С алгебраической точки зренияПравить

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

Вычисление МНПК с помощью сингулярного разложения описан в стандартных текстах. Мы можем решить уравнение

$XB\approx <!--\; \approx \; -->Y$$XB\approx Y$ 

относительно B, где X — матрица m-на-n, а Y — матрица m-на-k

То есть мы пытаемся найти матрицу B, минимизирующую матрицы ошибок R и F для X и Y соответственно. То есть

${\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{R,F}\Vert <!--\; \Vert \; -->[RF]{\Vert <!--\; \Vert \; -->}_F,(X+R)B=Y+F$$\mathrm {argmin} _{R,F}\|[R\;F]\|_{F},\qquad (X+R)B=Y+F$ ,

где $[RF]$$[R\;F]$  — расширенная матрица с R и F рядом и $\Vert <!--\; \Vert \; -->\cdot <!--\; \cdot \; -->{\Vert <!--\; \Vert \; -->}_F$$\|\cdot \|_{F}$  является нормой матрицы, квадратным корнем из суммы квадратов всех элементов матриц, что эквивалентно квадратному корню из суммы квадратов длин строк или столбцов матрицы.

Это можно переписать как

$[(X+R)(Y+F)]\left[\begin{array}{c}B\\ -<!--\; -\; -->E_k\end{array}\right]=0.$$[(X+R)\;(Y+F)]{\begin{bmatrix}B\\-E_{k}\end{bmatrix}}=0.$ 

Где $E_k$$E_k$  является $k\times <!--\; \times \; -->k$$k\times k$  единичной матрицей. Целью является нахождение матрицы $[RF]$$[R\;F]$ , которая уменьшает ранг $[XY]$$[X\;Y]$  на k. Определим $[U][\mathrm{\Sigma <!--\; \Sigma \; -->}][V]\ast <!--\; \ast \; -->$$[U][\Sigma ][V]*$  как сингулярное разложение расширенной матрицы $[XY]$$[X\;Y]$ .

$[X\;Y]=[U_{X}\;U_{Y}]{\begin{bmatrix}\Sigma _{X}&0\\0&\Sigma _{Y}\end{bmatrix}}{\begin{bmatrix}V_{XX}&V_{XY}\\V_{YX}&V_{YY}\end{bmatrix}}^{*}=[U_{X}\;U_{Y}]{\begin{bmatrix}\Sigma _{X}&0\\0&\Sigma _{Y}\end{bmatrix}}{\begin{bmatrix}V_{XX}^{*}&V_{YX}^{*}\\V_{XY}^{*}&V_{YY}^{*}\end{bmatrix}}$ ,

где V разбита на блоки, соответствующие формам матриц X и Y.

Если использовать теорему Экарта-Янга, аппроксимация, минимизирующая норму ошибки, это такая аппроксимация, что матрицы $U$$U$  и $V$$V$  не меняются, в то время как $k$$k$  наименьших сингулярных значений заменяются нулями. То есть мы хотим

$[(X+R)\;(Y+F)]=[U_{X}\;U_{Y}]{\begin{bmatrix}\Sigma _{X}&0\\0&0_{k\times k}\end{bmatrix}}{\begin{bmatrix}V_{XX}&V_{XY}\\V_{YX}&V_{YY}\end{bmatrix}}^{*}$ 

так что, ввиду линейности,

$[R\;F]=-[U_{X}\;U_{Y}]{\begin{bmatrix}0_{n\times n}&0\\0&\Sigma _{Y}\end{bmatrix}}{\begin{bmatrix}V_{XX}&V_{XY}\\V_{YX}&V_{YY}\end{bmatrix}}^{*}.$ 

Мы можем удалить блоки из матриц U и Σ, упростив выражение до

$[RF]=-<!--\; -\; -->U_Y{\mathrm{\Sigma <!--\; \Sigma \; -->}}_Y{\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]}^{\ast <!--\; \ast \; -->}=-<!--\; -\; -->[XY]\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]{\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]}^{\ast <!--\; \ast \; -->}.$$[R\;F]=-U_{Y}\Sigma _{Y}{\begin{bmatrix}V_{XY}\\V_{YY}\end{bmatrix}}^{*}=-[X\;Y]{\begin{bmatrix}V_{XY}\\V_{YY}\end{bmatrix}}{\begin{bmatrix}V_{XY}\\V_{YY}\end{bmatrix}}^{*}.$ 

Это даёт R и F, таки что

$[(X+R)(Y+F)]\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]=0.$$[(X+R)\;(Y+F)]{\begin{bmatrix}V_{XY}\\V_{YY}\end{bmatrix}}=0.$ 

Теперь, если $V_{YY}$$V_{YY}$  не вырождена, что не всегда верно (заметим, что поведение МНПК в случае вырожденности $V_{YY}$$V_{YY}$  не вполне понятно), мы можем умножить справа обе части на $-<!--\; -\; -->V_{YY}^{-<!--\; -\; -->1}$$-V_{YY}^{-1}$ , чтобы привести нижний блок правой матрицы к отрицательной единичной матрице, что даёт

$[(X+R)(Y+F)]\left[\begin{array}{c}-<!--\; -\; -->V_{XY}{V}_{YY}^{-<!--\; -\; -->1}\\ -<!--\; -\; -->V_{YY}{V}_{YY}^{-<!--\; -\; -->1}\end{array}\right]=[(X+R)(Y+F)]\left[\begin{array}{c}B\\ -<!--\; -\; -->E_k\end{array}\right]=0,$$[(X+R)\;(Y+F)]{\begin{bmatrix}-V_{XY}V_{YY}^{-1}\\-V_{YY}V_{YY}^{-1}\end{bmatrix}}=[(X+R)\;(Y+F)]{\begin{bmatrix}B\\-E_{k}\end{bmatrix}}=0,$ 

а тогда

$B=-<!--\; -\; -->V_{XY}{V}_{YY}^{-<!--\; -\; -->1}.$$B=-V_{XY}V_{YY}^{-1}.$ 

Имплементация в системе GNU Octave:

function B = tls(X,Y)

[m n]   = size(X);            % n является шириной матрицы X (X[m x n])
Z       = [X Y];              % Z является расширением X на Y.
[U S V] = svd(Z,0);           % находим [[Сингулярное разложение|SVD]] матрицы Z.
VXY     = V(1:n,1+n:end);     % Берём блок матрицы V, состоящий из первых n строк и n+1 последних столбцов
VYY     = V(1+n:end,1+n:end); % Берём нижний правы блок матрицы V.
B       = -VXY/VYY;

end

Метод решения задачи, описанный выше и требующий, чтобы матрица $V_{YY}$$V_{YY}$  не была вырожденной, может быть слегка расширен так называемым классическим МНПК алгоритмом.

ВычислениеПравить

Стандартная имплементация классического алгоритма МНПК доступна на Netlib, см. также статьи. Все современные имплементации, базирующиеся, например, на использовании обычного метода наименьших квадратов, аппроксимируют матрицу $B$$B$  (которая в литературе обозначается как $X$$X$ ), как это делают Ван Хуффель и Вандевалле. Стоит заметить, однако, что полученная матрица $B$$B$  во многих случаях не является решением МНПК.

Нелинейная модельПравить

Для нелинейных систем^[en] похожие рассуждения показывают, что нормальное уравнение для итерационного цикла может быть переписано как

${\mathbf{J}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{J}\mathbf{\Delta <!--\; \Delta \; -->}\mathit{\beta <!--\; \beta \; -->}\mathbf{=}{\mathbf{J}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-<!--\; -\; -->}\mathbf{1}}\mathbf{\Delta <!--\; \Delta \; -->}\mathbf{y}.$$\mathbf {J^{T}M^{-1}J\Delta {\boldsymbol {\beta }}=J^{T}M^{-1}\Delta y} .$ 

Геометрическая интерпретацияПравить

Основная статья: Приближение с помощью кривых § Алгебраическое приближение и геометрическое приближение кривыми

Подробное рассмотрение темы: Регрессия Деминга

Если независимые переменные ошибок не имеют, отклонения представляют «вертикальное» расстояние между точкой данных и аппроксимирующей кривой (или поверхностью). В методе наименьших полных квадратов отклонения представляют расстояние между точкой данных и аппроксимирующей кривой, измеряемое в некотором направлении. Фактически, если обе переменные измеряются в одинаковых единицах и ошибки обоих переменных те же самые, то отклонение представляет кратчайшее расстояние от точки данных до аппроксимирующее кривой, то есть вектор отклонения перпендикулярен касательной к кривой. По этой причине этот тип регрессии называют иногда двумерной евклидовой регрессией или ортогональной регрессией.

Масштабно-инвариантные методыПравить

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

Кратко, метод наименьших полных квадратов не имеет свойства инвариантности по единицам измерения, т.е. он не является масштабно инвариантным. Для полноценности модели мы требуем, чтобы это свойство выполнялось. Дальнейшее продвижение вперёд, это понимание, что отклонения (расстояния), измеряемые в других единицах, могут быть скомбинированы, если используется умножение, а не сложение. Рассмотрим аппроксимацию прямой, для каждой точки данных произведение горизонтального и вертикального отклонений равно удвоенной площади треугольника, образованного отрезками отклонений и аппроксимирующей прямой. Мы выбираем прямую, минимизирующую сумму этих площадей. Нобелевский лауреат Пол Самуэльсон доказал в 1942, что в двумерном случае эта прямая выражается исключительно в терминах отношений квадратических отклонений и корреляции коэффициентов, которые (1) удовлетворяют уравнению, если наблюдения находятся на прямой линии; (2) обнаруживают масштабную инвариантность, (3) обнаруживают инвариантность при обмене переменных. Эта прямая переоткрывалась в различных дисциплинах и известна как стандартизованная главная ось, приведённая главная ось, функциональная зависимость средних геометрических, регрессия наименьших квадратов, диагональная регрессия и прямая наименьших площадей. Тофаллис расширил этот подход для работы с несколькими переменными.

См. такжеПравить

• Регрессия Деминга
• Модель с ошибками в переменных^[en]
• Линейная регрессия
• Метод наименьших квадратов

ЛитератураПравить

Дополнительная литератураПравить

• Paige C. C., Strakoš Z.,. Core problems in linear algebraic systems // SIAM J. Matrix Anal. Appl. — 2006. — Т. 27. — С. 861—875.
• Jo S., Kim S. W. Consistent normalized least mean square filtering with noisy data matrix. — 2005. — Т. 53. — С. 2112—2123. — (IEEE Trans. Signal Processing).
• DeGroat R. D., Dowling E. M. The data least squares problem and channel equalization. — 1993. — Т. 41. — С. 407—411. — (IEEE Trans. Signal Processing).
• Abatzoglou T., Mendel J. Constrained total least squares. — 1987. — Т. 12. — С. 1485—1488. — (Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP’87)).
• de Groen P. arxiv.org An introduction to total least squares. — 1996. — С. 237—253. — (Nieuw Archief voor Wiskunde, Vierde serie, deel 14).
• Perpendicular Regression Of A Line at MathPages
• Amiri-Simkooei A. R., Jazaeri S. Weighted total least squares formulated by standard least squares theory // Journal of Geodetic Science. — 2012. — Т. 2 (2). — С. 113—124.