Возникла необходимость решения матричного уравнения следующего не вполне обычного вида:
\sum_{i} A_i*X*A_i^T = \sum_{i} B_i
Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что
1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь...
2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься.
Может быть, известны какие-то относительно простые подходы, оперирующие только матрицами NxN? Ткните в литературу...
Здравствуйте, kfmn, Вы писали:
K>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что K>1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь... K>2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься.
И в чем проблемма тут и есть N*N уравнений
Здравствуйте, kfmn, Вы писали:
K>Всем привет!
K>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что K>1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь... K>2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься.
Тут только путь 2 имхо. В X N^2 элементов. Вряд ли что-то будет проще, чем система линейных уравнений для всех них, коэфф-ты которой будут суммами по i попарных произведений эл-тов A_i.
Здравствуйте, kfmn, Вы писали:
K>Всем привет!
K>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что K>1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь... K>2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься.
K>Может быть, известны какие-то относительно простые подходы, оперирующие только матрицами NxN? Ткните в литературу...
У вас при n=2 будет уравнение вида:
A_1*X*A_1^T + A_2*X*A_2^T = B_1 + B_2,
где A_1, A_2, B_1, B_2, X — матрицы NxN
Вы это хотели сказать?
Здравствуйте, kov_serg, Вы писали:
_>Здравствуйте, kfmn, Вы писали:
K>>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>>Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что K>>1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь... K>>2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься. _>И в чем проблемма тут и есть N*N уравнений
Ну,я просто подумал, что стандартный метод Гаусса для решения уравнения с матрицей NxN потребует порядка N^3 операций. Соответственно, в моем случае это выльется уже в N^6.
А вдруг есть способы, скажем, свести это к набору из N (или даже N^2) уравнений с матрицами NxN относительно неких векторов, которые как-то связаны с искомой матрицей X, и восстановить ее по ним... Или еще что-нибудь, чтобы от 6-й степени перейти куда-то пониже...
Здравствуйте, kfmn, Вы писали:
K>Здравствуйте, Alex.a777, Вы писали:
AA>>У вас при n=2 будет уравнение вида: AA>>A_1*X*A_1^T + A_2*X*A_2^T = B_1 + B_2, AA>>где A_1, A_2, B_1, B_2, X — матрицы NxN AA>>Вы это хотели сказать?
K>Да, все матрицы — размерности NxN, а слагаемых в сумме — некое фиксированное количество, которое с N никак не связано, может быть и 2.
1) Ну метод Гаусса решает уравнение вида AX=B. Как планируете его применять к Вашей системе? Я так понимаю, итерационный подход поверх метода решения СЛАУ?
2) Касательно выбора методов, смотрите https://www.mathworks.com/help/matlab/ref/mldivide_full.png
Здравствуйте, Alex.a777, Вы писали:
AA>Здравствуйте, kfmn, Вы писали:
K>>Здравствуйте, Alex.a777, Вы писали:
AA>>>У вас при n=2 будет уравнение вида: AA>>>A_1*X*A_1^T + A_2*X*A_2^T = B_1 + B_2, AA>>>где A_1, A_2, B_1, B_2, X — матрицы NxN AA>>>Вы это хотели сказать?
K>>Да, все матрицы — размерности NxN, а слагаемых в сумме — некое фиксированное количество, которое с N никак не связано, может быть и 2.
AA>1) Ну метод Гаусса решает уравнение вида AX=B. Как планируете его применять к Вашей системе? Я так понимаю, итерационный подход поверх метода решения СЛАУ? AA>2) Касательно выбора методов, смотрите https://www.mathworks.com/help/matlab/ref/mldivide_full.png
Возможно, я недостаточно ясно выразился про решение с помощью вытягивания в вектор. Я имел в виду следующее: если расписать систему покомпонентно (в скалярной форме), то каждое из уравнений будет линейно относительно элементов матрицы X. Т.е. если все столбцы матрицы X составить в один длинный вектор Y размерностью N^2, то получим обычную СЛАУ типа CY=D, но с матрицей N^2xN^2. И решать ее можно, но хочется чего-то более быстрого.
Здравствуйте, kfmn, Вы писали:
K>Здравствуйте, Alex.a777, Вы писали:
AA>>Здравствуйте, kfmn, Вы писали:
K>>>Здравствуйте, Alex.a777, Вы писали:
AA>>>>У вас при n=2 будет уравнение вида: AA>>>>A_1*X*A_1^T + A_2*X*A_2^T = B_1 + B_2, AA>>>>где A_1, A_2, B_1, B_2, X — матрицы NxN AA>>>>Вы это хотели сказать?
K>>>Да, все матрицы — размерности NxN, а слагаемых в сумме — некое фиксированное количество, которое с N никак не связано, может быть и 2.
AA>>1) Ну метод Гаусса решает уравнение вида AX=B. Как планируете его применять к Вашей системе? Я так понимаю, итерационный подход поверх метода решения СЛАУ? AA>>2) Касательно выбора методов, смотрите https://www.mathworks.com/help/matlab/ref/mldivide_full.png
K>Возможно, я недостаточно ясно выразился про решение с помощью вытягивания в вектор. Я имел в виду следующее: если расписать систему покомпонентно (в скалярной форме), то каждое из уравнений будет линейно относительно элементов матрицы X. Т.е. если все столбцы матрицы X составить в один длинный вектор Y размерностью N^2, то получим обычную СЛАУ типа CY=D, но с матрицей N^2xN^2. И решать ее можно, но хочется чего-то более быстрого.
Вот Вы себе придумали работу! Не понял первоначального посыла, да))
Не претендую на оптимальность, но есть итерационный подход:
0) выбираем X0
1) решаем СЛАУ относительно X_i+1: A_1*X_i+1 = [summ(B_i,i=1,n) — summ(A_i*X_i*A_i^T, i=2,n)] * (A_1^T)^-1
2) пока ||X_i+1 — X_i||>eps X_i = X_i+1 и на шаг 1)
Здравствуйте, Alex.a777, Вы писали:
AA>Здравствуйте, kfmn, Вы писали:
K>>Здравствуйте, Alex.a777, Вы писали:
AA>>>Здравствуйте, kfmn, Вы писали:
K>>>>Здравствуйте, Alex.a777, Вы писали:
AA>>>>>У вас при n=2 будет уравнение вида: AA>>>>>A_1*X*A_1^T + A_2*X*A_2^T = B_1 + B_2, AA>>>>>где A_1, A_2, B_1, B_2, X — матрицы NxN AA>>>>>Вы это хотели сказать?
K>>>>Да, все матрицы — размерности NxN, а слагаемых в сумме — некое фиксированное количество, которое с N никак не связано, может быть и 2.
AA>>>1) Ну метод Гаусса решает уравнение вида AX=B. Как планируете его применять к Вашей системе? Я так понимаю, итерационный подход поверх метода решения СЛАУ? AA>>>2) Касательно выбора методов, смотрите https://www.mathworks.com/help/matlab/ref/mldivide_full.png
K>>Возможно, я недостаточно ясно выразился про решение с помощью вытягивания в вектор. Я имел в виду следующее: если расписать систему покомпонентно (в скалярной форме), то каждое из уравнений будет линейно относительно элементов матрицы X. Т.е. если все столбцы матрицы X составить в один длинный вектор Y размерностью N^2, то получим обычную СЛАУ типа CY=D, но с матрицей N^2xN^2. И решать ее можно, но хочется чего-то более быстрого.
AA>Вот Вы себе придумали работу! Не понял первоначального посыла, да)) AA>Не претендую на оптимальность, но есть итерационный подход: AA>0) выбираем X0 AA>1) решаем СЛАУ относительно X_i+1: A_1*X_i+1 = [summ(B_i,i=1,n) — summ(A_i*X_i*A_i^T, i=2,n)] * (A_1^T)^-1 AA>2) пока ||X_i+1 — X_i||>eps X_i = X_i+1 и на шаг 1)
Спасибо за идею, но, во-первых, матрицы A_i могут не быть обратимы, а во-вторых, будет ли такой процесс сходиться это большой вопрос...
Здравствуйте, Слава, Вы писали:
С>Здравствуйте, kfmn, Вы писали:
K>>Может быть, известны какие-то относительно простые подходы, оперирующие только матрицами NxN? Ткните в литературу...
С>Я бы предложил подход с другой стороны — попытаться использовать GPU для самим операций над матрицами. Машина железная, пусть считает.
Грубая сила, ага.. Но если с матрицами 20x20 это прокатит, решать систему 400x400 не слишком накладно, то вот с 1000x1000 (исходными) — уже нет. Никакая современная GPU матрицу 10^6x10^6 не сожрет.
У меня, конечно, матрицы не такие огромные, скорее, ближе к 1-му примеру, но интересно алгоритмическое ускорение.
Здравствуйте, kfmn, Вы писали:
K>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i
Интересно. Если я ничего не напутал A * X * A^T = (A * X) * A^T = (X^T * A^T)^T * A^T = A^T^T * (X^T * A^T)^T^T = A * X^T * A^T
Т.е. всегда два решения с точностью до транспонирования — если решение верно для X, то оно будет верно и для X^T
Не знаю правда чем это тебе поможет
но это не зря, хотя, может быть, невзначай
гÅрмония мира не знает границ — сейчас мы будем пить чай
Здравствуйте, ·, Вы писали:
·>Здравствуйте, kfmn, Вы писали:
K>>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i ·>Интересно. Если я ничего не напутал ·>A * X * A^T = (A * X) * A^T = (X^T * A^T)^T * A^T ====== A^T^T * (X^T * A^T)^T^T = A * X^T * A^T ·>Т.е. всегда два решения с точностью до транспонирования — если решение верно для X, то оно будет верно и для X^T ·>Не знаю правда чем это тебе поможет
Напутал: равенство, выделенное длинным. Ты там хотел, видимо, еще раз транспонировать, а внешнее транспонирование (обратно) потерял...
Хотя искомая матрица X в самом деле симметричная (по смыслу задачи). Это, вообще говоря, сокращает число скалярных уравнений почти вдвое, но их все равно много
Здравствуйте, kfmn, Вы писали:
K>>>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>>>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>·>Интересно. Если я ничего не напутал K>·>A * X * A^T = (A * X) * A^T = (X^T * A^T)^T * A^T ====== A^T^T * (X^T * A^T)^T^T = A * X^T * A^T K>·>Т.е. всегда два решения с точностью до транспонирования — если решение верно для X, то оно будет верно и для X^T K>·>Не знаю правда чем это тебе поможет
K>Напутал: равенство, выделенное длинным. Ты там хотел, видимо, еще раз транспонировать, а внешнее транспонирование (обратно) потерял...
Ах, да. Получается A * X * A^T = (A * X^T * A^T)^T, не очень-то полезно...
но это не зря, хотя, может быть, невзначай
гÅрмония мира не знает границ — сейчас мы будем пить чай
Здравствуйте, kfmn, Вы писали:
AA>>Вот Вы себе придумали работу! Не понял первоначального посыла, да)) AA>>Не претендую на оптимальность, но есть итерационный подход: AA>>0) выбираем X0 AA>>1) решаем СЛАУ относительно X_i+1: A_1*X_i+1 = [summ(B_i,i=1,n) — summ(A_i*X_i*A_i^T, i=2,n)] * (A_1^T)^-1 AA>>2) пока ||X_i+1 — X_i||>eps X_i = X_i+1 и на шаг 1)
K>Спасибо за идею, но, во-первых, матрицы A_i могут не быть обратимы, а во-вторых, будет ли такой процесс сходиться это большой вопрос...
1) Достаточно обратимости любой из матриц A_i^T. В качестве опорной же можно выбрать любую. Если менять опорную матрицу, то не факт что сходимость улучшиться, а ограничений больше, да
2) Сходимость, подозреваю, будет определяться условием типа критерия Скарбороу, т.е. в практическом смысле бесполезным
3) Как бы то ни было, самые самые солверы для плотных хорошо обусловленных матриц — это или многосеточные, или семейство итерационных методов Крылова
Здравствуйте, kfmn, Вы писали:
K>Здравствуйте, kov_serg, Вы писали:
_>>Здравствуйте, kfmn, Вы писали:
K>>>Возникла необходимость решения матричного уравнения следующего не вполне обычного вида: K>>>\sum_{i} A_i*X*A_i^T = \sum_{i} B_i K>>>Матрицы A_i, B_i — известные, X — искомая, все квадратные NxN. Понятно, что K>>>1) сумму B_i можно посчитать заранее, но, возможно, эта исходная форма в виде суммы способна чему-то помочь... K>>>2) матрицу X можно вытянуть в вектор, и тогда получается одна длинная линейная система с размером N^2. Но ее решать, соответственно, намного дольше, да и составлять заколебешься. _>>И в чем проблемма тут и есть N*N уравнений
K>Ну,я просто подумал, что стандартный метод Гаусса для решения уравнения с матрицей NxN потребует порядка N^3 операций. Соответственно, в моем случае это выльется уже в N^6.
Не используй метод гаусс. В большинстве случае сопряженные градиенты (например https://pdfs.semanticscholar.org/466d/addfb6340c28cb8da548007028c8cc5df687.pdf) сойдутся гораздо быстрее. K>А вдруг есть способы, скажем, свести это к набору из N (или даже N^2) уравнений с матрицами NxN относительно неких векторов, которые как-то связаны с искомой матрицей X, и восстановить ее по ним... Или еще что-нибудь, чтобы от 6-й степени перейти куда-то пониже...
Если норма одной матрицы сильно больше других можно использовать метод возмущений.
Здравствуйте, kfmn, Вы писали:
K>Грубая сила, ага.. Но если с матрицами 20x20 это прокатит, решать систему 400x400 не слишком накладно, то вот с 1000x1000 (исходными) — уже нет. Никакая современная GPU матрицу 10^6x10^6 не сожрет.
А CPU сожрёт что ли ? 1млн*1млн*8 это 8 Тб данных, в двойной точности.
1000*1000 это всего лишь 4 мегабайта float-ов или 8mb даблов.
Современные GPU имеют от 8 GB на борту. Ты туда своих матриц 1000*1000 можешь напихать по-самые помидоры 1000 матриц с даблами, размерностью 1000*1000 каждая. Мало что ли ?
Предельным размером видится что-то порядка 22000 * 22000 @ Single Precision для 8GB карты. Влезет 4 матрицы.
Только учти даблы существенно медленее в современных не-проф картах.
Например моя 1080ti выдаёт:
— 12869 GFlops в Single Precision
— 417 GFlops (всего лишь) в Double Precision.
Древний Radeon HD7990 не обрезали в двойной точности, и у него завидные 1821 GFlops, хотя в одинарной не впечатляющие 7296 GFlops. Правда памяти у него всего по 3GB на чип.
Для сравнения: 4 ядерный 8 поточный i7-7700K на 4.5ггц имеет:
— 460 GFlops одинарной.
— 230 GFlops двойной точности
GPU почти в 30 раз быстрее в Single Precision вычислениях.
Кроме того на таких объёмах наверняка будет играть роль скорость памяти. У GPU память быстрее раз в 10.
1000 матриц с даблами, размерностью 1000*1000 каждая. Мало что ли ? 
