Сингулярное разложение и "самое слабое звено"

[B@R5uk]

Пусть у матрицы A количество столбцов больше числа строк. Если её строки линейно зависимые, то сингулярное разложение позволяет это выявить эту линейную зависимость по равным нулю сингулярным числам и определить ранг матрицы. А есть ли какой-нибудь способ (кроме перебора в лоб и расчёта ранга на каждом шаге) определить какие именно строки надо выкинуть, чтобы новая матрица A' перестала быть вырожденной?

[zgemm]

Есть такой алгоритм - ранг-ревеалинг. По порогу можно остановиться на том наборе столбцов, на котором у вас будет еще подмножество псевдо-невырожденных (по порогу) столбцов.

[B@R5uk]

О! Спасибо за ключевое слово! Первый запрос в Гугле выдал интересную статью. Если я правильно понял содержание, то алгоритмы различных факторизаций, призванные выявить ранг матрицы, по сути своей заключаются в компромиссе между быстродействием (в том числе для важного случая непрерывной модификации исходной матрицы путём добавления/удаления строк/столбцов) и точностью оценки сингулярных значений. Потому что только в сингулярном разложении (по определению) сингулярные числа определены точно. Для других разложений имеется некая теорема Островского для оценки сингулярных чисел:

Теорема. Для матриц

$$A\in\mathbb{C}^{m\times n},\qquad\qquad X\in\mathbb{C}^{n\times n},\qquad\qquad Y\in\mathbb{C}^{m\times m}$$

где X и Y — невырожденные, выполнено

$$\sigma_k\lbig(Y^\ast AX\rbig)=\theta_k\sigma_k\lbig(A),\qquad\qquad 1\le k\le\min(m,\;n)$$

где

$$\sigma_n\lbig(X)\sigma_m\lbig(Y)\le\theta_k\le\sigma_1\lbig(X)\sigma_1\lbig(Y)$$

Другими словами, чем ближе матрицы X и Y к унитарным (или ортогональным в случае действительных чисел), тем ближе числа в диагональной матрице

$$D=Y^\ast AX$$

к реальным сингулярным числам матрицы A.

Мне же быстродействие пока не нужно. Я хочу въехать в суть вопроса и понять, можно ли прикрутить сингулярное разложение (функция svd в Matlab) к тому, что я хочу сделать. Эта функция там хорошо реализована (быстро и точно считается), даже вычисление ранга и псевдообратной матрицы через неё вычисляются (можно прямо открыть код скриптов и проверить).

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

...можно остановиться на том наборе столбцов...

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

[Alex Krylov]

Элиминирование по Гауссу должно дать нужный результат (функция rref() в Matlab, reduced row echelon form, Gauss-Jordan elimination).
См. Gilbert Strang, Introduction to Linear Algebra, Sixth Edition, стр. 93 (3.2 Computing the Nullspace by Elimination: A=CR)

[B@R5uk]

Alex Krylov, я так понимаю, надо посмотреть, какие из строк в результате работы функции rref оказались сплошными нулями, и их и выкинуть из исходной матрицы? Что-то эта функция никак не дискриминирует строки, всегда зануляет последнюю.

[Geen]

Напрягает формулировка про "самое слабое звено". Вот допустим, что есть три строки и третья является линейной комбинацией первых двух - и какую же строку нужно выкинуть? - правильный ответ: можно выкинуть любую....

[B@R5uk]

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

$$A=\left(\begin{matrix}9&8&0&1&0&0\\8&9&0&0&1&0\\1&-1&0&1&-1&0\end{matrix}\right)$$

Видно, что последняя строка — это первая минус вторая. Однако, угол между первыми двумя строками мал, в то время, как последняя почти ортогональна первым двум. Если выкинуть последнюю строку, то число обусловленности будет где-то 12.0416, а если первую или вторую — 6.0635, то есть почти в 2 раза меньше (и значительно лучше: в идеале должно быть 1).

Но мне бы для начала хотя бы понять, как по результатам svd (сингулярного разложения) определить, когда и какие строки нельзя выкидывать (а то ранг упадёт).

[Geen]

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

как по результатам svd (сингулярного разложения) определить

Ну так посчитайте разложение, и убедитесь что никак.

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

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

В любом векторном пр-ве (конечномерном) можно выбрать ортонормированный базис.

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

Если выкинуть последнюю строку, то число обусловленности будет где-то 12.0416, а если первую или вторую — 6.0635, то есть почти в 2 раза меньше (и значительно лучше: в идеале должно быть 1).

И что из этого? Как можно догадаться что за задачу Вы решаете и откуда берёте свои матрицы, и зачем выкидываете строки?
А то может быть третья строка (ввиду малости нормы) это просто случайный шум...

-- 08.05.2025, 21:28 --

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

число обусловленности будет

Кстати, а как Вы его считаете?

[B@R5uk]

Так, вроде чуть-чуть начинает прояснятся. Пусть имеется матрица A, в которой число строк меньше числа столбцов и известно, что одна из строк выражается как линейная комбинация остальных (то есть ранг r матрицы на 1 меньше числа строк):

$$A=\left(\begin{matrix}a_1\\a_2\\\ldots\\a_m\end{matrix}\right)\in\mathbb{R}^{m\times n},\qquad\qquad r+1=m<n$$

где a с индексом — строки матрицы. Находим для неё сингулярное разложение:

$$A=U\Sigma V^{T},\qquad U^TU=UU^T=I_m,\qquad V^TV=VV^T=I_n$$

$$\Sigma=\left(\begin{matrix}\sigma_1&&&&0&\vline&\\&\sigma_2&&&&\vline&\\&&\ddots&&&\vline&0_{m\times(n-m)}\\&&&\sigma_{m-1}&&\vline&\\0&&&&0&\vline&\end{matrix}\right)$$

где I с индексом — это единичная матрица размера, указанного индексом, а 0 с двумя индексами — прямоугольная матрица с нулями. Поскольку в матрице одна строка линейно выражается через другие, последнее её сингулярное значение нулевое:

$$\sigma_m=0$$

Перепишем сингулярное разложение в таком виде:

$$U^TA=\Sigma V^T=B$$

Поскольку все элементы на последней строке матрицы Σ нулевые, последняя строка матрицы B тоже будет содержать одни нули (умножение нулевой строки на произвольную матрицу даёт нулевую строку). В последней строке левой части будет стоять некоторая линейная комбинация строк матрицы A, коэффициентами которой будут элементы последнего столбца матрицы U. Эта линейная комбинация с необходимостью равна нулевой строке:

$$\left(\begin{matrix}u_{1m}\\u_{2m}\\\ldots\\u_{mm}\end{matrix}\right)^T\left(\begin{matrix}a_1\\a_2\\\ldots\\a_m\end{matrix}\right)=0_{1\times n}$$

Теперь, ВНИМАНИЕ!!! Самое интересное. Если какой-то из коэффициентов u с индексами будет нулевым и мы удалим из матрицы A строку, соответствующую этому коэффициенту, то в результате линейная комбинация потеряет одно нулевое слагаемое, и в ней ничего существенно не изменится: оставшиеся строки матрицы, просуммированные всё с теми же коэффициентами дадут всё ту же нулевую строку. То есть строки новой редуцированной матрицы будут всё так же линейно зависимы, её ранг будет на единицу меньше числа строк, которое будет на единицу меньше числа строк исходной матрицы.

Вывод: удаление строки, соответствующей нулевому коэффициенту нулевой линейной комбинации, приводит к уменьшению ранга матрицы A (что плохо: потеря информации). То есть, "сильное звено" svd таки позволяет выявить. Во всяком случае, это происходит в ситуации, когда ровно одна строка выражается через другие. Что происходит в общем случае, я пока не сообразил. Буду благодарен любым наставлениям.

(Оффтоп)

Geen в сообщении #1685403 писал(а):

а как Вы его считаете?

1.  
2. >> help cond 
3.  COND   Condition number with respect to inversion. 
4.     COND(X) returns the 2-norm condition number (the ratio of the 
5.     largest singular value of X to the smallest).  Large condition 
6.     numbers indicate a nearly singular matrix. 
7.   
8.     COND(X,P) returns the condition number of X in P-norm: 
9.   
10.        NORM(X,P) * NORM(INV(X),P).  
11.   
12.     where P = 1, 2, inf, or 'fro'.  
13.   
14.     Class support for input X: 
15.        float: double, single 
16.   
17.     See also rcond, condest, condeig, norm, normest. 
18.  
19.     Reference page in Help browser 
20.        doc cond 
21.  
22. >> 

Geen в сообщении #1685403 писал(а):

И что из этого?

Просто эвристический признак плохой обусловленности матрицы. Как правило, подтверждается вычислениями.

Geen в сообщении #1685403 писал(а):

Как можно догадаться что за задачу Вы решаете и откуда берёте свои матрицы, и зачем выкидываете строки?..

Извиняюсь, что непонятно объяснил ситуацию, и поясняю. Пока ничего не решаю, пытаюсь разобраться, как работает сингулярное разложение, какую информацию оно предоставляет о матрице (кроме уже озвученного ранга), можно ли его применить к задаче удаления линейно зависимых строк. Матрицы — любые модельные, на примере демонстрирующие какие-либо свойства. Шума нет никакого. Совсем нет, несмотря на то, что матрицы для тестирования гипотез я создаю случайным образом.

Выбрасывание строк призвано удалить из матрицы излишнюю информацию, которая в ней уже присутствует и без удаляемых строк. Хотелось бы сделать это таким образом, чтобы число обусловленности результирующей матрицы было по-лучше (ближе к 1).

[Geen]

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

help cond

Враньё, в некотором смысле - для Вашей матрицы обратная не определена. Т.е. это всегда равно бесконечности для любой неквадратной или вырожденной матрицы.

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

Просто эвристический признак плохой обусловленности матрицы.

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

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

Хотелось бы сделать это таким образом, чтобы число обусловленности результирующей матрицы было по-лучше (ближе к 1).

Для любого векторного (под)пространства всегда можно выбрать ортонормированный базис.

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

Выбрасывание строк призвано удалить из матрицы излишнюю информацию

В матрице нет никакой информации, тем более лишней. Информация может появиться в задаче.

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

Поскольку в матрице одна строка линейно выражается через другие, последнее её сингулярное значение нулевое

Вообще говоря, там нулевые идут до $n$...

[zgemm]

Alex Krylov в сообщении #1685391 писал(а):

Элиминирование по Гауссу должно дать нужный результат (функция rref() в Matlab, reduced row echelon form, Gauss-Jordan elimination).
См. Gilbert Strang, Introduction to Linear Algebra, Sixth Edition, стр. 93 (3.2 Computing the Nullspace by Elimination: A=CR)

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

-- 09.05.2025, 01:27 --

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

О! Спасибо за ключевое слово! Первый запрос в Гугле выдал интересную статью. Если я правильно понял содержание, то алгоритмы различных факторизаций, призванные выявить ранг матрицы, по сути своей заключаются в компромиссе между быстродействием (в том числе для важного случая непрерывной модификации исходной матрицы путём добавления/удаления строк/столбцов) и точностью оценки сингулярных значений.

Да, верно, та самая статья, Ник как раз очень понятно это объяснял. Я когда в 90-х об этом алгоритме читал, для меня было многое еще не понятным, а когда как-то в 2000 примерно услышал доклад Ника про РР, и все легло по полочкам.

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

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

Не, не так. Алгоритм находит набор столбцов с общей максимально низкой обусловленностью. Только вот столбцы портятся во время такой работы. Но это - не страшно, так как эти столбцы можно назад восстановить и такое восстановление численно устойчиво, хотя гарантировать то, что столбцы будут идентичны исходным уже нельзя. При очень больших размерностях (примерно от 10 миллионов строк для двойной точности) уже и тут могут накапливаться существенные погрешности из-за машинной арифметики, но против этого тоже есть свои приемы, но они уже не такие простые и больше граничат с какой-то полудоказанной магией, но, обычно хорошо работают.

[B@R5uk]

zgemm, спасибо. А какой быстрый метод факторизации вы посоветуете, если мне нужно постоянно пополнять матрицу строками и удалять их для поддержания низкого числа обусловленности? Ник в своей статье обсудил три с разновидностями, но он делал упор на корректность и точность, а не на быстродействие. У меня, правда, миллионы строк не планируются, максимум сотни.

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

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

А если во время работы поддерживать список элементов, на которых производилось pivoting, то можно же по ним восстановить, какие именно столбцы вошли в расчёт матрицы с минимальной обусловленностью? То есть не нужно будет делать вычисления в обратную сторону. Или я что-то недопонимаю?

(Оффтоп)

Geen в сообщении #1685409 писал(а):

Враньё...

Более 5 000 000 зарегистрированных пользователей Matlab по всему миру (включая учёных и инженеров) не заметили. Наверное, были слишком рассеянными, куда катится мир?

Geen в сообщении #1685409 писал(а):

...для Вашей матрицы обратная не определена. Т.е. это всегда равно бесконечности для любой неквадратной или вырожденной матрицы.

Ну как же так? Есть же такое понятие, как псевдообратная матрица. И теорема Мура-Пенроуза о её единственности. И даже то же SVD-разложение позволяет с лёгкостью её найти:

$$A=U\Sigma V^{T},\qquad\qquad A^+=V\Sigma^+U^{T}$$

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

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

Вывод: удаление строки, соответствующей нулевому коэффициенту нулевой линейной комбинации (ЛК), приводит к уменьшению ранга матрицы A (что плохо: потеря информации).

Если удаление такой строки — это не просто плохо, а ошибочно, то логично прийти к следующей эвристической гипотезе: удаление строк, соответствующих малым по модулю (но не ненулевым) коэффициентам — это плохо, а удаление строк, соответствующих большим по модулю коэффициентам ЛК — это, наоборот, хорошо.

Я решил проверить эту гипотезу численным экспериментом. Для этого неким образом я строю случайную матрицу, нахожу минимальный ненулевой ("плохой" выбор) и максимальный ("хороший" выбор) коэффициенты ЛК, запоминаю их индексы, перебором нахожу, что будет получаться при удалении каждой из строк, упорядочиваю результаты в порядке обусловленности урезанной матрицы и нахожу, на каких позициях в этом упорядоченном списке оказываются "плохой" и "хороший" выборы. Из статистики видно, что удаление строки максимального коэффициента с вероятностью более 90% даёт один из 3-х наилучших результатов:

Код программы (файл study_svd_weak_link_2.m):

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

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

clc
clearvars
format compact
%format long
format short

%   Количество экспериментов для набора статистики
num = 5000;
%   Количество столбцов матрицы
N = 25;
%   Количество строк с линейной зависимостью
M = 12;
%   Количество "свободных" строк (не участвующих в ЛЗ)
L = 3;

%   Гистограмма
hist_data = zeros (M + L, 2);
for k = 1 : num
    %   Построение случайной матрицы
    a = randn (M, N);
    [u s v] = svd (a, 'econ');
    %   С одной линейно зависимой строкой
    s (M, M) = 0;
    %   И несколькими строками, не входящими в линейную зависимость
    a = [u * s * v'; randn(L, N)];
    %   Сингулярное разложение
    [u s v] = svd (a, 'econ');
    %   Линейная комбинация строк матрицы, дающая нулевую строку
    lin_comb = u (:, end);
    %   Гипотетически, выбросить нужно строку, у которой
    %   коэффициент в линейной комбинации строк максимален по модулю
    row_best_hypothetical = find (abs (lin_comb) == max (abs (lin_comb)), 1);
    %   Предположительно, худший случай -- когда индекс в ЛК минимален по
    %   модулю, но ещё не равен нулю
    row_wrst_hypothetical = find (abs (lin_comb) == min (abs (lin_comb (1 : M))), 1);
    %   Информация об обусловленности результата: индекс удаляемой строки и
    %   величина, обратная числу обусловленности матрицы
    cond_info = zeros (M + L, 2);
    %   Перебор всех возможных строк
    for l = 1 : M + L
        %   Матрица без одной строки
        b = a ([1 : l - 1 l + 1 : end], :);
        %   Величина, обратная её числу обусловленности
        cond_info (l, 2) = 1 / cond (b);
        %   Индекс удаляемой строки
        cond_info (l, 1) = l;
    end
    %   Сортировка в порядке убывания числа обусловленности (возрастания
    %   обратной величины)
    cond_info = sortrows (cond_info, 2);
    %   Нахождение обратного преобразования индексов сортировки для
    %   определения ранга получающегося результата (то есть индекса в ряду
    %   значений от наихудшего к наилучшему)
    indices = [(1 : M + L)' cond_info(:, 1)];
    indices = sortrows (indices, 2);
    %   Вычисление гистограммы
    m = indices (row_best_hypothetical);
    hist_data (m, 1) = 1 + hist_data (m, 1);
    m = indices (row_wrst_hypothetical);
    hist_data (m, 2) = 1 + hist_data (m, 2);
end
%disp ([row_best_hypothetical row_wrst_hypothetical])
%disp (indices ([row_best_hypothetical row_wrst_hypothetical])')
%disp ([cond_info diag(s) lin_comb sqrt(sum (a .^ 2, 2))])

%   Отображение гистограммы
bar (100 * hist_data / num)
grid on
title ({'Ранжирование для "наихудшего" и "наилучшего"', 'выборов удаляемой строки матрицы'})
xlabel ({'Ранг получаемого результата', '(чем больше номер, тем ближе число обусловленности к 1)'})
ylabel ('Доля, %')
legend ('"Хороший" выбор', '"Плохой" выбор')
 

Понятно, что численный эксперимент — это не доказательство, что будет происходить в произвольном случае — ещё вопрос, а использованные мной в эксперименте матрицы имеют весьма специальный вид, и это будет влиять на результат. Но закономерность на лицо, поэтому хотелось бы знать: а существует ли какая-нибудь теория, объясняющая, откуда эта закономерность берётся?

[Geen]

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

Есть же такое понятие, как псевдообратная матрица
.

И какое это имеет отношение к тому хэлпу, который Вы цитировали?

[Евгений Машеров]

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

[zgemm]

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

zgemm, спасибо. А какой быстрый метод факторизации вы посоветуете, если мне нужно постоянно пополнять матрицу строками и удалять их для поддержания низкого числа обусловленности? Ник в своей статье обсудил три с разновидностями, но он делал упор на корректность и точность, а не на быстродействие. У меня, правда, миллионы строк не планируются, максимум сотни.

Тогда точно QR. Так называемый RRQR (rank-revealing QR) есть в lapack, который как я понимаю, можно из любых питонов, матлабов и октав вызывать, ну или нативно из С/С++/Fortran.

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

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