СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ В ЛИНЕЙНОЙ ЗАДАЧЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ

СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ В ЛИНЕЙНОЙ ЗАДАЧЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ

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

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

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

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

СИНГУЛяРНОЕ РАЗЛОЖЕНИЕ В ЛИНЕЙНОЙ ЗАДАчЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ

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

математики

Исполнил:

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

Владикавказ 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 заданных базисных функций (:

[pic].

Задача состоит в том, чтобы в уравнении подобрать такие b, чтобы

минимизировать суммы квадратов отклонений e=y–Xb, где X – есть так

называемая матрица плана, в которой строками являются n–мерный вектора с

компонентами, зависящими от xj: [pic] каждая строка соответствует

определенному значению xj. Коэффициенты можно найти решая нормальные

уравнения [pic], откуда [pic]. Покажем это. Возведем в квадрат выражение

для е:

[pic]

т. к. [pic].

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

Откуда и получаем [pic].

Следует отметить, что последнее выражение имеет в определенной степени

формальный характер, т. к. решение нормальных уравнений, как правило,

проводится без вычисления обратной матрицы (метод Крамера) такими методами

как метод Гаусса, Холесского и т. д.

Пример. Пусть заданы результаты четырех измерений (рис. 1): y=0 при

x=0; y=1 при x=1; y=2 при x=3; y=5 при x=4. Задача заключается в том, чтобы

провести через эти точки прямую [pic] таким образом, чтобы сумма квадратов

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

прямой [pic] по результатам измерений. Мы получаем переопределенную

систему:

[pic]

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

[pic]

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

вид [pic]

Таким образом, оптимальная прямая задается уравнением [pic] Метод

точечной квадратичной аппроксимации (метод наименьших квадратов) не

предполагает, что мы должны приближать экспериментальные данные лишь с

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

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

например, мы работаем с радиоактивным материалом. Тогда выходными данными у

являются показания счетчика Гейгера в различные моменты времени t. Пусть

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

знаем период полураспада каждого из них, но не знаем, в каких пропорциях

эти вещества смешаны. Если обозначить их количества через С и D, то

показания счетчика будут вести себя подобно сумме двух экспонент, а не как

прямая:

[pic]. (1)

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

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

[pic]

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

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

различные моменты времени [pic], и (1) выполняется лишь приближенно:

[pic]

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

относительно C и D практически невозможно. Но мы в состоянии получить

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

Ситуация будет совершенно иной, если нам известны количества веществ C

и D и нужно отыскать коэффициенты ( и (. Это нелинейная задача наименьших

квадратов, и решить ее существенно труднее. Мы по–прежнему будем

минимизировать сумму квадратов ошибок, но сейчас она уже не будет

многочленом второй степени относительно ( и (, так что приравнивание нулю

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

решений.

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

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

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

вектора невязок (( Ax-b ((.

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

A=HRKT (2)

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

[pic], (3)

где: R11 – kxk–матрица ранга k; K – ортогональная kxk–матрица.

Определим вектор

[pic] (4)

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

[pic]. (5)

Определим [pic] как единственное решение системы R11y1=g1. Тогда:

1. Все решения задачи о минимизации ((Ax-b(( имеют вид [pic], где y2

произвольно.

2. Любой такой вектор [pic] приводит к одному и тому же вектору невязки

[pic]. (6)

3. Для нормы r справедливо [pic]

4. Единственным решением минимальной длины является вектор [pic]

Доказательство. В выражении для квадрата нормы невязки заменим A на

HRKT в соответствии с (2) и умножая на ортогональную матрицу HT (умножение

на ортогональную матрицу не меняет евклидову норму вектора) получим

[pic] (7)

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

[pic].

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

[pic]

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

[pic]

Последнее выражение имеет минимальное значение [pic] при R11y1=g1, а в

этом уравнении единственным решением является [pic], так как ранг матрицы

R11 равен к. Общее решение y выражается формулой [pic], где y2 произвольно.

Для вектора [pic] имеем

[pic],

что устанавливает равенство (3). Среди векторов [pic] наименьшую длину

имеет тот, для которого y2=0. Отсюда следует, что решением наименьшей длины

будет вектор [pic]. Теорема доказана.

Всякое разложение матрицы А типа (2) мы будем называть ортогональным

разложением А. Заметим, что решение минимальной длины, множество всех

решений и минимальное значение для задачи минимизации ((Ax-b(( определяются

единственным образом. Они не зависят от конкретного ортогонального

разложения.

При проведении разложения необходимо приводить матрицы к диагональному

виду. Для этого обычно используются два преобразования: Гивенса и

Хаусхолдера, оставляющие нормы столбцов и строк матриц неизменными.

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

Лемма. Пусть дан 2–вектор [pic], причем [pic] либо [pic].Существует

ортогональная 2(2 матрица такая, что:

[pic] (8)

Доказательство. Положим:

[pic].

Далее прямая проверка.

Матрица преобразования представляет собой матрицу вращений

[pic]

или отражений

[pic]

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

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

преобразования представляет из себя следующее выражение: [pic],

(9)

или, если вектор v нормирован, т.е. используется вектор единичной длины

[pic], то [pic]. В обоих случаях H – симметричная и ортогональная матрица.

Покажем это:

[pic].

Отсюда следует: что [pic], т.е. симметричность и ортогональность. В

комплексном случае матрица [pic] эрмитова[1] и унитарна[2]. Предположим,

что дан вектор х размерности m, тогда существует матрица H такая, что

[pic], где

[pic]

а ( = +1, при положительной первой компоненте вектора х и = –1, при

отрицательной.

Доказательство. Положим [pic] действительная матрица. Любую

действительную матрицу можно привести в треугольному виду

[pic]

Далее принимаем во внимание то, что [pic] и получаем следующее:

[pic]

1.4. Сингулярное разложение матриц

Пусть X – матрица данных порядка Nxp, где N>p, и пусть r – ранг матрицы

X. Чаще всего r=p, но приводимый ниже результат охватывает общий случай, он

справедлив и при условии rk, что противоречит предположению rankA=k. Итак, QAP имеет форму,

указанную правой частью (4). Теорема доказана.

Подматрицу [R:T] из правой части (18) можно теперь преобразовать к

компактной форме, требуемой от матрицы R из теоремы 2. Это преобразование

описывает следующая лемма.

Лемма 1. Пусть [R:T] – к(к–матрица, причем R имеет ранг к. Существует

ортогональная n(n–матрица W такая, что

[pic]

где [pic] – нижняя треугольная матрица ранга к.

Доказательство леммы вытекает из теоремы 3, если отождествить величины

n, k, [R:T], W из формулировки леммы с соответствующими величинами m, n,

AT, QT теоремы 3.

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

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

m(m–матрица Н и ортогональная n(n–матрица К такие, что

[pic] (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 , можно написать

[pic] (24)

Здесь [pic] и [pic] – ортогональные к(к–матрицы, а [pic] –

невырожденная диагональная матрица, диагональные элементы которой

положительны и не возрастают. Из (24) следует, что матрицу R в уравнении

(19) можно записать в виде

[pic] (25)

где:

[pic] – ортогональная m(m–матрица;

[pic] – ортогональная n(n–матрица;

[pic] – ортогональная m(n–матрица;

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

[pic] (26)

заключаем из (24) – (26), что A=USVT, где U, S, V имеют свойства,

указанные в формулировке теоремы 5. Это завершает доказательство.

Заметим, что сингулярные числа матрицы А определены однозначно, в то

время, как в выборе ортогональных матриц U, V есть произвол. Пусть ( –

сингулярное число А, имеющее кратность l. Это значит, что для упорядоченных

сингулярных чисел найдется индекс I такой, что

[pic]

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

[pic]

Здесь Р – ортогональная l(l–матрица Если A=USVT – сингулярное

разложение А и si=…=si+l-1, то сингулярным разложением А будет также и

[pic], где [pic].

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

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

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

или выбранного алгоритма.

Например:

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

Корни этого уравнения есть 2+10-3 и 2-10-3. Однако изменение свободного

члена на 10-6 может вызвать изменение в корнях, равное 10-3.

Операции с матрицами, как правило, приводят к решению систем линейных

уравнений. Коэффициенты матрицы в правой части системы линейных уравнений

редко известны точно. Некоторые системы возникают из эксперимента, и тогда

коэффициенты подвержены ошибкам наблюдения. Коэффициенты других систем

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

этим необходимо знать, как влияют ошибки в коэффициентах матрицы на

решение. Именно для этого вводится понятие обусловленности матрицы.

По определению число обусловленности есть величина [pic]. Для более

подробного описания числа обусловленности нам понадобится понятие нормы в

пространстве векторов и матриц.

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

обозначаемый [pic], удовлетворяющий следующим условиям:

1) положительной определенности – [pic]

2) положительной однородности – [pic];

3) неравенству треугольника – [pic].

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

нормой вектора [pic] называется функционал [pic], удовлетворяющий условиям

1 – 3 для нормы вектора:

1) [pic];

2) [pic]

3) [pic]

4) мультипликативное неравенство – [pic]

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

. норма суммы модулей [pic]

. евклидова норма [pic]

. норма максимума модуля [pic]

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

. [pic]

. [pic]

. [pic]

Здесь [pic] являются сингулярными числами[3] матрицы А; это

положительные значения квадратных корней [pic] из собственных значений

[pic] матрицы АТА (которая при невырожденной матрице А положительно

определена[4], в противном случае положительно полуопределена

(неотрицательно определена[5]) и поэтому имеет только вещественные

собственные значения ( 0). Для вещественных симметричных матриц сингулярные

числа равны абсолютным величинам собственных значений: [pic].

Умножение вектора х на матрицу А приводит к новому вектору Ах, норма

которого может очень сильно отличаться от нормы вектора х.

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

[pic]

Максимум и минимум берутся по всем ненулевым векторам. Заметим, что

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

матрицы А,

[pic] (7)

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

Для матрицы А существует сингулярное разложение [pic], тогда [pic],

отсюда [pic]. Аналогично для обратной матрицы [pic] и [pic]. Отсюда

следует, что собственные числа матрицы [pic] – 1/[pic] есть величины,

обратные собственным числам матрицы [pic] – [pic]. При этом очевидно, что

[pic]. Из последнего выражения вместе с (7) следует [pic]. Таким образом

обусловленность матрицы равна произведению нормы матрицы на норму обратной

матрицы.

Рассмотрим систему уравнений Ax=b, и другую систему, полученную

изменением правой части: A(x+(x)=b+(b . Будем считать (b ошибкой в b, а (x

соответствующей ошибкой в x, хотя нам нет необходимости считать ошибки

малыми. Поскольку A((x)=(b, то определения M и m немедленно приводят к

неравенствам [pic] Следовательно , при m(0,

[pic]

Величина [pic] есть относительное изменение правой части, а величина

[pic] – относительная ошибка, вызванная этим изменением. Аналогичные

выкладки можно провести не только с элементами вектора правой части но и с

элементами самой матрицы А и найти зависимость между относительным

изменением элементов матрицы и относительной ошибкой вызванной этим

изменением. Отсюда следует, что число обусловленности выполняет роль

множителя в увеличении относительной ошибки.

Приведем некоторые свойства числа обусловленности. Ясно, что M(m и

поэтому cond(А)(1. Если Р – матрица перестановок[7], то компоненты вектора

Px лишь порядком отличаются от компонент вектора х. Отсюда следует, что

[pic] и cond(P)=1 . В частности cond(I)=1. Если А умножается на скаляр с,

то cond(cА)= cond(А). Если D – диагональная матрица, то [pic]

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

2.1. Алгоритмы

QR–алгоритм начинается с разложения матрицы по Грамму-Шмидту [pic],

затем меняются местами сомножители: [pic] Эта матрица подобна

первоначальной, [pic] Этот процесс продолжается, причем собственные

значения не изменяются:

[pic]

Эта формула описывает QR–алгоритм без сдвигов. Обычно время которое

тратится на такой процесс пропорционально кубу размерности матрицы – n3.

Страницы: 1, 2