7.6. ПОИСК ОБЩЕГО ПСЕВДОРЕШЕНИЯ


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

Напомним формулы метода окаймления

P =
A
 b 
 bT 
 c 
P–1 =
 A–1+ppT/d 
 –p/d 
 –pT/d 
 1/d 
p = A–1b,d=cbTp


Метод Гревилля работает с каймой в виде столбца (или строки) и обобщает предыдущий алгоритм на случай любых матриц

P =
A
 b 
P+ =
 A+aTp 
 aT 
p = A+b, d = b–Ap, a = dT/dTd,


старт с d=b, если d=0, то a=0. По ходу счета, если d=0, то a=pTA+/(1+pTp).

Смысл переменных обоих алгоритмов позволяет планировать обработку данных. В самом деле, малый делитель d означает близость к вырожденности матрицы P, этот признак позволяет сортировать столбцы (и строки, если нужно) с целью повысить эффективность процедуры. В методе Гревилля вектор Ap является ортогональной проекцией b на гиперплоскость, образованную ранее обработанными столбцами, входящими в A. Отношение норм векторов d и b дает нам в руки синус угла раствора столбца P по отношению к указанной гиперплоскости, он не превышает единицы. Следовательно, всегда можно назначить не относительный, а абсолютный порог грубости псевдоинверсии.

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

В таком случае рекурсивный алгоритм поиска общего псевдорешения системы уравнений идентификации

θ = θ0 + P+(R–Pθ0)


сводится к виду

θk = [θk–1pq; q], θ0 = [θ0k–1; q0], размерность k=1..n,


до тех пор, пока ранг наращиваемой левым верхним углом A матрицы P растет, имеем p=A–1b, d = cbTp, q=(rpTRk–1)/d, в противном случае нижние строки матрицы P, отвечающие плохо обусловленной части уравнений, игнорируются, а столбцы домолачиваются алгоритмом Гревилля p=A+b, q = (q0 + pTk–1 – θ0k–1))/ (1+pTp).

Обозначим текущую сумму квадратов невязок Δk = ||Y–Zθ(k)||2, рекурсивно вычисляемый вектор θ(k) дополнен до полной размерности нулями, она убывает до тех пор, пока вычисляемый ранг P нарастает

Δk = Δk–1q/d, Δ0 = ||Y||2.


В вычислительной математике сложились шаблоны – процедуры диагонализации и триангуляризации матриц. Именно универсальность является их слабой стороной. Учитывая специфику вырожденных задач, триангуляризация и сходные процессы должны зависеть от правой части решаемых уравнений. Обратим внимание на то, что падение оптимизируемого критерия на шаге рекурсивного алгоритма существенно зависит от отношения q/d. Этот факт заставляет по-новому взглянуть на распространенные алгоритмы выбраковки почти зависимых строк и столбцов P на основе всего лишь малой величины d. Путь безопасных вычислений и путь получения оценки с хорошими аппроксимирующими качествами едва ли не диаметрально противоположны. Именно малые по абсолютной величине делители d могут приводить к значительному снижению невязки. В данном случае алгоритм оказался поставщиком важной информации, помогающей выбрать опорный элемент.

При поиске более сильных в вычислительном отношении продолжений, рассмотрим формулы, построенные при помощи сингулярного разложения матрицы Z=UDVT системы уравнений идентификации.

В вырожденных задачах

D =
 σk 
0
0
0
Z+ = V
σk–10
00
UT,


где σk = diag1 σ2 …) – матрица ненулевых сингулярных чисел Z, выстроенных в порядке убывания; UTU=E, VVT=E. Разделим окаймляющие матрицы на блоки U=[Uk Un], V=[Vk Vn] в соответствии с делением D.

Если размер прямоугольной матрицы Z внушает сомнение, вместо нее можно использовать симметричное представление P=ZTZ=UDUT системы нормальных уравнений метода наименьших квадратов Pθ = R.

Формулу обобщенного псевдорешения можно записать в двух, казалось бы, эквивалентных вариантах

θ = θ0 + Z+(Y–Zθ0) или θ = Z+Y+(E–Z+Z)θ0.


В практике вычислительной математики двух равноценных путей не бывает, ее мощь, собственно, и состоит в выискивании рационального направления приложения усилий процессора. Если ранг матрицы Z невысок, первый путь связан с меньшим объемом вычислений. Второй путь годится для решения невырожденных и слабо вырожденных систем. Отсюда получаются две ветви алгоритма оценки параметров.

Размерность σk меньше n/2Размерность σk больше n/2
Общее псевдорешение Zθ = Y
θ = θ0 + Vkk–1 UkTY – VkTθ0)θ=V[σk–1UkTY;VnTθ0]
Общее псевдорешение Pθ = R
θ = θ0 + Ukk–1 UkTY – UkTθ0)θ=U[σk–1UkTY;UnTθ0]


Rambler's Top100