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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ Математический факультет Кафедра прикладной математики ДИПЛОМНЫЙ ПРОЕКТ сингулярное разложение в линейной задаче метода наименьших квадратов Заведующий кафедрой прикладной математики Исполнил Научный руководитель Владикавказ 2002 СОДЕРЖАНИЕ ВВЕДЕНИЕ 3 Глава 1. Метод наименьших квадратов 7 1.1. Задача наименьших квадратов 2. Ортогональное вращение Гивенса 3. Ортогональное преобразование Хаусхолдера 4. Сингулярное разложение матриц 11 1.5. QR разложение 6. Число обусловленности 20 глава 2. Реализация сингулярного разложения 1. Алгоритмы 2. Реализация разложения 3. Пример сингулярного разложения 29 глава 3. Использование сингулярного разложения в методе наименьших квадратов 33 ЗАКЛЮЧЕНИЕ 38 ЛИТЕРАТУРА 39 ПРИЛОЖЕНИЕ 1. Исходные тексты программы 40 ПРИЛОЖЕНИЕ 2. контрольный пример 45 ВВЕДЕНИЕ Метод наименьших квадратов обычно используется как составная часть некоторой более общей проблемы.

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

Пусть даны действительная mn матрица A ранга kminm, n и действительный m вектор b. Найти действительный n вектор x0, минимизирующий евклидову длину вектора невязки Ax b. Пусть y n мерный вектор фактических значений, x n мерный вектор значений независимой переменной, b коэффициенты в аппроксимации y линейной комбинацией n заданных базисных функций. Задача состоит в том, чтобы в уравнении подобрать такие b, чтобы минимизировать суммы квадратов отклонений ey Xb, где X есть так называемая матрица плана, в которой строками являются n мерный вектора с компонентами, зависящими от xj каждая строка соответствует определенному значению xj. Коэффициенты можно найти решая нормальные уравнения, откуда. Покажем это. Возведем в квадрат выражение для е т. к Это выражение имеет экстремум в точке, где 0 Откуда и получаем. Следует отметить, что последнее выражение имеет в определенной степени формальный характер, т. к. решение нормальных уравнений, как правило, проводится без вычисления обратной матрицы метод Крамера такими методами как метод Гаусса, Холесского и т. д. Пример.

Пусть заданы результаты четырех измерений рис. 1 y0 при x0 y1 при x1 y2 при x3 y5 при x4. Задача заключается в том, чтобы провести через эти точки прямую таким образом, чтобы сумма квадратов отклонений была минимальна.

Запишем уравнение, описывающее проведение прямой по результатам измерений.

Мы получаем переопределенную систему или Xby. Нам понадобится матрица XTX и обратная к ней Тогда решение bXTX-1XTy по методу наименьших квадратов будет иметь вид Таким образом, оптимальная прямая задается уравнением Метод точечной квадратичной аппроксимации метод наименьших квадратов не предполагает, что мы должны приближать экспериментальные данные лишь с помощью прямых линий.

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

Пусть, например, мы работаем с радиоактивным материалом. Тогда выходными данными у являются показания счетчика Гейгера в различные моменты времени t. Пусть наш материал представляет собой смесь двух радиоактивных веществ, и мы знаем период полураспада каждого из них, но не знаем, в каких пропорциях эти вещества смешаны. Если обозначить их количества через С и D, то показания счетчика будут вести себя подобно сумме двух экспонент, а не как прямая . 1 На практике, поскольку радиоактивность измеряется дискретно и через различные промежутки времени, показания счетчика не будут точно Рис. 1. Аппроксимация прямой линией. соответствовать 1. Вместо этого мы имеем серию показаний счетчика в различные моменты времени, и 1 выполняется лишь приближенно Если мы имеем более двух показаний, m 2, то точно разрешить эту систему относительно C и D практически невозможно.

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

Ситуация будет совершенно иной, если нам известны количества веществ C и D и нужно отыскать коэффициенты и. Это нелинейная задача наименьших квадратов, и решить ее существенно труднее. Мы по прежнему будем минимизировать сумму квадратов ошибок, но сейчас она уже не будет многочленом второй степени относительно и, так что приравнивание нулю производной не будет давать линейных уравнений для отыскания оптимальных решений. Глава 1. Метод наименьших квадратов 1. Задача наименьших квадратов Задача наименьших квадратов заключается в минимизация евклидовой длины вектора невязок Ax-b. Теорема 1. Пусть А mn матрица ранга k, представленная в виде AHRKT 2 где H ортогональная mm матрица R mn матрица вида , 3 где R11 kxk матрица ранга k K ортогональная kxk матрица. Определим вектор 4 и введем новую переменную . 5 Определим как единственное решение системы R11y1g1. Тогда 1. Все решения задачи о минимизации Ax-b имеют вид, где y2 произвольно. 2. Любой такой вектор приводит к одному и тому же вектору невязки . 6 3. Для нормы r справедливо 4. Единственным решением минимальной длины является вектор Доказательство. В выражении для квадрата нормы невязки заменим A на HRKT в соответствии с 2 и умножая на ортогональную матрицу HT умножение на ортогональную матрицу не меняет евклидову норму вектора получим 7 Далее из 3 и 5 следует, что. Из 4 следует Подставляя оба последних выражения в 7 получим Последнее выражение имеет минимальное значение при R11y1g1, а в этом уравнении единственным решением является, так как ранг матрицы R11 равен к. Общее решение y выражается формулой, где y2 произвольно.

Для вектора имеем, что устанавливает равенство 3. Среди векторов наименьшую длину имеет тот, для которого y20. Отсюда следует, что решением наименьшей длины будет вектор. Теорема доказана.

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

Они не зависят от конкретного ортогонального разложения.

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

Ортогональное вращение Гивенса

Ортогональное вращение Гивенса Лемма. Пусть дан 2 вектор, причем либо. Существует ортогональная 22 матрица такая, что 8 Доказательство.

Положим. Далее прямая проверка. Матрица преобразования представляет собой матрицу вращений или отражений 1.3.

Ортогональное преобразование Хаусхолдера

Теорема о сингулярном разложении утверждает, что 10 где V матрица поря... Далее выбираем ортогональную m-1m 1 матрицу P2 следующим образом. Пусть RT кк матрица, причем R имеет ранг к. Теорема 5. матрицы.

Реализация сингулярного разложения

Алгоритмы QR алгоритм начинается с разложения матрицы по Грамму-Шмидту... а также используется алгоритм со сдвигом. В общем случае, наддиагональный элемент матрицы As на s-ом шаге асимпт... Последовательность вычислений в этом случае описывается следующими соо... Если сдвиг ks выбрать близко к величине наименьшее собственное значени...

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

Реализация разложения. Таким образом, разложение производится в два этапа. Сначала матрица А ... Начальный угол можно выбрать так, чтобы преобразование было QR преобра... При таком выборе параметра s метод обладает глобальной и почти всегда ... 2.3.

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

Нули уже стоящие в первом столбце, не должны быть испорчены, длина пер... Фактически использование H3 попутно изменяет знак элемента 1.080 в мат... Последнее преобразование H3 в этом примере могло бы быть тождественным... . Пример сингулярного разложения.

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

Отсюда следует, что коэффициенты b можно получить решая уравнение UTyV... если все j, j1 n, являющиеся диагональными элементами отличны от нуля,... При этом получим следующие наборы коэффициентов для двойной и обычной ... Для значений t от 1900 до 1970 величина функции не превосходит 0.0017,... Это означает, что данные берутся в виде отклонений от среднего, которо...

ЗАКЛЮЧЕНИЕ

ЗАКЛЮЧЕНИЕ В работе описаны компьютерные методы решения задачи наименьших квадратов. Для использования данных методов составлена соответствующая программа на алгоритмическом языке 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.

Исходные тексты программы

контрольный пример Входные данные матрица изначально сингулярна первая... DO 40 I1,N SSUI,JYI 40 CONTINUE SSSIGMAJ DO 50 I1,N CICI SVI,J 50 CONT... SIGMA1 SIGMA1SIGMAJ CJ0. DO 30 J1,N IFSIGMAJ .GT. DO 321 I1,M 321 SSAJ,ICISS 322 Y0JSS write 6,570 WRITE 6,6 Y0I,I1,3 C ...