Найти псевдорешение системы уравнений методом наименьших квадратов

VMath

Инструменты сайта

Основное

Навигация

Информация

Действия

Содержание

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

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

Пусть из физических соображений можно считать (предполагать), что величины $ 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_ $ и $ x_j $ пара точек $ (x_,y_i) $ и $ (x_,y_j) $ определяет прямую. Но другая пара точек определяет другую прямую, и у нас нет оснований выбрать какую-нибудь одну из всех прямых.

Часто в задаче удаленность точки от прямой измеряют не расстоянием, а разностью ординат $ k\,x_i+b-y_i $, и выбирают прямую так, чтобы сумма квадратов всех таких разностей была минимальна. Коэффициенты $ k_0 $ и $ b_ <0>$ уравнения этой прямой дают некоторое решение стоящей перед нами задачи, которое отнюдь не является решением системы линейных уравнений $$ k\,x_1+b=y_1,\dots, k\,x_+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_x^ $$ Величина $ \varepsilon_i=f(x_i)-y_i $ называется $ i_<> $-й невязкой 1) . Уравнения $$ \left\<\begin a_0+a_1x_1+\dots+a_x_1^&=&y_1, \\ a_0+a_1x_2+\dots+a_x_2^&=&y_2, \\ \dots & & \dots \\ a_0+a_1x_m+\dots+a_x_m^&=&y_m \end \right. $$ называются условными. Матрица этой системы — матрица Вандермонда (она не обязательно квадратная).

Предположим что данные интерполяционной таблицы $$ \begin x & x_1 & x_2 & \dots & x_m \\ \hline y & y_1 & y_2 &\dots & y_m \end $$ не являются достоверными: величины $ x_<> $ нам известны практически без искажений (т.е. на входе процесса мы имеем абсолютно достоверные данные), а вот измерения величины $ y_<> $ подвержены случайным (несистематическим) погрешностям.

Задача. Построить полином $ f_<>(x) $ такой, чтобы величина $$ \sum_^m [f(x_j)-y_j]^2 $$ стала минимальной. Решение задачи в такой постановке известно как метод наименьшик квадратов 2) .

В случае $ \deg f_<> =m-1 $ мы возвращаемся к задаче интерполяции в ее классической постановке. Практический интерес, однако, представляет случай $ \deg f_<> n $: $$S=\left(\begin 1 &1&1&\ldots&1\\ x_1 &x_2&x_3&\ldots&x_\\ \vdots&& & &\vdots\\ x_1^ &x_2^&x_3^&\ldots&x_m^ \end\right) \cdot \left(\begin 1 &x_1&x_1^2&\ldots&x_1^\\ 1 &x_2&x_2^2&\ldots&x_2^\\ \ldots&& & &\ldots\\ 1 &x_m&x_m^2&\ldots&x_m^ \end\right)$$ $$\det S = \sum_<1\le j_1 0 $. По теореме Крамера система нормальных уравнений имеет единственное решение.

Осталось недоказанным утверждение, что полученное решение доставляет именно минимум сумме квадратов невязок. Этот факт следует из доказательства более общего утверждения — о псевдорешении системы линейных уравнений. Этот результат приводится ☟ НИЖЕ. ♦

Собственно минимальное значение величины cуммы квадратов невязок, а точнее усреднение по количеству узлов $$ \sigma=\frac<1>\sum_^m (f(x_j) -y_j)^2 $$ называется среднеквадратичным отклонением.

Показать, что линейный полином $ y=a_<0>+a_1x $, построенный по методу наименьших квадратов, определяет на плоскости $ (x_<>,y) $ прямую, проходящую через центроид

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

$$ \begin 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 $$

Решение. Имеем $$ 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 6 & 10.5 \\ 10.5 & 22.75 \end \right) \left( \begin a_0 \\ a_1 \end \right)= \left( \begin 9.23 \\ 19.06 \end \right), $$ получаем уравнение прямой в виде $$ y= 0.375 + 0.665\, x\ .$$

Вычислим и полиномы более высоких степеней. $$ f_2(x)=-1.568+3.579\, x-0.833\,x^2 \ , $$ $$ f_3(x)=2.217-5.942\,x+5.475\,x^2-1.201\, x^3 \ , $$ $$ f_4(x)= -4.473+17.101\,x-19.320\,x^2+9.205\, x^3-1.487\,x^4 \ , $$ $$ f_5(x)= 16.390-71.235\,x+111.233\,x^2-77.620\,x^3+25.067\,x^4-3.0347\, x^5 . $$

Среднеквадратичные отклонения: $$ \begin \deg & 1 & 2 & 3 & 4 & 5 \\ \hline \sigma & 0.717 & 0.448 & 0.204 &0.086 & 0 \end $$ ♦

Возникает естественный вопрос: полином какой степени следует разыскивать в МНК? При увеличении степени точность приближения, очевидно, увеличивается. Вместе с тем, увеличивается сложность решения системы нормальных уравнений и даже при небольших степенях $ n $ (меньших $ 10 $) мы столкнемся с проблемой чувствительности решения к точности представления входных данных.

Влияние систематических ошибок

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

$$ \begin x & 0.5 & 1 & 1.5 & 2 & 2.5 & 3 \\ \hline y & 0.35 & 0.80 & 1.70 & 1.85 & 2.51 & 2.02 \end $$ имеет вид (охра) $$ y=0.175+0.779\, x \, . $$ Теперь заменим значение $ y_5 $ на $ 0.2 $. Уравнение прямой принимает вид: $$ y=0.483+0.383\, x \, . $$ График (зеленый) существенно изменился. Почему это произошло? — Дело в том, что эффективность метода наименьших квадратов зависит от нескольких предположений относительно входных данных: в нашем случае — значений $ y $. Предполагается, что эти величины являлись результатами экспериментов, измерений, и, если они подвержены погрешностям, то эти погрешности носят характер несистематических флуктуаций вокруг истинных значений. Иными словами, изначально предполагается, что в действительности точки плоскости должны лежать на некоторой прямой. И только несовершенство наших методов измерений (наблюдений) демонстрирует смещение их с этой прямой. Ответ для исходной таблицы визуально подтвержает это предположение: экспериментальные точки концентрируются вокруг полученной прямой и величины невязок (отклонений по $ y $-координате) имеют «паритет» по знакам: примерно половина точек лежит выше прямой, а половина — ниже.

После замены значения $ y_5 $ на новое, значительно отличающееся от исходного, существенно меняется величина $ 5 $-й невязки $ \varepsilon_5= ax_5+b-y_5 $. А поскольку в минимизируемую функцию эта невязка входи еще и в квадрате, то понятно, что изначальная прямая просто не в состоянии правильно приблизить новую точку.

Эта проблема становится актуальной в тех случаях, когда в «истинно случайный» процесс привносятся намеренные коррективы. Данные начинают подвергаться существенным искажениям, возможно, даже имеющим «злой» умысел 3) .

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

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

Рассмотрим теперь обобщение задачи предыдущего пункта. В практических задачах часто бывает нужно найти решение, удовлетворяющее большому числу возможно противоречивых требований. Если такая задача сводится к системе линейных уравнений $$ \left\<\begin 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_x_1 +a_x_2+\ldots+a_x_n &=&b_m \end\right. \iff AX= <\mathcal B>$$ при числе уравнений $ m_<> $ большем числа неизвестных $ n_<> $, то такая переопределенная система, как правило, несовместна. В этом случае задача может быть решена только путем выбора некоторого компромисса — все требования могут быть удовлетворены не полностью, а лишь до некоторой степени.

Псевдорешением системы $ AX= <\mathcal B>$ называется столбец $ X\in \mathbb R^n $, обеспечивающий минимум величины $$ \sum_^m [a_x_1 +a_x_2+\ldots+a_x_n-b_i]^2 \ . $$

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

$$ AX= <\mathcal B>$$ и оно является решением системы $$ \left[A^<\top>A \right]X=A^ <\top> <\mathcal B>\ . $$ Это решение будет единственным тогда и только тогда, когда $ \operatorname 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 1 & 1 \\ 1 & -1 \\ 2 & 1 \end \right),\ \operatorname A =2,\ <\mathcal B>= \left( \begin 2 \\ 0 \\ 2 \end \right), \ A^<\top>A= \left( \begin 6 & 2 \\ 2 & 3 \end \right), \ A^ <\top> <\mathcal B>= \left( \begin 6 \\ 4 \end \right). $$

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

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

На дубовой колоде лежит мелкая монетка. К колоде

по очереди подходят четыре рыцаря и каждый наносит удар мечом, стараясь попасть по монетке. Все промахиваются. Расстроенные, рыцари уходят в харчевню пропивать злосчастную монетку. Укажите максимально правдоподобное ее расположение, имея перед глазами зарубки: $$ \begin 3\, x &- 2\, y&=& 6,\\ x &-3\,y&=&-3,\\ 11\,x& + 14\,y&=& 154, \\ 4\,x&+y&=&48. \end $$

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

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

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

Пример. Найти псевдообратную матрицу к матрице $$ A= \left( \begin 1 & 0 \\ 0 & 1 \\ 1 & 1 \end \right) \ . $$

Решение. $$ A^<\top>= \left( \begin 1 & 0 & 1 \\ 0 & 1 & 1 \end \right) \ \Rightarrow \ A^ <\top>\cdot A = \left( \begin 2 & 1 \\ 1 & 2 \end \right) \ \Rightarrow $$ $$ \Rightarrow \ (A^ <\top>\cdot A)^ <-1>= \left( \begin 2/3 & -1/3 \\ -1/3 & 2/3 \end \right) \ \Rightarrow $$ $$ \Rightarrow \quad A^ <+>= (A^ <\top>\cdot A)^ <-1>A^ <\top>= \left( \begin 2/3 & -1/3 & 1/3 \\ -1/3 & 2/3 & 1/3 \end \right) \ . $$ При этом $$ A^ <+>\cdot A = \left( \begin 1 & 0 \\ 0 & 1 \end \right),\quad A \cdot A^ <+>= \left( \begin 2/3 & -1/3 & 1/3 \\ -1/3 & 2/3 & 1/3 \\ 1/3 & 1/3 & 2/3 \end \right) \ , $$ т.е. матрица $ A^ <+>$ не будет правой обратной для матрицы $ A_<> $. ♦

Вычислить псевдообратную матрицу для $$ \mathbf\ \left( \begin 1 & 0 \\ 1 & 1 \\ 1 & 1 \end \right) \quad ; \quad \mathbf\ \left( \begin x_1 \\ x_2 \\ x_3 \end \right) \ . $$

Концепция псевдообратной матрицы естественным образом возникает из понятия псевдорешения системы линейных уравнений . Если $ A^ <+>$ существует, то псевдорешение (как правило, переопределенной и несовместной!) системы уравнений $ AX=\mathcal B_<> $ находится по формуле $ X= A^ <+>\mathcal B $ при любом столбце $ \mathcal B_<> $. Верно и обратное: если $ E_<[1]>, E_<[2]>,\dots, E_ <[m]>$ – столбцы единичной матрицы $ E_m $: $$ E_<[1]>=\left( \begin 1 \\ 0 \\ 0 \\ \vdots \\ 0 \end \right),\ E_<[2]>=\left( \begin 0 \\ 1 \\ 0 \\ \vdots \\ 0 \end \right),\dots, E_<[m]>=\left( \begin 0 \\ 0 \\ 0 \\ \vdots \\ 1 \end \right),\ $$ а псевдорешение системы уравнений $ AX=E_ <[j]>$ обозначить $ X_ $ (оно существует и единственно при условии $ \operatorname (A) = n $), то $$ A^<+>=\left[X_1,X_2,\dots,X_m \right] \ . $$

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

$$ \min_> \|AX-E_m\|^2 $$ где минимум разыскивается по всем вещественным матрицам $ X_<> $ порядка $ n\times m_<> $, а $ || \cdot || $ означает евклидову норму матрицы (норму Фробениуса) : $$ \|[h_]_\|^2=\sum_ h_^2 \ . $$ При сделанных предположениях решение задачи единственно.

С учетом этого результата понятно как распространить понятие псевдообратной матрицы на случай вещественной матрицы $ A_^<> $, у которой число строк меньше числа столбцов: $ m ☞ ЗДЕСЬ

Источники

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

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

Как показано выше, система линейных алгебраических уравнений с неизвестными

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

Для дальнейших рассуждений нам потребуется количественная характеристика для «измерения» и «сравнения» между собой матриц-столбцов. Для этого каждому столбцу поставим в соответствие неотрицательное действительное число.

Модулем (нормой) столбца называется неотрицательное число

Заметим, что величина характеризует погрешность решения системы. Если величина велика, то столбец — плохое приближенное решение. Если погрешность мала, то столбец — хорошее приближенное решение. Если же является решением системы (в привычном понимании), то .

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

Любая система имеет единственное псевдорешение

где матрица — псевдообратная для матрицы системы.

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

Алгоритм нахождения псевдорешения неоднородной системы

1. Найти псевдообратную матрицу одним из способов, рассмотренных ранее.

2. Найти псевдорешение .

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

Решение. 1. Найдем псевдообратную матрицу. При решении примера 5.6 матрица системы была приведена к простейшему виду

По матрицам и находим псевдообратную матрицу. Вычисляем произведения и разбиваем их на блоки

По формуле (4.21) вычисляем псевдообратную матрицу

2. Находим псевдорешение

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

Модуль этого столбца имеет выражение

Вместо минимизации модуля будем искать минимум его квадрата, чтобы избавиться от квадратного корня. Эти задачи, разумеется, эквивалентны, так как точки минимума неотрицательной функции и ее квадрата совпадают. Итак, будем искать минимум функции

зависящей от двух переменных.

Для нахождения стационарных точек приравниваем частные производные первого порядка нулю (применяем необходимое условие экстремума):

Упрощая уравнения, получаем систему которую решаем, например, по правилу Крамера:

Поскольку угловые миноры матрицы Гессе положительные: 0,» png;base64,iVBORw0KGgoAAAANSUhEUgAAAHIAAAAUBAMAAACjXnQJAAAALVBMVEVHcEwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACttl6nAAAADnRSTlMAgUEBwBDoZach0DFRkT91eccAAAG6SURBVDjLY2CgBLBuIFcnx0vs4hMYWAMYmLUN0MUnBcBY815g1cj8ysNtAnOVSA6aeIRjLozZ8W4CVp3v3nky8GQyrFNAFW9hCBOAejO1rwGrzkWaDAyMLxjugdQxHoYJcz5nYHwMYTIK6GH1KLMzkOD2YfArABkvVQwVZn/EwPgGwgwLYHmBXacSKChY08ABwipVChG2fcTA9BDCFGRgQvLoJiUQ0ATpTBFMB1pa5Q31VJUomGaB62RNZWDwA3kUEquFgiAgDdKZylCnwMAu4QRzQ9VNVJ2MwADQS2XglHLAdLAeyMS+ArjWVhSdYUBvsDxkYJ3ogJm07JKB9LlkmE4VkE6ORwzsEJ0SQMz0LICBEaLzINy1XK+BOnlbGfSgccBc1gl2IzBsn4P5mcYGxsZ+AjCdiBAyyQS6Nu4NzE5maAhxA+PzJcMcYPSkuABBXipMJwIAA6BvA9dSqD9ZxUSh4h4MbA3MeQUMbO/A4CGmTuYWkSQG7hUSiSjRycAQ1e4SwAB0JLMxBGDqZLA+BMwmRkXguOYtRogfBMYfkwCS4xxIyZnmB+BMbreUAhJ0RhmQWRgwVzMAAOFba2RqKNaHAAAAAElFTkSuQmCC» style=»vertical-align: middle;» /> 0″ png;base64,iVBORw0KGgoAAAANSUhEUgAAALoAAAAXBAMAAABOokjyAAAALVBMVEVHcEwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACttl6nAAAADnRSTlMAgb0BQZwQIVnQ8DHgccy0xpgAAAKxSURBVEjH1ZTPaxNBFMdf9mfCRkgCIlKQdEOrQQlLVrBKCS0NehDCUgwYDWXxog1SgkcP0mqFNodStQh6EDGQFE8R9OhFpQhCyU0UFJI0mFqzf4Nvdjeb/dWL0kOHMJudmffZ9+M7D+AwjOljBwgPfv3eODg6FQ/F99ni/v+zgQi1H5395b9eBk4F+rI0XOEves/g6iWcShHXDp01DO91feH0rvi8DN9im0N4srDojPo+MrliAee827yS/6I/81rZl65pJyFYhxdTg5W7kcknjiPJFHLnFaoDE4rLmt+Eih5YLaX40jcw4kAbTlgxP5NK684zOdy7pk5W4SE0yfs7y0/2JWQIlo3kfBNPJ3ASTkGqaS4wLeAkD51vA60Ko1HdQ2ZsgA8n4DwJ9LZKdf3pWZXk9rU6cKe35L4zSA+2pkeAFUWjIow4Y+xkEpDu4DMKjC3xb7NkEBXQH6JYUGH1tCXqrufOID2we/0nIZtR8bMGPhfX6dwWACnOLUNtS1EyRgi9BmtTEDy6QZZvEnobQh+99D8SW7eXs/jApFM75CLg3xrwV854lEPjEaDn0DN6G51O93UDeKM70BjQuxDoO9RyfNGk75G0o1ctYBX3XcYrnKmSc1WTjmiq76W3yM9up9MXjLxHSS1+q+w6G3dmJryHdHYcch0zMwGT7szMhIvOF/TMXEgACpXbkiVZJokPKc6qpusY3kIP1qqmnYCZcWk3pwCHmenZ4LOLsi4BXe9HXok4tmsAq6q7M8FcI/R0qHcQpZLiLowERZUa3mDGlAwwO0QTYU0ffRDG3UUVY5+AF2PDsOdHHzs9SGlaB248Slp6YMastlH58Rk/LxsD7qhNF11YRg2fXbZJacXdq2mShfcz1vuwE9DnVmyG3NWsf4eW/7Xny7auQWmaAodr/AUkDbKlexLEIQAAAABJRU5ErkJggg==» style=»vertical-align: middle;» />, то найденная стационарная точка является точкой минимума.

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

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

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

Решение. Нетрудно заметить, что эта система несовместна (сложив второе и третье уравнения, получим уравнение , которое противоречит первому уравнению системы). Составим матрицу системы

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

2. По формуле (5.28) получаем псевдорешение

Покажем, что это решение минимизирует погрешность . Для этого найдем минимум квадрата этой функции (см. пример 5.8)

2. Приравнивая к нулю частные производные по переменным , получаем после упрощений систему уравнений

Методом Гаусса находим общее решение этой системы (5.29)

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

Найдем точку минимума этой функции одной переменной:

Поскольку функция выпуклая, наименьшее по модулю решение (5.29) получается при . Это решение совпадает с найденным ранее псевдорешением .

МЕТОД НАИМЕНЬШИХ КВАДРАТОВ И ПРИБЛИЖЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ В ЭКОНОМИЧЕСКИХ ЗАДАЧАХ

Кирилл Сукин 5 лет назад Просмотров:

1 Глава 3 МЕТОД НАИМЕНЬШИХ КВАДРАТОВ И ПРИБЛИЖЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ В ЭКОНОМИЧЕСКИХ ЗАДАЧАХ 3 Метод наименьших квадратов Метод наименьших квадратов (часто называемый МНК) обычно упоминается в двух контекстах Во-первых, широко известно его применение в регрессионном анализе как метода построения моделей на основе зашумленных экспериментальных данных При этом помимо собственно построения модели обычно осуществляется оценка погрешности, с которой были вычислены её параметры, иногда решаются и некоторые другие задачи Во-вторых, МНК часто применяется просто как метод аппроксимации, без какой-либо привязки к статистике При решении различных экономических задач мы часто сталкиваемся с необходимостью обработки массивов экспериментальных данных В частности, с необходимостью решения системы линейных уравнений (СЛУ) AX = B, () где A = (a j ) матрица размера m, X вектор-столбец из неизвестных, те X = ( ), а B вектор-столбец из m элементов, те B = (b b m ) Рассматриваемая СЛУ, как правило, несовместна, если в ней число уравнений превышает число неизвестных Поэтому говорить о ее точном решении невозможно Невязкой системы AX = B по -ой проекции относительно вектора X называется разность δ = ( a+ K+ a) b ( =, m) Определение Вектором невязки системы () относительно вектора называется вектор δ = ( AX B) 35

2 Ошибкой S(X) системы относительно вектора X называется квадрат длины (или нормы) вектора невязки, те сумма квадратов невязок по всем проекциям, те S( X) = δ + K + δm Отметим, что в евклидовом m-мерном точечном пространстве ошибка S() представляет собой квадрат расстояния между точками AX и B, где B = (b b m ) и ( AX ) = ( a + K+ a, K, a + K + a m m ) Если X = (, K, ) точное решение системы (), то векторстолбец B есть линейная комбинация столбцов A,, A матрицы A, где = (,, m),, = (,, ) A a K a K A a K a m Если X точное решение системы (), то вектор-столбец B является линейной комбинацией столбцов A,, A матрицы A, где = (,, m),, = (,, ) A a K a K A a K a m В качестве приближенного решения системы () выбирается такая точка X, для которой ошибка S(X) имеет наименьшее значение Наилучшее приближенное решение X* системы () совпадает с точным решением системы: 36 ( A, A) + K + ( A, A) = ( A, B), KK KK KK ( A, A) + K + ( A, A) = ( A, B) В матричной форме эта система записывается в виде A AX ( ) = A B () Определение СЛУ () (или ( )) называется системой нормальных уравнений для системы () Теорема Если столбцы матрицы A линейно независимы, то система нормальных уравнений () имеет единственное решение X*, и это решение дает наименьшее значение функции S( X) = AX B суммы квадратов ошибок по проекциям СЛУ () Доказательство этой теоремы можно найти в [0]

3 δ кр Приближенным решением системы () называется произвольный -мерный вектор X, если при этом ошибка S(X) не превосходит заранее заданного критического значения δ кр Псевдорешением системы () называется вектор X, в котором ошибка S(X) системы принимает свое наименьшее значение Если ошибка S(X) при этом не превосходит заданной точности, то псевдорешение является наилучшим приближенным решением системы () Теорема гарантирует существование и единственность псевдорешения системы () в достаточно общем случае Тот факт, что ошибка S(X) представляет собой сумму квадратов ошибок по всем проекциям объясняет название метода наименьших квадратов 4B3 Применение метода наименьших квадратов На практике метод наименьших квадратов (МНК) часто используется при анализе статистических данных для выявления функциональной зависимости между несколькими наблюдаемыми величинами Ограничимся случаем двух величин и y, между которыми предполагается наличие функциональной связи y = P(), где P ( ) = p + p+ p + K+ p многочлен степени k 0 При этом мы считаем, что число наблюдений, те пар (, y)( =,, K, ) экспериментальных данных не менее, чем число k + неизвестных коэффициентов многочлена В этом случае при решении СЛУ p + p + + p = y K K K K k p0 + p + K + pk = y k 0 K k, используется метод наименьших квадратов Матрица коэффициентов системы при этом имеет вид k K k K A = K K K K k K k k 37

4 Пример Пусть произведено взвешиваний одного и того же драгоценного камня Найти вес камня, используя метод наименьших квадратов 38

5 Решение Исходная система имеет уравнений с одним неизвестным, где вес камня: = b, K, = b или в матричном виде b b X = K K b При система, как правило, несовместна и A = (,, K,) Система принимает вид = b + K + b Отсюда получаем b + K+ b = псевдорешение исходной системы Если между величинами существует линейная зависимость y = α + β и проведено измерений (, y ), где =,, то МНК позволяет получить приближенную линейную зависимость y = α * + β * Пусть известно, что y = α + β и проведено измерений (, y ), где =,, K, Тогда для определения параметров α и β имеем СЛУ: α + β = y, K K K, α + β = y или в матричном виде y α K K = β K (3) y 39

6 Система нормальных уравнений для (3) принимает вид α + β = y α + β = y Решая систему (4), найдем y y β* = Если ввести обозначения =, y = y средние значения наблюдаемых величин, то первое уравнение примет вид α + β = y С учетом обозначений средних величин можно теперь записать решение системы (4) в виде формул * y y β =, ( ) α* = y β * Таким образом, в случае линейной зависимости МНК дает линейное уравнение y = α * + β *, где коэффициенты α* и β* находятся по формулам (5) Отметим, что в случае квадратичной (параболической) зависимости y= α + β+ γ система нормальных уравнений () прини- мает вид α + β + γ = y, 3 α + β + γ = y, 3 4 α + β + γ = y (4) (5) (6) 40

7 5B33 Приближенное решение систем линейных уравнений и обусловленность матриц Когда мы занимаемся задачей сглаживания экспериментальных данных, мы имеем дело, как правило, с несовместной СЛУ, в которой число неизвестных меньше числа уравнений Для решения такой задачи часто применятся МНК В результате мы получаем наилучшее приближенное решение исходной несовместной СЛУ Следует отметить, что, имея дело с экспериментальными данными, мы имеем дело с приближенными числами Поэтому в такой задаче важно уметь оценить, как сильно изменится решение исходной СЛУ, если каким-то образом изменится вектор-столбец свободных членов или исходная матрица коэффициентов СЛУ Для такой оценки вводится число обусловленности квадратной невырожденной матрицы Отметим, что при использовании МНК решается система нормальных уравнений (), матрица которой квадратная, симметрическая и невырожденная Определение Нормой матрицы A называется такое число A, для которого верны следующие свойства: ) A 0( A 0) > ; ) A + B A + B ; 3) λ A = λ A ; 4) AB A B ; 5) A A ; где A, B матрицы, вектор, λ число Примерами норм матрицы A = ( a j ) являются следующие числа: A где, =, j a j ; A = a, j j ; A 3 A = sup, 0 A введенные ранее нормы (длины) векторов и A Норма матрицы A называется согласованной с нормой вектора, A если A = sup 0 Всем перечисленным требованиям удовлетворяет, например, aj, A где a j элементы матрицы A, или aj, или sup,, j где,, j 0 A введенные ранее нормы (длины) векторов и A 4

8 Для согласованной нормы и единичной матрицы E, очевидно, имеем E= для любого вектора линейного пространства L, и, следовательно, E = Рассмотрим СЛУ A= b, где A квадратная матрица порядка, de A 0, b 0 В этом случае система имеет единственное ненулевое решение На практике при решении СЛУ вычисления производятся с округлениями, те неточно Погрешность вычислений при этом часто можно интерпретировать как погрешность правой части Рассмотрим, наряду с СЛУ A= b, систему A( + ) = b+ δ b, где δ b 0 погрешность правой части, а погрешность решения Из двух последних равенств вытекает Aδ = δ b и = A δb В силу свойств нормы матрицы имеем, что b A и A δb Перемножая эти неравенства, получим b δ A A δ b Разделив последнее неравенство на b 0, имеем A A δb Определение Если A невырожденная квадратная матрица, то число ca ( ) = A A называется числом обусловленности матрицы A С учетом определения числа обусловленности последнее неравенство можно переписать в виде δb ca ( ) (7) b Неравенство (7) дает оценку относительной погрешности вектора решения в зависимости от вектора-столбца свободных членов b Отметим, что число обусловленности так же важно, если матрица A имеет погрешность, а вектор b известен точно В самом деле, решим систему ( A + δ A)( + ) = b Тогда + δ= ( A+ δa) b равенства второе, получим 4 δ= (( A+ δa) A ) b и из b = A b Вычитая из первого Введем обозначение B = A+ δ A и используем тождество A ( A B) B = A AB A BB = EB A E = B A

9 Далее имеем = ( B A ) b = ( A ( A B) B ) b = = A δ AA ( + δa) ) b= A δa ( + δ) Отсюда A δa + δ, те A δ A Домножив и разделив правую часть последнего неравенства на A, + получим δ A ca ( ) + A, (8) те неравенство, аналогичное неравенству (7) Для вычисления евклидовой нормы матрицы и ее числа обусловленности используют две следующие теоремы Теорема Евклидова норма матрицы A, согласованная с евклидовой нормой векторов, равна A = ma λ ( AA ) (9) Для симметрической матрицы формула (9) упрощается A = ma λ ( A) (0) Здесь через λ ( B) обозначены собственные значения матрицы B Теорема Если матрица A симметричная, то число обусловленности матрицы можно вычислить по формуле ma λ ca ( ) = () m λ j j Если матрица A не является симметричной, то число обусловленности равно ca ( ) = caa ( ) () 43

10 6BПримеры Пусть имеются следующие экспериментальные данные о цене нефти и индексе нефтяных компаний: 7,30 7,08 8,30 8,80 9,0 8,50 y Считая, что между величинами и y существует линейная зависимость, применяя метод наименьших квадратов (МНК), найти эту зависимость Решение Найдем необходимые для расчетов средние величины, yy,, Запишем результаты этих вычислений в виде таблицы y y 7, ,36 98,5984 7, ,70 90, , ,00 334, , ,00 353,44 5 9, ,00 368,64 6 8, ,00 34,5 09, ,06 988,509 s 8, ,5 33,40 Используя формулы (5), получим β*,35, α* 37,9 По экспериментальным данным,0,5,0,5 3,0 y 0,50 0,30 0,5 0,8 0, построить с помощью МНК линейную зависимость y от Сравнить полученную зависимость с зависимостью y = Решение Найдем необходимые для расчетов средние величины, yy,, Запишем результаты этих вычислений в виде таблицы 44

11 y y,0 0,50 0,50,0,5 0,30 0,45,5 3,0 0,5 0,50 4,0 4,5 0,8 0,45 6,5 5 3,0 0, 0,36 9,0 0,35,6,5 s 0,7 0,45 4,5 Используя формулы (5), получим β* 0,76; α* 0,6 Следовательно, НМК дает зависимость y = 0,6 0,76 Сравним значения S = ( f( ) y ) для найденной линейной и предложенной экспоненциальной зависимости Запишем промежуточные результаты этих вычислений в виде таблицы y ( y ) (0,6 0,76 y ),0 0,50 0,5 0 0,0096,5 0,30 0,3536 0, , ,0 0,5 0,5 0 0,0004 4,5 0,8 0,768 0, , ,0 0, 0,5 0, , ,35 0, ,00736 Итак, результаты вычислений дают S лин = 0,00736, а S эксп = 0,00908 Ответ: S лин = 0,00736, S эксп = 0,00908; экспоненциальная зависимость лучше, чем линейная отражает экспериментальные данные 3 Имеются следующие данные о расходах на производство и доходах от реализации в условных единицах,5,0,5 3,0 y,4, 3,6 5, 8,0 Найти квадратичную зависимость МНК y= α + β+ γ с помощью 45

12 Решение Найдем необходимые для расчетов суммы 3 4. y, y, y Результаты вычислений приведем в таблице y 3 4 y y,0,4,0,0,0,4,4,5,,5 3,375 5,065 3,3 4,95 3,0 3,6 4,0 8,0 6,0 7, 4,4 4,5 5, 6,5 5,65 39,065,75 3, ,0 8,0 9,0 7,0 8,0 4,0 7,0 0 0,3,5 55,0 4,5 48,65 4,65 В результате получим систему уравнений вида 5α + 0β +, 5γ = 0, 3, 0α +, 5β + 55γ = 48, 65,, 5α + 55β + 4,5γ = 4, 65 Решив систему методом исключения неизвестных, получим приближенный ответ α,9; β,696; γ,9 Ответ: α,9; β,696; γ,9 4 Найдите число обусловленности симметричной матрицы A и оцените относительную погрешность решения системы A( + ) = b+ δb, если 3 4 0, A = 3, b =, δb = 0, , Решение Находим собственные значения матрицы A: λ 3 A λe = λ 3 = ( λ )( λ 4λ 5) = λ 46

13 Корни характеристического уравнения равны λ = ; λ = + 9 6,36; λ 3 = 9,36 Далее, ma λ = + 9, m λ =, ca= ( ) + 9, b = 6, δ b = 0,5 Отсюда по формуле (7) имеем 0,5 ( + 9) = 0,05( + 9) 0,05 6,36 0,59 6 Ответ: ca= ( ) + 9, 0,59 5 Найдите число обусловленности матрицы A и оцените относительную погрешность решения системы A( + ) = b+ δb, если 0 0 0, A = 3, b =, δb = 0, , Решение 3 Вычисляем AA= Находим собственные зна- 3 чения матрицы A A: λ 3 3 AA λe= 3 3 λ 3 = λ + 7λ 9λ+ 3= 0 3 λ Корни характеристического уравнения равны λ = ; λ = ,45; λ 3 = 3 6 0,55 Далее, ma λ = 3 + 6, m λ = 3 6, 3+ 6 caa ( ) = = , ca ( ) = = 3 + 3,46 47

14 4 Далее, b = 4, δ b = Отсюда по формуле (7) имеем, 0 что 3,46 0, 0,35 Ответ: ca= ( ) 3 + 3,46, 0,35 6 Оцените относительную погрешность решения системы ( A + δ A)( + ) = b, если 3 0, 0 0,05 A = 3 и δ A = 0 0, 0, 3 3 0,05 0, 0, 4 Решение В примере 5 уже было вычислено число обусловленности ca ( ) = + 9 = A Находим теперь собственные значения матрицы δ A: 0, λ 0 0,05 δa λe = 0 0, λ 0, = 0, 05 0, 0, 4 λ = (0, λ)( λ 0, 6λ + 0, 0675) = 0 Решая полученное характеристическое уравнение, получим собственные значения матрицы δ A: λ = 0, ; λ = 0, 45; λ = 0,5 Используя формулу (8), получим оценку относительной погрешности решения системы: 0,45 ( + 9) = 0,45 + ( + 9) Ответ: 0,45 + Ряд наблюдений значений некоторого показателя, упорядоченный в зависимости от значений другого показателя, называют динамическим или временным рядом Другими словами динамический ряд есть последовательность пар (, y ), в которой значения показателя упорядочены по возрастанию 48

15 Как правило, одним из показателей является время, которое часто обозначается условными моментами времени,,3, Отрезок динамического ряда, на котором значения показателя y изменяются монотонно, называется трендом Заметим, что систему нормальных уравнений можно упростить и уменьшить абсолютные значения величин, если перенести начало координат в середину динамического ряда Если до переноса начала координат равно,,3,, то после переноса получим: ) для четного числа членов равно K, 5, 3. 3,5, K, ) для нечетного числа членов равно K, 3. 0. 3, K В этом случае формулы (5) упрощаются и коэффициенты линейного тренда могут быть вычислены по более простым формулам: y α = y ; β = (3) В случае параболического тренда в этом случае тоже имеются достаточно простые формулы для коэффициентов α, β и γ, а именно: y α = ( y y) 4 ( ) β = y, ( y y) γ = (4) 4 ( ) 7 По данным ежемесячных объемов выпуска продукции фирмы за 8 месяцев: y рассчитать: а) коэффициенты α и β линейного тренда y = α + β и прогноз на месяц вперед; б) коэффициенты α, β и γ параболического тренда y= α + β+ γ и прогноз на месяц вперед Решение Для расчета коэффициентов линейного и параболического трендов воспользуемся выражениями, полученными из системы, 49

16 нормальных уравнений Перенесем начало координат (ί’) Необходимые вычисления занесем в таблицу Линейный тренд y ( ) y Итого Вычислим коэффициенты линейного тренда по формулам: Α = Σ y / = 4 634/8 = 3079,5; y 8066 β = = 48,0 (‘) 68 Таким образом, величина среднего уровня ряда при = 0 составляет 3079,5, среднемесячное уменьшение выпуска продукции составляет 48,0 Уравнение линейного тренда: y = 3079,5 48,0 ‘ Прогноз на 9-й месяц составит: y = 3079,5 48,0 9 = 647,6 Параболический тренд y (‘) y ‘ 3 (‘) 4 (‘) y (‘) Итого Вычислим коэффициенты параболического тренда по формулам (334): α = 3077,05; β = 48,0; γ = 0,05 Прогноз на 9-й месяц: y = , ,05 8 = 653,47 50

17 Ответ: Уравнение параболического тренда: y= 3077,05 48,0 ‘ + 0,05( ‘) Прогноз на 9-й месяц: y = 653, 47 8 В таблице приведены данные о доходности и рыночном индексе за 0 лет Предполагая зависимость между доходностью r и рыночным индексом r линейной, найти эту зависимость r = α + βr вид I Год r 0,8, 6,, 0,4 4,9 7,9,8 0,0 4,9 r 6,6 0,,5 3,9 0,4 8,4,6 0, 0,7 4,5 I Решение Система нормальных уравнений линейной зависимости имеет 0α + 30β =, 8, 30α + 63, 3β =, 696 Решая систему, получим α = 0,4, β = 0,79 Ответ: r = 0,4 + 0,79r I 7BУпражнения и тестовые вопросы 3 Что такое вектор невязки для несовместной СЛУ относительно вектора X? Найти вектор невязки 3+ 4y =, для СЛУ 5+ 3y = 3, относительно вектора X = (,3; ) + 5y = 3 Как связаны между собой вектор невязки и ошибка несовместной СЛУ относительно вектора X? Пусть вектор невязки несовместной СЛУ относительно вектора X равен δ = (0,3;0,;0,;0,;0,) Чему равна ошибка S(X) этой СЛУ относительно вектора X? 33 Что такое система нормальных уравнений для несовместной СЛУ? Написать систему нормальных уравнений несовместной 3+ 5y = 0, СЛУ 5+ 3y =, + 7y = I 5

18 34 Что называется числом обусловленности квадратной невырожденной матрицы? Указать формулы вычисления числа обусловленности матрицы в евклидовой норме в случае симметричной матрицы и в общем случае 35 Найти наилучшее приближенное решение несовместной СЛУ с помощью метода наименьших квадратов (МНК) + 3y = 8, 4+ 5y =, 5+ 3y = 9, 3+ 7y = 4 В ответе с точностью до 0,0 указать решение, вектор невязки δ, и длину вектора невязки 36 Пусть цена на товар (усле) и y объем продаж товара y Предполагая, что между и y существует линейная зависимость y = α + β, найти эту зависимость с помощью МНК Вычислить вектор невязки δ и ошибку системы S 37 По экспериментальным данным с помощью МНК постройте линейную зависимость y = α + β Сравните эту зависимость с зависимостью y= ,3, сопоставляя их ошибки S лин и S кв y В результате исследования зависимости между сроком эксплуатации автомобиля и расходами на его ремонт получены следующие данные: 5, лет S, тыс руб Найдите а) линейную зависимость S = α + β стоимости ремонта автомобиля от срока его эксплуатации;

19 б) предполагаемую величину затрат на 0-ый год его эксплуатации 53