Решение

22

У меня есть матрицы A и G . A является разреженным и имеет размер n × n с очень большим n (может быть порядка нескольких миллионов). G является матрицей высотой n × m с довольно небольшим m ( 1 < m < 1000 ), и в каждом столбце может быть только один 1 запись с остальным 0 «с, таким образом, что G T G = Я . А огромен, так что это действительно трудно инвертировать, и я могу решить линейную систему, такую ​​как хAGAn×nnGn×mm1<m<100010GTG=IA A = b итеративно, используя метод подпространства Крылова, такой как B i C G S t a b ( l ) , но у меня нет A - 1 явно.Ax=bBiCGStab(l)A1

Я хочу решить систему вида: ( G T A - 1 G ) x = b , где x и b - векторы длины m . Один из способов сделать это - использовать итерационный алгоритм в итерационном алгоритме, чтобы найти для A - 1 для каждой итерации внешнего итерационного алгоритма. Это было бы чрезвычайно дорого в вычислительном отношении, однако. Мне было интересно, есть ли в вычислительном отношении более простой способ решения этой проблемы.(GTA1G)x=bxbmA1

КОСТИС
источник
Я только добавил к своему ответу замечание об использовании структуры 0-1.
Арнольд Ноймайер

Ответы:

19

Введем вектор y : = - A - 1 G x и решим большую связанную систему A y + G x = 0 , G T y = - b для ( y , x ) одновременно, используя итерационный метод. Если A симметричен (как представляется вероятным, хотя вы не указали это явно), то система является симметричной (но неопределенной, хотя и квазиопределенной, если Ay:=A1GxAy+Gx=0GTy=b(y,x)AAположительно определен), что может помочь вам выбрать подходящий метод. (релевантные ключевые слова: матрица ККТ, квазиопределенная матрица).

Редактировать: Как A является сложной симметричной, так и расширенная матрица, но нет квазиопределенности. Однако вы можете использовать й подпрограмму для вычисления А * х = ¯ ¯ х ; поэтому вы можете адаптировать метод, такой как QMR, ftp://ftp.math.ucla.edu/pub/camreport/cam92-19.pdf (разработанный для реальных систем, но вы можете легко переписать его для сложных систем, используя сопряженный в место транспонирования), чтобы решить вашу проблему.AAxAx=Ax¯¯¯¯¯¯¯¯¯¯

Edit2: На самом деле, (0,1) -структура G означает, что вы можете удалить x amd и компоненты G T y символически, таким образом получая меньшую систему для решения. Это значит возиться со структурой A и оплачивать только тогда, когда A задано явно в разреженном формате, а не как линейный оператор.GxGTyAA

Арнольд Ноймайер
источник
Спасибо! А является комплексным симметричным. Есть ли основания ожидать, что состояние расширенной матрицы будет хуже, чем состояние исходной матрицы A ? Если m мало, расширенная матрица лишь незначительно больше по размеру, чем A, поэтому я подозреваю, что итеративное решение этой расширенной системы не должно быть намного сложнее, чем решение системы с A?A
Костис
Число условий двух систем, как правило, совершенно не связано; это очень сильно зависит от того, что G есть. - Я добавил в свой ответ информацию о том, как использовать сложную симметрию. G
Арнольд Ноймайер
Привет, народ! Спасибо за все отклики; это место отличное! Дополнение к исходному вопросу. Предположим теперь, что у меня ( G T A - H B A - 1 G ) x = b , где G и A имеют то же значение, что и в исходном вопросе, но B является матрицей nxn с недостатком ранга ( того же размера, что и A), и весь G T A - H B A - 1 G имеет полный ранг. Как бы вы пошли о решении новой системы, поскольку теперь обратное значение B не существует, поэтому вы не можете иметь A B - 1 A H(GTAHBA1G)x=bGTAHBA1GAB1AH . Я не думаю, что это будет работать просто с псевдообратным B тоже.
Костис
1
Введите y : = A - 1 G x и z : = A - H B y и действуйте аналогично разработанному случаю. (Возможно, вам также нужно разложить B на матрицы полного ранга и ввести дополнительный промежуточный вектор.)y:=A1Gxz:=AHByB
Арнольд Ноймайер
Привет, Арнольд. Спасибо, это действительно работает! Я проверил его на нескольких очень маленьких примерах, и он отлично работает. Однако, у моего итеративного решателя есть огромные проблемы, инвертирующие расширенную матрицу. В то время как для решения системы вида A x = b с исходной матрицей A требуется всего около 80 итераций (несколько секунд) , система с расширенной матрицей (которая равна 2n + mx 2n + m или 2n-mx 2n- м с использованием подхода @ wolfgang-bangerth) для решения одной RHS требуется более десятков тысяч итераций (несколько часов). Существуют ли стратегии для ускорения конвергенции? возможно предварительный кондиционер? Ax=b
Костис
7

Following Arnold's reply, there is something you can do to simplify the problem. Specifically, rewrite the system as Ay+Gx=0,GTy=bAy+Gx=0,GTy=b. Then note that from the statement that GG is tall and narrow and each row has only one 1 and zeros otherwise, then the statement GTy=bGTy=b means that a subset of the elements of yy have a fixed value, namely the elements of bb.

Let us say that for simplicity that GG has mm columns and nn rows and that exactly the first mm rows have ones in them and that be reordering the elements of xx I can make it so that GG has the m×mm×m identity matrix at the top and a nm×mnm×m zero matrix at the bottom. Then I can partition y=(yc,yf)y=(yc,yf) into mm "constrained" and nmnm "free" elements so that yc=byc=b. I can also partition AA so that A=(AccAcfAfcAff)A=(AccAfcAcfAff). From the equation Ay+Gx=0Ay+Gx=0 I then get the following: Accyc+Acfyf+x=0,Afcyc+Affyf=0

Accyc+Acfyf+x=0,Afcyc+Affyf=0
and using what we know about ycyc we have from the second of these equations Affyf=Afcb
Affyf=Afcb
and consequently x=AccbAcfA1ffAfcb.
x=AccbAcfA1ffAfcb.
In other words, the only matrix you have to invert is the subset of AA whose rows and columns are not mentioned in GG (the null space of GG). This you can easily do: (i) compute z=Afcbz=Afcb; (ii) use whatever solver you have to solve Affh=zAffh=z; (iii) compute x=AccbAcfhx=AccbAcfh.

In other words, given the structure of GG, solving the linear system you have is really not more difficult than solving a single linear system with AA.

Wolfgang Bangerth
источник
0

But we know GG, GTGT and AA, so

GTA1Gx=bGTA1Gx=b

GGTA1Gx=GbGGTA1Gx=Gb

Since GTG=IGTG=I, then GT=G1GT=G1, so GGT=IGGT=I:

A1Gx=GbA1Gx=Gb

AA1Gx=AGbAA1Gx=AGb

Gx=AGbGx=AGb

GTGx=GTAGbGTGx=GTAGb

x=GTAGbx=GTAGb

Unless I've missed something, you don't need any iteration, or any solver to calculate x given GG, AA and bb.

Phil H
источник
3
GTGT being a left inverse of GG does not imply that it is also a right inverse. Consider G=e1G=e1, where GT=eT1GT=eT1 is a left inverse, but GGT=e1eT1IGGT=e1eT1I.
Jack Poulson
1
GCnCm, hence GTG=Im×m, but GGTIn×n. Rather it's a projector on a subspace.
Deathbreath