Сингулярное разложение в линейной задаче метода наименьших квадратов
МИНИСТЕРСТВО ОБРАЗОВАНИЯ
РОССИЙСКОЙ ФЕДЕРАЦИИ
Математический факультет
Кафедра прикладной математики
ДИПЛОМНЫЙ
ПРОЕКТ
сингулярное разложение в линейной
задаче метода наименьших квадратов
Заведующий
кафедрой прикладной
математики
Исполнил:
Научный
руководитель
Владикавказ 2002
СОДЕРЖАНИЕ
ВВЕДЕНИЕ................................................................................................................................................................................ 3
Глава 1. Метод наименьших квадратов...................................................................................................... 7
1.1. Задача наименьших квадратов................................................................................................ 7
1.2. Ортогональное вращение Гивенса........................................................................................... 9
1.3. Ортогональное преобразование Хаусхолдера.............................................................. 10
1.4. Сингулярное разложение матриц........................................................................................... 11
1.5. QR–разложение...................................................................................................................................... 15
1.6. Число обусловленности................................................................................................................. 20
глава 2. Реализация сингулярного разложения.............................................................................. 25
2.1. Алгоритмы............................................................................................................................................... 25
2.2. Реализация разложения................................................................................................................ 27
2.3. Пример сингулярного разложения........................................................................................ 29
глава 3. Использование сингулярного разложения
в методе наименьших квадратов 33
ЗАКЛЮЧЕНИЕ..................................................................................................................................................................... 38
ЛИТЕРАТУРА....................................................................................................................................................................... 39
ПРИЛОЖЕНИЕ 1. Исходные тексты программы...................................................................................... 40
ПРИЛОЖЕНИЕ 2. контрольный пример........................................................................................................... 45
Метод наименьших квадратов обычно используется
как составная часть некоторой более общей проблемы. Например, при необходимости
проведения аппроксимации наиболее часто употребляется именно метод наименьших
квадратов. На этом подходе основаны: регрессионный анализ в статистике,
оценивание параметров в технике и т.д.
Большое количество реальных задач сводится к
линейной задаче наименьших квадратов, которую можно сформулировать следующим
образом.
Пусть даны действительная m´n–матрица A ранга k£min(m,n) и действительный m–вектор b. Найти действительный n–вектор
x0, минимизирующий евклидову длину вектора невязки Ax–b.
Пусть y – n–мерный
вектор фактических значений, x – n–мерный вектор
значений независимой переменной, b – коэффициенты в аппроксимации y
линейной комбинацией n заданных базисных функций j:
.
Задача состоит в том, чтобы в уравнении
подобрать такие b, чтобы минимизировать суммы квадратов отклонений e=y–Xb,
где X – есть так называемая матрица плана, в
которой строками являются n–мерный вектора с
компонентами, зависящими от xj: каждая строка
соответствует определенному значению xj.
Коэффициенты можно найти решая нормальные уравнения , откуда . Покажем это. Возведем в
квадрат выражение для е:
т. к. .
Это выражение имеет экстремум в точке, где =0
Откуда и получаем .
Следует отметить, что последнее выражение
имеет в определенной степени формальный характер, т. к. решение нормальных
уравнений, как правило, проводится без вычисления обратной матрицы (метод
Крамера) такими методами как метод Гаусса, Холесского и т. д.
Пример. Пусть заданы результаты четырех измерений (рис. 1): y=0 при x=0; y=1 при x=1; y=2 при x=3; y=5 при x=4. Задача заключается в том, чтобы провести через эти точки прямую таким образом, чтобы сумма
квадратов отклонений была минимальна. Запишем уравнение, описывающее проведение
прямой по
результатам измерений. Мы получаем переопределенную систему:
или Xb=y. Нам понадобится
матрица XTX и обратная к ней:
Тогда решение b=(XTX)-1XTy по методу
наименьших квадратов будет иметь вид
Таким образом, оптимальная прямая задается
уравнением Метод
точечной квадратичной аппроксимации (метод наименьших квадратов) не
предполагает, что мы должны приближать экспериментальные данные лишь с помощью
прямых линий. Во многих экспериментах связи могут быть нелинейными, и было бы
глупо искать для этих задач линейные соотношения. Пусть, например, мы работаем
с радиоактивным материалом. Тогда выходными данными у являются показания
счетчика Гейгера в различные моменты времени t. Пусть наш материал
представляет собой смесь двух радиоактивных веществ, и мы знаем период
полураспада каждого из них, но не знаем, в каких пропорциях эти вещества
смешаны. Если обозначить их количества через С и D, то показания
счетчика будут вести себя подобно сумме двух экспонент, а не как прямая:
. (1)
На практике, поскольку радиоактивность
измеряется дискретно и через различные промежутки времени, показания счетчика
не будут точно
Рис. 1. Аппроксимация прямой линией.
соответствовать (1). Вместо этого мы имеем серию
показаний счетчика в
различные моменты времени , и (1) выполняется лишь приближенно:
Если мы имеем более двух показаний, m>2,
то точно разрешить эту систему относительно C и D практически
невозможно. Но мы в состоянии получить приближенное решение в смысле минимальных
квадратов.
Ситуация будет совершенно иной, если нам
известны количества веществ C и D и нужно отыскать коэффициенты l и m. Это нелинейная
задача наименьших квадратов, и решить ее существенно труднее. Мы
по–прежнему будем минимизировать сумму квадратов ошибок, но сейчас она уже не
будет многочленом второй степени относительно l и m, так что приравнивание нулю
производной не будет давать линейных уравнений для отыскания оптимальных
решений.
Задача наименьших квадратов заключается
в минимизация евклидовой длины вектора невязок || Ax-b ||.
Теорема 1.
Пусть А – m´n–матрица ранга k, представленная в виде
A=HRKT (2)
где H ортогональная m´m матрица;
R – m´n–матрица вида
, (3)
где: R11 – kxk–матрица ранга k; K – ортогональная
kxk–матрица. Определим вектор
(4)
и введем новую переменную
. (5)
Определим как единственное решение
системы R11y1=g1. Тогда:
1.
Все решения задачи о минимизации ||Ax-b|| имеют вид , где y2 произвольно.
2.
Любой такой вектор приводит к одному и тому
же вектору невязки . (6)
3.
Для нормы r справедливо
4.
Единственным решением минимальной длины является
вектор
Доказательство. В
выражении для квадрата нормы невязки заменим A на HRKT в соответствии с (2) и умножая на ортогональную матрицу HT (умножение на ортогональную матрицу не меняет евклидову норму
вектора) получим
(7)
Далее из (3) и (5) следует, что
.
Из (4) следует
Подставляя оба последних выражения в (7)
получим
Последнее выражение имеет минимальное значение
при R11y1=g1, а в этом уравнении единственным решением
является , так
как ранг матрицы R11 равен к. Общее решение y выражается формулой , где y2 произвольно. Для вектора имеем
,
что устанавливает равенство (3). Среди
векторов наименьшую
длину имеет тот, для которого y2=0. Отсюда следует, что решением наименьшей длины будет вектор . Теорема доказана.
Всякое разложение матрицы А типа (2) мы
будем называть ортогональным разложением А. Заметим, что решение
минимальной длины, множество всех решений и минимальное значение для задачи
минимизации ||Ax-b|| определяются единственным образом. Они не зависят от конкретного
ортогонального разложения.
При проведении разложения необходимо приводить
матрицы к диагональному виду. Для этого обычно используются два преобразования:
Гивенса и Хаусхолдера, оставляющие нормы столбцов и строк матриц неизменными.
Лемма. Пусть дан
2–вектор , причем либо .Существует ортогональная 2´2 матрица такая, что:
(8)
Доказательство.
Положим:
.
Далее прямая проверка.
Матрица преобразования представляет собой
матрицу вращений
или отражений
Применяется
для преобразования матриц к диагональному виду. Матрица преобразования
представляет из себя следующее выражение: , (9)
или, если вектор v нормирован,
т.е. используется вектор единичной длины , то . В обоих случаях H –
симметричная и ортогональная матрица. Покажем это:
.
Отсюда следует: что , т.е. симметричность и
ортогональность. В комплексном случае матрица эрмитова[1] и унитарна[2]. Предположим,
что дан вектор х размерности m, тогда существует матрица H
такая, что , где
а s = +1, при положительной первой
компоненте вектора х и = –1, при отрицательной.
Доказательство.
Положим действительная
матрица. Любую действительную матрицу можно привести в треугольному виду
Далее принимаем во внимание то, что и получаем следующее:
Пусть X – матрица данных порядка Nxp, где N>p, и пусть r – ранг матрицы X. Чаще
всего r=p, но приводимый ниже результат охватывает общий случай, он
справедлив и при условии r<p.
Теорема о сингулярном разложении утверждает,
что
(10)
где V – матрица порядка Nxr,
столбцы которой ортонормированы, т.е. ; U – матрица с ортонормированными
столбцами порядка pxr; таким образом, ; Г – диагональная матрица порядка rxr,
диагональные элементы которой , называемые сингулярными числами матрицы X,
положительны. Используя диагональные элементы матрицы Г, столбцы матрицы V, и столбцы матрицы U, сингулярное
разложение матрицы X, определяемое по (10), можно записать в виде:
(11)
Имеют место следующие фундаментальные соотношения.
·
Квадратная симметричная матрица XX' порядка NxN,
имеет r положительных и N–r нулевых собственных чисел.
Положительными собственными числами XX' являются , а соответствующими собственными
значениями – .
Таким образом, сингулярные значения – это положительные квадратные корни из
положительных собственных чисел матрицы XX', а столбцы матрицы V
– соответствующие собственные векторы.
·
Квадратная симметричная матрица X'X порядка pxp,
имеет r положительных и p–r нулевых собственных чисел.
Положительными собственными числами X'X являются , а соответствующими собственными
значениями – ,
таким образом, сингулярные значения – это положительные квадратные корни из
положительных собственных чисел матрицы X'X, а столбцы матрицы U
– соответствующие собственные векторы.
Положительные собственные числа матрицы X'X
и XX' совпадают и равны . Более того, если um –
собственный вектор матрицы X'X, а vm – собственный вектор матрицы XX', соответствующие одному и тому
же собственному числу ,
то um и vm связаны
следующим соотношением
(12)
Эти соотношения дают возможность вычислять , зная , и наоборот. В компактной форме эти
соотношения можно записать следующим образом:
. (13)
Исследование матрицы X'X в факторном
анализе называется R-модификацией, а XX' – Q–модификацией.
Соотношения (12)–(13) показывают, что результаты Q–модификации можно
получить по результатам R–модификации и наоборот.
Практическая последовательность нахождения
сингулярного разложения следующая.
1.
Вычисляется X'X или XX', в
зависимости от того, порядок какой матрицы меньше. Предположим, что в данном
случае это X'X.
2.
Вычисляются положительные собственные числа матрицы X'X и
соответствующие им собственные векторы .
3.
Находятся сингулярные числа .
4.
Вычисляются по соотношению (11).
Пусть в разложении (11) собственные числа
расположены в порядке убывания. Аппроксимационные свойства соотношения (11)
являются еще более фундаментальными, чем само соотношение. Эти свойства
вытекают из решения следующих двух задач.
Задача 1. Дана симметричная матрица S, порядка pxp и ранга r
с неотрицательными собственными значениями. Требуется найти симметричную
матрицу Т, размерности pxp, с неотрицательными собственным
значениями заданного ранга k, k<r, являющуюся наилучшей
аппроксимацией матрицы S в смысле наименьших квадратов.
Задача 2. Дана прямоугольная матрица X, порядка Nxp и ранга
r и число k<r. требуется найти матрицу W порядка pxp
и ранга k, наилучшим образом аппроксимирующую матрицу X в смысле
наименьших квадратов.
Решением этих двух задач являются матрицы:
(14)
представляющие собой суммы k первых членов в
соответствующем разложении. Матрицы T и W называются наилучшими в
смысле наименьших квадратов “матричными аппроксимациями меньшего ранга” для
матриц S и X соответственно. Свойство наилучшей аппроксимации в
смысле наименьших квадратов можно выразить следующим образом: матрица T
ближе всего к матрице S в том смысле, что сумма квадратов всех элементов
матрицы S–T минимальна. Аналогично матрица W ближе всего к
матрице X в том смысле, что минимальна сумма квадратов элементов матрицы
X–W. Мерой близости или качества аппроксимации считается относительная
величина , т.е.
сумма r–k наименьших собственных чисел матрицы X’X. Иногда мерой
качества аппроксимации считается относительная величина
(15)
или функция от нее.
Рассмотрим наиболее распространенный случай p=r.
Матрица S может быть ковариационной
матрицей p линейно независимых переменных. Матрица T также может
представлять собой ковариационную матрицу p переменных, но так как ранг
матрицы T k<p, то эти p переменных линейно зависят от k
переменных. Таким образом, p исходных переменных, ковариационная матрица
которых есть S, могут быть приближенно выражены через k
переменных.
Во второй задаче исходную матрицу X
порядка Nxp можно выразить как X=VГU’, где V –
матрица порядка Nxp c ортонормированными столбцами; Г –
диагональная матрица порядка pxp, а U – квадратная
ортогональная матрица порядка pxp.
Матричную аппроксимацию меньшего ранга W
можно представить в виде
где – состоит из первых k столбцов матрицы
V, – из
первых k строк или столбцов матрицы Г, а – из первых k столбцов матрицы U.
поскольку W»X, то
(16)
При умножении этой матрицы справа на получаем
(17)
Матрица порядка pxk определяет преобразование
строк матрицы X из евклидова p–мерного пространства в евклидово k–мерное
пространство; уравнение (16) показывает, что существует преобразование матрицы X
порядка Nxp в матрицу порядка Nxk. Матрица X
содержит N точек в p–мерном евклидовом пространстве, которые
приближенно могут быть спроектированы в k–мерное евклидово пространство.
матрица определяет
координаты этих точек в k–мерном евклидовом пространстве.
Теорема 2. Пусть А
– m´n–матрица. Существует ортогональная m´m–матрица
Q такая, что в матрице QA=R под главной диагональю стоят только
нулевые элементы.
Доказательство.
Выберем ортогональную m´m–матрицу Q в соответствии с преобразованием
Хаусхолдера (9), так, чтобы первый столбец Q1A имел нулевые компоненты со 2–ой по m–ю. Далее выбираем
ортогональную (m-1)´(m–1)–матрицу P2 следующим образом. Будучи применена к m–1 вектору,
составленному из компонент со 2–ой по m–ю второго столбца матрицы Q1A, она аннулирует компоненты с 3–ей
по m–ю этого вектора. Матрица преобразования
ортогональна, и Q2Q1A имеет в первых двух столбцах нули под главной диагональю. Продолжая
таким образом, можно построить произведение, состоящее максимум из n
ортогональных преобразований, которое трансформирует А к верхней
треугольной форме. Формальное доказательство можно получить методом конечной
индукции.
Полученное представление матрицы произведением
ортогональной и верхней треугольной матриц называется QR–разложением.
Теорема 3. Пусть А
– m´n–матрица ранга к, причем k<n£m. Существуют ортогональная m´m–матрица
Q и m´n–матрица перестановки P такие, что
, (18)
где R – верхняя треугольная к´к–матрица
ранга к.
Доказательство.
Выберем матрицу перестановки Р таким образом, чтобы первые к
столбцов матрицы AP, были линейно независимы. Согласно теореме 2, найдется
ортогональная m´m–матрица Q такая, что QAP будет
верхней треугольной. Поскольку первые к столбцов АР линейно
независимы, это будет верно для первых к столбцов QAP.
Все элементы матрицы QAP, стоящие на
пересечении строк с номерами к+1,...,m и столбцов с номерами к+1,...,n,
будут нулями. В противном случае rankQAP>k, что противоречит
предположению rankA=k. Итак, QAP имеет форму, указанную правой
частью (4). Теорема доказана.
Подматрицу [R:T] из правой части
(18) можно теперь преобразовать к компактной форме, требуемой от матрицы R
из теоремы 2. Это преобразование описывает следующая лемма.
Лемма 1. Пусть [R:T]
– к´к–матрица, причем R имеет ранг к.
Существует ортогональная n´n–матрица W
такая, что
где – нижняя треугольная матрица ранга к.
Доказательство леммы вытекает из теоремы 3,
если отождествить величины n, k, [R:T], W из формулировки
леммы с соответствующими величинами m, n, AT, QT теоремы 3.
Используя теорему 3 и лемму 1 можно доказать
следующую теорему.
Теорема 4. Пусть А
– m´n–матрица ранга к . Найдутся
ортогональная m´m–матрица Н и ортогональная n´n–матрица К такие, что
(19)
причем R11 – невырожденная треугольная к´к–матрица.
Заметим, что выбором Н и К в
уравнении (19) можно добиться, чтобы R11 была верхней или нижней треугольной.
В (19) матрица А представлена
произведением A=HRKT, где R – некоторая прямоугольная матрица, ненулевые компоненты
которой сосредоточены в невырожденной треугольной подматрице. Далее мы покажем,
что эту невырожденную подматрицу R можно упростить далее до
невырожденной диагональной матрицы. Это разложение тесно связано со
спектральным разложением симметричных неотрицательно определенных матриц ATA
и AAT (см. 11).
Теорема 5.
Пусть А – m´n–матрица ранга k. Тогда существуют
ортогональная m´m–матрица U, ортогональная n´n–матрица V и диагональная m´n–матрица S такие, что
UTAV=S, A=USVT (20)
Матрицу S можно выбрать так, чтобы ее
диагональные элементы составляли невозрастающую последовательность; все эти
элементы неотрицательны и ровно к из них строго положительны.
Диагональные элементы S называются сингулярными
числами А. Докажем сперва лемму для специального случая m=n=rankA.
Лемма 2. Пусть А
– n´n–матрица ранга n. Тогда существует
ортогональная n´n–матрица U, ортогональная n´n–матрица V и диагональная n´n–матрица S такие, что UTAV=S, A=USVT
и последовательные диагональные элементы S положительны и
не возрастают.
Доказательство леммы. Положительно определенная симметричная матрица ATA допускает спектральное разложение
ATA=VDVT, (21)
где V – ортогональная n´n–матрица, а D – диагональная матрица, причем диагональные
элементы D положительны и не возрастают. Определим S как
диагональную n´n–матрицу, диагональные элементы которой суть
положительные квадратные корни из соответствующих диагональных элементов D.
Таким образом
D=STS=S2, S-1DS-1=I. (22)
Определим матрицу
U=AVS-1 (23)
Из (21), (22), (23) и ортогональности V следует,
что
UTU=S-1VTATAVS-1=S-1DS-1=I
т.е. U ортогональна.
Из (23) и ортогональности V выводим USVT=AVS-1SVT=AVVT=A Лемма доказана.
Доказательство теоремы 5. Пусть A=HRKT, где H, R, KT имеют свойства, указанные в теореме 4. Так как R11 из (19) – невырожденная треугольная к´к–матрица, то согласно лемме 2 , можно написать
(24)
Здесь и – ортогональные к´к–матрицы, а –
невырожденная диагональная матрица, диагональные элементы которой положительны
и не возрастают. Из (24) следует, что матрицу R в уравнении (19) можно
записать в виде
(25)
где:
– ортогональная m´m–матрица;
– ортогональная n´n–матрица;
– ортогональная m´n–матрица;
Теперь, определяя U и V
формулами
(26)
заключаем из (24) – (26), что A=USVT, где U, S, V
имеют свойства, указанные в формулировке теоремы 5. Это завершает
доказательство.
Заметим, что сингулярные числа матрицы А
определены однозначно, в то время, как в выборе ортогональных матриц U, V
есть произвол. Пусть s – сингулярное число А, имеющее кратность l. Это значит,
что для упорядоченных сингулярных чисел найдется индекс I такой, что
Положим k=min(m,n), и пусть Q –
ортогональная к´к–матрица вида
Здесь Р – ортогональная l´l–матрица Если A=USVT – сингулярное разложение А и si=…=si+l-1, то сингулярным разложением А
будет также и ,
где .
Некоторые вычислительные задачи поразительно
чувствительны к изменению данных. Этот аспект численного анализа не зависит от
плавающей арифметики или выбранного алгоритма.
Например:
Найти корни полинома: (x-2)2=10-6
Корни этого уравнения есть 2+10-3 и
2-10-3. Однако изменение свободного члена на 10-6 может
вызвать изменение в корнях, равное 10-3.
Операции с матрицами, как правило, приводят к
решению систем линейных уравнений. Коэффициенты матрицы в правой части системы
линейных уравнений редко известны точно. Некоторые системы возникают из
эксперимента, и тогда коэффициенты подвержены ошибкам наблюдения. Коэффициенты
других систем записываются формулами, что влечет за собой ошибки округлений. В
связи с этим необходимо знать, как влияют ошибки в коэффициентах матрицы на
решение. Именно для этого вводится понятие обусловленности матрицы.
По определению число обусловленности есть
величина . Для
более подробного описания числа обусловленности нам понадобится понятие нормы в
пространстве векторов и матриц.
Нормой вектора x в пространстве векторов называется функционал, обозначаемый , удовлетворяющий следующим
условиям:
1)
положительной определенности –
2)
положительной однородности – ;
3)
неравенству треугольника – .
Нормой квадратной матрицы А в пространстве матриц, согласованной с нормой вектора называется функционал , удовлетворяющий условиям
1 – 3 для нормы вектора:
1)
;
2)
3)
4)
мультипликативное неравенство –
Наиболее употребимы следующие нормы для
векторов:
·
норма суммы модулей
·
евклидова норма
·
норма максимума модуля
Нормы матриц:
·
·
·
Здесь являются сингулярными числами[3] матрицы А;
это положительные значения квадратных корней из собственных значений матрицы АТА
(которая при невырожденной матрице А положительно определена[4], в противном
случае положительно полуопределена (неотрицательно определена[5]) и поэтому
имеет только вещественные собственные значения ³ 0). Для вещественных
симметричных матриц сингулярные числа равны абсолютным величинам собственных
значений: .
Умножение вектора х на матрицу А
приводит к новому вектору Ах, норма которого может очень сильно
отличаться от нормы вектора х.
Область изменений может быть задана двумя
числами
Максимум и минимум берутся по всем ненулевым
векторам. Заметим, что если А вырождена, то m=0. Отношение M/m
называется числом обусловленности матрицы А,
(7)
Рассмотрим норму обратной[6] матрицы .
Для матрицы А существует сингулярное
разложение ,
тогда , отсюда . Аналогично для обратной
матрицы и . Отсюда следует, что
собственные числа матрицы – 1/ есть величины, обратные собственным числам
матрицы – . При этом очевидно, что . Из последнего выражения
вместе с (7) следует .
Таким образом обусловленность матрицы равна произведению нормы матрицы на норму
обратной матрицы.
Рассмотрим систему уравнений Ax=b, и
другую систему, полученную изменением правой части: A(x+Dx)=b+Db . Будем считать Db ошибкой в b, а Dx соответствующей ошибкой в x, хотя нам нет необходимости считать
ошибки малыми. Поскольку A(Dx)=Db, то определения M
и m немедленно приводят к неравенствам Следовательно , при m¹0,
Величина есть относительное изменение правой части, а
величина –
относительная ошибка, вызванная этим изменением. Аналогичные выкладки можно
провести не только с элементами вектора правой части но и с элементами самой
матрицы А и найти зависимость между относительным изменением элементов матрицы
и относительной ошибкой вызванной этим изменением. Отсюда следует, что число
обусловленности выполняет роль множителя в увеличении относительной ошибки.
Приведем некоторые свойства числа
обусловленности. Ясно, что M³m и поэтому cond(А)³1. Если Р
– матрица перестановок[7],
то компоненты вектора Px лишь порядком отличаются от компонент вектора х.
Отсюда следует, что и
cond(P)=1 . В частности cond(I)=1. Если А умножается на
скаляр с, то cond(cА)= cond(А). Если D – диагональная
матрица, то
QR–алгоритм
начинается с разложения матрицы по Грамму-Шмидту , затем меняются местами сомножители: Эта матрица подобна
первоначальной, Этот
процесс продолжается, причем собственные значения не изменяются:
Эта формула описывает QR–алгоритм без
сдвигов. Обычно время которое тратится на такой процесс пропорционально кубу
размерности матрицы – n3. Необходимо процесс ускорить, для чего используется предварительное
приведение матрицы А к форме Хессенберга[8] а
также используется алгоритм со сдвигом. Форма Хессенберга представляет из себя
верхнюю треугольную матрицу (верхняя форма Хессенберга) у которой сохранена
одна диагональ ниже главной, а элементы ниже этой диагонали равны нулю. Если
матрица симметрична, то легко видеть, что матрица Хессенберга превращается в трехдиагональную
матрицу[9].
При использовании матрицы Хессенберга время процесса пропорционально n2, а при использовании трехдиагональной
матрицы – n.
Можно использовать другие соотношения
где Qs – унитарная, а Ls – нижняя треугольная матрица. Такой
алгоритм носит название QL–алгоритма.
В общем случае, когда все собственные значения
матрицы различны, последовательность матриц As имеет пределом нижнюю треугольную матрицу , диагональные элементы которой
представляют собой собственные значения матрицы А, расположенные в
порядке возрастания их модулей. Если матрица А имеет кратные собственные
значения, то предельная матрица не является треугольной, а содержит
диагональные блоки порядка p, соответствующие собственному числу кратности p.
В общем случае, наддиагональный элемент матрицы As на s-ом шаге асимптотически
равен , где kij – постоянная величина. Сходимость QL–алгоритма вообще говоря
недостаточна. Сходимость можно улучшить, если на каждом шаге вместо матрицы As использовать матрицу As-ksI (QL–алгоритм со сдвигом). Последовательность вычислений в этом
случае описывается следующими соотношениями:
которые определяют матрицу . При этом асимптотическое поведение
элемента определено
соотношением , а
не , как прежде.
Если сдвиг ks выбрать близко к величине (наименьшее собственное
значение), то в пределе внедиагональные элементы первой строки будут очень
быстро стремиться к нулю. Когда ими можно пренебречь, элемент с рабочей точностью равен , остальные являются
собственными значениями оставшейся матрицы n-1-го порядка. Тогда,
если QL–алгоритм выполнен без ускорения сходимости, то все равно , и поэтому автоматически
можно выделить величину сдвига ks.
Если матрица А эрмитова, то очевидно,
что и все матрицы Аs эрмитовы;
если А действительная и симметричная, то все Qs ортогональны и все Аs
действительны и симметричны.
Таким образом, разложение производится в два
этапа. Сначала матрица А посредством двух конечных последовательностей
преобразований Хаусхолдера где , приводится к верхней двухдиагональной форме следующего вида:
Далее реализуется итерационный процесс
приведения двухдиагональной матрицы J0 к диагональной форме, так что имеет место следующая
последовательность: где а Si и Ti – диагональные матрицы.
Матрицы Ti выбираются так, чтобы последовательность матриц сходилась к
двухдиагональной матрице. Матрицы же Si выбирают так, чтобы
все Ji сохраняли двухдиагональную форму. Переход осуществляется с помощью
плоских вращений (10) – преобразований Гивенса. Отсюда, где
а матрица вычисляется аналогично с заменой на .
Пусть начальный угол произволен, однако
следующие значения угла необходимо выбирать так, чтобы матрица Ji+1 имела
ту же форму, что и Ji. Таким образом не аннулирует ни одного
элемента матрицы, но добавляет элемент ; аннулирует но добавляет ; аннулирует но добавляет и т.д., наконец, аннулирует и ничего не добавляет.
Этот процесс часто называют процессом
преследования. Так как , то , и Mi+1 – трехдиагональная матрица, точно так же,
как и Mi. Начальный угол можно выбрать так, чтобы
преобразование было QR–преобразованием со сдвигом, равным s.
Обычный QR–алгоритм со сдвигом можно
записать в следующем виде:
где – верхняя треугольная матрица. Следовательно, . Параметр сдвига s
определяется собственным значением нижнего минора (размерности 2´2) матрицы Mi. При таком выборе параметра s метод обладает глобальной и почти
всегда кубичной сходимостью.
Проведем преобразование Хаусхолдера на матрице
,
К первой компоненте первого столбца прибавляем
норму первого столбца, получим . Пусть
Преобразованная матрица A2 вычисляется следующим образом. Для первого
столбца имеем:
так как
Таким образом, в первый столбец были введены
нули и его длина не изменилась. Получим второй столбец:
для третьего столбца:
окончательно,
Столбцы матрицы A2 получаются вычитанием кратных вектора v1 из столбцов A1. Эти кратные порождаются скалярными
произведениями, а не отдельными элементами, как в гауссовом исключении.
Прежде чем вводить дальнейшие нули под диагональю,
преобразованием вида A3=A2Q1, получают нули в первой строке. Нули уже
стоящие в первом столбце, не должны быть испорчены, длина первого столбца
должна быть сохранена; поэтому при внесении нулей в первую строку нельзя менять
первый элемент строки, изменяем второй элемент и зануляем третий. Для матрицы
большего размера на этом шаге было бы получено n–2 нуля.
Преобразование порождается первой строкой A2:
Строка матрицы A3 с номером i получается по формуле
.
Таким образом, из каждой строки A2 вычитается надлежащее кратное . Это дает матрицу
Поскольку первая компонента нулевая, то нули первого
столбца A2
сохраняются в A3,
Так как Q1 ортогональная,
то длина каждой строки в A3 равна длине соответствующей строки в A2.
Теперь можно добиться новых нулей под
диагональю, не испортив полученных ранее:
Поскольку ранг этой матрицы равен лишь 2, то
теперь третий столбец имеет на диагонали и под диагональю элементы порядка
ошибки округления. Эти элементы обозначены в матрице через 0.000, чтобы
отличить их от элементов, в точности равных нулю. Если бы матрица имела полный
ранг, то нужно было бы выполнить еще одно преобразование, чтобы получить нули в
третьем столбце:
Если бы не ошибки округлений, то в данном
примере третий диагональный элемент был бы точным нулем. Элементы под
диагональю во всех столбцах указаны как точные нули, потому что преобразования
так и строились, чтобы получить там нули. Последнее преобразование H3 в этом примере могло бы быть
тождественным, однако тогда оно не было бы хаусхолдеровым отражением.
Фактически использование H3 попутно изменяет знак элемента – 1.080 в матрице A4.
Получилась искомая двухдиагональная матрица, и
первый этап закончен. Прямое использование ортогональных преобразований не
позволяет получить какие–либо новые нули. Для общего порядка n нужно n
преобразований H и n–2 преобразований Q, чтобы достигнуть
данного места. Число преобразований не зависит от строчной размерности m,
но от m зависит работа, затрачиваемая на выполнение каждого
преобразования.
При использовании метода сингулярного разложения
(SVD – Singular Value Decomposition) мы проводим
разложение для матрицы плана. При этом основное уравнение y=Xb приобретает вид y=USVTb. Отсюда следует, что коэффициенты b можно получить решая
уравнение UTy=SVTb.
Т.е. если все sj, j=1,…,n,
являющиеся диагональными элементами S отличны от нуля, то последнее уравнение
разрешимо и
, где .
Однако такое решение не всегда желательно,
если некоторые sj малы. Для правильного использования метода SVD мы должны ввести
границу t отражающую точность входных данных и точность использованной плавающей
арифметики. Всякое sj, большее, чем t, приемлемо, и соответствующее вычисляется по (1.20). Любое sj, меньшее, чем t, рассматривается как пренебрежимо малое, и соответствующему может быть придано
произвольное значение. С этой произвольностью значений связана не единственность
набора коэффициентов, получаемых методом наименьших квадратов. Изменения
входных данных и ошибки округлений, меньшие, чем t, могут привести к совершенно другому набору коэффициентов,
определяемых этим методом. Поскольку обычно желательно, чтобы эти коэффициенты
были по возможности малы, то полагаем =0, если sj £t.
Отбрасывание чисел sj,
меньших, чем t, приводит к уменьшению числа обусловленности с до . Поскольку число обусловленности является
множителем в увеличении ошибки, то следствием будет более надежное определение
коэффициентов .
Продемонстрируем использование метода на
следующем примере:
t
|
Y
|
1900
|
75994575
|
1910
|
91972266
|
1920
|
105710620
|
1930
|
123203000
|
1940
|
131669275
|
1950
|
150697361
|
1960
|
179323175
|
1970
|
203211926
|