Quelle est la raison pour laquelle LAPACK utilise

La routine QR de LAPACK stocke Q en tant que réflecteurs domestiques. Il met à l'échelle le vecteur de réflexion $v$ avec $1/v_1$ , de sorte que le premier élément du résultat devient $1$ , il n'a donc pas besoin d'être stocké. Et il stocke un vecteur $\tau$ séparé , qui contient les facteurs d'échelle nécessaires. Donc, une matrice de réflecteur est comme ceci:

où $v$ n'est pas normalisé. Alors que, dans les manuels, la matrice du réflecteur est

où $v$ est normalisé.

Pourquoi LAPACK met-il à l'échelle $v$ avec $1/v_1$ , au lieu de le normaliser?

$\tau$$v_1$$H$$\tau$$2$$v$$\sqrt 2/\|v\|$

(La raison de ma question est que j'écris une routine QR et SVD, et j'aimerais connaître la raison de cette décision, si je dois la suivre ou non)

C'est la variante bloquée de Householder-QR qui anime cette conception. Si vous regardez dans le livre de Golub et Van Loan (Ch 5.2 ou plus), ils parlent de la façon dont les k-itérations de l'algorithme peuvent être bloquées ensemble en accumulant les réflecteurs individuels dans un réflecteur de rang k de la forme , où et sont des matrices "hautes et maigres" de taille . Cet algorithme fait plus de travail mais est plus rapide dans la pratique car il est riche en appels gemm (). Malheureusement, cela est un gaspillage de stockage en raison de la nécessité de représenter et indépendamment.$\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$$\mathbf W$$\mathbf Y$$n \times k$$\mathbf W$$\mathbf Y$

Dans un article ultérieur (cité ci-dessous), Van Loan décrit une structure de données "symétrisée" plus efficace, un réflecteur de bloc de la forme . Ici est toujours , mais l'exigence de flop / stockage pour former a été éliminée en introduisant , une petite matrice triangulaire supérieure . Bien que la nécessité de multiplier par introduise une petite quantité de travail supplémentaire, c'est généralement un gain net car .$\mathbf I + \mathbf Y \mathbf T \mathbf Y^{\mathrm T}$$\mathbf Y$$n \times k$$\mathbf W$$\mathbf T$$k \times k$$\mathbf T$$k << n$

Au sein de LAPACK, l'algorithme non bloqué n'est vraiment qu'un cas limitant de l'algorithme de bloc, jusqu'au choix des symboles (ce qui nous amène à , une petite version du Triangle ).$k \rightarrow 1$$\tau$$1\times1$$\mathbf T$

Référence: Schreiber, Robert et Charles Van Loan. "Une représentation WY efficace pour le stockage des produits des transformations des ménages." SIAM Journal on Scientific and Statistical Computing 10.1 (1989): 53-57.

Vous n'avez pas besoin de stocker , vous pouvez le recalculer à partir du reste du vecteur. (Vous pouvez recalculer partir des autres entrées également dans la version normalisée, mais il s'agit clairement d'un calcul instable en raison de ces soustractions.)$\tau$$v_1$

En fait, vous pouvez réutiliser la partie triangulaire inférieure de pour stocker , afin que la factorisation soit entièrement calculée en place. Lapack se soucie beaucoup de ces versions d'algorithmes en place.$R$$v_2,...v_n$

Ma suggestion est basée sur la documentation d'Intel MKL https://software.intel.com/en-us/mkl-developer-reference-c-geqrf . Il ressemble aux valeurs sur et au-dessus de la diagonale du magasin de sortie R, il ne reste donc que le triangle inférieur pour Q. Il semble naturel d'utiliser un stockage supplémentaire pour les facteurs d'échelle.