Résolution

22

J'ai des matrices et . est clairsemé et est avec très grand (peut être de l'ordre de plusieurs millions.) est une matrice haut avec plutôt petit ( ) et chaque colonne peut seulement un seul entrée avec le reste étant « s, de telle sorte que . est énorme, il est donc très difficile à inverser, et je peux résoudre un système linéaire tel que manière itérative en utilisant une méthode de sous-espace de Krylov telle que , mais je n'ai pasAGAn×nnGn×mm1<m<100010GTG=IAAx=bBiCGStab(l)A1 explicitement.

Je veux résoudre un système de la forme: , où et sont des vecteurs de longueur . Une façon de le faire est d'utiliser un algorithme itératif dans un algorithme itératif pour résoudre pour chaque itération de l'algorithme itératif externe. Ce serait cependant extrêmement coûteux en calcul. Je me demandais s'il y avait un moyen plus simple sur le plan informatique de résoudre ce problème.(GTA1G)x=bxbmA1

Costis
la source
Je viens d'ajouter à ma réponse une remarque sur l'exploitation de la structure 0-1.
Arnold Neumaier

Réponses:

19

Introduisez le vecteur et résolvez le grand système couplé A y + G x = 0 , G T y = - b pour ( y , x ) simultanément, en utilisant une méthode itérative. Si A est symétrique (comme il semble probable que vous ne le déclariez pas explicitement), alors le système est symétrique (mais indéfini, bien que quasi-défini si Ay:=A1GxAy+Gx=0GTy=b(y,x)AAest définie positive), ce qui pourrait vous aider à choisir une méthode appropriée. (mots clés pertinents: matrice KKT, matrice quasi-définie).

Edit: Comme est symétrique complexe, la matrice augmentée l'est aussi, mais il n'y a pas de quasi-infinité. Vous pouvez cependant utiliser la routine A x pour calculer A x = ¯ A ¯ x ; vous pouvez donc adapter une méthode telle que QMR ftp://ftp.math.ucla.edu/pub/camreport/cam92-19.pdf (conçue pour les systèmes réels, mais vous pouvez facilement la réécrire pour les systèmes complexes, en utilisant l'adjoint dans lieu de la transposition) pour résoudre votre problème.AAxAx=Ax¯¯

Edit2: En fait, la structure (0,1) de signifie que vous pouvez éliminer x et les composants de G T y symboliquement, se retrouvant ainsi avec un système plus petit à résoudre. Cela signifie jouer avec la structure de A , et ne paie que lorsque A est donné explicitement dans un format clairsemé plutôt que comme un opérateur linéaire.GxGTyAA

Arnold Neumaier
la source
Merci! A est symétrique complexe. Y a-t-il des raisons de s'attendre à ce que l'état de la matrice augmentée soit pire que celui de la matrice ? Si m est petit, la matrice augmentée n'est que légèrement plus grande que A, donc je soupçonne que la résolution itérative de ce système augmenté ne devrait pas être beaucoup plus difficile que la résolution d'un système avec A? A
Costis
Le numéro de condition des deux systèmes est généralement assez indépendant; cela dépend beaucoup de ce que est. - J'ai ajouté à ma réponse des informations sur la façon d'exploiter la symétrie complexe. G
Arnold Neumaier
Salut les gars! Merci pour toutes vos réponses; cet endroit est super! Une extension de la question d'origine: supposons maintenant que j'ai , où G et A ont la même signification que dans la question d'origine mais B est une matrice nxn déficiente en rang ( même taille que A) et l'ensemble G T A - H B A - 1 G est de rang complet. Comment allez-vous résoudre le nouveau système, puisque maintenant l'inverse de B n'existe pas, vous ne pouvez donc pas avoir A B - 1 A H(GTAHBA1G)x=bGTAHBA1GAB1AH. Je ne pense pas non plus que cela fonctionnerait simplement avec le pseudoinverse de B.
Costis
1
Introduisez et z : = A - H B y , et procédez par analogie avec le cas élaboré. (Probablement, vous devez également factoriser B dans les matrices de rang complet et introduire un vecteur intermédiaire supplémentaire.)y:=A1Gxz:=AHByB
Arnold Neumaier
Salut Arnold. Merci, cela fonctionne vraiment! Je l'ai testé avec de très petits exemples de test, et cela fonctionne très bien. Cependant, mon solveur itératif a d'énormes problèmes pour inverser la matrice augmentée. Alors qu'il ne faut que 80 itérations environ (quelques secondes) pour résoudre un système de la forme avec la matrice A d'origine, le système avec la matrice augmentée (qui est 2n + mx 2n + m ou 2n-mx 2n- m en utilisant l'approche de @ wolfgang-bangerth) prend des dizaines de milliers d'itérations (plusieurs heures) à résoudre pour un RHS. Existe-t-il des stratégies pour accélérer la convergence? peut-être un préconditionneur? Ax=b
Costis
7

Après la réponse d'Arnold, vous pouvez faire quelque chose pour simplifier le problème. Plus précisément, réécrivez le système sous la forme . Notez ensuite que d'après l'énoncé que G est grand et étroit et que chaque ligne n'a qu'un 1 et des zéros sinon, l'énoncé G T y = - b signifie qu'un sous-ensemble des éléments de y a une valeur fixe, à savoir les éléments de - b .Ay+Gx=0,GTy=bGGTy=byb

Disons que pour simplifier que a m colonnes et n lignes et qu'exactement les m premières lignes en contiennent et que réorganiser les éléments de x je puisse faire en sorte que G ait la matrice d'identité m × m en haut et une matrice n - m × m zéro en bas. Ensuite, je peux partitionner y = ( y c , y f ) en m éléments "contraints" et n - m "libres" afin queGmnmxGm×mnm×my=(yc,yf)mnm . Je peux également partitionner A pour que A = ( A c c A c f A f c A f f ) . De l'équation A y + G x = 0, j'obtiens alors: A c c y c + A c f y f + x = 0 ,yc=bAA=(AccAcfAfcAff)Ay+Gx=0 et en utilisant ce que nous savons sur y c, nous avons de la seconde de ces équations A f f y f = A f c b et par conséquent x = A c c b - A c f A - 1 f f A f c b . En d'autres termes, la seule matrice que vous devez inverser est le sous-ensemble de A

Accyc+Acfyf+x=0,Afcyc+Affyf=0
yc
Affyf=Afcb
x=AccbAcfAff1Afcb.
Adont les lignes et les colonnes ne sont pas mentionnées dans (l'espace nul de G ). Vous pouvez facilement le faire: (i) calculer z = A f c b ; (ii) utilisez le solveur que vous avez pour résoudre A f f h = z ; (iii) calculer x = A c c b - A c f h .GGz=AfcbAffh=zx=AccbAcfh

En d' autres termes, étant donné la structure de , la résolution du système linéaire que vous avez est vraiment pas plus difficile que la résolution d' un seul système linéaire avec A .GA

Wolfgang Bangerth
la source
0

Mais nous connaissons , G T et A , doncGGTA

GTA1Gx=b

GGTA1Gx=Gb

Puisque , alors G T = G - 1 , donc G G T = I :GTG=IGT=G1GGT=I

A1Gx=Gb

AA1Gx=AGb

Gx=AGb

GTGx=GTAGb

x=GTAGb

Sauf si j'ai raté quelque chose, vous n'avez besoin d'aucune itération, ni d'aucun solveur pour calculer x étant donné , A et b .GAb

Phil H
la source
3
étant un inverse gauche de G n'implique pas qu'il soit également un inverse droit. Considérons G = e 1 , où G T = e T 1 est un inversegauche, mais G G T = e 1 e T 1I . GTGG=e1GT=e1TGGT=e1e1TI
Jack Poulson
1
, donc G T G = I m × m , mais G G TI n × n . C'est plutôt un projecteur sur un sous-espace. GCnCmGTG=Im×mGGTIn×n
Deathbreath