Cálculo de variables de contraste polinomiales.

Dame una idea de cómo recodificar eficientemente una variable categórica (factor) en el conjunto de variables de contraste polinomiales ortogonales.

Para muchos tipos de variables de contraste (por ejemplo, desviación, simple, Helmert, etc.) el pase es:

1. Componga la matriz de coeficientes de contraste correspondiente al tipo.
2. Invierta o generalice-invierta para obtener la matriz de códigos.

Por ejemplo:

Suppose there is 3-group factor and we want to recode it into a set of deviation  contrast variables.
The last group is treated as reference. Then the contrast coefficients matrix L is

         Group1 Group2 Group3
   var1   2/3   -1/3   -1/3
   var2  -1/3    2/3   -1/3

and ginv(L) is then the sought-for coding matrix

         var1 var2
  Group1   1    0
  Group2   0    1
  Group3  -1   -1

(We might also use inv(L) instead if we add a row for constant, equal to 1/3, at the head of L.)

¿Existe la misma forma o similar para obtener variables de contraste polinomiales? En caso afirmativo, ¿cómo se vería la matriz C y cómo componerla? Si no, ¿cuál es la forma de calcular las variables de contraste polinomiales de manera eficiente (por ejemplo, por álgebra matricial)?

Como continuación de mi publicación anterior sobre este tema , quiero compartir una exploración tentativa (aunque incompleta) de las funciones detrás del álgebra lineal y las funciones R relacionadas. Se supone que esto es un trabajo en progreso.

Parte de la opacidad de las funciones tiene que ver con la forma "compacta" de la descomposición Householder . La idea detrás de la descomposición de Householder es reflejar los vectores a través de un hiperplano determinado por un vector unitario como en el diagrama a continuación, pero seleccionando este plano de manera intencional para proyectar cada vector de columna de la matriz original en el vector de unidad estándar. El vector normalizado norma-2 se puede usar para calcular las diferentes transformaciones de Householder .u A e 1 1 u I - $\mathbf {QR}$$\mathbf u$$\bf A$$\bf e_1$$1$$\bf u$$\mathbf{ I - 2\, uu^T x}$

La proyección resultante se puede expresar como

$\text{sign}(x_i=x_1)\times \lVert x \rVert \begin{bmatrix}1\\0\\0\\\vdots\\0\end{bmatrix}+\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}$

El vector representa la diferencia entre los vectores de columna en la matriz que queremos descomponer y los vectores correspondientes a la reflexión a través del subespacio o "espejo" determinado por .$\bf v$$\bf x$$\bf A$$\bf y$$\bf u$

El método utilizado por LAPACK libera la necesidad de almacenar la primera entrada en los reflectores Householder al convertirlos en 's. En lugar de normalizar el vector a con , es solo la entrada del primer puño la que se convierte en ; sin embargo, estos nuevos vectores, , todavía pueden usarse como vectores direccionales.$1$$\bf v$$\bf u$$\lVert u\rVert= 1$$1$$\bf w$

La belleza del método es que dado que en una descomposición es triangular superior, en realidad podemos aprovechar los elementos en debajo de la diagonal para completarlos con estos reflectores. Afortunadamente, las entradas principales en estos vectores son todas iguales a , evitando un problema en la diagonal "disputada" de la matriz: sabiendo que son todas , no necesitan ser incluidas, y pueden ceder la diagonal a las entradas de .$\bf R$$\bf QR$$0$$\bf R$$\bf w$$1$$1$$\bf R$

La matriz "compacta QR" en la función qr()$qrpuede entenderse aproximadamente como la adición de la matriz y la matriz triangular inferior de "almacenamiento" para los reflectores "modificados".$\mathbf R$

La proyección de Householder seguirá teniendo la forma , pero no trabajaremos con ( ), sino con un vector , de los cuales solo se garantiza que la primera entrada sea , y$\mathbf{ I - 2\, uu^T x}$$\bf u$$\lVert \bf x \rVert=1$$\bf w$$1$

$\mathbf{ I - 2\, uu^T x}=\mathbf{ I - 2\, \frac{w}{\lVert w \rVert} \frac{w^T}{\lVert w \rVert} x}=\mathbf{ I - 2\, \frac{w\,w^T}{\lVert w \rVert^2}\,x}\tag{1}$ .

Uno supondría que estaría bien almacenar estos reflectores debajo de la diagonal o excluyendo la primera entrada de , y llamarlo un día. Sin embargo, las cosas nunca son tan fáciles. En cambio, lo que se almacena debajo de la diagonal en es una combinación de y los coeficientes en la transformación Householder expresados ​​como (1), de modo que, definiendo como:$\bf w$$\bf R$$1$qr()$qr$\bf w$$\text{tau}$

$\Large \tau = \frac{w^T\,w}{2}= \frac{\lVert \bf w \rVert}{2}$ , los reflectores se pueden expresar como . Estos vectores "reflectores" son los que se almacenan justo debajo de en el llamado "compacto ".$\text{reflectors}=\large \bf w/\tau$$\bf R$$\bf QR$

Ahora estamos a un grado de distancia de los vectores , y la primera entrada ya no es , por lo tanto, la salida de deberá incluir la clave para restaurarlos ya que insistimos en excluir la primera entrada de los vectores "reflectores" a encajar todo en . Entonces, ¿estamos viendo los valores en la salida? Bueno, no, eso sería predecible. En cambio, en la salida de (donde se almacena esta clave) encontramos .$\bf w$$1$qr()qr()$qr$\tau$qr()$qraux$\rho=\frac{\sum \text{reflectors}^2}{2}= \frac{\bf w^Tw}{\tau^2} / 2$

Enmarcado en rojo a continuación, vemos los "reflectores" ( ), excluyendo su primera entrada.$\bf w/\tau$

Todo el código está aquí , pero como esta respuesta trata sobre la intersección de la codificación y el álgebra lineal, pegaré la salida para facilitar:

options(scipen=999)
set.seed(13)
(X = matrix(c(rnorm(16)), nrow=4, byrow=F))
           [,1]      [,2]       [,3]       [,4]
[1,]  0.5543269 1.1425261 -0.3653828 -1.3609845
[2,] -0.2802719 0.4155261  1.1051443 -1.8560272
[3,]  1.7751634 1.2295066 -1.0935940 -0.4398554
[4,]  0.1873201 0.2366797  0.4618709 -0.1939469

Ahora escribí la función de la House()siguiente manera:

   House = function(A){
    Q = diag(nrow(A))
    reflectors = matrix(0,nrow=nrow(A),ncol=ncol(A))
    for(r in 1:(nrow(A) - 1)){ 
        # We will apply Householder to progressively the columns in A, decreasing 1 element at a time.
        x = A[r:nrow(A), r] 
        # We now get the vector v, starting with first entry = norm-2 of x[i] times 1
        # The sign is to avoid computational issues
        first = (sign(x[1]) * sqrt(sum(x^2))) +  x[1]
        # We get the rest of v, which is x unchanged, since e1 = [1, 0, 0, ..., 0]
        # We go the the last column / row, hence the if statement:
        v = if(length(x) > 1){c(first, x[2:length(x)])}else{v = c(first)}
        # Now we make the first entry unitary:
        w = v/first
        # Tau will be used in the Householder transform, so here it goes:
        t = as.numeric(t(w)%*%w) / 2
        # And the "reflectors" are stored as in the R qr()$qr function:
        reflectors[r: nrow(A), r] = w/t
        # The Householder tranformation is:
        I = diag(length(r:nrow(A)))
        H.transf = I - 1/t * (w %*% t(w))
        H_i  = diag(nrow(A))
        H_i[r:nrow(A),r:ncol(A)] = H.transf
        # And we apply the Householder reflection - we left multiply the entire A or Q
        A = H_i %*% A
        Q = H_i %*% Q
    }
    DECOMPOSITION = list("Q"= t(Q), "R"= round(A,7), 
            "compact Q as in qr()$qr"=  
            ((A*upper.tri(A,diag=T))+(reflectors*lower.tri(reflectors,diag=F))), 
            "reflectors" = reflectors,
            "rho"=c(apply(reflectors[,1:(ncol(reflectors)- 1)], 2, 
                function(x) sum(x^2) / 2), A[nrow(A),ncol(A)]))
    return(DECOMPOSITION)
}

Comparemos la salida con las funciones integradas de R. Primero la función casera:

(H = House(X))
$Q
            [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

$R
          [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$`compact Q as in qr()$qr`
            [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$reflectors
            [,1]        [,2]      [,3] [,4]
[1,]  1.29329367  0.00000000 0.0000000    0
[2,] -0.14829152  1.73609434 0.0000000    0
[3,]  0.93923665 -0.67574886 1.7817597    0
[4,]  0.09911084  0.03909742 0.6235799    0

$rho
[1] 1.2932937 1.7360943 1.7817597 0.4754876

a las funciones R:

qr.Q(qr(X))
            [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

qr.R(qr(X))
          [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$qr
            [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$qraux
[1] 1.2932937 1.7360943 1.7817597 0.4754876