Calcul des variables de contraste polynomiales

Veuillez me donner une idée de la façon de recoder efficacement une variable catégorielle (facteur) dans l'ensemble de variables de contraste polynomiales orthogonales.

Pour de nombreux types de variables de contraste (par exemple, déviation, simple, Helmert, etc.), la réussite est la suivante:

1. Composez la matrice des coefficients de contraste correspondant au type.
2. Inversez ou généralisez-l'inverse pour obtenir la matrice de codes.

Par exemple:

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-t-il la même méthode ou une méthode similaire pour obtenir des variables de contraste polynomiales? Si oui, à quoi ressemblerait la matrice C et comment la composer? Si non, quelle est encore la façon de calculer efficacement les variables de contraste polynomiales (par exemple par algèbre matricielle).

Comme enchaînement à mon post précédent sur ce sujet, je veux partager une exploration provisoire (quoique incomplète) des fonctions derrière l'algèbre linéaire et les fonctions R connexes. C'est censé être un travail en cours.

Une partie de l'opacité des fonctions est liée à la forme "compacte" de la décomposition Householder . L'idée derrière la décomposition du Householder est de refléter les vecteurs à travers un hyperplan déterminé par un vecteur-unité comme dans le diagramme ci-dessous, mais en choisissant ce plan de manière ciblée afin de projeter chaque vecteur de colonne de la matrice d'origine sur le vecteur d'unité standard . Le vecteur normalisé-2 normalisé peut être utilisé pour calculer les différentes transformations 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 projection résultante peut être exprimée comme

$\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}$

Le vecteur représente la différence entre les vecteurs colonne dans la matrice que nous voulons décomposer et les vecteurs correspondant à la réflexion à travers le sous-espace ou "miroir" déterminée par .$\bf v$$\bf x$$\bf A$$\bf y$$\bf u$

La méthode utilisée par LAPACK libère le besoin de stockage de la première entrée dans les réflecteurs Householder en les transformant en . Au lieu de normaliser le vecteur en avec , c'est juste la première entrée qui est convertie en ; pourtant, ces nouveaux vecteurs - appelez-les peuvent toujours être utilisés comme vecteurs directionnels.$1$$\bf v$$\bf u$$\lVert u\rVert= 1$$1$$\bf w$

La beauté de la méthode est que, étant donné que dans une décomposition est triangulaire supérieur, nous pouvons réellement profiter des éléments dans sous la diagonale pour les remplir avec ces réflecteurs . Heureusement, les entrées principales de ces vecteurs sont toutes égales à , ce qui évite un problème dans la diagonale "contestée" de la matrice: sachant qu'elles sont toutes elles n'ont pas besoin d'être incluses et peuvent donner la diagonale aux entrées de .$\bf R$$\bf QR$$0$$\bf R$$\bf w$$1$$1$$\bf R$

La matrice "QR compacte" dans la fonction qr()$qrpeut être comprise en gros comme l'addition de la matrice et de la matrice triangulaire "de stockage" inférieure pour les réflecteurs "modifiés".$\mathbf R$

La projection Householder aura toujours la forme , mais nous ne travaillerons pas avec ( ), mais plutôt avec un vecteur , dont seule la première entrée est garantie à , et$\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}$ .

On pourrait supposer qu'il serait très bien de stocker ces réflecteurs sous la diagonale ou excluant la première entrée de , et l'appeler un jour. Cependant, les choses ne sont jamais aussi faciles. Au lieu de cela, ce qui est stocké en dessous de la diagonale dans est une combinaison de et des coefficients dans la transformation Householder exprimée comme (1), de sorte que, définissant comme:$\bf w$$\bf R$$1$qr()$qr$\bf w$$\text{tau}$

$\Large \tau = \frac{w^T\,w}{2}= \frac{\lVert \bf w \rVert}{2}$ , les réflecteurs peuvent être exprimés sous la forme . Ces vecteurs "réflecteurs" sont ceux stockés juste sous dans le soi-disant "compact ".$\text{reflectors}=\large \bf w/\tau$$\bf R$$\bf QR$

Maintenant, nous sommes à un degré des vecteurs , et la première entrée n'est plus , d'où la sortie de devra inclure la clé pour les restaurer car nous insistons pour exclure la première entrée des vecteurs "réflecteur" à s'adapter à tout . Voyons-nous donc les valeurs dans la sortie? Eh bien, non, ce serait prévisible. Au lieu de cela, dans la sortie de (où cette clé est stockée), nous trouvons .$\bf w$$1$qr()qr()$qr$\tau$qr()$qraux$\rho=\frac{\sum \text{reflectors}^2}{2}= \frac{\bf w^Tw}{\tau^2} / 2$

Ainsi encadrés en rouge ci-dessous, on voit les "réflecteurs" ( ), à l'exclusion de leur première entrée.$\bf w/\tau$

Tout le code est ici , mais comme cette réponse concerne l'intersection du codage et de l'algèbre linéaire, je vais coller la sortie pour plus de facilité:

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

Maintenant, j'ai écrit la fonction House()comme suit:

   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)
}

Comparons la sortie aux fonctions intégrées R. D'abord la fonction maison:

(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

aux fonctions 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