Сингулярное разложение в линейной задаче метода наименьших квадратов

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ

Математический факультет

Кафедра прикладной математики

ДИПЛОМНЫЙ ПРОЕКТ

сингулярное разложение в линейной задаче метода наименьших квадратов

Заведующий кафедрой прикладной

математики          

Исполнил:                                                           

Научный руководитель

Владикавказ 2002

СОДЕРЖАНИЕ

 TOC \o "1-2" ............................................................................................................................................................................. \h

Глава 1. Метод наименьших квадратов.................................................................................................. \h

1.1. Задача наименьших квадратов......................................................................................................... PAGEREF _Toc11128509 \h 7

1.2. Ортогональное вращение Гивенса................................................................................................... PAGEREF _Toc11128510 \h 9

1.3. Ортогональное преобразование Хаусхолдера.......................................................................... PAGEREF _Toc11128511 \h 10

1.4. Сингулярное разложение матриц................................................................................................... PAGEREF _Toc11128512 \h 11

1.5. QR–разложение........................................................................................................................................ PAGEREF _Toc11128513 \h 15

1.6. Число обусловленности....................................................................................................................... PAGEREF _Toc11128514 \h 20

глава 2. Реализация сингулярного разложения.......................................................................... \h

2.1. Алгоритмы.................................................................................................................................................. PAGEREF _Toc11128516 \h 25

2.2. Реализация разложения....................................................................................................................... PAGEREF _Toc11128517 \h 27

2.3. Пример сингулярного разложения.................................................................................................. PAGEREF _Toc11128518 \h 29

глава 3. Использование сингулярного разложения  в методе наименьших квадратов              \h

ЗАКЛЮЧЕНИЕ................................................................................................................................................................... \h

ЛИТЕРАТУРА..................................................................................................................................................................... \h

ПРИЛОЖЕНИЕ 1. Исходные тексты программы............................................................................... \h

ПРИЛОЖЕНИЕ 2. контрольный пример..................................................................................................... \h

ВВЕДЕНИЕ

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

Большое количество реальных задач сводится к линейной задаче наименьших квадратов, которую можно сформулировать следующим образом.

Пусть даны действительная m´n–матрица A ранга k£min(m,n) и действительный m–вектор b. Найти действительный n–вектор x₀, минимизирующий евклидову длину вектора невязки Ax–b.

Пусть y – n–мерный вектор фактических значений, x – n–мерный вектор значений независимой переменной, b – коэффициенты в аппроксимации y линейной комбинацией n заданных базисных функций j:

Задача состоит в том, чтобы в уравнении подобрать такие b, чтобы минимизировать суммы квадратов отклонений e=y–Xb, где X – есть так называемая матрица плана, в которой строками являются n–мерный вектора с компонентами, зависящими от x_j:  каждая строка соответствует определенному значению x_j. Коэффициенты можно найти решая нормальные уравнения е:

т. к.  

Это выражение имеет экстремум в точке, где

Откуда и получаем

Следует отметить, что последнее выражение имеет в определенной степени формальный характер, т. к. решение нормальных уравнений, как правило, проводится без вычисления обратной матрицы (метод Крамера) такими методами как метод Гаусса, Холесского и т. д.

Пример. Пусть заданы результаты четырех измерений (рис. 1): y=0 при x=0; y=1 при x=1; y=2 при x=3; y=5 при x=4. Задача заключается в том, чтобы провести через эти точки прямую  таким образом, чтобы сумма квадратов отклонений была минимальна. Запишем уравнение, описывающее проведение прямой  по результатам измерений. Мы получаем переопределенную систему:

или Xb=y. Нам понадобится матрица X^TX и обратная к ней:

Тогда решение b=(X^TX)^-1X^Ty по методу наименьших квадратов будет иметь вид

Таким образом, оптимальная прямая задается уравнением  Метод точечной квадратичной аппроксимации (метод наименьших квадратов) не предполагает, что мы должны приближать экспериментальные данные лишь с помощью прямых линий. Во многих экспериментах связи могут быть нелинейными, и было бы глупо искать для этих задач линейные соотношения. Пусть, например, мы работаем с радиоактивным материалом. Тогда выходными данными у являются показания счетчика Гейгера в различные моменты времени t. Пусть наш материал представляет собой смесь двух радиоактивных веществ, и мы знаем период полураспада каждого из них, но не знаем, в каких пропорциях эти вещества смешаны. Если обозначить их количества через С и D, то показания счетчика будут вести себя подобно сумме двух экспонент, а не как прямая:

                                                   (1)

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

Рис. 1. Аппроксимация прямой линией.

соответствовать (1). Вместо этого мы имеем серию показаний счетчика  в различные моменты времени

Если мы имеем более двух показаний, m>2, то точно разрешить эту систему относительно C и D практически невозможно. Но мы в состоянии получить приближенное решение в смысле минимальных квадратов.

Ситуация будет совершенно иной, если нам известны количества веществ C и D и нужно отыскать коэффициенты l и m. Это нелинейная задача наименьших квадратов, и решить ее существенно труднее. Мы по–прежнему будем минимизировать сумму квадратов ошибок, но сейчас она уже не будет многочленом второй степени относительно l и m, так что приравнивание нулю производной не будет давать линейных уравнений для отыскания оптимальных решений.

Глава 1. Метод наименьших квадратов

1.1. Задача наименьших квадратов

Задача наименьших квадратов заключается в минимизация евклидовой длины вектора невязок || Ax-b ||.

Теорема 1. Пусть А – m´n–матрица ранга k, представленная в виде 

A=HRK^T                                                     (2)

где H ортогональная m´m матрица; R – m´n–матрица вида

                                              (3)

где: R₁₁  –  kxk–матрица ранга k; K –  ортогональная  kxk–матрица. Определим вектор

                                          (4)

и введем новую переменную

                                      (5)

Определим  как единственное решение системы  R₁₁y₁=g₁. Тогда:

1.            ||Ax-b|| имеют вид y₂ произвольно.

2.              приводит к одному и тому же вектору невязки                                              (6)

3.            r справедливо

4.           

Доказательство. В выражении для квадрата нормы невязки заменим A на HRK^T в соответствии с (2) и умножая на ортогональную матрицу H^T  (умножение на ортогональную матрицу не меняет евклидову норму вектора) получим

                    (7)

Далее из (3) и (5) следует, что

Из (4) следует

Подставляя оба последних выражения в (7) получим

Последнее выражение имеет минимальное значение  при R₁₁y₁=g₁, а в этом уравнении единственным решением является R₁₁ равен к. Общее решение y выражается формулой y₂ произвольно. Для вектора  имеем

что устанавливает равенство (3). Среди векторов  наименьшую длину имеет тот, для которого y₂=0. Отсюда следует, что решением наименьшей длины будет вектор

Всякое разложение матрицы А типа (2) мы будем называть ортогональным разложением А. Заметим, что решение минимальной длины, множество всех решений и минимальное значение для задачи минимизации ||Ax-b|| определяются единственным образом. Они не зависят от конкретного ортогонального разложения.

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

1.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' совпадают и равны u_m – собственный вектор матрицы X'X, а v_m – собственный вектор матрицы XX', соответствующие одному и тому же собственному числу u_m и v_m связаны следующим соотношением

                        (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–мерном евклидовом пространстве.

1.5. QR–разложение

Теорема 2. Пусть А – m´n–матрица. Существует ортогональная m´m–матрица Q такая, что в матрице QA=R под главной диагональю стоят только нулевые элементы.

Доказательство. Выберем ортогональную m´m–матрицу Q в соответствии с преобразованием Хаусхолдера (9), так, чтобы первый столбец Q₁A имел нулевые компоненты со 2–ой по m–ю. Далее выбираем ортогональную (m-1)´(m–1)–матрицу P₂ следующим образом. Будучи применена к  m–1 вектору, составленному из компонент со 2–ой по m–ю второго столбца матрицы Q₁A, она аннулирует компоненты с 3–ей по m–ю этого вектора. Матрица преобразования

ортогональна, и Q₂Q₁A имеет в первых двух столбцах нули под главной диагональю. Продолжая таким образом, можно построить произведение, состоящее максимум из 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, A^T, Q^T теоремы 3.

Используя теорему 3 и лемму 1 можно доказать следующую теорему.

Теорема 4. Пусть А – m´n–матрица ранга к . Найдутся ортогональная m´m–матрица Н и ортогональная n´n–матрица К такие, что

                                        (19)

причем R₁₁ – невырожденная треугольная к´к–матрица.

Заметим, что выбором Н и К в уравнении (19) можно добиться, чтобы R₁₁ была верхней или нижней треугольной.

В (19) матрица А представлена произведением A=HRK^T, где R – некоторая прямоугольная матрица, ненулевые компоненты которой сосредоточены в невырожденной треугольной подматрице. Далее мы покажем, что эту невырожденную подматрицу R можно упростить далее до невырожденной диагональной матрицы. Это разложение тесно связано со спектральным разложением симметричных неотрицательно определенных матриц A^TA и AA^T (см. 11).

Теорема 5. Пусть А – m´n–матрица ранга k. Тогда существуют ортогональная m´m–матрица U, ортогональная n´n–матрица V и диагональная m´n–матрица S такие, что

U^TAV=S, A=USV^T                                      (20)

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

Диагональные элементы S называются сингулярными числами А. Докажем сперва лемму для специального случая m=n=rankA.

Лемма 2. Пусть А – n´n–матрица ранга n. Тогда существует ортогональная n´n–матрица U, ортогональная n´n–матрица V и диагональная n´n–матрица S такие, что U^TAV=S, A=USV^T и последовательные диагональные элементы S положительны и не возрастают.

Доказательство леммы. Положительно определенная симметричная матрица A^TA допускает спектральное разложение

A^TA=VDV^T,                                        (21)

где V – ортогональная n´n–матрица, а D – диагональная матрица, причем диагональные элементы D положительны и не возрастают. Определим S как диагональную n´n–матрицу, диагональные элементы которой суть положительные квадратные корни из соответствующих диагональных элементов D. Таким образом

D=S^TS=S², S^-1D^S-1=I.                                  (22)

Определим матрицу

U=AVS^-1                                                     (23)

Из (21), (22), (23) и ортогональности V следует, что

U^TU=S^-1V^TA^TAVS^-1=S^-1DS^-1=I т.е. U ортогональна. Из (23) и ортогональности V выводим USV^T=AVS^-1SV^T=AVV^T=A Лемма доказана.

Доказательство теоремы 5. Пусть A=HRK^T, где H, R, K^T имеют свойства, указанные в теореме 4. Так как R₁₁ из (19) – невырожденная треугольная к´к–матрица, то согласно лемме 2 , можно написать

                                                  (24)

Здесь  и  – ортогональные к´к–матрицы, а  – невырожденная диагональная матрица, диагональные элементы которой положительны и не возрастают. Из (24) следует, что матрицу R в уравнении (19) можно записать в виде

                                                  (25)

где:

       – ортогональная m´m–матрица;

       – ортогональная n´n–матрица;

           – ортогональная m´n–матрица;

Теперь, определяя U и V формулами

                                       (26)

заключаем из (24) – (26), что A=USV^T, где U, S, V имеют свойства, указанные в формулировке теоремы 5. Это завершает доказательство.

  Заметим, что сингулярные числа матрицы А определены однозначно, в то время, как в выборе ортогональных матриц U, V есть произвол. Пусть s – сингулярное число А, имеющее кратность l. Это значит, что для упорядоченных сингулярных чисел найдется индекс I такой, что

  Положим k=min(m,n), и пусть Q – ортогональная к´к–матрица вида

Здесь Р – ортогональная l´l–матрица Если A=USV^T – сингулярное разложение А и s_i=…=s_i+l-1, то сингулярным разложением А будет также и         

1.6. Число обусловленности

Некоторые вычислительные задачи поразительно чувствительны к изменению данных. Этот аспект численного анализа не зависит от плавающей арифметики или выбранного алгоритма.

Например:

Найти корни полинома: (x-2)²=10^-6

Корни этого уравнения есть 2+10^-3 и 2-10^-3. Однако изменение свободного члена на 10^-6 может вызвать изменение в корнях, равное 10^-3.

Операции с матрицами, как правило, приводят к решению систем линейных уравнений. Коэффициенты матрицы в правой части системы линейных уравнений редко известны точно. Некоторые системы возникают из эксперимента, и тогда коэффициенты подвержены ошибкам наблюдения. Коэффициенты других систем записываются формулами, что влечет за собой ошибки округлений. В связи с этим необходимо знать, как влияют ошибки в коэффициентах матрицы на решение. Именно для этого вводится понятие обусловленности матрицы.

По определению число обусловленности есть величина

Нормой вектора x в пространстве векторов  называется функционал, обозначаемый

1)          

2)          

3)          

Нормой квадратной матрицы А в пространстве матриц, согласованной с нормой вектора  называется функционал 

1)          

2)          

3)          

4)          

Наиболее употребимы следующие нормы для векторов:

·                         

·                                  

·                    

Нормы матриц:

·              

·              

·              

Здесь  являются сингулярными числами[3] матрицы А; это положительные значения квадратных корней  из собственных значений  матрицы А^ТА (которая при невырожденной матрице А положительно определена[4], в противном случае положительно полуопределена (неотрицательно определена[5]) и поэтому имеет только вещественные собственные значения ³ 0). Для вещественных симметричных матриц сингулярные числа равны абсолютным величинам собственных значений:

Умножение вектора х на матрицу А приводит к новому вектору Ах, норма которого может очень сильно отличаться от нормы вектора х.

Область изменений может быть задана двумя числами

Максимум и минимум берутся по всем ненулевым векторам. Заметим, что если А вырождена, то m=0. Отношение M/m называется числом обусловленности матрицы А,

                                (7)

Рассмотрим норму обратной[6] матрицы

Для матрицы А существует сингулярное разложение  и  – 1/ есть величины, обратные собственным числам матрицы  –

Рассмотрим систему уравнений 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–алгоритм без сдвигов. Обычно время которое тратится на такой процесс пропорционально кубу размерности матрицы – n³. Необходимо процесс ускорить, для чего используется предварительное приведение матрицы А к форме Хессенберга[8] а также используется алгоритм со сдвигом. Форма Хессенберга представляет из себя верхнюю треугольную матрицу (верхняя форма Хессенберга) у которой сохранена одна диагональ ниже главной, а элементы ниже этой диагонали равны нулю. Если матрица симметрична, то легко видеть, что матрица Хессенберга превращается в трехдиагональную матрицу[9]. При использовании матрицы Хессенберга время процесса пропорционально n², а при использовании трехдиагональной матрицы – n.

Можно использовать другие соотношения

где Q_s – унитарная, а L_s – нижняя треугольная матрица. Такой алгоритм носит название QL–алгоритма.

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

В общем случае, наддиагональный  элемент  матрицы A_s на s-ом шаге асимптотически равен k_ij – постоянная величина. Сходимость QL–алгоритма вообще говоря недостаточна. Сходимость можно улучшить, если на каждом шаге вместо матрицы A_s использовать матрицу A_s-k_sI (QL–алгоритм со сдвигом). Последовательность вычислений в этом случае описывается следующими соотношениями:

которые определяют матрицу  определено соотношением k_s выбрать близко к величине  (наименьшее собственное значение), то в пределе внедиагональные элементы первой строки будут очень быстро стремиться к нулю. Когда ими можно пренебречь, элемент  с рабочей точностью равен n-1-го порядка. Тогда, если QL–алгоритм выполнен без ускорения сходимости, то все равно k_s.

Если матрица А эрмитова, то очевидно, что и все матрицы А_s эрмитовы; если А действительная и симметричная, то все Q_s ортогональны и все А_s действительны и симметричны.

2.2. Реализация разложения

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

Далее реализуется итерационный процесс приведения двухдиагональной матрицы J₀ к диагональной форме, так что имеет место следующая последовательность:  где  а S_i и T_i  – диагональные матрицы.

Матрицы T_i выбираются так, чтобы последовательность матриц  сходилась к двухдиагональной матрице. Матрицы же S_i выбирают так, чтобы все J_i  сохраняли двухдиагональную форму. Переход  осуществляется с помощью плоских вращений (10) – преобразований Гивенса. Отсюда,  где

а матрица  вычисляется аналогично с заменой  на

Пусть начальный угол  произволен, однако следующие значения угла необходимо выбирать так, чтобы матрица J_i+1 имела ту же форму, что и J_i. Таким образом  не аннулирует  ни одного элемента матрицы, но добавляет элемент  аннулирует  но добавляет  аннулирует  но добавляет  и т.д., наконец,  аннулирует  и ничего не добавляет.

Этот процесс часто называют процессом преследования. Так как M_i+1 – трехдиагональная матрица, точно так же, как и M_i. Начальный угол  можно выбрать так, чтобы преобразование  было QR–преобразованием со сдвигом, равным s.

Обычный QR–алгоритм со сдвигом можно записать в следующем виде:

где  – верхняя треугольная матрица. Следовательно, s  определяется собственным значением нижнего минора (размерности 2´2) матрицы M_i. При таком выборе параметра s метод обладает глобальной и почти всегда кубичной сходимостью.

2.3. Пример сингулярного разложения

Проведем преобразование Хаусхолдера на матрице

К первой компоненте первого столбца прибавляем норму первого столбца, получим

Преобразованная матрица A₂ вычисляется следующим образом. Для первого столбца имеем:

так как

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

для третьего столбца:

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

Столбцы матрицы A₂ получаются вычитанием кратных вектора v₁ из столбцов A₁. Эти кратные порождаются скалярными произведениями, а не отдельными элементами, как в гауссовом исключении.

Прежде чем вводить дальнейшие нули под диагональю, преобразованием вида A₃=A₂Q₁, получают нули в первой строке. Нули уже стоящие в первом столбце, не должны быть испорчены, длина первого столбца должна быть сохранена; поэтому при внесении нулей в первую строку нельзя менять первый элемент строки, изменяем второй элемент и зануляем третий. Для матрицы большего размера на этом шаге было бы получено n–2 нуля.

Преобразование порождается первой строкой A₂:

Строка матрицы A₃ с номером i получается по формуле

Таким образом, из каждой строки A₂ вычитается надлежащее кратное

Поскольку первая компонента A₂ сохраняются в A₃, Так как Q₁ ортогональная, то длина каждой строки в A₃ равна длине соответствующей строки в A₂.

Теперь можно добиться новых нулей под диагональю, не испортив полученных ранее:

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

Если бы не ошибки округлений, то в данном примере третий диагональный элемент был бы точным нулем. Элементы под диагональю во всех столбцах указаны как точные нули, потому что преобразования так и строились, чтобы получить там нули. Последнее преобразование H₃ в этом примере могло бы быть тождественным, однако тогда оно не было бы хаусхолдеровым отражением. Фактически использование H₃ попутно изменяет знак  элемента – 1.080 в матрице A₄.

Получилась искомая двухдиагональная матрица, и первый этап закончен. Прямое использование ортогональных преобразований не позволяет получить какие–либо новые нули. Для общего порядка n нужно n преобразований H и n–2  преобразований Q, чтобы достигнуть данного места. Число преобразований не зависит от строчной размерности m, но от m зависит работа, затрачиваемая на выполнение каждого преобразования.

глава 3. Использование сингулярного разложения  в методе наименьших квадратов

При использовании метода сингулярного разложения (SVD – Singular Value Decomposition) мы проводим разложение для матрицы плана. При этом основное уравнение y=Xb приобретает вид y=USV^Tb. Отсюда следует, что коэффициенты b можно получить решая уравнение U^Ty=SV^Tb. Т.е. если все s_j, j=1,…,n, являющиеся диагональными элементами S отличны от нуля, то последнее уравнение разрешимо и

Однако такое решение не всегда желательно, если некоторые s_j малы. Для правильного использования метода SVD мы должны ввести границу t отражающую точность входных данных и точность использованной плавающей арифметики. Всякое s_j, большее, чем t, приемлемо, и соответствующее  вычисляется по (1.20). Любое s_j, меньшее, чем t, рассматривается как пренебрежимо малое, и соответствующему  может быть придано произвольное значение. С этой произвольностью значений связана не единственность набора коэффициентов, получаемых методом наименьших квадратов. Изменения входных данных и ошибки округлений, меньшие, чем t, могут привести к совершенно другому набору коэффициентов, определяемых этим методом. Поскольку обычно желательно, чтобы эти коэффициенты были по возможности малы, то полагаем s_j £t.

Отбрасывание чисел s_j, меньших, чем t, приводит к уменьшению числа обусловленности с  до

Продемонстрируем использование метода на следующем примере:

Следует определить значение Y при X =1980.

Если аппроксимировать эти данные квадратичным многочленом  и использовать двойную точность, то в результате получим следующие коэффициенты Y=227780000 , а для обычной точности Y=145210000. Ясно, что второй набор коэффициентов бесполезен. Исследуем достоверность результатов. Матрица плана для данного примера имеет размеры 8 ´ 3

Рис. 2. Численный пример

Ее сингулярные числа:

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

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

Теперь коэффициенты находятся в гораздо лучшем согласии друг с другом. Кроме того, коэффициенты стали существенно меньше, а это значит, что не будет столь большого, как прежде, взаимного уничтожения слагаемых при вычислении квадратичного многочлена. Прогнозное значение Y(1980) будет соответственно 212910000 и 214960000. Эффект обычной точности еще заметен, однако результаты уже не являются катастрофическими.

Можно также определить набор нулевых коэффициентов, соответствующих пренебрежимо малому сингулярному числу. Вот эти коэффициенты: t от 1900 до 1970 величина функции  не превосходит 0.0017, поэтому при любом a коэффициенты можно изменить a. Любой из четырех перечисленных нами наборов коэффициентов можно получить из другого подобным изменением.

Во–вторых, можно улучшить ситуацию заменой базиса. Модели

гораздо более удовлетворительны. Важно при этом то, что независимая переменная преобразуется из интервала [1900, 1970] в какой–нибудь более приемлемый интервал вроде [0, 70] или, еще лучше, [–3.5, 3.5]. Числа обусловленности при этом равны 5750 и 10.7 соответственно. последнее значение более чем приемлемо даже при счете с обычной точностью.

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

Во втором случае, после преобразования матрицы плана ее обусловленность сильно уменьшается, и, соответственно, повышается точность расчетов.

Данную программу можно использовать и при решении системы линейных уравнений вместо методов Гаусса, Жордана, Холесского и пр. В приложении 2 приведен пример расчета линейной системы, которая изначально не может быть решена этими методами вследствие вырожденности матрицы коэффициентов. Тем не менее, исследуемый метод дает нам правильное решение.

ЗАКЛЮЧЕНИЕ

В работе описаны компьютерные методы решения задачи наименьших квадратов. Для использования данных методов составлена соответствующая программа на алгоритмическом языке FORTRAN. Программа апробирована, результаты тестирования показывают работоспособность программы.

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

ЛИТЕРАТУРА

1.           

2.           

3.           

4.           

5.           

6.           

7.           

8.           

9.           

10.      

11.      

ПРИЛОЖЕНИЕ 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)

                CALL SVD(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) GO TO 565

                IF(ABS(W(L1))+ANORM.EQ.ANORM) GO TO 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 TO 565

                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   .7400000E 02

 -0.2400000E 02 -0.1000000E 02 -0.5700000E 02

 -0.8000000E 01 -0.4000000E 01 -0.1700000E 02

 -0.1400000E 02  0.1300000E 02  0.1000000E 01

Полученный результат

МАТРИЦА А

   .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

Видно, что правые части соответствуют начальным данным. Решение верно.