Можно ли частично обратить метод Гаусса?

[B@R5uk]

У меня имеется некоторая матрица

$$A\in\mathbb{R}^{m\times n},\quuad\qquad m>n$$

Из которой мне надо отобрать

$$m'<n$$

линейно независимых строк. Гарантировано, что этих строк ровно столько, и даже что их можно выбрать не единственным способом. Желательно, однако, чтобы выбираемые строки имели как можно меньший индекс, то есть предпочтение отдаётся самым верхним строкам матрицы A (до тех пор, пока они линейно независимы, зависимые — пропускаются).

Для реализации этого я решил использовать модифицированный метод Гаусса, который заключается в следующем. В процессе работы, строится матрица B, которая постепенно приближается к единичной матрице (с точностью до перестановки столбцов). Каждая строка характеризуется индексом некоторого ведущего элемента, который делается равным 1. Если в каком-то столбце есть ведущий элемент, то это будет единственный ненулевой элемент этого столбца. Чтобы достичь этого свойства матрицы B, на каждом шаге из матрицы A берётся очередная строка, для каждой её ненулевой координаты, находящейся в столбце с ведущим элементом матрицы B из взятой строки вычитаются строки матрицы B, помноженные на элементы строки из A, так, чтобы эти элементы обратились в 0 в столбцах с ведущими элементами B:

$$a'=a-\sum\limits_{k=1}^{r}a_{p(k)}b_k,\qquad\qquad b_{k,p(k)}=1$$

Здесь a — очередная строка матрицы A, a' — её модификация, a с индексом — элемент строки a, b с индексом — строки матрицы B, у которых в позиции $p(k)$ находится ведущий элемент (равный единице). Число r — это количество строк в матрице B на начало текущей итерации. При этом теряется хорошее свойство разреженности матрицы A, но тут уж ничего не поделаешь.

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

$$b_{r+1}=\frac{1}{a'_{max}}\;a',\qquad\qquad a'_{p(r+1)}=a'_{max}$$

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

$$b'_k=b_k-\frac{b_{k,p(r+1)}}{a'_{max}}\;a',\qquad\qquad 1\le k\le r$$

В результате получается что-то в этом духе:

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

    1.0000   -0.2441         0         0         0         0         0         0   -1.0000    0.2441
         0    0.0849         0         0         0         0         0    1.0000         0   -1.0849
         0    0.1387         0    1.0000         0         0         0         0         0   -1.1387
         0   -0.0563         0         0    1.0000         0         0         0   -1.0000    0.0563
         0   -0.4150         0         0         0         0    1.0000         0   -1.0000    0.4150
         0   -0.6874         0         0         0    1.0000         0         0         0   -0.3126
         0   -0.8487    1.0000         0         0         0         0         0   -1.0000    0.8487
 

В 1 и 3—8 столбцах находятся ведущие элементы, остальные столбцы (2, 9, 10) соответствуют свободным переменным.

По окончании всего этого производится проверка корректности и частенько оказывается, что одна строка была выбрана неудачно и её надо удалить и заменить какой-то другой строкой матрицы A. В связи с этим у меня вопросы:

1. Действительно ли то, что я делаю, является разновидностью метода Гаусса?
2. Можно ли удалить неудачную строку из матрицы B, не пересчитывая всю матрицу B с самого начала (другими словами, за время хотя бы $O(n^2)$, а не $O(n^3)$)?
3. Если да, то какую дополнительную информацию для этого надо хранить?
4. Существует ли какой-нибудь другой, возможно, более рациональный подход к этой задаче?

[dgwuqtj]

Да, это метод Гаусса. Обратить его, может, и можно, только при выполнении арифметических операций в обратном порядке будут накапливаться ошибки округления. Насколько это критично?

[B@R5uk]

dgwuqtj, спасибо. На счёт ошибок — не критично. Во-первых не так страшен чёрт, как его малюют, а во-вторых я пока исследую саму возможность проделывать описанное (и спрошенное). Когда формулы будут готовы в явном виде, можно будет и ошибки прикинуть.

Я так понимаю, чтобы была возможность быстро удалить строку, надо поддерживать в процессе вычислений некоторый набор "волшебных" столбцов. Каждый столбец соответствует строке, и для удаления любой из строк матрицы B нужно к ней добавить соответствующий этой строке столбец, помноженный на саму удаляемую строку:

$$B'=B+c_kb_k$$

Где k — индекс удаляемой строки, b и c с индексом — это, соответственно, удаляемая строка и её "волшебный" столбец. В матрице B' после этой операции удаляемая строка просто вычёркивается. Результат должен быть тот же самым, как если бы начали с матрицы A, выбрав те же самые строки (кроме удалённой), и выбрав в этих строках те же самые ведущие элементы. Они могут оказаться в этом случае не самыми лучшими, но цель будет достигнута.

Я верно рассуждаю?

[zgemm]

Первое, что приходит в голову, это или LU с полным пивотированием, или, как некоторые умники его называют - Адаптивная Корсс Аппроксимация.

Если на пальцах - надо делать так.

1. Вначале в каждой строке найти максимальный по модулю элемент, помножить это значение на какую-то весовую функцию, в которой учитывается близость этой строки к верху.
2. Выполнить один шаг LU для этой строки и выбрать столбец с максимумом по норме.
3. Найти дополнение по Шуру на все оставшееся.
4. Повторить для дополнения по Шуру все то же самое с первого шага и так делать пока не получилось столько строк, сколько вы хотите.

Все другие методы будут существенно более вычислительно сложными.

[B@R5uk]

zgemm, мне не понятен момент: что в описанном вами алгоритме позволяет удалять строки из результата после его завершения?

zgemm в сообщении #1687386 писал(а):

существенно более вычислительно сложными

Ну, метод Гаусса имеет кубическое время $O(n^3)$. Если удастся дополнить его структурой, которая позволит удалять строки за квадратичное время $O(n^2)$, то я буду доволен.

[zgemm]

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

zgemm, мне не понятен момент: что в описанном вами алгоритме позволяет удалять строки из результата после его завершения?

zgemm в сообщении #1687386 писал(а):

существенно более вычислительно сложными

Ну, метод Гаусса имеет кубическое время $O(n^3)$. Если удастся дополнить его структурой, которая позволит удалять строки за квадратичное время $O(n^2)$, то я буду доволен.

Гаусс и LU - это идентичные алгоритмы. Вам надо только правильно выбирать элемент, который вы будете исключать и я об этом написал. Чтобы после исключения одной пары строка-столбец двигаться дальше, вам надо все остальные элементы пересчитать, эта процедура конечно квадратична. Сделать что-то за меньшее время - практически нельзя, так как у LU с выбором это $3n^2$, а теоретически возможный минимум - как раз треть от этого.

Да почитайте хотя бы краткий курс численного анализа Тыртышникова, там не сильно заумно об этом написано.

[B@R5uk]

zgemm в сообщении #1687415 писал(а):

Вам надо только правильно выбирать элемент, который вы будете исключать и я об этом написал.

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

[B@R5uk]

Хорошо, кажется, я сообразил, что надо с Гауссом делать. Вот, мой модифицированный алгоритм преобразует урезанную матрицу A в матрицу B:

$$A_{red}\xrightarrow{algo}B=\left(\begin{array}{cc}I_{m'}&B_{side}\end{array}\right)\Pi$$

$$A_{red},B\in\mathbb{R}^{m'\times n},\qquad\qquad B_{size}\in\mathbb{R}^{m'\times(n-m')},\qquad\qquad\Pi\in\mathbb{R}^{n\times n}$$

где I — это единична матрица размером m', а Π — перестановочная матрица, символизирующая тот факт, что ведущие элементы в B выбираются где придётся. Необходимые мне "волшебные" столбцы можно извлечь из "псевдообратной" матрицы С, которая получится, если матрицу A дополнить единичной матрицей и заставить алгоритм работать над расширенной матрицей:

$$\left(\begin{array}{c|c}A_{red}&I_{m'}\end{array}\right)\xrightarrow{algo}\left(\begin{array}{c|c}B&C\end{array}\right)$$

Матрица C будет следить за историей модификации матрицы A: куда какая строка была добавлена с каким коэффициентом. В самом деле:

$$CA_{red}=B$$

Теперь, чтобы удалить k-ю строку матрицы A и всё её влияние из матрицы B мне нужно вычесть k-ю строку матрицы B из остальных с коэффициентами, с которыми она в них входит. Эти коэффициенты получаются из k-го столбца матрицы C:

$$B'=B-\frac{1}{c_{kk}}c_kb_k$$

Здесь $c_{kk}$ — это k-й диагональный элемент матрицы C, c с индексом — её k-й столбец и b с индексом — удаляемая срока матрицы B. После проделанной операции, матрица B потеряет все следы удаляемой строки, в том числе на её месте будут стоять нули, и они могут быть просто вычеркнуты.

За одно, этот подход явно сигнализирует о ситуации (существование которой верно заметил dgwuqtj), когда старый выбор ведущих элементов будет неудачным для новой матрицы B' без удалённой строки. А именно, когда k-й диагональный элемент матрицы C (на который осуществляется деление) окажется близок к нулю или даже просто заметно меньше остальных элементов k-го столбца:

$$\kappa_{cond}=\max\limits_{1\le l\le m'}^{}\left|\frac{c_{lk}}{c_{kk}}\right|>\kappa_{threshold}$$

Слишком большие величины добавляются к элементам матрицы B, что в силу конечной точности арифметики с плавающей запятой приведёт к потере информации. Пороговая величина может быть 2, 10, 32, 100 или что-нибудь в этом духе в зависимости от задачи. Что в этой ситуации делать, надо подумать. Ниже программа Matlab, тестирующая работоспособность выше сказанного.

Файл study_gauss_reversal.m:

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

%   dxdy.ru/topic160459.html
%   2025, May, B@R5uk

clc
clearvars
format compact
%format long
format short

num_cols = 7;
num_rows = 4;
num_rem = 2;

a = randn (num_rows, num_cols);
[b, c, p] = study_gauss_reversal_func (a);

aa = a ([1 : num_rem - 1 num_rem + 1 : end], :);
pp = p ([1 : num_rem - 1 num_rem + 1 : end]);

[bb1, cc, pp] = study_gauss_reversal_func (aa, pp);

b_row = b (num_rem, :) / c (num_rem, num_rem);
bb2 = b - c (:, num_rem) * b_row;
bb2 = bb2 ([1 : num_rem - 1 num_rem + 1 : end], :);

disp (max (abs (c (:, num_rem) / c (num_rem, num_rem))))
disp (' ')
disp (bb1 - bb2)
disp (' ')
disp (c * a - b)
 

Файл study_gauss_reversal_func.m:

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

function [b, c, p] = study_gauss_reversal_func (a, p)

num_rows = size (a, 1);

if 1 == nargin
    p = zeros (1, num_rows);
    p_flag = 1;
else
    p_flag = 0;
end

b = zeros (size (a));
c = zeros (num_rows);
for k = 1 : num_rows
    b_row_new = a (k, :);
    b_row_vals = b_row_new (p (1 : k - 1));
    b_row_new = b_row_new - b_row_vals * b (1 : k - 1, :);
    c_row_new =           - b_row_vals * c (1 : k - 1, :);
    c_row_new (k) = 1 + c_row_new (k);
    if p_flag
        p (k) = find (abs (b_row_new) == max (abs (b_row_new)), 1);
    end
    b_row_max = b_row_new (p (k));
    b_row_new = b_row_new / b_row_max;
    c_row_new = c_row_new / b_row_max;
    b_col = b (1 : k - 1, p (k));
    b (1 : k - 1, :) = b (1 : k - 1, :) - b_col * b_row_new;
    c (1 : k - 1, :) = c (1 : k - 1, :) - b_col * c_row_new;
    b (k, :) = b_row_new;
    c (k, :) = c_row_new;
end

end
 

[B@R5uk]

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

Что в этой ситуации делать, надо подумать.

Что-то я туплю. Строки матриц A и B соответствуют друг другу постольку-поскольку, то есть никак не соответствуют; матрица B получается из матрицы A умножением на C:

$$CA=B$$

По этой причине, при удалении k-ой строки $a_k$ матрицы A совсем не обязательно должна удалятся тоже k-я строка матрицы B, даже скорее наоборот. Строка $a_k$ входит во все строки матрицы B (или почти во все), и "степень этого вхождения" определяется коэффициентами k-го столбца матрицы C. Логично, что наибольший по модулю такой коэффициент будет указывать на строку $b_l$ матрицы B, которая "в наибольшей степени содержит" строку $a_k$. Её-то и надо удалять. Такой подход изящно решает проблему деления на нуль плохой обусловленности выбора, сделанного выше (в том числе возможного деления на нуль, если диагональный коэффициент матрицы C выше оказался таковым).

Итак, приведённый ранее модифицированный алгоритм Гаусса по набору строк матрицы A строит матрицы B и C:

$$\left(\begin{array}{c|c}A_{red}&I_{m'}\end{array}\right)\xrightarrow{algo}\left(\begin{array}{c|c}B&C\end{array}\right)$$

Затем он получает индекс k строки матрицы A, которую надо удалить. Чтобы определить, через строку с каким индексом l матрицы B будет производиться удаление, алгоритм должен найти максимальный (по модулю) элемент матрицы C в k-ом столбце:

$$l=\underset{1\le t\le m'}{\arg\max}\lBig|c_{tk}\rBig|$$

Затем строятся модифицированные матрицы B' и C':

$$B'=B-\frac{1}{c_{lk}}c_{\cdot k}b_{l\cdot}$$

$$C'=C-\frac{1}{c_{lk}}c_{\cdot k}c_{l\cdot}$$

Теперь в знаменателе находится максимально возможный элемент, и никакой потери информации не происходит. После этого из матрицы B' удаляется l-я строка (состоящая из одних нулей), а из матрицы C' — k-й столбец и l-я строка (они тоже будут нулевыми).

Демонстрационная программа. Файл study_gauss_reversal.m (Файл study_gauss_reversal_func.m остался прежним и может быть взят из предыдущего поста):

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

%   dxdy.ru/topic160459.html
%   2025, May, B@R5uk

clc
clearvars
format compact
%format long
format short

num_cols = 7;
num_rows = 4;
num_rem_a = 2;

a = randn (num_rows, num_cols);
[b, c, p] = study_gauss_reversal_func (a);

c_col_abs = abs (c (:, num_rem_a));
num_rem_b = find (max (c_col_abs) == c_col_abs, 1);

aa = a ([1 : num_rem_a - 1 num_rem_a + 1 : end], :);
pp = p ([1 : num_rem_b - 1 num_rem_b + 1 : end]);

[bb1, cc1, pp] = study_gauss_reversal_func (aa, pp);

bb2 = b - c (:, num_rem_a) * b (num_rem_b, :) / c (num_rem_b, num_rem_a);
cc2 = c - c (:, num_rem_a) * c (num_rem_b, :) / c (num_rem_b, num_rem_a);

bb2 = bb2 ([1 : num_rem_b - 1 num_rem_b + 1 : end], :);
cc2 = cc2 ([1 : num_rem_b - 1 num_rem_b + 1 : end], :);
cc2 = cc2 (:, [1 : num_rem_a - 1 num_rem_a + 1 : end]);

disp (bb1 - bb2)
disp (' ')
disp (cc1 - cc2)
 

ОК, задача решена. Хоть что-то доводится до конца. Теперь надо подумать, можно ли её прикрутить куда-нибудь, чтобы она делала что-то полезное.

[realeugene]

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

Желательно, однако, чтобы выбираемые строки имели как можно меньший индекс, то есть предпочтение отдаётся самым верхним строкам матрицы A (до тех пор, пока они линейно независимы, зависимые — пропускаются).

Сама идея не очень хорошая, если только поле не конечное. В процессе такого построения строится базис некоторого подпространства строк размерности меньшей, чем исходное число строк. При этом критерий включения или не включения очередной строки в качестве базисной то, что она не лежит или же уже лежит в подпространстве, натянутом на предыдущие строки. Для ответа на этот вопрос очередная строка проецируется в ортогональное дополнение предыдущих строк, и сравнивается длина этой проекции с нулём. Если эта проекция не нулевая, она принимается в качестве очередного базисного вектора. Но из-за вычислительных ошибок плавающей арифметики точное равенство нулю практически невозможно, более того, если эта проекция мала, то она может оказаться направленной куда угодно.

[B@R5uk]

realeugene в сообщении #1687569 писал(а):

Но из-за вычислительных ошибок плавающей арифметики точное равенство нулю практически невозможно...

Учтено (проверяется, что максимальная по модулю координата-элемент строки больше или меньше порога):

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

while num_left > k
    %   Берём очередную строку
    %   ...
    %   Проверяем, осталось ли что-нибудь от строки
    wl_max_val = max (abs (work_line (1 : end - 1)));
    %   Если нет, то переходим к проверке следующей
    if wl_max_val < 1e-8
        continue
    end
    %   ...
end
 

Коэффициент с -8 степенью десятки выбран от балды. В силу характера моих исходных данных он может быть любым в районе где-то от -4 степени до -13 степени, то есть наличие линейной зависимости прямо-таки прыгает в лицо когтями вперёд. В других задачах к вопросу придётся подойти более осторожно, но это не значит, что он не решаем. Просто потребуется критерий, учитывающий больше нюансов задачи.

В коде предыдущих постов я эту проверку не делал, так как эта проблема является оффтопиком к озвученному в теме вопросу.

[realeugene]

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

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

Собственно, насколько мне известно, именно поэтому разработчики пакетов линейной алгебры и наизобретали разные критерии наилучшего выбора следующей строки.