Реферат: СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ В ЛИНЕЙНОЙ ЗАДАЧЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ
МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙФЕДЕРАЦИИ
Математическийфакультет
Кафедраприкладной математики
ДИПЛОМНЫЙ ПРОЕКТ
сингулярное разложение в линейной задаче методанаименьших квадратов
Заведующий кафедройприкладной
математики
Исполнил:
Научный руководитель
Владикавказ 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,так что приравнивание нулю производной не будет давать линейных уравнений дляотыскания оптимальных решений.
/>/>/>Глава 1. Метод наименьшихквадратов/>1.1. Задача наименьшихквадратовЗадача наименьших квадратов заключается в минимизация евклидовой длины вектора невязок || 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|| определяются единственным образом. Они независят от конкретного ортогонального разложения.
При проведении разложения необходимо приводить матрицы кдиагональному виду. Для этого обычно используются два преобразования: Гивенса иХаусхолдера, оставляющие нормы столбцов и строк матриц неизменными.
/>1.2. Ортогональноевращение ГивенсаЛемма. Пусть дан 2–вектор />, причем /> либо />.Существуетортогональная 2´2 матрица такая,что:
/> (8)
Доказательство. Положим:
/>.
Далее прямая проверка.
Матрица преобразования представляет собой матрицу вращений
/>
или отражений
/>
/>1.3. Ортогональноепреобразование ХаусхолдераПрименяется для преобразованияматриц к диагональному виду. Матрица преобразования представляет из себяследующее выражение: />, (9)
или, если вектор v нормирован,т.е. используется вектор единичной длины />, то />. В обоих случаях H –симметричная и ортогональная матрица. Покажем это:
/>.
Отсюда следует: что />, т.е. симметричность иортогональность. В комплексном случае матрица /> эрмитова[1] и унитарна[2].Предположим, что дан вектор х размерности m, тогда существуетматрица H такая, что />, где
/>
а s = +1,при положительной первой компоненте вектора х и = –1, при отрицательной.
Доказательство. Положим /> действительная матрица.Любую действительную матрицу можно привести в треугольному виду
/>
Далее принимаем во внимание то, что /> и получаем следующее:
/>
/>/>/>/>1.4.Сингулярное разложение матрицПусть 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 – матрица порядка Nxpc ортонормированными столбцами; Г – диагональная матрица порядка 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–мерном евклидовом пространстве.
/>1.5. QR–разложениеТеорема 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, то сингулярным разложением Абудет также и />, где />.
/>/>/>1.6. Число обусловленностиНекоторые вычислительные задачи поразительно чувствительны кизменению данных. Этот аспект численного анализа не зависит от плавающейарифметики или выбранного алгоритма.
Например:
Найти корни полинома: (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 – диагональная матрица, то />
/>глава 2. Реализациясингулярного разложения/>2.1. Алгоритмы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действительны и симметричны.
/>2.2. РеализацияразложенияТаким образом, разложение /> производится в два этапа.Сначала матрица А посредством двух конечных последовательностейпреобразований Хаусхолдера где />, приводится к верхнейдвухдиагональной форме следующего вида:
/>
Далее реализуется итерационный процесс приведениядвухдиагональной матрицы J0кдиагональной форме, так что имеет место следующая последовательность: /> где /> а Siи Ti – диагональные матрицы.
Матрицы Tiвыбираются так, чтобы последовательность матриц /> сходилась кдвухдиагональной матрице. Матрицы же Si выбирают так, чтобывсе Ji сохраняли двухдиагональную форму. Переход /> осуществляетсяс помощью плоских вращений (10) – преобразований Гивенса. Отсюда, /> где
/>
а матрица /> вычисляется аналогично сзаменой /> на/>.
Пусть начальный угол /> произволен, однакоследующие значения угла необходимо выбирать так, чтобы матрица Ji+1 имела ту же форму,что и Ji. Таким образом /> неаннулирует ни одного элемента матрицы, но добавляет элемент />; /> аннулирует/> нодобавляет />;/> аннулирует/> нодобавляет /> ит.д., наконец, /> аннулирует/> и ничегоне добавляет.
Этот процесс часто называют процессом преследования. Так как/>, то />, и Mi+1 – трехдиагональнаяматрица, точно так же, как и Mi.Начальный угол /> можновыбрать так, чтобы преобразование /> было QR–преобразованиемсо сдвигом, равным s.
Обычный QR–алгоритм со сдвигом можно записать вследующем виде:
/>
где /> –верхняя треугольная матрица. Следовательно, />. Параметр сдвига s определяется собственным значением нижнего минора (размерности 2´2) матрицы Mi.При таком выборе параметраs метод обладает глобальной и почти всегдакубичной сходимостью.
/>2.3. Пример сингулярногоразложенияПроведем преобразование Хаусхолдера на матрице />,
К первой компоненте первого столбца прибавляем норму первогостолбца, получим />. Пусть />
Преобразованная матрица 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 зависит работа, затрачиваемая на выполнение каждогопреобразования.
/>глава 3. Использованиесингулярного разложения в методе наименьших квадратовПри использовании метода сингулярного разложения (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Следует определить значение Y при X =1980.
Если аппроксимировать эти данные квадратичным многочленом /> и использовать двойнуюточность, то в результате получим следующие коэффициенты />. При одинарной точностивычислений коэффициенты будут иметь значения />.У этих двух наборов коэффициентов не совпадают даже знаки. Если такую модельиспользовать для предсказания, то для коэффициентов, вычисленных с двойнойточностью, прогноз будет Y=227780000, а дляобычной точности Y=145210000. Ясно, что второйнабор коэффициентов бесполезен. Исследуем достоверность результатов. Матрицаплана для данного примера имеет размеры 8 ´3
/>
Рис. 2. Численныйпример
/>
Ее сингулярные числа: />.
Число обусловленности равно />,что говорит о близости базисных функций 1, t и t2к линейной зависимости. Для того, чтобы исправить ситуацию можно предпринятьодну из следующих мер.
Во–первых, можно выбрать границу для относительнойошибки, которая бы отражала точность данных и точность арифметики. Если взятьграницу в интервале />, то отбросимтретье сингулярное число. При этом получим следующие наборы коэффициентов длядвойной и обычной точности:
/>
Теперь коэффициенты находятся в гораздо лучшем согласии другс другом. Кроме того, коэффициенты стали существенно меньше, а это значит, чтоне будет столь большого, как прежде, взаимного уничтожения слагаемых привычислении квадратичного многочлена. Прогнозное значение Y(1980)будет соответственно 212910000 и 214960000. Эффект обычной точности ещезаметен, однако результаты уже не являются катастрофическими.
Можно также определить набор нулевых коэффициентов,соответствующих пренебрежимо малому сингулярному числу. Вот эти коэффициенты: />. Для значений t от1900 до 1970 величина функции /> непревосходит 0.0017, поэтомупри любом a коэффициенты можноизменить />, и при этом значения,выдаваемые моделью изменятся не более чем на 0.0017a. Любой из четырех перечисленных нами наборов коэффициентовможно получить из другого подобным изменением.
Во–вторых, можно улучшить ситуацию заменой базиса.Модели
/>
гораздо более удовлетворительны. Важно при этом то, чтонезависимая переменная преобразуется из интервала [1900, 1970] в какой–нибудьболее приемлемый интервал вроде [0, 70] или, еще лучше, [–3.5, 3.5]. Числаобусловленности при этом равны 5750 и 10.7 соответственно. последнее значениеболее чем приемлемо даже при счете с обычной точностью.
Удобнее всего воспользоваться стандартными способамистатистического анализа, т.е. матрицу плана преобразуем к стандартизованномуварианту Матрица стандартизованных данных есть матрица наблюдений с нулевым средними дисперсией 1. Это означает, что данные берутся в виде отклонений от среднего,которое мы считаем равным 0, вводим нормировку деля каждый член столбца матрицына корень квадратный из суммы квадратов отклонений.
Во втором случае, после преобразования матрицы плана ееобусловленность сильно уменьшается, и, соответственно, повышается точностьрасчетов.
Данную программу можно использовать и при решении системылинейных уравнений вместо методов Гаусса, Жордана, Холесского и пр. Вприложении 2 приведен пример расчета линейной системы, которая изначально неможет быть решена этими методами вследствие вырожденности матрицыкоэффициентов. Тем не менее, исследуемый метод дает нам правильное решение.
/>/>/>ЗАКЛЮЧЕНИЕВ работе описаны компьютерные методы решения задачинаименьших квадратов. Для использования данных методов составленасоответствующая программа на алгоритмическом языке FORTRAN.Программа апробирована, результаты тестирования показывают работоспособностьпрограммы.
Результаты данной разработки могут быть использованы в самыхразнообразных расчетах, где необходимо провести аппроксимацию данных заданнымифункциями.
/>/>/>ЛИТЕРАТУРА1. Беллман Р. Введение в теорию матриц. -М.: Наука, 1969, 368с.
2. Гантмахер Ф.Р. Теория матриц. -М.: Наука, 1988, 548с.
3. Ланкастер П. Теория матриц. -М.: Наука, 1982, 387с.
4. Лоусон Ч., Хенсон Р. Численное решение задач наименьших квадратов. М.:Статистика, 1979, 447с
5. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980
6. Мэйндоналд Дж. Вычислительные алгоритмы в прикладной статистике. М.:Финансы и статистика, 1988, 350с
7. Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980, 454с
8. Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке АЛГОЛ. Линейнаяалгебра, М.: Машиностроение, 1976, 390с
9. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. -М.:Физматгиз, 1963, 536с.
10. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математическихвычислений. М.: Мир, 1980, 279с
11. Харебов К.С. Компьютерные методы решения задачи наименьших квадратов ипроблемы собственных значений. Владикавказ.: Изд-во СОГУ, 1995, 76 с.
/>/>/>ПРИЛОЖЕНИЕ 1. Исходные текстыпрограммыREAL A(3,3), U(3,3), V(3,3),SIGMA(3), WORK(3),Y(3),C(3),Y0(3)
INTEGER I,IERR, J, M, N, NM
OPEN(6,FILE=«SVD.OUT»,STATUS=«UNKNOWN»,FORM=«FORMATTED»)
OPEN (5,FILE= «SVD.IN»,STATUS=«UNKNOWN»,FORM=«FORMATTED»)
140 FORMAT(3I5)
150 FORMAT(4E15.7)
READ(5,140) NM,M,N
DO 131 I=1,M
READ(5,150) (A(I,J),J=1,N)
131 CONTINUE
READ (5,150) (Y(I),I=1,M)
CALLSVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,IERR,WORK)
IF(IERR.NE.0) WRITE (6,2) IERR
2 FORMAT(15H TROUBLE.IERR=,I4)
WRITE(6,120)
120 FORMAT(/'МАТРИЦА А')
DO 121 I=1,M
WRITE(6,130) (A(I,J),J=1,N)
130 FORMAT(8E15.7)
121 CONTINUE
WRITE (6,160) (Y(I),I=1,N)
160 FORMAT(/'ПРАВЫЕ ЧАСТИ'/8E15.7)
210 FORMAT(/'СИНГУЛЯРНЫЕ ЧИСЛА')
WRITE(6,210)
DO 3 J=1,N
WRITE(6,6) SIGMA(J)
3 CONTINUE
SMA=SIGMA(1)
SMI=SIGMA(1)
DO 211 J=2,N
IF(SIGMA(J).GT.SMA) SMA=SIGMA(J)
IF(SIGMA(J).LT.SMI.AND.SIGMA(J).GT.0.)SMI=SIGMA(J)
211 CONTINUE
OBU=SMA/SMI
230 FORMAT(/'ЧИСЛО ОБУСЛОВЛЕННОСТИ=',E15.7)
WRITE(6,230) OBU
SIGMA1=0.
DO 30 J=1,N
IF(SIGMA(J) .GT. SIGMA1) SIGMA1=SIGMA(J)
C(J)=0.
30 CONTINUE
TAU=SIGMA1*0.1E-6
DO 60 J=1,N
IF(SIGMA(J).LE.TAU) GO TO 60
S=0.
DO 40 I=1,N
S=S+U(I,J)*Y(I)
40 CONTINUE
S=S/SIGMA(J)
DO 50 I=1,N
C(I)=C(I) + S*V(I,J)
50 CONTINUE
60 CONTINUE
write (6,560)
WRITE (6,6) (C(I),I=1,3)
DO 322 J=1,N
SS=0.
DO 321 I=1,M
321 SS=A(J,I)*C(I)+SS
322 Y0(J)=SS
write (6,570)
WRITE (6,6) (Y0(I),I=1,3)
C WRITE(6,7)
C DO 4 I=1,M
C WRITE(6,6) (U(I,J),J=1,N)
C4 CONTINUE
C WRITE(6,7)
C DO 5 I=1,N
C WRITE(6,6) (V(I,J),J=1,N)
C5 CONTINUE
6 FORMAT(3E15.7)
560 format(2x,'roots')
570 format(2x,'right')
7 FORMAT(1H )
STOP
E N D
SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)
REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)
LOGICAL MATU,MATV
IERR=0
DO 100 I=1,M
DO 100 J=1,N
U(I,J)=A(I,J)
100 CONTINUE
G=0.0
SCALE=0.0
ANORM=0.0
DO 300 I=1,N
L=I+1
RV1(I)=SCALE*G
G=0.0
S=0.0
SCALE=0.0
IF(I.GT.M) GO TO 210
DO 120 K=I,M
120 SCALE=SCALE+ABS(U(K,I))
IF(SCALE.EQ.0.0) GO TO 210
DO 130 K=I,M
U(K,I)=U(K,I)/SCALE
S=S+U(K,I)**2
130 CONTINUE
F=U(I,I)
G=-SIGN(SQRT(S),F)
H=F*G-S
U(I,I)=F-G
IF(I.EQ.N) GO TO 190
DO 150 J=L,N
S=0.0
DO 140 K=I,M
140 S=S+U(K,I)*U(K,J)
F=S/H
DO 150 K=I,M
U(K,J)=U(K,J)+F*U(K,I)
150 CONTINUE
190 DO 200 K=I,M
200 U(K,I)=SCALE*U(K,I)
210 W(I)=SCALE*G
G=0.0
S=0.0
SCALE=0.0
IF(I.GT.M.OR.I.EQ.N) GO TO 290
DO 220 K=L,N
220 SCALE=SCALE+ABS(U(I,K))
IF(SCALE.EQ.0.0) GO TO 290
DO 230 K=L,N
U(I,K)=U(I,K)/SCALE
S=S+U(I,K)**2
230 CONTINUE
F=U(I,L)
G=-SIGN(SQRT(S),F)
H=F*G-S
U(I,L)=F-G
DO 240 K=L,N
240 RV1(K)=U(I,K)/H
IF(I.EQ.M) GO TO 270
DO 260 J=L,M
S=0.0
DO 250 K=L,N
250 S=S+U(J,K)*U(I,K)
DO 260 K=L,N
U(J,K)=U(J,K)+S*RV1(K)
260 CONTINUE
270 DO 280 K=L,N
280 U(I,K)=SCALE*U(I,K)
290 ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))
300 CONTINUE
IF(.NOT.MATV) GO TO 410
DO 400 II=1,N
I=N+1-II
IF(I.EQ.N) GO TO 390
IF(G.EQ.0.0) GO TO 360
DO 320 J=L,N
320 V(J,I)=(U(I,J)/U(I,L))/G
DO 350 J=L,N
S=0.0
DO 340 K=L,N
340 S=S+U(I,K)*V(K,J)
DO 350 K=L,N
V(K,J)=V(K,J)+S*V(K,I)
350 CONTINUE
360 DO 380 J=L,N
V(I,J)=0.0
V(J,I)=0.0
380 CONTINUE
390 V(I,I)=1.0
G=RV1(I)
L=I
400 CONTINUE
410 IF(.NOT.MATU) GO TO 510
MN=N
IF(M.LT.N) MN=M
DO 500 II=1,MN
I=MN+1-II
L=I+1
G=W(I)
IF(I.EQ.N) GO TO 430
DO 420 J=L,N
420 U(I,J)=0.0
430 IF(G.EQ.0.0) GO TO 475
IF(I.EQ.MN) GO TO 460
DO 450 J=L,N
S=0.0
DO 440 K=L,M
440 S=S+U(K,I)*U(K,J)
F=(S/U(I,I))/G
DO 450 K=I,M
U(K,J)=U(K,J)+F*U(K,I)
450 CONTINUE
460 DO 470 J=I,M
470 U(J,I)=U(J,I)/G
GO TO 490
475 DO 480 J=I,M
480 U(J,I)=0.0
490 U(I,I)=U(I,I)+1.0
500 CONTINUE
510 DO 700 KK=1,N
K1=N-KK
K=K1+1
ITS=0
520 DO 530 LL=1,K
L1=K-LL
L=L1+1
IF(ABS(RV1(L))+ANORM.EQ.ANORM) GOTO 565
IF(ABS(W(L1))+ANORM.EQ.ANORM) GOTO 540
530 CONTINUE
540 C=0.0
S=1.0
DO 560 I=L,K
F=S*RV1(I)
RV1(I)=C*RV1(I)
IF(ABS(F)+ANORM.EQ.ANORM) GO TO565
G=W(I)
H=SQRT(F*F+G*G)
W(I)=H
C=G/H
S=-F/H
IF(.NOT.MATU) GO TO 560
DO 550 J=1,M
Y=U(J,L1)
Z=U(J,I)
U(J,L1)=Y*C+Z*S
U(J,I)=-Y*S+Z*C
550 CONTINUE
560 CONTINUE
565 Z=W(K)
IF(L.EQ.K) GO TO 650
IF(ITS.EQ.30) GO TO 1000
ITS=ITS+1
X=W(L)
Y=W(K1)
G=RV1(K1)
H=RV1(K)
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
G=SQRT(F*F+1.0)
F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C=1.0
S=1.0
DO 600 I1=L,K1
I=I1+1
G=RV1(I)
Y=W(I)
H=S*G
G=C*G
Z=SQRT(F*F+H*H)
RV1(I1)=Z
C=F/Z
S=H/Z
F=X*C+G*S
G=-X*S+G*C
H=Y*S
Y=Y*C
IF(.NOT.MATV) GO TO 575
DO 570 J=1,N
X=V(J,I1)
Z=V(J,I)
V(J,I1)=X*C+Z*S
V(J,I)=-X*S+Z*C
570 CONTINUE
575 Z=SQRT(F*F+H*H)
W(I1)=Z
IF(Z.EQ.0.0) GO TO 580
C=F/Z
S=H/Z
580 F=C*G+S*Y
X=-S*G+C*Y
IF(.NOT.MATU) GO TO 600
DO 590 J=1,M
Y=U(J,I1)
Z=U(J,I)
U(J,I1)=Y*C+Z*S
U(J,I)=-Y*S+Z*C
590 CONTINUE
600 CONTINUE
RV1(L)=0.0
RV1(K)=F
W(K)=X
GO TO 520
650 IF(Z.GE.0.0) GO TO 700
W(K)=-Z
IF(.NOT.MATV) GO TO 700
DO 690 J=1,N
690 V(J,K)=-V(J,K)
700 CONTINUE
GO TO 1001
1000 IERR=K
1001 RETURN
E N D
/>ПРИЛОЖЕНИЕ 2. контрольныйпримерВходные данные/>
(матрица изначально сингулярна – первая строка равна сумме второй итретьей с обратным знаком)
3 3 3
.3200000E 02 .1400000E 02 .7400000E02
-0.2400000E 02 -0.1000000E 02 -0.5700000E02
-0.8000000E 01 -0.4000000E 01 -0.1700000E02
-0.1400000E 02 0.1300000E 02 0.1000000E01
Полученный результатМАТРИЦА А
.3200000E+02 .1400000E+02 .7400000E+02
-.2400000E+02 -.1000000E+02 -.5700000E+02
-.8000000E+01 -.4000000E+01 -.1700000E+02
ПРАВЫЕ ЧАСТИ
-.1400000E+02 .1300000E+02 .1000000E+01
СИНГУЛЯРНЫЕ ЧИСЛА
.1048255E+03
.7310871E-06
.1271749E+01
ЧИСЛО ОБУСЛОВЛЕННОСТИ= .1433830E+09
Корни
.1215394E+01 .1821742E+01 -.1059419E+01
Правые корни после проверки
-.1400000E+02 .1300000E+02 .1000001E+01
Видно, что правые частисоответствуют начальным данным. Решение верно.