¿Cómo realizar la validación cruzada para PCA para determinar el número de componentes principales?

Estoy tratando de escribir mi propia función para el análisis de componentes principales, PCA (por supuesto, ya hay mucho escrito pero solo estoy interesado en implementar cosas por mí mismo). El principal problema que encontré es el paso de validación cruzada y el cálculo de la suma de cuadrados pronosticada (PRENSA). No importa qué validación cruzada use, es una pregunta principalmente sobre la teoría detrás, pero considere la validación cruzada de dejar uno fuera (LOOCV). De la teoría descubrí que para realizar LOOCV necesitas:

1. eliminar un objeto
2. escalar el resto
3. realizar PCA con cierto número de componentes
4. escalar el objeto eliminado de acuerdo con los parámetros obtenidos en (2)
5. predecir el objeto de acuerdo con el modelo de PCA
6. calcular PRENSA para este objeto
7. volver a realizar el mismo algoritmo a otros objetos
8. resumir todos los valores de PRENSA
9. lucro

Debido a que soy muy nuevo en el campo, para asegurarme de tener razón, comparo los resultados con la salida de algún software que tengo (también para escribir algún código sigo las instrucciones en el software). Obtengo los mismos resultados al calcular la suma residual de los cuadrados y $R^2$ , pero calcular la PRENSA es un problema.

¿Podría decirme si lo que implemento en el paso de validación cruzada es correcto o no?

case 'loocv'

% # n - number of objects
% # p - number of variables
% # vComponents - the number of components used in CV
dataSets = divideData(n,n); 
         % # it is just a variable responsible for creating datasets for CV 
         % #  (for LOOCV datasets will be equal to [1, 2, 3, ... , n]);'
tempPRESS = zeros(n,vComponents);

for j = 1:n
  Xmodel1 = X; % # X - n x p original matrix
  Xmodel1(dataSets{j},:) = []; % # delete the object to be predicted
  [Xmodel1,Xmodel1shift,Xmodel1div] = skScale(Xmodel1, 'Center', vCenter, 
                                              'Scaling', vScaling); 
          % # scale the data and extract the shift and scaling factor
  Xmodel2 = X(dataSets{j},:); % # the object to be predicted
  Xmodel2 = bsxfun(@minus,Xmodel2,Xmodel1shift); % # shift and scale the object
  Xmodel2 = bsxfun(@rdivide,Xmodel2,Xmodel1div);
  [Xscores2,Xloadings2] = myNipals(Xmodel1,0.00000001,vComponents); 
          % # the way to calculate the scores and loadings
                % # Xscores2 - n x vComponents matrix
                % # Xloadings2 - vComponents x p matrix
  for i = 1:vComponents
    tempPRESS(j,i) = sum(sum((Xmodel2* ...
       (eye(p) - transpose(Xloadings2(1:i,:))*Xloadings2(1:i,:))).^2));
  end
end
PRESS = sum(tempPRESS,1);

En el software ( PLS_Toolbox ) funciona así:

for i = 1:vComponents
    tempPCA = eye(p) - transpose(Xloadings2(1:i,:))*Xloadings2(1:i,:);
    for kk = 1:p
        tempRepmat(:,kk) = -(1/tempPCA(kk,kk))*tempPCA(:,kk);
          % # this I do not understand
        tempRepmat(kk,kk) = -1; 
          % # here is some normalization that I do not get
    end 
    tempPRESS(j,i) = sum(sum((Xmodel2*tempRepmat).^2)); 
end

Entonces, hacen una normalización adicional usando esta tempRepmatvariable: la única razón que encontré fue que aplican LOOCV para PCA robusto. Desafortunadamente, el equipo de soporte no quiso responder mi pregunta ya que solo tengo una versión demo de su software.

Lo que está haciendo está mal: ¡no tiene sentido calcular PRESS para PCA así! Específicamente, el problema radica en su paso # 5.

---

Enfoque ingenuo de PRENSA para PCA

$n$$d$$\mathbf x^{(i)} \in \mathbb R^d, \, i=1 \dots n$$\mathbf x^{(i)}$$\mathbf X^{(-i)}$$k$$\mathbf U^{(-i)}$$\left \|\mathbf x^{(i)} - \hat{\mathbf x}^{(i)}\right\|^2 = \left\|\mathbf x^{(i)} - \mathbf U^{(-i)} [\mathbf U^{(-i)}]^\top \mathbf x^{(i)}\right\|^2$$i$

$$\mathrm{PRESS} \overset{?}{=} \sum_{i=1}^n \left\|\mathbf x^{(i)} - \mathbf U^{(-i)} [\mathbf U^{(-i)}]^\top \mathbf x^{(i)}\right\|^2.$$

Por simplicidad, estoy ignorando los problemas de centrado y escalado aquí.

El enfoque ingenuo está mal

El problema anterior es que usamos para calcular la predicción , y eso es algo muy malo.$\mathbf x^{(i)}$$\hat{\mathbf x}^{(i)}$

Tenga en cuenta la diferencia crucial para un caso de regresión, donde la fórmula para el error de reconstrucción es básicamente la misma , pero la predicción se calcula usando las variables predictoras y no usando . Esto no es posible en PCA, porque en PCA no hay variables dependientes e independientes: todas las variables se tratan juntas.$\left\|\mathbf y^{(i)} - \hat{\mathbf y}^{(i)}\right\|^2$$\hat{\mathbf y}^{(i)}$$\mathbf y^{(i)}$

En la práctica, significa que PRESS, como se calculó anteriormente, puede disminuir al aumentar el número de componentes y nunca alcanzar un mínimo. Lo que llevaría a pensar que todos los componentes son significativos. O tal vez en algunos casos alcanza un mínimo, pero aún tiende a sobreajustar y sobreestimar la dimensionalidad óptima.$k$$d$

Un enfoque correcto

Hay varios enfoques posibles, ver Bro et al. (2008) Validación cruzada de modelos de componentes: una mirada crítica a los métodos actuales para una visión general y comparación. Un enfoque es dejar de lado una dimensión de un punto de datos a la vez (es decir, lugar de ), para que los datos de entrenamiento se conviertan en una matriz con un valor faltante , y luego predecir ("imputar") este valor faltante con PCA. (Por supuesto, uno puede mantener al azar una fracción mayor de elementos de la matriz, por ejemplo, 10%). El problema es que computar PCA con valores faltantes puede ser computacionalmente bastante lento (depende del algoritmo EM), pero debe repetirse aquí muchas veces. Actualización: ver http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/$x^{(i)}_j$$\mathbf x^{(i)}$ para una discusión agradable y la implementación de Python (PCA con valores faltantes se implementa a través de mínimos cuadrados alternos).

Un enfoque que me pareció mucho más práctico es dejar de lado un punto de datos a la vez, calcular PCA en los datos de entrenamiento (exactamente como se indica arriba), pero luego recorrer las dimensiones de , déjelos fuera uno a la vez y calcule un error de reconstrucción usando el resto. Esto puede ser bastante confuso al principio y las fórmulas tienden a volverse bastante desordenadas, pero la implementación es bastante sencilla. Permítanme primero dar la fórmula (algo aterradora) y luego explicarla brevemente:$\mathbf x^{(i)}$$\mathbf x^{(i)}$

$$\mathrm{PRESS_{PCA}} = \sum_{i=1}^n \sum_{j=1}^d \left|x^{(i)}_j - \left[\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j}\right]_j\right|^2.$$

Considere el bucle interno aquí. Dejamos un punto y calculamos componentes principales en los datos de entrenamiento, . Ahora mantenemos cada valor como prueba y utilizamos las dimensiones restantes para realizar la predicción . La predicción es la coordenada de "la proyección" (en el sentido de los mínimos cuadrados) de en el subespacio por . Para calcularlo, encuentre un punto en el espacio de la PC que esté más cerca de$\mathbf x^{(i)}$$k$$\mathbf U^{(-i)}$$x^{(i)}_j$$\mathbf x^{(i)}_{-j} \in \mathbb R^{d-1}$$\hat{x}^{(i)}_j$$j$$\mathbf x^{(i)}_{-j}$$\mathbf U^{(-i)}$$\hat{\mathbf z}$$\mathbb R^k$$\mathbf x^{(i)}_{-j}$ calculando donde es con la fila expulsado, y significa pseudoinverso. Ahora asigne al espacio original: y tome su coordenada . $\hat{\mathbf z} = \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j} \in \mathbb R^k$$\mathbf U^{(-i)}_{-j}$$\mathbf U^{(-i)}$$j$$[\cdot]^+$$\hat{\mathbf z}$$\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j}$$j$$[\cdot]_j$

Una aproximación al enfoque correcto

No entiendo la normalización adicional utilizada en PLS_Toolbox, pero aquí hay un enfoque que va en la misma dirección.

Hay otra forma de mapear en el espacio de componentes principales: , es decir, simplemente tome la transposición en lugar de pseudo-inversa. En otras palabras, la dimensión que queda fuera para la prueba no se cuenta en absoluto, y los pesos correspondientes también se eliminan simplemente. Creo que esto debería ser menos preciso, pero a menudo podría ser aceptable. Lo bueno es que la fórmula resultante ahora se puede vectorizar de la siguiente manera (omito el cálculo):$\mathbf x^{(i)}_{-j}$$\hat{\mathbf z}_\mathrm{approx} = \left [\mathbf U^{(-i)}_{-j}\right]^\top\mathbf x^{(i)}_{-j}$

$$\begin{align} \mathrm{PRESS_{PCA,\,approx}} &= \sum_{i=1}^n \sum_{j=1}^d \left|x^{(i)}_j - \left[\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^\top\mathbf x^{(i)}_{-j}\right]_j\right|^2 \\ &= \sum_{i=1}^n \left\|\left(\mathbf I - \mathbf U \mathbf U^\top + \mathrm{diag}\{\mathbf U \mathbf U^\top\}\right) \mathbf x^{(i)}\right\|^2, \end{align}$$

donde escribí como para compacidad, y significa poner todos los elementos no diagonales a cero. ¡Tenga en cuenta que esta fórmula se ve exactamente como la primera (PRENSA ingenua) con una pequeña corrección! Tenga en cuenta también que esta corrección solo depende de la diagonal de , como en el código PLS_Toolbox. Sin embargo, la fórmula sigue siendo diferente de lo que parece implementarse en PLS_Toolbox, y no puedo explicar esta diferencia.$\mathbf U^{(-i)}$$\mathbf U$$\mathrm{diag}\{\cdot\}$$\mathbf U \mathbf U^\top$

Actualización (febrero de 2018): anteriormente llamé a un procedimiento "correcto" y otro "aproximado", pero ya no estoy tan seguro de que esto sea significativo. Ambos procedimientos tienen sentido y creo que ninguno es más correcto. Realmente me gusta que el procedimiento "aproximado" tenga una fórmula más simple. Además, recuerdo que tenía un conjunto de datos donde el procedimiento "aproximado" arrojó resultados que parecían más significativos. Desafortunadamente, ya no recuerdo los detalles.

---

Ejemplos

Así es como se comparan estos métodos para dos conjuntos de datos conocidos: el conjunto de datos Iris y el conjunto de datos de vino. Tenga en cuenta que el método ingenuo produce una curva decreciente monotónicamente, mientras que otros dos métodos producen una curva con un mínimo. Tenga en cuenta además que en el caso de Iris, el método aproximado sugiere 1 PC como el número óptimo, pero el método pseudoinverso sugiere 2 PC. (Y mirando cualquier diagrama de dispersión de PCA para el conjunto de datos de Iris, parece que ambas primeras PC tienen alguna señal). Y en el caso del vino, el método pseudoinverso apunta claramente a 3 PC, mientras que el método aproximado no puede decidir entre 3 y 5.

---

Código de Matlab para realizar validación cruzada y trazar los resultados

function pca_loocv(X)

%// loop over data points 
for n=1:size(X,1)
    Xtrain = X([1:n-1 n+1:end],:);
    mu = mean(Xtrain);
    Xtrain = bsxfun(@minus, Xtrain, mu);
    [~,~,V] = svd(Xtrain, 'econ');
    Xtest = X(n,:);
    Xtest = bsxfun(@minus, Xtest, mu);

    %// loop over the number of PCs
    for j=1:min(size(V,2),25)
        P = V(:,1:j)*V(:,1:j)';        %//'
        err1 = Xtest * (eye(size(P)) - P);
        err2 = Xtest * (eye(size(P)) - P + diag(diag(P)));
        for k=1:size(Xtest,2)
            proj = Xtest(:,[1:k-1 k+1:end])*pinv(V([1:k-1 k+1:end],1:j))'*V(:,1:j)'; 
            err3(k) = Xtest(k) - proj(k);
        end

        error1(n,j) = sum(err1(:).^2);
        error2(n,j) = sum(err2(:).^2);
        error3(n,j) = sum(err3(:).^2);
    end    
end

error1 = sum(error1);
error2 = sum(error2);
error3 = sum(error3);
%// plotting code
figure
hold on
plot(error1, 'k.--')
plot(error2, 'r.-')
plot(error3, 'b.-')
legend({'Naive method', 'Approximate method', 'Pseudoinverse method'}, ...
    'Location', 'NorthWest')
legend boxoff
set(gca, 'XTick', 1:length(error1))
set(gca, 'YTick', [])
xlabel('Number of PCs')
ylabel('Cross-validation error')

Para agregar un punto aún más general a la buena respuesta de @ amoeba:

Una diferencia práctica y crucial entre los modelos supervisados ​​y no supervisados ​​es que, para los modelos no supervisados, debe pensar mucho más en qué consideraría equivalente y qué no.

Los modelos supervisados ​​siempre tienen su salida final de una manera en la que no necesita preocuparse demasiado por esto: por definición y construcción, afirma tener el mismo significado que , por lo que puede compararlo directamente.$\hat y$$\hat y$$y$

Para construir medidas de rendimiento significativas, debe pensar qué tipos de libertad del modelo no tienen sentido para su aplicación y cuáles no. Eso llevaría a una PRENSA en los puntajes, posiblemente (¿generalmente?) Después de algún tipo de rotación / volteo similar a Procrustes.

PRESIONE en x Mi conjetura es (ahora no tengo tiempo para averiguar qué hacen sus 2 líneas de código, pero ¿tal vez podría pasar por las líneas y echar un vistazo?):

Para obtener una medida que sea útil para determinar una buena complejidad del modelo a partir de una medida que proporcione una bondad de ajuste que generalmente aumentará hasta que se alcance el modelo de rango completo, debe penalizar por modelos demasiado complejos. Lo que a su vez significa que esta penalización es a) crucial yb) el ajuste de la penalización ajustará la complejidad elegida.

---

Nota al margen: me gustaría agregar que sería muy cuidadoso con este tipo de optimización automatizada de la complejidad del modelo. En mi experiencia, muchos de estos algoritmos solo producen pseudo-objetividad, y a menudo tienen el costo de funcionar bien solo para ciertos tipos de datos.