17. Метод наименьших квадратов в многомерном случае

Будем рассматривать здесь многомерную модель как модель со многими входами и одним выходом:

Рис. 18. Многомерная система

Сигналы на входах и выходах считаем центрированными, т.е.

M[x_j(t)] = 0, j = 1, 2, ..., k; M[y(t)] = 0.

Модель ищем в виде линейной комбинации входных сигналов

Помеху считаем аддитивной

Выбираем N >> k моментов времени t_i, i = 1, 2, ..., N и в каждый из этих моментов наблюдаем все входы и выход.

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

Введем обозначения:

x_j(t_i) = x_ij; y(t_i) = y_i; (t_i) = _i.

тогда

В матричной форме

Y = XB + E,

где

Критерий приближения

Так как E = Y - X B, то

J = (Y - XB)^T(Y - XB) = (Y^T – B^TX^T)(Y – XB) =

= Y^TY – B^TX^TY – Y^TXB + B^TX^TXB

Поскольку все слагаемые в последнем уравнении скаляры, то

B^TX^TY = (B^TX^TY)^T = Y^TXB

Тогда

J = Y^TY – 2B^TX^TY + B^TX^TXB

Вектор коэффициентов B можно вычислить из условия: dJ/dB = 0:

Здесь использовано правило для двух векторов y и z:

Далее воспользуемся правилом: если Q – квадратная матрица, а z – вектор соответствующего размера, то

Для квадратной матрицы X^TX и вектора B получаем:

Тогда

dJ/dB = -2X^TY + 2X^TXB = 0

Откуда получаем т.н. нормальные уравнения МНК:

X^TXB = X^TY.

Решение нормальных уравнений

B = (X^TX)^-1 X^TY

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

Запишем матрицы системы

В одномерном случае (k = 1) положим x_i1 = x_i, тогда

y_i = b₁ x_i + _i

что совпадает с ранее полученным ранее выражением.

18. Рекуррентный метод наименьших квадратов

Пусть по N наблюдениям вычислен вектор коэффициентов

B_N = (X^TX)^-1 X^TY = P_N X^TY

где P_N = (X^TX)^-1, и пусть в момент t_N+1 получена еще одна серия отсчетов

y_N+1, x_N+1,1, x_N+1,2, ..., x_N+1,k.

Необходимо найти способ откорректировать вектор B_N по информации, полученной на (N+1)-ом шаге, т.е. найти способ вычисления по формуле

B_N+1 = B_N + f(y_N+1, x_N+1,j), j = 1, 2, ..., k.

В этом состоит рекуррентный МНК. Для коррекции вектора B_N не требуется ни обращения матриц, ни повторного формирования и решения нормальных уравнений. Метод может использоваться как при оперативной, так и при ретроспективной идентификации.

Пусть

Сформируем клеточные (или блочные) матрицы

Тогда

Воспользуемся далее формулой из леммы об обращении матриц (см. Современные методы идентификации систем. Под ред. П. Эйкхоффа. М.: Мир, 1983, стр. 391).

Лемма. Если необходимые обратные матрицы существуют, тогда

(A + BCD)^1 = A^1  A^1B(C^1 + DA^1B) ^1DA^1

Доказательство.

Если лемма справедлива, то по определению обратных матриц произведение исходной матрицы A + BCD на ее обратную должно дать единичную матрицу I:

(A + BCD)[ A^1  A^1B(C^1 + DA^1B) ^1DA^1] = AA^1 + BCDA^1 

 AA^1B(C^1 + DA^1B) ^1DA^1  BCDA^1B(C^1 + DA^1B) ^1DA^1 = I +

+ BCDA^1  B(C^1 + DA^1B) ^1DA^1  BCDA^1B(C^1 + DA^1B) ^1DA^1 =

= I+ BCDA^1  BC[C^1(C^1 + DA^1B) ^1 + DA^1B(C^1 + DA^1B) ^1]DA^1 =

= I + BCDA^1  BC[(C^1 + DA^1B)(C^1 + DA^1B) ^1]DA^1 =

= I + BCDA^1  BCDA^1 = I, ч.т.д.

Пусть теперь C = 1, B = b – матрица-столбец, D = b^T, тогда можно записать

Применяя последнее выражение, когда A = X^TX; b = x, получим:

Отсюда:

Т.о. для вычисления B_N+1 необходимо знать B_N; P_N и новую информацию y_N+1; x. Вычисления состоят только в умножении и сложении матриц, т.е. отсутствуют операции обращения или решения системы уравнений.

Обратим внимание на то, что произведение

равно значению выходного сигнала, предсказываемое моделью (построенной за N шагов) по вектору новых значений x. Если это предсказываемое значение равно измеренному значению выхода, т.е. если

x^T B_N = y_N+1,

то B_N+1 = B_N и никакого пересчета вектора B_N не требуется.

Для применения формулы пересчета вектора B_N+1 на следующем (N+2)-ом шаге потребуется, очевидно, матрица P_N+1. Эту матрицу можно получить коррекцией матрицы P_N.

Так как B_N = P_N X^TY, то выражение

B_N+1 = (X^TX + xx^T)^-1 (X^TY + xy_N+1)

можно записать как

B_N+1 = P_N+1 (X^TY + xy_N+1),

где

Последнее выражение, полученное с помощью леммы об обращении матриц, позволяет получить откорректированную матрицу P_N.