В работе исследуется воздействие искажений матрицы линейной модели на статистические характеристики оценок по методу наименьших квадратов.
Ключевые слова: оценки МНК, линейная модель,матричные ошибки.
V.A.Galanina PhD, Associate Professor L.A. Reshetov PhD, Associate Professor M.V. Sokolovskay Senior lecturer A.E.Farafonova Master
St. Petersburg State University of Aerospace Instrumentation
Influence of matrix errors on parameter estimates by the least squares method
The paper investigates the effect of distorsions of the linear model matrix on the statistical characteristics of the least squares estimates.
Keywords: LSE,linear model,matrix errors.
Метод наименьших квадратов является наиболее эффективным и простым способом оценки неизвестных параметров матричных моделей.Достаточно указать ,что большинство алгоритмов распознавания изображений , существующих в настоящее время, используют градиентные процедуры минимизации ошибки , которые , как известно,
представляют собой одну из рекуррентных форм метода наименьших квадратов[1].Ошибки модели обычно задают в виде аддитивного случайного вектора и чаще , для упрощения , используют дискретный “белый”шум.Матрица модели при этом предполагается известной. В ряде практических задач [2] возможно появление матричных ошибок модели, что приводит к возникновению смещения оценок параметров и определение величины этих смещений является основным предметом исследования настоящей работы.
Введение. Входной сигнал в линейной модели определяют в виде суммы
y = Ax + ε , (1)
где y вектор входных данных (размер m×1) , x вектор неизвестных параметров (размер n×1), A известная матрица (размер m×n), ε вектор помехового воздействия (размер m×1). Часто предполагается , что вектор εимеет нулевое математическое ожидание (среднее значение) и ковариационную матрицу N(размер m×m). Методом наименьших квадратов принято называть следующую минимизационную процедуру
minx(y – Ax)TN-1(y – Ax) . (2)
Выполняя стандартные операции по поиску минимума квадратичной формы, получаем оценку вектора параметров x по методу наименьших квадратов(МНК)
x^ = (ATN-1A)-1ATN-1y , (3)
где AT транспонированная матрица , N-1 обратная матрица .Если помеха ε имеет многомерное нормальное распределение , то оценка (3) является также оценкой максимального правдоподобия [3].При реализации этой оценки необходимо знать ковариационную матрицу шума N. Учитывая ограниченную длину предъявляемых реализаций для большинства практических задач, наличие надёжной информации о ковариационной структуре помехи маловероятно.В подобных условиях допускают , что помеха ε “белый” шум, то есть её ковариационная матрица N = σ2I , где I единичная матрица.Оценка (3)при этом будет иметь вид
x^ = (ATA)-1ATy (4)
и имеет в статистике определённое название OLS(ordinary least square) ,а также принято её называть оценкой МНК. Мы в дальнейшем будем пользоваться этой терминалогией. Следствием этого предположения является то, что оценка (4) будет несмещённой , так как математическое ожидание M(x^) = x , а ковариационная матрицаоценки (4) равна
Kx^ = σ2(ATA)-1 . (5)
Дисперсия оценок зависит в первую очередь от дисперсии шума σ2 и углового расстояния между столбцами матрицы A .Лучшие результаты достигаются для ортогональных векторов , так как матрица ATA в этом случае диагональна.
Вычислим сумму квадратов остаточных разностей RSS(residual sum of squares). Её принимают как критерий качества процедуры минимизации квадратичной формы.Действительно, если в квадратичной форме J = (y – Ax)T(y – Ax) заменить вектор x его оценкой МНК (4), то квадратичная форма будет иметь вид
J = yTPA┴ y , (6)
где PA┴ проектор(ортогональный проектор) на подпространство векторов ортогональных столбцам матрицы A .Напомним , что проекционной матрицей(проектором ) P называют симметричную матрицу , обладающую следующим свойством P = P2 Если предполагаемая модель входного вектора верна , то . y = Ax + ε и RSS равно
J = εTPA┴ε . (7)
Нетрудно найти математическое ожидание RSS , используя известное свойство квадратичных форм εTPA┴ε = Tr(PA┴εεT), где Tr(.) след матрицы.Тогда математическое ожидание M(J) = Tr(PA┴N) или
M(J) = (m – n)σ2 . (8)
Увеличение числа столбцов n матрицы A повышает качество аппроксимации и уменьшает математическое ожидание RSS .При m = n ошибка аппроксимации M(J) равна нулю , так как аппроксимирующая кривая проходит через все отсчёты входных данных и задача аппроксимации переходит в задачу интерполяции.
До сих пор мы предполагали , что матрицы A , содержащиеся во входной реализации y = Ax + ε и в оценке x^ = (ATA)-1ATy , совпадают. Часто это не так.В силу различных причин входная реализация может иметь вид y = A~x + ε , где матрица A~ также случайна.Решение проблемы в общей постановке является сложной задачей и выходит за рамки данного исследования.В настоящей работе мы вводим следующие допущение:случайная матрица A~ и случайный вектор ε статиститчески независимы. Матрицу A~ мы можем представить в виде суммы A~ = A + δA~, а математическое ожидание случайной составляющей матрицы M(δA~) будем обозначать δA .
Основной результат. В начале более подробно рассмотрим влияние методических ошибок задания матрицы A , то есть полагаем , что A~ = A + δA и тогда y = (A + δA)x + ε. Математическое ожидание оценки МНК в этом случае равно
M(x^) = x + (ATA)-1ATδAx (9)
и смещение оценки , как мы видим , δx = M(x^) – x = (ATA)-1ATδAx .Определим также величину остаточных разностей RSS
J = (δAx+ε)TPA┴ (δAx+ε) . (10)
Найдём проекции векторов , образующих столбцы матрицы δA , на подпространство векторов, образующих столбцы матрицы A ,
QA= PAδA
и, так как PA= A(ATA)-1AT, то
QA= ASA , (11)
где SA = (ATA)-1ATδA квадратная матрица. Запишем ортогональное разложение матрицы δA в виде суммы δA = QA + QA┴ , где QA┴ = PA┴δA , а проекционная матрица PA┴ = I - PA осуществляет проектирование на подпространство векторов ортогональных столбцам матрицы A . Проекции векторов QA┴ = PA┴δA определены не единственным образом , так как любая система линейно независимых векторов , лежащая в соответствующей гиперплоскости, является подходящей.
Используя проекторы , найдём смещение оценки и математическое ожидание RSS в следующей форме
δx = SAx , (12)
M(J) = xT(QA┴)TQA┴x + (m – n)σ2 . (13)
Наблюдая полученные нами формулы (12) и (13) , мы можем заметить существенное различие воздействия аддитивных векторных помех ε и матричных помех δA .Влияние векторной помехи не связано с вектором искомых параметров x , а матричные ошибки пропорциональны элементам этого вектора.Важно также отметить , что смещение оценки определено проекцией QA= PAδA , а величина RSS проекцией QA┴ = PA┴δA .Следствием этого является то, что возможна ситуация , когда во входной реализации присутствует помеха , а оценка МНК остаётся несмещённой.Например, если δA ≠ 0 , но SA = 0.Тогда смещение δx = 0 , а величина M(J) принимает максимально возможное значение xTδATδAx + (m – n)σ2.
Представляет интерес определение предельно допустимого смещения оценки параметров,если предположить , что ковариационная матрица оценки Kx^ = σ2(ATA)-1. Сопоставим норму ковариационной матрицы Kx^ и квадрат нормы вектора смещения параметров δx .Известно ,что норма матрицы Kx^ равна
║Kx^║= σ2║(ATA)-1║ , (14)
а число обусловленности матрицы ATA
ν = ║ATA║║(ATA)-1║ . (15)
Тогда формула (14) для нормы матрицы будет иметь вид
║Kx^║= νσ2║ATA║-1 . (16)
Квадрат нормы вектора смещения δx , применяя эвклидову норму вектора, должен быть определён следующим образом
║δx║2 = δxTδx = xTSATSAx . (17)
Если мы хотим , чтобы ошибки оценки , вызванные смещением, не превышали ошибок от аддитивного шума, то мы должны обеспечить выполнение неравенства
║δx║2 ≤ ║Kx^║ . (18)
Учитывая формулы (16) и (17), получаем
xTSATSAx ≤ νσ2║ATA║-1 . (19)
Для достижения окончательного результата нам необходимо определить характер возмущения матрицы A . Если рассмотреть самый простой вид синхронных возмущений δA = αA, где α - скаляр , то матрица SA = αI . В итоге мы имеем
α2 ≤ νσ2║ATA║-1║x║-2. (20)
Из неравенства (20) следует , что допустимый уровень возмущений матрицы A тем больше ,чем больше зашумлённость входной реализации σ2,больше число обусловленности матрицы ATA и меньше её норма .Важно заметить ,что норма неизвестного вектора параметров ║x║тоже определяет допустимый уровень возмущений.
Это обстоятельство является главной особенностью воздействия матричных помех и, как следствие , не позволяет построить априорные оценки погрешности измерений.
Перейдём к анализу воздействия случайной составляющей матрицы δA~ на статистические характеристики оценки МНК при предположении, что математическое ожидание δA~ равно нулю M(δA~) = 0. Из допущения о независимости вероятностных характеристик случайного вектора ε и случайной матрицы δA~ следует , что ковариационная матрица оценки содержит два слагаемых: ковариационную матрицу σ2(ATA)-1 ,которая определена вектором ε (см.формулу5), и ковариационную матрицу M(SAxxTSAT), порождённую флуктуациями матрицы A .Здесь SA квадратная случайная матрица, а действие M(.) означает,как и прежде, операцию вычисления математического ожидания.Таким образом , мы можем записать ковариационную матрицу оценки в следующей форме
Kx^ = σ2(ATA)-1 + M(SAxxTSAT) . (21)
Из этой формулы следует , что межэлементная корреляция оценок отдельных составляющих вектора x, при наличии возмущений матрицы A , может иметь существенные вариации.Причиной этого является тот факт , что ранг матрицы M(SAxxTSAT) равен единице , а ранг матрицы (ATA)-1 равен n .Тогда возможные изменения баланса между флуктуациями аддитивного вектора ε и случайными изменениями матрицы δA~ приведут к значительным изменениям характера ковариационной матрицы.Наиболее наглядно это видно на примере синхронных возмущений δA~ = αA, где α случайное число .В результате формула для ковариационной матрицы будет иметь вид
Kx^ = σ2(ATA)-1 + M(α2)xxT . (22)
Если математическое ожидание случайной величины αравно нулю M(α) = 0 , то M(α2) это дисперсия α. Обозначая её Dα , получаем Kx^ = σ2(ATA)-1 + DαxxT . В этой формуле соотношение между параметрами σ2 и Dα регулирует вид ковариационной матрицы.Если ввести определение относительной величины u = Dα / σ2 , то мы можем сейчас представить формулу (22) в окончательной редакции
Kx^ = σ2[(ATA)-1 + uxxT] (23)
и использовать её для расчётов в поясняющем примере.
Пример.Пусть в пяти точках 0,1,2,3,4 координатной оси Оx задан набор входных данных, которыми заполняют вектор y . Мы хотим аппроксимировать эти данные прямой,используя метод наименьших квадратов.Матрица A имеет размер 5×2. Первый столбец должен содержать единицы и он имитирует смещение прямой , а во второй столбец помещены отсчёты линейно нарастающей функции.Таким образом, транспонированная матрица AT будет имееть такой вид
1
1
1
1
1
0
1
2
3
4
Допустим , что дисперсия помехи ε равна единице σ2 = 1 , а корреляция между элементами случайного вектора ε отсутствует.Воспользуемся формулой(23)для вычисления ковариационной матрицы оценки МНК параметров аппроксимирующей прямой.Если изменения матрицы A отсутствуют (u= 0) и помехи вызываются лишь вектором ε , то ковариационная матрица оценки Kx^ равна
0.6
-0.2
-0.2
0.1
На главной диагонали матрицы размещены значения дисперсий оценки соответствующих параметров прямой.Мы видим,что дисперсия оценки первого элемента вектора неизвестных параметров прямой x1, который определяет смещение этой прямой, равна 0.6, а дисперсия второго параметра x2 , задающего наклон прямой,равна 0.1 . Внедиагональные элементы матрицы характеризуют степень статистической связи оценок.После нормировки элементов ковариационной матрицы среднеквадратическими отклонениями, мы получаем корреляционную матрицу Rx^
1
-0.816
-0.816
1
Внедиагональные элементы матрицы -0.816 являются коэффициентом корреляции оценок параметров прямой.Он демонстрирует высокую степень связи оценок МНК этих параметров при отсутствии матричных возмущений , причём знак коэффициента корреляции отрицательный. Далее мы приведём результаты вычисления коэффициента корреляции оценок с учётом случайных изменений матрицы A .
Для конкретизации положим , что величина первого параметра в векторе x равна единице x1 = 1 и проведём вычисления для трёх значений параметра x2 .В нижеследующей таблице приведены результаты расчёта коэффициента корреляции оценок параметров для шести значений величины u= Dα / σ2 в диапазоне от 0 до 0,5 с шагом 0,1 .
u
0
0.1
0.2
0.3
0.4
0.5
x2 = 0.5
-0.816
-0.507
-0.287
-0.126
0
0.101
x2 = 1
-0.816
-0.267
0
0.167
0.283
0.369
x2 = 2
-0.816
0
0.236
0.37
0.46
0.526
Из данных таблицы следует, что присутствие даже малых матричных возмущений приводит к изменению величины и знака коэффициента корреляции оценок.Скорость изменений зависит от соотношения компонент вектора параметров.При увеличении нормы вектора x , если оставаться в рамках модели синхронных возмущений, модуль коэффициента корреляции будет стремиться к единице.
Выводы.Традиционная модель помехового воздействия с аддитивным вектором шума далеко не всегда соответствует реальной модели входных данных.Присутствие матричных возмущений может резко изменить статистические характеристики оценок МНК , а это не позволяет правильно спрогнозировать условия численного эксперимента.Проведённый нами математический анализ явления показывает пути построения простых правил обнаружения неадекватности модели. Платой за эту информацию будет ,как обычно, дополнительное увеличение объёма выборочных данных.
Библиографический список 1. Галанина В.А.,Решетов Л.А.,Соколовская М.В. Взаимосвязь байесовских оценок и оценок с линейными ограничениями //Сборник докладов конференции”Моделирование и ситуационное управление качеством сложных систем”СПб,ГУАП,2020,С.25-28 2. Альберт А. Регрессия,псевдоинверсия и рекуррентное оценивание.М.:Наука,1977.223c. 3. Галанина В.А.,Решетов Л.А.,Соколовская М.В. Рекуррентная форма информационнго байесовского решения для плохо обусловленных систем //Сборник статей “Волновая электроника и инфокоммуникационные системы”СПб,ГУАП,2020,ч.1,С.188-194
Сборник докладов “Моделирование и ситуационное управление качеством сложных систем” Вторая Всеросийская научная конференция 14-22 апреля 2021г.
Опубликовано в ноябре 2021г. стр.33-38
-----------------------------------------------
Решетов Леонид
lr11111@yandex.ru
