Диплом: Сингулярное разложение в линейной задаче метода наименьших квадратов - текст диплома. Скачать бесплатно.
Банк рефератов, курсовых и дипломных работ. Много и бесплатно. # | Правила оформления работ | Добавить в избранное
 
 
   
Меню Меню Меню Меню Меню
   
Napishem.com Napishem.com Napishem.com

Диплом

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

Банк рефератов / Математика

Рубрики  Рубрики реферат банка

закрыть
Категория: Дипломная работа
Язык диплома: Русский
Дата добавления:   
 
Скачать
Microsoft Word, 4101 kb, скачать бесплатно
Заказать
Узнать стоимость написания уникальной дипломной работы

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

21 МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ Математический факультет Кафедра прикладной математики ДИПЛОМНЫЙ ПРОЕКТ сингулярное разложение в линейной задаче метода наименьших квадратов Заведующий кафедрой прикладной математики Исполнил : Научный руководитель Владикавказ 200 2 СОДЕРЖАНИЕ ВВЕДЕНИЕ 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 – вектор x 0 , минимизирующий евклидову длину вектора невязки Ax – b . Пусть y – n – мерный вектор фактических значений , x – n – мерный вектор значений независимой переменной , b – коэффициенты в аппроксимации y линейной комбина цией n заданных базисных функций : . Задача состоит в том , чтобы в уравнении подобрать такие b , чтобы минимизировать суммы квад ратов отклонений e=y – Xb , где X – есть так называемая матрица плана , в которой строками являются n – мерный вектора с компонентами , зависящими от x j : каждая стро ка соответствует определенному значению x j . Коэффициенты можно найти решая нормальные уравнения , откуда . Покажем это . Возведем в квадрат выражение для е : т . к . . Это выражение имеет экстремум в точке , где =0 Откуда и получаем . Следует отметить , что последнее выражение имеет в определенной степени формальный характер , т . к . решение нормальных уравнений , как правило , проводится без вычисления обратной матрицы (метод Крамера ) таким и методами как метод Гаусса , Холесского и т . д. Пример. Пусть заданы результаты четырех измерений (рис . 1): y= 0 при x= 0 ; y= 1 при x= 1 ; y= 2 при x= 3 ; y= 5 при x= 4 . Задача заключается в том , чтобы провести через эти точки прямую таким образом , чтобы сумма квадратов отклонений была минимальна . Запишем уравнение , описывающее проведение прямой по результатам измерений . Мы получаем переопределенную систему : или Xb=y . Нам понадобится матрица X T X и обратная к ней : Тогда решение b =( X T X ) -1 X T y по методу наименьших квадратов будет иметь вид Таким образом , оптимальная прямая задается уравнением Метод точечной квадратичной аппроксимации (метод наименьших квадратов ) не предполагает , что мы должны приближать экспериментальные данные лишь с помощью прямых линий . Во многих экспериментах связи могут быть нелинейными , и было бы глупо искать для этих задач линейные соотношения . Пусть , например , мы раб о таем с радиоактивным материалом . Тогда выходными данными у являются показания счетчика Гейгера в различные моменты времени t . Пусть наш материал представляет собой смесь двух радиоактивных веществ , и мы знаем период полураспада каждого из них , но не знаем, в каких пропорциях эти вещества смешаны . Если обозначить их количества через С и D , то показания счетчика будут вести себя подобно сумме двух экспонент , а не как прямая : . (1) На практике , поскольку радиоактивность измеряется дискретно и через различные промежутки времени , показания счетчика не будут точно Рис . 1. Аппроксимация прямой лин ией. соответствовать (1). Вместо этого мы имеем серию показаний счетчика в различные моменты времени , и (1) выполняется лишь приближенно : Если мы имеем более двух показаний , m> 2, то точно разрешить эту систему относительно C и D практически невозможно . Но мы в состоянии получить приближенное решение в смыс ле минимальных квадратов . Ситуация будет совершенно иной , если нам известны количества веществ C и D и нужно отыскать коэффициенты и . Это нелинейная задача наименьших квадратов , и решить ее с ущественно труднее . Мы по– прежнему будем минимизировать сумму квадратов ошибок , но сейчас она уже не будет многочленом второй степени относительно и , так что приравнивание нулю производной не будет давать линейных уравнений для отыскания оптимальных решений . Глава 1. Метод наименьших квадратов 1.1. Задача наименьших квад ратов Задача наименьших квадратов заключается в минимизация евклидовой длины вектора невязок Ax-b . Теорема 1. Пусть А – m n – матрица ранга k , представленная в виде A = HRK T (2) где H ортогональная m m матрица ; R – m n – матрица вида , (3) где : R 11 – k x k – матрица ранга k ; K – ортогональная k x k – матрица . Определим вектор (4) и введем новую переменную . (5) Определим как единственное решение системы R 11 y 1 = g 1 . Тогда : 1. Все решения задачи о минимизации Ax-b имеют вид , где y 2 произвольно. 2. Любой такой вектор приводит к одному и тому же вектору невязки . (6) 3. Для нормы r справедливо 4. Единственным решением минимальной длины является вектор Доказательство . В выражении для квадрата нормы невязки заме ним A на HRK T в соответствии с (2) и умножая на ортогональную матрицу H T (умножение на ортогональную матрицу не меняет евклидову норму вектора ) получим (7) Далее из (3) и (5) следует , что . Из (4) следует Подставляя оба последних выражения в (7) получим Последнее выражение имеет минимальное значение при R 11 y 1 = g 1 , а в этом уравнении единственным реш ением является , так как ранг матрицы R 11 равен к . Общее решение y выражается формулой , где y 2 произвольно . Для вектора имеем , что устанавливает равенство (3). Среди векторов наименьшую длину имеет тот , для которого y 2 =0. Отсюда следует , что решением наименьшей длины будет вектор . Теорема доказана. Всякое разложение матрицы А типа (2) мы будем называть ортогональным разложением А . Заметим , что решение минимальной длины , множество всех решений и минимальное значение для задачи минимизации Ax-b определяются единственным образом . Они не зависят от конкретного ортогонального разложения . При проведении разложения необходимо приводить матрицы к диагональному виду . Дл я этого обычно используются два преобразования : Гивенса и Хаусхолдера , оставляющие нормы столбцов и строк матриц неизменными. 1.2. Ортогональное вращение Гивенса Лемма. Пусть дан 2 – вектор , причем либо .Существует ортогональная 2 2 матрица такая , что : (8) Доказательство. Положим : . Далее прямая проверка. Матрица преобразования представляет собой матрицу вращений или отражений 1.3. Ортогональное преобразование Хаусхолдера Применяется д ля преобразования матриц к диагональному виду . Матрица преобразования представляет из себя следующее выражение : , (9) или , если вектор v нормирован , т.е . и спользуется вектор единичной длины , то . В обоих случаях H – симметричная и ортогональная матрица . Покажем это : . Отсюда следует : что , т.е . симметричность и ортогональность . В комплексном случае матрица эрмитова Матрица А эрмитова ес ли она совпадает со своей комплексно сопряженной . и унитарна Матрица А унитарная если , где – сопряженная матрица. . Предположим , что дан вектор х размерности m , тогда существует матрица H такая , что , где а = + 1, при положительной первой компоненте вектора х и = – 1, при отрицательной. Доказательство. Положим действительная матрица . Любую действительную матрицу можно привести в треугольному виду Далее принимаем во внимание то , что и получаем следующее : 1.4. Сингулярное разложение матриц Пусть X – матрица данных порядка N x p , где N>p , и пусть r – ранг матрицы X . Чаще всего r=p , но приводимый ниже результат охватывает общий случай , он справедлив и при условии r

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 11 – невырожденная треугольная к к– матрица. Заметим , что выбором Н и К в уравнении (19) можно добиться , чтобы R 11 была верхней или нижней треугольной. В (19) матрица А представлена произведением A = HRK T , где R – некоторая прямоугольная матрица , ненулевые компоненты которой сосредоточены в невырожденной треугольной подматрице . Далее мы покажем , что эту невырожденную подматрицу R можно упростить далее до невырожденной диагональной матрицы . Это разложение тесно связано со спектральным разложением симметричных неотрицательно определенных матриц A T A и AA T (см . 11). Теорема 5. Пусть А – m n – матрица ранга k . Тогд а существуют ортогональная m m – матрица U , ортогональная n n – матрица V и диагональная m n – матрица S такие , что U T AV = S , A = USV T (20) Матрицу S можно выбрать так , чт обы ее диагональные элементы составляли невозрастающую последовательность ; все эти элементы неотрицательны и ровно к из них строго положительны. Диагональные элементы S называются сингулярными числами А . Докажем сперва лемму для специального случая m=n=ran kA . Лемма 2. Пусть А – n n – матрица ранга n . Тогда существует ортогональная n n – матрица U , ортогональная n n – матрица V и диагональная n n – матрица S такие , что U T AV = S , A = USV T и последовательные диагональные элементы S положительны и не возрастают. Доказательство леммы . Положительно определенная симметричная матрица A T A допускает спектральное разложение A T A = VDV T , (21) где V – ортогональн ая n n – матрица , а D – диагональная матрица , причем диагональные элементы D положительны и не возрастают . Определим S как диагональную n n – матрицу , диагональные элементы которой суть положительные квадратные корни из соответствующих диагональных элементов D . Таким образом D = S T S = S 2 , S -1 D S -1 = I . (22) Определим матрицу U = AVS -1 (23) Из (21), (22), (23) и ортогональности V следует , что U T U=S -1 V T A T AVS -1 = S -1 DS -1 =I т . е . U ортогональна . Из (23) и орт огональности V выводим USV T = AVS -1 SV T = AVV T = A Лемма доказана. Доказательство теоремы 5 . Пусть A = HRK T , где H , R , K T имеют свойства , указанные в теореме 4. Так как R 11 из (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 есть про извол . Пусть – сингулярное число А , имеющее кратность l . Это значит , что для упорядоченных сингулярных чисел найдется индекс I такой , что Положим k=min(m,n ), и пусть Q – ортогональная к к– матрица вида Здесь Р – ортогональная l l – матрица Если A = USV T – сингулярное разложение А и s i =… = s i + 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) мультипликативное неравенство – Наиболее употребимы следующие нормы для векторов : · норма суммы модулей · евклидова норма · норма максимума модуля Нормы матри ц : · · · Здесь являются сингулярными числами Сингулярным разложением произвольной m n – матрицы называется разложение вида , где U и V – ортогональные матрицы , а S – диагональная матрица с неотрицательными диагональными элементами . Диагональны е элементы S ( , i= 1,..., k , где k=min ( m,n )) называются сингулярными числами А . Это множество чисел однозначно определяется матрицей А . Число ненулевых сингуляр ных чисел равно рангу А . матрицы А ; это положительные значения квадратных корней из собственных значений матрицы А Т А (которая при невырожденной матрице А положительно определена Симметричная матрица положительно определена , если все ее собственные значения положительны . Положительно определенная матрица P облад ает также тем свойством , что для всех . , в противном случае положительно полуопределена (неотрицательно определена Симметричная матрица неотрицательно определена , если все ее собственные значения неотрицательны . Такая матрица P обладает также тем свойством , что для всех . Для произвольной m x n – матрицы А матрица симметрична и неотрицательно определена . Она положительно определена , если rankA=n . ) и поэтому имеет только вещественные собственные значения 0). Для вещественных симметричных матриц син гулярные числа равны абсолютным величинам собственных значений : . Умножение вектора х на матрицу А приводит к новому вектору Ах , норма которого может очень си льно отличаться от нормы вектора х . Область изменений может быть задана двумя числами Максимум и минимум берутся по всем ненулевым векторам . Заметим , что есл и А вырождена , то m= 0. Отношение M/m называется числом обусловленности матрицы А , (7) Рассмотрим норму обратной Обратной матрицей для квадратной невырожденной матрицы А называется такая матрица , для которой . м атрицы . Для матрицы А существует сингулярное разложение , тогда , отсюда . Аналогично для обратной матрицы и . Отсюда следу ет , что собственные числа матрицы – 1 / есть величины , обратн ые собственным числам матрицы – . При этом очевидно , что . Из последнего выражения вместе с (7) следует . Таким образо м обусловленность матрицы равна произведению нормы матрицы на норму обратной матрицы. Рассмотрим систему уравнений Ax=b , и другую систему , полученную изменением правой части : A(x+ x)=b+ b . Будем сч итать b ошибкой в b , а x соответствующей ошибкой в x , хотя нам нет необходимости считать ошибки малыми . Поскольку A( x)= b , то определения M и m немедленно приводят к неравенствам Следовательно , при m 0, Величина есть относительное изменение правой части , а величина – относительная ошибка , вызванная этим изменением . Аналогичные выкладки можно провести не только с элементами вектора правой части но и с элементами самой матрицы А и найти зависимость между относ ительным изменением элементов матрицы и относительной ошибкой вызванной этим изменением . Отсюда следует , что число обусловленности выполняет роль множителя в увеличении относительной ошибки. Приведем некоторые свойства числа обусловленности . Ясно , что M m и поэтому cond(А ) 1. Если Р – матрица перестановок Матрица перестановки - это квадратная матрица , столбцы которой получаются перестановкой столбцов единичной матрицы . Матрица перестановки ортого нальна. , то компоненты вектора Px лишь порядком отличаются от компонент вектора х . Отсюда следует , что и cond(P)= 1 . В частности cond(I)= 1. Если А умножаетс я на скаляр с , то cond(cА )= cond(А ) . Если D – диагональная матрица , то глава 2. Реализация сингулярного разложения 2.1. Алгоритмы QR – алгоритм начинается с разложения матрицы по Грамму-Шмидту , затем меняются местами сомножители : Эта м атрица подобна первоначальной , Этот процесс продолжается , причем собственные значения не изменяются : Эта формула описывает QR – алгоритм без сдвигов . Обычно время которое тратится на такой процесс пропорционально кубу размерности матрицы – n 3 . Необходимо процесс ускорить , для чего используется предварительное приведе ние матрицы А к форме Хессенберга Матрица А хессенбергова (верхняя хессенбергова ) если для j 1. Трехдиагональная матрица – это частный случай хесенберговой матрицы. . При использовании матрицы Хессенберга время п роцесса пропорционально n 2 , а при использовании трехдиагональной матрицы – n . Можно использовать другие соотношения где Q s – унитарная , а L s – нижняя треугол ьная матрица . Такой алгоритм носит название QL – алгоритма . В общем случае , когда все собственные значения матрицы различны , последовательность матриц A s имеет пределом нижнюю треугольную матрицу , диагональные элементы которой представляют собой собственные значения матрицы А , расположенные в порядке возрастания их модулей . Если матрица А имеет кратные собственные значения , то предельная матриц а не является треугольной , а содержит диагональные блоки порядка p , соответствующие собственному числу кратности p . В общем случае , наддиагональный элемент матрицы A s на s -ом шаге асимптотически равен , где k ij – пост оянная величина . Сходимость QL – алгоритма вообще говоря недостаточна . Сходимость можно улучшить , если на каждом шаге вместо матрицы A s использовать матрицу A s - k s I ( QL – алгоритм со сдвигом ). Последовательность вычислений в этом случае описывается следующими с оотношениями : которые определяют матрицу . При этом асимптотическое поведение элемента определено соотношением , а не , как прежде . Если сдвиг k s выбрать близко к величине (наименьшее собственное значение ), то в пределе внедиагональные э лементы первой строки будут очень быстро стремиться к нулю . Когда ими можно пренебречь , элемент с рабочей точностью равен , остальные являются собственными значениями оставшейся матрицы n- 1 - го порядка . Тогда , если QL – алгоритм выполнен без ускорения сходимости , то все равно , и поэтому автоматически можно выделить величину сдвига k s . Если матрица А эрмитова , то очевидно , что и все матрицы А s эрмитовы ; если А действительная и симметричная , то все Q s ортогональны и все А s действительны и симметричны. 2.2. Реализация разложения Таким образом , разложение производится в два этапа . Сначала матрица А посредством двух конечных последовательностей преобразований Хаусхолдера где , приводится к верхней двухдиагона льной форме следующего вида : Далее реализуется итерационный процесс приведения двухдиагональной матрицы J 0 к диагональной форме , так что имеет место следующ ая последовательность : где а 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 2 вычисляется следующим образом . Для первого столбца имеем : так как Таким образом , в первый столбец были введены нули и его длина не изменилась . Получим второй столбец : для третьего столбца : окончательно, Столбцы матрицы A 2 получаются вычитанием кратных вектора v 1 из столбцов A 1 . Эти кратные порождаются скалярными произведениями , а не отдельными элементами , как в гауссовом исключении. Прежде чем вводить дальнейшие нули под диагональю , преобразованием вида A 3 = A 2 Q 1 , получают нули в первой строке . Нули уже стоящие в первом столбце , не должны быть испорчены , длина первого столбца должна быть сохранена ; поэтому при внесении нулей в первую строку нельзя менять первый элемент строки , изменяем второй элемент и зануляем третий . Для матрицы большего размера на этом шаге было бы получено n – 2 нуля. Преобразование порождается первой строкой A 2 : Строка матрицы A 3 с номером i получается по формуле . Таким образом , из каждой строки A 2 вычитается на длежащее кратное . Это дает матрицу Поскольку первая компонен та нулевая , то нули первого столбца A 2 сохраняются в A 3 , Так как Q 1 ортогональная , то длина каждой строки в A 3 равна длине соответствующей строки в A 2 . Теперь можно добиться новых нулей под диагональю , не испортив полученных ранее : Поскольку ранг этой матрицы равен лишь 2, то теперь третий столбец имеет на диагонали и под диагональю элементы порядка ошибки округления . Эти элементы обозначены в матрице через 0.000, чтобы отличить их от элементов , в точности равных нулю . Если бы матрица имела полный ранг , то нужно было бы выполнить еще одно преобразование , чтобы полу ч ить нули в третьем столбце : Если бы не ошибки округлений , то в данном примере третий диагональный элемент был бы точным нулем . Элементы под диагональю во все х столбцах указаны как точные нули , потому что преобразования так и строились , чтобы получить там нули . Последнее преобразование H 3 в этом примере могло бы быть тождественным , однако тогда оно не было бы хаусхолдеровым отражением . Фактически использование H 3 попутно изменяет знак элемента – 1.080 в матрице A 4 . Получилась искомая двухдиагональная матрица , и первый этап закончен . Прямое использование ортогональных преобразований не позволяет получить какие– либо новые нули . Для общего порядка n нужно n преобр азований H и n – 2 преобразований Q , чтобы достигнуть данного места . Число преобразований не зависит от строчной размерности m , но от m зависит работа , затрачиваемая на выполнение каждого преобразования. глава 3. Использование сингул ярного разложения в методе наименьших квадратов При использовании метода сингулярного разложения ( SVD – S ingular V alue D ecomposition) мы проводим разложение для матрицы плана . При этом основное уравнение y = Xb приобретает вид y = U V T b . Отсюда следует , что коэффициенты b можно получить решая уравнение U T y = V T b . Т.е . если все j , j =1,…, n , являющиеся диагональными элементами отличны от нуля , то последнее уравнение разрешимо и , где . Однако такое решение не всегда желательно , если некоторые j малы . Для правильного использования метода SVD мы должны ввести границу отражающую точность входных данных и точность использ ованной плавающей арифметики . Всякое j , большее , чем , приемлемо , и соответствующее вычисляется по (1.20). Любое j , меньшее , чем , рассматривается как пренебрежимо малое , и соответствующему может быть придано произвольное значение . С этой произвольностью значений связана не единственность набора коэффициентов , получаемых методом наименьших квадратов . Изменения входных данных и о шибки округлений , меньшие , чем , могут привести к совершенно другому набору коэффициентов , определяемых этим методом . Поскольку обычно желательно , чтобы эти коэффициенты были по возможности малы , то полагаем =0, если j . Отбрасывание чисел j , меньших , чем , приводит к уменьшению числа обусловленности с до . Поскольку число обусловленности является множителем в увеличении ошибки , то следствием будет более надежное определение коэффициентов . Продемонстрируем использование метода на следующем примере : 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 и t 2 к линейной зависимости . Для того , чтобы исправить ситуацию можно предпринять одну из следующих мер . Во– первых , можно выбрать границу для относительной ошибки , которая бы отражала точность данных и точность арифметики . Если взять границу в интервале , то отбросим третье сингулярное число . Пр и этом получим следующие наборы коэффициентов для двойной и обычной точности : Теперь коэффициенты находятся в гораздо лучшем согласии друг с другом . Кроме то го , коэффициенты стали существенно меньше , а это значит , что не будет столь большого , как прежде , взаимного уничтожения слагаемых при вычислении квадратичного многочлена . Прогнозное значение Y (1980) будет соответственно 212910000 и 214960000. Эффект обычно й точности еще заметен , однако результаты уже не являются катастрофическими. Можно также определить набор нулевых коэффициентов , соответствующих пренебрежимо малому сингулярному числу . Вот эти коэффициенты : . Для значений t от 1900 до 1970 величина функции не превосходит 0.0017 , поэтому при любом коэффициенты можно изменить , и при этом значения , выдаваемые моделью изменятся не более чем на 0.0017 . Любой из чет ырех перечисленных нами наборов коэффициентов можно получить из другого подобным изменением. Во– вторых , можно улучшить ситуацию заменой базиса . Модели гораз до более удовлетворительны . Важно при этом то , что независимая переменная преобразуется из интервала [1900, 1970] в какой– нибудь более приемлемый интервал вроде [0, 70] или , еще лучше , [ – 3.5, 3.5]. Числа обусловленности при этом равны 5750 и 10.7 соответс т венно . последнее значение более чем приемлемо даже при счете с обычной точностью. Удобнее всего воспользоваться стандартными способами статистического анализа , т.е . матрицу плана преобразуем к стандартизованному варианту Матрица стандартизованных данных ес ть матрица наблюдений с нулевым средним и дисперсией 1. Это означает , что данные берутся в виде отклонений от среднего , которое мы считаем равным 0, вводим нормировку деля каждый член столбца матрицы на корень квадратный из суммы квадратов отклонений. Во в тором случае , после преобразования матрицы плана ее обусловленность сильно уменьшается , и , соответственно , повышается точность расчетов. Данную программу можно использовать и при решении системы линейных уравнений вместо методов Гаусса , Жордана , Холесского и пр . В приложении 2 приведен пример расчета линейной системы , которая изначально не может быть решена этими методами вследствие вырожденности матрицы коэффициентов . Тем не менее , исследуемый метод дает нам правильное решение. ЗАКЛЮЧЕНИЕ В работе описаны компьютерные методы решения задачи наименьших квадратов . Для использования данных методов составлена соответствующая программа на алгоритмическом языке FO RTRAN . Программа апробирована , результаты тестирования показывают работоспособность программы. Результаты данной разработки могут быть использованы в самых разнообразных расчетах , где необходимо провести аппроксимацию данных заданными функциями. ЛИТЕРАТУРА 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) 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) SIG MA(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) 2 90 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 T O 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=R V1(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 -.10 59419E+01 Правые корни после проверки -.1400000E+02 .1300000E+02 .1000001E+01 Видно , что правые части соответствуют начальным данным . Решение верно.

1Архитектура и строительство
2Астрономия, авиация, космонавтика
 
3Безопасность жизнедеятельности
4Биология
 
5Военная кафедра, гражданская оборона
 
6География, экономическая география
7Геология и геодезия
8Государственное регулирование и налоги
 
9Естествознание
 
10Журналистика
 
11Законодательство и право
12Адвокатура
13Административное право
14Арбитражное процессуальное право
15Банковское право
16Государство и право
17Гражданское право и процесс
18Жилищное право
19Законодательство зарубежных стран
20Земельное право
21Конституционное право
22Конституционное право зарубежных стран
23Международное право
24Муниципальное право
25Налоговое право
26Римское право
27Семейное право
28Таможенное право
29Трудовое право
30Уголовное право и процесс
31Финансовое право
32Хозяйственное право
33Экологическое право
34Юриспруденция
 
35Иностранные языки
36Информатика, информационные технологии
37Базы данных
38Компьютерные сети
39Программирование
40Искусство и культура
41Краеведение
42Культурология
43Музыка
44История
45Биографии
46Историческая личность
47Литература
 
48Маркетинг и реклама
49Математика
50Медицина и здоровье
51Менеджмент
52Антикризисное управление
53Делопроизводство и документооборот
54Логистика
 
55Педагогика
56Политология
57Правоохранительные органы
58Криминалистика и криминология
59Прочее
60Психология
61Юридическая психология
 
62Радиоэлектроника
63Религия
 
64Сельское хозяйство и землепользование
65Социология
66Страхование
 
67Технологии
68Материаловедение
69Машиностроение
70Металлургия
71Транспорт
72Туризм
 
73Физика
74Физкультура и спорт
75Философия
 
76Химия
 
77Экология, охрана природы
78Экономика и финансы
79Анализ хозяйственной деятельности
80Банковское дело и кредитование
81Биржевое дело
82Бухгалтерский учет и аудит
83История экономических учений
84Международные отношения
85Предпринимательство, бизнес, микроэкономика
86Финансы
87Ценные бумаги и фондовый рынок
88Экономика предприятия
89Экономико-математическое моделирование
90Экономическая теория

 Анекдоты - это почти как рефераты, только короткие и смешные Следующий
Живи так, как тебе хочется. Но не требуй от других жить так, как хочешь этого ты.
Anekdot.ru

Узнайте стоимость курсовой, диплома, реферата на заказ.

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

Смотрите также:


Банк рефератов - РефератБанк.ру
© РефератБанк, 2002 - 2016
Рейтинг@Mail.ru