Question Mystifié par qr.Q (): qu'est-ce qu'une matrice orthonormée en forme "compacte"?


R a un qr() fonction, qui effectue la décomposition QR en utilisant LINPACK ou LAPACK (dans mon expérience, ce dernier est 5% plus rapide). L'objet principal renvoyé est une matrice "qr" qui contient dans la matrice triangulaire supérieure R (c'est-à-dire R=qr[upper.tri(qr)]). Jusqu'ici tout va bien. La partie triangulaire inférieure de qr contient Q "sous forme compacte". On peut extraire Q de la décomposition qr en utilisant qr.Q(). Je voudrais trouver l'inverse de qr.Q(). En d'autres termes, j'ai Q et R, et je voudrais les placer dans un objet "qr". R est trivial mais Q ne l'est pas. Le but est de s'y appliquer qr.solve(), qui est beaucoup plus rapide que solve() sur les grands systèmes.


18
2018-06-13 05:37


origine


Réponses:


introduction

R utilise le LINPACK  dqrdc routine, par défaut, ou la LAPACK  DGEQP3 routine, si spécifié, pour calculer la décomposition QR. Les deux routines calculent la décomposition à l'aide de réflexions de maison. Une matrice mxn A est décomposée en une matrice orthogonale de taille économique mxn (Q) et une matrice triangulaire supérieure nxn (R) en tant que A = QR, où Q peut être calculé par le produit de t matrices de réflexion de m-1 et n: Q = H1H2... Ht.

Chaque matrice de réflexion Hje peut être représenté par un vecteur de longueur (m-i + 1). Par exemple, H1 nécessite un vecteur de longueur m pour le stockage compact. Toutes les entrées de ce vecteur, sauf une, sont placées dans la première colonne du triangle inférieur de la matrice d’entrée (la diagonale est utilisée par le facteur R). Par conséquent, chaque réflexion nécessite un scalaire de stockage supplémentaire, fourni par un vecteur auxiliaire (appelé $qraux dans le résultat de R qr).

La représentation compacte utilisée est différente entre les routines LINPACK et LAPACK.

La voie LINPACK

Une réflexion de maison est calculée comme Hje = I - vjevjeT/ pje, où je suis la matrice d'identité, pje est l'entrée correspondante dans $qrauxet vje est comme suit:

  • vje[1..i-1] = 0,
  • vje[i] = pje
  • vje[i + 1: m] = A [i + 1..m, i] (c'est-à-dire une colonne du triangle inférieur de A après l'appel qr)

Exemple LINPACK

Travaillons à travers le exemple de l'article de décomposition QR sur Wikipedia dans R.

La matrice en cours de décomposition est

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

Nous faisons la décomposition, et les parties les plus pertinentes du résultat sont présentées ci-dessous:

> Aqr = qr(A)
> Aqr
$qr
            [,1]         [,2] [,3]
[1,] -14.0000000  -21.0000000   14
[2,]   0.4285714 -175.0000000   70
[3,]  -0.2857143    0.1107692  -35

[snip...]

$qraux
[1]  1.857143  1.993846 35.000000

[snip...]

Cette décomposition a été faite (sous les couvertures) en calculant deux réflexions de Householder et en les multipliant par A pour obtenir R. Nous allons maintenant recréer les réflexions à partir des informations de $qr.

> p = Aqr$qraux   # for convenience
> v1 <- matrix(c(p[1], Aqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.8571429
[2,]  0.4285714
[3,] -0.2857143

> v2 <- matrix(c(0, p[2], Aqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.9938462
[3,] 0.1107692

> I = diag(3)   # identity matrix
> H1 = I - v1 %*% t(v1)/p[1]   # I - v1*v1^T/p[1]
> H2 = I - v2 %*% t(v2)/p[2]   # I - v2*v2^T/p[2]

> Q = H1 %*% H2
> Q
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

Maintenant, vérifions que le Q calculé ci-dessus est correct:

> qr.Q(Aqr)
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

Cela semble bon! Nous pouvons également vérifier que QR est égal à A.

> R = qr.R(Aqr)   # extract R from Aqr$qr
> Q %*% R
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

La voie LAPACK

Une réflexion de maison est calculée comme Hje = I - pjevjevjeT, où je suis la matrice d'identité, pje est l'entrée correspondante dans $qrauxet vje est comme suit:

  • vje[1..i-1] = 0,
  • vje[i] = 1
  • vje[i + 1: m] = A [i + 1..m, i] (c'est-à-dire une colonne du triangle inférieur de A après l'appel qr)

Il existe une autre tournure lors de l'utilisation de la routine LAPACK dans R: le pivot de colonne est utilisé, donc la décomposition résout un problème différent lié: AP = QR, où P est un matrice de permutation.

Exemple LAPACK

Cette section fait le même exemple qu'auparavant.

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> Bqr = qr(A, LAPACK=TRUE)
> Bqr
$qr
            [,1]        [,2]       [,3]
[1,] 176.2554964 -71.1694118   1.668033
[2,]  -0.7348557  35.4388886  -2.180855
[3,]  -0.1056080   0.6859203 -13.728129

[snip...]

$qraux
[1] 1.289353 1.360094 0.000000

$pivot
[1] 2 3 1

attr(,"useLAPACK")
[1] TRUE

[snip...]

Notez le $pivot champ; nous y reviendrons. Maintenant, nous générons Q à partir des informations du Aqr.

> p = Bqr$qraux   # for convenience
> v1 = matrix(c(1, Bqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.0000000
[2,] -0.7348557
[3,] -0.1056080


> v2 = matrix(c(0, 1, Bqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.0000000
[3,] 0.6859203


> H1 = I - p[1]*v1 %*% t(v1)   # I - p[1]*v1*v1^T
> H2 = I - p[2]*v2 %*% t(v2)   # I - p[2]*v2*v2^T
> Q = H1 %*% H2
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

Encore une fois, le Q calculé ci-dessus est en accord avec le Q. fourni par R

> qr.Q(Bqr)
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

Enfin, calculons QR.

> R = qr.R(Bqr)
> Q %*% R
     [,1] [,2] [,3]
[1,]  -51    4   12
[2,]  167  -68    6
[3,]   24  -41   -4

Remarquez la différence? QR est A avec ses colonnes permutées compte tenu de l'ordre dans Bqr$pivot au dessus.


19
2017-09-17 01:49



J'ai recherché ce même problème comme le demande l'OP et je ne pense pas que ce soit possible. Fondamentalement, la question OP est de savoir si le Q explicitement calculé permet de récupérer le H1 H2 ... Ht. Je ne pense pas que ce soit possible sans calculer le QR à partir de zéro, mais je serais également très intéressé de savoir s'il existe une telle solution.

J'ai un problème similaire à celui de l'OP, mais dans un contexte différent, mon algorithme itératif doit muter la matrice A en ajoutant des colonnes et / ou des lignes. La première fois, le QR est calculé en utilisant DGEQRF et donc le format LAPACK compact. Après que la matrice A soit mutée, par ex. Avec de nouvelles lignes, je peux rapidement construire un nouvel ensemble de réflecteurs ou de rotateurs qui annihileront les éléments non nuls de la diagonale la plus basse de mon R existant et construiront un nouveau R mais j'ai maintenant un ensemble de H1_old H2_old ... Hn_old et H1_new H2_new ... Hn_new (et de même tau) qui ne peuvent pas être mélangés en une seule représentation compacte QR. Les deux possibilités que j'ai sont, et peut-être l'OP a les deux mêmes possibilités:

  1. Conservez toujours Q et R séparés séparément, que vous les calculiez la première fois ou après chaque mise à jour au prix de flops supplémentaires, tout en maintenant la mémoire requise bien limitée.
  2. Respectez le format LAPACK compact, mais à chaque nouvelle mise à jour, conservez une liste de tous ces mini-ensembles de réflecteurs de mise à jour. Au moment de résoudre le système, on ferait un grand Q '* c, c'est-à-dire H1_u3 * H2_u3 * ... * Hn_u3 * H1_u2 * H2_u2 * ... * Hn_u2 * H1_u1 * H2_u1 ... * Hn_u1 * H1 * H2 * ... * Hn * c où ui est le numéro de mise à jour QR et qu'il y a beaucoup de multiplications à faire et de mémoire à suivre, mais certainement le moyen le plus rapide.

La longue réponse de David explique essentiellement ce qu'est le format compact QR, mais pas comment accéder à ce format QR compact ayant comme entrée les Q et R explicites calculés.


0
2018-03-27 14:06