Содержание к диссертации
Введение
Глава 1. Взвешенный метод полных наименьших квадратов 12
1.1. Формулировка взвешенного TLS - метода 12
1.2. SVD - анализ WTLS - метода 16
1.3. Геометрическая интерпретация WTLS - задачи 19
1.4. Обусловленность вычислительных WTLS задач 22
1.5. Статистические свойства WTLS-метода 31
1.6. Выводы 38
Глава 2. Вычислительные аспекты взвешенного метода полных наименьших квадратов 39
2.1. Обзор численных алгоритмов метода полных наименьших квадратов 39
2.1.1. Численные методы и подходы решения задач полных наименьших квадратов 39
2.1.2. Методы решения плохо обусловленных СЛАУ 45
2.1.3. Прямой рекуррентный метод 50
2.1.4. Алгоритмы вычисления собственных значений 66
2.2. Метод расширенной системы уравнений 74
2.3. Выводы 78
Глава 3. Математическое моделирование процесса растворимости химических веществ 80
3.1. Формулировка задачи математического моделирования процесса растворимости 80
3.2. Результаты вычислений 83
3.3. Выводы 87
Заключение 88
Список литературы 89
- Геометрическая интерпретация WTLS - задачи
- Статистические свойства WTLS-метода
- Численные методы и подходы решения задач полных наименьших квадратов
- Алгоритмы вычисления собственных значений
Введение к работе
Общая характеристика работы
Актуальность темы. Почти любая задача математического моделирования, в которой исходных данных достаточно для того, чтобы переопределить решение, требует применения того или иного метода аппроксимаций. По ряду вполне объективных причин наиболее часто в качестве критерия аппроксимации выбирают метод наименьших квадратов. Основной причиной широкого использования для решения многих практических задач математического моделирования критерия наименьших квадратов является его "грубость" по отношению к априорным предположениям, используемым при решении конкретных практических задач [29]. Кроме того, задачи наименьших квадратов часто возникают и как составная часть некоторой более обширной вычислительной проблемы [22]. Например, определение орбиты космического корабля нередко сводится математиками к решению многоточечной краевой задачи для обыкновенного дифференциального уравнения. При этом вычисление орбитальных парамет- ров обычно требует нелинейного оценивания в смысле наименьших квадратов; в последнем используют различные схемы линеаризации. К задаче наименьших квадратов можно подойти как к задаче отыскания для заданной точки функционального пространства ближайшей точки в заданном подпространстве.
В диссертационной работе рассматривается взвешенный вари- . ——— —> ант метода полных наименьших квадратов, который позволяет решать большой класс задач линейного параметрического оценивания с учетом неоднородных ошибок в исходных данных [16].
Задача линейного параметрического оценивания является достаточно общей для широкого класса научных дисциплин, таких, как теория сигналов, автоматическое управление, теория систем, а также часто возникает при решении различных технических, статистических, физических, экономических, медицинских и других проблем.
Впервые метод полных наименьших квадратов был использован для решения задач регрессионного анализа в математической статистике академиком Ю.В. Линником [21]. В дальнейшем исследованию и применению метода полных наименьших квадратов посвящены работы российских ученых В.В. Федорова, А.И. Жданова [8, 10], А.В. Крянева [20], а также известных зарубежных ученых Дж. Голуба (G. Golub) [54 - 56], Ван Лоана (С. Van Loan), Ван Хаффеля (S. Van Huftel), Д. Вандевейла (J. Vandewalle), М. Кендалла, А. Стьюарта [16]. К сожалению, все эти работы основаны на предположении однородности ошибок в исходных данных, в то время как для больший- ства возникающих в практике задач это предположение не выполняется. Это хорошо известно даже на примере использовании обычного метода наименьших квадратов [22, 41]. Однако проблемы, возникающие при применении метода полных наименьших квадратов в условиях неоднородных ошибок в исходных данных существенно сложнее соответствующих проблем, возникающих в обычном взвешенном методе наименьших квадратов.
Поэтому актуальной на сегодняшний день является разработка эффективных численных алгоритмов для взвешенного варианта метода полных наименьших квадратов, предназначенного для решения задач математического моделирования в условиях неоднородности ошибок в исходных данных.
Цель диссертационной работы заключается в исследовании взвешенного варианта метода полных наименьших квадратов, позволяющего решать задачи математического моделирования и линейного параметрического оценивания в условиях неоднородных ошибок в исходных данных и разработке эффективных численных алгоритмов решения этой задачи.
Для достижения поставленной в работе цели были решены сле дующие задачи: /{/ d' "
1. Исследована обусловленность вычислительной WTLS - задачи.
2. На основе сингулярного анализа сформулированной задачи был получен критериальный вид для задачи взвешенных полных наи меньших квадратов.
Исследованы условия существования и единственности решения задачи взвешенных полных наименьших квадратов. , ,*.
Разработан эффективный численный алгоритм для решения WTLS - задачи.
Разработан метод математического моделирования процессов растворимости химических веществ, описываемых нелинейными алгебраическими уравнениями в неявном виде.
Научная новизна заключается в следующем.
Исследован критерий для WTLS - задачи в виде отношения двух положительно определенных квадратичных форм.
Получены условия существования и единственности рассматриваемой WTLS - задачи.
Исследована обусловленность вычислительной WTLS - задачи.
Для уменьшения числа обусловленности вычислительной задачи WTLS - метода предложено ее преобразование к эквивалентной задаче решения некоторой совместной системы линейных алгебраических уравнений (СЛАУ).
В предположении о стохастическом характере возмущений в исходных данных доказана сильная состоятельность оценок решений, получаемых предложенным WTLS - методом. Этот асимптотический (в статистическом смысле) результат подтверждает обоснованность предложенного WTLS - метода.
6. С применением WTLS - метода решена задача определения па- раметров математических моделей процессов, описываемых нелинейными алгебраическими уравнениями в неявном виде. Теоретическая и практическая значимость значимость работы состоит в том, что разработанный взвешенный вариант метода полных наименьших квадратов позволяет решать большой класс задач параметрического оценивания и математического моделирования в условиях сильно неоднородных (по точности) экспериментальных данных. Особую теоретическую и практическую значимость полученные результаты имеют для решения задач регрессионного анализа с ошибками в независимых переменных, играющего основополагающую роль в теории идентификации систем, эконометрике и многих других задачах, связанных с обработкой данных. Полученные в работе теоретические результаты применены для решения задачи идентификации параметров математической модели растворимости химических веществ.
Методы исследований. При формулировке и доказательстве результатов в диссертационной работе используются положения линейной алгебры, вычислительной линейной алгебры и современного численного анализа, а также теория параметрической идентификации систем. При разработке программного обеспечения использовался пакет MATLAB (версия 5.3).
Апробация работы. Результаты диссертации докладывались и обсуждались на:
1) Международной КОНфбрбНЦИИ "Математическое моделирование
ММ-2001" (г. Самара, 2001 г.);
Всероссийской конференции по прикладной и промышленной математике (г. Самара, 2001 г.);
12-й научной межвузовской конференции "Математическое моделирование и краевые задачи" (г. Самара, 2002 г.);
3-й Международной конференции молодых ученых и студентов "Актуальные проблемы современной науки" (г. Самара, 2002г.)
Личный вклад. Основные теоретические положения разработаны совместно с научным руководителем проф. А.И. Ждановым. Доказательство всех утверждений и теорем, исследование приложений, анализ результатов и выводы из них выполнены автором самостоятельно.
Публикации. По теме диссертации опубликовано 5 работ [16 -19, 30].
Структура и объем работы. Диссертация состоит из введения, трех глав, заключения, списка литературы из 89 наименований источников отечественных и зарубежных авторов. Общий объем диссертации составляет 100 страниц.
На защиту выносятся следующие положения:
1. Условия существования и единственности для взвешенного метода полных наименьших квадратов.
Стохастическая формулировка WTLS - метода и доказательство сильной состоятельности оценок решений, получаемых предложенным методом.
Анализ обусловленности вычислительной задачи взвешенных полных наименьших квадратов.
Преобразование WTLS - задачи к эквивалентной задаче, состоящей в решении некоторой совместной СЛАУ, что позволяет значительно уменьшить число обусловленности исходной вычислительной WTLS - задачи.
Алгоритм, основанный на прямом проекционном методе, позволяющий наиболее эффективно (по быстродействию) решать расширенную СЛАУ для WTLS - задачи.
Решение задачи определения оценок параметров математических моделей для процессов растворимости химических веществ, описываемых нелинейными алгебраическими уравнениями в неявной форме.
Краткое содержание работы
Введение содержит обоснование актуальности рассмотренных в диссертации вопросов. Здесь же определяются цель исследования, научная новизна и практическое значение. Кратко изложено содержание диссертации.
Первая глава содержит математическую постановку задачи и на основе сингулярного анализа получены условия существования и единственности WTLS - задачи. Приведена стохастическая формулировка этой задачи и на ее основе показана сильная состоятельность оценок решений.
В главе II приводится обзор наиболее известных алгоритмов решения задач полных наименьших квадратов. Исследованы недостатки при использовании этих алгоритмов для решения WTLS-задач. Получен новый численный алгоритм решения WTLS-задач на основе их преобразования к расширенной СЛАУ, который позволяет существенно повысить численную устойчивость рассматриваемых в диссертации задач.
Геометрическая интерпретация WTLS - задачи
Как видно из главы 1, задача полных наименьших квадратов имеет две эквивалентные формулировки. Из-за этого и численные методы ее решения разбиваются на две группы. Первые связаны с использованием итерационного процесса для вычисления сингулярного вектора и минимального сингулярного числа матрицы (А, 6) одновременно. Для этого чаще всего используют различные модификации градиентного метода Ньютона в явной или неявной форме [82]. Вторая группа методов основана на решении СЛАУ [10] при заранее найденном минимальном сингулярном числе и. К ней можно отнести хорошо известные метод Холецкого, методы QR-разложения [22], а также предлагаемый в настоящей работе алгоритм на основе прямо-то проекционного метода. Ньютоновские методы.
Одна из модификаций метода Ньютона может быть представлена в следующей форме. Из результатов главы 1 следует, что при Л = 0"min решение х задачи полных наименьших квадратов удовлетворяют следующему уравнению: Уравнение (2.1) определяет систему из п + 1 нелинейных уравнений относительно ж и Л. Одним из способов решения (см. [82]) является исключение х для получения уравнений для Л = cr in: Данный итерационный процесс будет иметь квадратичную скорость сходимости. Существенным недостатком данного метода является то, что сходимость Л к cr in и ж к XTLS обеспечивается только в том случае, если начальное приближение Л 0 удовлетворяет условию В общем случае такое приближение найти достаточно сложно. Модифицированный метод отношения Рэлея. В [82] предложена одна из наиболее удачных модификаций метода отношения Рэлея. Различные методы решения задачи можно получить с помощью применения Ньютоновского метода к следующей системе уравнений
В [82] отмечается, что данный метод соответствует методу обратных итераций со сдвигом, равным отношению Рэлея. Для симметричной проблемы собственных чисел данный метод имеет кубическую скорость сходимости [24] и является лучшим в классе стандартных Ньютоновских методов, применяемых к задаче (2.3). Для Итерационный процесс данного метода строится следующим способом. Если х - текущее приближение решения и рк - соответствующее ему отношение Рэлея, то следующее приближение х к+1 и масштабный множитель /Зк получаются из решения следующей симметричной СЛАУ: Так как J положительно определенная матрица, то решение может быть получено с помощью блочного исключения Гаусса, Модифицированный метод отношения Рэлея можно описать следующими уравнениями относительно разности векторов (2.3) Итерационный процесс метода задается уравнениями (2.4)-(2.6). В качестве начального приближения для х обычно выбирается решение системы (А, Ъ) в смысле обычных наименьших квадратов. Преимущество данного метода перед методом Ньютона в том, что процесс может сходится к минимальному сингулярному числу даже в том случае, когда ро = р(х ) ф. (c in mm)- Численные эксперименты показывают [82], что можно значительно улучшить начальное приближение, использовав несколько первых шагов метода обратных итераций. Критерием завершения итерационного процесса в [82] рекомендуется выбирать либо достижение неравенства где Со 0 - константа, подбираемая эвристическим путем [82], єь так называемый машинный нуль (т.е. максимальное положительное число, для которого выполняется равенство fl(l + Ємаш)"= 1 ) Через fl(a) обозначается значение величины, точно равной а, вычисленное в арифметике с плавающей точкой. Формально алгоритм метода может быть записан следующим образом
Статистические свойства WTLS-метода
В этом разделе предлагается прямой рекуррентный метод, который является альтернативной формулировкой прямого проекционного метода решения систем линейных алгебраических уравнений, и показывается его эквивалентность методу оптимального исключения. Это позволяет рассматривать последний как прямой рекуррентный метод. Предлагаемый подход дает потенциальную возможность создавать оптимальные по требуемому объему оперативной памяти машинные алгоритмы решения систем линейных алгебраических уравнений со многими правыми частями, которые особенно эффективно могут использоваться при решении расширенных систем линейных уравнений, которые будут получены ниже, в WTLS-задачах.
Рассматривается задача решения СЛАУ с невырожденными матрицами, т.е. задача вычисления решения СЛАУ
Прямые методы решения СЛАУ оцениваются по трем важнейшим качествам: числу выполняемых арифметических операций, требованию к объему оперативной памяти и максимальной точности, достижимой при их помощи.
Наилучшими, по выше перечисленным показателям, методами решения невырожденных СЛАУ продолжают оставаться методы типа исключения Гаусса. Как известно [27], число операций этих методов оценивается величиной 2п3/3 + 0(п2), где в число операций включаются как умножения/деления, так и сложения/вычитания.
Гауссовское исключение существует во многих вариантах, которые алгебраически тождественны. Из всего множества вариантов исключения Гаусса метод оптимального исключения требует для своей реализации минимальный объем оперативной памяти. Точно такие же характеристики имеет метод окаймления [5], хотя он и не относится, вообще говоря, к методам исключения, так как по своей математической структуре принадлежит к классу методов малоранговой модификации. Однако, метод окаймления представляет собой лишь некоторое видоизменение вычислительной схемы метода оптимального исключения и не дает никаких преимуществ по сравнению с классической вычислительной схемой метода оптимального исключения. Поэтому в дальнейшем будет рассматриваться лишь последний. В этом методе объем требуемой оперативной памяти оценивается при четном п величиной - п2/4 машинных слов и при нечетном п данный - показатель будет близок к этому.
Как известно [6], в методе оптимального исключения попеременно выполняются операции прямого и обратного хода метода Гаусса. Это позволяет не вводить в оперативное запоминающее устройство очередную строку СЛАУ до тех пор, пока не преобразованы предыдущие строки, т.е. используется последовательный ввод строк матрицы и правых частей системы. Таким образом, метод оптимального исключения имеет безусловно рекуррентную структуру.
В настоящей работе предлагается прямой рекуррентный метод решения СЛАУ, основанный на использовании формулы малоранговой модификации обратной матрицы (формулы Шермана - Моррисо-на [1]). Показана эквивалентность этого метода методу оптимального исключения, хотя по своей математической структуре полученный рекуррентный метод не является методом_ исключения. Этот метод относится к классу прямых проекционных методов [2].
Указанная эквивалентность рассматриваемого в статье прямого рекуррентного метода, прямого проекционного метода и метода оптимального исключения дает возможность исследовать последний как конечный рекуррентный процесс, а также представить его в компактной рекуррентно - матричной форме. Это, как будет ниже показано, позволяет получить оптимальные по требуемому объему оперативной памяти машинные алгоритмы решения СЛАУ со многими правыми частями, которые особенно эффективны в алгоритмах итерационного уточнения (см., например, [25]).
Как уже указывалось, предлагаемый в работе прямой рекуррент ный метод решения СЛАУ основан на формуле Шермана - Моррисона [1], которая формулируется следующим образом.
Численные методы и подходы решения задач полных наименьших квадратов
Большинство подходов, применяемых к решению задачи полных наименьших квадратов, в той или иной форме связаны с решением СЛАУ вида (2.7). Как уже отмечалось, желание эффективно решить такую СЛАУ требует выбора методов, использующих ее специфику. Специфика нормальной смещенной системы вида (2.7) состоит в том, что: - матрица системы является симметричной и положительно определенной; - эта матрица часто бывает плохо обусловленной; - система может быть сведена к эквивалентной СЛАУ с расширенной матрицей значительно большей размерности, но более численно устойчивой. Известно, что методы решения систем линейных алгебраических уравнений подразделяются на прямые и итерационные. К прямым методам решения СЛАУ относятся методы, дающие точное (без учета погрешностей округления результатов арифметических операций) решение за конечное число действий. Остальные методы относят к итерационным. В случае плохой обусловленности СЛАУ итерационные методы не могут обеспечить регулярность вычислительного процесса. Поэтому в данной работе рассматривались только прямые методы решения СЛАУ. Метод Холецкого. Данный метод [27] (в литературе он еще называется метод квадратного корня [27]) применяется к СЛАУ с симметричной положительно определенной матрицей, а следовательно, может быть применен и к смещенной нормальной системе вида (2.7). Метод Холецкого в применении к задаче полных наименьших квадратов основан на разложении симметричной положительно определенной матрицы на верхнюю LT Є IRpXp и нижнюю L треугольные матрицы. Этот метод является модификацией метода исключения Гаусса [27]. Формально алгоритм вычисления матрицы L может быть записан следующим образом: Существует также модификация этого алгоритма, при которой допускается перестановка строк матрицы, что позволяет на каждом шаге осуществлять выбор максимального диагонального элемента 1ц для повышения численной устойчивости.
После того, как вычислено разложение (2.8), решение системы вида (2.7) может быть получено решением двух систем с треугольными матрицами: Решение систем (2.9) представляет собой прямую и обратную подстановки (то есть решение СЛАУ с треугольными матрицами), что требует порядка 0{р2) операций с плавающей точкой. Решение с помощью метода Холецкого задачи ПНК требует в общей сложности з"Р3 + \v2 + 0(1р) операций с плавающей точкой. В работе [27] показано, что если симметричная положительно неопределенная матрица С Є Ш?Хр удовлетворяет условию то разложение Холецкого LT L может быть вычислено с использованием машинной арифметики и будет удовлетворять соотношению Метод Холецкого является одним из наиболее распространенных методов для решения нормальных СЛАУ. Он характеризуется простотой реализации алгоритма, низкой вычислительной трудоемкостью, существуют эффективные реализации метода для больших и разреженных матриц [41]. Одним из важнейших достоинств этого метода является также низкий уровень требований к оперативной памяти. Однако данный метод имеет и существенный недостаток - на плохо обусловленных матрицах метод Холецкого зачастую приводит к неприемлемой погрешности в решении. Известно [6], что прямые алгоритмы решения СЛАУ в общем случае оцениваются по трем параметрам: количество арифметических операций, объем оперативной памяти, необходимой для численного решения задачи, и численная обусловленность алгоритма (максимальная достижимая точность решения). Также при использовании того или иного метода на практике в качестве критерия можно-учитывать удобство и простоту его реализации. Однако в настоящее время при выборе метода решения СЛАУ все больше и больше на первый план выходит возможность построения параллельного вычислительного процесса на базе этого метода. Особенно это характерно, когда исследователю приходится сталкиваться с реализацией алгоритмов решения задач большой размерности на векторных (параллельных) компьютерах. Наилучшими по этим показателям методами продолжают оставаться методы типа исключения Гаусса. Всевозможные модификации гауссовского исключения существуют во множестве вариантов. Однако все эти варианты являются "алгебраически тождественными", т.е имеют одинаковую вычислительную сложность и структуру вычислений. В [2] показано, что все они принадлежат подклассу так называемых проекционных ABS-алгоритмов. Из всех вариантов исключения Гаусса метод оптимального исключения [1], [5] требует для своей реализации на ЭВМ минимальный объем оперативной памяти. Данный метод эквивалентен методу окаймления [5], который представляет собой некоторое видоизменение вычислительной схемы метода оптимального исключения.
В классе методов оптимального исключения [5] на каждом итерационном шаге попеременно выполняются операции прямого и обратного хода метода Гаусса. Это свойство позволяет распараллелить процесс вычисления решения, что является одним из важных преимуществ. Эффективной модификацией метода оптимального исключения является прямой рекуррентный метод [14], относящийся к классу прямых проекционных методов. Алгебраическая эквивалентность прямо го проекционного метода методу оптимального исключения доказана в [2]. Отличие прямого рекуррентного метода от остальных алгоритмов, принадлежащих к тому же классу, в том, что для каждого шага известна структура матриц и векторов используемых при итерационном процессе. К тому же, как показано в [2], прямой рекуррентный метод по вычислительным характеристикам является одним из наиболее эффективных в своем классе.
Алгоритмы вычисления собственных значений
Для решения задачи полных наименьших квадратов необходимо вычислять минимальные сингулярные числа (или соответствующие собственные числа) симметричных матриц. В процессе выбора алгоритмов нахождения минимальных собственных чисел матрицы учитывалось несколько особенностей конкретной задачи: 1) матрица (А, Ь) может быть плохо обусловленной; 2) необходимо вычислить только минимальное собственное число, без нахождения соответствующего собственного вектора.
Одно из наиболее важных преимуществ разработанного алгоритма в том, что нелинейная векторная проблема оптимизации (1.2) разбивается на две подзадачи: скалярную нелинейную задачу нахождения минимального сингулярного числа, являющегося коэффициентом смещения для смещенной нормальной системы, и задачу решения СЛАУ. Скалярная нелинейная задача минимизации решается значительно проще, чем векторная как с точки зрения вычислительной сложности, так и с точки зрения численной устойчивости и сходимости методов.
В разработанных численных методах решения полных задач наименьших квадратов (1.25) - (1.28) подразумевается, что минимальное сингулярное число матрицы (Л,;Ъ) должно быть вычислено до начала применения основного алгоритма. И при проведении вычислительных исследований необходимо решить вопрос о поиске оптимального для задачи ПНК метода нахождения минимального собственного числа симметричной положительной матрицы. Существует большое количество методов нахождения собственных чисел матрицы. Большинство из них являются градиентными, реализующими одновременно и поиск собственного вектора. Однако градиентные методы имеют ряд недостатков.
Так, к примеру, в [26] предлагается квазиньютоновский алгоритм, который получает искомое минимальное собственное значение используется метод сопряженных градиентов. Хотя метод имеет кубическую скорость сходимости, наличие неизвестных априори параметров, сильно варьирующихся в зависимости от конкретных значений матрицы, очень затрудняет его использование. Метод обратных итераций со сдвигом Рэлея [24] и модифицированный метод отношения Рэлея, описанный в разделе 2.1 не обладают свойствами глобальной сходимости. Общим недостатком градиентных методов является высокая вычислительная трудоемкость поиска решения.
Существуют и не градиентные методы поиска сингулярного числа, в частности, метод Гивенса, метод Хаусхолдера. Они не связаны с поиском собственных векторов и, следовательно, требуют намного меньше вычислительных операций и более численно устойчивы. Оптимальным для задачи полных наименьших квадратов, при проведении вычислительных экспериментов, показал себя метод Хаусхолдера. Этот метод заключается в приведении матрицы к трех-диагональному виду с последующим нахождением собственного числа методом деления пополам [23]. Данный метод и использовался в настоящей работе при вычислениях значений минимальных и максимальных сингулярных чисел и поэтому будет описан более подробно.
Метод Хаусхолдера. Метод основан на приведении симметрич ной положительно определенной матрицы к трехдиагональному виду с помощью ортогональных подобных преобразований. Такое приведение может быть выполнено неитерационными методами и, следовательно, требует намного меньше вычислений. Хотя метод Хаусхолдера и метод Гивенса, основанный на плоских вращениях, по многим вычислительным параметрам эквивалентны [24], предпочтение можно отдать первому в силу того, что он не требует вычисления тригонометрических функций sin И COS.
Пусть дана исходная симметричная матрица А Є IRnXn. Задача состоит в том, чтобы найти унитарные матрицы Рг, такие, что причем матрица А является треугольной. В методе Хаусхолдера в качестве матриц преобразования Рг выбраны матрицы отражения и приведение состоит из N = п — 2 шагов, на r-м шаге с помощью матриц отражения обнуляются элементы в r-й строке и r-м столбце, причем сохраняются нули, полученные на предыдущих шагах. Пусть трехдиагональная матрица, Ьг_і Є ffi n , Вт-\ Є jg n-r)x(n-r) _ неК0Т0рЫе матрицы, получающиеся в результате преобразований. Матрица преобразований на r-м шаге Рг может быть
