¿Cuál es la razón por la que LAPACK usa

La rutina QR de LAPACK almacena Q como reflectores de cabeza de familia. Escala el vector de reflexión $v$ con $1/v_1$ , por lo que el primer elemento del resultado se convierte en $1$ , por lo que no tiene que almacenarse. Y almacena un vector $\tau$ separado , que contiene los factores de escala necesarios. Entonces una matriz reflectora es así:

donde $v$ no está normalizado. Mientras que, en los libros de texto, la matriz reflectora es

donde $v$ está normalizado.

¿Por qué LAPACK escala $v$ con $1/v_1$ , en lugar de normalizarlo?

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

(La razón de mi pregunta es que estoy escribiendo una rutina QR y SVD, y me gustaría saber la razón de esta decisión, si necesito seguirla o no)

Es la variante bloqueada de Householder-QR la que impulsa este diseño. Si miras en el libro de Golub y Van Loan (Capítulo 5.2 más o menos), hablan de cómo las iteraciones k del algoritmo pueden bloquearse juntas al acumular los reflectores individuales en un reflector de rango k de la forma , donde tanto como son matrices "altas y delgadas" con un tamaño . Este algoritmo hace más trabajo pero es más rápido en la práctica porque es rico en llamadas gemm (). Desafortunadamente, es un desperdicio de almacenamiento debido a la necesidad de representar y independiente.$\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$$\mathbf W$$\mathbf Y$$n \times k$$\mathbf W$$\mathbf Y$

En un artículo posterior (citado a continuación), Van Loan describe una estructura de datos "simétrica" ​​más eficiente, un reflector de bloque de la forma . Aquí todavía es , pero el requisito de flop / almacenamiento para formar se ha eliminado al introducir , una pequeña matriz triangular superior . Aunque la necesidad de multiplicar por introduce una pequeña cantidad de trabajo extra, generalmente es una ganancia neta porque .$\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$

Dentro de LAPACK, el algoritmo no bloqueado es realmente un caso limitante de del algoritmo de bloqueo, hasta la elección de símbolos (lo que nos lleva a , una pequeña versión de Triángulo ).$k \rightarrow 1$$\tau$$1\times1$$\mathbf T$

Cita: Schreiber, Robert y Charles Van Loan. "Una representación WY de almacenamiento eficiente para productos de transformaciones de Householder". Revista SIAM sobre Computación Científica y Estadística 10.1 (1989): 53-57.

No tiene que almacenar , puede volver a calcularlo desde el resto del vector. (Puede calcular desde las otras entradas también en la versión normalizada, pero es claramente un cálculo inestable debido a esas sustracciones).$\tau$$v_1$

En realidad, puede reutilizar la parte triangular inferior de para almacenar , de modo que la factorización se calcule completamente en el lugar. Lapack se preocupa mucho por estas versiones in situ de algoritmos.$R$$v_2,...v_n$

Mi sugerencia se basa en la documentación para Intel MKL https://software.intel.com/en-us/mkl-developer-reference-c-geqrf . Parece que los valores en y por encima de la diagonal de la salida almacenan R, por lo que solo queda un triángulo inferior para Q. Parece natural usar almacenamiento adicional para los factores de escala.