автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.11, диссертация на тему:Расширение возможностей компьютерно-алгебраической системы Maple для решения линейных задач метода наименьших квадратов

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М.В. ЛОМОНОСОВА

Факультет вычислительной математики и кибернетики

На правах рукописи

Матин Фар Машалла Набиолла

Расширение возможностей компьютерно -алгебраической системы Maple для решения линейных задач метода наименьших квадратов

Специальность 05.13.11 — математическое и программное обеспечение вычислительных машин, комплексов и компьютерных сетей

Автореферат

диссертации на соискание ученой степени кандидата физико-математических наук

Диссертация выполнена на кафедре общей математики факультета вычислительной математики и кибернетики Московского государственного университета им. М.В. Ломоносова

Научный руководитель — доктор физико-математических наук, профессор Х.Д. Икрамов

Официальные оппоненты:

— доктор физико-математических наук, профессор С.А Абрамов

— кандидат физико-математических наук М. В. Уфимцев

Ведущая организация — Институт вычислительной математики РАН

Защита диссертации состоится 28 мая 2004 г. в 11 часов на заседании диссертационного совета Д 501.001.44 в Московском государственном университете им. М.В. Ломоносова по адресу: 119992, ГСП-2, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет ВМиК, аудитория 685.

С диссертацией можно ознакомиться в библиотеке факультета ВМиК МГУ. Автореферат разослан "... "....... 2004 г.

Ученый секретарь диссертационного Совета

профессор

Н.П. Трифонов

Общая характеристика работы

Актуальность работы.

Если система линейных уравнений

Ах = Ь

(1)

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

т — Ь — Ах.

(2)

Способ минимизации, преимущественно используемый на практике, — это минимизация в смысле наименьших квадратов. Иначе говоря, ищутся векторы ху реализующие минимум евклидовой длины (2-нормы) невязки (2). Такие векторы называются псевдорешениями системы (1); псевдорешение наименьшей евклидовой длины называется нормальным псевдорешением. Сам процесс решения системы (1) в указанном смысле называют линейной задачей наименьших квадратов, а еще точнее, безусловной задачей, имея в виду, что на поиск вектора х могут быть наложены дополнительные ограничения, приводящие к различным вариантам условной задачи наименьших квадратов.

С задачей наименьших квадратов тесно связано псевдообращение матрицы по Муру-Пенроузу: если — псевдообратная для матрицы системы (1), то нормальное псевдорешение этой системы дается формулой

ха

= А+Ь.

(3)

Существует и ряд других типов обобщенного обращения матриц, представляющих интерес для приложений, например, обращение квадратной выро-

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

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

Ах = В, (4)

где В — матрица размера пхт. Иначе говоря, становится возможным одновременное решение m систем типа (1) с одной и той же матрицей А, но различными правыми частями. Это означает, в частности, что возможно вычисление псевдообратной матрицы для чего достаточно положить в (4) В = 1п. Однако такой рекомендации мы не найдем в описании процедуры Least Squares. Это не случайно и свидетельствует о двух обстоятельствах: 1) разработчики системы недооценивают роль псевдообращения в линейной алгебре вообще и в задаче наименьших квадратов, в частности; 2) они не осознают того факта, что матрица А+ есть рациональная функция элементов матрицы А и, следовательно, может быть вычислена точно

(в режиме точных вычислений).

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

Цель работы.

В настоящей работе автор ставил перед собой следующие цели:

1. Систематизировать существующие методы для задачи псевдообращения по Муру-Пенроузу; выделить среди них методы, являющиеся рациональными алгоритмами в разъясненном выше смысле; реализовать такие алгоритмы процедурами языка Maple и с помощью численных экспериментов определить наиболее эффективные алгоритмы.

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

3. Найти рациональные алгоритмы для экономного пересчета решения задачи наименьших квадратов при ее одноранговых модификациях. Реали-

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

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

Научная новизна работы.

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

Модифицированы существующие и в ряде случаев сконструированы новые рациональные алгоритмы для решения перечисленных выше задач. Эти алгоритмы реализованы как процедуры языка Maple. Для отладки и сравнения эффективности этих процедур построены оригинальные тестовые задачи; некоторые из них представляют самостоятельный интерес для теории матриц.

На защиту выносятся следующие результаты работы:

1. Алгоритмы и Мар1е-процедуры для точного вычисления псевдообратной матрицы Мура-Пенроуза и точного решения безусловной задачи наименьших квадратов.

2. Алгоритмы и Maple-процедуры для точного решения условной задачи наименьших квадратов с ограничениями в виде линейных уравнений и/или неравенств.

3. Алгоритмы и Maple-процедуры для экономного пересчета решения задачи наименьших квадратов при ее одноранговых модификациях.

4. Алгоритмы и Maple-процедуры для точного вычисления обобщенной обратной матрицы Дрейзина.

Практическая значимость работы.

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

Апробация работы.

Основные результаты диссертации докладывались на следующих семинарах и конференциях:

Объединенный научно-исследовательский семинар кафедр автоматизации систем вычислительных комплексов, алгоритмических языков и системного программирования; научный руководитель — профессор М.Р. Шура-Бура; факультет ВМиК МГУ;

"Компьютерная алгебра"; научный руководитель — профессор С.А. Абрамов; ВЦ РАН.

33-я и 34-я конференции математиков Ирана.

9-й семинар иранских студентов и аспирантов, обучающихся в университетах Европы (Бирмингэм, июнь 2002 г.)

Публикации.

Основные результаты работы отражены в публикациях [1,4,5,7,8].

Структура и объем работы.

Диссертация состоит из введения,трех глав основного текста, заключения, двух приложений и списка литературы.

Объем диссертации — 107 страниц.

Библиография включает в себя 25 наименований.

Краткое содержание работы

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

Глава 1. Псевдообращение матрицы и безусловная задача наименьших квадратов. Указанные две задачи неслучайно объединены в названии данной главы: в силу (3), знание псевдообратной матрицы А+ позволяет немедленное вычисление нормального псевдорешения системы

(1).

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

ся рациональными. Это методы Гревилля и скелетного разложения, описанные в энциклопедии Ф.Р. Гантмахера "Теория матриц", и метод (названный нами методом Эрмита), предложенный в статье Rao T.M. et al.// SIAM J. Numer. Anal. 1976. V. 13. P. 155-171. К этому следует добавить модификацию метода скелетного разложения, указанную в книге Ben-Israel A., Greville T.N.E. Generalized inverses: Theory and applications. New York Wiley, 1974. Напомним, что скелетным разложением т х га-матрицы А, имеющей ранг г, называется всякое ее представление вида

где матрицы В и С имеют соответственно р а з м е 7риы >Епс л и разложение (5) известно, то классическая формула сводит псевдообращение матрицы А (в общем случае, прямоугольной) к обычному обращению двух квадратных невырожденных матриц и нескольким матричным умножениям:

А+ = С'{СС')-\ВтВ)-хВ\ (6)

Заметим, что над комплексным полем в список операций, разрешенных в рациональном алгоритме, включается взятие сопряженного числа. Формула, найденная в упомянутой выше книге "Generalized inverses", предусматривает лишь одно обращение:

(7)

В первых четырех разделах гл. 1 дается описание указанных выше алгоритмов, а также способа вычисления матрицы посредством библиотечной программы LeastSquares системы Maple. Вплоть до седьмой версии Maple, А+ можно было находить лишь путем обращения к LeastSquares с запросом на решение п задач типа (1) с одной и той же матрицей А и правыми частями, равными соответственно 1-му, 2-му,... ,п-му столбцам единичной матрицы

Алгоритмы Гревилля, Эрмита и скелетного разложения (в двух модификациях, соответствующих формулам (6) и (7)) были реализованы четырьмя Мар1е-процедурами. Особенности этих реализаций и выводы из численных экспериментов с ними обсуждаются в разделе 1.3. Главный из этих выводов состоит в том, что наиболее эффективным методом псевдообращения является алгоритм скелетного разложения, основанный на формуле (7). Поэтому в последующих главах диссертации при необходимости вычисления псевдообратной матрицы это делается, как правило, посредством вызова процедуры ЯР-1, реализующей этот алгоритм-

Глава 2. Задачи наименьших квадратов с линейными ограничениями. Первая часть этой главы посвящена задаче наименьших квадратов с ограничениями в виде линейных уравнений или, кратко, задаче НКУ. Задача формулируется так: среди всех векторов я, точно удовлетворяющих (совместной) линейной системе

найти вектор, минимизирующий величину

11^-/112- (9)

После краткого обсуждения постановки задачи НКУ в разделе 2.1.1 мы описываем рациональные алгоритмы для ее решения в разделе 2.1.2. Здесь уместно отметить, что в гл. 20-22 известной книги Лоусон Ч., Хенсон Р. Численное решение задач метода наименьших квадратов. М.: Наука, 1986 (в дальнейшем ссылаемся на нее как на книгу ЛХ) приведены три алгоритма, используемых для решения задачи НКУ в стандартных библиотеках численного анализа. Один из этих алгоритмов (см. гл. 22 названной книги), итерационный по своему характеру, мы сразу исключаем из рассмотрения. Два других алгоритма, которые в своем первоначальном они-

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

В разделе 2.1.3 обсуждаются тестовые задачи, используемые для экспериментов с алгоритмами решения задачи НКУ. Эти задачи используются к тому же для проверки ряда последующих алгоритмов; они представляют и самостоятельный интерес для теории матриц. Задачи состоят в вычислении расстояния, измеряемого посредством нормы Фробениуса (в русской литературе называемой обычно евклидовой матричной нормой), от заданной матрицы А до линейных многообразий, образованных соответственно обобщенными стохастическими матрицами, обобщенными дважды стохастическими матрицами и обобщенными магическими квадратами, а также вычислении наилучших аппроксимаций матрицы А в этих матричных классах. Достоинство этих задач в том, что, с одной стороны, они допускают интерпретацию как задачи НКУ, где вектор неизвестных составляют элементы искомой матрицы-проекции В на нужное линейное многообразие; с другой стороны, они имеют аналитическое решение, обосновываемое в данном разделе. Например, решение первой задачи дается формулой

где ^ — матрица порядка N все элементы которой равны числу 1 /N, а

В разделе 2.1.4, завершающем первую часть главы 2, описаны результаты экспериментов с двумя процедурами, реализующими построенные рациональные алгоритмы для задачи НКУ. Выясняется, в частности, что процедура LSE2, основанная на модифицированном алгоритме из гл. 21 книги ЛХ, всегда эффективней процедуры LSE1, реализующей (модифицированный) алгоритм из гл. 20 той же книги.

Во второй части главы 2 исследуется задача наименьших квадратов с ограничениями в виде линейных неравенств (или задача НКН), Задача формулируется так: среди всех векторов х, удовлетворяющих системе линейных неравенств

Gx > h, (10)

найти вектор, минимизирующий норму невязки (9). Неравенства между вещественными векторами одинаковой размерности понимаются как покомпонентные неравенства.

Обсуждение постановки задачи НКН в разделе 2.2.1 выявляет центральную роль специального случая этой задачи, называемого задачей NNLS (Nonnegative Least Squares — наименьшие квадраты с условием неотрицательности): минимизировать функцию (9) при условии

х>0. (И)

Для решения задачи NNLS имеется эффективный алгоритм, также называемый NNLS. В разделе 2.2.2 приведено его описание, воспроизводящее описание в [ЛХ, гл. 23, § 3]. В последующем обсуждении алгоритма отмечается, что наибольшие трудозатраты в нем связаны с решением безусловной задачи наименьших квадратов, повторяемым при каждом исполнении шага 6. При этом матрицы задач, отвечающих двум последовательным исполнениям данного шага, получаются одна из другой либо введением дополнительного ненулевого столбца, либо, наоборот, удалением нескольких (в типичном случае, одного) ненулевых столбцов. Эта ситуация, подробно разбираемая в главе 3, позволяет применить алгоритмы экономного пересчета нормальных псевдорешений при одноранговых модификациях задачи. Эти алгоритмы вводятся в тело основного алгоритма NNLS; результаты, полученные в численных экспериментах с модифицированным NNLS, обсуждаются в разделе 2.2.4. Эксперименты используют тестовые задачи,

описываемые в разделе 2.2.3. Отметим, что, как и задачи раздела 2.1.3, эти тестовые задачи имеют интересное матричное истолкование.

Еще одним важным частным случаем задачи НКН является задача LDP (Least Distance Programming): минимизировать Ц1Ц2 при условии (10). Стандартный метод численного решения задачи НКН, описанный, например, в [ЛХ, гл. 23], состоит в сведении ее к задаче LDP, которая, в свою очередь, сводится к задаче NNLS. При этом процедура, используемая для перехода

(12)

предусматривает ортогональное разложение матрицы Е (см. (9)) и, следовательно, не является рациональной.

В разделах 2.2.5 и 2.2.6 мы предлагаем рациональный алгоритм для перехода от задачи НКН к задаче NNLS. Алгоритм состоит из следующих основных шагов:

1. Правая часть исходной задачи (см. (9)) проектируется на образ матрицы Е. Искомая проекция /в определяется путем решения линейной системы, матрица которой есть матрица Грама системы столбцов матрицы Е.

2. Замена переменных

преобразует первоначальную задачу к виду: минимизировать \\у\\ при условиях

yeimE, GE+y>h-GE+fE. (13)

Эта задача отличается от стандартной задачи LDP наличием дополнительного условия

3. Используя LU-разложение матрицы Et находим базис ядра транспонированной матрицы Матрица составленная по строкам из векторов

найденного базиса, позволяет заменить условие у G im Е в (13) системой линейных уравнений, которой должен удовлетворять вектор у. В результате задача приобретает вид: минимизировать ]|у|| при условиях

4. Заменяем каждое уравнение

системы Sy = 0 двумя линейными неравенствами anxi H-----h SinXn > О,

После такой замены задача минимизации, полученная на шаге 3, приобретает вид стандартной задачи LDP.

5. Переход

LDP NNLS.

выполняем в соответствии со стандартным алгоритмом (см. [ЛХ, гл. 23, § 4]). Мы напоминаем этот алгоритм в разделе 2.2.6.

В заключительном разделе 2.2.7 обсуждается работа всего комплекса Maple-процедур, реализующих процесс решения задачи НКН. Примечательно, что в качестве тестовой здесь используется уже знакомая нам задача о наилучшей аппроксимации посредством обобщенной стохастической матрицы. При этом система линейных уравнений, описывающая данный матричный класс, заменяется эквивалентной системой линейных неравенств по указанному выше правилу. В результате матричная задача, которая в параграфе 2.1 интерпретировалась как задача НКУ, превращается в задачу с линейными ограничениями-неравенствами.

Глава 3. Пересчет псевдообратных матриц и нормальных псевдорешений при одноранговых модификациях задачи. В этой главе исследуется ситуация, типичная для многих приложений метода наименьших квадратов: после того как решена задача (1) CWX п-матрщей А, в систему вводится новое уравнение

От+1,1®1 Н-----Ь «m+l,nSn = Ът+\

(или, наоборот, из системы (1) удаляется одно из уравнений) и требуется решить новую задачу наименьших квадратов

лишь одним уравнением отличающуюся от задачи (1).

В подобной ситуации у пользователя Maple (как и других систем компьютерной алгебры) нет иного выбора, как решать задачу (14) "заново"без какого-либо использования информации о решении предыдущей задачи и с выполнением примерно того же объема арифметической работы. Между тем, существует возможность экономичного (и рационального, в принятом в данной диссертации смысле) пересчета нормального псевдорешения "ста-рой"задачи в решение "новой", требующая на порядок меньшей работы. Она основана на формулах пересчета псевдообратной матрицы при одноранговых модификациях матрицы Л, которые в практике псевдообращения играют ту же роль, какую для обычного обращения выполняет известная формула Шермана-Моррисона

Здесь cud — векторы-столбцы той же размерности п, что и матрица А.

В первой части главы 3 рассматривается безусловная задача наименьших квадратов. Сформулировав задачу об экономичном пересчете в раз-

деле 3.1.1, мы переходим к ее подробному обсуждению в следующем разделе. Сами по себе формулы пересчета псевдообратной матрицы не являются новыми: они заимствованы нами из гл. 3 книги Campbell S.L., Meyer CD. Generalized inverses of linear transformations. London: Pitman, 1979. Новой является их интерпретация как рациональных формул, позволяющих точное перевычисление псевдообратных матриц. Заметим, что в выбранную схему одноранговой модификации укладываются (после надлежащей трансформации) задачи о пересчете при введении дополнительной строки (столбца) или, наоборот, удалении строки (столбца), с которых выше началось обсуждение данной главы. Собственно говоря, классический метод псевдообращения, принадлежащий Гревиллю (см. гл. 1), представляет собой последовательность подобных пересчетов.

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

Во второй части главы 3 техника пересчета применяется к задаче НКУ. Обсуждение в этой части основано на недавней публикации Zhou J. et al. // SIAM J. Matrix Anal. Appl. 2002. V. 24, N 1. P. 150-164. В этой статье рассматривается решение рекурсивной задачи наименьших квадратов

(15)

где Е — матрица размера тхпи Em+i (/m+l) получается из Ет (fm) приписыванием дополнительной строки (компоненты). Авторы предлагают модификацию формулы Гревилля для пересчета псевдообратных матриц, позволяющую существенную экономию памяти при решении большой серии задач (15), т.е. при больших т. Что еще более важно, эта модифи-

кация позволяет указать способ экономного пересчета нормальных псевдорешений в рекурсивной задаче, где к целевой функции (15) добавляется постоянная система линейных связей (8).

После обсуждения задачи о пересчете в разделе 3.2.1 мы приводим в следующем разделе модифицированные формулы Гревилля из названной выше статьи и интерпретируем их как рациональные алгоритмы. В применении к задаче НКУ такой алгоритм может быть описан следующим образом. Пусть хт — нормальное псевдорешение задачи (15), (8). Положим го — С+й. Определим для каждого т две п х п-матрицы

Эти две матричные последовательности пересчитываются рекурсивно, исходя из начальных значений

Формулы их пересчета различны для двух возможных случаев: (%n+lQm =

Например, в

первом случае вначале вычисляется а затем матрица

Рт+1 — (/„ — ^т+^т+О^т-

Матрица (}т :в данном случае не изменяется: фт+1 = Если фт+1 — последняя компонента вектора то нормальное псевдорешение перевычисляется теперь по формуле

Ят+1 = Хт + (Фт+1 ~ ет+1хт)Ьт+1-

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

формул больше, чем для стандартных формул Гревилля. Более того, при непродуманной реализации новых формул число операций может возрасти на порядок. В разделе 3.2.3 обсуждается правильный способ реализации. Здесь же описаны тестовые задачи (в одной из них снова используются обобщенные стохастические матрицы) и приведены результаты численных экспериментов с Maple-процедурами, реализующими пересчет в рекурсивной задаче НКУ.

К основному тексту диссертации примыкает Приложение 1. В нем рассмотрены две задачи теории матриц, в решении которых используется операция псевдообращения. Первая го этих задач — обращение вырожденной матрицы по Дрейзину. Пусть А — пХn-матрица индекса т. Напомним, что индексом квадратной матрицы называется наименьшее натуральное число к такое, что

гапЫ* = гапкЛ*+1.

Обратная матрица Дрейзина X = А° матрицы А определяется как (единственное) решение системы матричных уравнений

Несмотря на присутствие в системе (16) квадратичного уравнения ХАХ — X, матрицу .A® можно вычислить рационально. Более того, рациональные алгоритмы (не рассматривавшиеся прежде с такой точки зрения) можно найти среди существующих. Мы обнаружили два таких алгоритма. Один из них описан в параграфе 7.2 цитировавшейся выше книги Кэмп-белла и Мейера; мы обсуждаем этот алгоритм в разделе П1.1.2. Другой в виде формулы

(р —любое натуральное число, не меньшее индекса т матрицы А) предложен в недавней статье Wang G. // Linear Algebra Appl. 2002. V. 348. P.

265-272. Поскольку доказательство формулы (17) в этой статье содержит ошибку, мы приводим в разделе Ш.1.3 собственное доказательство.

Составной частью обоих алгоритмов является вычисление индекса матрицы А. Используя идею одной из работ В.Н. Кублановской, мы показываем в разделе Ш.1.4, что нет необходимости в явном построении степепей А, как это обычно делается (например, в алгоритме Кэмпбелла-Мейера). Этот усовершенствованный алгоритм для вычисления А" мы называем методом последовательной редукции.

В разделе П1.1.5 мы обсуждаем построение тестовых задач для алгоритмов обращения по Дрейзину, а в разделе Ш.1.6 — результаты численных экспериментов с Мар1е-реализациями всех описанных выше алгоритмов. В целом эти результаты свидетельствуют о значительном превосходстве метода последовательной редукции над остальными методами.

В качестве второго приложения процедуры псевдообращения мы рассматриваем задачу о проверке матричного свойства

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

Вх > 0.

В предположении, что строки матрицы В линейно независимы, мы описываем в параграфе П1.2 алгоритм для проверки К- коположительности, составной частью которого является псевдообращепие некоторой матрицы. С Мар1е-реализации этого алгоритма и началась работа над диссертацией. Приложение 2 — это совокупность 18 таблиц, иллюстрирующих численные эксперименты, описанные в трех главах основного текста и Приложении 1.

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

Основные результаты работы

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

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

Публикации по теме диссертации

1. Икрамов Х.Д., Матин Фар М. Одноранговые модификации и пересчет псевдообратных матриц // Вести. МГУ. Сер. 15. Вычисл. матем. и кибернетика, 2002. N4. С. 12-17.

2. Matin far M. On computer algebra procedures for the pseudoinversion of matrices // 9th Iranian Students Seminar in Europe (Birmingham 2002).

3. Ikramov Kh.D., Matin fax M. On computer algebra procedures for the pseudoinversion of matrices // Absracts and General Guide. 33rd Iranian Mathematics Conference. Ferdowsi University of Mashhad. Mashhad, Iran,

2002, p. 87.

4. Икрамов Х.Д., Матин Фар М. О компьютерно-алгебраических процедурах для псевдообратных матриц // Ж . вычисл. матем. и матем. физ.

2003. Т. 43. N2. С. 163-168.

5. Икрамов Х.Д., Матин Фар М. Пересчет нормальных псевдорешений при одноранговых модификациях матрицы // Ж. вычисл. матем. и матем. физ. 2003. Т. 43. N4. С. 493-505.

6. Ikramov Kh.D., Matin far M. Rank one modifications and updating pseudoinverses // Absracts and General Guide. 34rd Iranian Mathematics Conference. Shahrood University of Technology. Shahrood, Iran, 2003, p. 64.

7. Икрамов Х.Д., Матин Фар M. О компьютерно-алгебраических проце дурах для линейной задачи наименьших квадратов с линейными связями // Ж . вычисл. матем. и матем. физ. 2004. Т. 44. N2. С. 206-212.

8. Икрамов Х.Д., Матин Фар М. О компьютерно-алгебраических процедурах для линейной задачи наименьших квадратов с ограничениями в виде линейных неравенств // Ж . вычисл. матем. и матем. физ. 2004. Т. 44. N5. С. 771-776.

9. Икрамов Х.Д., Матин Фар М. О компьютерно-алгебраической реализации метода наименьших квадратов на неотрицательном ортанте (принята к печати. Зап. научн. сем. ПОМИ. 2004. Т. 309.)

10. Икрамов Х.Д., Матии Фар М., Чесноков А.А. Компьютерно-алгебраические процедуры для обращения вырожденной матрицы по Дрей-зину (принята к печати. Ж . вычисл. матем. и матем. физ. 2004. Т. 44. N7.)

11. Икрамов Х.Д., Матин Фар М. Пересчет нормальных псевдорешений в рекурсивной задаче наименьших квадратов с линейными связями (принята к печати. Ж . вычисл. матем. и матем. физ. 2004. Т. 44. N10.)

Издательство ООО "МАКС Пресс". Лицензия ИД № 00510 от 01.12.99 г. Подписано к печати 20.04.2004 г. Формат 60x90 1/16. Усл.печл. 1,25. Тираж 100 экз. Заказ 448. Тел. 939-3890,939-3891,928-1042. Тел./факс 939-3891. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В.Ломоносова.

»13i«8

Введение.

Глава 1. Псевдообращение матрицы и безусловная задача наименьших квадратов.

1.1. Рациональные алгоритмы.

1.1.1. Скелетное разложение.

1.1.2. Метод Эрмита.

1.1.3. Метод Гревилля.

1.2. Вычисление псевдообратной матрицы в системе Maple

1.3. Результаты численных экспериментов.

Глава 2. Задачи наименьших квадратов с линейными ограничениями

2.1. Ограничения в виде линейных уравнений.

2.1.1. Постановка задачи HKY.

2.1.2. Рациональные алгоритмы.

2.1.3. Тестовые задачи.

2.1.4. Результаты численных экспериментов.

2.2. Ограничения в виде линейных неравенств.

2.2.1. Постановка задачи НКН.

2.2.2. Алгоритм NNLS.

2.2.3. Тестовые задачи.

2.2.4. Результаты численных экстриментов.

2.2.5. Преобразование задачи НКН в задачу LDP.

2.2.6. От задачи LDP к задаче NNLS.

2.2.7. Результаты численных экспериментов.

Глава 3. Пересчёт псевдообратных матриц и нормальных псевдорешений при одноранговых модификациях задачи.

3.1. Безусловная ЗНК.

3.1.1. Постановка задачи.

3.1.2. Формулы пересчёта.

3.1.3. Результаты численных экспериментов.

3.2. Задача НКУ.

3.2.1. Постановка задачи.

3.2.2. Модифиицрованные формулы Гревилля.

3.2.3. Детали реализации и численные эксперименты.

Напомним постановку линейной задачи наименьших квадратов (ЗНК): для заданных матрицы А размера (шхп)и га-мерного вектора Ь найти вектор х размерности п, минимизирующий величину в смысле наименьших квадратов или, более коротко, псевдорешением этой системы. Если система (0.2) совместна, то ее псевдорешения совпадают с решениями в обычном смысле.

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

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

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

1. Для матрицы А находим то или иное ортогональное (или унитарное) разложение, т.е. представление типа х) = \\АХ-Ъ\\2.

Всякий такой вектор х называется решением системы

Ах = Ь

0.1)

0.2)

А = НШ

0.3) где Н и К — унитарные матрицы соответственно порядка т и п, а Я — матрица вида

Я = Rn 0 х о о с невырожденной (г х г)-подматрицей Яц. 2. Представим вектор д = Н*Ь в виде

9 = \

91 v 92 , где д\ — подвектор размерности г. Обозначим через у\ единственное решение системы линейных уравнений

RnVi = 9i

Если г = ?2, то вектор х = Kyi есть единственное и, значит, нормальное псевдорешение системы (0.2). В противном случае, линейное многообразие псевдорешений описывается формулой х = К где у2 — произвольный вектор размерности п —г. Выбор у2 = 0 определяет нормальное псевдорешение х = К ~ \ Vi о

Два наиболее популярных типа разложения (0.3) — это С^Я-разложение матрицы Л и ее сингулярное разложение. Первый тип разложения может быть осуществлен разными способами — посредством той ли иной разновидности процесса ортогонализации Грама-Шмидта или с помощью элементарных унитарных матриц. В любом варианте это конечный процесс, использующий арифметические операции (к числу которых в случае комплексной системы прибавляется операция сопряжения) и квадратичные радикалы. Матрица К в QR-разложении — это некоторая матрица- перестановка, а подматрица Яц — треугольпая. Сингулярное разложение есть итерационный процесс, в котором вычисляется диагональная матрица Яц, составленная из сингулярных чисел матрицы А.

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

Сопоставим задаче (0.2) систему линейных уравнений

А* Ах = А*Ь. (0.4)

Она называется системой нормальных уравнений для задачи (0.2). Хорошо известно, что система нормальных уравнений всегда совместна (хотя исходная система (0.2) может быть и несовместной ), а ее решения суть псевдорешения задачи (0.2).

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

Если г < п, то система нормальных уравнений недостаточна для отыскания нормального псевдорешения. К ней необходимо присоединить какую-либо процедуру, выявляющую вектор минимальной длины в линейном многообразии образованном псевдорешениями; можно, например, найти проекцию на С нулевого вектора. Такое проектирование тоже может быть осуществлено посредством рационального алгоритма.

Вместо этого двухступенчатого процесса (построение многообразия С с помощью системы нормальных уравнений -Ь последующее проектирование нуля на С) можно предложить прямое описание нормального псевдорешения как функции от входных данных задачи (0.2). Это описание основывается на использовании аппарата псевдообратных матриц.

Понятие о псевдообращении как о частичной компенсации обычной операции обращения для случаев, когда матрица А — квадратная, но вырожденная, или прямоугольная, было введено в линейную алгебру в 1920-х годах американским математиком Е. Муром. Изложенные не лучшим образом, идеи Мура не находили отклика среди алгебраистов вплоть до 1950-х годов, когда в более ясной форме они были повторены английским статистиком Пенроузом. В частности, Пенроуз сформулировал систему уравнений (называемых теперь уравнениями Пенроуза), единственным решением которых является матрица, псевдообратная для матрицы А : А, (0.5)

ХАХ = X, (АХ)* = АХ, (ХА)* = ХА.

0.6) (0.7) (0.8)

Эта матрица получила символ А+ и нередко называется (псевдо)обратной матрицей Мура-Пенроуза, чтобы отличить ее от других типов обобщенных обратных матриц.

Несмотря на присутствие в системе уравнений Пенроуза квадратичного уравнения (0.6), псевдообратная матрица А+ может быть вычислена посредством рационального алгоритма. Этот удивительный факт более подробно обсуждается в первой главе диссертации. Возвращаясь к задаче (0.2), укажем, что если матрица А+ известна, то нормальное псевдорешение х можно найти по формуле

Эта формула и дает рациональное решение задачи вычисления нормального псевдорешения. В случае квадратной невырожденной матрицы А псевдообратная матрица А+ совпадает с обычной обратной матрицей А"1, а формула (0.9) переходит в хорошо известное равенство

Важная роль, которую приобрела в современной теории матриц операция псевдообращения, до недавнего времени не находила отражения в системах компьютерной алгебры. Так, в системе Maple вплоть до ее седьмой версии не было процедур для прямого вычисления матрицы А+. Если пользователю матрица А+ все же была необходима, то приходилось многократно обращаться к библиотечной процедуре LeastSquares, решающей задачу наименьших квадратов (0.2): задавая в качестве Ь координатные х = А+Ъ.

0.9) х = A~lb. векторы n-мерного пространства ei,., еп, пользователь получал в качестве нормальных псевдорешений столбцы псевдообратной матрицы.

Лишь в седьмой версии Maple процедура LeastSquares была модифицирована таким образом, чтобы можно было решать более общие задачи наименьших квадратов

Ах = В (0.10) с (тхр)-матрицей В. Теперь псевдообращение матрицы А может быть достигнуто единственным обращением к LeastSquares с единичной матрицей 1т в качестве правой части В задачи (0.10).

Однако и теперь возможности Maple даже в отношении собственно задачи наименьших квадратов остаются весьма ограниченными. Например, помимо безусловной ЗНК, обсуждавшейся до сих пор, на практике очень часто возникает необходимость в минимизации функционала (0.1) при тех или иных ограничениях на искомый вектор х, Наиболее распространенными типами ограничений являются системы линейных уравнений, системы линейных неравенств, комбинации тех и других и квадратичные ограничения вида ЦжЦг < с или ||ж — < с. Во второй главе диссертации показано, что задачи наименьших квадратов с линейными ограничениями могут быть решены рационально, и описаны соответствующие алгоритмы. Ни один из них не содержится в Maple.

Не эффективна Maple и в такой ситуации, типичной для многих приложений: после того как задача (0.2) решена, в систему вводится новое уравнение и требуется решить новую ЗНК

Ах = Ъ. (0.11)

Другая типичная ситуация: задача (0.11) получается из (0.2) удалением одного уравнения. В обоих этих случаях у пользователя Maple нет альтернативы решению задачи (0.11) как новой посредством той же процедуры

LeastSquares без какого-либо использования информации о решении предыдущей задачи.

В третьей главе диссертации обсуждаются методы эффективного пересчета нормального пссвдорешения при одноранговой модификации матрицы задачи (0.2):

А A + cd\

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

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

В каждой главе диссертации и Приложении 1 мы стараемся придерживаться следующего общего плана:

1. Постановка задачи.

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

2. Описание рациональных алгоритмов, решающих данную задачу.

Если "готовых"рациональных алгоритмов среди имеющихся методов нет, мы указываем необходимые модификации, превращающие эти методы в рациональные.

3. Обсуждение особенностей реализации алгоритмов из п.2 в виде процедур языка Maple.

4. Описание тестовых задач.

5. Обсуждение результатов тестирования (включающее в себя сравнение эффективности различных алгоритмов).

Вычисления проводились на персональном компьютере Pentium III-650 с версиями б и 7 системы Maple.

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

Разделы "Псевдообратные операторы и матрицы"и "Решение линейных систем по методу наименьших квадратов "были впервые введены в программу обязательных университетских курсов линейной алгебры в книге В.В. Воеводина "Линейная алгебра"(1974 г.) и сопровождавшем ее "Задачнике по линейной алгебре"Х.Д. Икрамова. Комплекс Мар1е-процедур, разработанных в данной диссертации, существенно облегчает учащемуся овладение этими разделами, избавляя от рутинных вычислений и позволяя сосредоточиться на сути.

Заключение

С теоретической точки зрения, основным результатом настоящей диссертации является доказательство возможности получения точных решений для двух важных классов линейно-алгебраических задач с точно заданными рациональными коэффициентами. Это, во-первых, обобщенное обращение матриц, среди многочисленных разновидностей которого наиболее полезными являются псевдообращение по Муру-Пенроузу (глава 1) и обращение вырожденных квадратных матриц по Дрейзину (см. П1.1). Это, во-вторых, линейные задачи наименьших квадратов. Возможность точного решения безусловной задачи наименьших квадратов (ЗНК) была осознана уже довольно давно, что отразилось, в частности, во включении процедуры LeastSquares уже в ранние версии Maple. Такого осознания не было в отношении ЗНК с ограничениями в виде линейных уравнений и/или неравенств. В главе 2 мы показали возможность точного решения этой задачи. То же относится к вопросам экономного пересчета псевдорешений при одноранговых модификациях задачи наименьших квадратов (глава 3).

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

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

Упомянем еще, что самостоятельный интерес для теории матриц представляет решение некоторых метрических матричных задач, использованных в главе 2 в качестве тестовых примеров.

1. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1966.

2. Ben-Israel A., Greville T.N.E. Generalized Inverses: Theory and Applications. New York: Wiley-Interscience, 1974.

3. Rao T.M., Subramanian K., Krishnamurthy E.V. Residue arithmetic algorithms for exact computation of д-inverses of matrices // SIAM J. Numer. Anal. 1976. V. 13. P. 155-171.

4. Грегори P.Т., Кришнамурти E.B. Безошибочные вычисления. Методы и приложения. М.: Мир, 1988.

5. Boullion T.L., Odell P.L. Generalized Inverse Matrices. New York: Wiley-Interscience, 1971.

6. Greville T.N.E. Some applications of the pseudoinverse of a matrix // SIAM Review. I960. V. 2. P. 15-22.

7. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. M.-JL: Физматгиз, 1963.

8. Икрамов Х.Д., Матин фар М. О компьютерно-алгебраических процедурах для псевдообращения матриц // ЖВМ и МФ. 2003. Т. 43, N2. С. 163-168.

9. Лоусон Ч., Хенсон Р. Численное решение задач метода наименьших квадратов. М.: Наука, 1986.

10. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.

11. Khoury R.N. Closest matrices in the space of generalized doubly stochastic matrices //J. Math. Anal. Appl. 1998. V. 222, N 2. P. 562-568.

12. Икрамов Х.Д., Матин фар M. Одноранговые модификации и пересчет псевдообратных матриц // Вести. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2002. N4. С. 12-17.

13. Икрамов Х.Д., Матин фар М. Пересчет нормальных псевдорешений при одноранговых модификациях матрицы // ЖВМ и МФ.2003. Т.43. N4. С.493-505.

14. Икрамов Х.Д., Матин фар М. О компьютерно-алгебраических процедурах для линейной задачи наименьших квадратов с линейными связями // ЖВМ и МФ. 2004. Т. 44. N2. С. 206-212

15. Икрамов Х.Д., Матин фар М. О компьютерно-алгебраической реализации метода наименьших квадратов на неотрицательном ортанте // Зап. научн. сем. ПОМИ. 2004. Т. 309.

16. Каханер Д., Молер К., Нэш С. Численные методы и программное обеспечение. М.: Мир, 1998.

17. Campbell S.L., Meyer C.D. Generalized inverses of linear transformations. London: Pitman, 1979.

18. Zhou J., Zhu Y., Li X.R., You Z. Variants of the Greville formula with applications to exact recursive least squares // SIAM J. Matrix Anal. Appl. 2002. V. 24, N 1. P. 150-164.

19. Wilkinson J.H. Note on the practical significance of the Drazin inverse. Nat. Physics Lab. Report, 1979, N 13.

20. Campbell S.L., Meyer C.D., Rose N.J. Applications of the Drazin inverse to linear systems of differential equations // SIAM J. Appl. Math. 1976. 31. P. 411-425.

21. Wang G. The reverse order law for the Drazin inverses of multiple matrix products // Linear Algebra Appl. 2002. V. 348. P. 265-272.

22. Икрамов Х.Д. Задачник по линейной алгебре. М.: Наука, 1975.

23. Кублановская В.Н. Об одном способе решения полной проблемы собственных значений вырожденной матрицы // Журнал вычисл. матем. и матем. физ. 1966. Т. 6. N 4. С. 611-620.

24. Икрамов Х.Д, Савельева Н.В. Коположительные матрицы // Журнал вычисл. матем.и матем. физ. 1999. Т. 39. N8. С. 1253-1279.

25. Ikramov Kh.D., Savel'eva N.V. Conditionally definite matrices // J.Math. Sci. 2000. V. 98. N1. P. 1-50.