Введение в вычислительную математику

Численное решение переопределенных СЛАУ. Метод наименьших квадратов

Разбить на страницы
Показывать лекцию целиком

3.1. Пример использования метода наименьших квадратов (МНК)

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

Пусть f — линейная (или близкая к линейной) функция аргумента x: f(x) = u1x + u0. В точках xk известны значения функции f(xk). Тогда u0, u1 — коэффициенты, которые необходимо подобрать так, чтобы выполнялись условия u1xk + u0 = fk, k = 0,1,2,3,4, fk = f(xk).

Получим систему пяти уравнений относительно двух неизвестных. Это — переопределенная система. Она не имеет классического решения, так как в общем случае не существует прямой, проходящей через все 5 точек (это возможно только тогда, когда какие - либо три уравнения полученной системы линейными преобразованиями сводятся к двум другим — система линейно зависима).

Рассмотрим общий случай. Пусть коэффициенты {u0, u1} необходимо определить по результатам n + 1 измерения. Введем функцию, равную сумме квадратов невязок rk = u1xk + u0 - fk

$$\Phi (u_1,u_0) = \sum\limits_{k = 0}^n {r_k^2} = \sum\limits_{k = 0}^n {(u_1 x_k + u_0 - f_k)^2} .$$

Примем за обобщенное решение переопределенной СЛАУ такие {u0, u1} для которых $$\Phi (u_{0}, u_{1})$$ принимает наименьшие значение. Для определения обобщенного решения из условия минимума суммы квадратов невязки получаем систему двух уравнений, имеющую классическое решение:

$$$ \frac{\partial \Phi }{\partial u_0} = 0, \frac{\partial \Phi }{\partial u_1} = 0. $$$

Выбор функции $$\Phi (u_{0}, u_{1})$$ имеет некоторый произвол. Например, возможно каждому измерению придать некоторый вес bk. От набора таких весовых множителей зависело бы решение системы. В этом случае функция $$\Phi$$ будет

$$\Phi (u_0,u_1) = \sum\limits_{k = 0}^n {b_k{(u_1 x_k - u_0 - f_k)}^2}.$$

Если в качестве невязки выбрать rk = | u1xk + u0 - fk |, то получим задачу линейного программирования на отыскании минимума функции

$$\Phi (u_1 u_2) = \sum\limits_{k = 0}^3 \left|{u_1 x_k + u_0 - f_k}\right|.$$

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

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

$$r_k = \sum\limits_{j = 0}^p {u_j\varphi_j (x_k)} - f(x_k), k = 1, \ldots, n,$$

где $$\varphi_j (x)$$ — некоторые функции, образующие базис, например, тригонометрические: $$\varphi_j (x) = \sin (jx).$$ Выражение $$\sum\limits_{j = 0}^p {u_j \varphi_j (x)}$$ называется обобщенным полиномом. В приведенном выше примере в качестве базисных функций были выбраны степенные функции $$\varphi_j (x) = x^j .$$ Обобщенный полином превратился в алгебраический.

В случае выбора произвольной системы базисных функций переопределенная СЛАУ и функционал $$\Phi (u_{0}, \dots , u_{p})$$ будут

$$\begin{gather*} u_0 \varphi_0 (x_0) + \ldots + u_p \varphi_p(x_0) = f_0, \\ \ldots \\ u_0 \varphi_0 (x_n) + \ldots + u_p\varphi_p (x_n) = f_n,\\ \Phi (u_0,\ldots,u_n) = \sum\limits_{i = 0}^n (\sum\limits_{j = 0}^p u_j \varphi_j(x_i) - f_i)^2 \end{gather*}$$

Отыщем обобщенное решение методом наименьших квадратов. Приравнивая все частные производные по компонентам обобщенного решения к нулю $$$ \frac{\partial \Phi }{\partial u_k} = 0 $$$ (условия минимума) и изменяя порядок суммирования, получаем СЛАУ

$$\sum\limits_{j = 0}^p {\left({\sum\limits_{i = 0}^n{\varphi_j (x_i)\varphi_k (x_i)}}\right)u_j = \sum\limits_{i = 0}^n {f_i\varphi_k (x_i)}}, k = 0, \ldots, p,$$

или

$$\begin{gather*} (\varphi_0, \varphi_0) u_0 + (\varphi_0, \varphi_1)u_1 + \ldots + (\varphi_0, \varphi_p)u_p = (\varphi_0, f), \\ (\varphi_1, \varphi_0) u_0 + (\varphi_1, \varphi_1)u_1 + \ldots + (\varphi_1, \varphi_p)u_p = (\varphi_1, f), \\ (\varphi_p, \varphi_0) u_0 + (\varphi_p, \varphi_1)u_1 + \ldots + (\varphi_p, \varphi_p)u_p = (\varphi_p, f), \end{gather*}$$

Система метода наименьших квадратов имеет вид $$\mathbf{Du} = \mathbf{f}$$ с матрицей $$\mathbf{D},$$ элементами которой являются скалярные произведения $$(\varphi_i, \varphi_j) = \sum\limits_{i = 0}^n \varphi_j (x) \varphi_k (x_i)$$ Это — матрица Грама. Ее свойства известны из курса линейной алгебры, эта матрица симметричная и положительно определенная. Таким образом, решение исследуемой СЛАУ существует и единственно. В правой части системы стоят проекции свободного члена исходной задачи на подпространство базисных функций $$(\varphi,f) = \sum\limits_{i = 0}^n {\varphi_j(x_i)f_i}.$$

Здесь учтено, что $$$ \frac{\partial \Phi }{\partial u_k} = 2 \sum\limits_{i = 0}^n {\varphi_k(x_i)\left({\sum\limits_{j = 0}^p {u_j\varphi_j (x_i) - f_i}} \right)} $,$$ или, в развернутом виде

$$\begin{gather*} \sum\limits_{i = 0}^n {\varphi_0 (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_0(x_i) + \ldots + u_p \varphi_p(x_i) - f(x_i)} \right) = 0}, \\ \sum\limits_{i = 0}^n {\varphi_1 (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_1 (x_i) + \ldots + u_p \varphi_p (x_i) - f(x_i)} \right) = 0}, \\ \ldots, \\ \sum\limits_{i = 0}^n {\varphi_{s} (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_1 (x_i) + \ldots + u_p \varphi_p (x_i) - f(x_i)} \right) = 0.} \end{gather*}$$

Часто выбирают $$\varphi_k (x) = x^k,$$ в этом случае система уравнений принимает следующую форму:

$$\sum\limits_{j = 0}^p (\sum\limits_{i = 0}^n x_i^{j + k})u_j = \sum\limits_{i = 0}^n f_i x_i^k, k = 1, \ldots, p.$$

Эта система может быть легко выписана в компонентах:

$$\begin{gather*} u_0 + \left({\sum\limits_{i = 0}^n {x_i}} \right)u_1 + \ldots + \left({\sum\limits_{i = 0}^n {x_i^p}}\right)u_p = \sum\limits_{i = 0}^n{f(x_i)}, \\ \left({\sum\limits_{i = 0}^n {x_i}} \right)u_0 + \left({\sum\limits_{i = 0}^n {x_i^2}} \right)u_1 + \ldots + \left({\sum\limits_{i = 0}^n {x_i^{p + 1}}} \right) = \sum\limits_{i = 0}^n {x_i f(x_i)}, \\ \ldots \\ \left({\sum\limits_{i = 0}^n {x^p_i}} \right)u_0 + \left({\sum\limits_{i = 0}^n {x^{p + 1}_i}}\right) + \ldots + \left({\sum\limits_{i = 0}^n{x^{2p + 1}_i}}\right)u_p = \sum\limits_{i = 0}^n {x^p_if(x_i)}. \end{gather*}$$

В случае использования ортонормированных систем базисных функций $$\varphi_j (x),$$ т.е., при выполнении условия $$(\varphi_k, \varphi_j) = \delta_{kj}$$ решение принимает простой вид $$$ u_0 = \frac{(\varphi_0, f)}{(\varphi_0,\varphi_0)} = {\left\| {\varphi_0} \right\|}^{- 2} (\varphi_0, f); u_1 = {\left\| {\varphi_1} \right\|}^{- 2} (\varphi_1, f); \ldots u_p = {\left\| {\varphi_p} \right\|}^{- 2} (\varphi_p, f) $.$$ Поскольку система функции ортонормирована, то $$\left\| {\varphi_j}\right\| = 1$$ и $$u_k = (\varphi_k, f), k = 1, \ldots, n.$$ Эти коэффициенты называются коэффициентами Фурье, а обобщенный многочлен с этими коэффициентами — обобщенным многочленом Фурье. В частности, в качестве базисных можно использовать ортогональную систему функций на $$[- \pi, \pi ]: \{\sin{kx},\cos{kx}\}$$ k = 1, ..., n. Такие представления называют отрезками тригонометрических рядов Фурье или конечными рядами Фурье .

Докажем теорему о методе наименьших квадратов, обобщающую изложенную информацию.

Запишем переопределенную СЛАУ

$$\begin{gather*} a_{11} u_1 + \ldots + a_{1p} u_p = f_1, \\ \ldots \\ a_{n1} u_1 + \ldots + a_{np} u_p = f_n, n > p. \end{gather*}$$

$$\mathbf{u} = {\{u_1, \ldots,u_p \}}^T \in L^p, \mathbf{f} = {\{f_1, \ldots,f_n \}}^T \in L^n,$$

где линейные нормированные пространства Lp и Ln имеют размерности p и n соответственно. Перепишем (3.2) в матричной форме:

$$\mathbf{Au} = \mathbf{f},$$

где

$$\mathbf{A} = \left( \begin{array}{ccc} {a_{11}} \ldots {a_{1p}} \\ \ldots \ldots \ldots \\ {a_{n1}} \ldots {a_{np}}} \end{array} \right).$$

Наряду с основным скалярным умножением в Ln

$${(\mathbf{x},\mathbf{y})}^n = \sum\limits_{k = 1}^n{x_k y_k}$$

введем скалярное умножение с весовой матрицей $$\mathbf{B}$$:

$${\left[{\mathbf{x},\mathbf{y}}\right]}^n = {(\mathbf{Bx},\mathbf{y})}^n, \mathbf{B} = {\mathbf{B}}^* > 0, \mathbf{x}, \mathbf{y} \in L^n$$

Оба этих умножения удовлетворяют аксиомам скалярного умножения элементов линейного пространства. Матрица $$\mathbf{B}$$ является весовой и определяет вклад невязки каждого слагаемого суммы (3.1). Система (3.2) не имеет классического решения. Определим обобщенное решение этой системы как элемент линейного пространства $$\mathbf{v},$$ придающий наименьшее значение квадратичной форме:

$$\Phi (\mathbf{u}) = {\left[\mathbf{Au} - \mathbf{f}, \mathbf{Au} - \mathbf{f}\right]}^n.$$

Теорема. Пусть столбцы матрицы $$\mathbf{A}$$ линейно независимы, т.е. ранг $$\mathbf{A}$$ равен p. Тогда существует единственный элемент p - мерного евклидова пространства $$\mathbf{v} \in L^P,$$ являющийся обобщенным решением системы (3.2), и решением СЛАУ

$${\mathbf{A}}^*{\mathbf{BAu}} = {\mathbf{A}}^*{\mathbf{Bf}},$$

состоящей из p скалярных уравнений относительно неизвестных $$\left\{{\mathbf{u}_k}\right\}_{k = 1}^{k = p},$$ доставляющий минимум квадратичной форме

$$\Phi (\mathbf{u}) = {\left[\mathbf{Au} - \mathbf{f}, \mathbf{Au} - \mathbf{f}\right]}^n.$$

Доказательство.

Покажем, что решение СЛАУ $${\mathbf{A}}^*{\mathbf{BAu}} = {\mathbf{A}}^*{\mathbf{Bf}}$$ существует и единственно. Введем обозначение $${\mathbf{q}}_k \in L^n$$ для вектора — k столбца матрицы системы $$\mathbf{A}$$:

$${\mathbf{q}}_k = {\{a_{1k}, \ldots, a_{nk} \}}^T; k = 1, \ldots, p.$$

Несложно показать, что матрица $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ системы (3.5) есть квадратная матрица p x p. Элемент dij этой матрицы, стоящий на пересечении i строки и j столбца, есть $$d_{ij} = {({\mathbf{q}}_i, {\mathbf{Bq}}_j)}^n = {({\mathbf{Bq}}_i,{\mathbf{q}}_j)}^n = {\left[{{\mathbf{q}}_i,{\mathbf{q}}_j}\right]}^n,$$ в силу коммутативности скалярного произведения dij = dij, что означает самосопряженность матрицы $$\mathbf{D}: \mathbf{D} = {\mathbf{D}}^*.$$

Покажем, что матрица $$\mathbf{D}$$ невырождена и положительно определена. Напомним, что $${(\mathbf{f},\mathbf{A\xi })}^n = {({\mathbf{A}}^*{\mathbf{f}},{\mathbf{\xi }})}^P; {\mathbf{f}}_1 \in L^n, \mathbf{\xi } \in L^P.$$ Это равенство проверяется непосредственно, если записать обе его части в развернутом виде. Невырожденность матрицы $$\mathbf{D}$$ следует из того, что ранг матрицы $$\mathbf{A}$$ равен p.

В таком случае $$0 < {\left[{\mathbf{Ax},\mathbf{x}}\right]}^n = {(\mathbf{BAx},\mathbf{Ax})}^n = {({\mathbf{A}}^*{\mathbf{BAx}},{\mathbf{x}})}^p = {(\mathbf{Dx},\mathbf{x})}^p.$$ Поскольку $$\mathbf{D}$$ невырождена и положительно определена, то (3.5) имеет единственное решение $$\mathbf{v} \in L^p.$$ Теперь покажем, что $$\mathbf{v}$$ — единственное обобщенное решение системы. Для любого вектора $$\mathbf{\Delta } \ne 0$$ выполнено $$\Phi (\mathbf{v} + \mathbf{\Delta }) > \Phi (\mathbf{v})$$:

$$\Phi (\mathbf{v} + \mathbf{\Delta }) = {\left[{\mathbf{A}(\mathbf{v} + \mathbf{\Delta }) - \mathbf{f},\mathbf{A}(\mathbf{v} + \mathbf{\Delta }) - \mathbf{f}}\right]}^n = \\ = {\left[{\mathbf{Av} - \mathbf{f},\mathbf{Av} - \mathbf{f}}\right]}^n - 2{\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n + {\left[{\mathbf{A\Delta },\mathbf{A\Delta }}\right]}^n = \\ = \Phi (\mathbf{v}) + {2\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n + {(\mathbf{D\Delta },\mathbf{\Delta })}^p = \Phi (\mathbf{v}) + {(\mathbf{D\Delta },\mathbf{\Delta })}^p > \Phi (\mathbf{v}),$$

что и требовалось доказать.

При доказательстве использовалось

$${\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n = {(\mathbf{B}(\mathbf{Av} - \mathbf{f}),\mathbf{A\Delta })}^n = {({\mathbf{A}}^*{\mathbf{BAv}} - {\mathbf{A}}^*{\mathbf{Bf}},\mathbf{\Delta })}^p = 0,$$

поскольку $${\mathbf{A}}^*{\mathbf{BAv}} = {\mathbf{A}}^*{\mathbf{Bf}}.$$

Так как матрица $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ — симметричная и положительно определенная, то для численного решения полученной СЛАУ можно воспользоваться итерационными методами.

Если система векторов $$\left\{{{\mathbf{q}}_k}\right\}_{k = 1}^p$$ оказывается ортонормированной, т.е. $$[{\mathbf{q}}_i,{\mathbf{q}}_j ] = \delta_{ij} i,j = 1, \ldots, p,$$ то матрица $$\mathbf{D}$$ оказывается единичной. Ее элементы и есть скалярные произведения $$[{\mathbf{q}}_i,{\mathbf{q}}_j ].$$ В этом случае решением системы будет

$$\mathbf{u} = \mathbf{A^*Bf}.$$

Следует отметить, что если базисные функции $$\varphi_i (x), i = 0, \ldots, p$$ не выбираются специальным образом, то при достаточно больших p $$(p \ge 5)$$ полученная СЛАУ оказывается плохо обусловленной. Строки матрицы $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ могут оказаться почти линейно зависимыми. Простейшим примером такого почти линейно зависимого базиса является система функций xi, i = 1, ..., p при больших p. В этом случае желательно использовать ортогональные функциональные базисы, однако такой выбор не всегда возможен и удобен.

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

$$f(x) = \sum\limits_{j = 0}^p {u_j x^j}.$$

Скалярные произведения на отрезке [0, 1], записанные в интегральной форме (т.е. при $$n \to \infty$$ ), будут иметь вид

$$$ (\varphi_i,\varphi_j) = \int\limits_0^1 x^i x^j dx = \int\limits_0^1 x^{i + j}dx = \frac{1}{i + j + 1}. $$$

В таком случае СЛАУ после применения МНК, т.е. минимизации функционала $$\Phi (u) = \int\limits_0^1 {(F(x) - f(x))}^2dx,$$ где F(x) — заданная функция, будет:

$$\begin{gather*} u_0 + \frac{1}{2}u_1 + \ldots + \frac{1}{p + 1}u_p = \int\limits_0^1{F(x)dx}, \\ \frac{1}{2}u_0 + \frac{1}{3}u_1 + \ldots + \frac{1}{p + 2}u_p = \int\limits_0^1 xF(x)dx, \\ \ldots \\ \frac{1}{p + 1}u_0 + \frac{1}{p + 2}u_1 + \ldots + \frac{1}{2p + 1}u_p = \int\limits_0^1 pF(x)dx, \end{gather*}$$

или

$${\mathbf{H}}_{p + 1}\mathbf{u} = \mathbf{F^{\prime}},$$

где

$$\begin{gather*} \mathbf{u} = {\{u_0, \ldots, u_p \}}^T,\quad \mathbf{F^{\prime}} = {\left\{{\int\limits_0^1 {F(x)dx}, \ldots, \int\limits_0^1 {x^p F(x)dx}} \right\}}^T, \\ {\mathbf{H}}_{p + 1} = {\left\{\frac{1}{i + j - 1}\right\}}^{p + 1}_{i,j} = 1 \end{gather*}$$

Матрица $$\mathbf{H}_{p + 1}$$ называется матрицей Гильберта. Это классический пример плохо обусловленной матрицы. Число обусловленности очень быстро растет с ростом p. Так при p = 1 $$\mu = \left\|\mathbf{H}\right\| \cdot \left\|{\mathbf{H}}^{- 1}\right\| \approx 20,$$ при $$p = 9 \mu \approx 10^{13}.$$ Если получим СЛАУ для дискретной системы точек, т.е. для $$$ (\varphi_j,\varphi_k) = \sum\limits_{i = 0}^n x^{k + j}_i, (\varphi_j,f) = \sum\limits_{i = 0}^n x^k_i f(x_i)x_i = \frac{i}{n} $,$$ то ее матрица будет асимптотически приближаться к матрице Гильберта $$\mathbf{H}_{p + 1}$$ при $$n \rightarrow \infty.$$

3.2. Понятие о методах решения плохо обусловленных СЛАУ

Улучшить качество численного решения СЛАУ метода наименьших квадратов возможно, если использовать различные преобразования матрицы $$\mathbf{A}.$$

Большинство прямых методов решения линейных систем основано либо на замене исходной системы $$\mathbf{Au} = \mathbf{f}$$ ( $$\mathbf{A}$$ — квадратная матрица) на эквивалентную $$\mathbf{CAu} = \mathbf{Cf},$$ либо на представлении матрицы $$\mathbf{A}$$ в виде произведения других матриц, таких, чтобы новая система либо решалась более просто, либо была лучше (по крайней мере, не хуже) обусловлена, чем исходная.

Подход, использующий спектральную эквивалентность матриц $$\mathbf{A}$$ и $${\mathbf{C}}^{- 1}$$ (в смысле границ спектра собственных значений), основан на умножении $$\mathbf{A}$$ на близкую в некотором смысле матрицу $$\mathbf{C}.$$ Последняя матрица выбирается таким образом, чтобы произведение было близким к единичной матрице (при этом $${\mathbf{C}} \ne {\mathbf{A}}^{- 1},$$ так как обращение плохо обусловленной матрицы приводит к накоплению вычислительных ошибок. Число обусловленности матрицы $$\mathbf{CA}$$ $$\mu = \lambda _{max}/\lambda _{min}$$ будет близко к единице. Метод энергетически эквивалентных операторов оказался эффективным при численном решении сеточных уравнений при разностной аппроксимации уравнений в частных производных эллиптического типа.

Идея предобусловливания СЛАУ состоит в том, чтобы вместо исходной системы $$\mathbf{Au} = \mathbf{f}$$ решать систему $$\mathbf{A^{\prime}u^{\prime}} = \mathbf{f^{\prime}},$$ где $$\mathbf{A^{\prime}} = \mathbf{C^{- 1}AC^{- 1}}, \mathbf{u^{\prime}} = \mathbf{Cu}, \mathbf{f^{\prime}} = \mathbf{C^{- 1}f},$$ матрица $$\mathbf{C}$$ выбирается так, чтобы она была симметричной положительной, хорошо обусловленной.

Определение. Матрица $$\mathbf{Q}$$ с вещественными элементами qij является ортогональной, если

$$\mathbf{Q}^* = \mathbf{Q}^{- 1}.$$

Пусть матрица $$\mathbf{C}$$ невырождена. Тогда она представима в виде

$$\mathbf{C} = \mathbf{QR},$$

где $$\mathbf{Q}$$ — ортогональная, а $$\mathbf{R}$$ — верхняя треугольная матрицы. В качестве матрицы $$\mathbf{Q}$$ часто используется симметричная ортогональная матрица $$\mathbf{H}$$:

$$\mathbf{H} = \mathbf{E} - 2{\mathbf{ww}}^T,$$

где $$\mathbf{w}$$ - произвольный вектор - столбец, такой, что $$({\mathbf{w}}^T{\mathbf{w}}) = 1.$$ Матрица $$(\mathbf{w,wT})$$ есть произведение вектора - столбца $$\mathbf{w}$$ на вектор - строку $$\mathbf{wT}$$ (преобразование Хаусхолдера или метод отражений ).

Заметим, что симметричность матрицы $$\mathbf{H}$$ (матрицы отражений) показывается непосредственной проверкой; ортогональность $$\mathbf{H}$$ можно показать следующим образом:

$${\mathbf{HH}}^T = (\mathbf{E} - 2{\mathbf{ww}}^T)(\mathbf{E} - 2{\mathbf{ww}}^T) = \\ = \mathbf{E} - 2{\mathbf{ww}}^T - 2{\mathbf{ww}}^T + 4{\mathbf{ww}}^T \cdot {\mathbf{ww}}^T = \\ = \mathbf{E} - 4{\mathbf{ww}}^T - 4\mathbf{w}({\mathbf{w}}^T{\mathbf{w}}){\mathbf{w}}^T = \mathbf{E}, \\ $$

так как $${\mathbf{w}}^T{\mathbf{w}} = 1.$$

Матрица СЛАУ представляется в виде $$\mathbf{A} = {\mathbf{H}}^T{\mathbf{R}},$$ после чего решается эквивалентная система уравнений $$\mathbf{Ru} = \mathbf{Hf}.$$

3.3. Задачи

  • Для функции $$f(x) = \sqrt {x}$$ на отрезке [0,1] построить многочлен F(x) = u0 + u1x среднеквадратичного приближения со скалярными произведениями:

    $$(\varphi_i,\varphi_j) = \int\limits_0^1 {x^ix^j}{dx},\quad f,\varphi_i) = \int\limits_0^1 f(x)x^idx.$$

    Решение. Введем базисные функции $$\varphi_0 (x) = 1, \varphi_1 (x) = x,$$ и вычислим скалярные произведения

    $$$ (\varphi_0,\varphi_0) = \int\limits_0^1 {1^2 dx} = 1, (\varphi_1,\varphi_1) = \int\limits_0^1 {x^2 dx =} \frac{1}{3},$$$

    $$$ (\varphi_0,\varphi_1) = \int\limits_0^1 {xdx}} = \frac{1}{2}, (f,\varphi_0) = \int\limits_0^1 {\sqrt{1} x dx} = \frac{2}{3},$$$

    $$$ (f,\varphi_1) = \int\limits_0^1 {\sqrt{x} x dx} = \frac{2}{5},$$$

    Для вычисления коэффициентов получим СЛАУ

    $$$ u_0 + \frac{1}{2}u_1 = \frac{2}{3},$$$

    $$$ \frac{1}{2}u_0 + \frac{1} 3u_1 = \frac{2}{5}, $$$

    откуда $$$ u_0 = \frac{4}{{15}}, c_1 = \frac{4}{5}, F(x) = \frac{4}{15} + \frac{4}{5}x.$$$

  • Получить СЛАУ
    c11x + c12y = f1
    c22x + c21y = f2,

    если она задана в форме метода наименьших квадратов: $$\mathbf{A^*BAu} = \mathbf{A^*Bf},$$ где

    $$\mathbf{A} = \left( \begin{array}{cc} {a_{11}} {a_{12}} \\ {a_{21}} {a_{22}} \\ {a_{31}} {a_{32}} \end{array} \right), \mathbf{B} = \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right), \mathbf{u} = {(u_1,u_2)}^T, \mathbf{f} = {(f_1,f_2,f_3)}^T .$$

    Решить эту систему для случая a11 = a32 = 1, a21 = a12 = 2, a31 = 0, a22 = 1, f1 = 1, f2 = 2, f3 = 1.

    Решение. Проведем необходимые вычисления:

    $$\begin{gather*} \left( \begin{array}{ccc} {a_{11}} {a_{21}} {a_{31}} \\ {a_{12}} {a_{22}} {a_{32}} \end{array} \right) \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right) \left( \begin{array}{cc} {a_{11}} {a_{12}} \\ {a_{21}} {a_{22}} \\ {a_{31}} {a_{32}} \end{array} \right) = \\ = \left( \begin{array}{ccc} {a_{11}} {a_{21}} {a_{31}} \\ {a_{12}} {a_{22}} {a_{32}} \end{array} \right) \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right) \left( \begin{array}{l} {f_1} \\ {f_2} \\ {f_3} \end{array} \right), \\ \left( \begin{array}{cc} {a_{11}^2 + a_{21}^2 + a_{31}^2} {a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}} \\ {a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}} {a_{12}^2 + a_{22}^2 + a_{32}^2} \end{array} \right) = \\ = \left( \begin{array}{l} {u_1} \\ {u_2} \end{array} \right) = \left( \begin{array}{l} {a_{11}f_1 + a_{21}f_2 + a_{31}f_3} \\ {a_{12}f_1 + a_{22}f_2 + a_{32}f_3} \\ \end{array} \right), \\ (a_{11}^2 + a_{21}^2 + a_{31}^2)u_1 + (a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32})u_2 = a_{11}f_1 + a_{21}f_2 + a_{31}f_3, \\ (a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}) u_1 + (a_{12}^2 + a_{22}^2 + a_{32}^2)u_2 = a_{12}f_1 + a_{22}f_2 + a_{32}f_3. \end{gather*}$$

    Решение в числах предлагается найти читателям.

  • С помощью метода наименьших квадратов найти коэффициенты полинома второй степени f(x) = u0 + u1x + u2x2, если таблица измерений задана $$\left. {\left\{{x_i, f_i}\right\}}\right|_{i = 0}^n .$$

    Решение. Переопределенная система уравнений:

    $$u_0 + u_1 x_k + u_2 x_k^2 = f_k, k = 0, 1, \ldots, n,$$

    определяет функционал

    $$\Phi (u_0, \ldots, u_n) = \sum\limits_{k = 0}^n{r_k^2} = \sum\limits_{k = 0}^n {(u_0 + u_1 x_k + u_2 x_k^2 - f_k)}^2.$$

    Условия минимума последнего

    $$$ \frac{\partial \Phi }{\partial u_k} = 0. $$$

    После проведения алгебраических преобразований получим для коэффициентов систему линейных уравнений

    $$\begin{gather*} u_0 + \left({\sum\limits_{k = 0}^n {x_k}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^2}}\right)u_2 = \sum\limits_{k = 0}^n {f_k}, \\ \left({\sum\limits_{k = 0}^n {x_k}} \right)u_0 + \left({\sum\limits_{k = 0}^n {x_k^2}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^3}} \right)u_2 = \sum\limits_{k = 0}^n {f_k} x_k, \\ \left({\sum\limits_{k = 0}^n {x_k^2}} \right)u_0 + \left({\sum\limits_{k = 0}^n {x_k^3}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^4}} \right)u_2 = \sum\limits_{k = 0}^n{f_k} x_k ^2. \end{gather*}$$
  • 3.4. Задачи для самостоятельного решения

  • Задана таблица {xk, fk}, k = 0, 1, ..., n, n > 1.

    Найти линейную функцию $$f(x) = \alpha x + \beta,$$ минимизирующую функционал $$\Phi (\alpha,\beta) = \sum\limits_{k = 0}^n {(f_k - \alpha x_k - \beta)^2} .$$

    Используя этот результат, решить переопределенную СЛАУ

    x + y = 3,  x + 3y = 7,  2x - y = 0.2, 3x + y = 5.
  • Пусть $$\mathbf{B}$$ — квадратная матрица размером n x n, $$\mathbf{u}$$ — n -мерный вектор, $$\Phi (x)$$ — функционал, $$\Phi (x) = \left({\mathbf{Bu} - \mathbf{u}x,\mathbf{Bu} - \mathbf{u}x}\right).$$ Доказать, что $$\Phi (x)$$ достигает минимума при $$x = (\mathbf{Bu},\mathbf{u})/(\mathbf{u},\mathbf{u}).$$
  • Барометрическое давление изменяется с высотой по закону p = aebh.

    Определить коэффициенты a, b по результатам наблюдений, приведенных в таблице ( h — высота в метрах над уровнем моря; p — давление в мм рт. ст.):

    h 0 270 840 1452 2116 3203
    p 760 737 686 636 584 508
  • Скорость корабля связана с мощностью его двигателя формулой P = a0 + a3v3 ( P — мощность в лошадиных силах, v — скорость в узлах).

    Определить a0, a3 по табличным данным.

    v 6 8 10 12 13
    P 423 805 1378 2357 2893
  • Найти приближение многочленом третьей степени методом наименьших квадратов для функции $$f(x) = \sin (\pi x)$$ по значениям в точках x0 = - 1, x1 = - 0,5, x2 = 0, x4 = 1.
  • Вычислить матрицы $$\mathbf{A^*BA}$$ и $$\mathbf{A^*B}$$ и найти решение системы уравнений вида $$\mathbf{A^*BAu} = \mathbf{A^*Bf},$$ если заданы матрицы $$\mathbf{A}, \mathbf{B}$$ и вектор $$\mathbf{f}.$$

    $$\mathbf{A} = \left( \begin{array}{cc} 1 1 \\ 2 {- 1} \\ 1 3 \\ 3 1 \end{array} \right), \mathbf{B} = \left( \begin{array}{cc} 1 {0,5} \\ {0,5} 1 \end{array} \right), \mathbf{f} = {(3, 0,2, 7, 5)}^T.$$

  • Функция $$f(x) = \sqrt{1 + {\sin}^2 (x - 1)}$$ приближенно заменяется тригонометрическим многочленом P(x) = a0 + a1sin x + b1cos x + a2sin 2x + b2cos 2x по десяти точкам x0, ..., x9 с помощью метода наименьших квадратов.

    Опишите алгоритм вычисления коэффициентов a0, a1, b1, a2, b2.

  • Доказать, что прямая, проведенная по методу наименьших квадратов, проходит через точку с координатами $$$ {x^{\prime}} = \frac{\sum\limits_i {y_ix_i}}{\sum\limits_i {y_i}}, y^{\prime} = \frac{\sum\limits_i {x_iy_i}}{\sum\limits_i {x_i}} $.$$
  • Страницы:

    3.1. Пример использования метода наименьших квадратов (МНК)

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

    Пусть f — линейная (или близкая к линейной) функция аргумента x: f(x) = u1x + u0. В точках xk известны значения функции f(xk). Тогда u0, u1 — коэффициенты, которые необходимо подобрать так, чтобы выполнялись условия u1xk + u0 = fk, k = 0,1,2,3,4, fk = f(xk).

    Получим систему пяти уравнений относительно двух неизвестных. Это — переопределенная система. Она не имеет классического решения, так как в общем случае не существует прямой, проходящей через все 5 точек (это возможно только тогда, когда какие - либо три уравнения полученной системы линейными преобразованиями сводятся к двум другим — система линейно зависима).

    Рассмотрим общий случай. Пусть коэффициенты {u0, u1} необходимо определить по результатам n + 1 измерения. Введем функцию, равную сумме квадратов невязок rk = u1xk + u0 - fk

    $$\Phi (u_1,u_0) = \sum\limits_{k = 0}^n {r_k^2} = \sum\limits_{k = 0}^n {(u_1 x_k + u_0 - f_k)^2} .$$

    Примем за обобщенное решение переопределенной СЛАУ такие {u0, u1} для которых $$\Phi (u_{0}, u_{1})$$ принимает наименьшие значение. Для определения обобщенного решения из условия минимума суммы квадратов невязки получаем систему двух уравнений, имеющую классическое решение:

    $$$ \frac{\partial \Phi }{\partial u_0} = 0, \frac{\partial \Phi }{\partial u_1} = 0. $$$

    Выбор функции $$\Phi (u_{0}, u_{1})$$ имеет некоторый произвол. Например, возможно каждому измерению придать некоторый вес bk. От набора таких весовых множителей зависело бы решение системы. В этом случае функция $$\Phi$$ будет

    $$\Phi (u_0,u_1) = \sum\limits_{k = 0}^n {b_k{(u_1 x_k - u_0 - f_k)}^2}.$$

    Если в качестве невязки выбрать rk = | u1xk + u0 - fk |, то получим задачу линейного программирования на отыскании минимума функции

    $$\Phi (u_1 u_2) = \sum\limits_{k = 0}^3 \left|{u_1 x_k + u_0 - f_k}\right|.$$

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

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

    $$r_k = \sum\limits_{j = 0}^p {u_j\varphi_j (x_k)} - f(x_k), k = 1, \ldots, n,$$

    где $$\varphi_j (x)$$ — некоторые функции, образующие базис, например, тригонометрические: $$\varphi_j (x) = \sin (jx).$$ Выражение $$\sum\limits_{j = 0}^p {u_j \varphi_j (x)}$$ называется обобщенным полиномом. В приведенном выше примере в качестве базисных функций были выбраны степенные функции $$\varphi_j (x) = x^j .$$ Обобщенный полином превратился в алгебраический.

    В случае выбора произвольной системы базисных функций переопределенная СЛАУ и функционал $$\Phi (u_{0}, \dots , u_{p})$$ будут

    $$\begin{gather*} u_0 \varphi_0 (x_0) + \ldots + u_p \varphi_p(x_0) = f_0, \\ \ldots \\ u_0 \varphi_0 (x_n) + \ldots + u_p\varphi_p (x_n) = f_n,\\ \Phi (u_0,\ldots,u_n) = \sum\limits_{i = 0}^n (\sum\limits_{j = 0}^p u_j \varphi_j(x_i) - f_i)^2 \end{gather*}$$

    Отыщем обобщенное решение методом наименьших квадратов. Приравнивая все частные производные по компонентам обобщенного решения к нулю $$$ \frac{\partial \Phi }{\partial u_k} = 0 $$$ (условия минимума) и изменяя порядок суммирования, получаем СЛАУ

    $$\sum\limits_{j = 0}^p {\left({\sum\limits_{i = 0}^n{\varphi_j (x_i)\varphi_k (x_i)}}\right)u_j = \sum\limits_{i = 0}^n {f_i\varphi_k (x_i)}}, k = 0, \ldots, p,$$

    или

    $$\begin{gather*} (\varphi_0, \varphi_0) u_0 + (\varphi_0, \varphi_1)u_1 + \ldots + (\varphi_0, \varphi_p)u_p = (\varphi_0, f), \\ (\varphi_1, \varphi_0) u_0 + (\varphi_1, \varphi_1)u_1 + \ldots + (\varphi_1, \varphi_p)u_p = (\varphi_1, f), \\ (\varphi_p, \varphi_0) u_0 + (\varphi_p, \varphi_1)u_1 + \ldots + (\varphi_p, \varphi_p)u_p = (\varphi_p, f), \end{gather*}$$

    Система метода наименьших квадратов имеет вид $$\mathbf{Du} = \mathbf{f}$$ с матрицей $$\mathbf{D},$$ элементами которой являются скалярные произведения $$(\varphi_i, \varphi_j) = \sum\limits_{i = 0}^n \varphi_j (x) \varphi_k (x_i)$$ Это — матрица Грама. Ее свойства известны из курса линейной алгебры, эта матрица симметричная и положительно определенная. Таким образом, решение исследуемой СЛАУ существует и единственно. В правой части системы стоят проекции свободного члена исходной задачи на подпространство базисных функций $$(\varphi,f) = \sum\limits_{i = 0}^n {\varphi_j(x_i)f_i}.$$

    Здесь учтено, что $$$ \frac{\partial \Phi }{\partial u_k} = 2 \sum\limits_{i = 0}^n {\varphi_k(x_i)\left({\sum\limits_{j = 0}^p {u_j\varphi_j (x_i) - f_i}} \right)} $,$$ или, в развернутом виде

    $$\begin{gather*} \sum\limits_{i = 0}^n {\varphi_0 (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_0(x_i) + \ldots + u_p \varphi_p(x_i) - f(x_i)} \right) = 0}, \\ \sum\limits_{i = 0}^n {\varphi_1 (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_1 (x_i) + \ldots + u_p \varphi_p (x_i) - f(x_i)} \right) = 0}, \\ \ldots, \\ \sum\limits_{i = 0}^n {\varphi_{s} (x_i)\left({u_0 \varphi_0 (x_i) + u_1 \varphi_1 (x_i) + \ldots + u_p \varphi_p (x_i) - f(x_i)} \right) = 0.} \end{gather*}$$

    Часто выбирают $$\varphi_k (x) = x^k,$$ в этом случае система уравнений принимает следующую форму:

    $$\sum\limits_{j = 0}^p (\sum\limits_{i = 0}^n x_i^{j + k})u_j = \sum\limits_{i = 0}^n f_i x_i^k, k = 1, \ldots, p.$$

    Эта система может быть легко выписана в компонентах:

    $$\begin{gather*} u_0 + \left({\sum\limits_{i = 0}^n {x_i}} \right)u_1 + \ldots + \left({\sum\limits_{i = 0}^n {x_i^p}}\right)u_p = \sum\limits_{i = 0}^n{f(x_i)}, \\ \left({\sum\limits_{i = 0}^n {x_i}} \right)u_0 + \left({\sum\limits_{i = 0}^n {x_i^2}} \right)u_1 + \ldots + \left({\sum\limits_{i = 0}^n {x_i^{p + 1}}} \right) = \sum\limits_{i = 0}^n {x_i f(x_i)}, \\ \ldots \\ \left({\sum\limits_{i = 0}^n {x^p_i}} \right)u_0 + \left({\sum\limits_{i = 0}^n {x^{p + 1}_i}}\right) + \ldots + \left({\sum\limits_{i = 0}^n{x^{2p + 1}_i}}\right)u_p = \sum\limits_{i = 0}^n {x^p_if(x_i)}. \end{gather*}$$

    В случае использования ортонормированных систем базисных функций $$\varphi_j (x),$$ т.е., при выполнении условия $$(\varphi_k, \varphi_j) = \delta_{kj}$$ решение принимает простой вид $$$ u_0 = \frac{(\varphi_0, f)}{(\varphi_0,\varphi_0)} = {\left\| {\varphi_0} \right\|}^{- 2} (\varphi_0, f); u_1 = {\left\| {\varphi_1} \right\|}^{- 2} (\varphi_1, f); \ldots u_p = {\left\| {\varphi_p} \right\|}^{- 2} (\varphi_p, f) $.$$ Поскольку система функции ортонормирована, то $$\left\| {\varphi_j}\right\| = 1$$ и $$u_k = (\varphi_k, f), k = 1, \ldots, n.$$ Эти коэффициенты называются коэффициентами Фурье, а обобщенный многочлен с этими коэффициентами — обобщенным многочленом Фурье. В частности, в качестве базисных можно использовать ортогональную систему функций на $$[- \pi, \pi ]: \{\sin{kx},\cos{kx}\}$$ k = 1, ..., n. Такие представления называют отрезками тригонометрических рядов Фурье или конечными рядами Фурье .

    Докажем теорему о методе наименьших квадратов, обобщающую изложенную информацию.

    Запишем переопределенную СЛАУ

    $$\begin{gather*} a_{11} u_1 + \ldots + a_{1p} u_p = f_1, \\ \ldots \\ a_{n1} u_1 + \ldots + a_{np} u_p = f_n, n > p. \end{gather*}$$

    $$\mathbf{u} = {\{u_1, \ldots,u_p \}}^T \in L^p, \mathbf{f} = {\{f_1, \ldots,f_n \}}^T \in L^n,$$

    где линейные нормированные пространства Lp и Ln имеют размерности p и n соответственно. Перепишем (3.2) в матричной форме:

    $$\mathbf{Au} = \mathbf{f},$$

    где

    $$\mathbf{A} = \left( \begin{array}{ccc} {a_{11}} \ldots {a_{1p}} \\ \ldots \ldots \ldots \\ {a_{n1}} \ldots {a_{np}}} \end{array} \right).$$

    Наряду с основным скалярным умножением в Ln

    $${(\mathbf{x},\mathbf{y})}^n = \sum\limits_{k = 1}^n{x_k y_k}$$

    введем скалярное умножение с весовой матрицей $$\mathbf{B}$$:

    $${\left[{\mathbf{x},\mathbf{y}}\right]}^n = {(\mathbf{Bx},\mathbf{y})}^n, \mathbf{B} = {\mathbf{B}}^* > 0, \mathbf{x}, \mathbf{y} \in L^n$$

    Оба этих умножения удовлетворяют аксиомам скалярного умножения элементов линейного пространства. Матрица $$\mathbf{B}$$ является весовой и определяет вклад невязки каждого слагаемого суммы (3.1). Система (3.2) не имеет классического решения. Определим обобщенное решение этой системы как элемент линейного пространства $$\mathbf{v},$$ придающий наименьшее значение квадратичной форме:

    $$\Phi (\mathbf{u}) = {\left[\mathbf{Au} - \mathbf{f}, \mathbf{Au} - \mathbf{f}\right]}^n.$$

    Теорема. Пусть столбцы матрицы $$\mathbf{A}$$ линейно независимы, т.е. ранг $$\mathbf{A}$$ равен p. Тогда существует единственный элемент p - мерного евклидова пространства $$\mathbf{v} \in L^P,$$ являющийся обобщенным решением системы (3.2), и решением СЛАУ

    $${\mathbf{A}}^*{\mathbf{BAu}} = {\mathbf{A}}^*{\mathbf{Bf}},$$

    состоящей из p скалярных уравнений относительно неизвестных $$\left\{{\mathbf{u}_k}\right\}_{k = 1}^{k = p},$$ доставляющий минимум квадратичной форме

    $$\Phi (\mathbf{u}) = {\left[\mathbf{Au} - \mathbf{f}, \mathbf{Au} - \mathbf{f}\right]}^n.$$

    Доказательство.

    Покажем, что решение СЛАУ $${\mathbf{A}}^*{\mathbf{BAu}} = {\mathbf{A}}^*{\mathbf{Bf}}$$ существует и единственно. Введем обозначение $${\mathbf{q}}_k \in L^n$$ для вектора — k столбца матрицы системы $$\mathbf{A}$$:

    $${\mathbf{q}}_k = {\{a_{1k}, \ldots, a_{nk} \}}^T; k = 1, \ldots, p.$$

    Несложно показать, что матрица $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ системы (3.5) есть квадратная матрица p x p. Элемент dij этой матрицы, стоящий на пересечении i строки и j столбца, есть $$d_{ij} = {({\mathbf{q}}_i, {\mathbf{Bq}}_j)}^n = {({\mathbf{Bq}}_i,{\mathbf{q}}_j)}^n = {\left[{{\mathbf{q}}_i,{\mathbf{q}}_j}\right]}^n,$$ в силу коммутативности скалярного произведения dij = dij, что означает самосопряженность матрицы $$\mathbf{D}: \mathbf{D} = {\mathbf{D}}^*.$$

    Покажем, что матрица $$\mathbf{D}$$ невырождена и положительно определена. Напомним, что $${(\mathbf{f},\mathbf{A\xi })}^n = {({\mathbf{A}}^*{\mathbf{f}},{\mathbf{\xi }})}^P; {\mathbf{f}}_1 \in L^n, \mathbf{\xi } \in L^P.$$ Это равенство проверяется непосредственно, если записать обе его части в развернутом виде. Невырожденность матрицы $$\mathbf{D}$$ следует из того, что ранг матрицы $$\mathbf{A}$$ равен p.

    В таком случае $$0 < {\left[{\mathbf{Ax},\mathbf{x}}\right]}^n = {(\mathbf{BAx},\mathbf{Ax})}^n = {({\mathbf{A}}^*{\mathbf{BAx}},{\mathbf{x}})}^p = {(\mathbf{Dx},\mathbf{x})}^p.$$ Поскольку $$\mathbf{D}$$ невырождена и положительно определена, то (3.5) имеет единственное решение $$\mathbf{v} \in L^p.$$ Теперь покажем, что $$\mathbf{v}$$ — единственное обобщенное решение системы. Для любого вектора $$\mathbf{\Delta } \ne 0$$ выполнено $$\Phi (\mathbf{v} + \mathbf{\Delta }) > \Phi (\mathbf{v})$$:

    $$\Phi (\mathbf{v} + \mathbf{\Delta }) = {\left[{\mathbf{A}(\mathbf{v} + \mathbf{\Delta }) - \mathbf{f},\mathbf{A}(\mathbf{v} + \mathbf{\Delta }) - \mathbf{f}}\right]}^n = \\ = {\left[{\mathbf{Av} - \mathbf{f},\mathbf{Av} - \mathbf{f}}\right]}^n - 2{\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n + {\left[{\mathbf{A\Delta },\mathbf{A\Delta }}\right]}^n = \\ = \Phi (\mathbf{v}) + {2\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n + {(\mathbf{D\Delta },\mathbf{\Delta })}^p = \Phi (\mathbf{v}) + {(\mathbf{D\Delta },\mathbf{\Delta })}^p > \Phi (\mathbf{v}),$$

    что и требовалось доказать.

    При доказательстве использовалось

    $${\left[{\mathbf{Av} - \mathbf{f},\mathbf{A\Delta }}\right]}^n = {(\mathbf{B}(\mathbf{Av} - \mathbf{f}),\mathbf{A\Delta })}^n = {({\mathbf{A}}^*{\mathbf{BAv}} - {\mathbf{A}}^*{\mathbf{Bf}},\mathbf{\Delta })}^p = 0,$$

    поскольку $${\mathbf{A}}^*{\mathbf{BAv}} = {\mathbf{A}}^*{\mathbf{Bf}}.$$

    Так как матрица $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ — симметричная и положительно определенная, то для численного решения полученной СЛАУ можно воспользоваться итерационными методами.

    Если система векторов $$\left\{{{\mathbf{q}}_k}\right\}_{k = 1}^p$$ оказывается ортонормированной, т.е. $$[{\mathbf{q}}_i,{\mathbf{q}}_j ] = \delta_{ij} i,j = 1, \ldots, p,$$ то матрица $$\mathbf{D}$$ оказывается единичной. Ее элементы и есть скалярные произведения $$[{\mathbf{q}}_i,{\mathbf{q}}_j ].$$ В этом случае решением системы будет

    $$\mathbf{u} = \mathbf{A^*Bf}.$$

    Следует отметить, что если базисные функции $$\varphi_i (x), i = 0, \ldots, p$$ не выбираются специальным образом, то при достаточно больших p $$(p \ge 5)$$ полученная СЛАУ оказывается плохо обусловленной. Строки матрицы $$\mathbf{D} = {\mathbf{A}}^*{\mathbf{BA}}$$ могут оказаться почти линейно зависимыми. Простейшим примером такого почти линейно зависимого базиса является система функций xi, i = 1, ..., p при больших p. В этом случае желательно использовать ортогональные функциональные базисы, однако такой выбор не всегда возможен и удобен.

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

    $$f(x) = \sum\limits_{j = 0}^p {u_j x^j}.$$

    Скалярные произведения на отрезке [0, 1], записанные в интегральной форме (т.е. при $$n \to \infty$$ ), будут иметь вид

    $$$ (\varphi_i,\varphi_j) = \int\limits_0^1 x^i x^j dx = \int\limits_0^1 x^{i + j}dx = \frac{1}{i + j + 1}. $$$

    В таком случае СЛАУ после применения МНК, т.е. минимизации функционала $$\Phi (u) = \int\limits_0^1 {(F(x) - f(x))}^2dx,$$ где F(x) — заданная функция, будет:

    $$\begin{gather*} u_0 + \frac{1}{2}u_1 + \ldots + \frac{1}{p + 1}u_p = \int\limits_0^1{F(x)dx}, \\ \frac{1}{2}u_0 + \frac{1}{3}u_1 + \ldots + \frac{1}{p + 2}u_p = \int\limits_0^1 xF(x)dx, \\ \ldots \\ \frac{1}{p + 1}u_0 + \frac{1}{p + 2}u_1 + \ldots + \frac{1}{2p + 1}u_p = \int\limits_0^1 pF(x)dx, \end{gather*}$$

    или

    $${\mathbf{H}}_{p + 1}\mathbf{u} = \mathbf{F^{\prime}},$$

    где

    $$\begin{gather*} \mathbf{u} = {\{u_0, \ldots, u_p \}}^T,\quad \mathbf{F^{\prime}} = {\left\{{\int\limits_0^1 {F(x)dx}, \ldots, \int\limits_0^1 {x^p F(x)dx}} \right\}}^T, \\ {\mathbf{H}}_{p + 1} = {\left\{\frac{1}{i + j - 1}\right\}}^{p + 1}_{i,j} = 1 \end{gather*}$$

    Матрица $$\mathbf{H}_{p + 1}$$ называется матрицей Гильберта. Это классический пример плохо обусловленной матрицы. Число обусловленности очень быстро растет с ростом p. Так при p = 1 $$\mu = \left\|\mathbf{H}\right\| \cdot \left\|{\mathbf{H}}^{- 1}\right\| \approx 20,$$ при $$p = 9 \mu \approx 10^{13}.$$ Если получим СЛАУ для дискретной системы точек, т.е. для $$$ (\varphi_j,\varphi_k) = \sum\limits_{i = 0}^n x^{k + j}_i, (\varphi_j,f) = \sum\limits_{i = 0}^n x^k_i f(x_i)x_i = \frac{i}{n} $,$$ то ее матрица будет асимптотически приближаться к матрице Гильберта $$\mathbf{H}_{p + 1}$$ при $$n \rightarrow \infty.$$

    3.2. Понятие о методах решения плохо обусловленных СЛАУ

    Улучшить качество численного решения СЛАУ метода наименьших квадратов возможно, если использовать различные преобразования матрицы $$\mathbf{A}.$$

    Большинство прямых методов решения линейных систем основано либо на замене исходной системы $$\mathbf{Au} = \mathbf{f}$$ ( $$\mathbf{A}$$ — квадратная матрица) на эквивалентную $$\mathbf{CAu} = \mathbf{Cf},$$ либо на представлении матрицы $$\mathbf{A}$$ в виде произведения других матриц, таких, чтобы новая система либо решалась более просто, либо была лучше (по крайней мере, не хуже) обусловлена, чем исходная.

    Подход, использующий спектральную эквивалентность матриц $$\mathbf{A}$$ и $${\mathbf{C}}^{- 1}$$ (в смысле границ спектра собственных значений), основан на умножении $$\mathbf{A}$$ на близкую в некотором смысле матрицу $$\mathbf{C}.$$ Последняя матрица выбирается таким образом, чтобы произведение было близким к единичной матрице (при этом $${\mathbf{C}} \ne {\mathbf{A}}^{- 1},$$ так как обращение плохо обусловленной матрицы приводит к накоплению вычислительных ошибок. Число обусловленности матрицы $$\mathbf{CA}$$ $$\mu = \lambda _{max}/\lambda _{min}$$ будет близко к единице. Метод энергетически эквивалентных операторов оказался эффективным при численном решении сеточных уравнений при разностной аппроксимации уравнений в частных производных эллиптического типа.

    Идея предобусловливания СЛАУ состоит в том, чтобы вместо исходной системы $$\mathbf{Au} = \mathbf{f}$$ решать систему $$\mathbf{A^{\prime}u^{\prime}} = \mathbf{f^{\prime}},$$ где $$\mathbf{A^{\prime}} = \mathbf{C^{- 1}AC^{- 1}}, \mathbf{u^{\prime}} = \mathbf{Cu}, \mathbf{f^{\prime}} = \mathbf{C^{- 1}f},$$ матрица $$\mathbf{C}$$ выбирается так, чтобы она была симметричной положительной, хорошо обусловленной.

    Определение. Матрица $$\mathbf{Q}$$ с вещественными элементами qij является ортогональной, если

    $$\mathbf{Q}^* = \mathbf{Q}^{- 1}.$$

    Пусть матрица $$\mathbf{C}$$ невырождена. Тогда она представима в виде

    $$\mathbf{C} = \mathbf{QR},$$

    где $$\mathbf{Q}$$ — ортогональная, а $$\mathbf{R}$$ — верхняя треугольная матрицы. В качестве матрицы $$\mathbf{Q}$$ часто используется симметричная ортогональная матрица $$\mathbf{H}$$:

    $$\mathbf{H} = \mathbf{E} - 2{\mathbf{ww}}^T,$$

    где $$\mathbf{w}$$ - произвольный вектор - столбец, такой, что $$({\mathbf{w}}^T{\mathbf{w}}) = 1.$$ Матрица $$(\mathbf{w,wT})$$ есть произведение вектора - столбца $$\mathbf{w}$$ на вектор - строку $$\mathbf{wT}$$ (преобразование Хаусхолдера или метод отражений ).

    Заметим, что симметричность матрицы $$\mathbf{H}$$ (матрицы отражений) показывается непосредственной проверкой; ортогональность $$\mathbf{H}$$ можно показать следующим образом:

    $${\mathbf{HH}}^T = (\mathbf{E} - 2{\mathbf{ww}}^T)(\mathbf{E} - 2{\mathbf{ww}}^T) = \\ = \mathbf{E} - 2{\mathbf{ww}}^T - 2{\mathbf{ww}}^T + 4{\mathbf{ww}}^T \cdot {\mathbf{ww}}^T = \\ = \mathbf{E} - 4{\mathbf{ww}}^T - 4\mathbf{w}({\mathbf{w}}^T{\mathbf{w}}){\mathbf{w}}^T = \mathbf{E}, \\ $$

    так как $${\mathbf{w}}^T{\mathbf{w}} = 1.$$

    Матрица СЛАУ представляется в виде $$\mathbf{A} = {\mathbf{H}}^T{\mathbf{R}},$$ после чего решается эквивалентная система уравнений $$\mathbf{Ru} = \mathbf{Hf}.$$

    3.3. Задачи

  • Для функции $$f(x) = \sqrt {x}$$ на отрезке [0,1] построить многочлен F(x) = u0 + u1x среднеквадратичного приближения со скалярными произведениями:

    $$(\varphi_i,\varphi_j) = \int\limits_0^1 {x^ix^j}{dx},\quad f,\varphi_i) = \int\limits_0^1 f(x)x^idx.$$

    Решение. Введем базисные функции $$\varphi_0 (x) = 1, \varphi_1 (x) = x,$$ и вычислим скалярные произведения

    $$$ (\varphi_0,\varphi_0) = \int\limits_0^1 {1^2 dx} = 1, (\varphi_1,\varphi_1) = \int\limits_0^1 {x^2 dx =} \frac{1}{3},$$$

    $$$ (\varphi_0,\varphi_1) = \int\limits_0^1 {xdx}} = \frac{1}{2}, (f,\varphi_0) = \int\limits_0^1 {\sqrt{1} x dx} = \frac{2}{3},$$$

    $$$ (f,\varphi_1) = \int\limits_0^1 {\sqrt{x} x dx} = \frac{2}{5},$$$

    Для вычисления коэффициентов получим СЛАУ

    $$$ u_0 + \frac{1}{2}u_1 = \frac{2}{3},$$$

    $$$ \frac{1}{2}u_0 + \frac{1} 3u_1 = \frac{2}{5}, $$$

    откуда $$$ u_0 = \frac{4}{{15}}, c_1 = \frac{4}{5}, F(x) = \frac{4}{15} + \frac{4}{5}x.$$$

  • Получить СЛАУ
    c11x + c12y = f1
    c22x + c21y = f2,

    если она задана в форме метода наименьших квадратов: $$\mathbf{A^*BAu} = \mathbf{A^*Bf},$$ где

    $$\mathbf{A} = \left( \begin{array}{cc} {a_{11}} {a_{12}} \\ {a_{21}} {a_{22}} \\ {a_{31}} {a_{32}} \end{array} \right), \mathbf{B} = \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right), \mathbf{u} = {(u_1,u_2)}^T, \mathbf{f} = {(f_1,f_2,f_3)}^T .$$

    Решить эту систему для случая a11 = a32 = 1, a21 = a12 = 2, a31 = 0, a22 = 1, f1 = 1, f2 = 2, f3 = 1.

    Решение. Проведем необходимые вычисления:

    $$\begin{gather*} \left( \begin{array}{ccc} {a_{11}} {a_{21}} {a_{31}} \\ {a_{12}} {a_{22}} {a_{32}} \end{array} \right) \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right) \left( \begin{array}{cc} {a_{11}} {a_{12}} \\ {a_{21}} {a_{22}} \\ {a_{31}} {a_{32}} \end{array} \right) = \\ = \left( \begin{array}{ccc} {a_{11}} {a_{21}} {a_{31}} \\ {a_{12}} {a_{22}} {a_{32}} \end{array} \right) \left( \begin{array}{cc} 1 0 \\ 0 1 \end{array} \right) \left( \begin{array}{l} {f_1} \\ {f_2} \\ {f_3} \end{array} \right), \\ \left( \begin{array}{cc} {a_{11}^2 + a_{21}^2 + a_{31}^2} {a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}} \\ {a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}} {a_{12}^2 + a_{22}^2 + a_{32}^2} \end{array} \right) = \\ = \left( \begin{array}{l} {u_1} \\ {u_2} \end{array} \right) = \left( \begin{array}{l} {a_{11}f_1 + a_{21}f_2 + a_{31}f_3} \\ {a_{12}f_1 + a_{22}f_2 + a_{32}f_3} \\ \end{array} \right), \\ (a_{11}^2 + a_{21}^2 + a_{31}^2)u_1 + (a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32})u_2 = a_{11}f_1 + a_{21}f_2 + a_{31}f_3, \\ (a_{11}a_{12} + a_{21}a_{22} + a_{31}a_{32}) u_1 + (a_{12}^2 + a_{22}^2 + a_{32}^2)u_2 = a_{12}f_1 + a_{22}f_2 + a_{32}f_3. \end{gather*}$$

    Решение в числах предлагается найти читателям.

  • С помощью метода наименьших квадратов найти коэффициенты полинома второй степени f(x) = u0 + u1x + u2x2, если таблица измерений задана $$\left. {\left\{{x_i, f_i}\right\}}\right|_{i = 0}^n .$$

    Решение. Переопределенная система уравнений:

    $$u_0 + u_1 x_k + u_2 x_k^2 = f_k, k = 0, 1, \ldots, n,$$

    определяет функционал

    $$\Phi (u_0, \ldots, u_n) = \sum\limits_{k = 0}^n{r_k^2} = \sum\limits_{k = 0}^n {(u_0 + u_1 x_k + u_2 x_k^2 - f_k)}^2.$$

    Условия минимума последнего

    $$$ \frac{\partial \Phi }{\partial u_k} = 0. $$$

    После проведения алгебраических преобразований получим для коэффициентов систему линейных уравнений

    $$\begin{gather*} u_0 + \left({\sum\limits_{k = 0}^n {x_k}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^2}}\right)u_2 = \sum\limits_{k = 0}^n {f_k}, \\ \left({\sum\limits_{k = 0}^n {x_k}} \right)u_0 + \left({\sum\limits_{k = 0}^n {x_k^2}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^3}} \right)u_2 = \sum\limits_{k = 0}^n {f_k} x_k, \\ \left({\sum\limits_{k = 0}^n {x_k^2}} \right)u_0 + \left({\sum\limits_{k = 0}^n {x_k^3}}\right)u_1 + \left({\sum\limits_{k = 0}^n {x_k^4}} \right)u_2 = \sum\limits_{k = 0}^n{f_k} x_k ^2. \end{gather*}$$
  • 3.4. Задачи для самостоятельного решения

  • Задана таблица {xk, fk}, k = 0, 1, ..., n, n > 1.

    Найти линейную функцию $$f(x) = \alpha x + \beta,$$ минимизирующую функционал $$\Phi (\alpha,\beta) = \sum\limits_{k = 0}^n {(f_k - \alpha x_k - \beta)^2} .$$

    Используя этот результат, решить переопределенную СЛАУ

    x + y = 3,  x + 3y = 7,  2x - y = 0.2, 3x + y = 5.
  • Пусть $$\mathbf{B}$$ — квадратная матрица размером n x n, $$\mathbf{u}$$ — n -мерный вектор, $$\Phi (x)$$ — функционал, $$\Phi (x) = \left({\mathbf{Bu} - \mathbf{u}x,\mathbf{Bu} - \mathbf{u}x}\right).$$ Доказать, что $$\Phi (x)$$ достигает минимума при $$x = (\mathbf{Bu},\mathbf{u})/(\mathbf{u},\mathbf{u}).$$
  • Барометрическое давление изменяется с высотой по закону p = aebh.

    Определить коэффициенты a, b по результатам наблюдений, приведенных в таблице ( h — высота в метрах над уровнем моря; p — давление в мм рт. ст.):

    h 0 270 840 1452 2116 3203
    p 760 737 686 636 584 508
  • Скорость корабля связана с мощностью его двигателя формулой P = a0 + a3v3 ( P — мощность в лошадиных силах, v — скорость в узлах).

    Определить a0, a3 по табличным данным.

    v 6 8 10 12 13
    P 423 805 1378 2357 2893
  • Найти приближение многочленом третьей степени методом наименьших квадратов для функции $$f(x) = \sin (\pi x)$$ по значениям в точках x0 = - 1, x1 = - 0,5, x2 = 0, x4 = 1.
  • Вычислить матрицы $$\mathbf{A^*BA}$$ и $$\mathbf{A^*B}$$ и найти решение системы уравнений вида $$\mathbf{A^*BAu} = \mathbf{A^*Bf},$$ если заданы матрицы $$\mathbf{A}, \mathbf{B}$$ и вектор $$\mathbf{f}.$$

    $$\mathbf{A} = \left( \begin{array}{cc} 1 1 \\ 2 {- 1} \\ 1 3 \\ 3 1 \end{array} \right), \mathbf{B} = \left( \begin{array}{cc} 1 {0,5} \\ {0,5} 1 \end{array} \right), \mathbf{f} = {(3, 0,2, 7, 5)}^T.$$

  • Функция $$f(x) = \sqrt{1 + {\sin}^2 (x - 1)}$$ приближенно заменяется тригонометрическим многочленом P(x) = a0 + a1sin x + b1cos x + a2sin 2x + b2cos 2x по десяти точкам x0, ..., x9 с помощью метода наименьших квадратов.

    Опишите алгоритм вычисления коэффициентов a0, a1, b1, a2, b2.

  • Доказать, что прямая, проведенная по методу наименьших квадратов, проходит через точку с координатами $$$ {x^{\prime}} = \frac{\sum\limits_i {y_ix_i}}{\sum\limits_i {y_i}}, y^{\prime} = \frac{\sum\limits_i {x_iy_i}}{\sum\limits_i {x_i}} $.$$
  • Вернуться к учебному плану