- •2. Применение алгоритма
- •3. Список формальных параметров
- •2. Применение алгоритма
- •3. Список формальных параметров
- •4. Программы на языке алгол
- •5. Организация процедур и обозначения
- •6. Оценка точности решения
- •7. Примеры использования процедур и результаты их проверки
- •1. Теоретические предпосылки
- •2. Применение алгоритма
- •3. Список формальных параметров
- •4. Алгол- Программа
- •Integer p, 0, /, /;
Из книги алгоритмов линейной алгебры формулы частично вписаны
Алгоритм 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 e0 (е0 — машинная точность), можно гарантировать, что все б# будут меньше, чем е0 ||J(0)||oo. Элементами в матрице J, которые не превышают указанной величины, можно пренебречь. Следовательно, матрица J распадается на две части, которые можно рассматривать независимо.