Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:
Коршунов лабораторные / ЛинАлгбр алгол-пргр.rtf
Скачиваний:
17
Добавлен:
26.04.2015
Размер:
392.72 Кб
Скачать

Из книги алгоритмов линейной алгебры формулы частично вписаны

Алгоритм 1.9 . Решение систем уравнений и проблема наименьших квадратов

1. Теоретические предпосылки

Пусть А — матрица, имеет п строк и т столбцов, причем т <=n. Тогда для произвольного вектора b найдется единственный вектор х, минимизирующий евклидову норму || b - Ах || = min ||b - Аξ||, если только столбцы исходной матрицы А линейно независимы.

Справедливы следующие утверждения:

1. Вектор L -Ах ортогонален столбцам матрицы А, т. е. выполнено условие Ат (b - Ах) =0*.

2. Вектор х может быть определен из выражения х = А+ b, где матрицу А+, равную ТА) -1АТ, называют псевдообратной.

3. Вектор b - Ах можно представить как некоторую проекцию на вектор b: b - Ах = Нb, где Н = I- А (АТА)-1АТ — эрмитова матрица, удовлетворяющая условиям Н2 = Н и НА = 0.

Обычно вычислительная процедура начинается с решения так называемого нормального уравнения АтАх == АТb, (1)

которое следует из утверждения 1.

Вместо уравн. (1) можно использовать уравнение (AS)T Ах = (AS)T b, где S - невырожденная матрица размера т x т. Для удобства рекомендуется выбирать матрицу S таким образом, чтобы матрица U = AS образовывала ортонормированный базис в пространстве, порожденном столбцами матрицы А. Матрица UTA полученной системы линейных уравнений UTAx = UTb будет иметь то же число обусловленности **, что и исходная матрица А, и, вообще говоря, лучше обусловлена, чем матрица АТА в уравнении (1).

Применение процесса ортогонализации матрицы А для получения решения методом наименьших квадратов описано в [38].

Ниже приведена процедура с незначительными отличиями от описанного алгоритма. Она основана на представлении исходной матрицы в виде А = UDR, где матрица U состоит из т взаимно ортогональных ненулевых столбцов; R — невырожденная верхняя треугольная матрица размера mxn, a D = (UTU)~l. Такое разложение существует, если столбцы матрицы А линейно независимы.

_________________________________________________________________

Так как рассматриваем ниже процедуру для случая действительных матриц АS, то знак ** эрмитовой сопряженности заменен на знак транспонирования (Прим. ред. пер.)

Число обусловленности прямоугольной матрицы определяют : cond(A) = max || Ax ||/min || Ах||

||х||=1 ||х||=1

Для евклидовой нормы достаточно показать, что (UTA)T U та= а ТA . Учитывая, что S неособая матрица, получим следующие соотношения (UTA)T UTA,= ATUUTA = ATUUTUS-1 = ATUS-1 = ATA. Если матрица U не нормирована, то минимум числа обусловленности [5] остается неизменным.

=============================================================================

Условие UT (b - Ax) = 0 определяет невырожденную треугольную систему уравнений Rx = UTb, коя решается методом обратной подстановки. Заметим, что матрица r = UтА и вектор UTb получены с помощью одного и того же преобразования UT, примененного к матрице А и вектору b соответственно.

Вычислительный процесс приведения к форме Rx = UTb осуществляется за m шагов. На i-ом шаге будет вычислена матрица A(t) по известной матрице А(i-1) путем ортогонализации относительно t-гo столбца столбцов, следующих за ним. Выбрав начальную матрицу А(0) как объединение исходной матрицы А и совокупности векторов правых частей А(0) - [А | b1 |bb2 | . . . | bk ], получим в процессе ортогонализации на i-1-м шаге матрицу А(i-1), все столбцы которой ортогональны любому из i - 1 первых столбцов. Кроме того, первые i - 1 столбцов взаимно ортогональны. Для получения матрицы А(1) необходимо провести ортогонализацию столбцов матрицы, следующих за i-ым, относительно i-го столбца.

Обозначим столбцы матриц А(i-1) ** и А(i) через ak (i-1) и ak (i) соответственно. На i-ом шаге преобразуются лишь столбцы ak (i) для k> i в соответствии с выражением ak (i) = ak (i-1) - α ak (i-1) (2)

где коэффициент α необходимо выбирать так, чтобы векторы ak (i), и ak (i-1) были ортогональны.

Из условия ортогональности получаем выражение для коэффициента α (в фигурных скобках). Подставляя в (2) получаем окончательно: ak (i) = ak (i-1) - α ak (i-1) = ak (i-1) - {(ai (i-1) [ ai (i-1) ]T ) / ([ ai (i-1) ]T ai (i-1) )} ak (i-1) (3)

Это соотношение (3) справедливо и при k <= i. Если k < i, то из-за ортогональности векторов ai (i-1) и ak (i-1) преобразовывать столбцы с номерами k< I не нужно. А при k = i правая часть обращается в ноль. В ходе вычислений, конечно, преобразуются только столбцы, следующие за i-ым.

Таким образом, реализуемый процесс вычислений, по существу, представляет процедуру исключения Гаусса. Однако одновременно это и процесс ортогонализации, так что после т шагов первые т столбцов матрицы Аш образуют ортогональный базис U такой, что U = AS, где S — невырожденная верхняя треугольная матрица с единичной диагональю. (По существу, это процедура ортогонализации Шмидта, если не учитывать, что в последнем случае ортогонализация вектора выполняется лишь 1 раз.)

Остальные столбцы матрицы А(т), обусловленные наличием правых частей b1 . . ., bk будучи ортогональны столбцам матрицы U, одновременно ортогональны столбцам матрицы А и поэтому образуют разности b1 — Ах1, b2 — Ах2, . . bk — Ахk.

Более того, оператор U преобразуeт матрицу [А | b 1 | b2|… | bk] автоматически в процессе исключения, причем i-ая строка матрицы [UTA|UTb1 | UTb2. |Uтbk] есть строка, исключаемая на i шаге алгоритма. Это следует из того, что, во-первых, i-ая строка есть совокупность скалярных произведений i-го столбца матрицы U и столбцов матрицы [Aj |bj | b2 | . . .| bk] и равна скалярным произведениям i-го столбца матрицы А(i) на столбцы матрицы А(0). Поскольку соответствующие столбцы матриц А(0) и A(i) отличаются, возможно, только вектором, принадлежащим пространству первых t — 1 столбцов А(i), которые ортогональны i-му столбцу A(i), то элемент i-ой строки матрицы [UTA|UT b1 |UTb2 |..| UTbk] есть скалярное произведение i-го и j-го столбцов матрицы А(i). Однако это и есть та строка, которая исключается на i-ом шаге.

В результате, запоминая исключенные строки, сформируем систему уравнений с невырожденной верхней треугольной матрицей R: Rx = UTb = v, из которой и получаем решение х.

Применяя метод исключения Гаусса к нормальной системе уравнений, получаем эту же самую систему, ибо АТА= RTDR. Если матрица U ортонормирована, то R представляет собой верхнюю треугольную матрицу в разложении Холецкого для АТА.

2. Применение алгоритма

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

Процедуру ortholin 1 можно использовать для вычисления решения X на основе метода наименьших квадратов для системы уравнений АХ = В, где А —• матрица размера п X t, X — матрица т X k, В — матрица п X k. Если k =• = 1, то В есть вектор-столбец, и специально для этого случая процедура ortholin 2. Процедура ortho 1 это обращение матриц.

Учитывая рекомендации работ [12, 33], в эти три процедуры включено итерационное уточнение полученного решения. В отличие от них в процедуре ortho 2 итерационное уточнение не реализовано, зато в ней предусмотрена существенная экономия памяти за счет того, что массив исходной матрицы А использован для размещения промежуточных результатов, а затем и элементов обратной матрицы.

3. Список формальных параметров

Входные параметры процедуры ortholin /:

а — массив размера [1 : n, 1 : т] для размещения элементов матрицы системы линейных уравнений; этот массив должен быть описан во внешнем блоке^

п — число строк матрицы А;

т — число столбцов матрицы А, ^ п);

b — массив размера [1 : n, 1 : k] для размещения магрицы В правых частей системы линейных уравнений; массив должен быть описан во внешнем блоке;

k — число столбцов матрицы В;

eps — максимальная относительная ошибка округления.

Выходные параметры процедуры ortholin /:

х — массив для размещения матрицы решений размера (т X k)\ массив должен быть описан во внешнем блоке;

stop — метка во внешнем блоке для выхода из процедуры, если итерационное уточнение оказывается неэффективным.

Нет стр 116-121 с алгол-пр.

Алгоритм 1.10 стр128

Обычный QL-алгоритм со сдвигом можно записать в следующем виде: М— sI = TsRs; RsTs -f- si = ms, (4)

где Tj Ts = I и Rs — верхняя треугольная матрица. Следовательно, Ms = = T^MTS. В работе [26] показано, что не обязательно выполнять вычисления согласно соотношениям (4); этот сдвиг можно осуществить неявным образом. Пусть Т — произвольная матрица, та- (k= 1,2, . . ., /г), т. е. элементы первого столбца Ts равны элементам первого столбца Т и ТТТ <* I.

Тогда справедлива следующая теорема [26]: если поддиагональные элементы матрицы М ненулевые и М = Т^МТ, то матрицу М можно вычислить следующим образом: М = DMSD, где D — диагональная матрица, элементы которой суть ±

Рис. 2. Изменение структуры двухдиагональной матрицы при выполнении п преобразований вращения

Таким образом, матрицу Т2 в последовательности (3) надо выбирать так, чтобы ее первый столбец был пропорционален столбцам М — si; то же справедливо и для первого столбца матрицы произведения Т = Т2Т3 . . .Тя. Следовательно, если поддиагональные элементы М ненулевые, условия теоремы Франциса выполнены, и матрица Т совпадает с матрицей Ts (с точностью до знака элементов в столбце). Итак, преобразование (2) эквивалентно Q/^-преобразованию JTJ с заданным сдвигом s.

Параметр сдвига s определяется собственным значением нижнего минора (размера 2X2) матрицы М. При таком выборе параметра s метод обладает глобальной и почти всегда кубичной сходимостью [99 ].

Проверка сходимости алгоритма

Если элемент еп двухдиагональной матрицы удовлетворяет условию п\ ^ ^ б, где 6 — заданная точность, то | qn \ можно считать сингулярным числом, и тогда порядок исследуемой матрицы уменьшается на единицу. Если же выполнено условие \*k\ ^ б для k =f= n, то матрица распадается на две, и сингулярные числа вычисляются независимо.

Если qk = 0, то, по крайней мере, одно из сингулярных чисел должно быть равно нулю. В этом случае при отсутствии ошибок округления матрица распадается при выполнении преобразования с нулевым сдвигом. Теперь допустим, что на каком-то шаге |^/?|^б. В этом случае матрица J будет подвергнута дополнительным преобразованиям, проводимым над строками (k, k + 1), (k, k + 2), . . ., (k, n), которые состоят в следующем:

элемент {J} k, k+i будет аннулирован, а элементы {J}^, &+2> PU+i, k добавлены в результате левого преобразования матрицы J с помощью плоского вращения;

элемент {J}&, &+2 аннулирован, но добавлены элементы {J}^, k+z> и т. д., наконец,

элемент {J}&, n аннулирован и добавлен последний элемент (J}n,

стр 129

Полученная таким образом матрица имеет вид

*

Заметим, что в силу ортогональности

Таким образом, выбирая б = ||J(0)||oo e00 — машинная точность), можно гарантировать, что все б# будут меньше, чем е0 ||J(0)||oo. Элементами в матрице J, которые не превышают указанной величины, можно пренебречь. Следовательно, матрица J распадается на две части, которые можно рассматривать независимо.

Соседние файлы в папке Коршунов лабораторные