Quand les transformations orthogonales surpassent-elles l'élimination gaussienne?

22

Comme nous le savons, les méthodes de transformations orthogonales (rotations de Givens et réflexions de Housholder) pour les systèmes d'équations linéaires sont plus coûteuses que l'élimination gaussienne, mais ont théoriquement de plus bonnes propriétés de stabilité dans le sens où elles ne modifient pas le numéro de condition du système. Bien que je ne connaisse qu'un exemple académique d'une matrice gâchée par l'élimination gaussienne avec pivotement partiel. Et il est communément admis qu'il est très peu probable de rencontrer ce type de comportement dans la pratique (voir ces notes de cours [pdf] ).

Alors, où chercher la réponse sur le sujet? Implémentations parallèles? Mise à jour? ..

linear-algebra reference-request faleichik
la source

24

Précision

Trefethen et Schreiber ont écrit un excellent article, Stabilité moyenne des cas d'élimination gaussienne , qui traite du côté exactitude de votre question. Voici quelques-unes de ses conclusions:

"Pour la factorisation QR avec ou sans pivot de colonne, l'élément maximal moyen de la matrice résiduelle est , tandis que pour l'élimination gaussienne, il est . Cette comparaison révèle que l'élimination gaussienne est légèrement instable , mais l'instabilité ne serait détectable que pour de très gros problèmes de matrice résolus en basse précision. Pour la plupart des problèmes pratiques, l'élimination gaussienne est très stable en moyenne. "(Emphasis mine) $O(n^{1/2})$ $O(n)$
"Après les premières étapes de l'élimination gaussienne, les éléments de matrice restants sont approximativement normalement distribués, qu'ils aient commencé de cette façon ou non."

Il y a beaucoup plus dans le document que je ne peux pas saisir ici, y compris la discussion de la matrice du pire cas que vous avez mentionnée, donc je vous recommande fortement de le lire.

Performance

Pour les matrices réelles carrées, la LU avec pivotement partiel nécessite environ flops, tandis que la QR basée sur les ménages nécessite environ flops. Ainsi, pour des matrices carrées raisonnablement grandes, la factorisation QR ne sera environ deux fois plus chère que la factorisation LU. $2/3 n^3$ $4/3 n^3$

Pour les matrices , où , LU avec pivotement partiel nécessite flops, contre QR (qui est toujours le double de celui de la factorisation LU) . Cependant , il est surprenant que les applications produisent des matrices maigres très hautes ( ), et Demmel et al. avoir un bon article, Factorisation QR parallèle et séquentielle évitant la communication , qui (dans la section 4) discute d'un algorithme intelligent qui n'exige que messages soient envoyés lorsque processeurs sont utilisés, par rapport aux messages des approches traditionnelles . La dépense est que $m \times n$ $m \ge n$ $mn^2 - n^3/3$ $2mn^2 - 2n^3/3$ $m \gg n$ $\log p$ $p$ $n \log p$ $O(n^3 \log p)$ des flops supplémentaires sont effectués, mais pour les très petits cela est souvent préféré au coût de latence de l'envoi de plus de messages (au moins lorsqu'une seule factorisation QR doit être effectuée). $n$

Jack Poulson
la source

10

Je suis surpris que personne n'ait mentionné les problèmes de moindres carrés linéaires , qui se produisent fréquemment dans le calcul scientifique. Si vous souhaitez utiliser l'élimination gaussienne, vous devez former et résoudre les équations normales, qui ressemblent à:

A^{T} A x = A^{T} b,

$A^{T}Ax = A^{T}b,$

où est une matrice de points de données correspondant aux observations de variables indépendantes, est un vecteur de paramètres à trouver et est un vecteur de points de données correspondant aux observations d'une variable dépendante. $A$ $x$ $b$

Comme Jack Poulson le souligne fréquemment, le nombre de conditions de est le carré du nombre de conditions de , de sorte que les équations normales peuvent être désastreusement mal conditionnées. Dans de tels cas, bien que les approches basées sur QR et SVD soient plus lentes, elles donnent des résultats beaucoup plus précis. $A^{T}A$ $A$

Geoff Oxberry
la source

2

Upvoted, mais QR devrait en fait être à égalité avec LU si l' on considère l'inutile

opérations nécessaires pour former

(QR ne nécessite que

plus de flops que LU). L'approche SVD devrait cependant être plus lente (on peut penser à son coût à environ

).

n^{3}

$n^3$

A^{H} A

$A^H A$

2 / 3 n^{3}

$2/3 n^3$

6 n^{3}

$6n^3$

Jack Poulson

1

En plus de la stabilité garantie par l'utilisation de transformations orthogonales, le grand avantage de SVD est que la décomposition fournit son propre contrôle de condition, car le rapport de la plus grande à la plus petite valeur singulière est précisément le nombre de conditions (2 normes). Pour les autres décompositions, l'utilisation d'un estimateur de conditions (par exemple Hager-Higham) est, bien que moins coûteuse que la décomposition proprement dite, quelque peu «clouée».

JM

1

@JackPoulson Juste par curiosité, avez-vous une référence pour votre compte de flop pour SVD? D'après ce que je peux dire d'un rapide coup d'œil dans Golub & Van Loan (p. 254 3e édition), la constante semble plus élevée pour utiliser le SVD dans la résolution des problèmes des moindres carrés, mais je peux me tromper. Merci d'avance.

OscarB

1

@OscarB: C'était un nombre très approximatif du haut de ma tête, ce qui est inférieur à la formation de la SVD complète (car nous pouvons éviter les coûts de retransformation).

de travail est nécessaire pour la réduction de la forme bidiagonale ( par exemple,

), une certaine quantité de travail, dit

, est nécessaire pour la SVD bidiagonale (

), et ensuite

8 / 3 n^{3}

$8/3 n^3$

A = F B G^{H}

$A=FBG^H$

C

$C$

B = U Σ V^{H}

$B=U\Sigma V^H$

, ce qui devrait nécessiter un travail

. Ainsi, tout dépend de la taille du

... si MRRR fonctionne jamais ici, ce sera

, mais jusque-là, il est cubique et dépend du problème.

x := (G (V (i n v (Σ) (U^{H} (F^{H} b)))))

$x := (G (V (\mathrm{inv}(\Sigma) (U^H (F^H b)))))$

O (n^{2})

$O(n^2)$

C

$C$

O (n^{2})

$O(n^2)$

Jack Poulson

1

@JM Notez cependant que le numéro de condition du problème des moindres carrés n'est pas le numéro de condition "classique"

d'une matrice; c'est une quantité plus compliquée.

\frac{σ_{1}}{σ_{n}}

$\frac{\sigma_1}{\sigma_n}$

Federico Poloni

3

Comment mesurez-vous les performances? La vitesse? Précision? La stabilité? Un test rapide dans Matlab donne les informations suivantes:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

Ainsi, la résolution d'un système unique avec une décomposition LU est environ trois fois plus rapide que la résolution avec une décomposition QR, au prix d'un demi-chiffre décimal de précision (cet exemple!).

Pedro
la source

Tous les mérites que vous avez suggérés sont les bienvenus.

faleichik

3

L'article que vous citez défend l'élimination gaussienne en disant que même s'il est numériquement instable, il a tendance à bien fonctionner sur des matrices aléatoires et puisque la plupart des matrices auxquelles on peut penser sont comme des matrices aléatoires, nous devrions être d'accord. Cette même affirmation peut être dite de nombreuses méthodes numériquement instables.

Considérez l'espace de toutes les matrices. Ces méthodes fonctionnent très bien presque partout. C'est 99,999 ...% de toutes les matrices que l'on pourrait créer n'auront aucun problème avec les méthodes instables. Il n'y a qu'une très petite fraction de matrices pour lesquelles GE et d'autres auront du mal.

Les problèmes qui préoccupent les chercheurs se situent généralement dans cette petite fraction.

Nous ne construisons pas de matrices au hasard. Nous construisons des matrices avec des propriétés très spéciales qui correspondent à des systèmes très spéciaux, non aléatoires. Ces matrices sont souvent mal conditionnées.

Géométriquement, vous pouvez considérer l'espace linéaire de toutes les matrices. Il y a un sous-espace zéro volume / mesure de matrices singulières qui traverse cet espace. De nombreux problèmes que nous construisons sont regroupés autour de ce sous-espace. Ils ne sont pas distribués au hasard.

À titre d'exemple, considérons l'équation ou la dispersion de la chaleur. Ces systèmes ont tendance à supprimer les informations du système (tous les états initiaux gravitent vers un seul état final) et, par conséquent, les matrices qui décrivent ces équations sont extrêmement singulières. Ce processus est très peu probable dans une situation aléatoire pourtant omniprésente dans les systèmes physiques.

MRocklin
la source

2

Si le système linéaire est initialement mal conditionné, peu importe la méthode que vous utilisez: les décompositions LU et QR donneront des résultats inexacts. QR ne peut gagner que dans les cas où le processus d'élimination gaussien "gâte" une bonne matrice. Le principal problème est que l'on ne connaît pas de cas pratiques d'un tel comportement.

faleichik

Pour la plupart des applications scientifiques, nous obtenons généralement des matrices clairsemées, symétriques, définies positives et / ou dominantes en diagonale. À quelques exceptions près, il existe une structure dans la matrice qui nous permet d'exploiter certaines techniques par rapport à l'élimination gaussienne traditionnelle.

Paul

@Paul: D'autre part, l'élimination gaussienne dense est l'endroit où la plupart du temps est passé dans la méthode multifrontale pour les matrices non symétriques clairsemées.

Jack Poulson

6

@Paul Il n'est tout simplement pas vrai que "la plupart des applications produisent des matrices SPD / dominantes en diagonale". Oui, il existe généralement une structure exploitable, mais des problèmes non symétriques et indéfinis sont extrêmement courants.

Jed Brown

4

"En cinquante ans d'informatique, aucun problème de matrice qui excite une instabilité explosive n'est apparu dans des circonstances naturelles." - LN Trefethen et D. Bau Ils donnent une analyse probabiliste intéressante dans leur livre.

JM

Quand les transformations orthogonales surpassent-elles l'élimination gaussienne?

Réponses:

Précision

Performance