7.7. ПОИСК ВЗВЕШЕННОГО ПСЕВДОРЕШЕНИЯ


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

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

Вместо «решения» нерешаемых уравнений их следует разделить на информативную первую и малоинформативную вторую части

Pθ = R,
P=
 P1 
 P2 
R=
 R1 
 R2 


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

θ = θ0 + W(P1W)+(R1 – P1θ0).


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

Попытка найти псевдорешение в лоб, по формуле, приводит к необходимости псевдообращать произведение P1W, что плохо по двум причинам.

Во-первых, умножение на W нарушает симметрию матрицы P=ZZT и неоправданно сужает выбор численных средств.

Во-вторых, умножение P1 на матрицу весов способно ухудшить и без того плохую обусловленность задачи: за качество взвешенных оценок приходится платить серьезным увеличением нагрузки на вычислитель.

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

Предпосылки к построению численного алгоритма. На первом этапе цель вычислений может состоять в формировании наращиваемого левым углом невырожденного квадратного блока A симметричной матрицы P, сопровождаемом перестановкой строк P, R (и столбцов P, иначе она потеряет симметрию). Перестановка исключает из рассмотрения плохо обусловленные уравнения, группируемые своими коэффициентами в нижней части матриц, так что

P =
AB
C
,
R =
 R1 
 R2 


Наращиваемый блок нужно представить в виде произведения двух сомножителей A=LLT, где L – нижнетреугольная матрица. Если преобразование доведено до конца, в том смысле, что блок A охватывает собой всю матрицу P, то алгоритм заканчивается стандартно решением двух систем с треугольными матрицами

LX=R, LTθ=X,


что соответствует формуле θ = P–1R.

Интрига вычислительного метода продолжает развиваться дальше, если размер хорошо обусловленного, подчеркнем, блока A уступит размерности задачи. Наступает черед псевдообращения произведения P1W, где согласно принятым обозначениям P1 = (A|B), матрицу весов также сепарируем на блоки W=diag(W1,W2), так что

P1W = (LLT|B)W = L(LTW1|L–1BW2).


произведение LTW1 – невырожденная верхнетреугольная матрица.

Предположим, что мы отыскали ортогональную матрицу H, например, методом Хаусхолдера, такую, что умножение на нее аннулирует остаточный блок L–1B W2 , т. е. (P1WH)H–1 = L(G|0)H–1, причем G – теперь уже невырожденная нижнетреугольная матрица. Тогда псевдообращение произведения P1W = (P1WH)H–1 сведется к обычной инверсии и перестановкам невырожденных компонент L, G, H–1, (P1W)+ = H(G–1|0)TL–1. Искомое решение предстанет в виде

θ = θ0 + Pw+(R1 – P1θ0) = θ0 + WH[G–1L–1(R1 – P1θ0);0].


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

Метод Холецкого, нацеленный на такого рода нужды, порождает не один, а семейство алгоритмов, объединенных общей идеей, так что алгоритм предстоит не столько выбрать, сколько под конкретную потребность создать. В нашем случае итогом вычислений должна стать, в более строгом толковании целей, не треугольная, а трапециевидная матрица разложения P, несущая в себе помимо L информацию о L–1B, где L–1 в явном виде не вычисляется. Для пущей экономии в таких случаях матрицу P наращивают дополнительно вектором, подлежащим умножению на L–1. Обычно это R, но здесь еще и Pθ0 или разность R – Pθ0.

Распространенная модификация метода Холецкого формирует нижнетреугольную матрицу L на месте P, столбец за столбцом, см. рис. 7.3.

Рис. 7.3Рис. 7.4


Такой близорукий подход не позволяет планировать перестановки, столь важные при решении вырожденных задач. Далее будет рассмотрен иной вариант, который занят итерационным уточнением всей диагонали L. Тогда на каждом шаге рекурсии мы сможем, очевидно, опереться на любую строку и столбец P. Но и это еще не предел изворотливости, к которой нас настоятельно подталкивает сложность проблемы. Столь же итерационно, забегая вперед, можно рассчитывать добавочный вектор в P, помещаемый обычно в нижнюю строку. Секрет прост, его элементы, образующие L–1R, дают прогноз падения невязки МНК.

Горизонтальные стрелки рис. 7.4, обозначающие заявленные прогностические действия вдоль диагонали, следует отразить вертикально вниз на дополнительную строку и тогда картина вычислений будет полна.

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

Lkk = (Pkk – Σj=1:k–1Lkj2)½, Lik=(Pik – Σj=1:k–1LijLkj)/Lkk,


где k=2..n, старт k=1 с L11 = √P11 и Li1 = Pi1/L11, i=k+1..n+1.

Рис. 7.5 поясняет содержимое трех блоков трапециевидной структуры, которая имеет высоту матрицы P, расширенной строкой RT.


Рис. 7.5


Наращиваемая вниз лента столбцов позволяет итерационно вычислять квадраты положительных диагональных элементов, пусть на первом шаге Lii2 = Pii (совпадают с диагональю) и далее Lii2=Lii2 – Lik–12 для всех i=k..n. Забегающий вперед расчет синхронизован с коррекцией строки RT так, что элементы Ri = (Ri – Rk–1Li k–1)/Lii.

Вычисленные впрок, они дают прогноз падения невязки МНК, помогающий выбрать и переставить в позицию k ведущие строку и столбец P и элемент R.

Для продолжения итераций нужно сбрасывать делители при всех несостоявшихся конкурентах элементов Lkk и Rk, т. е. Rj = RjLjj для j=k+1..n. Этой репродуктивной фазы можно избежать, более аккуратно обращаясь с памятью вычислителя.

Начальное значение максимума квадрата невязки Δk = ||Y–Zθ(k)||2 несложно посчитать как Δ0 = ||Y||2, его можно вынести как дополнительный элемент расширенной диагонали P и уменьшать на выбранную величину Δk = Δk–1 – Rk2.

Заключительная фаза. Если процесс триангуляризации благополучно завершился вычислением разложением P=LLT, по месту хранения участвующего в итерациях вектора R оказывается произведение L–1R, так что совсем несложно будет решить систему LTθ = L–1R и получить ответ невырожденной задачи. Другое дело, если по той или иной причине триангуляризация останавливается, выдавая обширную информацию в блоках трапециевидной матрицы T=(LT|L–1B). Иногда это выгодно делать, чтобы повысить значимость вектора притяжения. Смысл последующих численных операций состоит в том, чтобы, обрезав ортогональными преобразованиями с матрицей H взвешенную при помощи W трапецию TW=(LTW1|L–1BW2), превратить ее в нижний треугольник G в TWH = (G|0).

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

H = (E – 2U1U1T/(U1TU1)) (E – 2U2U2T/(U2TU2)) . ... . (E – 2UkUkT/(UkTUk))


где Ui – опорные векторы преобразований Хаусхолдера, алгоритм составления которых подробно описан в [17 , с. 58].

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

θ = θ0 + Pw+(R1 – P1θ0) = θ0 + WH[G–1L–1(R1 – P1θ0); 0].


Учитывая, что (P1W)+ = (LTW)+ = (TW)+L–1, итоги можно подвести иначе, по формуле

θ = θ0 + Pw+(R1 – P1θ0) = θ0 + W(TW)+[L–1(R1 – P1θ0; 0],


привлекая для псевдоинверсии взвешенной трапециевидной матрицы более простой в алгоритмическом отношении метод Гревилля, а именно, ту его ветвь, которая обрабатывает заведомо зависимые столбцы прямоугольного «туловища», поскольку инверсия «угла» не составляет проблем.

Рассматриваемый подход вместе с контролируемым при помощи весового коэффициента прямым методом накопления данных

Pk = wPk–1 + (1–w)ZkZkT, Rk = wRk–1 + (1–w)Zkyk,


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

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

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

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

Процессор вычислительной машины не следует перегружать операциями с несоразмерными числами. За стремительное уменьшение квадратичного критерия мы расплачиваемся риском потерять точность вычислений. На этом перепутье, к сожалению, нет единственно верной стратегии, здесь возникает вопрос назначения ведущей строки, который сам по себе может составить предмет отдельного исследования. Итерация дает прогноз на один шаг вперед, а локально хороший путь может оказаться ущербным с точки зрения совокупного эффекта. Сходные проблемы в смежной области интегрирования дифференциальных уравнений привели к возникновению численных методов типа алгоритмов Рунге-Кутта.

Rambler's Top100