r linear-algebra lapack blas

Mystified by qr.Q(): ¿qué es una matriz ortonormal en forma "compacta"?



linear-algebra lapack (2)

R tiene una función qr() , que realiza la descomposición QR usando LINPACK o LAPACK (en mi experiencia, esta última es 5% más rápida). El objeto principal devuelto es una matriz "qr" que contiene en la matriz triangular superior R (es decir, R=qr[upper.tri(qr)] ). Hasta ahora tan bueno. La parte triangular inferior de qr contiene Q "en forma compacta". Uno puede extraer Q de la descomposición de qr usando qr.Q() . Me gustaría encontrar el inverso de qr.Q() . En otras palabras, tengo Q y R, y me gustaría ponerlas en un objeto "qr". R es trivial pero Q no lo es. El objetivo es aplicarlo a qr.solve() , que es mucho más rápido que solve() en sistemas grandes.

Introducción

R utiliza la rutina LINPACK dqrdc , por defecto, o la rutina LAPACK DGEQP3 , cuando se especifica, para calcular la descomposición de QR. Ambas rutinas computan la descomposición usando reflexiones de los dueños de casa. Una matriz mxn A se descompone en una matriz ortogonal de tamaño económico mxn (Q) y una matriz triangular superior nxn (R) como A = QR, donde Q puede calcularse mediante el producto de t matrices de reflexión de cabeza de familia, siendo t la menor de m-1 yn: Q = H 1 H 2 ... H t .

Cada matriz de reflexión H i puede representarse por un vector de longitud (m-i + 1). Por ejemplo, H1 requiere un vector de longitud-m para un almacenamiento compacto. Todas las entradas de este vector, excepto una, se colocan en la primera columna del triángulo inferior de la matriz de entrada (la diagonal es utilizada por el factor R). Por lo tanto, cada reflexión necesita un escalar más de almacenamiento, y esto es proporcionado por un vector auxiliar (llamado $qraux en el resultado de la qr de R).

La representación compacta utilizada es diferente entre las rutinas LINPACK y LAPACK.

El modo LINPACK

Una reflexión del $qraux se calcula como H i = I - v i v i T / p i , donde I es la matriz de identidad, p i es la entrada correspondiente en $qraux y v i es la siguiente:

  • v i [1..i-1] = 0,
  • v i [i] = p i
  • v i [i + 1: m] = A [i + 1..m, i] (es decir, una columna del triángulo inferior de A después de llamar qr )

Ejemplo de LINPACK

Veamos el ejemplo del artículo de descomposición de QR en Wikipedia en R.

La matriz que se descompone es

> 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

Hacemos la descomposición, y las partes más relevantes del resultado se muestran a continuación:

> 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...]

Esta descomposición se realizó (debajo de las coberturas) al calcular dos reflexiones de los propietarios de viviendas y multiplicarlas por A para obtener R. Ahora recrearemos las reflexiones de la información en $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

Ahora verifiquemos que la Q calculada arriba es correcta:

> 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

¡Se ve bien! También podemos verificar que QR es igual a 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

El modo LAPACK

Una reflexión del $qraux se calcula como H i = I - p i v i v i T , donde I es la matriz de identidad, p i es la entrada correspondiente en $qraux y v i es la siguiente:

  • v i [1..i-1] = 0,
  • v i [i] = 1
  • v i [i + 1: m] = A [i + 1..m, i] (es decir, una columna del triángulo inferior de A después de llamar qr )

Hay otro giro cuando se usa la rutina LAPACK en R: se usa la rotación de columnas, por lo que la descomposición está resolviendo un problema relacionado diferente: AP = QR, donde P es una matriz de permutación .

Ejemplo de LAPACK

Esta sección hace el mismo ejemplo que antes.

> 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...]

Note el campo $pivot ; Volveremos a eso. Ahora generamos Q a partir de la información del 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

Una vez más, la Q calculada anteriormente concuerda con la Q proporcionada por la 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

Por último, vamos a calcular QR.

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

¿Nota la diferencia? QR es A con sus columnas permutadas dado el orden en Bqr$pivot anterior.

He investigado para este mismo problema, como pregunta el OP, y no creo que sea posible. Básicamente, la pregunta del OP es si tener la Q calculada explícitamente, se puede recuperar el H1 H2 ... Ht. No creo que esto sea posible sin calcular el QR desde cero, pero también estaría muy interesado en saber si existe tal solución.

Tengo un problema similar al OP, pero en un contexto diferente, mi algoritmo iterativo necesita mutar la matriz A agregando columnas y / o filas. La primera vez, el QR se calcula utilizando DGEQRF y, por lo tanto, el formato compacto LAPACK. Después de que la matriz A sea mutada, por ejemplo, con nuevas filas puedo construir rápidamente un nuevo conjunto de reflectores o rotadores que aniquilarán los elementos que no sean cero de la diagonal más baja de mi R existente y construiré una nueva R, pero ahora tengo un conjunto de H1_old H2_old ... Hn_old y H1_new H2_new ... Hn_new (y de manera similar, tau) que no se pueden mezclar en una sola representación compacta de QR. Las dos posibilidades que tengo son, y tal vez el OP tiene las mismas dos posibilidades:

  1. Siempre mantenga la Q y la R separadas de manera explícita, ya sea cuando se calculan la primera vez o después de cada actualización a costa de fracasos adicionales pero manteniendo la memoria requerida bien delimitada.
  2. Manténgase en el formato compacto de LAPACK, pero cada vez que llegue una nueva actualización, mantenga una lista de todos estos mini juegos de reflectores de actualización. En el punto de resolver el sistema, se haría una gran Q ''* c, es decir, H1_u3 * H2_u3 * ... * Hn_u3 * H1_u2 * H2_u2 * ... * Hn_u2 * H1_u1 * H2_u1 ... * Hn_u1 * H1 * H2 * ... * Hn * c donde ui es el número de actualización de QR y esto es potencialmente una gran cantidad de multiplicaciones que hacer y memoria para realizar un seguimiento de, pero definitivamente la forma más rápida.

La respuesta larga de David explica básicamente qué es el formato QR compacto, pero no cómo obtener este formato QR compacto con la Q y R calculadas explícitas como entrada.