УДК 519.612.4
Д.В. Бутенина
кандидат ф-м.наук, доцент
В.А. Галанина
кандидат технических наук, доцент
Л.А. Решетов
кандидат технических наук, доцент
М.В. Соколовская
cт. преподаватель
Санкт- Петербургский государственный университет аэрокосмического приборостроения
Байесовское решение рандомизированной системы линейных уравнений
Рассмотрено решение системы линейных уравнений со случайной матрицей. В пределах гауссовской модели случайных параметров найдены байесовские оценки матрицы системы уравнений. Выполнено численное моделирование совместных оценок в частном случае одностолбцовой матрицы системы.
Ключевые слова: система линейных уравнений, случайная матрица, байесовская оценка.
D.V. Butenina
PhD, Associate Professor
V.A .Galanina
PhD, Associate Professor
L.A. Reshetov
PhD, Associate Professor
M.V. Sokolovskay
Senior lecturer
St. Petersburg State University of Aerospace Instrumentation
BAYESIAN SOLUTION OF A RANDOMIZED SYSTEM OF LINEAR EQATIONS
The solution of a system of linear equations with a random matrix is considered. Within Gaussian model of random parameters, Bayesian estimates of the matrix of a system of equatins are found. Numerical simulation of joint estimates in a particular case of a one-colomn matrix of thr system is performed.
Keywords: system of linear equations, random matrix, Bayesian estimate.
В большинстве задач математической физики допускается существование линейной модели анализируемого поля. В рамках этой модели поле представляется в виде суммы, возможно бесконечной, произведений некоторых функций (спектральные функции) и чисел, которые принято называть спектральными коэффициентами. Желательно, чтобы система функций, по которой производится разложение поля, была бы ортонормированной и полной. Кроме того, выбранная система функций должна быть согласована с особенностями поставленной физической задачи. Например, в задаче рассеяния внешнего поля малыми многослойными частицами для частиц, состоящих из вложенных сфер, наиболее адекватным является использование системы сферических функций, а для софокусных сфероидов целесообразно использовать систему сфероидальных функций [1]. К сожалению, часто при решении обратной задачи, то есть при восстановлении спектральных коэффициентов по результатам измерения поля на некотором множестве дискретных точек, вид спектральных функций неизвестен заранее. Такая задача, конечно же, требует матричной формулировки и формально описывается как решение системы линейных уравнений со случайной правой частью и случайной матрицей системы.
Понятно, что поиск большого числа неизвестных возможен только при наличии априорных данных и должен основываться на процедурах регуляризации. В настоящей работе мы используем один из вариантов статистической регуляризации, а именно байесовский подход. Основным результатом работы мы считаем получение в явном виде совместных оценок векторных и матричных параметров, определяющих систему линейных уравнений, а также обсуждение свойств этих оценок.
Рассмотрим решение системы рандомизированных (стохастических) линейных уравнений
Ax= y, A є Hmn , xє Rn , y є Rm , (1)
у которых матрица A и вектор y являются случайными величинами, Hmn – пространство прямоугольных вещественнозначных матриц размерности m×n, Rn – n-мерное евклидово пространство. В настоящей статье мы полагаем, что m≥ n и матрица A имеет полный столбцовый ранг, то есть столбцы матрицы A образуют линейно независимую систему векторов.
Для матрицы A = (a1│a2│…│an) є Hmn, состоящей из столбцов ai, введём вектор - столбец
A(v)= vec (A),
который получается, если каждый последующий столбец матрицы A расположить ниже предыдущих.
Всюду в дальнейшем предполагается, что случайные векторы xє Rn и A(v) є Rmn определены априорными распределениями вероятности p(x) и p(A). Обозначим через x0 и A0априорные значения x и A,соответственно.Решается задача: на основе реализации случайного вектора y построить такие оценки x и A, которые минимизировали бы невязку
‹ Ky-1(Ax– y),(Ax– y) ›, (2)
где Ky = K(y│x, A) – ковариационная матрица помех y , и наименее уклонялись бы от априорных значений x0 и A0.Здесь ‹ c,b › = cTb– скалярное произведение векторов c и b, cT- означает операцию транспонирования вектора.
При наличии априорных данных для получения оптимальных оценок параметров x и A стохастического уравнения (1) обычно используют байесовское правило. Тогда оценки x и A могут быть получены, если максимизировать апостериорную плотность вероятности, вычисляемую по формуле
p(x,A│y) = p(y│x,A)p(x,A) / p(y) . (3)
Мы ограничиваем себя случаем статистической независимости случайных величин x и A
p(x,A│y) = p(y│x,A)p(x)p(A) / p(y) . (4)
В практических задачах помехи y и параметры x , A часто имеют гауссовское распределение и для получения байесовских оценок достаточно провести минимизацию по x и A следующего выражения
T(x,A) = ‹ Ky-1(Ax – y),(Ax – y) › + ‹ Kx-1(x – x0),(x – x0) › + ‹ KA-1(A(v) – A0(v)),(A(v) – A0(v) ›, (5)
где Ky = K(y│x,A), Kx = K(x), KA = K(A(v)) – ковариационные матрицы, которые полагаются в данной работе положительно определёнными. Подобные оценки принято называть оценками по методу наименьших квадратов. Необходимые и достаточные условия минимума T(x,A) определены уравнениями
ƌT(x,A)/ƌx = 0 , ƌT(x,A)/ƌA(v) = 0
и приведены в статье [2] в следующем виде
(A*)T Ky-1A*x* - (A*)T Ky-1y + Kx-1(x* – x0) = 0, (6)
(Ky-1A*x*(x*)T)(v) – (Ky-1y(x*)T)(v) + KA-1((A*)(v) – (A0)(v)) = 0. (7)
Мы вводим обозначения x*и A*для значений параметров x и A , которые удовлетворяют уравнениям (6) и(7) и являются ,таким образом , байесовскими оценками.
1) Оценка параметра x
Если норма матрицы (A*)TKy-1A*для любого случайного значения матрицы A*ограничена, то из уравнения (6) получаем оценку
x* = ((A*)TKy-1A* + Kx-1)-1((A*)T Ky-1y + Kx-1x0) . (8)
Если же матрица A нам известна, то, заменяя в выражении (8) матрицу A*на матрицу A, мы получим хорошо известную оценку для коэффициентов регрессии [3]
x1* = (ATKy-1A + Kx-1)-1(AT Ky-1y + Kx-1x0). (9)
Ковариационная матрица оценки x1*имеет вид
K(x1*) = (ATKy-1A + Kx-1)-1ATKy-1A(ATKy-1A + Kx-1)-1 .
Если матрица ATKy-1Aневырожденная, то
K(x1*) = (I + (ATKy-1A )-1Kx-1)-1(ATKy-1A)-1(I + (ATKy-1A)-1Kx-1)-1 , (10)
где I– единичная матрица соответствующего размера. Благодаря матричному множителю
(I + (ATKy-1A)-1Kx-1)-1 достигается уменьшение нормы ковариационной матрицы K(x1*) оценки x1*при вводе априорных данных x0. Основным недостатком оценки (9) является её смещение при фиксированном значении случайной величины x. Действительно, условное математическое оценки x1*равно
M(x1*) = (ATKy-1A + Kx-1)-1(AT Ky-1Ax + Kx-1x0)
и M(x1*) ≠ x, если x ≠ x0.По этой причине значительно более распространенной оценкой является несмещённая оценка , которую можно получить при Kx-1= 0
x2* = (ATKy-1A )-1AT Ky-1y , (11)
хотя она требует достаточно хорошей обусловленности матрицы ATKy-1A.
Возвращаясь к оценке общего вида (8), можно заметить, что присутствие случайной матрицы A*приводит к появлению дополнительного смещения оценки x*из-за несоответствия истинного значения матрицы A и её оценки A*. Кроме того , число обусловленности матрицы (A*)TKy-1A*становится случайной величиной и это обстоятельство предъявляет определённые требования к качеству оценок матрицы A*и априорных данных A0 .
2) Оценка параметра A.
Для получения из уравнения (7) оценки матрицы A в явном виде необходимо принять дополнительные предположения относительно ковариационных матриц Ky и KA. В настоящей работе мы полагаем, что ковариационная матрица KA= K1× K2это результат кронекерова (прямого)произведения матриц K1є Hnn и K2є Hmm .Матрица K1описывает межстолбцовое вероятностное взаимодействие (ковариация столбцов) матрицы A ,а матрица K2– матрица ковариации строк. Отметим некоторые свойства кронекерова произведения матриц .
- Матрица A×B является обратимой тогда и только тогда, когда A и B являются обратимыми, и тогда
(A×B)-1= A-1× B-1
- Пусть A,B и C это матрицы размерности m×n, n×p , p×q ,соответственно. Тогда [4]
(ABC)(v)= (CT×A)B(v).
Используя эти свойства прямого произведения матриц можно доказать следующий результат.
Теорема. Пусть ковариационные матрицы Ky, Kx, KA являются положительно определёнными матрицами и случайные параметры x , A стохастического уравнения Ax = y, а также вектор правой части уравнения y, имеют гауссовское распределение. Тогда, если ковариационная матрица KA= K1× K2, где K1є Hnn и K2є Hmm, и при этом справедливо равенство K2= αKy, α > 0, то оценка матрицы A по методу наименьших квадратов будет иметь вид
A* = ( y(x*)T + α-1A0K1-1 )( x*(x*)T + α-1K1-1)-1, (12)
где x*- оценка случайного параметра x, A0- априорное значение случайного параметра A.
Доказательство. Применяя первое свойство кронекерова произведения матриц, мы можем найти матрицу KA-1= K1-1× K2-1. По определению ковариационная матрица K1 является симметричной матрицей и, следовательно, KA-1= (K1T)-1× K2-1. Тогда справедливо равенство KA-1((A*)(v)– (A0)(v)) = ((K1T)-1× K2-1)((A*)(v)– (A0)(v)) и, учитывая второе свойство кронекерова произведения, получаем
((K1T)-1× K2-1)((A*)(v) – (A0)(v)) = (K2-1(A* – A0)K1-1)(v) . (13)
Очевидно, что благодаря этому соотношению уравнение (7) примет следующий вид
(Ky-1A*x*(x*)T)(v) – (Ky-1y(x*)T)(v) + (K2-1(A* – A0)K1-1)(v) = 0. (14)
Так как выражения A + B = 0 и A(v)+ B(v)= 0 для матриц A,B одной размерности эквивалентны , то формулу (14) можно переписать в форме матричного уравнения
Ky-1A*x*(x*)T – Ky-1y(x*)T + K2-1(A* – A0)K1-1 = 0.
По условию теоремы K2= αKy. Учитывая это равенство, получаем
Ky-1A*x*(x*)T – Ky-1y(x*)T + α-1Ky-1(A*– A0)K1-1 = 0. (15)
Умножая левую и правую части уравнения (15) на Kyи выделяя слагаемые, содержащие матрицу A*, находим
A*( x*(x*)T + α-1K1-1) = y(x*)T + α-1A0K1-1 .
Так как матрица K1обратима и α > 0, то для любых конечных значений параметра α обратима матрица x*(x*)T + α-1K1-1. В конечном итоге получаем искомый результат
A*= ( y(x*)T + α-1A0K1-1)( x*(x*)T + α-1K1-1)-1 . ■
Наблюдая совместные оценки параметров (8) и (12), можно заметить, что они нелинейные и смещённые. Если допустить, что значение параметра x нам известно, то оценка параметра A
A1* = ( yxT + α-1A0K1-1 )( xxT + α-1K1-1)-1 (16)
будет линейной, но почти всегда смещённой оценкой. Действительно, если вычислить условное математическое ожидание
M(A1*) = (AxxT+ α-1A0K1-1)(xxT + α-1K1-1)-1, (17)
которое мы получаем путём осреднения оценки по распределению аддитивной помехи, то видно, что M(A1*) ≠ A , если A ≠ A0.В общем случае размерность вектора оцениваемых параметров существенно превышает размерность вектора правой части уравнения y. Для получения приемлемых качественных показателей оценки далее мы рассматриваем самый простой случай одностолбцовой матрицы A (n = 1). При этом вектор xи матрица K1являются скалярными величинами, а число оцениваемых параметров лишь на единицу превышает размерность вектора y.
Полагая, что оценка вектора x* (в данном случае x – число) никогда не равна нулю x*≠ 0, оценка вектора A* допускает простое квазиоптимальное представление
A*= y(x*)-1, (18)
которое справедливо асимптотически при отсутствии шума.
Подобная оценка не содержит априорных данных A0и, в ряде случаев, может быть несмещённой. Именно эта оценка (18), благодаря её простоте, и, сопутствующая ей оценка
x* = ((A*)TKy-1A* )-1(A*)T Ky-1y , A* ≠ 0 , (19)
являются предметом дальнейшего исследования, хотя мы понимаем, что отказ от априорных данных, пусть и частичный, приводит к повышенному влиянию аддитивных шумов. Часто ковариационная матрица шума имеет вид Ky= DyI , где Dy – дисперсия. Тогда
x*= ((A*)TA*)-1(A*)T y . (20)
Найдём область применения совместной оценки x*= ((A*)TA*)-1(A*)T y и матрицы A*= y(x*)-1 .Исходными данными для этих оценок являются вектор правой части yи априорное значение матрицы A0.Подставим A0 в формулу (20)
x*= ((A0)TA0)-1(A0)T y .
Если x и A это выборочные значения искомых параметров, то условное математическое ожидание оценки равно
M(x*) = ((A0)TA0)-1(A0)TAx. (21)
Запишем матрицу A в виде A= A0+ δA и тогда смещение оценки c= M(x*) – x будет определяться выражением
c= (│δA│/│A0│) x cosφ, (22)
где φ - угол между векторами δA и A0, │.│ - длина вектора (эвклидова норма вектора). Согласно (22), если A ≠ A0, смещение оценки x*равно нулю для ортогональных векторов δA и A0.Относительное смещение Δ = c/x = = (│δA│/│A0│) cosφ определяет область применения оценки x*через введение границы величины модуля │Δ│.
Исследование смещения оценки A*= y(x*)-1 представляет собой значительно более сложную задачу. Здесь мы примем следующее допущение: будем считать, что аддитивный шум очень мал и им можно пренебречь. Тогда условное математическое ожидание оценки A*, после осреднения по распределению аддитивного шума, будет равно
M(A*) = A(1+Δ)-1. (23)
Следует отметить, что величина смещения оценки M(A*) – A существенно зависит от вида выбранной векторной нормы.
Результаты численного эксперимента.
В качестве матрицы A мы использовали одностолбцовую матрицу, каждый элемент которой представляет собой отсчёты некоторой функции угла θ и 0○ ≤ θ ≤ 90○ .
Предполагается, что априорное значение матрицы A0описывает точки, распределённые на окружности, то есть каждый элемент этой матрицы равен единице. Таким образом, можно говорить об оценке отсчётов функций, геометрически мало отличающихся от окружности. Мы выбрали в качестве альтернативной фигуры эллипс и моделировали его отсчёты по формуле
r= b( 1- e2sin2θ )-1/2,
где e= (1 – b2/a2)1/2 – эксцентриситет эллипса, a и b – большая и малая полуоси эллипса, соответственно. При моделировании мы проводили вычисления для четырёх значений угла θ, который изменялся каждый раз на 30○, а размерность матрицы A была, следовательно,4×1. В таблице 1 мы приводим результаты вычисления элементов матрицы A для трёх отношений a/b, если малая полуось эллипса принимает значение b= 0/8.
Таблица 1 Значения элементов матрицы A для эллиптического фронта волны.
θ
0○
30○
60○
90○
A0
1
1
1
1
a/b = 1.25
0.8
0.839
0.936
1
cosφ = - 0.8
a/b = 1.55
0.8
0.866
1.067
1.24
cosφ = - 0.04
a/b = 1.85
0.8
0.882
1.168
1.48
cosφ = 0.295
Вычисления выполнялись при различных значениях параметра x, но ниже мы даём результаты счёта лишь для одного значения x= 5. В таблице 2 приведены отсчёты относительного смещения Δ оценки (20) как функции параметра a/b.
Таблица 2 Относительное смещение оценки параметра x .
a/b
1.25
1.35
1.45
1.55
1.65
1.75
1.85
│Δ│
0.106
0.072
0.04
0.007
0.024
0.054
0.082
cosφ
- 0.803
- 0.545
- 0.263
- 0.04
0.116
0.222
0.295
Расчёты показывают, что наименьшее смещение оценки параметра x достигается при a/b = 1.57, что соответствует ортогональному положению векторов A0и δA.Отсюда
следует, что для любого вектора A, который можно представить как сумму A0+ [I - A0(A0T A0)-1 A0T] z, где вектор z произволен, рассматриваемые нами оценки A и x будут несмещёнными. Согласно данным таблицы 2, относительное смещение оценки параметра x не превышает десяти процентов для эллипсов с отношением полуосей a/b, лежащем в диапазоне от 1.25 до 1.9.
Смещение оценки A*= y(x*)-1, в соответствии с формулой (23), полностью определено относительным смещением оценки x*. В таблице 3 приведены результаты оценивания вектора A для двух значений параметра a/b =1.85 и a/b= 1.25 при b = 0.8 .
Таблица 3 Смещение оценки векторного параметра A .
A │a/b =1.85
0.8
0.882
1.168
1.48
A*│a/b =1.85
0.739
0.815
1.079
1.367
1+ Δ = 1.082
A │a/b =1.25
0.8
0.839
0.936
1
A*│a/b =1.25
0.895
0.938
1.048
1.119
1+ Δ = 0.894
Мы видим, что при положительном значении относительного смещения параметра x отсчёты оценочного вектора A*уменьшаются по отношению к отсчётам истинного вектора A ,что приводит к смещению оценки. Если знак Δ отрицательный, то отсчёты вектора A*увеличиваются.
Выводы.
1 Решение системы линейных уравнений со случайной матрицей и случайной правой частью возможно, если имеется дополнительная априорная информация о вероятностных свойствах случайных параметров. Применение байесовского подхода в гауссовском предположении позволило получить оценки матрицы системы и вектора неизвестных в явном виде. Так как число оцениваемых параметров больше числа измерительных отсчётов, эти оценки следует рассматривать лишь как уточнение известных априорных данных.
2 В простейшем случае одностолбцовой матрицы системы допустимые значения матрицы, которые мы можем оценить с приемлимой точностью, концентрируются вокруг векторов, представляющих собой сумму вектора априорных данных A0и любого вектора из ортогонального к A0подпространства.
3 Численное моделироеание показало, что для A0, содержащего отсчёты на круге, вектор A, содержащий отсчёты на эллипсе, будет воспроизводится с десятипроцентной точностью, если отношение большой и малой полуосей эллипса находится в диапазоне от 1.25 до 1.9
Библиографический список
1. Фарафонов В.Г., Устимов В.И. Рассеяние света малыми многослойными частицами: обобщённый метод разделения переменных. Оптика и спектроскопия, 2018, том 124, вып.2, 255-263
2. Жуковский Е.Л., Мелешко В.И. О помехоустойчивых решениях в линейной условной алгебре. Доклады Академии наук СССР, 1978, том 241, 5, 1009-1012
3. Альберт А. Регрессия, псевдоинверсия и рекуррентное оценивание. Москва, Наука 1977, 223
Schott J.R. Matrix Analysis for Statistics. N.Y., John Wiley,
