УказательРазделыОбозначенияАвторО проекте


§

Вспомогательная страница к разделу ИНТЕРПОЛЯЦИЯ.


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

Пусть из физических соображений можно считать (предполагать), что величины x_{} и y_{} связаны линейной зависимостью вида y=kx+b, а неизвестные коэффициенты k_{} и b_{} должны быть оценены экспериментально. Экспериментальные данные представляют собой m>1 точек на координатной плоскости (x_1,y_1), \dots, (x_m,y_m). Если бы эти опыты производились без погрешностей, то подстановка данных в уравнение приводила бы нас к системе из m_{} линейных уравнений для двух неизвестных k_{} и b_{}:

y_1=k\,x_1+b, \dots, y_m=k\,x_m+b \ .

Из любой пары уравнений этой системы можно было бы однозначно определить коэффициенты k_{} и b_{}.

Однако, в реальной жизни опытов без погрешностей не бывает

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

А.М.Б.Розен. Физики шутят. М.Мир.1993.

Будем предполагать, что величины x_{1},\dots,x_m известны точно, а им соответствующие y_1,\dots,y_m — приближенно. Если m>2, то при любых различных x_{i} и x_j пара точек (x_{i},y_i) и (x_{j},y_j) определяет прямую. Но другая пара точек определяет другую прямую, и у нас нет оснований выбрать какую-нибудь одну из всех прямых.

Часто в задаче удаленность точки от прямой измеряют не расстоянием, а разностью ординат k\,x_i+b-y_i, и выбирают прямую так, чтобы сумма квадратов всех таких разностей была минимальна. Коэффициенты k_0 и b_{0} уравнения этой прямой дают некоторое решение стоящей перед нами задачи, которое отнюдь не является решением системы линейных уравнений

k\,x_1+b=y_1,\dots, k\,x_{m}+b=y_m

(вообще говоря, несовместной).

Рассмотрим теперь обобщение предложенной задачи. Пусть искомая зависимость между величинами y_{} и x_{} полиномиальная:

y_1=f(x_1),\dots , y_m=f(x_m), \quad npu \quad f(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1}

Величина \varepsilon_i=f(x_i)-y_i называется i_{}невязкой1). Уравнения

\left\{\begin{array}{ccl} a_0+a_1x_1+\dots+a_{n-1}x_1^{n-1}&=&y_1, \\ a_0+a_1x_2+\dots+a_{n-1}x_2^{n-1}&=&y_2, \\ \dots & & \dots \\ a_0+a_1x_m+\dots+a_{n-1}x_m^{n-1}&=&y_m \end{array} \right.

называются условными.

Предположим что данные интерполяционной таблицы

\begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_m \\ \hline y & y_1 & y_2 &\dots & y_m \end{array}

не являются достоверными: величины x_{} нам известны практически без искажений (т.е. на входе процесса мы имеем абсолютно достоверные данные), а вот измерения величины y_{} подвержены случайным погрешностям.

Задача. Построить полином f_{}(x) такой, чтобы величина

\sum_{j=1}^m [f(x_j)-y_j]^2

стала минимальной.

В случае \deg f_{} =m-1 мы возвращаемся к задаче интерполяции в ее классической постановке. Практический интерес, однако, представляет случай \deg f_{} < m-1, т.е. случай когда экспериментальных данных больше (обычно — много больше) чем значений параметров (коэффициентов полинома f_{}), которые требуется определить.

Так, в случае \deg f_{}=1 речь идет о построении прямой y=ax+b на плоскости (x,y), обеспечивающей

\min (\varepsilon_1^2+\varepsilon_2^2+\dots + \varepsilon_m^2) \, .

Здесь \varepsilon_j = a\,x_j+b-y_j.

Т

Теорема. Если m\ge n_{} и узлы интерполяции x_{1},\dots,x_m все различны, то существует единственный набор коэффициентов a_{0},\dots,a_{n-1}, обеспечивающий минимальное значение для

\sum_{j=1}^m (a_0+a_1x_j+\dots+a_{n-1}x_j^{n-1} -y_j)^2 \ .

Этот набор определяется как решение системы нормальных уравнений

\underbrace{ \left(\begin{array}{llllll} s_0 &s_1&s_2&\ldots&s_{n-2}& s_{n-1}\\ s_1 &s_2&s_3&\ldots&s_{n-1}& s_{n}\\ s_2 &s_3&s_4&\ldots&s_{n}& s_{n+1}\\ \vdots& & & && \vdots\\ s_{n-1} &s_n&s_{n+1}&\ldots &s_{2n-3}&s_{2n-2} \end{array}\right)}_{\displaystyle S_{n\times n}} \left(\begin{array}{l} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{array} \right)= \left(\begin{array}{l} t_0 \\ t_1 \\ t_2 \\ \vdots \\ t_{n-1} \end{array} \right)

при s_{k} = x_1^k+\dots+x_m^k, \ t_{k} = x_1^ky_1+\dots+x_m^k y_m.

§

Если одно из значений x_{j} равно 0_{} , то полагаем 0^{0} = 1, так что s_{0}=m при любых \{x_{1},\dots, x_m\} \subset \mathbb R.

Доказательство. Рассмотрим сумму как полином от неопределенных коэффициентов \{a_{j}\}_{j=0}^{n-1}:

F(a_0,\dots,a_{n-1})= \sum_{i=1}^m [f(x_i)-y_i]^2=
=\sum_{i=1}^m [(a_0+a_1x_i+\dots+a_{n-1}x_i^{n-1})-y_i]^2 \ .

На основании теоремы из пункта ЭКСТРЕМУМЫ ПОЛИНОМА такая функция может принимать экстремальные значения только на вещественных решениях системы уравнений

\partial F / \partial a_0=0, \dots, \partial F / \partial a_{n-1}=0 \ .

Распишем выражение для \partial F / \partial a_j:

\partial F / \partial a_j =\sum_{i=1}^m 2\,[f(x_i)-y_i] \frac{\partial (a_0+a_1x_i+\dots+a_{n-1}x_i^{n-1})}{\partial a_j}
=2\, \sum_{i=1}^m [f(x_i)-y_i] x_i^{j}=
=2\, \sum_{i=1}^m \left[\left(a_0x_i^{j}+a_1x_i^{j+1}+\dots+a_{n-1}x_i^{j+n-1} \right)-y_ix_i^{j} \right]=
= 2\left[a_0\, \sum_{i=1}^m x_i^{j}+a_1\, \sum_{i=1}^m x_i^{j+1}+ \dots+a_{n-1}\, \sum_{i=1}^m x_i^{j+n-1}- \sum_{i=1}^m y_ix_i^{j} \right]=
=2 [a_0\, s_j+a_1\, s_{j+1}+\dots+a_{n-1}\,s_{j+n-1}-t_j ]

Таким образом, условия \left\{\partial F / \partial a_j=0 \right\}_{j=0}^n можно переписать в виде системы нормальных уравнений.

Покажем теперь, что матрица этой системы имеет ненулевой определитель. Действительно, матрица S_{}ганкелева. При m=n_{} ее определитель вычисляется с помощью теоремы Бине-Коши:

\det S = \prod_{1\le j<k\le n} (x_k-x_j)^2 \ .

Воспользуемся той же теоремой и для случая m>n:

S=\left(\begin{array}{ccccc} 1 &1&1&\ldots&1\\ x_1 &x_2&x_3&\ldots&x_{m}\\ \ldots&& & &\ldots\\ x_1^{n-1} &x_2^{n-1}&x_3^{n-1}&\ldots&x_m^{n-1} \end{array}\right) \cdot \left(\begin{array}{ccccc} 1 &x_1&x_1^2&\ldots&x_1^{n-1}\\ 1 &x_2&x_2^2&\ldots&x_2^{n-1}\\ \ldots&& & &\ldots\\ 1 &x_m&x_m^2&\ldots&x_m^{n-1} \end{array}\right)
\det S = \sum_{1\le j_1< j_2 <\dots < j_n \le m} \left|\begin{array}{cccc} 1 &1&\ldots&1\\ x_{j_1} &x_{j_2}&\ldots&x_{j_n}\\ \ldots&& &\ldots\\ x_{j_1}^{n-1} &x_{j_2}^{n-1}&\ldots&x_{j_n}^{n-1} \end{array}\right| \cdot \left|\begin{array}{cccc} 1 &x_{j_1}&\ldots&x_{j_1}^{n-1}\\ 1 &x_{j_2}&\ldots&x_{j_2}^{n-1}\\ \ldots&& &\ldots\\ 1 &x_{j_n}&\ldots&x_{j_n}^{n-1} \end{array}\right|=
=\sum_{1\le j_1< j_2 <\dots < j_n \le m} \ \prod_{1\le L < K \le n} (x_{j_K}-x_{j_L})^2 \ .

Каждое слагаемое неотрицательно и отлично от нуля поскольку, по предположению, все \{x_j\}_{j=1}^m различны. Поэтому \det S >0. По теореме Крамера система нормальных уравнений имеет единственное решение.

?

Показать, что линейный полином y=a_{0}+a_1x, построенный по методу наименьших квадратов, определяет на плоскости (x_{},y) прямую, проходящую через центр тяжести системы точек (x_{1},y_1),\dots,(x_m,y_m).

П

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

\begin{array}{c|cccccc} x & 0.5 & 1 & 1.5 & 2 & 2.5 & 3 \\ \hline y & 0.35 & 0.80 & 1.70 & 1.85 & 3.51 & 1.02 \end{array}

Решение. Имеем

s_0=6, s_1=0.5 + 1 + 1.5 + 2 + 2.5 + 3=10.5,
s_2=0.5^2 + 1^2 + 1.5^2 + 2^2 + 2.5^2 + 3^2=22.75,
t_0=0.35 + 0.80 + 1.70 + 1.85 + 3.51 + 1.02=9.23,
t_1 =0.5\cdot 0.35 + 1 \cdot 0.80 + 1.5 \cdot 1.70 + 2 \cdot 1.85 +
+2.5 \cdot 3.51 + 3 \cdot 1.02=19.06 .

Решаем систему нормальных уравнений

\left( \begin{array}{cc} 6 & 10.5 \\ 10.5 & 22.75 \end{array} \right) \left( \begin{array}{c} a_0 \\ a_1 \end{array} \right)= \left( \begin{array}{c} 9.23 \\ 19.06 \end{array} \right),

получаем уравнение прямой в виде

y= 0.375 + 0.664\, x\ .

Псевдорешение системы линейных уравнений

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

\left\{\begin{array}{ccc} a_{11}x_1 +a_{12}x_2+\ldots+a_{1n}x_n &=&b_1\\ a_{21}x_1 +a_{22}x_2+\ldots+a_{2n}x_n &=&b_2\\ \ldots& & \ldots \\ a_{m1}x_1 +a_{m2}x_2+\ldots+a_{mn}x_n &=&b_m \end{array}\right. \iff AX={\mathcal B}

при числе уравнений m_{} большем числа неизвестных n_{}, то такая переопределенная система, как правило, несовместна. В этом случае задача может быть решена только путем выбора некоторого компромисса — все требования могут быть удовлетворены не полностью, а лишь до некоторой степени.

Псевдорешением системы AX={\mathcal B} называется столбец X\in \mathbb R^n, обеспечивающий минимум величины

\sum_{i=1}^m [a_{i1}x_1 +a_{i2}x_2+\ldots+a_{in}x_n-b_i]^2 \ .
§

Такому определению можно также соотнести вероятностную интерпретацию. Пусть для определения неизвестных величин x_{1},\dots,x_n проводятся m_{} экспериментов, описываемых линейными уравнениями:

a_{i1}x_1 +a_{i2}x_2+\ldots+a_{in}x_n =b_i \quad npu \quad i\in \{1,\dots,m\} \ .

При этом величины \{ a_{ij} \}, i \in \{1,\dots, m\}, j\in \{1,\dots, n\} — известные постоянные, не подверженные сопутствующим экспериментам (наблюдениям) погрешностям, а вот величины \{b_{i}\}_{i=1}^m этим погрешностям подвержены. Формально каждое из равенств следует рассматривать как приближенное. Понятно, что при таких обстоятельствах не имеет смысла гоняться за точным решением системы AX={\mathcal B} (его может и не существовать вовсе!). Искать следует приближенное решение, оптимальное в некотором смысле. Оказывается, что именно выбор критерия в виде минимума квадратов разностей левых и правых частей уравнения обеспечивает то, что псевдорешение дает максимально правдоподобные значения неизвестных величин x_1,\dots,x_{n}.

T

Теорема. Существует псевдорешение системы

AX={\mathcal B}

и оно является решением системы

\left[A^{\top}A \right]X=A^{\top} {\mathcal B} \ .

Это решение будет единственным тогда и только тогда, когда \operatorname{rank} A =n.

Система \left[A^{\top}A \right]X=A^{\top} {\mathcal B} называется нормальной системой по отношению к системе AX={\mathcal B}. Формально она получается домножением системы AX={\mathcal B} слева на матрицу A^{\top}. Заметим также, что если m=n_{} и \det A \ne 0, то псевдорешение системы совпадает с решением в традиционном смысле.

Доказательство ЗДЕСЬ.

Если нормальная система имеет бесконечное количество решений, то обычно в качестве псевдорешения берут какое-то одно из них — как правило то, у которого минимальна сумма квадратов компонент («длина»).

П

Пример. Найти псевдорешение системы

x_1+x_2 = 2, \ x_1-x_2 = 0,\ 2\, x_1+x_2 = 2 \ .

Решение. Имеем:

A=\left( \begin{array}{rr} 1 & 1 \\ 1 & -1 \\ 2 & 1 \end{array} \right),\ \operatorname{rank} A =2,\ {\mathcal B} = \left( \begin{array}{r} 2 \\ 0 \\ 2 \end{array} \right), \ A^{\top}A= \left( \begin{array}{rr} 6 & 2 \\ 2 & 3 \end{array} \right), \ A^{\top} {\mathcal B} = \left( \begin{array}{r} 6 \\ 4 \end{array} \right).

Ответ. x_1=5/7, x_2 = 6/7.

?

Показать, что матрица A^{\top}A всегда симметрична.

?

На дубовой колоде лежит мелкая монетка. К колоде по очереди подходят четыре рыцаря и каждый наносит удар мечом, стараясь попасть по монетке. Все промахиваются. Расстроенные, рыцари уходят в харчевню пропивать злосчастную монетку. Укажите максимально правдоподобное ее расположение, имея перед глазами зарубки:

\begin{array}{rrcr} 3\, x &- 2\, y&=& 6,\\ x &-3\,y&=&-3,\\ 11\,x& + 14\,y&=& 154, \\ 4\,x&+y&=&48. \end{array}







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

Псевдообратная матрица

!

Эта матрица определяется не только для квадратной матрицы.

Пусть сначала матрица A_{} порядка m\times n_{}вещественная и m \ge n_{} (число строк не меньше числа столбцов). Если \operatorname{rank} (A) = n (столбцы матрицы линейно независимы), то псевдообратная к матрице A_{} определяется как матрица

A^{+}=(A^{\top}A)^{-1} A^{\top} \ .

Эта матрица имеет порядок n \times m_{}. Матрица (A^{\top}A)^{-1} существует ввиду того факта, что при условии \operatorname{rank} (A) = n будет выполнено \det (A^{\top} A) > 0 (см. упражнение в пункте ТЕОРЕМА БИНЕ-КОШИ или же пункт СВОЙСТВА ОПРЕДЕЛИТЕЛЯ ГРАМА ). Очевидно, что A^{+} \cdot A = E_{n}, т.е. псевдообратная матрица является левой обратной для матрицы A_{}. В случае m=n_{} псевдообратная матрица совпадает с обратной матрицей: A^{+}=A^{-1}.

П

Пример. Найти псевдообратную матрицу к матрице

A= \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{array} \right) \ .

Решение.

A^{\top}= \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 1 \end{array} \right) \quad \Rightarrow \quad A^{\top} \cdot A = \left( \begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \quad \Rightarrow \quad (A^{\top} \cdot A)^{-1} = \left( \begin{array}{cc} 2/3 & -1/3 \\ -1/3 & 2/3 \end{array} \right) \quad \Rightarrow
\Rightarrow \quad A^{+} = (A^{\top} \cdot A)^{-1} A^{\top} = \left( \begin{array}{rrr} 2/3 & -1/3 & 1/3 \\ -1/3 & 2/3 & 1/3 \end{array} \right) \ .

При этом

A^{+} \cdot A = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right),\quad A \cdot A^{+} = \left( \begin{array}{rrr} 2/3 & -1/3 & 1/3 \\ -1/3 & 2/3 & 1/3 \\ 1/3 & 1/3 & 2/3 \end{array} \right) \ ,

т.е. матрица A^{+} не будет правой обратной для матрицы A_{}.

?

Вычислить псевдообратную матрицу для

\mathbf{a)}\ \left( \begin{array}{cc} 1 & 0 \\ 1 & 1 \\ 1 & 1 \end{array} \right) \quad ; \quad \mathbf{b)}\ \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) \ .

Концепция псевдообратной матрицы естественным образом возникает из понятия псевдорешения системы линейных уравнений. Если A^{+} существует, то псевдорешение (как правило, переопределенной и несовместной!) системы уравнений AX=\mathcal B_{} находится по формуле X= A^{+} \mathcal B при любом столбце \mathcal B_{}. Верно и обратное: если E_{[1]}, E_{[2]},\dots, E_{[m]} – столбцы единичной матрицы E_m:

E_{[1]}=\left( \begin{array}{c} 1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right),\ E_{[2]}=\left( \begin{array}{c} 0 \\ 1 \\ 0 \\ \vdots \\ 0 \end{array} \right),\dots, E_{[m]}=\left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \vdots \\ 1 \end{array} \right),

а псевдорешение системы уравнений AX=E_{[j]} обозначить X_{j} (оно существует и единственно при условии \operatorname{rank} (A) = n), то

A^{+}=\left[X_1,X_2,\dots,X_m \right] \ .
Т

Теорема. Пусть A_{} вещественная матрица порядка m\times n_{}, m \ge n_{} и \operatorname{rank} (A) = n. Тогда псевдообратная матрица A^{+} является решением задачи минимизации

\min ||AX-E_m||^2

где минимум разыскивается по всем вещественным матрицам X_{} порядка n\times m_{}, а || \cdot || означает евклидову норму матрицы :

||[h_{jk}]_{j,k}||^2=\sum_{j,k} h_{jk}^2 \ .

При сделанных предположениях решение задачи единственно.

§

Образно говоря, если уж невозможно найти обратную матрицу для матрицы A_{m\times n}^{}, давайте найдем хотя бы такую матрицу X_{n\times m}, чтобы отклонение произведения A\cdot X от единичной матрицы E_m, вычисленное как квадрат евклидова расстояния между матрицами A\cdot X и E_m, было бы минимальным.

С учетом этого результата понятно как распространить понятие псевдообратной матрицы на случай вещественной матрицы A_{m\times n}^{}, у которой число строк меньше числа столбцов: m < n_{}. Будем искать эту матрицу как решение задачи минимизации

\min ||YA-E_n||^2

где минимум разыскивается по всем вещественным матрицам Y_{} порядка n\times m_{}. Пусть \operatorname{rank} (A) = m, т.е. строки матрицы линейно независимы. Тогда псевдообратная к матрице A_{} определяется как матрица

A^{+}= A^{\top} (A\cdot A^{\top})^{-1} \ .

Очевидно, что в этом случае A\cdot A^{+}=E_m.

Задачи

Источники

[1]. Беклемишев Д.В. Дополнительные главы линейной алгебры. М.Наука.1983, с.187-234

1) Residual (англ.)

2018/12/10 16:17 редактировал au