авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 8 |

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

-- [ Страница 5 ] --

4/5 3/5 12/5 12/ 4 В любом варианте QT AQ2 =, где =, что и должно быть.

Теперь найдем + и A+. Имеем (из примеров 10.1, 10.2 и 10.3): + = = 1/5 | 0. Поэтому ищем A+ = Q2+QT. В варианте 3/5 4/ A+ = Q2 + QT = 1 1/5 0 = 3/25 4/25.

4/5 3/ 10 Теоретические основы В варианте 3/5 4/ A+ = Q2 + QT = 1 1/5 0 = 3/25 4/25.

4/5 3/ Здесь также имеем A+ A = 3/25 4/25 = 1 = I.

Это произошло в силу того, что rank A = r = n = 1. Как мы знаем, в этом случае A+ = (AT A)1AT = 25 3 4 = 3/25 4/25.

Вывод из примера 10.4. Пользоваться сингулярным разложением A = Q1 QT, чтобы находить A+, непросто. Для этого нужно:

1. Найти собственные векторы x1, x2,..., xn матрицы AT A.

Тогда Q2 = x1 x2... xn.

2. Найти собственные значения матрицы AT A: 1 2... n 0.

3. Оценить ранг r матрицы A, т. е. разграничить {1 2.

..

r 0} и {r+1 = r+2 =... = n = 0} и образовать µi = i, D = diag (µ1, µ2,..., µr ).

4. Найти yi = Axi/µi, где µi = i, i = 1, 2,..., r.

5. Доопределить каким-либо образом систему {yi } до ортонормирован ного базиса в Rm.

6. Образовать Q1 = y1 y2... ym.

В качестве альтернативного решения для нахождения A+ рассмотрим на следующем примере так называемое усеченное LU -разложение матрицы.

Пример 10.5. Дана матрица [13] 1 A=2 5, rank A = 2 = n.

3 Выполним полное LU -разложение и затем усеченное LU -разложение:

12 1 0 1 = 2 1 1 2 = LU, A=2 5=2 37 311 U L U L 10.5 Отыскание псевдообратной матрицы rank L = rank U = rank A = 2.

В общем случае: пусть rank A = r, A = LU, тогда:

1. L всегда квадратная с единичной диагональю, т. е. L = L(m, m) и det L = 1, имеет (m r) нулевых нижних строк.

2. U Отбросим последние (m r) строк в U и последние (m r) столбцов в L (в этом и заключается переход к усеченному разложению), тогда получим:

A = LU, L = L(m, r), U = U (r, n), rank L = r = rank U = rank A.

Свойства L и U (так как A = LU и rank A = rank L = r):

1 U : R(U T ) = R(AT ) совпадение пространств строк, 2 L: R(L) = R(A) совпадение пространств столбцов.

Теорема 10.10 ( [13]).

A+ = U T (U U T )1(LT L)1LT.

произвольный вектор в Rm. Рассмотрим Доказательство. Пусть b вектор y = [ U T (U U T )1(LT L)1LT ] b.

Вектор y R(U T ) = R(AT ). Умножим y слева на A = LU :

U T (U U T )1(LT L)1LT b = L(LT L)1LT b.

Ay = LU По определению, L(LT L)1LT, есть матрица проектирования на R(L), т. е.

на R(A), поэтому U T (U U T )1(LT L)1LT A b = p, проекция b на R(A).

где p Таким образом, вектор y удовлетворяет следующим условиям:

1. y (AT ).

проекция любого вектора b на R(A).

2. Ay = p, где p 10 Теоретические основы 3. y = U T (U U T )1(LT L)1LT b.

Отсюда U T (U U T )1(LT L)1LT = A+. Приведем заключительные примеры для лучшего понимания структуры и свойств псевдообратной матрицы A+.

Пример 10.6. Если A = A(1, 1), то 0, если A = 0, A+ = 1/A, если A = 0.

Пример 10.7. Если A = diag (1, 2,..., n ), то 0, если j = 0, A+ = diag (+, +,..., +), где + = n 1 2 j 1/j, если j = 0.

Пример 10.8. Матрицы 40 A1 = и A2 = 0 почти одинаковы. Однако, их псевдообратные матрицы 1/4 0 1/4 A+ = A+ = и 1 0 сильно отличаются.

Эти примеры, имеющие общее объяснение в теореме о сингулярном раз ложении, говорят о том, что операция псевдообращения разрывная. Дей ствительно, функция 1/, если = 0, + = 0, если = имеет разрыв в точке = 0.

Отсутствие непрерывности операции псевдообращения приводит к серьез ным вычислительным трудностям и, возможно, к большим вычислительным ошибкам при решении задачи МНК, особенно, если r n или r m.

К счастью, на практике задача МНК обычно возникает при m n = r, т. е. с матрицами A = A(m, n) максимального столбцового ранга для пере определенных систем Ax = z. В этом случае A+ = (AT A)1AT и МНК решение x единственно и равно x0.

10.6 Основные теоремы по МНК и псевдоинверсии 10.6 Основные теоремы по МНК и псевдоинверсии Пусть z Rm и A = A(m, n). Тогда:

Теорема 10.11 ( [1]).

(1) найдется x0 Rn единственный, если это вектор с минимальной нор мой, минимизирующий z Ax 2, (2) вектор x0 является единственным вектором из R(AT ), удовлетворяю проекция вектора z на R(A).

щим уравнению A = z, где z = p x Теорема 10.12 ( [1]). Среди всех векторов x, минимизирующих zAx, вектор x0, имеющий минимальную норму, является единственным вектором вида x0 = AT y, удовлетворяющим уравнению AT A0 = AT z.

x Иными словами, вектор x0 может быть получен с помощью любого век тора y0, удовлетворяющего уравнению AT AAT y = AT z, по формуле x0 = AT y.

Теорема 10.13 ( [1]). Для всякой матрицы A = A(m, n) существует псевдообратная матрица A+ = A+ (n, m), такая что для произвольного век тора z Rm x 0 = A+ z является вектором с минимальной нормой среди всех векторов x, миними зирующих z Ax.

Теорема 10.14 ( [1]). Для любой матрицы A = A(m, n) псевдообрат ная матрица A+ обладает свойствами:

(1) A+ = (AT A)+AT, (6) (AAT )+ = (AT )+A+, (2) (AT )+ = (A+)T, (7) AA+ A = A, (3) A+ = AT (AAT )+, (8) A+ AA+ = A+, (4) (A+)+ = A, (9) AA+ = (AA+)T, (5) (AT A)+ = A+ (AT )+, (10) A+ A = (A+A)T.

10 Теоретические основы Основное правило обращения произведения матриц, (BA)1 = A1B 1, не выполняется для псевдообратных матриц, то есть (BA)+ = A+B +.

Теорема 10.15 ( [1]). Для любой симметрической матрицы A = = A(n, n) с действительными элементами предельная матрица PA = lim(A + I)1 A = lim A(A + I) 0 существует. Она является матрицей проектирования на R(A) = R(AT ). Это означает, что для любого вектора z Rn вектор z = PA z является проекцией z на R(A) = R(AT ) на пространство строк или столб цов данной симметрической матрицы.

Теорема 10.16 ( [1]).

(1) Для произвольных z Rm и A = A(m, n) вектор x минимизирует кри терий z Ax тогда и только тогда, когда x имеет вид:

x = A+ z + (I A+ A)y для некоторого y Rn.

(2) Вектор x, минимизирующий z Ax 2, является единственным тогда и только тогда, когда A+ A = I. Последнее условие выполняется тогда и только тогда, когда только нулевой вектор составляет ядро (нуль пространство) матрицы A, т. е. при r = n.

(3) Уравнение Ax = z имеет решение тогда и только тогда, когда A+ Az = = z. Это выполняется в том и только том случае, когда z R(A).

Вектор x является решением уравнения Ax = z тогда и только тогда, когда он задается в виде x = A+ z + (I A+ A)y для произвольного y Rn. Это решение единственно (и тогда оно равно A+z), если и только если AA+z = z A+A = I.

и 10.7 Вычисление матриц проектирования 10.7 Вычисление матриц проектирования Дана матрица A. Возьмем B = AT. Имеем R(AT ) пространство строк матрицы A, R(AT ) пространство столбцов матрицы AT.

R(B) Определение 10.8. Матрица проектирования P на R(B) есть такая матрица, которая обладает свойстом:

(b P b) R(B), где P b = проекция вектора b на R(B), (b P b) = b b перпендикуляр к R(B).

Запишем = b P b = (I P )b, b матрица проектирования на R (B) = N (B T ) = N (A).

где I P Из определения псевдообратной матрицы (см. определение 10.5) следует, что матрица проектирования P вектора b на R(A) в общем случае такова:

P b = A0 = AA+b P = AA+.

P b = p, x = Если взять B = AT и проектировать вектор b на R(B), то мы должны взять матрицу проектирования в виде P = BB +. Но B + = (AT )+ = (A+)T, следовательно P = AT (A+)T.

Если A = A(m, n), то P имеет размеры: (n m)(n m)T = (n n). Если B = B(n, m), то B + = B + (m, n), тогда P имеет размеры: (n m)(m n) = = (n n).

Выводы:

1. Если ищут матрицу вида P = I AT (AAT )1A, то это есть матрица проектирования любого вектора на ядро, т. е. на нуль-пространство N (A) матрицы A. Но в таком виде ее можно определить, только если (AAT )1 существует, т. е. если rank A = m (A имеет полный строчный ранг).

2. Если rank A = r m, что возможно иногда при m n, то (AAT ) не существует. В этом случае для этой же матрицы P справедливо 10 Теоретические основы наиболее общее выражение, а именно: P = I AT (A+)T, где также можно иметь в виду, что всегда (A+)T = (AT )+. Это означает, что нужно уметь отыскивать A+. Для этого есть различные, не очень простые, вычислительные методы (см. книгу [1], а также книгу [13]).

Пример 10.9. Определение проектора PA = AA+. Дано:

1 1 0 A = A1 = 2 2 1 3 = A(n, m) = A(3, 4), 1 1 2 n = 3, m = 4, r = rank(A) = 2 n = 3.

Делаем LU -разложение:

1 1 0 A= 2 1 0 0 1 1.

1 2 1 0 0 L U Находим усеченное LU -разложение:

1 1 0 A = LU = 2 1.

0 0 1 1 U L Отсюда следует:

10 14 1 10 14.

A+ = U + L+ = U T (U U T )1 (LT L)1LT = 5 8 U+ L+ 5 22 Найдем матрицы проектирования как на пространство строк матрицы A, т. е. на R(AT ), так и на пространство столбцов матрицы A, т. е. на R(A).

1. Матрица проектирования на R(AT ):

1 2 10 10 1 2 1 PA = AT (A+)T = 14 14 8 22 = 0 1 2 22 22 41 10.7 Вычисление матриц проектирования 60 60 30 30 2 2 1 1 60 60 30 30 1 2 2 1 1 =.

= 150 30 30 90 60 5 1 1 3 2 30 30 60 90 1 1 2 3 2 1 1 2 3 1 I PA = I AT (A+)T =.

5 1 1 2 1 1 2 PA проектирует на R(AT ), Здесь:

I PA проектирует на N (A).

Найдем базис пространства N (A). Он образован следующими двумя векторами (так как n r = 2):

1 b1 = b2 = и 1.

0 Любой вектор в N (A) задается в виде:

1 y1 + y2 1, где y1, y2 числа.

0 2. В то же время матрица проектирования на R(A):

10 14 1 1 0 1 2 2 1 3 1 10 14 22 = + PA = AA = 150 5 8 1 1 2 5 22 25 50 25 5 10 1 50 130 10 = 10 26 2.

= 150 25 10 145 5 2 25 10 + 4 2.

I PA = I AA = 5 2 10 Теоретические основы PA проектирует на R(A), Здесь:

I PA проектирует на N (AT ) = R(A).

Найдем базис пространства R(A):

1 v1 = 2, v2 = 1.

1 Эти два вектора оказались взаимно ортогональны:

v1 v2 = 1 2 1 1 = 2 + 2 = 0.

T Найдем базис пространства N (AT ). Это пространство определяется как совокупность векторов y, таких что: y T A = 0. Или иначе: y T v1 = 0 и y T v2 = 0. Отсюда найдем:

y = 2, если выбрать y3 = 1.

Действительно, 1 5 2 1 2 = 0 и 5 2 1 1 = 0.

1 Базис пространства N (AT ) состоит из одного вектора. Возьмем в каче стве базисного вектора v3 = 2 y, т. е. вектор 2. v3 = 1.

0. В качестве произвольного вектора для проектирования возьмем (рис. 10.3) 6 z = 0 = 6 0.

0 10.7 Вычисление матриц проектирования e N (AT ) v3 y 1 v2 R(A) e v1 = z перпендикуляр на R(A) z e Рис. 10.3. Проектирование вектора на линейное подпространство Найдем его проекцию на R(A):

5 10 5 1 10 26 2 6 0 = 2 = v1 (совпало с v1).

z = PA z = 5 2 29 Найдем его проекцию на N (AT ):

6 1 z = z z = 0 2 = 2 = y (совпало с y).

0 Пример 10.10.

1 2 1 2 A = A2 = = A(n, m) = A(4, 3), 0 1 n = 4, m = 3, r = rank(A) = 2 n.

По сравнению с примером 10.9 данная матрица совпадает с транспони рованной матрицей примера 10.9. Используем свойство (A+)T = (AT )+ для 10 Теоретические основы нахождения псевдообратной матрицы для нашего примера:

10 10 A+ = 14 14 8 22.

22 22 41 Отсюда, найдем матрицу проектирования на пространство строк матрицы A, т. е. на R(AT ):

10 14 1 1 0 1 10 14 = T +T PA = A (A ) = 2 2 1 150 5 8 1 1 2 5 22 25 50 25 5 10 1 50 130 10 = 10 26 2.

= 150 25 10 145 5 2 Следовательно, 25 10 I PA = I AT (A+)T = 4 2.

5 2 Здесь:

PA проектирует на пространство строк матрицы A, R(AT ), I PA проектирует на нуль-пространство матрицы A, N (A).

Геометрическая интерпретация. Для примера 10.9 построить картинку невозможно (так как там b R4), а для примера 10.10 можно (в этом слу чае b R3 ). Существует связь матриц двух последних примеров, а именно:

A1 = AT. Соответственно, имеем R(AT ) = R(A2 ) R(AT ) = R(A1) 1 N (A1) = N (AT ) N (A2) = N (AT ) 2 Из примера 10.9 имеем матрицу P1, которая проектирует на N (A1), и мат рицу I P1, которая проектирует на R(AT ). Из примера 10.10 имеем матрицу P2, которая проектирует на N (AT ), и матрицу I P2, которая проектирует на R(A1) (см. рис. 10.2 на стр. 213).

10.8 Рекурсия в задаче МНК 10.8 Рекурсия в задаче МНК Постановка вопроса. Даны матрица A = A(m, n) и вектор z Rm.

Требуется найти единственный вектор x0 с минимальной нормой, миними зирующий z Ax 2.

Можно ли искать его последовательно?

1 Известно, что A0 = z, где z R(A), z z R(A), x0 R(AT ).

x Нормальное псевдорешение x0 несовместной системы Ax = z удовле творяет равенству x 0 = A+ z, где A+ псевдообратная матрица.

A 2 Произвольно расщепим матрицу A на блоки, A = и, соответ A z ственно, вектор z на подвекторы, z =, так чтобы A1 = A1 (k, n), z A2 = A2(s, n), k + s = m, z1 R, z2 Rs.

k 3 На первом шаге найдем нормальное псевдорешение только первой сис темы Ax1 = z1. Оно равно x0 = A+z1. Проекция вектора z1 на R(A1):

z1 = A1x0 = (A1A+)z1.

4 Рассмотрим систему A1 z x=, (10.13) A2 z которая отличается от исходной системы A1 z x= (10.14) A2 z тем, что вместо z1 используется z1. Для нее нормальное псевдорешение обозначим x0. Оно равно + A1 z x0 =.

A2 z Вопрос: верно ли равенство x0 = x0? Напомним, что + A1 z x0 =.

A2 z 10 Теоретические основы Решение вопроса:

Для системы (10.13) запишем нормальные уравнения:

A1 z AT AT AT AT x=, 1 2 1 A2 z AT A1 + AT A2 x = AT z1 + AT z2.

1 (10.15) 1 2 Затем для системы (10.14) запишем нормальные уравнения:

A1 z AT AT AT AT x=, 1 2 1 A2 z AT A1 + AT A2 x = AT z1 + AT z2. (10.16) 1 2 1 Сравним в правой части AT z1 и AT z1. Имеем 1 AT z1 = AT (A1A+ )z1.

1 1 Более простое доказательство:

AT (AA+) = AT (AA+)T = (AA+A)T = AT.

Лемма 10.2. Докажем, что AT = AT (AA+).

Доказательство. Имеем (см. стр. 220):

AT = U T L T, A+ = U + L +, AA+ = LU U + L+, A = LU, AT (AA+) = U T LT · LU · U +L+ = U T LT = AT.

Лемма доказана. Поэтому AT z1 = AT z1.

1 Докажем то же самое другим способом:

AT (1) = AT (z z ) = AT z AT z, где z N (AT ), т. е. AT z = 0.

1z 1 1 1 Поэтому AT z1 = AT z1. Таким образом, правые части уравнений (10.15) 1 и (10.16) совпадают. Поэтому решения уравнений (10.15) и (10.16) оди наковы, то есть x0 = x0.

10.9 Основные свойства симметрических / эрмитовых матриц 10.9 Основные свойства симметрических / эрмитовых матриц Свойства симметрических матриц, A = AT, в теории МНК служат осно вой многих результатов. В алгебре эти свойства обобщены на случай эрми товых матриц, т. е. матриц A, в которых понятие симметричности в отно шении комплекснозначных элементов расширено добавлением требования комплексной сопряженности aij = aji, что выражается записью A = A, где A сопряженно транспонированная матрица к A.

Напомним некоторые определения, а затем приведем четыре основные свойства эрмитовых матриц.

Определение 10.9. Число называется собственным значением (nn)-матрицы A с соответствующим собственным вектором x, если Ax = = x, т. е. есть любой из n корней уравнения det(A I) = 0.

характеристическое уравнение для матрицы A.

Это Определение 10.10. Ортогональной матрицей называется квадрат ная вещественная матрица Q, столбцы которой ортонормированы, то есть выполнено условие QT Q = I.

Упражнение 10.1. Матрица QT ортогональна тогда и только тогда, когда ортогональна матрица Q.

Определение 10.11. Унитарной матрицей называется квадратная матрица U, столбцы которой ортонормированы с добавлением требования комплексной сопряженности элементов столбцов: U U = I.

Упражнение 10.2. Матрица U унитарна тогда и только тогда, когда унитарна матрица U.

Теорема 10.17. Если A = A, то для всех комплексных векторов x число xAx вещественно.

Доказательство. Вычислите (xAx) как эрмитову (1 1)-матрицу и примите во внимание, что в любой эрмитовой матрице диагональные эле менты должны быть вещественными.

Теорема 10.18. Каждое собственное значение эрмитовой матрицы вещественно.

10 Теоретические основы Доказательство. Умножьте Ax = x слева на x и используйте преды дущую теорему 10.17.

Теорема 10.19 ( [13]). Собственные векторы эрмитовой матрицы, отвечающие различным собственным значениям, взаимно ортогональны.

Доказательство. Пусть = µ и Ax = x, Ay = µy. Возьмем (1 n) матрицу xA = x, сопряженно транспонированную к Ax = x, заменим A = A и учтем, что = (по теореме 10.18), т. е. запишем x A = x.

Умножая это равенство справа на y, а равенство Ay = µy слева на x, полу чим x Ay = x y, xAy = µxy.

Следовательно, xy = µxy, а так как = µ, то xy = 0, т. е. вектор x ортогонален к y. Этот результат является основным по своей важности. Разумеется, любые кратные собственным векторам x/ и y/, равным образом оста ются собственными векторами. При выборе = x, = y и т. п.

мы имеем возможность нормировать все собственные векторы к единичной длине и, таким образом, далее использовать ортонормированные собствен ные векторы {u1, u2,..., un} эрмитовой матрицы A. Запишем их в виде мат рицы U = [ u1 | u2 |... | un ], которая, по определению, является унитарной:

U U = I. Имеем Aui = i ui, по определению собственных значений i, i = 1, 2,..., n. Эту систему уравнений перепишем в виде одного уравнения, если введем обозначение = diag {1, 2,..., n }:

AU = U, что равносильно соотношениям U AU =, A = U U.

Тем самым приходим к следующей теореме о спектральном разложении эрмитовой матрицы.

Теорема 10.20. Если A = A, то существует диагонализирующая матрица T, которая является также и унитарной, T T = I, такая что T AT =, A = T T, (10.17) где = diag {1, 2,..., n} диагональная матрица, составленная из соб ственных значений матрицы A. Если эти значения простые (среди i нет 10.9 Основные свойства симметрических / эрмитовых матриц кратных), то T = U, где U = [ u1 | u2 |... | un ] матрица, составленная из соответствующих собственных векторов матрицы A. В этом случае спек тральное разложение матрицы A имеет вид A = U U = 1 u1u + 2 u2 u +... + n un u.

1 2 n При наличии кратных собственных значений любая эрмитова матрица A по-прежнему имеет полный набор ортонормированных собственных векто ров и, следовательно, может быть диагонализирована с помощью некоторой унитарной матрицы T, т. е. (10.17) справедливо в общем случае. Каждая эрмитова матрица с k различными собственными значениями имеет свое спектральное разложение A = T T = 1 P1 + 2 P2 +... + k Pk, где Pi есть проекция матрицы A на собственное подпространство, соответ ствующее значению i.

Доказательство. Полное доказательство можно найти в [13] или [23].

Это очень важный результат, имеющий и обратную формулировку:

только эрмитовы матрицы обладают одновременно и вещественными соб ственными значениями, и ортонормированными собственными векторами.

Если T 1AT равняется некоторой вещественной диагональной матрице и матрица T унитарна, T 1 = T, то матрица A обязательно является эрмито вой: A = (T T ) = T T = A. Теорема 10.20 имеет множество примене ний при обработке экспериментальных данных, не только в регрессион ном анализе, но также в области так называемого многомерного анализа.

Интересны примеры, когда по экспериментальной корреляционной матрице R нужно провести анализ основных компонент и факторный анализ [13]. Для специалистов-прикладников эти примеры очень поучительны. Они показы вают, как два специалиста, вводя различное число p факторов и диагональ ную компоненту D в разложении R = F F T + D, а также любую ортогональ ную (p p)-матрицу Q в представлении F = F Q, могут получать совер шенно различные интерпретации F и F (матрицы факторных коэффициен тов) одних и тех же экспериментальных данных. Мы эти частные вопросы не рассматриваем и за деталями отсылаем к [13] и специальной литературе по факторному анализу. Когда A = AT, (10.17) принимает следующий вид:

T T AT =, A = T T T.

Оценивание по методу наименьших квадратов 11.1 Модели, регрессии и оценки Изучение объектов, явлений или процессов реального мира заключается в построении их моделей, т. е. тех или иных форм описания для выявле ния существенных закономерностей. Когда тип, т. е. форма модели выбрана, требуется определить ее наилучшее параметрическое наполнение. При этом критерием естественно считать соответствие между откликами объекта и модели, когда они погружены в одни и те же или одинаковые внешние усло вия. Иными словами, та модель будет наилучшей, чей отклик, по данным проведенного наблюдения, менее всего уклоняется от отклика реального объекта в разнообразных, но одинаковых для модели и объекта условиях.

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

Регрессионный анализ в статистике, регрессионное моделирование, иден тификация моделей и оценивание состояния объекта в теории систем, адап тивная обработка сигналов это примеры очень близких или тесно связан ных задач, имеющих общую математическую основу, один и тот же матема тический метод решения [21, 26, 27, 32, 33, 34, 50, 53, 62, 64, 66, 69, 82, 88].

Если модель линейная по параметрам, а критерий ее качества есть квад рат расхождения (невязки) между откликами модели и объекта, этот метод и есть знаменитый линейный Метод Наименьших Квадратов (МНК) [47, 53], возникновение которого связывается с работами Гаусса и Лежандра в начале 19 столетия, см. подробнее стр. 199. При отказе от линейности модели по ее 11.1 Модели, регрессии и оценки параметрам получаем нелинейную задачу о наименьших квадратах [22, 35].

В разделах 11 и 12 этой книги мы ограничиваемся обыкновенной линей ной статической детерминистской задачей НК [15]. Она хорошо разработана, но имеет обобщения по разным направлениям (для них она является базо вой). Основные обобщения:

1. Динамическая задача: оцениваемый вектор не постоянен, а эволюци онирует во времени, подчиняясь какому-либо уравнению состояния [61, 115]. Это обобщение присутствует в разд. 13, 14.

2. Стохастическая задача: оцениваемый вектор постоянен, но случаен, и он измеряется на фоне случайных помех [61, 115]. Это обобщение начина ется в подразд. 11.9. Если задача, к тому же, динамическая, то допуска ются случайные возмущения в уравнении состояния или/и в уравнении наблюдений [61, 115]. Это обобщение присутствует в разд. 13, 14.

3. Полная МНК-задача: имеются возмущения в матрице наблюдений, т. е.

в матрице регрессоров. Это обобщение здесь отсутствует. См. Total Least Squares в [18, 38, 100].

4. Неоднородность погрешностей в исходных данных [49]. Это обобщение присутствует в разд. 13, 14.

5. Зависимость ковариации погрешностей от неизвестного параметра (Generalized Least Squares) [106]. Это обобщение здесь отсутствует. Воз можен переход к адаптивному оцениванию [73].

6. Зависимость от неизвестного параметра не только ковариации погреш ностей, но и других элементов стохастической модели источника дан ных. Это обобщение здесь отсутствует. Возможен переход к адаптив ному оцениванию [73].

7. Улучшение вычислительных схем [15, 89, 93, 95, 103, 104, 107, 108, 119, 123, 125, 135, 129]. Это обобщение присутствует в разд. 13, 14 и состав ляет их основное содержание.

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

11.2 Линейная задача наименьших квадратов Во многих приложениях, связанных с обработкой экспериментальных данных, необходимо отыскивать такой вектор x Rn, линейные комбина матрица размера (m n), ции компонент которого, Ax, где A = A(m, n) как можно более близки или, еще лучше, равны данным значениям, обра зующим вектор z Rm, т. е. Ax z. Если мерой близости двух векторов считать квадрат евклидовой нормы разностного вектора, в данном случае, вектора v zAx, то указанная задача есть линейная задача о наименьших квадратах (см. также подразд. 10.4).

Возможность сделать равным нулю вектор невязок v = zAx существует тогда и только тогда, когда z R(A), где R(A) пространство столбцов матрицы A. В этом случае имеем совместную систему уравнений Ax = z.

Однако, z вектор наблюдений, то есть экспериментальных данных и A матрица, которую задают до того, как получат z и которую в различных приложениях называют либо матрицей регрессоров, либо матрицей наблю дений, либо матрицей плана эксперимента. Совсем не обязательно, что усло вие z R(A) будет выполнено, например, из-за случайных погрешностей v во время регистрации экспериментальных данных. Тогда z = Ax + v, (11.1) и решение по методу наименьших квадратов (для краткости, МНК решение) есть вектор x, доставляющий минимум функционалу качества (см.

рис. 11.1):

m T v(j)2 = v T v.

J(x) = (z Ax) (z Ax) = (11.2) j= Требуя минимум этого критерия, для искомого x получаем так называе мые нормальные уравнения (см. подразд. 10.4, стр. 209):

AT A = AT z.

x (11.3) 11.2 Линейная задача наименьших квадратов вход A выход z Объект + Модель v = z Ax Ax x min v x Рис. 11.1. Линейная задача наименьших квадратов Их решение всегда существует (обе части равенства (11.3) принадлежат од ному и тому же пространству R(AT ) столбцов матрицы AT ), но может быть не единственным (если rank A n). В последнем случае из всех x выбирают то единственное, x0, которое имеет минимальную норму x0. Этот вектор называют нормальным псевдорешением. Известно (см. подразд. 10.4), что x0 = A+z, (11.4) где A+ псевдообратная матрица к A. Как уже отмечалось, в каче стве определения A+ применяют различные формулировки. Здесь для этого используем геометрический подход: A+ есть такая матрица в выраже нии (11.4), что для любого z Rm вектор x0 Rn удовлетворяет двум условиям:

A0 = p, p R(A), z p R(A).

x (11.5) x0 R(AT ).

(11.6) Условие (11.5) требует, чтобы x0 отвечал совместной системе A0 = p, где x проекция вектора z на R(A), а условие (11.6) требует, чтобы этот x p T был взят из пространства R(A ) строк матрицы A. Условие (11.5), таким образом, выбирает x0 = x, чтобы минимизировать функционал (11.2), а условие (11.6) среди всех таких x выбирает единственный x0 с минимальной нормой.

Часто матрицу A выбирают так, чтобы она имела полный столбцовый ранг, rank A = n. В этом случае m n, x единственно и равно x0, A+ = 1 T T = (A A) A и x0 = (AT A)1AT z.

(11.8) 11 Оценивание по методу наименьших квадратов Однако иногда такое условие не выполняется, и тогда x0 = A+ z, где A+ псевдообратная матрица (см. подразд. 10.4–10.5, стр. 209–212).

11.3 Статистическая интерпретация Предположим, что вектор ошибок v в уравнении (11.1) образован из слу чайных величин с нулевым средним и известной матрицей ковариации E vv T = Pv, E {v} = 0, (11.9) где E {·} оператор математического ожидания (среднего) над ·, и Pv ПО (положительно определенная) матрица. Найдем квадратно-корневое разло жение Pv = SS T (например, разложение Холесского). Если теперь умножить вектор z (11.1) на S 1, то данные z = S 1z получают представление z = Ax + v (11.10) с матрицей A = S 1A и ошибками v = S 1v. Этот вектор ошибок всегда имеет единичную ковариацию:

E v v T = E S 1vv T S T = S 1 E vv T S T = S 1SS T S T = Im, единичная матрица размера (m m). Вследствие этого данные где Im z называют нормализованными экспериментальными данными. Значение представления (11.10) заключается в том, что оно демонстрирует, как скон струировать вектор некоррелированных между собой измерений с единич ной дисперсией из вектора, элементы которого произвольно взаимно корре лированы (декоррелировать и нормализовать его). Ниже предполагаем, что данные z (11.1) уже декоррелированы и нормализованы, так что E vv T = Im, E {v} = 0, (11.11) единичная матрица размера (m m). При этом из (11.3) находим где Im AT A = AT z = AT Ax + AT v, x AT A( x) = AT v.

x Отсюда, если det(AT A) = 0, имеем E {} = x, x (11.12) (AT A) E ( x)( x)T (AT A) = AT E vv T A = AT A.

x x (11.13) 11.4 Включение априорных статистических данных Соотношение (11.12) выражает собой свойство несмещенности решения (оценки) x относительно неизвестного (постоянного) вектора x, измеряемого в виде экспериментальных данных z, (11.1) или (11.10). Соотношение (11.13) дает выражение для ковариации оценки x в виде Px = E ( x)( x)T = (AT A)1.

x x (11.14) при определении x по нормализованным экспериментальным данным.

Обратная матрица Px от ковариации Px называется информационной матрицей. Ее обозначение будет x или просто. При использовании нор мализованных данных она равна AT A, а в более общем случае (11.9) она равна = AT Pv A.

11.4 Включение априорных статистических данных Предположим, что в добавление к линейной системе (11.1) мы имеем априорную несмещенную оценку неизвестного вектора x в виде x и соот Это означает, что ветствующую априорную информационную матрицу.

E {} = x и x 1 = E ( x)( x)T = P, x x (11.15) где P ковариация оценки x. Найдем какой-нибудь квадратный корень 1/ из матрицы, например, по одному из разложений Холесского (см. разд. 6):

= 1/2T /2 = RT R, x где 1/2 = RT. Образуем вектор v = (1/2)T ( x) = R( x). Этот вектор x имеет смысл нормализованной ошибки для априорной оценки x вектора x.

Действительно, его ковариация равна единичной матрице размера (n n):

E v v T = R E ( x)( x)T RT = T /21 1/2 = In.

x x Так как о векторе x, кроме экспериментальных данных z, (11.1), известна априорная оценка x с ковариацией P = 1, эту информацию целесообразно включить в контекст задачи о наименьших квадратах, рассматривая моди фицированный функционал качества J1(x) = v T v + v T v вместо (11.2). Он соединяет в себе квадрат нормы нормализованной ошибки (невязки) апри орной оценки x v = R( x) = T /2( x), x 11 Оценивание по методу наименьших квадратов + v = z Ax A z Объект Ax Модель 2 + v 2) min( v x x Rx Априорная модель x + x v = z Rx z = R R Рис. 11.2. Включение априорных данных в линейную задачу НК с квадратом нормы нормализованной ошибки (невязки) экспериментальных данных v = z Ax. Так как x J1(x) = ( x)T ( x) + (z Ax)T (z Ax) = x (11.16) z z = ( Rx)T ( Rx) + (z Ax)T (z Ax), x где z = R, то J1(x) может быть интерпретирован просто как критерий каче ства метода наименьших квадратов применительно к расширенной системе (рис. 11.2) z v R = x+, (11.17) z v A включающей, помимо текущих экспериментальных данных z, дополнитель ные экспериментальные данные z, соответствующие имеющейся в наличии априорной информации x,.

Обозначим через x МНК-решение расширенной системы (11.17), достав ляющее минимум функционалу (11.16). Из этого критерия для x получаем, аналогично (11.3), нормальные уравнения x x ( + AT A) = + AT z. (11.18) Простая модификация данной в п. 11.3 статистической интерпретации при водит к выражению x x ( x) = ( x) + AT v. (11.19) 11.5 Включение предшествующего МНК-решения Так как E { x} = 0 и E {v} = 0, то и E { x} = 0, то есть x есть также x x несмещенная оценка: E {} = x. Если ошибка априорной оценки и ошибка x измерения взаимно некоррелированы, E ( x)v T = 0, то после возве x дения в квадрат обеих частей (11.19) и осреднения получим ковариацию Px = E ( x)( x)T = ( + AT A)1 = 1, x x (11.20) где через обозначена информационная матрица апостериорной оценки x, = + AT A.

Замечание 11.1. Матрица не обязана быть невырожденной, хотя в (11.15) это формально предполагалось. В действительности, априорное знание некоторых (или всех) компонент вектора x может быть исчезающе мало, так что соответствующие строки и столбцы в информационной мат рице заполняются исчезающе малыми числами или даже нулями. При этом соотношение (11.15) сохраняет силу в пределе, в том смысле, что в кова риационной матрице P соответствующие диагональные элементы стремятся к +, в то время как другие остаются ограниченными. Заметьте, что (11.18) как условие минимума функционала (11.16) получается при произвольной неотрицательно определенной матрице. Если = 0, то (11.18) сводится к (11.3) и (11.20) сводится к (11.14). Таким образом, информационная и ко вариационная матрицы взаимно обратны не только формально в силу опре деления (11.15), но и по существу как меры достоверности/недостоверности априорной оценки x вектора x.

11.5 Включение предшествующего МНК-решения Расширенная система (11.17) показала, что априорные статистические сведения о векторе x, поступающие в виде несмещенной оценки x и ее кова, могут быть интерпретированы как воображаемые добавочные риации P x результаты z = R некоего эксперимента. Это наводит на мысль, что и чисто алгебраическую задачу отыскания МНК-решения системы уравнений можно решать последовательно: предварительно найти x как МНК-решение части системы, а затем включить это x в полную систему, чтобы найти ее МНК-решение x. Такое разделение системы на две части априорную и можно назвать ее расщеплением (см. формальное решение текущую рассматриваемого здесь вопроса о расщеплении в подразд. 10.8).

11 Оценивание по методу наименьших квадратов Пусть система уравнений произвольно расщеплена на эти две подсистемы:

f w R = x+. (11.21) z v A МНК-решение x этой полной системы, доставляющее минимум функцио налу J1 (x) = w w + v T v, есть решение нормальных уранений T f R R T AT R T AT x=. (11.22) z A Допустим, найдено МНК-решение x для подсистемы f = Rx + w из крите рия минимума функционала J(x) = wT w. Как отмечено в (11.5)–(11.6), оно удовлетворяет двум условиям:

x z R(R), f z R(R), R = z, x R(RT ).

Разностный вектор r = f z ортогонален пространству столбцов R(R) мат рицы R и, следовательно, лежит в левом нуль-пространстве N (RT ), опреде ляемом как все векторы y, удовлетворяющие уравнению RT y = 0. Поэтому z RT f = RT ( + r) = RT z + RT r = RT z.

Следовательно, уравнения (11.22) совпадают с уравнениями z R R T AT R T AT x=, z A которые, в свою очередь совпадают с уранениями (11.18), так как RT R =.

Тем самым доказано, что МНК-решение x данной системы (11.21) совпадает с МНК-решением системы (11.17), отличающейся от (11.21) тем, что в нее x вместо f включен вектор z, равный проекции f на R(R), z = R, где x МНК-решение левой подсистемы в (11.21).

11.6 Рекурсия МНК в стандартной информационной форме Интерпретация априорных статистических данных как дополнительных наблюдений, или, что равносильно, учет имеющегося МНК-решения подси стемы после добавления в систему новой порции уравнений, является крае угольным камнем рекурсии для МНК. Это дает возможность обрабатывать 11.6 Рекурсия МНК в стандартной информационной форме экспериментальные данные по мере их поступления, то есть решать задачу о наименьших квадратах по мере поступления новых уравнений. Это рав носильно также тому, что исходную большую совокупность эксперименталь ных данных можно расщеплять произвольно на порции и последовательно включать их в обработку. Результат такой последовательной обработки, как доказано выше (подразд. 10.8 и 11.5), будет равносилен результату обра ботки всей совокупности экспериментальных данных целиком.

Результаты при статистической интерпретации рекурсивны потому, что текущие величины оценка x и ковариация Px становятся априорными и комбинируются с новыми данными, чтобы образовать обновленные оценку и ковариацию. При этом существенно, что результаты (11.18) и (11.20) не зависят (теоретически) от того, какой квадратный корень R в разложении = RT R использован. Эта свобода позволяет выбирать R из соображений большей вычислительной точности. Кроме того, если лишь окончательная оценка (МНК-решение) необходима, то лучше не находить промежуточных AT Aj и сумму оценок, а просто накапливать информационную матрицу j T Aj zj, и лишь в нужный момент (например, в самом конце) вычислить решение.

Информационную форму последовательного МНК запишем, вводя мат рицу ( | d).

I. Инициализация. Устанавливают начальные значения x0, 0 :

d0 = 0 x0, d := 0 d0.

Эти начальные значения берут из априорных данных: x0 = x, 0 =.

II. Обработка наблюдений. Вводят очередную порцию наблюдений z = = Ax + v:

d := d + AT A z. (11.23) В общем случае ненормализованных статистических данных z вме сто (11.23) используют алгоритм:

+ AT R d := d Az. (11.24) III. Выдача результата. После последовательной обработки всех порций наблюдений или в нужный момент, когда 1 существует, вычисляют x = 1d, Px = 1.

11 Оценивание по методу наименьших квадратов 11.7 Рекурсия МНК в стандартной ковариационной форме Пусть априорная и апостериорная оценки, x и x, характеризуются невы рожденными информационными матрицами и, соответственно. Тогда существуют обратные к ним ковариационные матрицы, (11.15) и (11.20).

Разрешим (11.18) относительно x:

x x = ( + AT A)1 + ( + AT A)1AT z.

Обозначим L = ( + AT A)1, K = ( + AT A)1AT (11.25) и преобразуем:

x = L + Kz KA + KA = x + K(z A), x x x x так как L+KA = In. Для определения K воспользуемся в (11.25) следующей важной леммой.

Лемма 11.1.

(1 121 21)1 = 1 + 112 (2 21 112 )121 1, (11.26) 2 1 1 1 где предполагается, что все матрицы согласованы по размерам и тре буемые обращения матриц существуют.

Упражнение 11.1. Докажите лемму 11.1, рассматривая блочные мат рицы 1 12 P1 P =, 21 2 P21 P и выписывая поблочно равенства P = I и P = I.

Применим лемму 11.1 при 1 =, 12 = AT, 21 = T = A, 1 = I.

12 Получим ( + AT A)1 = 1 1AT (A1AT + I)1A1.

Обозначим: P = 1, P = 1. Имеем из подразд. 11.5 выражение (11.23):

= + AT A. Следовательно, P = P P AT (AP AT + I)1AP. Так как K = P AT, то K = P AT P AT (AP AT + I)1AP AT = = P AT (AP AT + I)1[AP AT + I AP AT ] = P AT (AP AT + I)1.

11.7 Рекурсия МНК в стандартной ковариационной форме Таким образом, при статистической интерпретации (см. подразд. 11.4) по лучаем возможность уточнять априорную несмещенную оценку x и умень шать ее ковариацию P, благодаря включению в процесс обработки вновь поступивших результатов наблюдений z = Ax +v и применяя к ним следую щий алгоритм, известный как стандартный алгоритм Калмана (рис. 11.3).

A z Объект + Априорная z = A x Обновление модель модели x, P Экстраполяция K модели + x := x, P := P Апостериорная модель x, P Рис. 11.3. Рекурсия МНК в стандартной ковариационной форме (схема Калмана).

Этап 1 обновление модели по наблюдениям. Этап 2 экстраполяция модели меж ду наблюдениями. Априорная модель дает предсказание z для отклика z объекта.

Разность имеет смысл обновляющего процесса. Весовая матрица K, умноженная на, обеспечивает обновление априорной модели по наблюдению z, т. е. переход к апостериорной модели I. Инициализация. Начальные значения x0, P0 :

x := x0, P := P0. (11.27) II. Обработка наблюдений. Очередная порция наблюдений z = Ax + v :

K = P AT (AP AT + I)1, (11.28) P = P KAP, x = x + K(z A).

x 11 Оценивание по методу наименьших квадратов III. Экстраполяция. Распространение оценки x и ее ковариации P между наблюдениями, т. е. к моменту повторения этапа II со следующей пор цией наблюдений :

P := P, x := x.

(11.29) Замечание 11.2. Первая и вторая формулы в (11.29) выражают пра вило экстраполяции только для статической МНК-задачи оценивания, т. е.

для случая, когда оцениваемый вектор x не изменяется в процессе наблю дения: x = const. Они заменяются другими в случае динамической МНК задачи оценивания первыми двумя формулами в системе уравнений (13.3) (см. ниже стр. 273).

Равным образом данный алгоритм пригоден и без статистической интер претации (см. подразд. 11.5), когда алгебраическая задача отыскания МНК решения переопределенной системы решается последовательно. Такое реше ние может стартовать с условно пустой системы уравнений. Практически, это должно отвечать условию = 0, которое легко реализовать в инфор мационной форме (подразд. 11.5). В ковариационной форме данное условие можно реализовать лишь приближенно, например, полагая P0 = 2I, где 0. При таком выборе величина x0 практически не имеет влияния на дальнейший процесс, так что она может быть взята равной нулевому век тору, x0 = 0. После такой инициализации уравнения исходной переопреде ленной системы могут вводиться в этап обработки измерений последователь ными порциями, и в порциях может содержаться любое, не обязательно одно и то же, число уравнений, выражаемое в числе строк матрицы A.

Как следует из подразд. 11.5, от указанного числа МНК-решение всей алгебраической системы уравнений не зависит. В связи с этим, с вычисли тельной точки зрения удобным оказывается добавление в очередной пор ции лишь по одному уравнению. Тогда матрица A содержит всего одну строку, которую теперь обозначим как транспонированный вектор-столбец a, aT = (a1, a2,..., an ). В этом случае z = aT x + v, (11.30) и обработка наблюдений (11.28) принимает особенно простой вид:

= aT P a+1, P = P KaT P, x = x +K(z aT x). (11.31) K = P a/, Это алгоритм скалярной (последовательной) обработки. В нем умножение на обратную матрицу (AP AT + I)1 заменяется делением на скалярную величину.

11.7 Рекурсия МНК в стандартной ковариационной форме Упражнение 11.2. С применением леммы 11.1 к выражению (11.24) докажите, что в общем случае ненормализованных экспериментальных дан ных (11.1) при их статистической интерпретации с ковариацией ошибок v, равной R, матрица Калмана K в алгоритме (11.28) определяется выраже нием K = P AT (AP AT + R)1. (11.32) Упражнение 11.3. Статистическая интерпретация алгоритма (11.31) дана в подразд. 11.3 для случая, когда экспериментальные данные (11.1) нормализованы, то есть характеризуются математическими ожидани ями (11.11). Это выражается добавлением +1 в выражении для, (11.31), причем эта +1 есть не что иное как E vv T = 1 для ошибки v в (11.30).

Докажите, что в более общем случае, когда матрица R в (11.32) является диагональной, R = diag (r1, r2,..., rm ), (11.33) и только в этом случае матричный алгоритм (11.28) с матрицей K из (11.32) эквивалентен m-кратному повторению скалярного алгоритма вида (11.31) с = aT P a + r, где aT i-я строка матрицы A, r = ri и z i-й элемент вектора z, (11.1), при i = 1, 2,..., m.

Таким образом, скаляризованный алгоритм Калмана (11.31) имеет сле дующее общее представление:

= aT P a+r, P = P KaT P, x = x +K(z aT x), (11.34) K = P a/, если наблюдения (11.1) в их статистической интерпретации составлены из m отдельных, независимых друг от друга, скалярных данных вида (11.30), каждое с ковариацией r = ri, i = 1, 2,..., m.

Еще раз отметим, что в применении к решению переопределенной сис темы алгебраических уравнений в (11.34) следует считать r = 1, то есть использовать (11.31).

Замечание 11.3. Условие (11.33) не является ни в коей мере ограничительным для использования (11.34). Используя разложение Холес ского без квадратных корней (см. разд. 6), любую R 0 можно предста вить в виде R = U DU T или R = LDLT и затем перейти к измерениям z = U 1 z или z = L1z, чтобы диагонализировать матрицу ковариаций ошибок наблюдений.

11 Оценивание по методу наименьших квадратов Упражнение 11.4. Докажите, что в алгоритме (11.28) с определением K по выражению (11.32) вычисление P может быть представлено в так назы ваемой симметричной форме Джозефа:

P = (I KA)P (I KA)T + KRK T, которую иначе называют стабилизированным алгоритмом Калмана, так как она предотвращает возможную потерю положительной определенности матрицы P, присущую стандартному алгоритму (11.28) с P = P KAP.

При скалярной обработке наблюдений в алгоритме (11.34) это выражение для P, соответственно, заменяется на скаляризованный алгоритм Джозефа P = (I KT )P (I K T ) + rKK T.

11.8 Ковариационный алгоритм Поттера для МНК Вместо матриц P и P, по своей природе положительно определенных, далее оперируем с квадратными корнями из них, соответственно, S и S, отвечающими равенствам S S T и S S T. Перепишем выражение для P в (11.34) в виде S S T = S(In f f T /)S T, f = S T a, = r + f T f, где n размерность векторов x, x, и потребуем так выбрать число, чтобы обеспечить справедливость следующего разложения:

In f f T / = (In f f T )(In f f T ).

Отсюда для получаем квадратное уравнение и из двух его решений выби раем = (1/)/(1 + r/), поскольку выбор знака ” + ” обеспечивает меньший уровень относительных ошибок при этих вычислениях. Обозначим = 1/(1 + r/), тогда = /. В итоге вместо (11.34) получаем следующий ряд формул:

f = ST a, T = f f +r, = 1/(1 + r/), (11.35) K = Sf /, = S Kf, T S T x = x + K(z a x), 11.9 Полная статистическая интерпретация рекурсии в МНК который и составляет алгоритм Поттера. Он численно более устойчив, чем стандартный ковариационный алгоритм Калмана (11.34), но ему эквивален тен. В целом, для него характерно следующее.

Вычисление S в (11.35) равносильно счету с двойной точностью матрицы P в (11.34) при использовании той же разрядности чисел в компьютере, или, иначе, равносильная точность вычислений матрицы P может быть достигнута значительно быстрее. Для матрицы P теперь отсутствует опас ность потери положительной определенности, присущая операции вычита ния в (11.34), поскольку здесь вычисляют S, а P = S S T и P 0 при det(S) = 0. Недостатком алгоритма (11.35) является наличие операции извлечения квадратного корня, отдельной для каждого скалярного наблю дения z = aT x+v, и возможная потеря специального (желательно, треуголь ного) вида матрицы S в общем случае. Действительно, для экономии памяти и объема вычислений обычно стремятся иметь матрицы S и S в виде тре угольных матриц (обе нижние треугольные или обе верхние треуголь ные), что соответствует разложениям Холесского: P = S S T и P = S S T.

Однако, если стартовать с матрицы S треугольного вида, выполняя иници ализацию в соответствии с (11.27), то из-за вычитания в (11.35) матрица S в общем случае не остается треугольной. Например, пусть для S и S выбрана верхняя треугольная форма. Тогда только при a = (1, 0,..., 0)T, где скаляр, для S в (11.35) будет сохранена та же верхняя треугольная форма, благодаря чему выполнение этапа экстраполяции, согласно (11.29), сводится к простому присваиванию: S := S. Если же выбранная для S треугольная форма будет утрачена, то этап экстраполяции матрицы потребует предва рительной триангуляризации матрицы S, то есть операции S := triang (S).

Триангуляризация triang (·) должна быть выполнена ортогональным преоб разованием матрицы (·), например, преобразованиями Хаусхолдера, Гивенса или же Грама–Шмидта, которые рассматривались выше (разд. 7).

11.9 Полная статистическая интерпретация рекурсии в МНК Простая статистическая интерпретация МНК-решения, данная в п. 11.3, позволила объяснить в п. 11.4 идею включения априорных данных в про цесс решения задачи о наименьших квадратах. Она была простой в том смысле, что предполагала случайную природу лишь для ошибки v в экспери ментальных данных z, (11.1), при необходимости находить апостериорную 11 Оценивание по методу наименьших квадратов v x z Ax A + Рис. 11.4. Представление экспериментальных данных z в виде результата измерения неизвестного вектора x с помощью матрицы A в присутствии случайных ошибок v оценку x для некоторого (постоянного) неизвестного вектора x с учетом некоторой имеющейся (априорной) оценки x (см. рис. 11.4). Для той интер претации было достаточно использовать первые два момента случайной ошибки v: математическое ожидание E {v} = 0 и ковариацию E vv T = Pv.

Ниже, описывая случайные векторы, мы также ограничиваемся первыми двумя моментами распределения вероятностей, то есть фактически прини маем гипотезу о нормальном (гауссовом) распределении. Отличие будет в том, что с этого момента мы вводим статистическое описание не только для v, но и для оцениваемого вектора x. При этом надо иметь в виду, что такое описание для x должно всегда рассматриваться как условное распределе ние. Именно это обстоятельство делает приводимую ниже статистическую интерпретацию МНК полной, т. е. учитывающей рекурсию процесса обнов ления оценок в два этапа.

1. Экстраполяция оценок по времени между моментами измерений;

озна чает обновление не по измерению, а по времени;

его иногда называют темпоральным обновлением ;

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

2. Обновление оценок по измерениям;

означает переход от текущих апри орных оценок (текущего априорного условного распределения) к сле дующим апостериорным оценкам (новому апостериорному условному распределению вероятностей) для оцениваемого вектора x.

Предположим, что ошибка наблюдения v в уравнении (11.1) есть случай ный гауссов вектор. Мы предполагаем, что это уравнение нормализовано, т. е. записано по типу уравнения (11.10), поэтому E vv T = Im, E {v} = 0, (11.36) где E {·} оператор математического ожидания, 0 нулевой m-вектор, единичная матрица (m m)-матрица.


Im 11.9 Полная статистическая интерпретация рекурсии в МНК Замечание 11.4. Предположение о предварительной нормализации не нарушает общности и введено только для упрощения записей. В любой момент оно может быть отозвано путем умножения нормализованных дан ных на корень квадратный L из матрицы R, где LLT = R (см. п. 11.3) и R положительно определенная матрица (R 0). Последнее, т. е. R 0, означает, что не найдется такого невырожденного преобразования вектора наблюдений z, (11.1), после которого результат содержал бы элементы, сво бодные от случайных ошибок. Иными словами, все экспериментальные дан ные содержат случайные погрешности. Случай, когда некоторые измерения являются точными, может быть рассмотрен особо [115].

Таким образом, плотность распределения вероятностей вектора v опре делена выражением 1/ exp T R1, fv () = (2)m|R | (11.37) где согласно (11.36), R = E vv T = Im.

Допустим вместе с этим, что еще до проведения эксперимента по схеме z = Ax + v (11.38) и независимо от него вы располагаете априорной информацией о неизвест ном векторе x. Это предварительное знание, явившееся результатом ваших предыдущих экспериментов либо просто кем-то другим полученное и вам сообщенное, выражается величиной x. Сама по себе величина x тоже может быть не вполне надежной, то есть вам следует ее рассматривать как слу чайную, а то конкретное ее значение, которое вам сообщено, вы обознача ете. Это имеет смысл реализованного значения случайного вектора x.

которое вы принимаете за центр распределения вектора x, Кроме этого, вам нужно иметь сведения о разбросе величины x относительно центра, т. е.

. Имея и P и продолжая работать лишь с нормаль нужна ковариация P ными распределениями, вы выражаете эту априорную информацию следу ющей условной (априорно принятой) плотностью распределения вероятно стей вектора x:

1/ exp ( )T P 1 ( ).

fx| (|) = (2)n|P | (11.39) x Здесь вы формально предполагаете, что P обратима (т. е. P 0), чтобы иметь возможность записать это выражение.

11 Оценивание по методу наименьших квадратов Теперь вы проводите свое наблюдение по схеме (11.38) и получаете кон кретный его результат, реализованное значение случайного вектора z.

Вы хотите использовать результат, чтобы выработать наилучшую оценку x для вектора x. Как для этого скомбинировать априорное знание (, P ) с полученным результатом ? Ответ на этот вопрос кроется в апостериорном распределении, т. е. в плотности fx|z, (, ) распределения вероятностей x неизвестного вектора x при условии, что известны: (1) результат для век тора наблюдений z и (2) сообщенное значение априорного данного вектора x с ковариацией P, что отражено принятием априорной плотности (11.39).

Если удастся найти эту fx|z, (, ), то по ней можно будет принять ре x шение, поскольку именно плотность распределения для любой случайной величины служит наиболее полной характеристикой.

Найдем fx|z, (, ). По формуле Байеса для условных плотностей x fx,z, (,, ) x fx|z, (, ) =. (11.40) x fz,(, ) x Преобразуем (11.40) по известным из теории вероятностей законам, опус кая, для простоты промежуточных записей, аргументы. Получаем x fz|x, · fx| · fx fz|x, (, )fx| ( ) x x x fx|z, (, ) = =. (11.41) x fz| · fx fz| ( ) x x Здесь плотность fx| ( ) известна как априорная и задана выраже x нием (11.39). Остается найти другие две. Для их нахождения опираемся на (11.38). Имеем z = Ax + v = A + v = A + v, x= x= x= так как v от не зависит. Отсюда и из (11.37) получаем 1/ exp ( A)T R1( A). (11.42) fz|x, (, ) = (2)m|R | x Для другой плотности, находящейся в знаменателе (11.41), снова рассмат риваем выражение (11.38), но теперь при условии x =. Найдем совместную плотность векторов x и v при этом условии:

x fx,v, (,, ) fv|x, (, )fx| ( )fx() x x fx,v| (, ) = = = fv ()·fx| (|).

x x fx() fx() (11.43) 11.9 Полная статистическая интерпретация рекурсии в МНК Последнее равенство получено вследствие принятого свойства ошибок v в эксперименте (11.38): они не зависят ни от x, ни от x. В выражении (11.43) перемножаются две гауссовы плотности: (11.37) и (11.39), поэтому резуль тат тоже гауссова плотность. Вектор z, (11.38), образован как взвешен ная сумма двух векторов, x и v, совместная плотность которых гауссова.

Поэтому плотность в знаменателе (11.41) гауссова. Чтобы ее записать, достаточно найти по (11.38) первые два условные момента. Первый момент, учитывая (11.39), равен E z x = = A E x x = + 0 = A, так как E v x = = 0. Находим второй момент, используя (11.39) и (11.37):

(z A)(z A)T = AP AT + R, x= R = Im E (последнее в силу нормализации, см. замечание 11.4). Теперь имеем воз можность записать 1/ fz|( ) = (2)m AP AT + Im exp{·}, (11.44) x где {·} = 2 ( A)T (AP AT + Im )1( A).

Все три плотности для (11.41) найдены: (11.42), (11.39) и (11.44). Под ставляя их в (11.41), получаем 1/ = (2)n R · P exp ( + ), fx|z, (, ) (11.45) x AP AT + R T 1 = ( A) R ( A), T = P 1, (11.46) T.

AP AT + R A A = Замечание 11.5. Здесь и далее для удобства (прослеживания место положения матрицы ковариаций ошибки v) вместо Im, (11.36), сохранено общее обозначение R, как и в (11.37), которое в любой момент можно заме нить на Im, если будет принято предположение о нормализации наблюдений.

11 Оценивание по методу наименьших квадратов Очевидно, (11.45), (11.46) определяют нормальную плотность распреде ления, однако, ее явное определение требует двух действий: (1) приведе ние суммы трех квадратичных форм (11.46) к одной квадратичной форме и (2) приведение отношения трех определителей к одному определителю.

Проведем эти преобразования. При первом действии активно используем лемму 11.1 (см. п. 11.7), беря ее в виде:

1 P 1 + AT R1 H = P P AT AP H T + R AP (11.47) с обозначением = P 1 для априорной информационной матрицы. Пере пишем (11.47) в нескольких эквивалентных формах. Имеем 1 P 1 + AT R1A AT = P AT P AT AP AT + R AP AT = 1 = P AT AP AT + R AP AT + R AP AT = P AT AP AT + R R.

Отсюда 1 P 1 P 1 + AT R1A AT R1 = AT AP AT + R. (11.48) Умножая (11.47) слева и справа на P 1, получаем 1 P 1 P 1 + AT R1 A P 1 = P 1 AT AP AT + R A. (11.49) Умножая (11.47) слева на A и справа на AT и обозначая C = AP AT + R, находим A P 1 + AT R1 A AT = AP AT AP AT C 1AP AT = = (C R) (C R)C 1(C R) = R R AP AT + R R.

Отсюда 1 R1 A P 1 + AT R1A AT R1 = R1 AP AT + R. (11.50) Вычисляя квадратичную форму в (11.45) как ( + ), введем промежу точные обозначения:

, P = P 1 + AT R1A T a = A R, (11.51) b = P 1, c = a + b.

11.9 Полная статистическая интерпретация рекурсии в МНК Теперь после раскрытия скобок в (11.46), приведения подобных членов в ( + ) и при подстановках (11.48), (11.49) и (11.50) получим T + = T P 1 2 T c + cT P c = P c P 1 P c. (11.52) Кроме того, P c = P AT R1 + P 1 = = P AT R1 + P 1 + AT R1A AT R1A = (11.53) P 1 + AT R1A + AT R1 A =P = = + P AT R1 A.

Из (11.48) и (11.51) видно, что P AT R1 = K = P AT AP AT + R (11.54) есть матрица Калмана, определенная выражением (11.32) (см. п. 11.7).

В свою очередь, величина (11.53), входящая в квадратичную форму (11.52) плотности распределения (11.45), есть не что иное как среднее значение век тора x, обусловленное двумя событиями: z = и x = (первое результа том измерения и второе результатом априорной оценки). Обозначим это условное среднее как x, тогда (11.53) перепишется в виде:

x = + K( A). (11.55) Кроме того, из (11.52) видно, что матрица P в (11.51) есть ковариация апо стериорного распределения (11.45), и по лемме (11.47) она дается выраже нием P = P P AP AT + R AP = P KAP. (11.56) Сравнивая (11.54), (11.55) и (11.56) с подразд. 11.7, убеждаемся в том, что исходя из другой чисто статистической задачи оценивания вектора x по зашумленным наблюдениям его компонент через матрицу A (см. рис. 11.4), мы получаем полное алгебраическое совпадение с рекурсией для МНК в стандартной ковариационной форме, если при этом принято во внимание, что z = и x =, а в качестве искомой наилучшей оценки x для вектора x мы приняли апостериорное среднее значение (11.55), что соответствует критерию максимума апостериорной вероятности (МАВ) (11.40), (11.41), (11.45).

11 Оценивание по методу наименьших квадратов v текущие данные z= x Ax x max fx|z, (|, ) + A x задача решение + x= A x + A K априорные данные апостериорные данные P AP KAP P A K + Рис. 11.5. Статистическая задача получения оценки x по критерию МАВ и ее решение (полная статистическая интерпретация МНК-решения). Между моментами получе ния текущих данных апостериорные данные занимают место априорных данных, и процесс решения повторяется Таким образом, при нормальных законах распределения и независимых ошибках наблюдения v МНК-решение x интерпретируется как МАВ-оценка вектора x (рис. 11.5). Именно МАВ-оценка отвечает на поставленный выше вопрос (см. между (11.39) и (11.40)): как наилучшим образом скомбиниро вать априорные данные (, P ) с текущими данными.

Для решения этой задачи, как видно из изложенного, пришлось вычис лить апостериорную плотность (11.41) и затем найти максимизирующий ее аргумент, он оказался равен (11.55). Однако потребовались сложные вычисления, которые пока еще не коснулись отмеченного выше приведения отношения определителей в (11.45) к одному определителю. Из условия нор мировки плотности видно, что он должен оказаться равен |P |, но это еще предстоит доказать. Отложим это доказательство на окончание этого под разд. 11.9 (см. стр. 260), а сейчас рассмотрим, что будет, если изменить кри терий, то есть ради простоты максимизировать по не (11.41), а только одну из плотностей в числителе (11.41), именно, (11.42), что означает критерий максимального правдоподобия, МП. Очевидно, такое решение получается моментально: максимум (11.42) по совпадает с минимумом квадратичной формы в показателе экспоненты. При R = I (что уже обсуждалось как при 11.9 Полная статистическая интерпретация рекурсии в МНК нятое неограничительное условие нормализации зашумленных наблюдений) имеем x = arg min ( A)T ( A) = xМНК = xМП, то есть x есть МНК-решение, xМНК, системы (11.1), переписанной здесь в виде = A + v с R = I, и одновременно, это x есть оценка максимума правдоподобия, xМП.


Тем самым, получена полная статистическая интерпретация МНК решения с точки зрения теории оценок:

1. Если в МНК-решении учтена априорная информация, то полученное решение совпадает с оценкой МАВ, xМАВ (11.55). Эту оценку называют также байесовской, поскольку ее вывод основан на формуле Байеса (11.40). Любой алгоритм (т. е. инструмент) вычисления оценки назы вают оценивателем. Байесовский оцениватель в линейной задаче оце нивания с нормальными законами распределения случайных величин дается формулами (11.55) и (11.56). Он совпадает с формулами филь тра Калмана на этапе обновления оценок по измерениям, и он имеет наглядную схемную интерпретацию на рис. 11.5.

2. Если в МНК-решении не учтена априорная информация, то получен ное решение имеет смысл оценки МП, xМП. Пренебрежение априорной информацией может быть вызвано тем, что она очень мала (недосто верна). Этот факт выражается простой записью: P = 2I при 0, как уже отмечалось на стр. 248.

3. Поскольку недостоверность или скудность априорной информации вы ражается большим ростом диагональных элементов матрицы P по срав нению с внедиагональными элементами: P 1 0 при P 1 = 2 I, 0, то P 1 AT A, что формально видно из (11.51) при R = I. Оцени ватель, который пренебрегает априорной информацией, называют фи шеровским оценивателем в связи с именем R. A. Fisher известного статистика. В этом оценивателе xМАВ AT A AT = xМНК, (11.57) что следует из преобразований выражения (11.53) для (11.55) при использовании (11.51) в пределе при P 1 0 (и с неограничительным условием нормализации наблюдений R = I), если, конечно, матрица 11 Оценивание по методу наименьших квадратов AT A обратима. Если же информационная матрица AT A получен ных измерений z не обратима, тогда вместо (11.57) предельное соотношение записывается в виде xМАВ A+ + I A+A = xМНК, (11.58) что вытекает из теоремы 10.16 (см. подразд. 10.1) и означает сходимость xМАВ к общему решению нормальных уравнений МНК.

Замечание 11.6. В (11.58) первое слагаемое есть нормальное псев дорешение при минимизации z Ax 2, где z =, а второе проекция на N (A).

любого Теперь, в завершение необходимых действий по приведению (11.45) к стандартному виду плотности для нормального закона распределения, кроме (11.52), докажем, что R ·P =P, (11.59) AT + R AP где P определено формулой (11.56). Для этого введем вспомогательную мат рицу P AT P P =.

AP AT + R AP Выполним ее блочное U L-разложение, где U верхняя треугольная мат рица с положительно определенными блоками на диагонали, а L нижняя треугольная матрица с единичной диагональю. Получаем P KAP P AT I P =.

T T K I 0 AP A + R Это позволяет найти определитель для P как произведение определителей диагональных блоков. С учетом (11.56), имеем P = P · AP AT + R. (11.60) С другой стороны, воспользуемся следующей леммой обращения симметрич ной положительно определенной матрицы [115].

Лемма 11.2. При симметричных X1 0, X2 0 и произвольной T X21 = X21 справедливо равенство 1 1 X1 I X1 X X1 X = T X21 X2 0 I X2 X12X1 X I.

T X12X1 I 11.10 Основные результаты Доказательство. Проведите самостоятельно прямой проверкой. Беря определитель в этом выражении, получаем X1 X12 = |X1 |1 · X2 X12X1 X T.

T X12 X Применим этот результат к P, обозначая X1 = P, X2 = AP AT + R, X12 = = P A. Находим P = P · R. (11.61) Сопоставляя (11.60) и (11.61), получаем (11.59). 11.10 Основные результаты 1. Критерий наименьших квадратов min (z Ax)T (z Ax) (11.62) x для решения алгебраических линейных систем общего вида z = Ax + v (11.63) не является изначально статистическим критерием. Решение x, отвеча ющее этому критерию, дается нормальными уравнениями AT A x = AT z (11.64) и имеет вид суммы взаимно ортогональных векторов:

x = A+ z + I A+ A y (11.65) для некоторого y. Это решение единственно тогда и только тогда, когда A+A = I, что означает, что только нулевой вектор составляет ядро N (A) матрицы A (т. е. A имеет полный столбцовый ранг). Нормаль ное псевдорешение x0 есть единственный вектор из (11.65), имеющий минимальную норму, т. е. x0 = A+ z.

2. Статистическая интерпретация задачи (11.62) означает следующее:

(а) x есть постоянный вектор (случайный или нет), который необхо димо оценить по результатам наблюдения z в виде (11.63).

11 Оценивание по методу наименьших квадратов (б) A есть матрица наблюдений, которая показывает, какие линейные комбинации элементов вектора x включены, по условиям задачи, в вектор наблюдения z.

(в) При этом v есть случайная ошибка наблюдения, ковариация кото рой R задана как положительно определенная матрица. Однако, если задача (11.62), (11.63) чисто алгебраическая, R = I.

(г) При статистической интерпретации, матрица AT A в (11.64) назы вается информационной матрицей. Она обратима тогда и только тогда, когда rank A = n (матрица A имеет полный столбцовый ранг). Обратная к матрица называется ковариационной матри цей ошибки оценивания вектора x, P = 1, и она характеризует степень неопределенности в решении задачи оценки (в скалярном случае она равна дисперсии).

3. Статистическая задача оценивания, указанная в предыдущем пункте, приводит к таким же (алгебраически эквивалентным) результатам, как решение алгебраической задачи (11.62), (11.63).

4. Решение задачи (11.62), (11.63), так же как и решение ее статистиче ского эквивалента, может быть получено в рекуррентной форме, удоб ной с вычислительной точки зрения. Возможны две базовые эквива лентные рекуррентные формы (взаимно инверсные): информационная схема оценивания (неявная относительно апостериорной оценки x) и ковариационная схема оценивания (явная относительно x).

5. Оцениватель максимального правдоподобия, выбирает оценку xМП так, чтобы максимизировать функцию правдоподобия (11.42). При гауссо вых независимых ошибках она является решением тех же самых нор мальных уравнений AT R1A = AT R1, x (11.66) которые отвечают критерию взвешенных наименьших квадратов min (z Ax)T R1 (z Ax), z =. (11.67) x По смыслу, xМП есть значение неизвестного вектора x, которое с наи большей вероятностью обеспечивает то событие, что при случайном наблюдении (11.63) результат окажется равным (МП обеспечивает x наибольшую вероятность события z = ). Эта оценка не учитывает никакой априорной информации о векторе x. С другой стороны, оценка, 11.10 Основные результаты удовлетворяющая взвешенному критерию (11.67) и даваемая, соот ветственно, решением взвешенных нормальных уравнений (11.66), называется марковской оценкой в связи с теоремой Гаусса–Маркова, которая утверждает, что эта оценка xМП, удовлетворяющая уравне ниям (11.66), эффективная в классе всех линейных несмещенных оценок, см., например, [49, 17].

6. Оцениватель максимума апостериорной вероятности выбирает оценку xМАВ, в отличие от xМП, принципиально в рекуррентной форме, т. е. учитывает априорную информацию и наилучшим образом комбини рует ее с текущими наблюдениями. Оценка xМАВ вычисляется в явном виде в алгоритме Калмана (11.54), (11.55), (11.56) (см. подразд. 11.7) или в неявном виде в информационном алгоритме (см. подразд. 11.6).

При недостаточной (т. е. крайне неопределенной, скудной) априорной информации, оценка xМАВ будет такой же, как xМП, или, что то же самое, xМНК в задаче (11.67) с решением из нормальной системы (11.66).

7. Рекурсия (рекуррентная схема вычислений) в алгоритмах МНК и оптимального оценивания является краеугольным камнем теории и практики эффективных вычислительных алгоритмов, рассматривае мых ниже в разд. 13 и 14.

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

Одновременное решение нормальных уравнений 12.1 Метод нормальных уравнений Нормальные уравнения (11.3) показаны выше (см. стр. 238). Довольно часто их используют для отыскания МНК-решения x [97]. Этот подход про тивоположен последовательным методам, берущим начало от идеи расщеп ления исходной переопределенной системы на априорную и текущую части (см. подразд. 11.5). Поэтому иногда его называют одновременным реше нием всех нормальных уравнений, характеризующих всю исходную пере определенную систему Ax z (см. стр. 238 и выражение (11.1)).

Алгоритм использует полный вектор Метод Нормальных Уравнений.

наблюдений z Rm и всю матрицу плана эксперимента A Rmn, имеющую rank(A) = n. По матрице A и вектору z алгоритм вычисляет решение x задачи наименьших квадратов, доставляющее min z Ax. Алгоритм:

Вычисляет нижний треугольник матрицы = AT A.

1.

Вычисляет d = AT z.

2.

Вычисляет разложение Холесского = SS T.

3.

Решает сначала систему Sy = d и затем систему S T x = y.

4. Здесь в качестве разложения Холесского может быть взято либо нижнее треугольное разложение (S = L), либо верхнее треугольное разложение (при S = U ) см. подразд. 6.2, стр. 91. Оба эти варианта требуют многократного применения операции извлечения квадратного корня. Если в п. 3 алгоритма вычислять разложение Холесского без операции квадратного корня: либо = SDS T c S = L, либо S = U, то п. 4 заменится на последовательное решение трех систем: Sw = d Dy = w S T x = y.

12.2 Формирование матрицы A 12.2 Формирование матрицы A Используйте следующие четыре варианта формирования матрицы A для оценки скорости и точности двух методов решения переопределенных систем линейных алгебраических уравнений, как указано в пп. В и Г задания (см.

Задание ниже в подразд. 12.3).

Для n = 2 до 12 с шагом 5 выполнять Для m = n + 1 до 26 с шагом 3 выполнять Для i = 1 до m выполнять Для j = 1 до n выполнять (i 1)j Вариант 1: aij = sin, m Вариант 2: aij =, (i 1)j 1 + m Вариант 3: aij = 1/(i + j), Вариант 4: aij = 100(RAN 1/2), где RAN равномерно распределенная в интервале [0, 1] случайная вели чина, получаемая независимо для каждого элемента aij матрицы A.

12.3 Задание на лабораторный проект № А. Спроектировать и отладить подпрограмму решения несовместной сис темы Ax = z, A = A(m, n), m n, rank(A) = n, в смысле наимень ших квадратов при помощи заданного метода ортогонального приведения.

Обосновать проект и дать набор инструкций для пользователей подпро граммы. Сделать подсчет операций (отдельно сложения, умножения, деле ния и извлечения квадратного корня) в зависимости от m и n, где m число строк матрицы A, а n число столбцов. Рекомендуется в качестве основы вашего проекта использовать ту программу, которая была вами написана и отлажена в рамках лабораторного проекта № 6 для решения совместной системы уравнений Ax = f с квадратной матрицей A методом ортогональ ного приведения. Для этого в указанной программе достаточно осуществить незначительные изменения.

12 Одновременное решение нормальных уравнений Б. Повторить п. А задания на основе метода нормальных уравнений A A = AT z, которым удовлетворяет искомое решение x в смысле наи T x меньших квадратов, называемое нормальным псевдорешением несовместной системы Ax z. Для этого применить вашу программу решения системы уравнений с симметричной положительно определенной матрицей методом квадратного корня (разложение Холесского) из лабораторной работы № 5.

В. Спроектировать и провести вычислительный эксперимент для срав нения скорости выполнения двух программ по пп. А. и Б. Использовать четыре различных варианта генерации n векторов длины m для формиро вания матрицы A (см. подразд. 12.2) при 2 n 12, n + 1 m 26.

Результаты представить в виде таблиц и графиков, которые иллюстрируют поведение каждого метода на каждом варианте генерации матрицы A. Дать обобщенную (по вариантам матрицы) картину зависимости времени выпол нения от значений параметров m и n матрицы. Проанализировать соотно шение между фактическим временем выполнения и числом операций, рас считанным по пп. А и Б. Правые части уравнений формировать, как указано в следующем пункте задания.

Г. Подобно п. В, сравнить точность нахождения нормального псевдореше ния переопределенной системы линейных уравнений для методов из пп. А и Б. Для этого также четырьмя способами сгенерировать матрицу A (см.

подразд. 12.2), выбрать (принудительно задать) точное нормальное псев дорешение x и образовать вектор z = Ax. К элементам этого вектора добавить случайные числа vi, чтобы образовать правые части zi = zi + avi, i = 1, 2,..., m. Написать подпрограмму генерации псевдослучайных чисел vi так, чтобы каждое vi, имело стандартное нормальное распределение (с нулевым средним значением и единичной дисперсией). При этом любые слу чайные величины vi, vj (i = j) должны моделироваться как попарно неза висимые. Предусмотреть множитель-переключатель a, чтобы по желанию включать или отключать добавление случайных чисел vi, или же регулиро вать их уровень.

В качестве точного нормального псевдорешения взять вектор x = = [1, 2,..., m]T. Для оценки точности оценивания использовать норму век тора = max(|xi|).

x i При использовании программы, где выполняется ортогональное приведе ние QA = B, для проверки правильности метода убедиться в справедливо 12.3 Задание на лабораторный проект № сти равенства AT A B T B = 0. Для этого использовать норму матрицы n |aij |.

A = max i j= Д. Представить обобщенную аттестацию двух подпрограмм и соответ ствующих методов, основанную на проведенных наблюдениях. Обсудить любые сравнительные достоинства и недостатки, поддающиеся количествен ной оценке, и предложить план дальнейших вычислительных эксперимен тов, которые могли бы помочь уточнить различия между рассмотренными выше методами решения переопределенных систем уравнений.

Е. Решить следующую прикладную задачу [82]. Для i = 1, 2,..., m (m кратно четырем) имеем yi = x1wi + x2wi1, wi = sin(2i/m), di = 2 cos(2i/m).

Найти оптимальное значение x = (1, x2)T вектора коэффициентов x = x = (x1, x2)T, доставляющее минимум средней квадратической ошибке m (yi di )2.

J(x) = m i= Решение выполнить двумя способами: аналитически и численно. Аналити ческое решение должно включать:

1) эквивалентную постановку задачи в виде переопределенной системы Ax z, 2) решение для нее нормальных уравнений, дающее T cosec(2/m), x = 2 ctg(2/m) 3) представление критерия качества для общего случая в виде J(x) = Jmin + (x x)T (x x), информационная матрица, x = (AT A)1AT z где = AT A нор мальное псевдорешение, 4) доказательство того, что в данном конкретном случае Jmin = min(J(x)) = J() = 0, x x 12 Одновременное решение нормальных уравнений 5) вычисление собственных значений (1, 2 ) матрицы, дающее 1 1 cos(2/m), 1 = 2 = 1 + cos(2/m), 2 6) вычисление соответствующих собственных векторов (v1, v2) матрицы, 7) представление критерия качества для общего случая в виде J(x) = Jmin + eT QQT e, где e = x x отклонение x от оптимального значения x, Q матрица ортонормированных собственных векторов матрицы, диагональная матрица = diag [1, 2 ] составлена из собственных значений матрицы, так что = QQ1 = QQT.

Изобразить на экране в системе координат [x1, x2] = xT линии постоянных уровней критерия J(x) = const для шести значений const = {1, 2, 3, 4, 5, 6} в окрестности точки минимума критерия для одного из значений m = = 4, 8, 12, 16, 20, 24, 28 или 32 (по выбору). Объяснить геометрический смысл матриц Q и в последнем представлении критерия.

Численное решение должно включать вычисление решения x с помощью двух методов (по пп. А и Б) и сравнительную оценку точности двух решений для нормы вектора x (в зависимости от m = 4, 8, 12, 16, 20, 24, 28, 32).

Эту зависимость необходимо представить таблицей и графиком.

12.4 Варианты задания на лабораторный проект № По теме Одновременное решение нормальных уравнений студентам предлагается выполнение лабораторной работы проекта № 8.

Задание на этот проект содержит 28 вариантов, которые аналогичны вариантам, приведенным в табл. 7.2 (см. стр. 139) для проекта № 6. Все варианты различаются по следующим признакам:

четыре варианта заполнения треугольной матрицы R;

• три вида ортогональных преобразований:

• 1) отражения Хаусхолдера, 2) вращения Гивенса, 3) ортогонализация Грама–Шмидта, 12.4 Варианты задания на лабораторный проект № Таблица 12.1. Варианты задания на лабораторный проект № Вариант Отражения Вращения Ортогонализация заполнения Хаусхолдера Гивенса Грама-Шмидта матрицы R a b a b c d e Rne 1 2 3 4 5 6 Rnw 8 9 10 11 12 13 0 Rse 15 16 17 18 19 20 0 Rsw 22 23 24 25 26 27 столбцово-ориентированный алгоритм;

a строчно-ориентированный алгоритм;

b классическая схема;

c модифицированая схема;

d модифицированая схема с выбором ведущего вектора.

e две разновидности алгоритма ортогонализации по методу Хаусхолдера • и по методу Гивенса:

1) столбцово-ориентированный алгоритм, 2) строчно-ориентированный алгоритм, три разновидности ортогонализации по методу Грпма–Шмидта:

• 1) классическая схема, 2) модифицированная схема, 3) модифицированная схема с выбором ведузего вектора.

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

Устойчивые алгоритмы фильтрации 13.1 Фильтрация Калмана в историческом аспекте Вскоре после открытия Калманом своего ставшего знаменитым нового подхода к задачам фильтрации и предсказания [101] выяснилось, что это изящное решение, которое теперь принято называть обыкновенным или стандартным фильтром Калмана (the conventional Kalman lter, CKF), работает хорошо лишь в исключительных так называемых хорошо обу словленных задачах и расходится в большинстве практических задач.

Явление расходимости теоретического алгоритма CKF, детально исследо ванное 10 лет спустя [98], породило поиски альтернатив, алгебраически экви валентных CKF, но в вычислениях значительно более устойчивых.

В 1963 году Поттер нашел первое решение [125] в связи с разработкой аме риканского лунохода (Lunar Excursion Module, LEM) для Программы Апол лон. По природе, это квадратно-корневой алгоритм, названный квадратно корневым фильтром Поттера, PSRF (Potter Square Root Filter). Будучи приспособлен лишь для ограниченных случаев (некоррелированность ска лярных измерений и отсутствие шума в уравнении состояния), PSRF по родил две ключевые идеи: (1) факторизация (разложение на множители) ковариационных матриц, чтобы устранить опасность потери положительной определенности, и (2) декорреляция и последующая скаляризация вектор ных измерений, чтобы устранить операцию вычисления обратной матрицы в алгоритме CKF.

Любой квадратно-корневой алгоритм фильтрации численно более устой чив, чем CKF, так как число обусловленности квадратного корня из мат рицы есть квадратный корень из числа обусловленности соответствующей матрицы. Этот факт широко известен и вошел в учебные курсы численных методов алгебры (разложение Холесского) cм. выше разд. 6 и [74].

13.1 Фильтрация Калмана в историческом аспекте Другая идея, также вошедшая в учебные курсы, ортогонализация мат риц (cм. выше разд. 7 и [74]) развивалась параллельно для численного решения плохо обусловленных и переопределенных систем по методу наи меньших квадратов (МНК) [112]. Важный промежуточный результат тех лет опубликование обзора [104] и монографии Бирмана [15].

Кроме исходных аэрокосмических приложений, эти и другие современные методы калмановской фильтрации вошли в другие области: эконометрика [93], метеорология [110], спутниковая геодезия [91], телекоммуникационные сети [121], высоконадежная спутниковая навигация [111], обработка изобра жений в реальном времени [124] и многое другое.



Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 8 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.