Дж. Деммель - Вычислительная линейная алгебра, страница 14
Описание файла
PDF-файл из архива "Дж. Деммель - Вычислительная линейная алгебра", который расположен в категории "". Всё это находится в предмете "квантовые вычисления" из 7 семестр, которые можно найти в файловом архиве МГУ им. Ломоносова. Не смотря на прямую связь этого архива с МГУ им. Ломоносова, его также можно найти и в других разделах. .
Просмотр PDF-файла онлайн
Текст 14 страницы из PDF
Старая гипотеза, согласно которой дор < п, недавно была опровергнута [99, 122). Все еще остается открытой задача построения хорошей верхней оценки для дср (который,как по-прежнему полагают, является величиной порядка 0(п)). Дополнительные 0(п ) сравнений, производимых в методе СЕСР при выборе главных элементов (0(пз) сравнений на шаг, в то время как в методе СЕРР достаточно 0(п) сравнений на шаг), делают этот метод значительно более медленным, чем метод СЕРР, особенно на высокопроизводительных компьютерах, где операции с плавающей точкой выполняются почти так же быстро, как сравнение.
Поэтому обращение к методу СЕСР редко бывает оправданным (см., однако, разд. 2.4.4, 2.5.1 и 5.4.3). Пример 2.3. Рисунки 2.1 и 2.2 иллюстрируют оценки обратной ошибки, о которых шла речь выше. Для обоих рисунков генерировались пять случайных матриц А каждой указанной размерности; элементы матриц суть независимые нормально распределенные величины со средним 0 и средним квадратичным отклонением 1. (Тестирование на подобных случайных матрицах может подчас неверно ориентировать в отношении поведения метода для некоторых реальньгх задач; все же такое тестирование достаточно информативно.) Для каждой матрицы генерировался случайный вектор Ь с аналогичными характеристиками.
Для решения системы Ах = Ь использовались оба метода СЕРР и СЕСР. На рис. 2.1 изображены коэффициенты роста дрр и дор. Как и следовало ожидать, оба медленно растут с ростом размерности. Рис. 2.2 показывает поведение двух наших верхних оценок для обратной ошибки: Зпзвдрр (или Зпзвдор) и Зпв ))',, —. Показана также подлинная обратная ошибка, вычисляемая, бО Глава 2. Решение линейных уравнений КоэФфициенты роста: метод ОЕРР = О, метод ОЕСР =е В О 0 10 20 ЗО 40 ВО ВО 70 ВО 90 100 Размерность матрицы Рис.
2.1. Коэффициенты роста для случайных матриц, о = дее, + = дое. Обратная ошибка в гауссовом исключении с частичным выбором главного элемента 10 10 10-12 толк 10 0 10 20 30 40 60 60 70 60 90 100 Размерностьматрицы и Обратная ошибка в гауссовом исключении с полным выбором главного элемента 10 !Оик том4 10 0 10 20 90 40 60 60 70 60 90 100 Размерностьматрицы Рис. 2.2. Обратная ошибка в гауссовом исключении для случайных матриц, х = Зпзед, + = ЗпЦ )Ц . (Ьтщ', ДАО, „о = ЗАх — Ьб„/НАО„'б~хб, ).
61 2.4. Анализ ошибок как это описано в теореме 2.2. Машинное эпсилон указано сплошной горизонтальной прямой на уровне е = 2 зз 1.1 10 'е. Обе оценки действительно являются верхними границами для подлинной обратной ошибки, но переоценивают ее на несколько порядков. По поводу Ма11аЬ-программы, производящей эти графики, см. НОМЕРАСЕ/Маь1аЬ/ргкоь.ш. О 2.4.3. Оценка чисел обусловленности Чтобы вычислить практичную оценку ошибки, основанную на неравенстве типа (2.5), нужно уметь оценивать число ~)А г ~!. Этого достаточно, чтобы суметь оценить н число обусловленности к(А) = 11А (( 'йА!), поскольку бАй' находится легко.
Одна из возможностей состоит в том, чтобы вычислить матрицу А ' в явном виде, а затем определить ее норму. Однако это стоило бы 2нз операций, больше, чем 1~из операций, уже выполненных для гауссова исключения. (Отметим следствие из сказанного: нет никакого выигрыша в решении системы Ах = Ь путем вычисления матрицы А г с последующим умножением ее на 6. Это верно даже в том случае, когда при одной и той же А имеется много различных векторов Ь. См, вопрос 2.2.) Фактом является то, что большинство пользователей не станут вычислять оценки ошибки, если они обходятся так дорого. Поэтому вместо вычисления А ' мы разработаем гораздо более дешевый алгоритм оценивании числа 11А г11.
Такой алгоритм называется оцени1иком обусловленности и должен обладать следующими свойствами: 1. Он должен использовать 0(пз) операций, если заданы треугольные множители 1 и У матрицы А. Для достаточно больших п эта цена пренебрежимо мала по сравнению с з из операций, необходимых для метода СЕРР. 2. Он должен давать оценку, почти всегда не превосходящую удесятиренного значения 11А ' ~!. Это все, что требуется от оценки ошибки, которая должна говорить нам, сколько верных знаков в решении задачи мы имеем. (Ошибка, соответствующая множителю 10, равносильна одному потерянному десятичному разряду ') Имеется много таких оценщиков (их обзор дан в [146]). Мы выбрали для описания метод, применимый не только в решении системы Ах = Ь, но н в ряде других задач; вследствие такой универсальности, он несколько медленней алгоритмов, специализированных для задачи Ах = Ь (хотя и для этой задачи он достаточно быстр).
Подобно большинству других оценщиков, наш алгоритм дает нишснюю границу для 11А г11, а не оценку сверху. Эмпирически эта оценка отличается от 11А '11 множителем, почти всегда не меньшим, чем 1/10 н обычно находящимся в пределах от — до —. В случае матриц, использованных 1 1 для рис. 2.1 и 2.2 (их числа обусловленности варьируются от 10 до 10з), наш оценщик давал несколько верных десятичных разрядов числа )~А г11 в 83% 1 Как уже отмечалось, никому не удалось построить оценщик, который бы аппроксимировал число 11А г1~ с некоторой гарантированной точностью и был бы в то же время значительно дешевле явного вычисления матрицы А г.
Была выдвинута гипотеза, что такого оценщика не существует, но она не была доказана. Глава 2. Решение линейеых уравнений ][В[)» = п»ах = шах~~ [ЬО[. [[Вх][» хФО )[х[)» Легко показать, что максимум по ненулевым векторам х достигается здесь для х = е е — — [О,..., О, 1, О,..., О]» . (Если шах, 2, [Ь,Д достигается при у = »о, то единственная ненулевая компонента вектора х имеет номер»а.) Поиск по всем е, (у = 1,..., п) означает вычисления всех столбцов матрицы В = А ', что слишком дорого. Поскольку [~В[[» = шахг,1г<» ~[Вх[[», мы можем вместо перебора е. применить метод наискорейшего падвема к функции /(х) = [~Вх[~» на множестве [[х[[» < 1.
Последнее, очевидно, является выпуклым множеством, а /(х) — выпуклой функцией. Действительно, для 0 < с» < 1 имеем /(с»х+ (1 — с»)у) = [[с»Вх+ (1 — с»)Ву[[» < с»[]Вх]~» + (1 — »»)][Ву~[» = о/(х) + (1 — о)У(у). Максимизация /(х) методом скорейшего подъема означает движение х в направлении градиента»7/(х) (если он существует), пока /(х) возрастает. Из выпуклости /(х) следует, что /(у) > /(х) + с7/(х) (у — х) (если с7/(х) существует). Чтобы вычислить»7/, предпалоэ»сиэл, что в сумме /(х) = ~, [ 2 Ьйхд[ все слагаемые 2,'1 Ьбх отличны от нуля (это верно почти всегда). Положим ~» = гйп2'1Ь„хд, тогда»,» = х1 и /(х) = 2 '» 2 ~»Ьихд.
Поэтому, = 2,'»»,»Ь»ь и С~У = (.'В = (В О . Итак, вычисление градиента»7/(х) состоит из трех шагов: ю = Вх, », = Бкп(ю) и»7/ = »,тВ. Алгоритм 2.5. Оценщик обусловленности Хэйдэ»сера (Науег) выдает ниэ»с- нюю границу []ю)[» длл [)В)[». /'наприл»ер, хн=-„''/ /'х'= 7/ / выбрать произвольный век»пор х с нормой [[х[~» =1 гереа» ю =Вх, », =зйп(ю), г =Вт», »/[ф[, < гхх»йеп ге»игп [[ю[[» ейе х = ед где [гу[ = [[г[[, е иЮ епд тереа» экспериментов; наихудшая оценка составляла 0.43 от правильного результата. Такой точности более чем достаточно для того, чтобы иметь возможность оценить число верных знаков в полученном решении.
Алгоритм оценивает 1-норму [~В[~» матрицы В; предполагается, что мы можем вычислять произведения Вх и Вгу для любых векторов х и у. Мы применим алгоритм к В = А», поэтому потребуется вычислять векторы А»х и А»у, т. е. решать линейные системы. При имеющемся ЕС-разложении матрицы А решение таких систем обходится лишь в О(пг) операций. Алгоритм был разработан в [138, 146, 148], а его новейшая версия описана в [147]. Напомним, что ][В[~» определяется как 63 2.4. Анализ ошибок Теорема 2.6. 1. Если на выходе алгоритма получено значение ЙшЙд, то 11ш))д — — 11Вх)(д есть локальный максимум функции !)Вх11д. 2.
В пРотивном слУчае, ()ВедЙ (значение в конЦе цикла) > 11Вх)( (значение в ддачале цикла). Тем самым значение функции 1'(х) увеличено. Доказательство. 1. В данном случае йг11, < гтх. Вблизи х функция ((х) = 11Вх~(д .д,;Ьбх. линейно зависит от х, поэтому ((у) = д(х)+\уд(х) (у — х) = ,((х) + г~(у — х), где гт = ду,((х). Чтобы показать, что х — локальный максимум, проверим, что г~(у — х) < 0 при ~)у)~д — — 1. Имеем г (у — х) = г у — г х = ~гд у; — г х < гэ 1г1) )уд! — г х < 11г11, . ()у11д — г~х = )(г)), — г~х < О, что и требовалось. 2.
В данном слУчае ЙЦ, > гтх. Положим х = е в18п(Яд), где (Уд) выбРано так, что 11гд.11 = 11г11 . Тогда ,((х) > ((х) + д7 ( (У вЂ” х) = ((х) + г (х — х) = 1(х) + г х — г х = 1(х) + 1гд/ — г х > ((х), где последнее неравенство выполняется вследствие выбора ~. В (147, 148) описаны результаты тестирования несколько улучшенной версии алгоритма на многочисленных случайных матрицах порядков 10, 25 и 50 н с числами обусловленности к = 10,10з,10е, 10э.