How to retro-engineer hyperPCA and obtain the components from the coefficients and original input?

Dear all, I am using the hyperpca matlab function to obtain the principal components for some hyperspectral images. The function is called like this:

[components,coefficients,variance]=hyperpca(hypercube,numcoefficients);

Where hypercube is the hyperspectral image passed as an m*n*p int16 matrix.

If I understand correctly PCA is a function reprojecting the input data (hypercube) into the components space (components) and therefore a linear function.

My expectations would be that multiplying the input data by the coefficients would return the principal components. I proceeded as follows:

reshaped_inputs=reshape(hypercube,size(hypercube,1)*size(hypercube,2),size(hypercube,3));

reshaped_new_components=single(reshaped_inputs)*coefficients; % need to cast the inputs from int16 to single otherwise mtimes won't work

new_components=reshape(reshaped_new_components,size(hypercube,1),size(hypercube,2),size(coefficients,2));

I would then expect that new_components is == components, but instead the difference between those two matrices is far from zero.

Why is that? Is there some sort of scaling/centering that I'm missing?

Thanks for any inputs!

Answer 1

[Mahesh Belagal Sarvisetty]

By default, hyperpca mean-centering the input data in order to compute principal components. So, you can set MeanCentered value as false to match new_components with the components from the hyperpca function or you can mean center the input data before multiplying with the coefficients.

Here is quick snippet with MeanCentered value false

% Read Hypercube

hCube = hypercube("paviaU.dat");

% Perform PCA

[princComp, coeff] = hyperpca(hCube,20,"MeanCentered",false);

% Extract and reshape the DataCube

dataCube = (hCube.DataCube);

sz = size(dataCube);

dataCube = reshape(dataCube,[sz(1)*sz(2) sz(3)]);

% Generate the new principal component

newPrincComp = dataCube * coeff;

sz = size(princComp);

newPrincComp = reshape(newPrincComp,[sz(1) sz(2) sz(3)]);

% Calculate max absolute difference

diff = imabsdiff(princComp,newPrincComp);

max(diff(:))

[Federico]

Dear @Mahesh Belagal Sarvisetty, thank you for your answer!

I followed your example, applied it to my datacube, but, unfortunately, there's still something off when I attempt it.

The maximum difference I see between the original principal components and the back-calculated ones is quite huge (2.0621e+04).

I don't know what could be the cause, my hypothesis are:

• Could too many zeroes on the image border affect the PCA computation? The Pavia image had some zero pixels, but not that many...
• I have to convert my datacube from int16 to single. I've checked some numbers manually (since many functions, such as the simple A-B do not work on int16 data) and I can't find any difference between the int16 and the single version of the datacube.
• I'm asking something that's not really possible. For example, the same toolbox has the inverseProjection function which should reconstruct the full hyperspectral datacube from the PCA components and coefficients. Even on the Indian Pines example data the minimum difference between the reconstructed cube and the old datacube is >960

% MathWorks Example

hcube=hypercube('indian_pines.dat');

[pcDataCube,coeff]=hyperpca(hcube,10);

reconstructedData=inverseProjection(pcDataCube,coeff);

% Difference calculation

diff=imabsdiff(hcube.DataCube,reconstructedData);

disp(['Minimum difference: ',sprintf('%.2f',min(diff(:)))]);

disp(['Average difference: ',sprintf('%.2f',mean(diff(:)))]);

disp(['Maximum difference: ',sprintf('%.2f',max(diff(:)))]);

% Displayed output:

'Minimum difference: 961.11'

'Average difference: 2503.44'

'Maximum difference: 9315.66'

Any idea or comment?

[Mahesh Belagal Sarvisetty]

inverseProjection input 'coeff' contains coefficients for only 10 principal components. In order to reconstruct original data we have use coefficients of all principal components and also set the MeanCentered value as false.

Take a look at this example, reconstructedData is same as original datacube

hcube = hypercube('indian_pines.dat');
[pcDataCube,coeff] = hyperpca(hcube,220,"MeanCentered",false);
reconstructedData = inverseProjection(pcDataCube,coeff);

% Difference calculation
diff = imabsdiff(hcube.DataCube,reconstructedData);
disp(['Minimum difference: ',sprintf('%.2f',min(diff(:)))]);
disp(['Average difference: ',sprintf('%.2f',mean(diff(:)))]);
disp(['Maximum difference: ',sprintf('%.2f',max(diff(:)))]);