МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ
Математический факультет
Кафедра прикладной математики
ДИПЛОМНЫЙ ПРОЕКТ
сингулярное разложение в линейной задаче метода наименьших квадратов
Заведующий кафедрой прикладной
математики
Исполнил:
Научный руководитель
Владикавказ 2002
СОДЕРЖАНИЕ
TOC \o “1-2” ВВЕДЕНИЕ PAGEREF _Toc11128507 \h 3
Глава 1. Метод наименьших квадратов PAGEREF _Toc11128508 \h 7
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. Реализация сингулярного разложения PAGEREF _Toc11128515 \h
25
2.1. Алгоритмы PAGEREF _Toc11128516 \h 25
2.2. Реализация разложения PAGEREF _Toc11128517 \h 27
2.3. Пример сингулярного разложения PAGEREF _Toc11128518 \h 29
глава 3. Использование сингулярного разложения в методе наименьших
квадратов PAGEREF _Toc11128519 \h 33
ЗАКЛЮЧЕНИЕ PAGEREF _Toc11128520 \h 38
ЛИТЕРАТУРА PAGEREF _Toc11128521 \h 39
ПРИЛОЖЕНИЕ 1. Исходные тексты программы PAGEREF _Toc11128522 \h 40
ПРИЛОЖЕНИЕ 2. контрольный пример PAGEREF _Toc11128523 \h 45
ВВЕДЕНИЕ
Метод наименьших квадратов обычно используется как составная часть
некоторой более общей проблемы. Например, при необходимости проведения
аппроксимации наиболее часто употребляется именно метод наименьших
квадратов. На этом подходе основаны: регрессионный анализ в статистике,
оценивание параметров в технике и т.д.
Большое количество реальных задач сводится к линейной задаче наименьших
квадратов, которую можно сформулировать следующим образом.
Пусть даны действительная m(n–матрица A ранга k(min(m,n) и
действительный m–вектор b. Найти действительный n–вектор x0,
минимизирующий евклидову длину вектора невязки Ax–b.
Пусть y – n–мерный вектор фактических значений, x – n–мерный вектор
значений независимой переменной, b – коэффициенты в аппроксимации y
линейной комбинацией n заданных базисных функций (:
.
. Покажем это. Возведем в квадрат выражение для е:
.
=0
.
Следует отметить, что последнее выражение имеет в определенной степени
формальный характер, т. к. решение нормальных уравнений, как правило,
проводится без вычисления обратной матрицы (метод Крамера) такими
методами как метод Гаусса, Холесского и т. д.
по результатам измерений. Мы получаем переопределенную систему:
или Xb=y. Нам понадобится матрица XTX и обратная к ней:
Метод точечной квадратичной аппроксимации (метод наименьших квадратов)
не предполагает, что мы должны приближать экспериментальные данные лишь
с помощью прямых линий. Во многих экспериментах связи могут быть
нелинейными, и было бы глупо искать для этих задач линейные соотношения.
Пусть, например, мы работаем с радиоактивным материалом. Тогда выходными
данными у являются показания счетчика Гейгера в различные моменты
времени t. Пусть наш материал представляет собой смесь двух
радиоактивных веществ, и мы знаем период полураспада каждого из них, но
не знаем, в каких пропорциях эти вещества смешаны. Если обозначить их
количества через С и D, то показания счетчика будут вести себя подобно
сумме двух экспонент, а не как прямая:
. (1)
На практике, поскольку радиоактивность измеряется дискретно и через
различные промежутки времени, показания счетчика не будут точно
Рис. 1. Аппроксимация прямой линией.
, и (1) выполняется лишь приближенно:
Если мы имеем более двух показаний, m>2, то точно разрешить эту систему
относительно C и D практически невозможно. Но мы в состоянии получить
приближенное решение в смысле минимальных квадратов.
Ситуация будет совершенно иной, если нам известны количества веществ C и
D и нужно отыскать коэффициенты ( и (. Это нелинейная задача наименьших
квадратов, и решить ее существенно труднее. Мы по–прежнему будем
минимизировать сумму квадратов ошибок, но сейчас она уже не будет
многочленом второй степени относительно ( и (, так что приравнивание
нулю производной не будет давать линейных уравнений для отыскания
оптимальных решений.
Глава 1. Метод наименьших квадратов
1.1. Задача наименьших квадратов
Задача наименьших квадратов заключается в минимизация евклидовой длины
вектора невязок (( Ax-b ((.
Теорема 1. Пусть А – m(n–матрица ранга k, представленная в виде
A=HRKT (2)
где H ортогональная m(m матрица; R – m(n–матрица вида
, (3)
где: R11 – kxk–матрица ранга k; K – ортогональная kxk–матрица.
Определим вектор
(4)
и введем новую переменную
. (5)
как единственное решение системы R11y1=g1. Тогда:
, где y2 произвольно.
. (6)
Доказательство. В выражении для квадрата нормы невязки заменим A на HRKT
в соответствии с (2) и умножая на ортогональную матрицу HT (умножение
на ортогональную матрицу не меняет евклидову норму вектора) получим
(7)
Далее из (3) и (5) следует, что
.
Из (4) следует
Подставляя оба последних выражения в (7) получим
имеем
,
. Теорема доказана.
Всякое разложение матрицы А типа (2) мы будем называть ортогональным
разложением А. Заметим, что решение минимальной длины, множество всех
решений и минимальное значение для задачи минимизации ((Ax-b((
определяются единственным образом. Они не зависят от конкретного
ортогонального разложения.
При проведении разложения необходимо приводить матрицы к диагональному
виду. Для этого обычно используются два преобразования: Гивенса и
Хаусхолдера, оставляющие нормы столбцов и строк матриц неизменными.
1.2. Ортогональное вращение Гивенса
.Существует ортогональная 2(2 матрица такая, что:
(8)
Доказательство. Положим:
.
Далее прямая проверка.
Матрица преобразования представляет собой матрицу вращений
или отражений
1.3. Ортогональное преобразование Хаусхолдера
, (9)
. В обоих случаях H – симметричная и ортогональная матрица. Покажем
это:
.
, где
а ( = +1, при положительной первой компоненте вектора х и = –1, при
отрицательной.
действительная матрица. Любую действительную матрицу можно привести в
треугольному виду
и получаем следующее:
1.4. Сингулярное разложение матриц
Пусть X – матрица данных порядка Nxp, где N>p, и пусть r – ранг матрицы
X. Чаще всего r=p, но приводимый ниже результат охватывает общий случай,
он справедлив и при условии r
имеет форму, указанную правой частью (4). Теорема доказана.
Подматрицу [R:T] из правой части (18) можно теперь преобразовать к
компактной форме, требуемой от матрицы R из теоремы 2. Это
преобразование описывает следующая лемма.
Лемма 1. Пусть [R:T] – к(к–матрица, причем R имеет ранг к. Существует
ортогональная n(n–матрица W такая, что
– нижняя треугольная матрица ранга к.
Доказательство леммы вытекает из теоремы 3, если отождествить величины
n, k, [R:T], W из формулировки леммы с соответствующими величинами m, n,
AT, QT теоремы 3.
Используя теорему 3 и лемму 1 можно доказать следующую теорему.
Теорема 4. Пусть А – m(n–матрица ранга к . Найдутся ортогональная
m(m–матрица Н и ортогональная n(n–матрица К такие, что
(19)
причем R11 – невырожденная треугольная к(к–матрица.
Заметим, что выбором Н и К в уравнении (19) можно добиться, чтобы R11
была верхней или нижней треугольной.
В (19) матрица А представлена произведением A=HRKT, где R – некоторая
прямоугольная матрица, ненулевые компоненты которой сосредоточены в
невырожденной треугольной подматрице. Далее мы покажем, что эту
невырожденную подматрицу R можно упростить далее до невырожденной
диагональной матрицы. Это разложение тесно связано со спектральным
разложением симметричных неотрицательно определенных матриц ATA и AAT
(см. 11).
Теорема 5. Пусть А – m(n–матрица ранга k. Тогда существуют ортогональная
m(m–матрица U, ортогональная n(n–матрица V и диагональная m(n–матрица S
такие, что
UTAV=S, A=USVT (20)
Матрицу S можно выбрать так, чтобы ее диагональные элементы составляли
невозрастающую последовательность; все эти элементы неотрицательны и
ровно к из них строго положительны.
Диагональные элементы S называются сингулярными числами А. Докажем
сперва лемму для специального случая m=n=rankA.
Лемма 2. Пусть А – n(n–матрица ранга n. Тогда существует ортогональная
n(n–матрица U, ортогональная n(n–матрица V и диагональная n(n–матрица S
такие, что UTAV=S, A=USVT и последовательные диагональные элементы S
положительны и не возрастают.
Доказательство леммы. Положительно определенная симметричная матрица ATA
допускает спектральное разложение
ATA=VDVT, (21)
где V – ортогональная n(n–матрица, а D – диагональная матрица, причем
диагональные элементы D положительны и не возрастают. Определим S как
диагональную n(n–матрицу, диагональные элементы которой суть
положительные квадратные корни из соответствующих диагональных элементов
D. Таким образом
D=STS=S2, S-1DS-1=I. (22)
Определим матрицу
U=AVS-1 (23)
Из (21), (22), (23) и ортогональности V следует, что
UTU=S-1VTATAVS-1=S-1DS-1=I т.е. U ортогональна. Из (23) и
ортогональности V выводим USVT=AVS-1SVT=AVVT=A Лемма доказана.
Доказательство теоремы 5. Пусть A=HRKT, где H, R, KT имеют свойства,
указанные в теореме 4. Так как R11 из (19) – невырожденная треугольная
к(к–матрица, то согласно лемме 2 , можно написать
(24)
– невырожденная диагональная матрица, диагональные элементы которой
положительны и не возрастают. Из (24) следует, что матрицу R в уравнении
(19) можно записать в виде
(25)
где:
– ортогональная m(m–матрица;
– ортогональная n(n–матрица;
– ортогональная m(n–матрица;
Теперь, определяя U и V формулами
(26)
заключаем из (24) – (26), что A=USVT, где U, S, V имеют свойства,
указанные в формулировке теоремы 5. Это завершает доказательство.
Заметим, что сингулярные числа матрицы А определены однозначно, в то
время, как в выборе ортогональных матриц U, V есть произвол. Пусть ( –
сингулярное число А, имеющее кратность l. Это значит, что для
упорядоченных сингулярных чисел найдется индекс I такой, что
Положим k=min(m,n), и пусть Q – ортогональная к(к–матрица вида
.
1.6. Число обусловленности
Некоторые вычислительные задачи поразительно чувствительны к изменению
данных. Этот аспект численного анализа не зависит от плавающей
арифметики или выбранного алгоритма.
Например:
Найти корни полинома: (x-2)2=10-6
Корни этого уравнения есть 2+10-3 и 2-10-3. Однако изменение свободного
члена на 10-6 может вызвать изменение в корнях, равное 10-3.
Операции с матрицами, как правило, приводят к решению систем линейных
уравнений. Коэффициенты матрицы в правой части системы линейных
уравнений редко известны точно. Некоторые системы возникают из
эксперимента, и тогда коэффициенты подвержены ошибкам наблюдения.
Коэффициенты других систем записываются формулами, что влечет за собой
ошибки округлений. В связи с этим необходимо знать, как влияют ошибки в
коэффициентах матрицы на решение. Именно для этого вводится понятие
обусловленности матрицы.
. Для более подробного описания числа обусловленности нам понадобится
понятие нормы в пространстве векторов и матриц.
, удовлетворяющий следующим условиям:
;
.
, удовлетворяющий условиям 1 – 3 для нормы вектора:
;
Наиболее употребимы следующие нормы для векторов:
Нормы матриц:
.
Умножение вектора х на матрицу А приводит к новому вектору Ах, норма
которого может очень сильно отличаться от нормы вектора х.
Область изменений может быть задана двумя числами
Максимум и минимум берутся по всем ненулевым векторам. Заметим, что если
А вырождена, то m=0. Отношение M/m называется числом обусловленности
матрицы А,
(7)
.
.
0
b
d
f
h
j
l
®
°
a
ae
ae
e
e
h
h
?e
46hjlnpI?
T
V
?
?
?
’
?
1/4
i
?
o
oe
o
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
/o{o{o{ouopoioipopo
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
:u?e?eu?eaeYOYae?ae1/2??ae?eae?ae ’?ae?aeu?aeoaehOae
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
?
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
Wc h
h
wbwdwOwOw-x
x”x(x:x
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
?e–eaeaNaEeaEUeae h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
?
Y Y
Y"Y,Y0Y8Y
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
Wc h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
WcH* h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
ju
h
h
h
h
Wc h
h
h
h
jU
h
h
j„
h
h
ju
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
ji
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
. Таким образом обусловленность матрицы равна произведению нормы
матрицы на норму обратной матрицы.
Следовательно , при m(0,
– относительная ошибка, вызванная этим изменением. Аналогичные
выкладки можно провести не только с элементами вектора правой части но и
с элементами самой матрицы А и найти зависимость между относительным
изменением элементов матрицы и относительной ошибкой вызванной этим
изменением. Отсюда следует, что число обусловленности выполняет роль
множителя в увеличении относительной ошибки.
глава 2. Реализация сингулярного разложения
2.1. Алгоритмы
Этот процесс продолжается, причем собственные значения не изменяются:
Эта формула описывает QR–алгоритм без сдвигов. Обычно время которое
тратится на такой процесс пропорционально кубу размерности матрицы – n3.
Необходимо процесс ускорить, для чего используется предварительное
приведение матрицы А к форме Хессенберга а также используется алгоритм
со сдвигом. Форма Хессенберга представляет из себя верхнюю треугольную
матрицу (верхняя форма Хессенберга) у которой сохранена одна диагональ
ниже главной, а элементы ниже этой диагонали равны нулю. Если матрица
симметрична, то легко видеть, что матрица Хессенберга превращается в
трехдиагональную матрицу. При использовании матрицы Хессенберга время
процесса пропорционально n2, а при использовании трехдиагональной
матрицы – n.
Можно использовать другие соотношения
где Qs – унитарная, а Ls – нижняя треугольная матрица. Такой алгоритм
носит название QL–алгоритма.
кратности p.
, где kij – постоянная величина. Сходимость QL–алгоритма вообще говоря
недостаточна. Сходимость можно улучшить, если на каждом шаге вместо
матрицы As использовать матрицу As-ksI (QL–алгоритм со сдвигом).
Последовательность вычислений в этом случае описывается следующими
соотношениями:
, и поэтому автоматически можно выделить величину сдвига ks.
Если матрица А эрмитова, то очевидно, что и все матрицы Аs эрмитовы;
если А действительная и симметричная, то все Qs ортогональны и все Аs
действительны и симметричны.
2.2. Реализация разложения
, приводится к верхней двухдиагональной форме следующего вида:
а Si и Ti – диагональные матрицы.
где
.
и ничего не добавляет.
было QR–преобразованием со сдвигом, равным s.
Обычный QR–алгоритм со сдвигом можно записать в следующем виде:
. Параметр сдвига s определяется собственным значением нижнего минора
(размерности 2(2) матрицы Mi. При таком выборе параметра s метод
обладает глобальной и почти всегда кубичной сходимостью.
2.3. Пример сингулярного разложения
,
Преобразованная матрица A2 вычисляется следующим образом. Для первого
столбца имеем:
Таким образом, в первый столбец были введены нули и его длина не
изменилась. Получим второй столбец:
для третьего столбца:
окончательно,
Столбцы матрицы A2 получаются вычитанием кратных вектора v1 из столбцов
A1. Эти кратные порождаются скалярными произведениями, а не отдельными
элементами, как в гауссовом исключении.
Прежде чем вводить дальнейшие нули под диагональю, преобразованием вида
A3=A2Q1, получают нули в первой строке. Нули уже стоящие в первом
столбце, не должны быть испорчены, длина первого столбца должна быть
сохранена; поэтому при внесении нулей в первую строку нельзя менять
первый элемент строки, изменяем второй элемент и зануляем третий. Для
матрицы большего размера на этом шаге было бы получено n–2 нуля.
Преобразование порождается первой строкой A2:
Строка матрицы A3 с номером i получается по формуле
.
. Это дает матрицу
нулевая, то нули первого столбца A2 сохраняются в A3, Так как Q1
ортогональная, то длина каждой строки в A3 равна длине соответствующей
строки в A2.
Теперь можно добиться новых нулей под диагональю, не испортив полученных
ранее:
Поскольку ранг этой матрицы равен лишь 2, то теперь третий столбец имеет
на диагонали и под диагональю элементы порядка ошибки округления. Эти
элементы обозначены в матрице через 0.000, чтобы отличить их от
элементов, в точности равных нулю. Если бы матрица имела полный ранг, то
нужно было бы выполнить еще одно преобразование, чтобы получить нули в
третьем столбце:
Если бы не ошибки округлений, то в данном примере третий диагональный
элемент был бы точным нулем. Элементы под диагональю во всех столбцах
указаны как точные нули, потому что преобразования так и строились,
чтобы получить там нули. Последнее преобразование H3 в этом примере
могло бы быть тождественным, однако тогда оно не было бы хаусхолдеровым
отражением. Фактически использование H3 попутно изменяет знак элемента
– 1.080 в матрице A4.
Получилась искомая двухдиагональная матрица, и первый этап закончен.
Прямое использование ортогональных преобразований не позволяет получить
какие–либо новые нули. Для общего порядка n нужно n преобразований H и
n–2 преобразований Q, чтобы достигнуть данного места. Число
преобразований не зависит от строчной размерности m, но от m зависит
работа, затрачиваемая на выполнение каждого преобразования.
глава 3. Использование сингулярного разложения в методе наименьших
квадратов
При использовании метода сингулярного разложения (SVD – Singular Value
Decomposition) мы проводим разложение для матрицы плана. При этом
основное уравнение y=Xb приобретает вид y=U(VTb. Отсюда следует, что
коэффициенты b можно получить решая уравнение UTy=(VTb. Т.е. если все
(j, j=1,…,n, являющиеся диагональными элементами ( отличны от нуля, то
последнее уравнение разрешимо и
.
=0, если (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 и t2 к линейной
зависимости. Для того, чтобы исправить ситуацию можно предпринять одну
из следующих мер.
, то отбросим третье сингулярное число. При этом получим следующие
наборы коэффициентов для двойной и обычной точности:
Теперь коэффициенты находятся в гораздо лучшем согласии друг с другом.
Кроме того, коэффициенты стали существенно меньше, а это значит, что не
будет столь большого, как прежде, взаимного уничтожения слагаемых при
вычислении квадратичного многочлена. Прогнозное значение Y(1980) будет
соответственно 212910000 и 214960000. Эффект обычной точности еще
заметен, однако результаты уже не являются катастрофическими.
, и при этом значения, выдаваемые моделью изменятся не более чем на
0.0017(. Любой из четырех перечисленных нами наборов коэффициентов можно
получить из другого подобным изменением.
Во–вторых, можно улучшить ситуацию заменой базиса. Модели
гораздо более удовлетворительны. Важно при этом то, что независимая
переменная преобразуется из интервала [1900, 1970] в какой–нибудь более
приемлемый интервал вроде [0, 70] или, еще лучше, [–3.5, 3.5]. Числа
обусловленности при этом равны 5750 и 10.7 соответственно. последнее
значение более чем приемлемо даже при счете с обычной точностью.
Удобнее всего воспользоваться стандартными способами статистического
анализа, т.е. матрицу плана преобразуем к стандартизованному варианту
Матрица стандартизованных данных есть матрица наблюдений с нулевым
средним и дисперсией 1. Это означает, что данные берутся в виде
отклонений от среднего, которое мы считаем равным 0, вводим нормировку
деля каждый член столбца матрицы на корень квадратный из суммы квадратов
отклонений.
Во втором случае, после преобразования матрицы плана ее обусловленность
сильно уменьшается, и, соответственно, повышается точность расчетов.
Данную программу можно использовать и при решении системы линейных
уравнений вместо методов Гаусса, Жордана, Холесского и пр. В приложении
2 приведен пример расчета линейной системы, которая изначально не может
быть решена этими методами вследствие вырожденности матрицы
коэффициентов. Тем не менее, исследуемый метод дает нам правильное
решение.
ЗАКЛЮЧЕНИЕ
В работе описаны компьютерные методы решения задачи наименьших
квадратов. Для использования данных методов составлена соответствующая
программа на алгоритмическом языке FORTRAN. Программа апробирована,
результаты тестирования показывают работоспособность программы.
Результаты данной разработки могут быть использованы в самых
разнообразных расчетах, где необходимо провести аппроксимацию данных
заданными функциями.
ЛИТЕРАТУРА
Беллман Р. Введение в теорию матриц. -М.: Наука, 1969, 368с.
Гантмахер Ф.Р. Теория матриц. -М.: Наука, 1988, 548с.
Ланкастер П. Теория матриц. -М.: Наука, 1982, 387с.
Лоусон Ч., Хенсон Р. Численное решение задач наименьших квадратов. М.:
Статистика, 1979, 447с
Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980
Мэйндоналд Дж. Вычислительные алгоритмы в прикладной статистике. М.:
Финансы и статистика, 1988, 350с
Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980, 454с
Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке АЛГОЛ. Линейная
алгебра, М.: Машиностроение, 1976, 390с
Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. -М.:
Физматгиз, 1963, 536с.
Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических
вычислений. М.: Мир, 1980, 279с
Харебов К.С. Компьютерные методы решения задачи наименьших квадратов и
проблемы собственных значений. Владикавказ.: Изд-во СОГУ, 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) 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
Видно, что правые части соответствуют начальным данным. Решение верно.
.
– сопряженная матрица.
, i=1,...,k, где k=min(m,n)) называются сингулярными числами А. Это
множество чисел однозначно определяется матрицей А. Число ненулевых
сингулярных чисел равно рангу А.
.
симметрична и неотрицательно определена. Она положительно определена,
если rankA=n.
.
Матрица перестановки - это квадратная матрица, столбцы которой
получаются перестановкой столбцов единичной матрицы. Матрица
перестановки ортогональна.
для j1. Трехдиагональная матрица – это частный случай
хесенберговой матрицы.
PAGE
PAGE 2
2
/
1
2
2
2
1
2
2
/
1
2
2
2
1
1
)
(
=
,
)
(
v
v
v
s
v
v
v
c
?
?
?
Нашли опечатку? Выделите и нажмите CTRL+Enter