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

Метод наименьших квадратов (МНК, англ. Ordinary Least Squares, OLS ) - математический метод, применяемый для решения различных задач, основанный на минимизации суммы квадратов отклонений некоторых функций от искомых переменных. Он может использоваться для «решения» переопределенных систем уравнений (когда количество уравнений превышает количество неизвестных), для поиска решения в случае обычных (не переопределенных) нелинейных систем уравнений, для аппроксимации точечных значений некоторой функции. МНК является одним из базовых методов регрессионного анализа для оценки неизвестных параметров регрессионных моделей по выборочным данным.

Энциклопедичный YouTube

1 / 5

✪ Метод наименьших квадратов. Тема

✪ Метод наименьших квадратов, урок 1/2. Линейная функция

✪ Эконометрика. Лекция 5 .Метод наименьших квадратов

✪ Митин И. В. - Обработка результатов физ. эксперимента - Метод наименьших квадратов (Лекция 4)

✪ Эконометрика: Суть метода наименьших квадратов #2

Субтитры

История

До начала XIX в. учёные не имели определённых правил для решения системы уравнений , в которой число неизвестных меньше, чем число уравнений; до этого времени употреблялись частные приёмы, зависевшие от вида уравнений и от остроумия вычислителей, и потому разные вычислители, исходя из тех же данных наблюдений, приходили к различным выводам. Гауссу (1795) принадлежит первое применение метода, а Лежандр (1805) независимо открыл и опубликовал его под современным названием (фр. Méthode des moindres quarrés ) . Лаплас связал метод с теорией вероятностей , а американский математик Эдрейн (1808) рассмотрел его теоретико-вероятностные приложения . Метод распространён и усовершенствован дальнейшими изысканиями Энке , Бесселя , Ганзена и других.

Сущность метода наименьших квадратов

Пусть x {\displaystyle x} - набор n {\displaystyle n} неизвестных переменных (параметров), f i (x) {\displaystyle f_{i}(x)} , , m > n {\displaystyle m>n} - совокупность функций от этого набора переменных. Задача заключается в подборе таких значений x {\displaystyle x} , чтобы значения этих функций были максимально близки к некоторым значениям y i {\displaystyle y_{i}} . По существу речь идет о «решении» переопределенной системы уравнений f i (x) = y i {\displaystyle f_{i}(x)=y_{i}} , i = 1 , … , m {\displaystyle i=1,\ldots ,m} в указанном смысле максимальной близости левой и правой частей системы. Сущность МНК заключается в выборе в качестве «меры близости» суммы квадратов отклонений левых и правых частей | f i (x) − y i | {\displaystyle |f_{i}(x)-y_{i}|} . Таким образом, сущность МНК может быть выражена следующим образом:

В случае, если система уравнений имеет решение, то минимум суммы квадратов будет равен нулю и могут быть найдены точные решения системы уравнений аналитически или, например, различными численными методами оптимизации. Если система переопределена, то есть, говоря нестрого, количество независимых уравнений больше количества искомых переменных, то система не имеет точного решения и метод наименьших квадратов позволяет найти некоторый «оптимальный» вектор x {\displaystyle x} в смысле максимальной близости векторов y {\displaystyle y} и f (x) {\displaystyle f(x)} или максимальной близости вектора отклонений e {\displaystyle e} к нулю (близость понимается в смысле евклидова расстояния).

Пример - система линейных уравнений

В частности, метод наименьших квадратов может использоваться для «решения» системы линейных уравнений

где A {\displaystyle A} прямоугольная матрица размера m × n , m > n {\displaystyle m\times n,m>n} (т.е. число строк матрицы A больше количества искомых переменных).

Такая система уравнений в общем случае не имеет решения. Поэтому эту систему можно «решить» только в смысле выбора такого вектора x {\displaystyle x} , чтобы минимизировать «расстояние» между векторами A x {\displaystyle Ax} и b {\displaystyle b} . Для этого можно применить критерий минимизации суммы квадратов разностей левой и правой частей уравнений системы, то есть (A x − b) T (A x − b) → min x {\displaystyle (Ax-b)^{T}(Ax-b)\rightarrow \min _{x}} . Нетрудно показать, что решение этой задачи минимизации приводит к решению следующей системы уравнений

МНК в регрессионном анализе (аппроксимация данных)

Пусть имеется n {\displaystyle n} значений некоторой переменной y {\displaystyle y} (это могут быть результаты наблюдений, экспериментов и т. д.) и соответствующих переменных x {\displaystyle x} . Задача заключается в том, чтобы взаимосвязь между y {\displaystyle y} и x {\displaystyle x} аппроксимировать некоторой функцией , известной с точностью до некоторых неизвестных параметров b {\displaystyle b} , то есть фактически найти наилучшие значения параметров b {\displaystyle b} , максимально приближающие значения f (x , b) {\displaystyle f(x,b)} к фактическим значениям y {\displaystyle y} . Фактически это сводится к случаю «решения» переопределенной системы уравнений относительно b {\displaystyle b} :

F (x t , b) = y t , t = 1 , … , n {\displaystyle f(x_{t},b)=y_{t},t=1,\ldots ,n} .

В регрессионном анализе и в частности в эконометрике используются вероятностные модели зависимости между переменными

Y t = f (x t , b) + ε t {\displaystyle y_{t}=f(x_{t},b)+\varepsilon _{t}} ,

где ε t {\displaystyle \varepsilon _{t}} - так называемые случайные ошибки модели.

Соответственно, отклонения наблюдаемых значений y {\displaystyle y} от модельных f (x , b) {\displaystyle f(x,b)} предполагается уже в самой модели. Сущность МНК (обычного, классического) заключается в том, чтобы найти такие параметры b {\displaystyle b} , при которых сумма квадратов отклонений (ошибок, для регрессионных моделей их часто называют остатками регрессии) e t {\displaystyle e_{t}} будет минимальной:

где R S S {\displaystyle RSS} - англ. Residual Sum of Squares определяется как:

В общем случае решение этой задачи может осуществляться численными методами оптимизации (минимизации). В этом случае говорят о нелинейном МНК (NLS или NLLS - англ. Non-Linear Least Squares ). Во многих случаях можно получить аналитическое решение. Для решения задачи минимизации необходимо найти стационарные точки функции R S S (b) {\displaystyle RSS(b)} , продифференцировав её по неизвестным параметрам b {\displaystyle b} , приравняв производные к нулю и решив полученную систему уравнений:

МНК в случае линейной регрессии

Пусть регрессионная зависимость является линейной:

Пусть y - вектор-столбец наблюдений объясняемой переменной, а X {\displaystyle X} - это (n × k) {\displaystyle ({n\times k})} -матрица наблюдений факторов (строки матрицы - векторы значений факторов в данном наблюдении, по столбцам - вектор значений данного фактора во всех наблюдениях). Матричное представление линейной модели имеет вид:

Тогда вектор оценок объясняемой переменной и вектор остатков регрессии будут равны

соответственно сумма квадратов остатков регрессии будет равна

Дифференцируя эту функцию по вектору параметров b {\displaystyle b} и приравняв производные к нулю, получим систему уравнений (в матричной форме):

В расшифрованной матричной форме эта система уравнений выглядит следующим образом:

(∑ x t 1 2 ∑ x t 1 x t 2 ∑ x t 1 x t 3 … ∑ x t 1 x t k ∑ x t 2 x t 1 ∑ x t 2 2 ∑ x t 2 x t 3 … ∑ x t 2 x t k ∑ x t 3 x t 1 ∑ x t 3 x t 2 ∑ x t 3 2 … ∑ x t 3 x t k ⋮ ⋮ ⋮ ⋱ ⋮ ∑ x t k x t 1 ∑ x t k x t 2 ∑ x t k x t 3 … ∑ x t k 2) (b 1 b 2 b 3 ⋮ b k) = (∑ x t 1 y t ∑ x t 2 y t ∑ x t 3 y t ⋮ ∑ x t k y t) , {\displaystyle {\begin{pmatrix}\sum x_{t1}^{2}&\sum x_{t1}x_{t2}&\sum x_{t1}x_{t3}&\ldots &\sum x_{t1}x_{tk}\\\sum x_{t2}x_{t1}&\sum x_{t2}^{2}&\sum x_{t2}x_{t3}&\ldots &\sum x_{t2}x_{tk}\\\sum x_{t3}x_{t1}&\sum x_{t3}x_{t2}&\sum x_{t3}^{2}&\ldots &\sum x_{t3}x_{tk}\\\vdots &\vdots &\vdots &\ddots &\vdots \\\sum x_{tk}x_{t1}&\sum x_{tk}x_{t2}&\sum x_{tk}x_{t3}&\ldots &\sum x_{tk}^{2}\\\end{pmatrix}}{\begin{pmatrix}b_{1}\\b_{2}\\b_{3}\\\vdots \\b_{k}\\\end{pmatrix}}={\begin{pmatrix}\sum x_{t1}y_{t}\\\sum x_{t2}y_{t}\\\sum x_{t3}y_{t}\\\vdots \\\sum x_{tk}y_{t}\\\end{pmatrix}},} где все суммы берутся по всем допустимым значениям t {\displaystyle t} .

Если в модель включена константа (как обычно), то x t 1 = 1 {\displaystyle x_{t1}=1} при всех t {\displaystyle t} , поэтому в левом верхнем углу матрицы системы уравнений находится количество наблюдений n {\displaystyle n} , а в остальных элементах первой строки и первого столбца - просто суммы значений переменных: ∑ x t j {\displaystyle \sum x_{tj}} и первый элемент правой части системы - ∑ y t {\displaystyle \sum y_{t}} .

Решение этой системы уравнений и дает общую формулу МНК-оценок для линейной модели:

Для аналитических целей оказывается полезным последнее представление этой формулы (в системе уравнений при делении на n, вместо сумм фигурируют средние арифметические). Если в регрессионной модели данные центрированы , то в этом представлении первая матрица имеет смысл выборочной ковариационной матрицы факторов, а вторая - вектор ковариаций факторов с зависимой переменной. Если кроме того данные ещё и нормированы на СКО (то есть в конечном итоге стандартизированы ), то первая матрица имеет смысл выборочной корреляционной матрицы факторов, второй вектор - вектора выборочных корреляций факторов с зависимой переменной.

Немаловажное свойство МНК-оценок для моделей с константой - линия построенной регрессии проходит через центр тяжести выборочных данных, то есть выполняется равенство:

В частности, в крайнем случае, когда единственным регрессором является константа, получаем, что МНК-оценка единственного параметра (собственно константы) равна среднему значению объясняемой переменной. То есть среднее арифметическое, известное своими хорошими свойствами из законов больших чисел, также является МНК-оценкой - удовлетворяет критерию минимума суммы квадратов отклонений от неё.

Простейшие частные случаи

В случае парной линейной регрессии y t = a + b x t + ε t {\displaystyle y_{t}=a+bx_{t}+\varepsilon _{t}} , когда оценивается линейная зависимость одной переменной от другой, формулы расчета упрощаются (можно обойтись без матричной алгебры). Система уравнений имеет вид:

Отсюда несложно найти оценки коэффициентов:

Несмотря на то что в общем случае модели с константой предпочтительней, в некоторых случаях из теоретических соображений известно, что константа a {\displaystyle a} должна быть равна нулю. Например, в физике зависимость между напряжением и силой тока имеет вид U = I ⋅ R {\displaystyle U=I\cdot R} ; замеряя напряжение и силу тока, необходимо оценить сопротивление. В таком случае речь идёт о модели y = b x {\displaystyle y=bx} . В этом случае вместо системы уравнений имеем единственное уравнение

(∑ x t 2) b = ∑ x t y t {\displaystyle \left(\sum x_{t}^{2}\right)b=\sum x_{t}y_{t}} .

Следовательно, формула оценки единственного коэффициента имеет вид

B ^ = ∑ t = 1 n x t y t ∑ t = 1 n x t 2 = x y ¯ x 2 ¯ {\displaystyle {\hat {b}}={\frac {\sum _{t=1}^{n}x_{t}y_{t}}{\sum _{t=1}^{n}x_{t}^{2}}}={\frac {\overline {xy}}{\overline {x^{2}}}}} .

Случай полиномиальной модели

Если данные аппроксимируются полиномиальной функцией регрессии одной переменной f (x) = b 0 + ∑ i = 1 k b i x i {\displaystyle f(x)=b_{0}+\sum \limits _{i=1}^{k}b_{i}x^{i}} , то, воспринимая степени x i {\displaystyle x^{i}} как независимые факторы для каждого i {\displaystyle i} можно оценить параметры модели исходя из общей формулы оценки параметров линейной модели. Для этого в общую формулу достаточно учесть, что при такой интерпретации x t i x t j = x t i x t j = x t i + j {\displaystyle x_{ti}x_{tj}=x_{t}^{i}x_{t}^{j}=x_{t}^{i+j}} и x t j y t = x t j y t {\displaystyle x_{tj}y_{t}=x_{t}^{j}y_{t}} . Следовательно, матричные уравнения в данном случае примут вид:

(n ∑ n x t … ∑ n x t k ∑ n x t ∑ n x t 2 … ∑ n x t k + 1 ⋮ ⋮ ⋱ ⋮ ∑ n x t k ∑ n x t k + 1 … ∑ n x t 2 k) [ b 0 b 1 ⋮ b k ] = [ ∑ n y t ∑ n x t y t ⋮ ∑ n x t k y t ] . {\displaystyle {\begin{pmatrix}n&\sum \limits _{n}x_{t}&\ldots &\sum \limits _{n}x_{t}^{k}\\\sum \limits _{n}x_{t}&\sum \limits _{n}x_{t}^{2}&\ldots &\sum \limits _{n}x_{t}^{k+1}\\\vdots &\vdots &\ddots &\vdots \\\sum \limits _{n}x_{t}^{k}&\sum \limits _{n}x_{t}^{k+1}&\ldots &\sum \limits _{n}x_{t}^{2k}\end{pmatrix}}{\begin{bmatrix}b_{0}\\b_{1}\\\vdots \\b_{k}\end{bmatrix}}={\begin{bmatrix}\sum \limits _{n}y_{t}\\\sum \limits _{n}x_{t}y_{t}\\\vdots \\\sum \limits _{n}x_{t}^{k}y_{t}\end{bmatrix}}.}

Статистические свойства МНК-оценок

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

1. математическое ожидание случайных ошибок равно нулю, и
2. факторы и случайные ошибки - независимые случайные величины .

Второе условие - условие экзогенности факторов - принципиальное. Если это свойство не выполнено, то можно считать, что практически любые оценки будут крайне неудовлетворительными: они не будут даже состоятельными (то есть даже очень большой объём данных не позволяет получить качественные оценки в этом случае). В классическом случае делается более сильное предположение о детерминированности факторов, в отличие от случайной ошибки, что автоматически означает выполнение условия экзогенности. В общем случае для состоятельности оценок достаточно выполнения условия экзогенности вместе со сходимостью матрицы V x {\displaystyle V_{x}} к некоторой невырожденной матрице при увеличении объёма выборки до бесконечности.

Для того, чтобы кроме состоятельности и несмещенности , оценки (обычного) МНК были ещё и эффективными (наилучшими в классе линейных несмещенных оценок) необходимо выполнение дополнительных свойств случайной ошибки:

Данные предположения можно сформулировать для ковариационной матрицы вектора случайных ошибок V (ε) = σ 2 I {\displaystyle V(\varepsilon)=\sigma ^{2}I} .

Линейная модель, удовлетворяющая таким условиям, называется классической . МНК-оценки для классической линейной регрессии являются несмещёнными , состоятельными и наиболее эффективными оценками в классе всех линейных несмещённых оценок (в англоязычной литературе иногда употребляют аббревиатуру BLUE (Best Linear Unbiased Estimator ) - наилучшая линейная несмещённая оценка; в отечественной литературе чаще приводится теорема Гаусса - Маркова). Как нетрудно показать, ковариационная матрица вектора оценок коэффициентов будет равна:

V (b ^ O L S) = σ 2 (X T X) − 1 {\displaystyle V({\hat {b}}_{OLS})=\sigma ^{2}(X^{T}X)^{-1}} .

Эффективность означает, что эта ковариационная матрица является «минимальной» (любая линейная комбинация коэффициентов, и в частности сами коэффициенты, имеют минимальную дисперсию), то есть в классе линейных несмещенных оценок оценки МНК-наилучшие. Диагональные элементы этой матрицы - дисперсии оценок коэффициентов - важные параметры качества полученных оценок. Однако рассчитать ковариационную матрицу невозможно, поскольку дисперсия случайных ошибок неизвестна. Можно доказать, что несмещённой и состоятельной (для классической линейной модели) оценкой дисперсии случайных ошибок является величина:

S 2 = R S S / (n − k) {\displaystyle s^{2}=RSS/(n-k)} .

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

Необходимо отметить, что если классические предположения не выполнены, МНК-оценки параметров не являются наиболее эффективными и, где W {\displaystyle W} - некоторая симметрическая положительно определенная весовая матрица. Обычный МНК является частным случаем данного подхода, когда весовая матрица пропорциональна единичной матрице. Как известно, для симметрических матриц (или операторов) существует разложение W = P T P {\displaystyle W=P^{T}P} . Следовательно, указанный функционал можно представить следующим образом e T P T P e = (P e) T P e = e ∗ T e ∗ {\displaystyle e^{T}P^{T}Pe=(Pe)^{T}Pe=e_{*}^{T}e_{*}} , то есть этот функционал можно представить как сумму квадратов некоторых преобразованных «остатков». Таким образом, можно выделить класс методов наименьших квадратов - LS-методы (Least Squares).

Доказано (теорема Айткена), что для обобщенной линейной регрессионной модели (в которой на ковариационную матрицу случайных ошибок не налагается никаких ограничений) наиболее эффективными (в классе линейных несмещенных оценок) являются оценки т. н. обобщенного МНК (ОМНК, GLS - Generalized Least Squares) - LS-метода с весовой матрицей, равной обратной ковариационной матрице случайных ошибок: W = V ε − 1 {\displaystyle W=V_{\varepsilon }^{-1}} .

Можно показать, что формула ОМНК-оценок параметров линейной модели имеет вид

B ^ G L S = (X T V − 1 X) − 1 X T V − 1 y {\displaystyle {\hat {b}}_{GLS}=(X^{T}V^{-1}X)^{-1}X^{T}V^{-1}y} .

Ковариационная матрица этих оценок соответственно будет равна

V (b ^ G L S) = (X T V − 1 X) − 1 {\displaystyle V({\hat {b}}_{GLS})=(X^{T}V^{-1}X)^{-1}} .

Фактически сущность ОМНК заключается в определенном (линейном) преобразовании (P) исходных данных и применении обычного МНК к преобразованным данным. Цель этого преобразования - для преобразованных данных случайные ошибки уже удовлетворяют классическим предположениям.

Взвешенный МНК

В случае диагональной весовой матрицы (а значит и ковариационной матрицы случайных ошибок) имеем так называемый взвешенный МНК (WLS - Weighted Least Squares). В данном случае минимизируется взвешенная сумма квадратов остатков модели, то есть каждое наблюдение получает «вес», обратно пропорциональный дисперсии случайной ошибки в данном наблюдении: e T W e = ∑ t = 1 n e t 2 σ t 2 {\displaystyle e^{T}We=\sum _{t=1}^{n}{\frac {e_{t}^{2}}{\sigma _{t}^{2}}}} . Фактически данные преобразуются взвешиванием наблюдений (делением на величину, пропорциональную предполагаемому стандартному отклонению случайных ошибок), а к взвешенным данным применяется обычный МНК.

ISBN 978-5-7749-0473-0 .

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

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

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

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

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

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

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

иллюстрация

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Иллюстрация

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

UPDATE: тут у меня косяк, расстояние до прямой должно измеряться по вертикали, а не ортогональной проекцией. комментатор прав.

Иллюстрация

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

Иллюстрация

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

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

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

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

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

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

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

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

Изначальный коммит доступен . Для минимизации внешних зависимостей я взял код своего софтверного рендерера, уже на хабре. Для решения линейной системы я пользуюсь 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 &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 имеет нулевые компоненты напротив. Это значит, я вешаю пружинку на каждое ребро нашей треугольной сетки: все рёбра стараются получить одну и ту же вершину в качестве отправной и финальной точки.

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

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

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

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 единиц. То есть, мы повесили более сильную пружинку на крайние вершины, решение предпочтёт сильнее растянуть другие. Вот результат:

Давайте вдвое усилим пружинки между вершинами:
nlCoefficient(face[ j ], 2); nlCoefficient(face[(j+1)%3], -2);

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

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

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

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

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

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

Разрежу картинку пополам:

И выделю руками стул:

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

For (int i=0; i

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

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

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

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

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

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

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

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

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

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

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

И его фрагмент, чтобы было виднее:

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

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

3. Аппроксимация функций с помощью метода

наименьших квадратов

Метод наименьших квадратов применяется при обработке результатов эксперимента для аппроксимации (приближения) экспериментальных данных аналитической формулой. Конкретный вид формулы выбирается, как правило, из физических соображений. Такими формулами могут быть:

и другие.

Сущность метода наименьших квадратов состоит в следующем. Пусть результаты измерений представлены таблицей:

где f - известная функция, a 0 , a 1 , …, a m - неизвестные постоянные параметры, значения которых надо найти. В методе наименьших квадратов приближение функции (3.1) к экспериментальной зависимости считается наилучшим, если выполняется условие

то есть сумм a квадратов отклонений искомой аналитической функции от экспериментальной зависимости должна быть минимальна .

Заметим, что функция Q называется невязкой.

Так как невязка

то она имеет минимум. Необходимым условием минимума функции нескольких переменных является равенство нулю всех частных производных этой функции по параметрам. Таким образом, отыскание наилучших значений параметров аппроксимирующей функции (3.1), то есть таких их значений, при которых Q = Q (a 0 , a 1 , …, a m ) минимальна, сводится к решению системы уравнений:

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

Нахождение параметров линейной функции

Пусть экспериментальные данные надо представить линейной функцией:

Требуется подобрать такие значения a и b , для которых функция

будет минимальной. Необходимые условия минимума функции (3.4) сводятся к системе уравнений:

После преобразований получаем систему двух линейных уравнений с двумя неизвестными:

решая которую , находим искомые значения параметров a и b .

Нахождение параметров квадратичной функции

Если аппроксимирующей функцией является квадратичная зависимость

то её параметры a , b , c находят из условия минимума функции:

Условия минимума функции (3.6) сводятся к системе уравнений:

После преобразований получаем систему трёх линейных уравнений с тремя неизвестными:

при решении которой находим искомые значения параметров a , b и c .

Пример . Пусть в результате эксперимента получена следующая таблица значений x и y :

Требуется аппроксимировать экспериментальные данные линейной и квадратичной функциями.

Решение. Отыскание параметров аппроксимирующих функций сводится к решению систем линейных уравнений (3.5) и (3.7). Для решения задачи воспользуемся процессором электронных таблиц Excel .

1. Сначала сцепим листы 1 и 2. Занесём экспериментальные значения x i и y i в столбцы А и В, начиная со второй строки (в первой строке поместим заголовки столбцов). Затем для этих столбцов вычислим суммы и поместим их в десятой строке.

В столбцах C – G разместим соответственно вычисление и суммирование

2. Расцепим листы.Дальнейшие вычисления проведём аналогичным образом для линейной зависимости на Листе 1и для квадратичной зависимости на Листе 2.

3. Под полученной таблицей сформируем матрицу коэффициентов и вектор-столбец свободных членов. Решим систему линейных уравнений по следующему алгоритму:

Для вычисления обратной матрицы и перемножения матриц воспользуемся Мастером функций и функциями МОБР и МУМНОЖ .

4. В блоке ячеек H2: H 9 на основе полученных коэффициентов вычислим значенияаппроксимирующего полинома y i выч ., в блоке I 2: I 9 – отклонения D y i = y i эксп . - y i выч .,в столбце J – невязку:

Полученные таблицы и построенные с помощью Мастера диаграмм графики приведёны на рисунках6, 7, 8.

Рис. 6. Таблица вычисления коэффициентов линейной функции,

аппроксимирующей экспериментальные данные.

Рис. 7. Таблица вычисления коэффициентов квадратичной функции,

аппроксимирующей экспериментальные данные.

Рис. 8. Графическое представление результатов аппроксимации

экспериментальных данных линейной и квадратичной функциями.

Ответ. Аппроксимировали экспериментальные данные линейной зависимостью y = 0,07881 x + 0,442262 c невязкой Q = 0,165167 и квадратичной зависимостью y = 3,115476 x 2 – 5,2175 x + 2,529631 c невязкой Q = 0,002103 .

Задания. Аппроксимировать функцию, заданную таблично, линейной и квадратичной функциями.