Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 19
Текст из файла (страница 19)
е. критерий (7.2) охватывает оба случая. Однако сам метод наименьших квадратов (МНК), выраженный критерием (7.2), создан Лежандром в1805 году как алгебраическая процедура именно для несовместных систем иподтвержден как статистическая процедура Гауссом в 1809 году. МНК какалгебраическая процедура проиллюстрирован выше с помощью рис. 7.1(a),а как статистическая процедура — с помощью рис. 7.1(b). Замечательно, чтообе процедуры имеют одинаковые решения, T т.
е. алгебраически эти решения= I (см. рис. 7.1(b)) совпадают,эквивалентны и при E {v} = 0 и E vvпоэтому можно говорить о едином МНК-решении x̂.МНК-решение x̂ всегда существует как решение нормальных уравненийAT Ax̂ = AT z,(7.3)x̂ = A+ z + (I − A+ A)y(7.4)выражается формулойчерез произвольный вектор y ∈ Rn , где A+ — псевдообратная матрицадля матрицы A, и единственно тогда и только тогда, когда A+ A = I, чторавносильно условию, что только нулевой вектор составляет ядро (нульпространство) матрицы A, т. е. при rank A = n.Условие rank A = n, называемое условием полного столбцового рангаматрицы A, обусловливает случай m ≥ n, что при m > n означает переопределенную систему полного ранга в (7.1).
Этот типичный для практики1097 Ортогональные преобразованияслучай ниже и рассматривается, при этом из (7.3), (7.4) следует x̂ = A+z иA+ = (AT A)−1AT .Замечание 7.1. Слагаемое x̂0 , A+z в (7.4) есть единственное МНКрешение с минимальной нормой, называемое нормальным псевдорешением.Оно ортогонально второму слагаемому в (7.4), т. е. A+z ⊥ (I − A+ A)y, илежит в пространстве строк матрицы A, т. е. x̂0 ∈ R(AT ) (подробнее см.подразд. 10.4).Таким образом, типичный для практики случай имеет формальное решение x̂ = x̂0 = (AT A)−1AT z, и вычислительная задача наименьших квадратовзаключается в его эффективном отыскании. Об этом подробнее см. разд.
11.7.3Ортогональные матрицы и наименьшие квадратыВ рассматриваемой задаче о наименьших квадратах имеем критерийJ(x) = kz − Axk2,A(m, n),m ≥ n,rank A = n.(7.5)Пусть T , T (m, m), есть матрица некоторого ортогонального преобразования. В силу свойства C (см. подразд. 7.1) запишемJ(x) = kT (z − Ax)k2 = k(T z) − (T A)xk2.(7.6)При таком представлении видно, что минимум критерия J(x), равныйJ(x̂), не зависит от T . Этом фактом можно воспользоваться, т.
е. матрицуT можно выбрать так, что (T A) приобретает привлекательную для вычислений форму. Действительно, в подразд. 7.4 и 7.7 мы покажем, как можновыбрать T , чтобы преобразованная матрица имела вид R }nTA =(7.7)0 }m− nс верхнетреугольным блоком R, rank R = n.Если соответственно этому вектор T z разбить на блоки, т.
е. записатьz1 } nTz =,(7.8)z2 } m − nто J(x) от (7.6) приводится к видуJ(x) = kz1 − Rxk2 + kz2 k2.110(7.9)7.4 Преобразование ХаусхолдераПриведение критерия наименьших квадратов к виду (7.9) позволяетзаключить, что искомый вектор x̂, минимизирующий этот критерий, должен удовлетворять уравнениюRx̂ = z1 ,(7.10)которое легко решается обратной подстановкой (см. подразд. 7.6), и крометого,min J(x) = J(x̂) = kz2 k2.(7.11)В вычислительном отношении эти результаты гораздо более элегантны,чем неразумная трата сил на решение нормальных уравнений (7.3).
Однаковажнее всего то, что решение, использующее ортогональные преобразования(соотношения (7.7), (7.8), (7.10) и (7.11)), менее чувствительны к погрешностям, вызванным ошибками округления в компьютере. Это видно хотя быиз того, что выражение (7.7) влечет равенствоRT R = (T A)T (T A) = AT A,которое означает, что R является квадратным корнем из матрицы (AT A)системы нормальных уравнений (7.3). Следовательно, при решении системы(7.10) вдвое более эффективно используется разрядная сетка компьютера,чем при решении системы (7.3)1.7.4Преобразование ХаусхолдераПреобразования Хаусхолдера суть матричные представления, которыесоответствуют геометрическому понятию отражения [97, 15]. Пусть заданнекоторый ненулевой вектор u, который мы называем направляющим вектором.
Подпространство, ортогональное ему, есть гиперплоскость U⊥. Есливзять произвольный вектор y, то можно отразить его от U⊥ , в точностисоблюдая законы обычного оптического отражения от плоского зеркала(рис. 7.2).Обозначим отраженный вектор yr . Поскольку положение гиперплоскостиU⊥ не зависит от длины направляющеговектора, пронормируем его, т. е.образуем орт û = u/kuk. Проекция (y u) вектора y на прямую, задаваемуюнаправлением u, равна (y T û)û. Следовательно,y = (y u) + v, v ⊥ u, v ∈ U⊥.(7.12)1Представление в компьютере квадрата a2 любого действительного числа a требует удвоенной разрядности мантиссы, т. е.
счет по уравнению (7.10) равносилен счету с удвоенной разрядностью мантиссычисел по уравнению (7.3).1117 Ортогональные преобразованияyr − yU⊥vyry−(y u)(y u)u0Рис. 7.2. Геометрия преобразования Хаусхолдера. Задача 1 (прямая): даны векторыu и y, найти вектор yr , отраженный от гиперплоскости U⊥Отраженный вектор yr , как видно из рис. 7.2, имеет разложениеyr = −(y u) + v, v ⊥ u, v ∈ U⊥(7.13)с той же составляющей v, которая ортогональна вектору u, но с проекцией−(y u), которая (в силу знака −) направлена противоположно проекции(y û) вектора y на направление u.
Исключая v из (7.12) и (7.13), находимyr = y − 2(y u) = (I − βuuT )y = Tuy,(7.14)где β , 2/kuk2 = 2/uT u. Матрица Хаусхолдера Tu , (I − βuuT ), в вычислениях явно не участвующая, имеет фундаментальное значение для приложений в силу своих замечательных свойств.Свойство 1. Tu = TuT , т. е. Tu — симметрическая матрица.Свойство 2. Tu2 = I, т. е. Tu — идемпотентная матрица. Это легкопродемострировать алгебраически разложением матрицы Tu2 или же геометрически по рис.
7.2 как двукратное отражение вектора y относительно U⊥.Свойство 3. Если u(j) = 0, то (Tuy)(j) = y(j), т. е. если j-я компонентавектора u — нулевая, то Tu оставляет j-ю компоненту вектора y неизменной.Свойство 4. Если u ⊥ y, то Tuy = y.Свойство 5.Tuy = y − γu,112γ , 2y T u/uT u = βy T u.(7.15)7.4 Преобразование ХаусхолдераСвойство 5 — важное с практической точки зрения. Формирование матрицы Tu в качестве множителя для y потребовало бы на порядок большевычислений, чем того требует прямое вычисление Tuy по свойству 5. Этотакже означает, что не нужно тратить память для хранения Tu, что наиболее существенно проявляется при больших m.Триангуляризация матрицы методом Хаусхолдера. Обратимся косновному применению ортогональных преобразований.
Для этого решимзадачу, обратную к той, что рассмотрена выше, а именно: дан вектор y идано желаемое расположение отраженного вектора yr , — найти направление u такое, что Tuy = (s, 0, . . . , 0)T (рис. 7.3). Из свойства C, подразд. 7.1,норма (евклидова длина) вектора y не изменяется при ортогональном преобразовании, следовательно, определим ее какσ , kTuyk = |s| = (y T y)1/2.(7.16)Направление u может быть получено из свойства 5 (уравнение (7.15)), т.е.u = const · (y − se1 ).(7.17)Этот результат приводит к следующему свойству.Свойство 6.
Пусть s = − sgn [y(1)]σ, где sgn [·] — функция знака,(1,x ≥ 0,sgn [x] =−1, x < 0,и элементы вектора u определены выражением (7.17), т. е. u(1) = y(1) − s,u(i) = y(i), i > 1. Тогда Tuy = se1 и β , 2/uT u = −1/(su(1)).Замечание 7.2. Геометрический смысл выражения (7.17) ясен изрис.
7.3, где видно, что вектор yr ортогонален гиперплоскости U⊥ и параллелен (коллинеарен) вектору u.Непосредственное вычисление uT u показывает, что uT u = −2su(1), приэтом знак для s выбран противоположно знаку первого элемента y(1), т. е.так, чтобы максимизировать |u(1)| и тем уменьшить относительную погрешность вычисления разности u(1) = y(1)−s. Если свойство 6 применить к матрице A, взяв в качестве y ее первый столбец, то это даст первый шаг, которыйпослужит основой приведения матрицы к верхнетреугольному виду. Повторение таких действий шаг за шагом позволит осуществлять верхнюю триангуляризацию любой заданной матрицы A.1137 Ортогональные преобразованияe2aU⊥yuaryre10Рис. 7.3. Геометрия преобразования Хаусхолдера. Задача 2 (обратная): даны векторыy и yr, найти вектор u, задающий отражающую гиперплоскость U⊥ ; здесь yr = se1 =T= s 0···0Лемма 7.1.
Пусть дана матрица A(m, n). Тогда существует ортогональное преобразование Хаусхолдера Tu такое, что1n−1z}|{ z }| {1{s à .Tu A =m−10 (7.18)Замечание 7.3. Скаляр s и матрица à в (7.18) вычисляются непосредственно по данным в матрице A; s — по выражению (7.16) и свойству 6,а à — по свойству 5, (7.15). Первый столбец, который уместно назватьведущим столбцом в преобразовании Хаусхолдера, используют как вектор yв задаче 2 (см. рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
