Обратная матрица с сужением: как красиво посчитать?

[B@R5uk]

У меня есть симметричная положительно определённая матрица G (полученная умножением транспонированной матрицы A на себя, но не суть):

$$G=G^T,\qquad G\in\mathbb{R}^{n\times n}$$

а так же матрица C небольшого количества ортонормированных столбцов, образующих вместе со столбцами матрицы D ортогональную матрицу B базисных столбцов всего пространства (I — единичная):

$$B^TB=BB^T=I_n,\qquad B\in\mathbb{R}^{n\times n},\qquad B=\bigl(C\;\;D\bigr),\qquad C\in\mathbb{R}^{n\times p},\qquad D\in\mathbb{R}^{n\times(n-p)},\qquad p\ll n$$

То есть, данная по условию матрица C задаёт одно подпространство, а матрица D — другое, ортогональное первому, но в котором надо работать.

Вопрос такой: как красиво посчитать матрицу:

$$H=D\left(D^TGD\right)^{-1}D^T$$

Дело в том, что исходными данными для расчёта искомой матрицы H являются матрицы C и G. Но в явную формулу для H входит промежуточная величина — матрица D, которая рассчитывается по матрице C, причём неоднозначно. Если взять произвольную обратимую матрицу Q и с помощью неё от D перейти к новой D', то формула не изменится:

$$D'=DQ$$

$$H=D\left(D^TGD\right)^{-1}D^T=DQQ^{-1}\left(D^TGD\right)^{-1}Q^{-T}Q^TD^T=$$ $$=D'(Q^TD^TGDQ)(DQ)^T=D'\left(D'^TGD'\right)^{-1}D'^T$$

То есть матрица H инвариантна относительно выбора D, причём даже не обязательно, чтобы последняя была ортонормированная. Хотелось бы от этой избыточности в расчёте избавиться. Так же факт того, что число столбцов p матрицы C много меньше размерности пространства n, намекает на то, что использование матрицы D в расчёте H довольно нерационально.

Я уже долго пытаюсь придумать, как выразить H только через G и C, но ничего вразумительного пока не получается. Один тривиальный случай, когда матрица G является единичной решается легко:

$$D\left(D^TI_nD\right)^{-1}D^T=D\left(D^TD\right)^{-1}D^T=D\left(I_{n-p}\right)^{-1}D^T=DI_{n-p}D^T=DD^T=I_n-CC^T$$

Последнее равенство следует из разложения для B:

$$I_n=BB^T=\bigl(C\;\;D\bigr)\bigl(C\;\;D\bigr)^T=CC^T+DD^T$$

Ещё я пробовал посмотреть что получится для H при попытке представить её в виде спектрального разложения, но что-то ничего путного не приходит в голову. Ну, кроме того, что у неё есть кратное (кратности p) нулевое собственное значение и соответствующее ему подпространство, задаваемое столбцами матрицы C.

Следующее, что мне приходит в голову — это попытаться прикрутить формулу Шермана-Моррисона для модифицированной обратной матрицы. Но эта формула подразумевает, что обе матрицы — исходная и обратная — являются невырожденными, что в моём случае неверно: уменьшение размерности G при переходе к H означает зануление пачки строк и столбцов. Это означает, что эту формулу надо сначала допилить до чего-то псевдообратного в духе Мура-Пенроуза.

Подскажите, пожалуйста, что ещё можно придумать, чтобы выразить H через C и G? Может я пропустил/не заметил что-то очевидное, например в спектральном разложении?

[B@R5uk]

Неужели никто ничего путного не сможет посоветовать?

[magnetic_balls]

$D^+=D^T$

$$D(D^TGD)^+D^T=DD^+G^+(D^T)^+D^T=DD^TG^+DD^T=(I-CC^T)G^+(I-CC^T)$$

[B@R5uk]

magnetic_balls, спасибо за идею! К сожалению, ваша формула не работает. Простая и быстрая проверка:

Используется синтаксис Matlab M

%   dxdy.ru/topic162142.html
clc
clear all

n = 6;
m = 2;

g = randn (n);
g = g + g';

b = randn (n);
[b c d] = svd (b);
b = b * d;
c = b (:, 1 : m);
d = b (:, m + 1 : end);

h1 = d * (d' * g * d) ^ -1 * d';
h2 = eye (n) - c * c';
h2 = h2 * g ^ -1 * h2;

disp (h1 - h2)
 

Результат значительно больше машинного эпсилон:

Используется синтаксис Text

   -0.0829    0.1539    0.0427    0.0403   -0.1003    0.0144
    0.1539   -0.3350   -0.1000   -0.0735    0.2078   -0.0404
    0.0427   -0.1000   -0.0307   -0.0202    0.0608   -0.0132
    0.0403   -0.0735   -0.0202   -0.0196    0.0481   -0.0066
   -0.1003    0.2078    0.0608    0.0481   -0.1307    0.0234
    0.0144   -0.0404   -0.0132   -0.0066    0.0234   -0.0063
 

Я повникал в ваш вывод и понял, что затык вот в этой формуле, что вы используете:

$$(AB)^+=B^+A^+$$

Она иногда работает, а иногда нет:

Используется синтаксис Matlab M

%   dxdy.ru/topic162142.html
clc
clear all

n = 4;
m = 5;
k = 3;
l = 6;

a1 = randn (n, k);
b1 = randn (k, m);
c1 = pinv (a1 * b1);
d1 = pinv (b1) * pinv (a1);

a2 = randn (n, l);
b2 = randn (l, m);
c2 = pinv (a2 * b2);
d2 = pinv (b2) * pinv (a2);

disp (c1 - d1)
disp (' ')
disp (c2 - d2)
 

Результат:

Используется синтаксис Text

  1.0e-014 *
    0.1110    0.0527    0.1499   -0.0944
   -0.0222    0.0437    0.0555    0.0111
   -0.0173    0.0042   -0.0045    0.0250
    0.0500    0.0389    0.1277   -0.0833
    0.0069    0.0087    0.0271    0.0028
 
   -0.0809   -0.1437    0.0328   -0.0998
   -0.4690   -1.0215   -0.1138   -0.5219
    0.1982    0.5721    0.2747    0.1782
   -0.4041   -1.1803   -0.5826   -0.3591
   -0.5877   -1.4426   -0.4047   -0.6051
 

В первом случае разность порядка машинного эпсилон, то есть формула работает, когда же внутренний размер матриц произведения больше внешних, то всё ломается. В любом случае, magnetic_balls, ещё раз большое СПАСИБО за идею! Теперь мне хотя бы есть в каком направлении копать.

И сразу вопрос ко всем: посоветуйте, пожалуйста, литературу, в которой обсуждается псевдообращение произведения двух матриц. Что-то я на это правило нигде ещё не натыкался.

-- 15.01.2026, 14:02 --

Наткнулся вот на этот вопрос-ответ на MSE. На первый взгляд результаты от туда можно прикрутить к моей задаче. Надо только вникнуть, что за уличная магия там творится.

[B@R5uk]

Обозначив проектор на подпространство ${}_{{}_.}\operatorname{span}(C)$ как

$$P_C=CC^T$$

и воспользовавшись формулами по ссылке выше (не особо вникая что там происходит), у меня получилось что-то такое:

$$H=G^{-1}-G^{-1}P_C\left(G^{-1/2}P_C\right)^+G^{-1/2}$$

При этом требуется, чтобы G была положительно определена. При экспериментальной проверке эта формула даже работает:

Используется синтаксис Matlab M

%   dxdy.ru/topic162142.html
clc
clear all

n = 6;
m = 5;

g = randn (n);
g = g + g';
[v g] = eig (g);
g = abs (diag (g));
g_ir = v * diag (1 ./ sqrt (g)) * v';
g_i  = v * diag (1 ./       g ) * v';
g    = v * diag (           g ) * v';

b = randn (n);
[b c d] = svd (b);
b = b * d;
c = b (:, 1 : m);
d = b (:, m + 1 : end);

h1 = d * (d' * g * d) ^ -1 * d';
p_c = c * c';
h2 = g_i - g_i * p_c * pinv (g_ir * p_c) * g_ir;

disp (h1 - h2)
 

К сожалению, необходимость вычислять корень из матрицы, а потом ещё и псевдообратную от произведения сводит все усилия по снижению вычислительных трудозатрат на ноль. Ну, хотя бы получилась интересная и необычная формула. Может быть её ещё как-то можно вывернуть?

-- 15.01.2026, 20:32 --

Численный эксперимент показывает, что

$$P_C\left(G^{-1/2}P_C\right)^+$$

можно упростить до

$$\left(G^{-1/2}P_C\right)^+$$

потому что эти два выражения равны. Пока не понимаю, почему это происходит и всегда ли оно так. Прямое какое-то волшебство.

[B@R5uk]

B@R5uk в сообщении #1714903 писал(а):

Пока не понимаю, почему это происходит и всегда ли оно так.

Кажись, въехал. Псевдообращение произведения AP (где P — это ортогональный проектор на некое подпространство) можно выразить через его "экономное" сингулярное разложение (матрица Σ с тильдой "урезана" до ненулевых значений на диагонали):

$$AP=\tilde{U}\tilde{\Sigma}\tilde{V}^T$$

$$(AP)^+=\tilde{V}\tilde{\Sigma}^{-1}\tilde{U}^T$$

Умножение на P справа соответствует проектированию всех строк матрицы на подпространство этого проектора (умножение слева проектирует туда столбцы). Поэтому матрица V с тильдой в разложении выше является как раз базисом подпространства проектора:

$$P=\tilde{V}\tilde{V}^T$$

Это верно, даже если проектор до этого задавался каким-то другим базисом С, и справедливо равенство:

$$P=CC^T=\tilde{V}\tilde{V}^T$$

Но это не значит, что матрицы C и V с тильдой равны (хотя они связаны некоторой ортогональной матрицей). При построении псевдообращения произведения A на P выше подпространства столбцов и строк меняются местами (что в буквальном смысле очевидно из формулы: базисы этих подпространств — матрицы U и V — меняются местами). Поэтому умножение на P слева приводит к повторному проектированию уже спроектированных векторов и ничего не меняет:

$$P(AP)^+=(AP)^+$$

Это верно всегда, даже если матрица A при умножении на P теряет в ранге.

-- 15.01.2026, 21:37 --

В результате, моя искомая формула редуцируется до этого:

$$H=G^{-1}-G^{-1}\left(G^{-1/2}CC^T\right)^+G^{-1/2}$$

-- 15.01.2026, 21:39 --

А в силу симметричности матриц G и H её ещё можно записать так:

$$H=-G^{-1/2}\left(CC^TG^{-1/2}\right)^+G^{-1}+G^{-1}$$

[magnetic_balls]

Проверьте, пожалуйста, $H=((1-CC^T)G(1-CC^T))^+$. Если не ошибся опять, напишу вывод.

[B@R5uk]

Да, это верно. Более того, выражение под скобками

$$(1-CC^T)G(1-CC^T)$$

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

$$H=\left((1-CC^T)G(1-CC^T)+CC^T\right)^{-1}-CC^T$$

Подробнее я расписал здесь. Там же мне эту формулу и подсказали.

B@R5uk в сообщении #1714877 писал(а):

И сразу вопрос ко всем: посоветуйте, пожалуйста, литературу, в которой обсуждается псевдообращение произведения двух матриц.

Шикарно про это дело расписано в книге Stephen L. Campbell, Carl D. Meyer - Generalized Inverses of Linear Transformations в первой же главе. Не знаю, если в русскоязычной литературе что-нибудь на столько же годное и понятное.

(Оффтоп)

Если нет, то у меня возникают серьёзные опасения за отечественную науку. На одних хороших семинаристах и лекторах далеко не уедешь, нужны и хорошие учебники.

magnetic_balls, спасибо за помощь и подсказки!

[magnetic_balls]

Всё равно выложу свой вывод(можно ли оформить $\TeX$ лучше?)

(Оффтоп)

$$D(D^TGD)^+D^T=B(0I)^T((0I)B^TGB(0I)^T)^+(0I)B^T=B(0I)^TX^+(0I)B^T=$$ $$=B\begin{pmatrix} 0 & 0 \\ 0 & X^+ \end{pmatrix}B^T=B\begin{pmatrix} 0 & 0 \\ 0 & X \end{pmatrix}^+B^T=B(\begin{pmatrix} 0 & 0 \\ 0 & I \end{pmatrix}X\begin{pmatrix} 0 & 0 \\ 0 & I \end{pmatrix})^+B^T=$$ $$=B((0I)^T(0I)B^TGB(0I)^T(0I))^+B^T=(B^T^+(0I)^TD^TGD(0I)B^+)^+=(DD^TGDD^T)^+$$

Воспользовался верностью $(AQ)^+= Q^+A^+$ для ортогональной $Q$.