Метод наименьших полных квадратов (Bymk; ugnbyu,on] hklud] tfg;jgmkf)

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

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

Линейная модель

Основы

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

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

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

\mathbf {r=y-X{\boldsymbol {\beta }}} .

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

\mathbf {X^{T}WX{\boldsymbol {\beta }}=X^{T}Wy}

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

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

S=\mathbf {r_{x}^{T}M_{x}^{-1}r_{x}+r_{y}^{T}M_{y}^{-1}r_{y}}

,

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

\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 ограничениях. Задача решается с помощью множителей Лагранжа. После некоторых алгебраических преобразований^[3] получим

\mathbf {X^{T}M^{-1}X\Delta {\boldsymbol {\beta }}=X^{T}M^{-1}\Delta y} ,

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

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

\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},\beta )=\alpha +\beta x_{i}\!

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

M_{ii}=\sigma _{y,i}^{2}+\beta ^{2}\sigma _{x,i}^{2}

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

M_{ii}=\sigma _{y,i}^{2}+\left({\frac {dy}{dx}}\right)_{i}^{2}\sigma _{x,i}^{2}

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

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

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

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

XB\approx Y

относительно B, где X — матрица m-на-n, а Y — матрица m-на-k^{[note 2]}

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

\mathrm {argmin} _{R,F}\|[R\;F]\|_{F},\qquad (X+R)B=Y+F

,

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

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

[(X+R)\;(Y+F)]{\begin{bmatrix}B\\-E_{k}\end{bmatrix}}=0.

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

[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)]{\begin{bmatrix}V_{XY}\\V_{YY}\end{bmatrix}}=0.

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

[(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}.

Имплементация в системе 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}$ не была вырожденной, может быть слегка расширен так называемым классическим МНПК алгоритмом^[7].

Вычисление

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

Нелинейная модель

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

\mathbf {J^{T}M^{-1}J\Delta {\boldsymbol {\beta }}=J^{T}M^{-1}\Delta y} .

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

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

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

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

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

См. также

Примечания

↑ Альтернативная форма — $\mathbf {X^{T}WX{\boldsymbol {\Delta }}{\boldsymbol {\beta }}=X^{T}W{\boldsymbol {\Delta }}y}$ , где ${\boldsymbol {\Delta }}{\boldsymbol {\beta }}$ является сдвигом параметров от начальной оценки ${\boldsymbol {\beta }}$ , а ${\boldsymbol {\Delta }}\mathbf {y}$ — разница между y и значением, вычисленным по начальной оценке ${\boldsymbol {\beta }}$
↑ Здесь используется выражение XB ≈ Y, чтобы отразить употреблённые ранее выражения. В литературе чаще используется выражение AX ≈ B, т.е. с буквой X для обозначения n-х-k матрицы неизвестных регриссионных коэффициентов.

↑ Markovsky, Van Huffel, 2007, с. 2283–2302, 2007.
↑ Deming, 1943.
↑ Gans, 1992.
↑ Golub, Van Loan, 1980, с. 883–893.
↑ Golub, Van Loan, 1996, с. 596.
↑ Bjõrck, 1996.
↑ Van Huffel, Vandewalle, 1991.
↑ Van Huffel, 1988.
↑ Van Huffel, 1989, с. 111–119.
↑ Plešinger, 2008, с. 748–770.
↑ Stein.
↑ Samuelson, 1942, с. 80–83.
↑ Ricker, 1975, с. 1494–1498.
↑ Warton, Wright, Falster, Westoby, 2006, с. 259–291.
↑ Draper, Smith, 1998, с. 92–96.
↑ Tofallis, 2002.

Литература

Van Huffel S., Vandewalle J. The Total Least Squares Problems: Computational Aspects and Analysis. — Philadelphia PA: SIAM Publications, 1991. — Т. 9. — (Frontiers in applied mathematics). — ISBN 0-89871-271-0.
Golub G. H., Van Loan C. F. An analysis of the total least squares problem // SIAM J. on Numer. Anal.. — 1980. — Т. 17. — С. 883—893.

Gene H. Golub, Charles F. Van Loan. Matrix Computations. — 3rd. — The Johns Hopkins University Press, 1996.
Ake Bjõrck. Numerical Methods for Least Squares Problems. — SIAM (Society for Industrial and Applied Mathematics), 1996. — ISBN 978-0898713602.
Van Huffel S. Documented Fortran 77 programs of the extended classical total least squares algorithm, the partial singular value decomposition algorithm and the partial total least squares algorithm, Internal Report ESAT-KUL 88/1 ESAT Lab., Dept. of Electrical Engineering,. — Katholieke Universiteit Leuven, 1988.
Van Huffel S. The extended classical total least squares algorithm // J. Comput. Appl. Math.,. — 1989. — С. 111—119,.
Plešinger M. The Total Least Squares Problem and Reduction of Data in AX ≈ B. Doctoral Thesis. — TU of Liberec and Institute of Computer Science, AS CR Prague, 2008. — (Ph.D. Thesis). Архивная копия от 24 июля 2012 на Wayback Machine
Hnětynková I., Plešinger M., Sima D. M., Strakoš Z., Van Huffel S. [1]. — The total least squares problem in AX ≈ B. A new classification with the relationship to the classical works.: SIMAX, 2011. — Т. 32. — С. 748—770.
Yaakov J. Stein. Two Dimensional Euclidean Regression.
Paul A. Samuelson. A Note on Alternative Regressions // Econometrica. — The Econometric Society, 1942. — Т. 10, вып. 1. — С. 80—83. — doi:10.2307/1907024. — JSTOR 1907024.
Ricker W. E. A note concerning Professor Jolicoeur's Comments // Journal of the Fisheries Research Board of Canada. — 1975. — Т. 32. — С. 1494—1498. — doi:10.1139/f75-172.
David I. Warton, Ian J. Wright, Daniel S. Falster, Mark Westoby. Bivariate line-fitting methods for allometry // Biological Reviews. — Wiley, 2006. — Т. 81, вып. 2. — С. 259—291. — doi:10.1017/S1464793106007007.
Draper N.R., Smith H. Applied Regression Analysis. — 3rd edition. — 1998. — С. 92—96. — (Wiley series in probability and statistics). — ISBN 0-471-17982-8.
Chris Tofallis. Model Fitting for Multiple Variables by Minimising the Geometric Mean Deviation // Total Least Squares and Errors-in-Variables Modeling: Analysis, Algorithms and Applications / Sabine Van Huffel, P. Lemmerling. — Dordrecht [u.a.]: Kluwer Academic Publ., 2002. — ISBN 978-1402004766.
Markovsky I., Van Huffel S. Overview of total least squares methods // Signal Processing. — 2007. — Т. 87.
W.E. Deming. Statistical Adjustment of Data. — New York: John Wiley & Sons, 1943.
Peter Gans. Data Fitting in the Chemical Sciences. — Wiley, 1992. — ISBN 9780471934127.

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

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.

[2] Альтернативная форма — $\mathbf {X^{T}WX{\boldsymbol {\Delta }}{\boldsymbol {\beta }}=X^{T}W{\boldsymbol {\Delta }}y}$ , где ${\boldsymbol {\Delta }}{\boldsymbol {\beta }}$ является сдвигом параметров от начальной оценки ${\boldsymbol {\beta }}$ , а ${\boldsymbol {\Delta }}\mathbf {y}$ — разница между y и значением, вычисленным по начальной оценке ${\boldsymbol {\beta }}$

[7] Здесь используется выражение XB ≈ Y, чтобы отразить употреблённые ранее выражения. В литературе чаще используется выражение AX ≈ B, т.е. с буквой X для обозначения n-х-k матрицы неизвестных регриссионных коэффициентов.

[_575aca08a360360c-1] Markovsky, Van Huffel, 2007, с. 2283–2302, 2007.

[_e6e6d2fa0b36b3f2-3] Deming, 1943.

[_dd8644baa04751cd-4] Gans, 1992.

[_40bee0a270f8522f-5] Golub, Van Loan, 1980, с. 883–893.

[_4d029fc1a3a9ebb4-6] Golub, Van Loan, 1996, с. 596.

[_53986ea8c14e3a28-8] Bjõrck, 1996.

[_0c1fd00e3d923831-9] Van Huffel, Vandewalle, 1991.

[_f9c361cdd5e0dfee-10] Van Huffel, 1988.

[_405cabfa4d04af06-11] Van Huffel, 1989, с. 111–119.

[_2d3d73096ddfbef7-12] Plešinger, 2008, с. 748–770.

[_7fab2e9b08ab4546-13] Stein.

[_c345b852ed122880-14] Samuelson, 1942, с. 80–83.

[_172907286966d6f4-15] Ricker, 1975, с. 1494–1498.

[_8428c39132605e92-16] Warton, Wright, Falster, Westoby, 2006, с. 259–291.

[_9a2ae07adb87e7e0-17] Draper, Smith, 1998, с. 92–96.

[_167f66c826adfd3d-18] Tofallis, 2002.

[1]

[note 1]

[2]

[3]

[4]

[5]

[note 2]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]