Симметричная матрица -> внешнее произведение

[madschumacher]

Добрый день.

Очередные глупые вопросы лезут в голову. У меня имеется симметричная действительная матрица квадратная $\mathcal{A}$ размера $N \times N$. Допустим, я хочу найти наилучшее её представление в виде внешнего произведения вектора $\mathbf{u}$ длины $N$, скажем в смысле
$\sum_{i,j = 1}^{N} \left( \underbrace{(\mathcal{A} - \mathbf{u} \mathbf{u}^\dagger)_{ij}}_{\mathcal{A}_{ij} - u_i u_j} \right)^2 \rightarrow \min$.
Есть ли какие-то библиотечные алгоритмы, которые позволяют такое делать?

Судя по поиску, если бы меня интересовало представление в виде внешнего произведения двух векторов, то ответом было бы SVD.

В этом же случае я только смог свести это всё к уравнению
$\mathcal{A} \mathbf{u} = \underbrace{\mathcal{U}}_{\mathbf{u} \mathbf{u}^\dagger} \mathbf{u}$
и для его решения я могу представить себе итерационные методы, типа
$\mathcal{A} \mathbf{u}_{n+1} = \mathcal{U}_{n} \mathbf{u}_{n}$
где нижний индекс нумерует итерации.

[amon]

madschumacher в сообщении #1441418 писал(а):

$\sum_{i,j = 1}^{N} \left( \underbrace{(\mathcal{A} - \mathbf{u} \mathbf{u}^\dagger)_{ij}}_{\mathcal{A}_{ij} - u_i u_j} \right)^2 \rightarrow \min$.

А-ля метод мин. квадратов - не? Дифференцируем по $u_i,$ приравниваем производную к нулю, получаем $N$ уравнений (линейных). Только Вам, небось, нужно, что бы "мало" отличались не матрицы, а их собственные значения, а в таком случае эта норма не катит.

[madschumacher]

amon в сообщении #1441424 писал(а):

А-ля метод мин. квадратов - не?

Так там и получается нелинейное уравнение в итоге:

madschumacher в сообщении #1441418 писал(а):

$\mathcal{A} \mathbf{u} = \underbrace{\mathcal{U}}_{\mathbf{u} \mathbf{u}^\dagger} \mathbf{u}$

amon в сообщении #1441424 писал(а):

Только Вам, небось, нужно, что бы "мало" отличались не матрицы, а их собственные значения, а в таком случае эта норма не катит.

Да нет, вроде именно что матрицы, я пытаюсь очистить изображение особого сорта от артефактов...

[amon]

madschumacher в сообщении #1441425 писал(а):

Так там и получается нелинейное уравнение в итоге

Да, чего-то я нынче совсем не в форме.

[Xaositect]

SVD от симметрической матрицы будет симметричным и называется оно обычно Eigenvalue decomposition ($A = UDU^T$, $U$ ортогональная). Потом берем вектор соотв. максимальному сз и делаем $u \lambda u^T = (\sqrt{\lambda}u)(\sqrt{\lambda}u)^T$. LAPACK умеет.

Что, собственно, получается и у Вас: $Au = uu^T u = u \|u\|^2$, т.е. $u$-собственный вектор.

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

По-моему, Вам нужен собственный вектор, соответствующий максимальному собственному значению. Любой алгоритм, начиная со степенного метода, которым можно даже вручную побаловаться...

[arseniiv]

madschumacher
Кстати outer product обычно не переводится как «внешнее произведение», внешнее чаще соответствует exterior product.

[madschumacher]

Спасибо большое, Xaositect, Евгений Машеров, действительно работает. Для моих целей оказалось, что даже можно несколько собственных векторов брать.

arseniiv, спасибо, буду иметь в виду. А как тогда outer product будет по-русски корректно называться?

[arseniiv]

А вот тут я поискал и не нашёл (или начисто забыл). Можно попытаться говорить «тензорное произведение», в конце концов это частный случай именно его ($u\otimes v^\flat$, где $v^\flat$ — опускание (с помощью скалярного произведения) индекса у вектора, превращающее его в 1-форму — в индексной записи это всё будет $u^i v_j$, ну и кстати как вы наверно в курсе, в бра-кет-обозначениях то же самое есть $\lvert u\rangle\langle v\rvert$, а такие произведения «без упрощения» (когда не умножается непосредственно бра на кет) как раз тензорные). В присутствии скалярного произведения все тензорные произведения двух множителей, «использующие» два вектора — кроме первого это $u\otimes v, u^\flat\otimes v, u^\flat\otimes v^\flat$ и ещё четыре с обратным порядком — будут приводимыми одно к другому). Может, кто-то и более точное название слышал.

-- Вт фев 25, 2020 21:14:04 --

А, ну и да, вы же говорите о приближении оператора оператором ранга 1, и можно говорить «оператор ранга 1».

-- Вт фев 25, 2020 21:19:04 --

(Подобным образом вообще говорят о приближении произвольных тензоров тензорами ранга $m$ — только не того ранга, который связан с валентностью, а того, который равен наименьшему числу слагаемых вида $v_1\otimes\ldots\otimes v_n$, которые надо взять, чтобы в сумме получить такой тензор; для линейных операторов это совпадёт с их рангом.)

[ilghiz]

arseniiv в сообщении #1441508 писал(а):

А вот тут я поискал и не нашёл (или начисто забыл). Можно попытаться говорить «тензорное произведение», в конце концов это частный случай именно его ($u\otimes v^\flat$, где $v^\flat$ — опускание (с помощью скалярного произведения) индекса у вектора, превращающее его в 1-форму — в индексной записи это всё будет $u^i v_j$, ну и кстати как вы наверно в курсе, в бра-кет-обозначениях то же самое есть $\lvert u\rangle\langle v\rvert$, а такие произведения «без упрощения» (когда не умножается непосредственно бра на кет) как раз тензорные). В присутствии скалярного произведения все тензорные произведения двух множителей, «использующие» два вектора — кроме первого это $u\otimes v, u^\flat\otimes v, u^\flat\otimes v^\flat$ и ещё четыре с обратным порядком — будут приводимыми одно к другому). Может, кто-то и более точное название слышал.

Все тензорные алгоритмы все равно через SVD или задачи на собственные значения решаются, там по сути ничего не меняется.

2 TC: А матрица у Вас теоретическая, или практическая, и какие размеры ее Вас интересуют? От этого сильно метод решения может измениться. Так, например, для на попробовать с матрицей порядка $5 \times 5$, проще всего эту матрицу несколько раз в квадрад возвести и первый столбец взять, а в случае с $5 * 10^6 \times 5* 10^6 $ Ланцоша или Давидсона использовать, а между - может и лапаком.

[madschumacher]

ilghiz в сообщении #1441550 писал(а):

А матрица у Вас теоретическая, или практическая, и какие размеры ее Вас интересуют?

Порядка $10^3 \times 10^3$. Через Numpy все замечательно работает.

[arseniiv]

ilghiz в сообщении #1441550 писал(а):

Все тензорные алгоритмы все равно через SVD или задачи на собственные значения решаются, там по сути ничего не меняется.

Я вообще ничего не говорил о том как решать задачу, так как это действительно выше уже сказали, так что не понятно, к чему возражать. (Хотя кстати говоря SVD и eigendecomposition не определены для тензоров произвольного ранга вот прям так.) Отвечал лишь на вопрос, как называть выражения типа $\lvert u\rangle\langle v\rvert$.

[ilghiz]

Numpy - это фактически lapack. Для таких размерностей итерировать может быть дороже, чем в лоб через $QTQ^T, ~~Q^TQ=I$, и $T$ - трехдиагональная и далее QR алгоритмом, как собственно в лапаке и нумпи реализовано.

[madschumacher]

ilghiz в сообщении #1441555 писал(а):

Для таких размерностей итерировать может быть дороже, чем в лоб через $QTQ^T, ~~Q^TQ=I$, и $T$ - трехдиагональная и далее QR

Собственно, на решение уходит меньше минуты вычислительного времени (поиск собственных векторов - - не самая затратная операция кода). Так что оптимизировать код смысла мне пока нет.