11.2.4. Ортогонально-треугольное разложение

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

или

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

где - верхняя треугольная матрица размера -матрица размера с ортонормированиями столбцами (т. е. Поскольку

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

Нормальные уравнения принимают здесь вид

или, поскольку матрица невырождена,

После того как найдена матрица

решить треугольную систему уравнений (11.16) не составляет труда. Вектор остатков равен

и из уравнения (11.8) вытекает, что

Другое разложение мы получим, если рассмотрим ортогональную матрицу

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

(иногда это соотношение записывают в виде который и дает этому методу название -алгоритм) и

Мы снова находим

и решаем уравнение Используя выражение (11.17), получаем остатки в виде

Правда, как замечено Джентлменом, такой способ вычисления остатков может оказаться неустойчивым в вычислительном отношении. Если матрица X запоминается (как это делается в случае итерационного уточнения в § 11.6), то лучше всего получать по формуле В то же время включение вектора в (11.9) оказывается полезным в двух отношениях. Во-первых,

и, во-вторых, поскольку то

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

Сущность процедуры, лежащей в основе соотношения (11.19), заключается в ортогональном преобразовании (т. е. вращении) матрицы Такое преобразование можно осуществить последовательно с помощью серии "элементарных" вращений, подобных приведенным ниже преобразованиям Хаусхольдера и Гивенса. При этом матрицу можно найти, если необходимо, применяя эти преобразования к единичной матрице или (в случае преобразований Хаусхольдера) с помощью запоминания и перемножения элементарных вращений. Матрицу проще всего найти по формуле (в разд. 11.2.3 мы уже видели, что обращение треугольной матрицы не составляет труда).

Теперь подробно рассмотрим численные методы отыскания приведенных выше матриц (или ).

(а) Процесс ортогонализации Грамма — Шмидта

Вычисления квадратных корней можно избежать, если использовать разложение

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

Используя процесс ортогонализации, мы можем, таким образом, преобразовать матрицу X к матрице с ненормированными столбцами и решить нормальные уравнения, используя метод, указанный в (11.12). Если вместо X использовать расширенную матрицу то получаемый при этом дополнительный элемент матрицы равен ее (см. упр. 6 в конце главы).

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

В результате [Jordan (1968)] получается расширенная матрица

так что при этом оценка наименьших квадратов для х и остаток вычисляются одновременно. Матрицу можно получить как

Для преобразования матрицы X к матрице или матрицы имеются два основных алгоритма: классический алгоритм Грама — Шмидта и модифицированный алгоритм Грама-Шмидта (МАГШ). Классический алгоритм для преобразования матрицы X с геометрической точки зрения представляет следующее: на шаге столбец делают ортогональным каждому из предварительно ортогонализованных столбцов, и такую процедуру выполняют последовательно для столбцов с номерами В модифицированном алгоритме на шаге делают ортогональным столбцу столбцы с номерами и эта процедура выполняется последовательно для Для более подробного ознакомления с этими алгоритмами можно рекомендовать работы Bjorck (1967а, Ь, 1968), Jordan (1968), Golub (1969) и в особенности Clayton (1971) и Farebrother (1974). Экспериментальное сравнение этих алгоритмов [Rice (1966), Jordan (1968), Wampler (1970, с. 556—557)], а также теоретический анализ [Bjorck (1967а)] указывают на то, что модифицированный алгоритм является более точным и более устойчивым, чем классический. Если матрица X совсем плохо обусловлена (§ 11.4), то при использовании вычисляемые столбцы матриц или быстро теряют ортогональность. Это означает, что не следует использовать без дополнительной переортогонализации; последняя же приводит к значительному возрастанию объема вычислений. В то же время при использовании никакой переортогонализации не требуется.

(b) Преобразования Хаусхольдера

Преобразованием Хаусхольдера называют всякую квадратную матрицу вида — где Здесь и так что матрица симметрична и ортогональна. Особенно нас интересуют преобразования вида

где Такое преобразование можно записать также в форме

где - преобразование Хаусхольдера размера Укажем кратко, как матрицу из (11.18) можно представить в виде

произведения преобразований Хаусхольдера, имеющих вид (11.24).

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

где - вектор размера элемент которого равен 1, а все остальные элементы равны нулю. Геометрически это представляет собой ортогональное преобразование, состоящее во вращении столбцов матрицы X до совпадения с "первой" координатной осью (представляемой вектором с). Следующий шаг состоит во вращении путем использования соответствующего преобразования остальных столбцов вокруг этой первой оси до тех пор, пока не попадет в плоскость, образуемую первой и второй осями. При этом становится линейной комбинацией с и Поскольку то , и поэтому

Этот процесс можно продолжить, и на шаге мы будем иметь

Окончательно

Матрица являясь произведением ортогональных матриц, также ортогональна.

Детали этого алгоритма приводят Golub (1965) и Businger, Golub (1965), которые представили также программу на языке АЛГОЛ-60. Эти авторы привели также две небольшие модификации этого алгоритма, позволяющие повысить эффективность вычислений. Дело в том, что, во-первых, каждое из преобразований Хаусхольдера требует вычисления двух квадратных корней (подробности см. в работе Golub, Styan (1973, с. 255)). Однако ерли записать в виде (это проекционная

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

где —матрица перестановок, переставляющая столбцы матрицы X в соответствующем порядке. Поскольку то нормальные уравнения после выполнения перестановки столбцов матрицы принимает вид

или, поскольку матрица невырождена,

Далее мы решаем уравнение относительно х и получаем

Что дает сравнение метода Хаусхольдера с процессом Грама — Шмидта? Выше уже упоминалось о том, что при применении столбцы матрицы могут оказаться далеко не ортогональными, особенно если матрица X плохо обусловлена. Если мы хотим получить матрицу, "сравнимую" в смысле "степени ортогональности" с матрицей получаемой преобразованием Хаусхольдера, то эти столбцы необходимо переортогонализовать [Wilkinson (1965, с. 244)]. Преимуществом по отношению к преобразованию Хаусхольдера является возможность относительно более простого программирования этого алгоритма. Кроме того, как представляется из экспериментов [Jordan (1968), Wampeer (1970)], МАГШ является несколько более точной процедурой, нежели преобразование Хаусхольдера. В то же время требует и несколько большего объема вычислений.

(с) Преобразования Гивенса

Еще одним типом ортогональных преобразований являются преобразования Гивенса, задаваемые матрицей

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

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

преобразовалась в пару

где

и

Если то и преобразование просто переставляет строки и изменяет знак:

Если то изменяется и знак первой строки в (11.30). Если то и никаких изменений не происходит.

Таким образом, последовательным применением преобразований Гивенса мы можем аннулировать (обратить в нуль) все элементы матрицы X, расположенные ниже диагонали, так что матрицу можно представить в виде

Такую редукцию можно выполнить двумя способами [Wilkinson (1965, с. 239—240)].

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

(2) На шаге преобразовать строку в парах с каждой из выше расположенных строк, начиная с первой, так, чтобы первые строк матрицы X приняли вид верхней треугольной матрицы.

Метод (1) требует запоминания в ЭВМ всей матрицы X, тогда как в методе (2) каждая из строк матрицы X нужна только в момент ее использования.

Преобразования (11.28) и (11.29) можно применить также и к матрице При этом мы получим векторы в (11.19). Например, если преобразованы строки с номерами то

После приведения матрицы к виду (11.19) мы можем при желании продолжить вращения пар элементов вектора и последовательно привести этот вектор к виду . В этом случае остаточная сумма квадратов ее равна просто и вектор получается умножением на первого столбца матрицы (ср. с выражением Мы можем всегда обеспечить положительность значения изменяя в случае необходимости его знак дополнительным "поворотом" на угол

Редукция матрицы X к верхней треугольной форме требует

примерно умножений и извлечений квадратного корня в случае преобразования Хаусхольдера и примерно умножений и пр извлечений квадратного корня в случае преобразования Гивенса. По этой причине предпочтение обычно отдают методу Хаусхольдера. В то же время преобразование Гивенса имеет два явных преимущества. Во-первых, редукцию матрицы X можно производить строка за строкой (метод (2), указанный выше). Во-вторых, уже имеющиеся в матрице X нули можно легко использовать для уменьшения количества арифметических операций. Как указывает Gentelman (1973) (его статья и служит основой этого раздела), эти преимущества важны по трем соображениям. Во-первых, в регрессионных и особенно в факторных планах размеры матрицы X могут быть столь значительными, что ее невозможно хранить в оперативной памяти, а значит, приходится как-то генерировать ее или вызывать из внешней памяти. Поскольку строк в матрице X обычно гораздо больше, чем столбцов, то построчная обработка уменьшает требуемый объем оперативной памяти и (или) количество обращений к внешним запоминающим устройствам. Кроме того, часто оказывается более естественным генерировать или вызывать матрицу X построчно. Например, в задачах регрессии каждая строка соответствует наблюдению над моделью и может возникнуть потребность в обновлении матрицы ( путем добавления дополнительных наблюдений (строк). В факторных планах также обычно оказывается более удобным образовывать именно строки, а не столбцы [Fowlkes (1969)]. Во-вторых, матрицы планов X часто содержат большое число нулей. Джентлмен сообщает, что применение преобразований Гивенса и простое использование этих нулей (не приспособленное специально к какой-нибудь специфической форме матрицы) приводили при анализе соответствующих несбалансированных планов к уменьшению вычислительных затрат вплоть до 70%. В-третьих, часто даже уже после того, как произведен анализ модели и матрицы редуцированы к матрицам и бывает необходимо привлечь дополнительные наблюдения или ввести те или иные ограничения (разд. 11.5.4). И здесь независимо от того, каким образом была произведена редукция—методом Хаусхольдера или методом Гивенса, мы можем произвести перетриангуляцию новой матрицы с помощью преобразований Гивенса, используя уже отмеченные их преимущества (§ 11.8).

Джентлмен предложил также модификацию указанного метода, подобную модификации, приведенной для метода Грама-Шмидта. Эта модификация позволяет избежать извлечения квадратных корней и сокращает на четверть необходимое число умножений. Она весьма успешно конкурирует с методом Хаусхольдера [Gentleman (1974а, Ь)].