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

Метод наименьших квадратов(МНК) - математический метод, по нахождению приближающей функции - $f(x)$ по набору данных(точек) $\{(x_{1},y_{1}),(x_{2},y_{2}),\ldots ,(x_{n},y_{n})\}$, которая минимизирует сумму квадратов отклонений точек от найденной функции. 

Приближающая функция - $f(x)$

Сущность МНК

Требуется по значениям $x_{1},...,x_{n}$ и $y_{1},...,y_{n}$как можно точнее оценить $f(x)$, точнее означает с минимальными ошибками, с минимальным разбросом(дисперсией), однако зависимость $f(x)$ не должна быть просто построенной по точкам $(x_{i},y_{i})$.
Из Модели Регрессии следует минимизировать дисперсию ошибок:

$D[{\boldsymbol {\varepsilon }}]={\frac {1}{n}}\sum _{i=1}^{n}(y_{i}-f({\vec {x}}_{i},{\vec {c}}))^{2}$

где ${\vec {c}}=\{c_{0},..,c_{m}\}$ — вектор неизвестных параметров.

Далее будем рассматривать только случай однопараметрического (однофакторного) МНК, когда приближающая функция $f({\vec {x}}_{i},{\vec {c}})$ зависит только от одного параметра $x$, (т.е приближающая функция зависит только от одной переменной), поэтому считаем что дисперсия имеет вид:

$D[{\boldsymbol {\varepsilon }}]={\frac {1}{n}}\sum _{i=1}^{n}(y_{i}-f(x_{i},{\vec {c}}))^{2}$

Минимум $D[{\boldsymbol {\varepsilon }}]$ будем искать как:

${\frac {\partial D[{\boldsymbol {\varepsilon }}]}{\partial {\vec {c}}}}=0$

${\frac {\partial D[{\boldsymbol {\varepsilon }}]}{\partial {\vec {c}}}}={\frac {\partial D[{\boldsymbol {\varepsilon }}]}{\partial c_{0}}}{\vec {e}}_{0}+\ldots +{\frac {\partial D[{\boldsymbol {\varepsilon }}]}{\partial c_{m}}}{\vec {e}}_{m}=0$

где $\{{\vec {e}}_{k}\}_{k=0}^{m}=\{{\vec {e}}_{0},..,{\vec {e}}_{m}\}$ — базис из линейно-независимых функций.

Минимизируя дисперсию ошибок $D[{\boldsymbol {\varepsilon }}]$ по неизвестным параметрам $\{c_{0},..,c_{m}\}$, базисом может являться система линейно-независимых функций $\{{\vec {e}}_{k}\}_{k=0}^{m}=\{\varphi _{k}\}_{k=0}^{m}=\{\varphi _{0}(x),..,\varphi _{m}(x)\}$, при условии $m<n$, тогда приближающая функция есть разложение по этому базису: $f(x,{\vec {c}})=\sum _{k=0}^{m}c_{k}\varphi _{k}(x)$, а функция дисперсии ошибки будет иметь вид:

$D[{\boldsymbol {\varepsilon }}]={\frac {1}{n}}\sum _{i=1}^{n}(y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i}))^{2}$

Общий случай вычисления коэффициентов приближающей функции

---

Для этого дифференцируем $D[{\boldsymbol {\varepsilon }}]$(переобозначение: $D[{\boldsymbol {\varepsilon }}]\equiv D$) отдельно по каждому из параметров $c_{j}$ из $\{c_{k}\}_{k=0}^{m}$, где $j$ пробегает значения от $0$ до $m$:
${\frac {\partial D}{\partial c_{j}}}=0$, где $j=0..m$
Получаем систему уравнений относительно параметров $\{c_{k}\}_{k=0}^{m}$:
${\begin{cases}{\frac {\partial D}{\partial c_{0}}}=0\\{\frac {\partial D}{\partial c_{1}}}=0\\\ldots \\{\frac {\partial D}{\partial c_{m}}}=0\end{cases}}$

${\frac {\partial D}{\partial c_{j}}}={\frac {\partial }{\partial c_{j}}}\sum _{i=1}^{n}(y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i}))^{2}=-2\sum _{i=1}^{n}(y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i}))\sum _{k=0}^{m}{\frac {\partial c_{k}}{\partial c_{j}}}\varphi _{k}(x_{i})=0$

${\frac {\partial c_{k}}{\partial c_{j}}}={\begin{cases}1,j=k\\0,j\neq k\end{cases}}\equiv \delta _{jk}$

$-2\sum _{i=1}^{n}(y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i}))\sum _{k=0}^{m}\delta _{jk}\varphi _{k}(x_{i})=0$

Используя:

$\sum _{k=0}^{m}\delta _{jk}\varphi _{k}(x)=\varphi _{j}(x)$

приходим к:

$\sum _{i=1}^{n}[y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})]\varphi _{j}(x_{i})=0$,   где $j=0\ldots m$

Распишем предыдущую формулу в виде системы $m+1$ уравнений:

${\begin{cases}\sum _{i=1}^{n}[y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})]\varphi _{0}(x_{i})=0\\\sum _{i=1}^{n}[y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})]\varphi _{1}(x_{i})=0\\\dots \\\sum _{i=1}^{n}[y_{i}-\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})]\varphi _{m}(x_{i})=0\\\end{cases}}$

Процесс приведения системы уравнений к матричному виду:

Раскроем скобки и перенесём $y_{i}\varphi _{0}(x_{i})$ вправо: ${\begin{cases}\sum _{i=1}^{n}\varphi _{0}(x_{i})\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})=\sum _{i=1}^{n}y_{i}\varphi _{0}(x_{i})\\\sum _{i=1}^{n}\varphi _{1}(x_{i})\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})=\sum _{i=1}^{n}y_{i}\varphi _{1}(x_{i})\\\dots \\\sum _{i=1}^{n}\varphi _{m}(x_{i})\sum _{k=0}^{m}c_{k}\varphi _{k}(x_{i})=\sum _{i=1}^{n}y_{i}\varphi _{m}(x_{i})\\\end{cases}}$ В левой части, раскроем сумму по $k$ и перемножим: ${\begin{cases}\sum _{i=1}^{n}[c_{0}\varphi _{0}^{2}(x_{i})+c_{1}\varphi _{0}(x_{i})\varphi _{1}(x_{i})+...+c_{m}\varphi _{0}(x_{i})\varphi _{m}(x_{i})]=\sum _{i=1}^{n}y_{i}\varphi _{0}(x_{i})\\\sum _{i=1}^{n}[c_{0}\varphi _{1}(x_{i})\varphi _{0}(x_{i})+c_{1}\varphi _{1}^{2}(x_{i})+...+c_{m}\varphi _{1}(x_{i})\varphi _{m}(x_{i})]=\sum _{i=1}^{n}y_{i}\varphi _{1}(x_{i})\\\dots \\\sum _{i=1}^{n}[c_{0}\varphi _{m}(x_{i})\varphi _{0}(x_{i})+c_{1}\varphi _{m}(x_{i})\varphi _{1}(x_{i})+...+c_{m}\varphi _{m}^{2}(x_{i})]=\sum _{i=1}^{n}y_{i}\varphi _{m}(x_{i})\\\end{cases}}$

Используя следующие переобозначения:

$\langle \varphi _{j},\varphi _{k}\rangle =\sum _{i=1}^{n}\varphi _{j}(x_{i})\varphi _{k}(x_{i})$,  

$\langle y,\varphi _{j}\rangle =\sum _{i=1}^{n}y_{i}\varphi _{j}(x_{i})$

Запишем предыдущую систему уравнений в матричном виде: $A{\vec {c}}={\vec {b}}$,

где $A={\begin{pmatrix}\langle \varphi _{0},\varphi _{0}\rangle &\langle \varphi _{1},\varphi _{0}\rangle &\cdots &\langle \varphi _{m},\varphi _{0}\rangle \\\langle \varphi _{0},\varphi _{1}\rangle &\langle \varphi _{1},\varphi _{1}\rangle &\cdots &\langle \varphi _{m},\varphi _{1}\rangle \\\vdots &\vdots &\ddots &\vdots \\\langle \varphi _{0},\varphi _{m}\rangle &\langle \varphi _{1},\varphi _{m}\rangle &\cdots &\langle \varphi _{m},\varphi _{m}\rangle \\\end{pmatrix}}$, ${\vec {c}}={\begin{pmatrix}c_{0}\\c_{1}\\\vdots \\c_{m}\end{pmatrix}}$, ${\vec {b}}={\begin{pmatrix}\langle y,\varphi _{0}\rangle \\\langle y,\varphi _{1}\rangle \\\vdots \\\langle y,\varphi _{m}\rangle \\\end{pmatrix}}$

Матрица $A$ - называется матрицей Грама.

решаем систему матричным методом и находим вектор искомых коэффициентов приближающей функции: ${\vec {c}}=\{c_{0}\;,c_{1},\;\ldots ,\;c_{m}\;\}$${\vec {c}}={\begin{pmatrix}c_{0}\\c_{1}\\\vdots \\c_{m}\end{pmatrix}}=A^{-1}b={\begin{pmatrix}\langle \varphi _{0},\varphi _{0}\rangle &\langle \varphi _{1},\varphi _{0}\rangle &\cdots &\langle \varphi _{m},\varphi _{0}\rangle \\\langle \varphi _{0},\varphi _{1}\rangle &\langle \varphi _{1},\varphi _{1}\rangle &\cdots &\langle \varphi _{m},\varphi _{1}\rangle \\\vdots &\vdots &\ddots &\vdots \\\langle \varphi _{0},\varphi _{m}\rangle &\langle \varphi _{1},\varphi _{m}\rangle &\cdots &\langle \varphi _{m},\varphi _{m}\rangle \\\end{pmatrix}}^{-1}{\begin{pmatrix}\langle y,\varphi _{0}\rangle \\\langle y,\varphi _{1}\rangle \\\vdots \\\langle y,\varphi _{m}\rangle \\\end{pmatrix}}$

В результате, мы нашли приближающую функцию $f(x)$:

$f(x)={\begin{pmatrix}c_{0}\\c_{1}\\\vdots \\c_{m}\end{pmatrix}}{\begin{pmatrix}\varphi _{0}(x),&\varphi _{1}(x),&\ldots ,&\varphi _{m}(x)\end{pmatrix}}=c_{0}\varphi _{0}(x)+c_{1}\varphi _{1}(x)+...+c_{m}\varphi _{m}(x)$

Основные случаи

Степенное разложение(общий случай)

Часто в виде системы линейно-независимых функций $\{\varphi _{k}(x)\}_{k=0}^{m}$ выбирают набор степенных функций $\{1,x,..,x^{m}\}$ и приближающую функцию ищут в виде:

$f(x,{\vec {c}})=c_{0}+c_{1}x+c_{2}x^{2}+...+c_{m}x^{m}=\sum _{k=0}^{m}c_{k}x^{k}$

Тогда система уравнений $A{\vec {c}}={\vec {b}}$ будет иметь следующий вид:

(далее производится суммирование $\sum _{i=1}^{n}$ которое просто обозначаем $\sum$ )

где $A=$ ${\begin{pmatrix}\sum 1&\sum x_{i}&\cdots &\sum x_{i}^{m}\\\sum x_{i}&\sum x_{i}^{2}&\cdots &\sum x_{i}^{m+1}\\\vdots &\vdots &\ddots &\vdots \\\sum x_{i}^{m}&\sum x_{i}^{m+1}&\cdots &\sum x_{i}^{2m}\end{pmatrix}}$,    $b={\begin{pmatrix}\sum y_{i}\\\sum y_{i}x_{i}\\\vdots \\\sum y_{i}x_{i}^{m}\end{pmatrix}}$

Вектор коэффициентов: ${\vec {c}}={\begin{pmatrix}c_{0}\\c_{1}\\\vdots \\c_{m}\end{pmatrix}}=A^{-1}b={\begin{pmatrix}\sum 1&\sum x_{i}&\cdots &\sum x_{i}^{m}\\\sum x_{i}&\sum x_{i}^{2}&\cdots &\sum x_{i}^{m+1}\\\vdots &\vdots &\ddots &\vdots \\\sum x_{i}^{m}&\sum x_{i}^{m+1}&\cdots &\sum x_{i}^{2m}\end{pmatrix}}^{-1}{\begin{pmatrix}\sum y_{i}\\\sum y_{i}x_{i}\\\vdots \\\sum y_{i}x_{i}^{m}\end{pmatrix}}$

$f(x)={\begin{pmatrix}c_{0}\\c_{1}\\\vdots \\c_{m}\end{pmatrix}}\{1,x,\ldots ,x^{m}\}=c_{0}+c_{1}x+c_{2}x^{2}+...+c_{m}x^{m}$

На практике применима до степени $m\leq 5$ из-за быстрого нарастания ошибки при вычислении матрицы $A^{-1}$.

Приближение линейной зависимостью

Набор данных(точки) $\{x_{i},y_{i}\}_{i=1}^{n}$, приближается линейной зависимостью(прямой):

$f(x)=c_{0}+c_{1}x$

Базисные функции линейной зависимости: $\varphi _{0}(x)=1,\ \varphi _{1}(x)=x$

(далее производится суммирование $\sum _{i=1}^{n}$)

$A={\begin{pmatrix}\langle \varphi _{0},\varphi _{0}\rangle &\langle \varphi _{1},\varphi _{0}\rangle \\\langle \varphi _{0},\varphi _{1}\rangle &\langle \varphi _{1},\varphi _{1}\rangle \\\end{pmatrix}}={\begin{pmatrix}\sum 1&\sum x_{i}\\\sum x_{i}&\sum x_{i}^{2}\\\end{pmatrix}}$,

$b={\begin{pmatrix}\langle y,\varphi _{0}\rangle \\\langle y,\varphi _{1}\rangle \end{pmatrix}}={\begin{pmatrix}\sum y_{i}\\\sum y_{i}x_{i}\end{pmatrix}}$

$\sum _{i=1}^{n}1=n$

${\vec {c}}=A^{-1}{\vec {b}}={\begin{pmatrix}n&\sum x_{i}\\\sum x_{i}&\sum x_{i}^{2}\\\end{pmatrix}}^{-1}\cdot {\begin{pmatrix}\sum y_{i}\\\sum y_{i}x_{i}\end{pmatrix}}$

$A^{-1}={\frac {1}{n\sum x_{i}^{2}-(\sum x_{i})^{2}}}{\begin{pmatrix}\sum x_{i}^{2}&-\sum x_{i}\\-\sum x_{i}&n\\\end{pmatrix}}$

${\vec {c}}={\begin{pmatrix}c_{0}\\c_{1}\end{pmatrix}}=A^{-1}{\vec {b}}={\frac {1}{n\sum x_{i}^{2}-(\sum x_{i})^{2}}}{\begin{pmatrix}\sum x_{i}^{2}&-\sum x_{i}\\-\sum x_{i}&n\\\end{pmatrix}}\cdot {\begin{pmatrix}\sum y_{i}\\\sum y_{i}x_{i}\end{pmatrix}}$

${\vec {c}}={\begin{pmatrix}c_{0}\\c_{1}\end{pmatrix}}={\frac {1}{n\sum x_{i}^{2}-(\sum x_{i})^{2}}}{\begin{pmatrix}\sum x_{i}^{2}\sum y_{i}-\sum x_{i}\sum y_{i}x_{i}\\n\sum y_{i}x_{i}-\sum x_{i}\sum y_{i}\end{pmatrix}}$

Таким образом вычислили коэффициенты и нашли приближающую функцию:

$c_{0}={\frac {\sum x_{i}^{2}\sum y_{i}-\sum x_{i}\sum y_{i}x_{i}}{n\sum x_{i}^{2}-(\sum x_{i})^{2}}}$

$c_{1}={\frac {n\sum y_{i}x_{i}-\sum x_{i}\sum y_{i}}{n\sum x_{i}^{2}-(\sum x_{i})^{2}}}$

$f(x)=c_{0}+c_{1}x$

Ортогональность линейно-независимой системы функций

Чтобы упростить решение системы уравнений $A{\vec {c}}={\vec {b}}$ находя ${\vec {c}}$ матричным методом ${\vec {c}}=A^{-1}{\vec {b}}$, нужно матрицу $A$ привести к диагональному виду, для этого необходимо чтобы базисные функции $\{\varphi _{k}(x)\}_{k=0}^{m}$ по которым разложена приближающая функция $f(x)=\sum _{k=0}^{m}c_{k}\varphi _{k}(x)$ были ортогональны на определённом интервале, например на $[x_{1},x_{n}]$.

То есть их скалярное произведение: $\langle \varphi _{j},\varphi _{k}\rangle =\sum _{i=1}^{n}\varphi _{j}(x_{i})\varphi _{k}(x_{i})$

обладало бы свойством ортогональности: $\langle \varphi _{j},\varphi _{k}\rangle ={\begin{cases}\neq 0,\ j=k\\\ \ 0,\ j\neq k\end{cases}}$ на интервале $[x_{1},x_{n}]$

Тогда матрица Грама станет диагональной:

$A={\begin{pmatrix}\langle \varphi _{0},\varphi _{0}\rangle &0&\cdots &0\\0&\langle \varphi _{1},\varphi _{1}\rangle &\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &\langle \varphi _{m},\varphi _{m}\rangle \\\end{pmatrix}}$

а набор коэффициентов ${\vec {c}}=\{c_{k}\}_{k=0}^{m}$ может быть легко вычислен:

${\vec {c}}=A^{-1}{\vec {b}}={\begin{pmatrix}{\frac {1}{\langle \varphi _{0},\varphi _{0}\rangle }}&0&\cdots &0\\0&{\frac {1}{\langle \varphi _{1},\varphi _{1}\rangle }}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &{\frac {1}{\langle \varphi _{m},\varphi _{m}\rangle }}\\\end{pmatrix}}\cdot {\begin{pmatrix}\langle y,\varphi _{0}\rangle \\\langle y,\varphi _{1}\rangle \\\vdots \\\langle y,\varphi _{m}\rangle \\\end{pmatrix}}={\begin{pmatrix}{\frac {\langle y,\varphi _{0}\rangle }{\langle \varphi _{0},\varphi _{0}\rangle }}\\{\frac {\langle y,\varphi _{1}\rangle }{\langle \varphi _{1},\varphi _{1}\rangle }}\\\vdots \\{\frac {\langle y,\varphi _{m}\rangle }{\langle \varphi _{m},\varphi _{m}\rangle }}\\\end{pmatrix}}$

В данном случае, коэффициенты $\{c_{k}\}_{k=0}^{m}$ приближающей функции $f(x)$ называется коэффициентами Фурье: 

$c_{k}={\frac {\langle y,\varphi _{k}\rangle }{\langle \varphi _{k},\varphi _{k}\rangle }}$

А сама приближающая функция $f(x)$ - обобщённым многочленом Фурье:

$f(x)=\sum _{k=0}^{m}{\frac {\langle y,\varphi _{k}\rangle }{\langle \varphi _{k},\varphi _{k}\rangle }}\varphi _{k}(x)$

Однако, ортогональные функции, назовём их $\{e_{k}\}_{k=0}^{m}$ (чтобы отличать от первоначально заданного линейно-независимых функций $\{\varphi _{k}\}_{k=0}^{m}$) заранее неизвестны и определяются относительно того или иного скалярного произведения $\langle \varphi _{k},\varphi _{j}\rangle$.

Находятся $\{e_{k}\}_{k=0}^{m}$ методом ортогонализации линейной-независимого набора векторов $\{\varphi _{k}\}_{k=0}^{m}$ (процесс Грама-Шмидта):

$e_{0}=\varphi _{0}$

$e_{1}=\varphi _{1}-{\frac {\langle \varphi _{1},e_{0}\rangle }{\|e_{0}\|^{2}}}e_{0}$

$\ldots$

$e_{m}=\varphi _{m}-\sum _{k=0}^{m-1}{\frac {\langle \varphi _{m},e_{k}\rangle }{\|e_{k}\|^{2}}}e_{k}$

где $\|e_{k}\|^{2}=\langle e_{k},e_{k}\rangle$

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

$f(x)=\sum _{k=0}^{m}c_{k}e_{k}=\sum _{k=0}^{m}{\frac {\langle y,e_{k}\rangle }{\|e_{k}\|^{2}}}e_{k}$

где $c_{k}={\frac {\langle y,e_{k}\rangle }{\langle e_{k},e_{k}\rangle }}={\frac {\langle y,e_{k}\rangle }{\|e_{k}\|^{2}}}$

Однако, такой метод нахождение базисных функции $\{e_{k}\}_{k=0}^{m}$ является численно неустойчивым, так как происходит накопление ошибки включающей в себя вычисления предыдущих функций.

Для того чтобы избежать накопления ошибки при вычислениях базисных функций $\{e_{k}\}_{k=0}^{m}$, нужно воспользоваться рекуррентным соотношением ортогональных многочленов:

${e_{k+1}\ =\ (\varphi _{1}-\beta _{k})\ e_{k}\ -\ \alpha _{k}\ e_{k-1}},$

где

$\beta _{n}={\frac {\langle \varphi _{1}e_{k},e_{k}\rangle }{\langle e_{k},e_{k}\rangle }},\qquad \alpha _{k}={\frac {\langle \varphi _{1}e_{k},e_{k-1}\rangle }{\langle e_{k-1},e_{k-1}\rangle }}$.

Для случая когда приближение производится степенными функциями, т.е. $\{\varphi _{k}\}_{k=0}^{m}=\{x^{k}\}_{k=0}^{m}$ рекуррентное соотношение имеет вид:

$e_{k+1}=xe_{k}-{\frac {\langle xe_{k},e_{k}\rangle }{\langle e_{k},e_{k}\rangle }}e_{k}-{\frac {\langle xe_{k},e_{k-1}\rangle }{\langle e_{k-1},e_{k-1}\rangle }}e_{k-1}$

где

$k=1,2,..;\ \ \ \ e_{0}=1,\ \ \ \ e_{1}=x-{\frac {\langle x,e_{0}\rangle }{\langle e_{0},e_{0}\rangle }}e_{0}$

Иногда приведённые выше многочлены носят названия многочленов Чебышева, но не стоит их путать с классическими многочленами Чебышева.

Рекуррентное соотношение может быть упрощено, если область определения функции $y$, $[x_{1},x_{n}]$ отобразить на интервал симметричный относительно нуля $[-1,1]$:

$t:[x_{1},x_{n}]\rightarrow [-1,1]$

Примером отображения $t$ , для равноотстоящих точек, может быть функция:

$t(x)={\frac {2x-(x_{n}+x_{1})}{x_{n}-x_{1}}}$

Тогда рекуррентное соотношение будет иметь более простой вид :

$e_{k+1}(t)\ =\ t\ e_{k}(t)\ -\ \alpha _{k}\ e_{k-1}(t)$

где

$e_{0}=1;\ \ \ \ e_{1}=t;\ \ \ \ \ \alpha _{k}={\frac {\langle te_{k},e_{k-1}\rangle }{\langle e_{k-1},e_{k-1}\rangle }};$

Используемые материалы

1. "Вывод рекуррентного соотношения ортогональных многочленов из процесса ортогонализации Грама-Шмидта, а также схема применения полученного рекуррентного соотношения" Сухопаров С.Ю.[1]

2. http://solidbase.karelia.ru/edu/meth_calc/files/09.shtm «Аппроксимация функций методом наименьших квадратов»