5.11. Реализация метода Гаусса

Расчет направления задаваемого уравнением может быть выполнен несколькими способами. Здесь можно применить любой метод, пригодный для решения задач множественной линейной регрессии. Однако линейные задачи предъявляют в некотором смысле более жесткие требования к методам решения. Итоговое решение в них предполагается получить за один шаг, и потому матрица должна быть вычислена с высокой точностью. Для решения же нелинейных задач требуется несколько итераций. Небольшие ошибки на каждой итерации можно допустить при условии, что являются допустимыми выбираемые направления. Другими словами, в нелинейных задачах от матрицы в принципе требуется только положительная определенность. Следует отметить, однако, что значительные ошибки при расчете могут привести к резкому увеличению числа необходимых итераций.

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

Численные методы расчета вектора направления распадаются на два класса. Первый из них охватывает такие методы решения системы нормальных уравнений, которые не учитывают конкретной структуры этих уравнений. Очевидно, эти методы применимы независимо от того, являются ли регрессионные уравнения линейными по структуре или нет. Второй класс составляют методы, которые основываются на линейности регрессионной структуры показываются иногда неприменимыми при наличи» априорного распределения вероятности. Методы можно подразделить также на те, в которых

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

(а) Методы решения системы нормальных уравнений. В простейшем из них уравнение разрешается относительно с использованием обычных методов решения систем уравнений. Наиболее быстрым методом оказывается метод разложения Холецкого (см. раздел А.5), однако его рекомендуется применять лишь в тех случаях, когда известно, что положительно-определенная и хорошо обусловленная матрица. В общем случае мы рекомендуем употребить одну из разновидностей метода выбора направления (раздел 5.7) или метод Марквардта (раздел 5.8), каждый из которых обеспечивает улучшение обусловленности матрицы

Регрессионные методы. Сущность метода Дженриха и Сэмпсона [115] состоит в применении принципов шаговой регрессии (см. приложение А.3) к регрессионной задаче из раздела 5.10. Как и раньше, составляются нормальные уравнения, но компоненты вектора которые не могут привести к существенному уменьшению значения целевой функции, приравниваются нулю. Этот метод является методом выбора направлений, в котором направления совпадают с осями координат. Можно высказать предположение, что в большинстве нелинейных задач метод шаговой регрессии с последовательным включением и исключением окажется более эффективным, чем обычный регрессионный анализ.

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

Здесь мы упомянем два таких метода. В первом из них, методе Голуба [89], производится разложение (по методу Холецкого) матрицы непосредственно из при использовании преобразования Хаусхолдера. Другой метод, представляющий собой модифицированную ортогонализацию Грама — Шмидта, был предложен Лонгли в 1967 г. [138]. Он позволяет получить значительно более высокую точность, чем при непосредственном решении системы нормальных уравнений. В работе Голуба [90] указывалось, что этот метод несколько точнее процедуры Хаусхолдера.

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

Матрица В имеет строк и I столбцов, последние мы обозначим Матрица имеет размеры — это вектор-столбец с элементами. Система нормальных уравнений может теперь быть записана в матричном виде:

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

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

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

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

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

Используя выражения можно записать решение уравнения в виде

Вычисление особенно упрощается в том случае, когда А — легко обратимая матрица, например верхняя треугольная Матрица. Такую матрицу можно получить с помощью следующей процедуры, известной под названием модифицированной ортогонализации Грама — Шмидта:

1) образуется матрица Ниже будет описано преобразование этой матрицы, имеющей строк и столбцов. Обозначим через столбец матрицы

2) полагается ;

3) пусть ;

4) заменяется на заметим, что теперь нормированы в том смысле, что

5) пусть для ;

6) C заменяется на для Заметим, что тем самым векторы удается сделать ортогональными по отношению к без потери ранее установленного свойства их ортогональности к

7) если то процедура завершается, в противном случае заменяется на и переходят к операции 3.

Из примечаний, сделанных к операциям 4 и 6, ясно, что первые I столбцов итоговой матрицы С образуют ортонормированный (относительно базис для В и что последний столбец С содержит компоненты ортогональные к векторам в этом базисе. Другими словами, теперь мы имеем

Возвращаясь к нашему алгоритму, легко проверить, что

и

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

и

Из сопоставления уравнений следует:

Аналогично из уравнений следует, что

Таким образом, матрица А и вектор полностью определены и система уравнений теперь может быть разрешена относительно путем последовательных подстановок:

Может оказаться так, что векторы не будут линейно-независимыми. Допустим, ранг матрицы В (т. е. число линейно-независимых есть Тогда на итерациях процедуры ортогонализации окажется, что следовательно, и операции 4—6 выполнить не удастся. Наиболее простое решение в этом случае — оставить матрицу С неизменной и задать . Вследствие этого строк матрицы А и соответствующие элементы будут нулевыми, и, как это следует из равенства соответствующие окажутся неопределенными. Этим можно приписать произвольные значения, например, нулевые, а остальные значения вычислять с помощью

Ниже приводится простой числовой пример. Пусть

Выполняя последовательно операции в соответствии с описанным алгоритмом, получим:

4) заменяем первый столбец С на 10,056478, 0,169435, 0,423587,

6) второй и третий столбцы матрицы С заменяются на

4) второй столбец С заменяется на

6) третий столбец С заменяется на

Теперь мы имеем следовательно,

Легко убедиться, что этот вектор удовлетворяет нормальным уравнениям

В некоторых случаях нормальные уравнения имеют вид

где это заданная положительно-определенная матрица, заданный вектор.

Можно привести следующие примеры:

1. В методе Марквардта матрица состоит из диагональных элементов умноженных на скалярную величину X, а

2. Если в имеет априорное нормальное распределение вероятностей с ковариационной матрицей и средним то

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

Присоединяя к уравнения модели I дополнительных фиктивных уравнений, можно привести уравнения к нормальному регрессионному виду Пусть это такая матрица, что Определим следующие матрицы:

Легко проверить, что и Следовательно, нахождение линейной регрессии с заменой и на и эквивалентно решению уравнения