Журнал вычислительной математики и математической физики, 2023, T. 63, № 9, стр. 1415-1427

1. ВВЕДЕНИЕ

Матричный оператор ${\text{vec}}\left( . \right)$, векторизующий квадратную матрицу, располагая ее столбцы друг под другом, широко используется для облегчения математических вычислений с матрицами [1, 2]. В частности, линейное матричное уравнение с матрицей неизвестных $P$ с помощью ${\text{vec}}\left( . \right)$ преобразуется в систему линейных алгебраических уравнений, к которой применим весь известный арсенал аналитических и численных методов поиска ее решения. В [1, 2] также вводятся два варианта данной операции для случая, когда такая матрица является симметричной $P = {{P}^{{\text{T}}}}$: ${\text{vech}}\left( . \right)$ и ${\text{vecd}}\left( . \right)$. Обе операции возвращают вектор-столбец, который содержит в себе только уникальные элементы матрицы $P$, т.е. ее ниже-треугольную часть. Для ${\text{vech}}\left( . \right)$ и ${\text{vecd}}\left( . \right)$ различается лишь порядок расположения таких элементов, при этом операции обладают сходными свойствами, ввиду чего в дальнейшем тексте работы без потери общности будет упоминаться только операция ${\text{vech}}\left( . \right)$. Применение таких операций позволяет существенно сократить размерность получаемой после векторизации системы уравнений, а также в принципе гарантировать возможность получения решения системы уравнений (см. лемму 7.1 в [1]). Однако существуют задачи, в которых в матричных уравнениях матрица неизвестных параметров обладает более сложной симметрией по сравнению с уже упомянутым случаем $P = {{P}^{{\text{T}}}}$. Например, если такая матрица является симметричной и состоит из блоков равной размерности, каждый из которых также является симметричной квадратной матрицей – $P \otimes P$, где $P = {{P}^{{\text{T}}}}$. В этом случае при применении операции векторизации размерность решаемой задачи можно было бы сократить еще сильнее, поскольку повторяющихся элементов в такой матрице больше, чем в обычной симметричной. Однако ни одна из имеющихся операций ${\text{vec}}\left( . \right)$, ${\text{vech}}\left( . \right)$ не позволяет этого сделать. Ниже рассмотрим один из классов таких задач, в которых в качестве матрицы неизвестных возникает $P \otimes P$ при $P = {{P}^{{\text{T}}}}$.

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

При применении многих подобных методов первичной задачей является приведение исходной системы к линейному регрессионному уравнению (LRE) относительно неизвестных параметров [6]:

где y ∈ ${{{\text{R}}}^{{N \times 1}}}$ – вектор измеримых выходных величин, ω ∈ ${{{\text{R}}}^{{N \times M}}}$ – измеримый регрессор, θ ∈ ${{{\text{R}}}^{{M \times 1}}}$ – вектор неизвестных параметров. Разумеется, желательно, чтобы размерность θ была минимально возможной для решения поставленной задачи. Для улучшения динамических свойств регрессора ω возможно применение различных схем предобработки, таких как DRE [7], MRE [8] для получения матричного квадратного регрессора Ω ∈ ${{{\text{R}}}^{{M \times M}}}$, а также DREM [6] для получения скалярного регрессора Δ ∈ R. Независимо от размерности регрессора, для идентификации неизвестных параметров θ обычно используются классические алгоритмы идентификации, такие как, gradient descent или least squares.

В целом методы теории адаптивного управления обладают рядом недостатков, одним из основных среди которых является отсутствие гарантий на качество получаемых в системе переходных процессов [9]. Например, асимптотическая или экспоненциальная сходимость не гарантирует отсутствие перерегулирования и/или колебательности. Поэтому становится актуальным применение методов теории адаптивного управления для задач, традиционно решаемых в рамках теории оптимального управления.

Теория оптимального управления позволяет гарантировать наилучшее возможное качество управления с точки зрения целевой функции, выбранной пользователем [10]. Для линейных систем классическим решением задачи оптимального управления является синтез регулятора, минимизирующего линейно-квадратичный (LQ) функционал от состояний и управлений на бесконечном времени [11, 12]. Параметры такого регулятора вычисляются на основе матрицы P, являющейся решением алгебраического уравнения Риккати (ARE).

Для синтеза LQ регулятора в адаптивной постановке одним из подходов является приведение ARE к виду (1). Для этого осуществляют [13] подстановку в него вместо неизвестных матриц объекта A и B их динамических оценок $\hat {A}{\text{,}}\;\hat {B}$ или некоторых измеримых сигналов. На основе идентификации неизвестных параметров полученной таким образом регрессии (1) и определяют искомую матрицу P. Однако задача идентификации матрицы $P$ осложняется тем обстоятельством, что после приведения ARE к (1) вектор неизвестных параметров θ, кроме искомой матрицы $P$, включает в себя также перепараметризованную часть, определяемую выражением ${\text{vec}}\left( {P \otimes P} \right)$. Как хорошо известно, решение задачи идентификации неизвестных параметров при наличии перепараметризации обладает рядом существенных недостатков (неудовлетворительное качество оценок, потеря возбуждения регрессором и пр.), поэтому актуальным является уменьшение числа перепараметризованных параметров или их полное исключение. Как уже упоминалось выше, единственным существующим на сегодняшний день подходом к уменьшению числа параметров в перепараметризованной части уравнения ARE, записанного в виде (1), является использование вместо ${\text{vec}}\left( {P \otimes P} \right)$ процедуры частичной векторизации ${\text{vech}}\left( {P \otimes P} \right)$. Однако в силу симметрии матрицы $P$ и симметрии произведения $P \otimes P$ перепараметризованная часть в данной задаче содержит повторы и в нижне-треугольной части, а поэтому использование процедуры ${\text{vech}}\left( {P \otimes P} \right)$ не является оптимальным и возвращает вектор с повторяющимися элементами. В литературе авторам не удалось обнаружить алгоритма операции векторизации $P \otimes P$, результатом которого является вектор, содержащий только уникальные элементы. Также не были обнаружены оценки на число уникальных элементов в произведении $P \otimes P$ для матрицы $P = {{P}^{{\text{T}}}}$ произвольно выбранного размера. Необходимо упомянуть, что в работах [13, 14] декларативно упоминается необходимость наличия такой операции, но конструктивный алгоритм ее выполнения не приводится, а полученная в [13] формула для оценки числа уникальных параметров при некоторых значениях размера матрицы $P = {{P}^{{\text{T}}}}$ выдает нецелочисленный ответ. Некоторые дополнительные пояснения по результатам работам [13, 14] приведены в разд. 2 настоящей статьи.

Таким образом, в развитие результатов работ [1, 2] для случая более сложной симметрии матрицы неизвестных параметров, целью настоящей работы является разработка алгоритма векторизации произведения ($P \otimes P$) с целью получения вектора, включающего в себя только уникальные элементы такого произведения. Дополнительно будет получена зависимость числа уникальных элементов рассматриваемого произведения от размерности исходной матрицы P. Основным результатом работы является операция векторизации произведения $P \otimes P$ с исключением повторений, которая может быть использована как в задаче адаптивного LQ синтеза, так и в других задачах, в которых приходится оперировать выражением вида $P \otimes P{\text{,}}$ $P = {{P}^{{\text{T}}}}$.

2. ПОСТАНОВКА ЗАДАЧИ

Следуя работам [13, 14], рассмотрим классическое ARE, к решению которого сводится задача конструирования оптимального LQ регулятора:

Здесь $A \in {{R}^{{n \times n}}}$ – матрица состояний, $B \in {{R}^{{n \times m}}}$ – матрица коэффициентов усиления, $R = {{R}^{{\text{T}}}} > 0{\text{,}}$ $Q = {{Q}^{{\text{T}}}} > 0$ – матрицы соответствующей размерности, являющиеся параметрами целевого функционала LQ оптимизации, $P = {{P}^{{\text{T}}}} > 0$ – матрица соответствующей размерности, являющаяся решением ARE.

Применив ${\text{vec}}\left( . \right)$ к выражению (3), получим

Использовав второе свойство операции векторизации, имеем

Так как $P = {{P}^{{\text{T}}}}{\text{,}}$ $Q = {{Q}^{{\text{T}}}}{\text{,}}$ $B{{R}^{{ - 1}}}{{B}^{{\text{T}}}} = {{\left( {B{{R}^{{ - 1}}}{{B}^{{\text{T}}}}} \right)}^{{\text{T}}}}$, то ${\text{vec}}\left( {P \otimes P} \right){\text{,}}$ ${\text{ve}}{{{\text{c}}}^{{\text{T}}}}\left( {B{{R}^{{ - 1}}}{{B}^{{\text{T}}}}} \right){\text{,}}$ $vec\left( Q \right)$ содержат повторы, а поэтому для исключения неуникальных элементов рационально перейти с учетом свойств, описанных в определении 2, от $vec\left( . \right)$ к ${\text{vech}}\left( . \right)$.

Утверждение 1. Уравнение Риккати (4) относительно матрицы $P = {{P}^{{\text{T}}}} > 0$ может быть записано в виде линейного уравнения с перепараметризацией:

в котором матрицы $L \in {{R}^{{k \times {{n}^{2}}}}}{\text{,}}$ $D \in {{R}^{{{{n}^{2}} \times k}}}{\text{,}}$ ${{L}_{1}} \in {{R}^{{r \times {{n}^{4}}}}}{\text{,}}$ ${{D}_{1}} \in {{R}^{{{{n}^{4}} \times r}}}$ задают операции исключения и возвращения над $Q$ и $P \otimes P$ и по определению 2 удовлетворяют условиям $L = {{\left( {{{D}^{{\text{T}}}}D} \right)}^{{ - 1}}}{{D}^{{\text{T}}}}{\text{,}}$ ${{L}_{1}} = {{\left( {D_{1}^{{\text{T}}}{{D}_{1}}} \right)}^{{ - 1}}}D_{1}^{{\text{T}}}{\text{,}}$ $r = 0.5{{n}^{2}}\left( {{{n}^{2}} + 1} \right),$ $k = 0.5n\left( {n + 1} \right)$.

Доказательство утверждения непосредственно следует из применения операции частичной векторизации к ${\text{vec}}\left( {P \otimes P} \right)$, ${\text{vec}}\left( P \right)$, ${\text{ve}}{{{\text{c}}}^{{\text{T}}}}\left( {B{{R}^{{ - 1}}}{{B}^{{\text{T}}}}} \right)$, ${\text{vec}}\left( Q \right)$ в уравнении (4) с учетом свойств матриц возвращения и исключения. Рассмотрим иллюстративный пример.

Пример 2. Пусть n = 2, m = 1, а матрицы ARE определены следующим образом:

Тогда ARE в форме (4) примет вид:

где Очевидно, что второе и третье уравнения (6) совпадают, как и то, что ${{P}_{{{\text{kron}}}}} = {\text{vec}}\left( {P \otimes P} \right)$ содержит множество повторов и всего шесть уникальных элементов.

Теперь рассмотрим запись уравнения (3) с использованием частичной векторизации (5):

В данном случае

а вектор ${{P}_{{{\text{kron}}}}} = {{\left[ {\begin{array}{*{20}{c}} {p_{1}^{2}}&{{{p}_{1}}{{p}_{2}}}&{{{p}_{1}}{{p}_{2}}}&{p_{2}^{2}}&{{{p}_{1}}{{p}_{4}}}&{p_{2}^{2}}&{{{p}_{2}}{{p}_{4}}}&{{{p}_{1}}{{p}_{4}}}&{{{p}_{2}}{{p}_{4}}}&{p_{4}^{2}} \end{array}} \right]}^{{\text{T}}}}.$

Из сравнения (7) и (6) следует, что число уравнений (по определению k из утверждения 1) сократилось до трех (повторяющееся уравнение было отброшено), а длина вектора ${{P}_{{{\text{kron}}}}}$ (по определению r из утверждения 1) – до 10. Тем не менее в ${{P}_{{{\text{kron}}}}}$ все также присутствует четыре повтора, число которых, очевидно, будет существенно возрастать с ростом размерности n.

Таким образом, пример 2 демонстрирует, что операция $\operatorname{vech} \left( . \right)$ не позволяет сформировать из ${{P}_{{{\text{kron}}}}}$ вектор уникальных параметров минимальной длины, поскольку неуникальные элементы в $P \otimes P$ находятся не только в верхнетреугольной части, но и в нижнетреугольной.

В работах [13, 14] выражения, сходные с (5), рассматриваются как линейная регрессия, вектор параметров которой идентифицируется on-line. В такой задаче вектор параметров должен иметь минимально возможную перепараметризацию – фактически, состоять только из уникальных составляющих. Однако работы [13, 14] ограничены наличием только процедуры ${\text{vech}}\left( . \right)$, что приводит к проблемам, продемонстрированным в примере 2. Поэтому в [13] сделано допущение о том, что желаемая процедура векторизации существует и уже применена к (5), не поясняя, как именно ее выполнить. Кроме того, они дают оценку на количество уникальных элементов в $P \otimes P$, $P = {{P}^{{\text{T}}}}$: $\tfrac{{{{n}^{4}} + {{n}^{3}} + 7{{n}^{2}} + 6n}}{8}$, которая, как будет продемонстрированно ниже, является неверной.

Поэтому актуальной является задача разработки операции, выполняющей преобразование произведения $P \otimes P$ в вектор, наполненный только уникальными элементами.

Цель. В настоящей работе требуется, во-первых, построить алгоритм формирования матрицы исключения $\bar {L} \in {{R}^{{r \times {{n}^{4}}}}}$, выполняющей преобразование ${\text{vec}}\left( {P \otimes P} \right)$ в вектор, состоящий из уникальных элементов (оптимизируя тем самым размерность $r$), а во-вторых, получить аналитическую формулу зависимости числа уникальных элементов в произведении $P \otimes P$ от размерности матрицы $P$.

3. ОСНОВНОЙ РЕЗУЛЬТАТ

Чтобы сформулировать алгоритм исключения из $P \otimes P$ повторяющихся элементов, рассмотрим подробней представленный на фиг. 1 результат умножения по Кронекеру $P$ саму на себя в случае n = 2.

Таким образом, поскольку результат произведения $P \otimes P$ при $P = {{P}^{{\text{T}}}}$ является симметричной матрицей, то элементы верхнетреугольной формы (1,2), (1,3), (1,4), (2,3), (2,4), (3,4) считаются неуникальными. Кроме того, если результат произведения $P \otimes P$ разбить на блоки размером $n \times n$, то в каждом из полных блоков нижнетреугольной части в силу $P = {{P}^{{\text{T}}}}$ наблюдается симметрия относительно главной диагонали блока, а тогда элемент (3,2) считается неуникальным. Дополнительно, чем правее и ниже расположен блок в нижнетреугольной части, тем больше в нем повторов элементов с блоками, стоящими левее и/или выше – межблочная симметрия (элементы (4,3) и (3,1) считаются неуникальными). На фиг. 1 в красные рамки взяты элементы, находящиеся ниже главной диагонали $P \otimes P$, но при этом являющиеся неуникальными при обходе нижнетреугольной части по блокам сверху вниз и слева направо.

На основе описанной и отмеченной на фиг. 1 естественной симметрии в матрице $P \otimes P$ формализуем в виде следующего утверждения алгоритм формирования матрицы исключения $\bar {L} \in {{R}^{{r \times {{n}^{4}}}}}$ и возвращения $\bar {D} \in {{R}^{{{{n}^{4}} \times r}}}$.

Утверждение 2. Матрица исключения $\bar {L} \in {{R}^{{r \times {{n}^{4}}}}}$, обеспечивающая $\bar {L}{\text{vec}}\left( {P \otimes P} \right) \in {{R}^{r}}{\text{,}}$ $r = \left( {\tfrac{{{{n}^{4}}}}{{12}} + \tfrac{{{{n}^{3}}}}{3} + \tfrac{{5{{n}^{2}}}}{{12}} + \tfrac{n}{6}} \right)$ может быть сформирована следующим конструктивным алгоритмом:

где $p$ – номер вектор-строки с числом строк $n$ в блочной матрице ${{E}_{{{{\wp }_{l}}}}}$, q – номер вектор-столбца с числом столбцов $n$ в блочной матрице ${{E}_{{{{\wp }_{l}}}}}$, $i$ – номер строки в блоке $\left( {p,q} \right)$,$j$ – номер столбца в блоке $\left( {p,q} \right)$, ${{u}_{l}} \in {{R}^{r}}$ – вектор, заполненный нулями и содержащий единицу на l-й позиции, ${{E}_{\wp }}$– множество мощностью $r$ всех возможных матриц с учетом $\left( {i,j,p,q} \right) \in \wp $, заполненных нулями и содержащих единицы в позициях $\left( {i,j,p,q} \right){\text{,}}\;\left( {j,i,p,q} \right)$, $\left( {i,j,q,p} \right){\text{,}}\;\left( {j,i,q,p} \right)$, $\left( {p,q,i,j} \right){\text{,}}\;\left( {q,p,i,j} \right)$, $\left( {p,q,j,i} \right){\text{,}}\;\left( {q,p,j,i} \right)$.

Доказательство. Поскольку алгоритм (8) является конструктивным, то справедливость утверждения 2 следует из отсутствия логических противоречий в рассуждениях, использованных для его формирования. Далее кратко остановимся на основных рассуждениях, на основе которых и был предложен алгоритм (8).

В первую очередь, матрицу $P \otimes P \in {{R}^{{{{n}^{2}} \times {{n}^{2}}}}}$ разобьем на ${{n}^{2}}$ блоков размерностью $n \times n$. Пусть каждый блок имеет уникальный составной индекс $\left( {p{\text{, }}q} \right)$, где $p$ – номер вектор-строки с числом строк $n$ в блочной матрице $P \otimes P$, $q$ – номер вектор-столбца с числом столбцов $n$ в блочной матрице $P \otimes P$. Для обращения к каждому элементу блока также введем уникальный составной индекс $\left( {i{\text{, }}j} \right)$, где $i$ – номер строки в блоке $\left( {p{\text{, }}q} \right)$, $j$ – номер столбца в блоке $\left( {p{\text{, }}q} \right)$. Таким образом, каждый элемент матрицы $P \otimes P$ может быть однозначно определен уникальным составным индексом $\left( {i{\text{, }}j{\text{, }}p{\text{, }}q} \right)$.

Теперь, пользуясь симметрией, отмеченной на фиг. 1, определим позиции уникальных элементов в матрице $P \otimes P$, которые находятся на или ниже главной диагонали.

Ввиду общей симметричности матрицы $P \otimes P$, такие элементы должны находиться в блоках, находящихся ниже главной диагонали $P \otimes P$, или включающих ее. Таким образом, при учете p = 1, …, n, q = 1, …, n, необходимым условием уникальности является $p \geqslant q$.

Учитывая, что каждый блок матрицы $P \otimes P$ сам представляет собой симметричную матрицу, необходимо считать уникальными только элементы блока, которые находятся на или ниже его главной диагонали. Таким образом, при учете i = 1, …, n, j = 1, …, n, необходимым условием является $i \geqslant j$.

Принимая во внимание блочную структуру матрицы $P \otimes P$ и тот факт, что $P = {{P}^{{\text{T}}}}$, для элемента с индексом $\left( {i{\text{, }}j{\text{, }}p{\text{, }}q} \right)$ обязательно будет симметричным элемент $\left( {p,q,i,j} \right)$. Ведь элемент $\left( {i{\text{, }}j{\text{, }}p{\text{, }}q} \right)$ – это фактически, следуя определению операции произведения Кронекера, произведение $\left( {i{\text{, }}j} \right)$ элемента матрицы P на ее элемент $\left( {p{\text{, }}q} \right)$, а $\left( {p,q,i,j} \right)$ – это произведение $\left( {p{\text{, }}q} \right)$ элемента матрицы P на элемент $\left( {i{\text{, }}j} \right)$. Поэтому необходимым условием оптимальности является $p \leqslant i,$ $q \leqslant j$.

Таким образом, объединив сформулированные условия, может быть сформировано множество индексов уникальных элементов матрицы $P \otimes P$: $\wp {\text{ :}} = \left\{ {\left( {i,j,p,q} \right)} \right.:{\text{ }}i \geqslant j > 0,{\text{ }}p \geqslant q > 0,$ $\left. {p \leqslant i \leqslant n,{\text{ }}q \leqslant j \leqslant n} \right\}.$ Примем пока, что таких уникальных элементов r, а обоснование формулы для подсчета их числа приведем ниже.

Для формирования матрицы $\bar {D}$ путем, аналогичным (2) [1], необходимо для каждого l-го (l = 1, …, r) уникального элемента матрицы $P \otimes P$ (l-го члена множества $\wp $, который обозначим ${{\wp }_{l}}$) сформировать матрицу ${{E}_{{{{\wp }_{l}}}}} \in {{R}^{{{{n}^{2}} \times {{n}^{2}}}}}$, которая содержит единицы на всех позициях, соответствующих позициям, в которых в исходной матрице $P \otimes P$ встречается элемент ${{\wp }_{l}} = \left( {{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}}} \right)$.

Таких позиций может быть максимум восемь: сама позиция элемента $\left( {{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}}} \right)$; симметричный элемент относительно главной диагонали блока $\left( {p{\text{, }}q} \right)$ – $\left( {{{j}^{{\left( l \right)}}},\;{{i}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}}} \right)$; как уже упоминалось выше, элемент $\left( {{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}}} \right)$; также элемент, симметричный $\left( {{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}}} \right)$ относительно главной диагонали блока (i, j) – $\left( {{{q}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}}} \right)$, а также четыре элемента, симметричных уже упомянутым четырем относительно главной диагонали матрицы $P \otimes P$ – $\left( {{{i}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}}} \right){\text{, }}\left( {{{j}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}}} \right)$, $\left( {{{p}^{{\left( l \right)}}},{{q}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}}} \right){\text{, }}\left( {{{q}^{{\left( l \right)}}},{{p}^{{\left( l \right)}}},{{j}^{{\left( l \right)}}},{{i}^{{\left( l \right)}}}} \right)$. Для некоторых элементов в матрице $P \otimes P$ часть из этих восьми индексов могут указывать на один и тот же элемент (например, для элемента (1, 1, 1, 1)). В этом нет противоречия, поскольку это будет лишь означать, что установить единицу в некоторый элемент матрицы ${{E}_{{{{\wp }_{l}}}}}$ необходимо несколько раз.

Транспонированная векторизация матрицы ${{E}_{{{{\wp }_{l}}}}}$ представляет собой вектор размерностью $1 \times {{n}^{4}}$, который должен занять l-й столбец в матрице $\bar {D}$, что сходно с (2). Заявленная в утверждении процедура формирования вектора ${{u}_{l}} \in {{R}^{r}}$ и формула

позволяют сформировать именно такую матрицу.

От полученной матрицы дублирования с помощью псевдообращения Мура-Пенроуза возможно получить матрицу исключения $\bar {L}$. При этом существование $\bar {L}$ гарантируется взаимной ортогональностью ${\text{vec}}\left( {{{E}_{{{{\wp }_{l}}}}}} \right)$ $\forall l = \overline {1,r} $, в свою очередь обеспечиваемой $\left( {i,j,p,q} \right) \in \wp $, т.е. тем фактом, что элементы сформированного множества $\wp $ уникальны.

Для подсчета числа уникальных элементов в матрице $P \otimes P$ воспользуемся следующим подходом. В первую очередь, число уникальных элементов матрицы $P \otimes P$ не больше числа уникальных элементов в ее нижне-треугольной части, определяемого по формуле 0.5(n⁴ + n²). Из этой величины необходимо вычесть число элементов, расположенных выше главной диагонали своего блока для всех блоков, расположенных полностью ниже главной диагонали матрицы $P \otimes P$. Таких блоков ровно n² – 0.5(n² + n), поскольку всего их n², а расположенных на или выше главной диагонали – 0.5(n² + n). В каждом из таких блоков по n² – 0.5(n² + n) элементов выше главной диагонали. Поэтому из 0.5(n⁴ + n²) необходимо вычесть (n² – 0.5(n² + n))².

Наконец, необходимо подсчитать число неуникальных элементов в нижне-треугольных частях всех блоков матрицы $P \otimes P$, расположенных на или ниже главной диагонали $P \otimes P$. Рассмотрение матрицы $P \otimes P$ необходимо вести по блочным вектор-столбцам q = 1, …, n. В первом таком столбце по мере увеличения номера вектор-строки p = 1, …, n, количество неуникальных элементов будет возрастать в виде ряда, члены которого представляют собой сумму элементов арифметической прогрессии 0, 1, 2, 3, 4, …, т.е. 0, 1, 3, 6, 10, … ввиду упомянутого выше совпадения элементов (i, j, p, q) и (p, q, i,  j) для всех допустимых i,  j, p, q. Сумма такой последовательности вычисляется как (n – 1)n(n + 1)/6. По мере роста q число неуникальных элементов в каждом последующем столбце будет возрастать, поскольку при формировании множества $\wp $ были выставлены требования $p \leqslant i{\text{,}}$ $q \leqslant j$. То есть при возрастании q на единицу в каждом из (n – q) блоков соответствующего вектор-столбца рассматривается на один столбец элементов меньше, что, при учете, что $p \geqslant q$ означает отбрасывание еще (nq – (q – 1)q/2) неуникальных элементов. При этом для оставшихся в рассмотрении элементов каждого блока действует все то же правило суммы элементов арифметической прогрессии. Таким образом, общее количество неуникальных элементов в нижне-треугольных частях всех блоков матрицы $P \otimes P$, расположенных на или ниже главной диагонали $P \otimes P$, рассчитывается по формуле:

А общее число уникальных элементов в матрице $P \otimes P$ вычисляется по формуле

что завершает доказательство утверждения.

На основе конструктивного алгоритма формирования матрицы исключения $\bar {L}$ введем определение операции исключения повторений из ${\text{vec}}\left( {P \otimes P} \right)$.

Определение 5. Операция исключения повторений $\operatorname{vecu} \left( . \right)$ из вектора ${\text{vec}}\left( {M \otimes M} \right) \in {{R}^{{{{n}^{4}}}}}{\text{,}}$ $M = {{M}^{{\text{T}}}} \in {{R}^{{n \times n}}}$, задается следующим образом:

где $\bar {L}$ и $r$ определены в утверждении 2, а $\bar {L}$ обладает следующими свойствами:

1) $\bar {L}$ имеет полный строчный ранг r;

2) ${{\bar {L}}^{\dag }} = {{\bar {L}}^{{\text{T}}}}$, где ${{\bar {L}}^{\dag }}$ – псевдоинверсия Мура–Пенроуза матрицы $\bar {L}$.

Существует матрица $\bar {D}$ такая, что:

3) $\bar {L}{{\bar {L}}^{{\text{T}}}} = \bar {L}\bar {D} = {{I}_{r}}$;

4) ${\text{vec}}\left( M \right) = \bar {D}\bar {L}\operatorname{vec} \left( M \right) = \bar {D}\operatorname{vecu} \left( M \right).$

Справедливость определения 5 непосредственно следует из доказательства утверждения 2.

Замечание 1. Результатом выполнения операции $\operatorname{vecu} \left( . \right)$ является вектор уникальных элементов $P \otimes P$, а также матрицы исключения и дублирования, поэтому введение отдельных дополнительных массивов для длительного хранения индексов уникальных элементов в $P \otimes P$ (которые упоминались в доказательстве) не требуется.

Ниже приведен код на языке Matlab Script Language, реализующий в соответствии с конструктивным алгоритмом (8) формирование матриц дополнения $\bar {D}$ и исключения $\bar {L}$.

%порядок задачи

n=2;

%число уникальных элементов в матрице P kron P

sum=n*(n-1)*(n+1)/6;

for q=1:1:(n-1)

sum=sum+(n-q)*(n*q-q*(q-1)/2)+(n-q)*(n-q+1)*(n-q-1)/6;

end

result_unique = (n^4+n^2)/2 - sum - (n*n-n*(n+1)/2)^2;

%код для формирования матрицы возвращения

D=zeros(n^4,result_unique);

l=1;

for q=1:1:n

for p=q:1:n

for j=q:1:n

if j>p

k=j;

else

k=p;

end

for i=k:1:n

E=zeros(n^2,n^2);

E((p-1)*n+i,(q-1)*n+j)=1;

E((p-1)*n+j,(q-1)*n+i)=1;

E((q-1)*n+i,(p-1)*n+j)=1;

E((q-1)*n+j,(p-1)*n+i)=1;

E((i-1)*n+p,(j-1)*n+q)=1;

E((i-1)*n+q,(j-1)*n+p)=1;

E((j-1)*n+p,(i-1)*n+q)=1;

E((j-1)*n+q,(i-1)*n+p)=1;

E=reshape(E,n^2*n^2,1);

D(:,l)=E;

l=l+1;

end

end

end

end

%код для формирования матрицы исключения

L=pinv(D);

Оценим сложность предложенного конструктивного алгоритма. Необходимо отметить, что для подсчета числа уникальных элементов result_unique она является линейной O(n). Для построения матрицы D циклы for организованы таким образом, что будет осуществлен обход только уникальных r (см. (9)) штук элементов $P \otimes P$. При этом для каждого такого элемента будет осуществляться векторизация матрицы E со сложностью O(n⁴). То есть сложность этой части алгоритма составляет O(n⁴ ⋅ r). Сложность вычисления матрицы L совпадает со сложностью получения псевдообратной матрицы Мура–Пенроуза, которая для данного случая будет равна O((n⁴)² ⋅ r) = O(n⁸ ⋅ r). Наибольшую сложность имеет последняя часть алгоритма. Она же и определяет его сложность в целом.

4. ЧИСЛЕННЫЙ ПРИМЕР

Рассмотрим два иллюстративных примера (n = 2, n = 3), поясняющих предложенный конструктивный алгоритм формирования матрицы исключения.

Пример 3. Продолжим рассмотрение ARE из примера 2. На фиг. 2 приведено расположение уникальных элементов в произведении $P \otimes P$ в соответствии с определением множества $\wp $.

Таким образом, p = 1, 2; q = 1, 2. Согласно сформированным правилам, множество $\wp = \left\{ {\left( {1,1,1,1} \right),\left( {2,1,1,1} \right),\left( {2,2,1,1} \right),\left( {2,1,2,1} \right),\left( {2,2,2,1} \right),\left( {2,2,2,2} \right)} \right\}$, а согласно формуле расчета (9), количество уникальных элементов $r = \left( {\tfrac{{{{2}^{4}}}}{{12}} + \tfrac{{{{2}^{3}}}}{3} + \tfrac{{5 \cdot {{2}^{2}}}}{{12}} + \tfrac{2}{6}} \right) = 6$.

На основе определения множества $\wp $ составим шесть матриц ${{E}_{{{{\wp }_{l}}}}}$:

По формуле (8) сформируем матрицу:

Все ее столбцы уникальны, поэтому матрица $\bar {L}$ существует и имеет вид

Тогда, пользуясь четвертым свойством операции из определения 5, по аналогии с (5) получим:

Откуда, подставив значения матриц, в рассматриваемом случае получим:

где $V = \left[ {\begin{array}{*{20}{c}} {b_{1}^{2}}&{2{{b}_{1}}{{b}_{2}}}&0&{b_{2}^{2}}&0&0 \\ 0&{b_{1}^{2}}&{{{b}_{1}}{{b}_{2}}}&{{{b}_{1}}{{b}_{2}}}&{b_{2}^{2}}&0 \\ 0&0&0&{b_{1}^{2}}&{2{{b}_{1}}{{b}_{2}}}&{b_{2}^{2}} \end{array}} \right],$ а вектор ${{P}_{{{\text{kron}}}}} = {{\left[ {\begin{array}{*{20}{c}} {p_{1}^{2}}&{{{p}_{1}}{{p}_{2}}}&{{{p}_{1}}{{p}_{4}}}&{p_{2}^{2}}&{{{p}_{2}}{{p}_{4}}}&{p_{4}^{2}} \end{array}} \right]}^{{\text{T}}}}$.

Таким образом, получена минимально возможная параметризация исходной задачи. При использовании операции ${\text{vec}}\left( {P \otimes P} \right)$ вектор параметров P_kron содержал 16 элементов, после применения операции ${\text{vech}}\left( {P \otimes P} \right)$ – 10 (см. пример 2), а после применения разработанной в работе операции ${\text{vecu}}\left( {P \otimes P} \right)$ – 6. Данный пример демонстрирует практическую ценность предлагаемой операции. В задаче идентификации неизвестных параметров регрессии типа (5) общее число неизвестных параметров было сокращено с 10 до 6, что позволяет избежать излишней перепараметризации и максимально сократить размерность решаемой задачи.

Также на примере (10) наглядно видно, что результат выполнения операции $\operatorname{vecu} \left( {P \otimes P} \right)$ единожды используется для получения параметризации типа (10) для дальнейшего решения ARE. Длительное хранение матриц исключения и дублирования не требуется. Хранится лишь вектор уникальных элементов $P \otimes P$, который по количеству элементов строго меньше общего числа элементов $P \otimes P$.

Пример 4. Пусть n = 3. В данном примере ограничимся изображением с отметками уникальных элементов в матрице $P \otimes P$ и подсчетом их количества по формуле (9). На фиг. 3 приведено расположение уникальных элементов в матрице $P \otimes P$ в соответствии с определением множества $\wp $.

Количество уникальных элементов $r = \left( {\tfrac{{{{3}^{4}}}}{{12}} + \tfrac{{{{3}^{3}}}}{3} + \tfrac{{5 \cdot {{3}^{2}}}}{{12}} + \tfrac{3}{6}} \right) = 20$, что подтверждается фиг. 3. При использовании операции ${\text{vec}}\left( {P \otimes P} \right)$ вектор параметров P_kron содержал 81 элемент, после применения операции ${\text{vech}}\left( {P \otimes P} \right)$ содержал 45 элементов, а после применения разработанной в работе операции ${\text{vecu}}\left( {P \otimes P} \right)$содержал 20 элементов.

Несложно убедиться, что для n = 4 по сравнению с операцией ${\text{vec}}\left( {P \otimes P} \right)$ при использовании $\operatorname{vecu} \left( {P \otimes P} \right)$ число параметров сокращается с 256 до 50, для n = 5 – с 625 до 105 и т.д.

5. ЗАКЛЮЧЕНИЕ

В данном исследовании разработан конструктивный алгоритм выполнения операции произведения Кронекера симметричной матрицы на саму себя $P \otimes P$ с последующей векторизацией результата таким образом, чтобы получаемый вектор включал в себя только уникальные элементы – пары произведений элементов исходной матрицы. Данный алгоритм основан на определении операции Кронекера и симметричности исходной матрицы. Предложен алгоритм поиска индексов уникальных элементов в $P \otimes P$, а также формулы для получения матриц исключения и дублирования, которые базируются на классической работе [1] по частичной векторизации матриц. Кроме того, была получена зависимость числа таких уникальных элементов от размерности исходной матрицы.

Такой подход позволяет получать минимальную по числу параметров параметризацию уравнения Риккати в виде (1), привнося тем самым необходимое обоснование в работы [13, 14], где подобная процедура была заявлена декларативно.

Предложенная операция может применяться не только в задачах, связанных с решением ARE, но и в любых других, где возникает необходимость векторизации выражения вида $P \otimes P$, при условии P= P^T.

Актуальной представляется задача расширения полученного результата на случай векторизации произведений вида $A \otimes B$, где $A \in {{R}^{{n \times n}}}$ и $B \in {{R}^{{n \times n}}}$ – произвольные квадратные матрицы.

В дальнейшем авторы планируют воспользоваться предложенной операцией $\operatorname{vecu} \left( . \right)$ и полученными результатами по минимальной параметризации уравнения Риккати для построения адаптивного LQ регулятора.