Журнал вычислительной математики и математической физики, 2019, T. 59, № 11, стр. 1823-1835

ВВЕДЕНИЕ

Решение многих прикладных задач математического моделирования, теории автоматического управления и теории систем [1]–[6], а также, например, регрессионного анализа часто приводит к переопределенным, недоопределенным или вырожденным СЛАУ. Кроме того, на практике часто возникают обратные задачи с плохо обусловленными матрицами, для решения которых, как правило, невозможно непосредственно вычислить обратную матрицу, применить метод Гаусса или другие классические методы линейной алгебры.

В работах [1]–[6] была разработана и применена теория вложения систем, основывающаяся на анализе алгебраических особенностей оператора системы и позволяющая синтезировать полный класс решений задачи управления данной системой. Для анализа алгебраических особенностей систем в работах [7] и [8] была разработана технология канонизации, которая позволяет строить полный класс решений любых СЛАУ: переопределенных, недоопределенных, вырожденных или плохо обусловленных. Технология канонизации успешно применяется для решения задач теории автоматического управления, например, в работах [9]–[11].

Похожие идеи излагаются в работах Г. Стрэнга [12] и [13], где вводятся понятия правого и левого нуль-пространств, а также ортогональных им пространств строк и столбцов матрицы, и устанавливается связь между LU-разложением матрицы и базисами данных подпространств. Отдельно отметим, что канонизацию матрицы можно рассматривать как обобщение сингулярного разложения, на что было указано в работе [6].

Для получения канонизации матрицы в работе [6] предлагается планшетный метод, основанный на методе Гаусса–Жордана и пригодный как для ручного счета, так и для программной реализации на ЭВМ. Существует программная реализация планшетного метода канонизации [14], главным достоинством которой является возможность работы с СЛАУ, имеющими дробно-рациональные коэффициенты, что позволяет решать задачи синтеза систем управления, сформулированные в форме матричных передаточных функций. В качестве основных недостатков планшетного метода можно отметить невозможность получения канонизаторов и делителей нуля с заданными свойствами, например, являющихся унитарными, а также теоретическую неустойчивость метода Гаусса.

Из работ по вычислительной линейной алгебре (например, [15] и [16]) известно, что наиболее эффективными являются алгоритмы решения СЛАУ, основанные на матричных разложениях. В частности, наиболее быстрыми методами для квадратных матриц являются те, что основаны на LU-разложении, для матриц произвольного размера – на QR-разложении, максимальная скорость может быть достигнута для эрмитовых положительно-определенных матриц при помощи разложения Холецкого, а для решения вырожденных и плохо обусловленных систем наиболее подходящими являются QR-, LQ- и сингулярное разложения. Следовательно, актуальной задачей является разработка эффективного алгоритма канонизации, основанного на данных разложениях и при этом сохраняющего возможность работы с СЛАУ, имеющими дробно-рациональные коэффициенты.

Проблема вычисления и оценки числа обусловленности задачи поднимается во многих работах, например, в [17] и [18]. Прямое вычисление числа обусловленности матрицы осложнено необходимостью ее обращения, в то время как эффективные методы решения СЛАУ не выполняют его явно. Поэтому практический интерес представляют методы вычисления и оценки числа обусловленности, не требующие обращения матрицы системы.

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

В разд. 1 настоящей работы приводятся общие сведения о технологии канонизации матриц (см., например, [6]–[8]). В разд. 2 осуществляется сведение задачи канонизации матрицы к задаче вычисления одного из матричных разложений, разрабатывается алгоритм канонизации, основанный на данных разложениях, и предлагается метод оценки числа обусловленности задачи канонизации, основанный на вычислении норм матриц, получаемых в результате канонизации, который не требует обращения исходной матрицы. В разд. 3 приводятся результаты вычисления канонизации и числа обусловленности тестовых матриц, выполненных в среде MATLAB с помощью предлагаемого алгоритма. В разд. 4 обсуждаются области применения предлагаемого алгоритма.

1. ОБЩИЕ СВЕДЕНИЯ О КАНОНИЗАЦИИ МАТРИЦ

Пусть рассматривается левосторонняя СЛАУ общего вида

Здесь и далее A_m × n – числовая или дробно-рациональная (полиномиальная) матрица коэффициентов размера m × n, B_m × p – числовая или дробно-рациональная (полиномиальная) матрица правых частей размера m × p, X_n × p – матрица с неизвестными числовыми или дробно-рациональными (полиномиальными) элементами размера n × p, m – число строк матриц A_m × n и B_m × p, n – число столбцов матрицы A_m × n и строк матрицы X_n × p, p – число столбцов матриц X_n × p и B_m × p.

В случае, если матрица A_m × n является квадратной и невырожденной, то есть m = n и rank(A_m × m) = m, то формально решением СЛАУ вида (1.1) является следующее выражение:

где (A_m × m)^–1 – обратная матрица. По теореме Кронекера–Капелли решение данной системы существует, если выполняется условие

Однако, если матрица ${{A}_{{m \times n}}}$ не является квадратной или невырожденной, то есть $m \ne n$ или rank(A_m × n) < min(m, n), то выражение (A_m × n)^–1 не имеет смысла. В этом случае, в зависимости от выполнения условий m < n или m > n, говорят о том, что система (1.1) является недоопределенной или переопределенной, соответственно, записать аналитическое решение в явном виде (1.2) оказывается невозможным.

В работах [6]–[8] предлагается обобщение подхода, описываемого выражениями (1.2) и (1.3). Данный метод авторы называют методом канонизации. Он заключается в нахождении в общем случае пятерки матриц

Здесь и далее матрица $\bar {A}_{{(m - r) \times m}}^{L}$ называется левым, а $\bar {A}_{{n \times (n - r)}}^{R}$ – правым делителями нуля максимального ранга,$\tilde {A}_{{r \times m}}^{L}$ – левым канонизатором (делителем единицы), $\tilde {A}_{{n \times r}}^{R}$ – правым канонизатором (делителем единицы) и ${{\tilde {A}}_{{n \times m}}}$ – сводным канонизатором, r – ранг матрицы ${{A}_{{m \times n}}}$.

Известны соотношения (строго обоснованные в работах [6] и [8]), связывающие вышеописанные матрицы

Здесь ${{I}_{{r \times r}}}$ – единичная матрица размера $r \times r$, а ${{O}_{{(m - r) \times n}}}$ и ${{O}_{{m \times (n - r)}}}$ – нулевые матрицы соответствующих размеров. Учитывая основную теорему линейной алгебры (см., например, [12] и [13]), из соотношений (1.5)–(1.8) следует, что левый делитель нуля $\bar {A}_{{(m - r) \times m}}^{L}$ составлен из строк, являющихся базисом левого нуль-пространства матрицы ${{A}_{{m \times n}}}$, правый делитель нуля $\bar {A}_{{n \times (n - r)}}^{R}$ – из столбцов, являющихся базисом (правого) нуль-пространства матрицы ${{A}_{{m \times n}}}$, строки левого канонизатора $\tilde {A}_{{r \times m}}^{L}$ непосредственно связаны с базисом пространства строк матрицы ${{A}_{{m \times n}}}$, а столбцы правого канонизатора $\tilde {A}_{{n \times r}}^{R}$ – с базисом пространства столбцов матрицы ${{A}_{{m \times n}}}$ [6].

Используя соотношения (1.5) и (1.8), можно показать, что сводный канонизатор ${{\tilde {A}}_{{m \times n}}}$ является обобщением понятия обратной матрицы ${{({{A}_{{m \times n}}})}^{{ - 1}}}$. Кроме того, в работе [8] доказываются теоремы, утверждающие, что на основе условия (1.3) и канонизации (1.4) можно сформулировать условие разрешимости СЛАУ (1.1)

и в совокупности с соотношениями (1.5) и (1.8) записать полный класс решений СЛАУ (1.1) в виде Здесь${{\eta }_{{(n - r) \times p}}}$ – матрица размера $(n - r) \times p$ с произвольными элементами ${{n}_{{i,j}}}$, $i = \overline {1,(n - r)} $, $j = \overline {1,p} $. Обозначением ${{\{ {{X}_{{n \times p}}}\} }_{\eta }}$ указывают на множественность возможных решений ${{X}_{{n \times p}}}$, зависящих от матричного параметра ${{\eta }_{{(n - r) \times p}}}$.

Таким образом, условие (1.9) и выражение (1.10) дают полное решение СЛАУ (1.1) в аналитическом виде. Следует отметить, что аналогичные результаты можно получить для правосторонних и двусторонних СЛАУ (см., например, [6] и [8]).

2. СВЕДЕНИЕ ЗАДАЧИ КАНОНИЗАЦИИ К ЗАДАЧЕ МАТРИЧНОГО РАЗЛОЖЕНИЯ

АЛГОРИТМ

Шаг 1. Определить размер матрицы ${{A}_{{m \times n}}}$.

Шаг 2. Если $m > n$, то перейти к шагу 3, если $m < n$, то перейти к шагу 5, иначе перейти к шагу 7.

Шаг 3. Вычислить матрицы ${{Q}_{{m \times m}}}$ и ${{R}_{{m \times n}}}$, осуществив QR-разложение (2.16) матрицы ${{A}_{{m \times n}}}$.

Шаг 4. Вычислить канонизацию (1.4), сформировав матрицы $\bar {A}_{{n \times (n - r)}}^{R}$, $\bar {A}_{{(m - r) \times m}}^{L}$, $\tilde {A}_{{n \times r}}^{R}$, $\tilde {A}_{{r \times m}}^{L}$ и ${{\tilde {A}}_{{n \times m}}}$ в соответствии с формулами (2.18)–(2.22). Перейти к шагу 11.

Шаг 5. Вычислить матрицы ${{L}_{{m \times n}}}$ и ${{Q}_{{n \times n}}}$, осуществив LQ-разложение (2.27) матрицы ${{A}_{{m \times n}}}$.

Шаг 6. Вычислить канонизацию (1.4), сформировав матрицы $\bar {A}_{{n \times (n - r)}}^{R}$, $\bar {A}_{{(m - r) \times m}}^{L}$, $\tilde {A}_{{n \times r}}^{R}$, $\tilde {A}_{{r \times m}}^{L}$ и ${{\tilde {A}}_{{n \times m}}}$ в соответствии с формулами (2.28)–(2.32). Перейти к шагу 11.

Шаг 7. Вычислить матрицы ${{P}_{{m \times m}}}$, ${{L}_{{m \times m}}}$ и ${{U}_{{m \times m}}}$, осуществив LU-разложение (2.33) матрицы ${{A}_{{m \times n}}}$.

Шаг 8. Вычислить канонизацию (1.4), сформировав матрицы $\bar {A}_{{n \times (n - r)}}^{R}$, $\bar {A}_{{(m - r) \times m}}^{L}$, $\tilde {A}_{{n \times r}}^{R}$, $\tilde {A}_{{r \times m}}^{L}$ и ${{\tilde {A}}_{{n \times m}}}$ в соответствии с формулами (2.35)–(2.39). Перейти к шагу 11.

Шаг 9. Вычислить матрицы ${{U}_{{m \times m}}}$, ${{{\Sigma }}_{{m \times n}}}$ и ${{V}_{{n \times n}}}$, осуществив сингулярное разложение (2.5) матрицы ${{A}_{{m \times n}}}$.

Шаг 10. Вычислить канонизацию (1.4), сформировав матрицы $\bar {A}_{{n \times (n - r)}}^{R}$, $\bar {A}_{{(m - r) \times m}}^{L}$, $\tilde {A}_{{n \times r}}^{R}$, $\tilde {A}_{{r \times m}}^{L}$ и ${{\tilde {A}}_{{n \times m}}}$ в соответствии с формулами (2.11)–(2.15). Перейти к шагу 13.

Шаг 11. Вычислить число обусловленности $\kappa \left( A \right)$.

Шаг 12. Если ${{\left( {\kappa \left( A \right)} \right)}^{{ - 1}}} < u \cdot \left\| A \right\| \cdot \max \left( {m,n} \right)$, то перейти к шагу 9, иначе перейти к шагу 13.

Шаг 13. Конец алгоритма.

Следует отметить, что в большинстве пакетов прикладных программ (MATLAB, Wolfram Mathematica, Mathcad и т.д.) для вычисления матричных разложений, используемых предлагаемым алгоритмом, существуют эффективные вычислительные процедуры. Отметим, что шаг 11 предлагаемого алгоритма предполагает вычисление числа обусловленности $\kappa \left( A \right)$ матрицы ${{A}_{{m \times n}}}$, а шаг 12 – сравнение ${{\left( {\kappa \left( A \right)} \right)}^{{ - 1}}}$ с допустимым уровнем точности $u \cdot \left\| A \right\| \cdot \max \left( {m,n} \right)$, зависящим от машинной точности $u$, нормы исходной матрицы $\left\| A \right\|$ и максимальной размерности матрицы $\max \left( {m,n} \right)$ (см., например, [15] и [16]).

Далее предлагается способ вычисления числа обусловленности, возникающий в результате вычисления канонизации (1.4), позволяющий не осуществлять явного обращения исходной матрицы. Данный способ может быть особенно полезным в случаях, когда исходная матрица является прямоугольной или вырожденной.

3. РЕЗУЛЬТАТЫ

Для иллюстрации работы предлагаемого алгоритма осуществим канонизацию выборки объемом $N = 100\,000$ в общем случае прямоугольных матриц. Количества строк и столбцов матриц являются независимыми дискретными случайными величинами, равномерно распределенными в интервале $\left[ {2,10} \right]$. Элементы матриц также являются независимыми дискретными случайными величинами, равномерно распределенными в интервале $\left[ { - 10,10} \right]$. На фиг. 1 представлены распределения чисел обусловленности задач обращения и вычисления канонизации матрицы для данной выборки, вычисляемые по формулам (2.40), (2.41) и (2.43) соответственно; на оси абсцисс отложено число обусловленности, на оси ординат – количество матриц с соответствующим числом обусловленности. По гистограммам видно, что предлагаемая оценка числа обусловленности (2.43) практически идентична оценкам (2.40), (2.41). Отметим, что количество матриц с большим числом обусловленности ($\kappa \left( A \right) > {{10}^{3}}$) достаточно невелико, однако не равно нулю. На фиг. 2 показаны значения относительной ошибки канонизации, определяемые по следующей формуле:

Видно, что относительная ошибка $\delta {{\tilde {A}}_{{n \times m}}}$ не превышает $u \cdot \max \left( {m,n} \right) \cdot \kappa \left( {{{A}_{{m \times n}}}} \right)$ при любых размерах и значениях числа обусловленности матрицы ${{A}_{{m \times n}}}$.

Далее рассмотрим процесс канонизации плохо обусловленной обратной по отношению к матрице Гильберта матрицы размера $5 \times 5$

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

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

Поскольку матрица $T$ является невырожденной, матрицы ${{\bar {T}}^{L}}$ и ${{\bar {T}}^{R}}$ не существуют. Числа обусловленности, вычисленные по формулам (2.40) и (2.42), совпадают и равны $4.7661 \times {{10}^{5}}$; максимальная ошибка составляет $u \cdot \max (m,n) \cdot \kappa ({{A}_{{m \times n}}}) = 2.9104 \times {{10}^{{ - 10}}}$, а относительная ошибка равна $\delta {{\tilde {A}}_{{n \times m}}} = 6.5157 \times {{10}^{{ - 12}}}$.

Теперь рассмотрим процесс канонизации прямоугольной матрицы

Осуществляя ее канонизацию в соответствии с предлагаемым алгоритмом, получаем следующие матрицы: Поскольку все строки исходной матрицы $A$ линейно независимы, ее левый делитель нуля (матрица ${{\bar {A}}^{L}}$) не существует.

Числа обусловленности, вычисляемые по формулам (2.41) и (2.43), оказались приблизительно равны – разница между ними оказалась равной $8.8818 \times {{10}^{{ - 16}}}$, кроме того, поскольку канонизация данной матрицы осуществлялась методом LQ-разложения, оценка ее числа обусловленности является точной, то есть совпадает с числом обусловленности, полученным по формуле (2.42). Максимальная ошибка составляет $u \cdot \max (m,n) \cdot \kappa ({{A}_{{m \times n}}}) = 4.4409 \times {{10}^{{ - 15}}}$, а относительная ошибка равна $\delta {{\tilde {A}}_{{n \times m}}} = 7.2075 \times {{10}^{{ - 16}}}$.

4. ОБСУЖДЕНИЕ

Определим области применения предлагаемого алгоритма канонизации. Отметим, что алгоритмы матричных разложений являются устойчивыми и могут быть применены к матрицам любого размера. Кроме того, из работы [14] известно, что для решения плохо обусловленных СЛАУ целесообразно применять алгоритм сингулярного разложения с последующей аппроксимацией плохо обусловленной матрицы ${{{\Sigma }}_{{m \times n}}}$ хорошо обусловленной матрицей более низкого ранга. Таким образом, с помощью предлагаемого алгоритма канонизации можно решать переопределенные, недоопределенные, вырожденные, а также плохо обусловленные СЛАУ. Однако следует отметить, что вычислительная сложность такого разложения, например, методом Голуба–Райнша (Golub-Reinsch method) [19] составляет порядка $4{{m}^{2}}n + 8m{{n}^{2}} + 9{{n}^{3}}$, а методом R-SVD – порядка $4{{m}^{2}}n + 22{{n}^{3}}$ операций с плавающей точкой, что существенно медленнее стандартных алгоритмов вычисления QR- или LU-разложений. Так, например, алгоритм QR-разложения, основанный на преобразованиях Хаусхолдера, также является устойчивым и требует лишь порядка $2m{{n}^{3}} - 2{\text{/}}3{{n}^{3}}$ операций, а для квадратных матриц LU-разложение позволяет выполнить канонизацию существенно быстрее, поскольку требует лишь порядка $2{\text{/}}3{{m}^{3}}$ операций с плавающей точкой [15]. Таким образом, для вычисления канонизации представляется целесообразным использовать сингулярное разложение только в случаях, когда требуется повышенная точность результата канонизации или матрица является плохо обусловленной. Во всех остальных случаях желательно использовать более эффективные с вычислительной точки зрения алгоритмы, например, QR- и LU-разложения.

Следует отметить, что хотя LU-разложение может применяться и к прямоугольным матрицам, более целесообразным представляется в таких случаях использовать QR- или LQ-разложения, поскольку исходная прямоугольная матрица по определению является вырожденной и требует более точных методов вычисления базисов фундаментальных подпространств матрицы ${{A}_{{m \times n}}}$.

Отдельно укажем на тот факт, что формулы (2.18)–(2.22) и (2.28)–(2.32) напрямую не требуют обратимости матриц ${{R}_{{m \times n}}}$ и ${{L}_{{m \times n}}}$ соответственно, что позволяет осуществлять на их основе канонизацию матриц произвольного ранга и размера. Кроме того, хотя в формулах (2.18), (2.20), (2.29) и (2.31) формально присутствует обращение матриц ${{R}_{{r \times r}}}$ и ${{L}_{{r \times r}}}$, оно, с учетом того, что они являются треугольными, на практике сводится к применению эффективных алгоритмов прямой и обратной подстановок. Как и в случае с QR-разложением, в формулах (2.35)–(2.39) формально присутствуют операции обращения матриц ${{L}_{{m \times m}}}$ и ${{U}_{{r \times r}}}$, однако обе матрицы являются треугольными обратимыми и, следовательно, для вычисления можно также применять эффективные алгоритмы прямой и обратной подстановок.

Отметим, что на практике также полезным свойством канонизации на основе QR- или LQ-разложений является возможность получения ортонормированного левого или правого делителя нуля и канонизатора, позволяющих более тонко контролировать норму решения (1.10), варьируя произвольную матрицу ${{\eta }_{{\left( {n - r} \right) \times p}}}$. Также заметим, что формула оценки числа обусловленности (2.43) с вычислительной точки зрения существенно проще формулы (2.42), поскольку не требует перемножения матриц $\tilde {A}_{{n \times r}}^{R}$ и $\tilde {A}_{{r \times m}}^{L}$. Равенство в (2.43) достигается в случае, если один из канонизаторов, правый или левый, является ортонормированным. Учитывая формулы (2.21) и (2.30), можно видеть, что данное условие выполняется при вычислении канонизации методами QR- или LQ-разложений, в результате которых оказывается унитарной матрица $\tilde {A}_{{r \times m}}^{L}$ или $\tilde {A}_{{n \times r}}^{R}$ соответственно. Если для оценки обусловленности задачи канонизации используется формула (2.42), то непосредственно перед ее вычислением возникает необходимость вычисления сводного канонизатора ${{\tilde {A}}_{{n \times m}}}$ путем перемножения правого и левого канонизаторов, что несколько снижает ее вычислительную эффективность перед оценкой (2.43). Таким образом, алгоритм канонизации матрицы, основанный на QR- или LQ-разложении, позволяет использовать более эффективную с вычислительной точки зрения оценку обусловленности (2.43) без потери точности в сравнении с оценкой (2.42).

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

В данной работе предложены новые методы нахождения канонизации матрицы, основанные на ее LU-, QR- или LQ-разложении, кроме того, представлены формулы, позволяющие осуществлять оценку обусловленности задачи канонизации, не осуществляя явного обращения исходной матрицы. Также на основе предлагаемых методов разработан и реализован на языке программирования MATLAB эффективный алгоритм вычисления канонизации матриц, позволяющий существенно повысить как устойчивость, так и быстродействие вычислительного процесса. Показано, что наиболее эффективным методом канонизации квадратных, возможно вырожденных матриц, является метод, основанный на LU-разложении произвольных прямоугольных матриц – метод, основанный на QR- или LQ-разложении, а плохо обусловленных матриц – метод, основанный на сингулярном разложении. Показаны точность и эффективность оценки обусловленности матрицы, получаемой на основе вычисления норм ее левого и правого канонизаторов. Доказано, что наиболее эффективную с вычислительной точки зрения оценку числа обусловленности, не теряющую при этом своей точности, можно получить, вычисляя канонизацию методами, основанными на QR- или LQ-разложениях.

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