§ 11. Погрешность приближенного решения системы уравнений и обусловленность матриц. Регуляризация

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

в действительности должна была решаться некоторая система

Пусть известны оценки и . Займемся оценкой погрешности решения.

Сначала выделим главный член погрешности. Будем обозначать решения (1) и (2) через X и и разность — через . Подставив выражения в (2), будем иметь

Вычитая из этого равенства (1), получим

откуда

b

Если и малы, то следует ожидать и малости . Тогда слагаемое имеет более высокий порядок малости; отбрасывая слагаемое, получаем

Отсюда следует оценка погрешности

(4)

Строгая оценка погрешности получается следующим образом. Вследствие (3) выполняется неравенство

Предположим, что . Перенеся последнее слагаемое в левую часть и поделив неравенство на коэффициент при , получим оценку

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

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

Соответственно этому в качестве меры обусловленности системы принимается число

Отсюда получаем оценку относительной погрешности решения через меру обусловленности системы и относительную погрешность правой части:

Так как , то

и

Иногда удобнее иметь более грубую характеристику свойств системы только через свойства матрицы А. Эту характеристику называют мерой (или числом) обусловленности матрицы А. Согласно этому определению и (6), имеем оценку

связывающую относительные погрешности правой части и решения только через свойства матрицы системы. Так как

то

Поскольку любая норма матрицы не меньше ее наибольшего по модулю собственного значения, то поскольку собственные значения матриц А и взаимно обратны, то

Таким образом,

В частности, при имеем и

Следовательно, в случае нормы

Рассмотрим вопрос о погрешности решения вследствие округления в ЭВМ правой части. Пусть, как обычно, t -двоичная разрядность чисел в ЭВМ. Каждый элемент правой части округляется с относительной погрешностью , т.е. с абсолютной погрешностью, равной , поэтому

Следовательно,

При практической работе вопрос: о строгой оценке погрешности полученного приближенного решения системы линейных уравнений с помощью полученных неравенств или каким-либо иным способом возникает редко. Однако информация о порядке погрешности решения часто полезна для получения качественных выводов о том, с какой точностью разумно решать задачу. Соотношения (4), (5) оценивают сверху погрешность решения, являющуюся следствием погрешности исходных данных. Из равенства (3) видно, что оценки (4), (5) довольно точны, поэтому обычно не имеет смысла стремиться получать решение задачи с погрешностью, существенно меньшей чем .

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

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

При реально заданной правой части решением будет

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

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

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

Первый способ. Задавая некоторое , находим решение системы

записывается в виде

Так как

то наличие малого параметра а несущественно изменит слагаемые, соответствующие большим . В то же время при имеем

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

Второй способ заключается в следующем. Будем решать систему уравнений каким-либо итерационным способом. Рассмотрим случай итераций по формуле

при некотором начальном приближении . Пусть

Подставим эти выражения в (7) и, приравнивая коэффициенты при , получим соотношения

Последовательно выражая каждое через предыдущие, имеем

Если , то при

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

В других случаях решение задачи находят, минимизируя некоторый функционал, близкий к , например функционал с малым .

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

при малых значение опять-таки подбирается экспериментально из сравнения результатов расчетов при различных .

Другая группа методов основана на представлении матрицы системы А в виде

где G и Р — ортогональные матрицы, а — двухдиагональная матрица, у которой могут быть отличными от нуля элементы при и .

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

Задача 1. Пусть — матрица размерности с элементами

1. Вычислить матрицу и доказать утверждение: при матрицы в некотором смысле хорошо обусловлены, а при и больших плохо обусловлены.

2. Выписать явно решение системы через правую часть.

3. Выписать явно через правую часть вектор , минимизирующий функционал

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

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

Если , то и все собственные значения матрицы А равны между собой по модулю, т.е. . За исключением этого частного случая , поэтому и . Таким образом, в случае симметричной матрицы применение симметризации увеличивает число обусловленности. Число обусловленности является непрерывной функцией от матрицы, поэтому для матриц, близких к симметричным, применение симметризации также увеличивает число обусловленности матриц. Таким образом, имеется непустой класс несимметричных матриц по отношению к которым алгоритм симметризации увеличивает число обусловленности.

Задача 2. Существуют ли несимметричные матрицы, для которых .

Рассмотрим еще один метод решения плохо обусловленных систем линейных алгебраических уравнений. Пусть

Относительно А будем считать, что в спектре матрицы есть как собственные числа порядка 1, так и собственные числа, близкие (или даже равные) к нулю. Это как раз и означает, что матрица А плохо обусловлена.

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

Назовем решением X уравнения (8) вектор, который минимизирует функционал невязки, а именно,

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

Уравнение (10), в отличие от (8), всегда имеет решение. Действительно, непосредственной проверкой убеждаемся, что . Необходимым и достаточным условием существования решения линейной системы уравнений (10) является ортогональность правой части ядру матрицы системы, т.е. вектор должен быть ортогонален ядру , которое, как мы отметили выше, совпадает с ядром . Но из вида правой части видно, что она действительно ортогональна ядру А. Таким образом, система (10) всегда имеет решение. В общем случае таких решений может быть несколько.

Описываемый ниже метод заключается в минимизации функционала методом оптимального покоординатного спуска; на каждом шаге выбирается координата, спуск по которой будет оптимальным в смысле минимизации . В качестве координатных (базисных) векторов можно выбрать любую ортонормированную систему.

Пусть — ортонормированная система векторов в (не обязательно базис) и по крайней мере для одного из векторов. Обозначим через W линейную оболочку векторов . Будем искать вектор, минимизирующий функционал невязки на подпространстве W. Для этого рассмотрим следующий итерационный метод. Положим

Если приближение уже найдено, то следующее приближение будем искать в виде , где

(11)

Наряду с приближениями введем невязки

Выпишем условия минимума функционала по . Имеем

Заметим, что при поиске минимума достаточно рассматривать только те , для которых так как в противном случае значение функционала не меняется. Функция , как функция переменной , является многочленом второй степени, причем коэффициент при положителен в силу замечания выше. Отсюда следует, что минимум по при фиксированном j существует и единствен, если Таким образом, из (13) следует, что удовлетворяет уравнению

откуда

При таком выборе

и

Суммируя вышесказанное, получаем следующий алгоритм:

1) Вычисляем векторы , и их нормы в дальнейшем рассматриваем только те векторы , для которых . Не уменьшая общности, будем считать, что число таких векторов равно .

2) Выбираем на основе априорной информации; в частности, можно взять .

3) Если найден, то вычисляем по формулам

(14)

4) Следующее приближение вычисляем по формуле

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

Этап 3 итерационного метода требует , а этап 4 — арифметических операций.

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

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

Исследуем сходимость итерационного метода. Имеет место

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

Доказательство. Положим

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

Пусть для определенности . Тогда справедлива цепочка неравенств

откуда и следует непрерывность рассматриваемого функционала.

Предположим, что утверждение леммы неверно. Тогда существует последовательность такая, что , где при . Так как в конечномерном пространстве сфера компактна, то существует сходящаяся подпоследовательность. Для простоты изложения предположим, что сходится сама последовательность

В силу непрерывности функционала имеем .

Следовательно, при имеем

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

Теорема. Последовательность приближений , получаемая в ходе итерационного метода (14), (15), является фундаментальной и сходится к некоторому вектору, минимизирующему функционал невязки на подпространстве W, со скоростью геометрической прогрессии. А именно, существует такое, что

Постоянная q при этом зависит от выбора базиса и оператора А.

Доказательство. Так как , то существует постоянная такая, что

Из (14), (15) следует, что невязки удовлетворяют соотношению

Положим . Тогда (17) примет вид

Поскольку выбиралось из условия минимума , то

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

Из (14) и (15) имеем

откуда, учитывая (16), получаем оценку

Применяя к полученному неравенству оценку (18), имеем

откуда следует цепочка соотношений

Таким образом, последовательность является фундаментальной и имеет предел . В силу построения, минимизирует функционал на подпространстве W. Теорема доказана.

Описанный выше метод решения систем уравнений с плохо обусловленными матрицами особенно эффективен в случае, когда априорная информация представлена в виде каких-либо сведений о структурных особенностях искомого решения; например, когда известны базисные функции , и решение представимо в виде разложения по небольшому числу этих функций. Такая ситуация часто имеет место в задачах цифровой обработки сигналов.

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

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