Skip to content

Problems with the FastICA filter

Description

There is a FastICAFilter in the DimensionalityReduction module. The filter is not tested (there is actually a test that is compiled but not executed). The DimensionalityReduction application uses this filter, but it also is untested.

In MR !336 (merged) we tried to add this test to CTest, but it produces different results depending on the number of threads used, which makes Dashboard testing hard. One solution would be to force the number of threads to 1 in the test, however, I don't think the result of FastICA should depend on the number of thread used and that this filter does not do what it is supposed to do.

Note that the result of the test also changes if we use float instead of double as pixel type.

The FastICA algorithm used is the one described in [1]. It is implemented as an ImageToImageFilter and uses the same framework as the other Dimensionality Reduction filters (pca, mnf ...) i.e. :

  • They can work both in "Forward" and "Backward" mode, forward meaning applying the reduction transformation and backward applying its inverse. This is a template parameter.

  • The GenerateData() method of the filter only apply the transformation, i.e. the result of the preprocessing PCA (whitening) followed by the result of the ICA. this transformation can be supplied to the filter via Setters, but if it is not, it will be computed in the UpdateOutputInformation() method using the FastICA algorithm (the filter also computes the PCA whitening transform if needed). This means that most of the processing of this filter is done in UpdateOutputInformation().

  • The FastICA filter uses (in the UpdateOutputInformation() method!) a FastICAInternalOptimizerVectorImageFilter which is a "kind of" Persistent filter (meaning it is not persistent but maybe it should be), this filter takes as input the output of the PCAFilter, and computes output pixels but also accumulators and is plugged to a StreamingStatisticsVectorImageFilter. This pipeline is executed at each iteration of the FastICA algorithm.

I found many problems in this filter, but nothing that really explains why the result depends on the number of thread :

  • The "persistent computation" of the FastICAInternalOptimizerVectorImageFilter i.e. computing quantities (m_Den and m_Beta) based on the accumulators is done in the AfterThreadedGenerateData() method. AfterThreadedGenerateData() is called every time GenerateData() is called, i.e. once for each piece of stream. So instead of computing m_Den and m_Beta using the whole image, it is computed using only the last stream of the input image. This means that the filter does not support streaming ... Note that it supports threading, and this does not explain the bug.

  • I think the formula used for the transformation matrix update step is wrong. Currently it is :

W(band, bd) -= m_Mu * ( estimator->GetMean()[bd]
                              - optimizer->GetBeta() * W(band, bd) / optimizer->GetDen() );

but according to [1] it should be :

W(band, bd) -= m_Mu * ( estimator->GetMean()[bd]
                              - optimizer->GetBeta() * W(band, bd) ) / optimizer->GetDen() ;

(note the parenthesis)

  • At the first step of the FastICA algorithm, we use the input Image as input of the filter instead of the PCA output. This is a problem because FastICA expects whitened input (uncorrelated, zero mean and unit variance).

  • The PCA is computed with the option useNormalization to true, The input image is normalized in mean and variance before the PCA. While removing the mean is required, I'm not sure about the variance normalization before the PCA for ICA in general and it doesn't work in OTB case. I'm not sure what the PCAFilter computes when the variance normalization is on, but it is not a whitening transformation (i.e. D^-1/2 * E, where D is the eigenvalue matrix and E the (transposed) eigenvector matrix). If we set UseVarianceForNormalization to false the data is correctly whitened.

  • The non linear function f which is stored as a function pointer can be set in the filter by the user. By default it is std::tanh. The problem is that we also need its derivative, which is hard-coded as 1-f^2, so the filter will only work for non linear function that verifies f'=1-f^2, i.e tanh.

  • The dimensionality reduction is done at the end of the algorithm, by selecting the five first independent components. But these components are in no specific order, so this is equivalent to taking five components at random (ICA is not really a dimensionality reduction algorithm ...). Instead the reduction should be done in the PCA step. This would also reduce computation time!

[1] A. Hyvarinen. Fast and robust fixed-point algorithms for independent component analysis. Neural Networks, IEEE Transactions on, 10(3):626–634, 1999

Edited by Cédric Traizet