Ajouter 'octave/fpica.m'
This commit is contained in:
parent
d7ca70f114
commit
ed7d6a9a40
|
|
@ -0,0 +1,902 @@
|
|||
function [A, W] = fpica(X, whiteningMatrix, dewhiteningMatrix, approach, ...
|
||||
numOfIC, g, finetune, a1, a2, myy, stabilization, ...
|
||||
epsilon, maxNumIterations, maxFinetune, initState, ...
|
||||
guess, sampleSize, displayMode, displayInterval, ...
|
||||
s_verbose, plottype);
|
||||
|
||||
%FPICA - Fixed point ICA. Main algorithm of FASTICA.
|
||||
%
|
||||
% [A, W] = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach,
|
||||
% numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon,
|
||||
% maxNumIterations, maxFinetune, initState, guess, sampleSize,
|
||||
% displayMode, displayInterval, verbose);
|
||||
%
|
||||
% Perform independent component analysis using Hyvarinen's fixed point
|
||||
% algorithm. Outputs an estimate of the mixing matrix A and its inverse W.
|
||||
%
|
||||
% whitesig :the whitened data as row vectors
|
||||
% whiteningMatrix :can be obtained with function whitenv
|
||||
% dewhiteningMatrix :can be obtained with function whitenv
|
||||
% approach [ 'symm' | 'defl' ] :the approach used (deflation or symmetric)
|
||||
% numOfIC [ 0 - Dim of whitesig ] :number of independent components estimated
|
||||
% g [ 'pow3' | 'tanh' | :the nonlinearity used
|
||||
% 'gaus' | 'skew' ]
|
||||
% finetune [same as g + 'off'] :the nonlinearity used in finetuning.
|
||||
% a1 :parameter for tuning 'tanh'
|
||||
% a2 :parameter for tuning 'gaus'
|
||||
% mu :step size in stabilized algorithm
|
||||
% stabilization [ 'on' | 'off' ] :if mu < 1 then automatically on
|
||||
% epsilon :stopping criterion
|
||||
% maxNumIterations :maximum number of iterations
|
||||
% maxFinetune :maximum number of iteretions for finetuning
|
||||
% initState [ 'rand' | 'guess' ] :initial guess or random initial state. See below
|
||||
% guess :initial guess for A. Ignored if initState = 'rand'
|
||||
% sampleSize [ 0 - 1 ] :percentage of the samples used in one iteration
|
||||
% displayMode [ 'signals' | 'basis' | :plot running estimate
|
||||
% 'filters' | 'off' ]
|
||||
% displayInterval :number of iterations we take between plots
|
||||
% verbose [ 'on' | 'off' ] :report progress in text format
|
||||
%
|
||||
% EXAMPLE
|
||||
% [E, D] = pcamat(vectors);
|
||||
% [nv, wm, dwm] = whitenv(vectors, E, D);
|
||||
% [A, W] = fpica(nv, wm, dwm);
|
||||
%
|
||||
%
|
||||
% This function is needed by FASTICA and FASTICAG
|
||||
%
|
||||
% See also FASTICA, FASTICAG, WHITENV, PCAMAT
|
||||
|
||||
% 15.1.2001
|
||||
% Hugo Gävert
|
||||
|
||||
% Added by DR
|
||||
% plottype [ 'dispsig' | 'classic' | '' | 'complot'
|
||||
% 'scatter' | 'compare' | 'sum' | 'sumerror' ]
|
||||
%
|
||||
% :mode of the output plot
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Global variable for stopping the ICA calculations from the GUI
|
||||
%global g_FastICA_interrupt=[];
|
||||
%if isempty(g_FastICA_interrupt)
|
||||
% clear global g_FastICA_interrupt;
|
||||
% interruptible = 0;
|
||||
%else
|
||||
% interruptible = 1;
|
||||
%end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Default values
|
||||
|
||||
if nargin < 3, error('Not enough arguments!'); end
|
||||
[vectorSize, numSamples] = size(X);
|
||||
if nargin < 20, s_verbose = 'on'; end
|
||||
if nargin < 19, displayInterval = 1; end
|
||||
if nargin < 18, displayMode = 'on'; end
|
||||
if nargin < 17, sampleSize = 1; end
|
||||
if nargin < 16, guess = 1; end
|
||||
if nargin < 15, initState = 'rand'; end
|
||||
if nargin < 14, maxFinetune = 100; end
|
||||
if nargin < 13, maxNumIterations = 1000; end
|
||||
if nargin < 12, epsilon = 0.0001; end
|
||||
if nargin < 11, stabilization = 'on'; end
|
||||
if nargin < 10, myy = 1; end
|
||||
if nargin < 9, a2 = 1; end
|
||||
if nargin < 8, a1 = 1; end
|
||||
if nargin < 7, finetune = 'off'; end
|
||||
if nargin < 6, g = 'pow3'; end
|
||||
if nargin < 5, numOfIC = vectorSize; end % vectorSize = Dim
|
||||
if nargin < 4, approach = 'defl'; end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the data
|
||||
|
||||
if imag(X) ~= 0
|
||||
error('Input has an imaginary part.');
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for verbose
|
||||
|
||||
str_param = lower(s_verbose);
|
||||
|
||||
if strcmp (str_param, 'on'),
|
||||
b_verbose = 1;
|
||||
elseif strcmp (str_param, 'off'),
|
||||
b_verbose = 0;
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for approach
|
||||
|
||||
str_param = lower(approach);
|
||||
if strcmp (str_param, 'symm'),
|
||||
approachMode = 1;
|
||||
elseif strcmp (str_param, 'defl'),
|
||||
approachMode = 2;
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''approach''\n', approach));
|
||||
end
|
||||
if b_verbose, fprintf('Used approach [ %s ].\n', approach); end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for numOfIC
|
||||
|
||||
if vectorSize < numOfIC
|
||||
error('Must have numOfIC <= Dimension!');
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the sampleSize
|
||||
if sampleSize > 1
|
||||
sampleSize = 1;
|
||||
if b_verbose
|
||||
fprintf('Warning: Setting ''sampleSize'' to 1.\n');
|
||||
end
|
||||
elseif sampleSize < 1
|
||||
if (sampleSize * numSamples) < 1000
|
||||
sampleSize = min(1000/numSamples, 1);
|
||||
if b_verbose
|
||||
fprintf('Warning: Setting ''sampleSize'' to %0.3f (%d samples).\n', ...
|
||||
sampleSize, floor(sampleSize * numSamples));
|
||||
end
|
||||
end
|
||||
end
|
||||
if b_verbose
|
||||
if b_verbose & (sampleSize < 1)
|
||||
fprintf('Using about %0.0f%% of the samples in random order in every step.\n',sampleSize*100);
|
||||
end
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for nonlinearity.
|
||||
|
||||
str_param = lower(g);
|
||||
if strcmp (str_param, 'pow3'),
|
||||
gOrig = 10;
|
||||
elseif strcmp (str_param, 'tanh'),
|
||||
gOrig = 20;
|
||||
elseif strcmp (str_param, 'gaus') | strcmp (str_param, 'gauss'),
|
||||
gOrig = 30;
|
||||
elseif strcmp (str_param, 'skew'),
|
||||
gOrig = 40;
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''g''\n', g));
|
||||
end
|
||||
|
||||
if sampleSize ~= 1
|
||||
gOrig = gOrig + 2;
|
||||
end
|
||||
if myy ~= 1
|
||||
gOrig = gOrig + 1;
|
||||
end
|
||||
|
||||
if b_verbose,
|
||||
fprintf('Used nonlinearity [ %s ].\n', g);
|
||||
end
|
||||
|
||||
finetuningEnabled = 1;
|
||||
str_param = lower(finetune);
|
||||
if strcmp (str_param, 'pow3'),
|
||||
gFine = 10 + 1;
|
||||
elseif strcmp(str_param, 'tanh'),
|
||||
gFine = 20 + 1;
|
||||
elseif strcmp(str_param, 'gaus'),
|
||||
gFine = 30 + 1;
|
||||
elseif strcmp(str_param, 'gauss'),
|
||||
gFine = 30 + 1;
|
||||
elseif strcmp(str_param, 'skew'),
|
||||
gFine = 40 + 1;
|
||||
elseif strcmp(str_param, 'off'),
|
||||
if (myy ~= 1),
|
||||
gFine = gOrig;
|
||||
else
|
||||
gFine = gOrig + 1;
|
||||
end
|
||||
finetuningEnabled = 0;
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''finetune''\n', ...
|
||||
finetune));
|
||||
end
|
||||
|
||||
if b_verbose & finetuningEnabled
|
||||
fprintf('Finetuning enabled (nonlinearity: [ %s ]).\n', finetune);
|
||||
end
|
||||
|
||||
str_param = lower(stabilization);
|
||||
if strcmp (str_param, 'on'),
|
||||
stabilizationEnabled = 1;
|
||||
elseif strcmp (str_param, 'off'),
|
||||
if myy ~= 1,
|
||||
stabilizationEnabled = 1;
|
||||
else
|
||||
stabilizationEnabled = 0;
|
||||
end
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''stabilization''\n', ...
|
||||
stabilization));
|
||||
end
|
||||
|
||||
if b_verbose & stabilizationEnabled
|
||||
fprintf('Using stabilized algorithm.\n');
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Some other parameters
|
||||
myyOrig = myy;
|
||||
% When we start fine-tuning we'll set myy = myyK * myy
|
||||
myyK = 0.01;
|
||||
% How many times do we try for convergence until we give up.
|
||||
failureLimit = 5;
|
||||
|
||||
|
||||
usedNlinearity = gOrig;
|
||||
stroke = 0;
|
||||
notFine = 1;
|
||||
long = 0;
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for initial state.
|
||||
|
||||
str_param = lower(initState);
|
||||
if strcmp (str_param, 'rand'),
|
||||
initialStateMode = 0;
|
||||
elseif strcmp (str_param,'guess'),
|
||||
if size(guess,1) ~= size(whiteningMatrix,2),
|
||||
initialStateMode = 0;
|
||||
if b_verbose
|
||||
fprintf('Warning: size of initial guess is incorrect. Using random initial guess.\n');
|
||||
end
|
||||
else
|
||||
initialStateMode = 1;
|
||||
if size(guess,2) < numOfIC
|
||||
if b_verbose
|
||||
fprintf('Warning: initial guess only for first %d components. Using random initial guess for others.\n', size(guess,2));
|
||||
end
|
||||
guess(:, size(guess, 2) + 1:numOfIC) = ...
|
||||
rand(vectorSize,numOfIC-size(guess,2))-.5;
|
||||
elseif size(guess,2)>numOfIC
|
||||
guess=guess(:,1:numOfIC);
|
||||
fprintf('Warning: Initial guess too large. The excess column are dropped.\n');
|
||||
end
|
||||
if b_verbose, fprintf('Using initial guess.\n'); end
|
||||
end
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''initState''\n', initState));
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Checking the value for display mode.
|
||||
|
||||
str_param = lower(displayMode);
|
||||
if strcmp (str_param, 'off') | strcmp (str_param, 'none'),
|
||||
usedDisplay = 0;
|
||||
elseif strcmp (str_param, 'on') | strcmp (str_param, 'signals');
|
||||
usedDisplay = 1;
|
||||
if (b_verbose & (numSamples > 10000))
|
||||
fprintf('Warning: Data vectors are very long. Plotting may take long time.\n');
|
||||
end
|
||||
if (b_verbose & (numOfIC > 25))
|
||||
fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
|
||||
end
|
||||
elseif strcmp (str_param, 'basis'),
|
||||
usedDisplay = 2;
|
||||
if (b_verbose & (numOfIC > 25))
|
||||
fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
|
||||
end
|
||||
elseif strcmp (str_param, 'filters'),
|
||||
usedDisplay = 3;
|
||||
if (b_verbose & (vectorSize > 25))
|
||||
fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
|
||||
end
|
||||
else
|
||||
error(sprintf('Illegal value [ %s ] for parameter: ''displayMode''\n', displayMode));
|
||||
end
|
||||
|
||||
% The displayInterval can't be less than 1...
|
||||
if displayInterval < 1
|
||||
displayInterval = 1;
|
||||
end
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
if b_verbose, fprintf('Starting ICA calculation...\n'); end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% SYMMETRIC APPROACH
|
||||
if approachMode == 1,
|
||||
|
||||
% set some parameters more...
|
||||
usedNlinearity = gOrig;
|
||||
stroke = 0;
|
||||
notFine = 1;
|
||||
long = 0;
|
||||
|
||||
A = zeros(vectorSize, numOfIC); % Dewhitened basis vectors.
|
||||
if initialStateMode == 0
|
||||
% Take random orthonormal initial vectors.
|
||||
B = orth(rand(vectorSize, numOfIC) - .5);
|
||||
elseif initialStateMode == 1
|
||||
% Use the given initial vector as the initial state
|
||||
B = whiteningMatrix * guess;
|
||||
end
|
||||
|
||||
BOld = zeros(size(B));
|
||||
BOld2 = zeros(size(B));
|
||||
|
||||
% This is the actual fixed-point iteration loop.
|
||||
for round = 1:maxNumIterations + 1,
|
||||
if round == maxNumIterations + 1,
|
||||
if b_verbose
|
||||
fprintf('No convergence after %d steps\n', maxNumIterations);
|
||||
fprintf('Note that the plots are probably wrong.\n');
|
||||
end
|
||||
A=[];
|
||||
W=[];
|
||||
return;
|
||||
end
|
||||
|
||||
%Aapo removed this:
|
||||
% if (interruptible & g_FastICA_interrupt)
|
||||
% if b_verbose
|
||||
% fprintf('\n\nCalculation interrupted by the user\n');
|
||||
% end
|
||||
% if ~isempty(B)
|
||||
% W = B' * whiteningMatrix;
|
||||
% A = dewhiteningMatrix * B;
|
||||
% else
|
||||
% W = [];
|
||||
%% A = [];
|
||||
% end
|
||||
% return;
|
||||
% end
|
||||
|
||||
|
||||
% Symmetric orthogonalization.
|
||||
B = B * real(inv(B' * B)^(1/2));
|
||||
|
||||
% Test for termination condition. Note that we consider opposite
|
||||
% directions here as well.
|
||||
minAbsCos = min(abs(diag(B' * BOld)));
|
||||
minAbsCos2 = min(abs(diag(B' * BOld2)));
|
||||
|
||||
if (1 - minAbsCos < epsilon)
|
||||
if finetuningEnabled & notFine
|
||||
if b_verbose, fprintf('Initial convergence, fine-tuning: \n'); end;
|
||||
notFine = 0;
|
||||
usedNlinearity = gFine;
|
||||
myy = myyK * myyOrig;
|
||||
BOld = zeros(size(B));
|
||||
BOld2 = zeros(size(B));
|
||||
|
||||
else
|
||||
if b_verbose, fprintf('Convergence after %d steps\n', round); end
|
||||
|
||||
% Calculate the de-whitened vectors.
|
||||
A = dewhiteningMatrix * B;
|
||||
break;
|
||||
end
|
||||
elseif stabilizationEnabled
|
||||
if (~stroke) & (1 - minAbsCos2 < epsilon)
|
||||
if b_verbose, fprintf('Stroke!\n'); end;
|
||||
stroke = myy;
|
||||
myy = .5*myy;
|
||||
% "mod" doesn't exist in Octave, use this function explicitely:
|
||||
if usedNlinearity-2.*floor(usedNlinearity./2) == 0
|
||||
usedNlinearity = usedNlinearity + 1;
|
||||
end
|
||||
elseif stroke
|
||||
myy = stroke;
|
||||
stroke = 0;
|
||||
% "mod" doesn't exist in Octave, use this function explicitely:
|
||||
if (myy == 1) & (usedNlinearity-2.*floor(usedNlinearity./2) ~= 0)
|
||||
usedNlinearity = usedNlinearity - 1;
|
||||
end
|
||||
elseif (~long) & (round>maxNumIterations/2)
|
||||
if b_verbose, fprintf('Taking long (reducing step size)\n'); end;
|
||||
long = 1;
|
||||
myy = .5*myy;
|
||||
if usedNlinearity-2.*floor(usedNlinearity./2) == 0
|
||||
usedNlinearity = usedNlinearity + 1;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
BOld2 = BOld;
|
||||
BOld = B;
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Show the progress...
|
||||
if b_verbose
|
||||
if round == 1
|
||||
fprintf('Step no. %d\n', round);
|
||||
else
|
||||
fprintf('Step no. %d, change in value of estimate: %.3f \n', round, 1 - minAbsCos);
|
||||
end
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Also plot the current state...
|
||||
switch usedDisplay
|
||||
case 1
|
||||
if rem(round, displayInterval) == 0,
|
||||
% There was and may still be other displaymodes...
|
||||
% 1D signals
|
||||
icaplot(plottype,(X'*B)');
|
||||
end
|
||||
case 2
|
||||
if rem(round, displayInterval) == 0,
|
||||
% ... and now there are :-)
|
||||
% 1D basis
|
||||
A = dewhiteningMatrix * B;
|
||||
icaplot(plottype,A');
|
||||
end
|
||||
case 3
|
||||
if rem(round, displayInterval) == 0,
|
||||
% ... and now there are :-)
|
||||
% 1D filters
|
||||
W = B' * whiteningMatrix;
|
||||
icaplot(plottype,W);
|
||||
end
|
||||
otherwise
|
||||
end
|
||||
|
||||
%drawnow;
|
||||
|
||||
switch usedNlinearity
|
||||
% pow3
|
||||
case 10
|
||||
B = (X * (( X' * B) .^ 3)) / numSamples - 3 * B;
|
||||
case 11
|
||||
% optimoitu - epsilonin kokoisia eroja
|
||||
% tämä on optimoitu koodi, katso vanha koodi esim.
|
||||
% aikaisemmista versioista kuten 2.0 beta3
|
||||
Y = X' * B;
|
||||
Gpow3 = Y .^ 3;
|
||||
Beta = sum(Y .* Gpow3);
|
||||
D = diag(1 ./ (Beta - 3 * numSamples));
|
||||
B = B + myy * B * (Y' * Gpow3 - diag(Beta)) * D;
|
||||
case 12
|
||||
Xsub=X(:, getSamples(numSamples, sampleSize));
|
||||
B = (Xsub * (( Xsub' * B) .^ 3)) / size(Xsub,2) - 3 * B;
|
||||
case 13
|
||||
% Optimoitu
|
||||
Ysub=X(:, getSamples(numSamples, sampleSize))' * B;
|
||||
Gpow3 = Ysub .^ 3;
|
||||
Beta = sum(Ysub .* Gpow3);
|
||||
D = diag(1 ./ (Beta - 3 * size(Ysub', 2)));
|
||||
B = B + myy * B * (Ysub' * Gpow3 - diag(Beta)) * D;
|
||||
|
||||
% tanh
|
||||
case 20
|
||||
hypTan = tanh(a1 * X' * B);
|
||||
B = X * hypTan / numSamples - ...
|
||||
ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / numSamples * ...
|
||||
a1;
|
||||
case 21
|
||||
% optimoitu - epsilonin kokoisia
|
||||
Y = X' * B;
|
||||
hypTan = tanh(a1 * Y);
|
||||
Beta = sum(Y .* hypTan);
|
||||
D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
|
||||
B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
|
||||
case 22
|
||||
Xsub=X(:, getSamples(numSamples, sampleSize));
|
||||
hypTan = tanh(a1 * Xsub' * B);
|
||||
B = Xsub * hypTan / size(Xsub, 2) - ...
|
||||
ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / size(Xsub, 2) * a1;
|
||||
case 23
|
||||
% Optimoitu
|
||||
Y = X(:, getSamples(numSamples, sampleSize))' * B;
|
||||
hypTan = tanh(a1 * Y);
|
||||
Beta = sum(Y .* hypTan);
|
||||
D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
|
||||
B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
|
||||
|
||||
% gauss
|
||||
case 30
|
||||
U = X' * B;
|
||||
Usquared=U .^ 2;
|
||||
ex = exp(-a2 * Usquared / 2);
|
||||
gauss = U .* ex;
|
||||
dGauss = (1 - a2 * Usquared) .*ex;
|
||||
B = X * gauss / numSamples - ...
|
||||
ones(size(B,1),1) * sum(dGauss)...
|
||||
.* B / numSamples ;
|
||||
case 31
|
||||
% optimoitu
|
||||
Y = X' * B;
|
||||
ex = exp(-a2 * (Y .^ 2) / 2);
|
||||
gauss = Y .* ex;
|
||||
Beta = sum(Y .* gauss);
|
||||
D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
|
||||
B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
|
||||
case 32
|
||||
Xsub=X(:, getSamples(numSamples, sampleSize));
|
||||
U = Xsub' * B;
|
||||
Usquared=U .^ 2;
|
||||
ex = exp(-a2 * Usquared / 2);
|
||||
gauss = U .* ex;
|
||||
dGauss = (1 - a2 * Usquared) .*ex;
|
||||
B = Xsub * gauss / size(Xsub,2) - ...
|
||||
ones(size(B,1),1) * sum(dGauss)...
|
||||
.* B / size(Xsub,2) ;
|
||||
case 33
|
||||
% Optimoitu
|
||||
Y = X(:, getSamples(numSamples, sampleSize))' * B;
|
||||
ex = exp(-a2 * (Y .^ 2) / 2);
|
||||
gauss = Y .* ex;
|
||||
Beta = sum(Y .* gauss);
|
||||
D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
|
||||
B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
|
||||
|
||||
% skew
|
||||
case 40
|
||||
B = (X * ((X' * B) .^ 2)) / numSamples;
|
||||
case 41
|
||||
% Optimoitu
|
||||
Y = X' * B;
|
||||
Gskew = Y .^ 2;
|
||||
Beta = sum(Y .* Gskew);
|
||||
D = diag(1 ./ (Beta));
|
||||
B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
|
||||
case 42
|
||||
Xsub=X(:, getSamples(numSamples, sampleSize));
|
||||
B = (Xsub * ((Xsub' * B) .^ 2)) / size(Xsub,2);
|
||||
case 43
|
||||
% Uusi optimoitu
|
||||
Y = X(:, getSamples(numSamples, sampleSize))' * B;
|
||||
Gskew = Y .^ 2;
|
||||
Beta = sum(Y .* Gskew);
|
||||
D = diag(1 ./ (Beta));
|
||||
B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
|
||||
|
||||
otherwise
|
||||
error('Code for desired nonlinearity not found!');
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
% Calculate ICA filters.
|
||||
W = B' * whiteningMatrix;
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Also plot the last one...
|
||||
switch usedDisplay
|
||||
case 1
|
||||
% There was and may still be other displaymodes...
|
||||
% 1D signals
|
||||
icaplot(plottype,(X'*B)');
|
||||
%drawnow;
|
||||
case 2
|
||||
% ... and now there are :-)
|
||||
% 1D basis
|
||||
icaplot(plottype,A');
|
||||
%drawnow;
|
||||
case 3
|
||||
% ... and now there are :-)
|
||||
% 1D filters
|
||||
icaplot(plottype,W);
|
||||
%drawnow;
|
||||
otherwise
|
||||
end
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% DEFLATION APPROACH
|
||||
if approachMode == 2
|
||||
|
||||
B = zeros(vectorSize);
|
||||
|
||||
% The search for a basis vector is repeated numOfIC times.
|
||||
round = 1;
|
||||
|
||||
numFailures = 0;
|
||||
|
||||
while round <= numOfIC,
|
||||
myy = myyOrig;
|
||||
usedNlinearity = gOrig;
|
||||
stroke = 0;
|
||||
notFine = 1;
|
||||
long = 0;
|
||||
endFinetuning = 0;
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Show the progress...
|
||||
if b_verbose, fprintf('IC %d ', round); end
|
||||
|
||||
% Take a random initial vector of lenght 1 and orthogonalize it
|
||||
% with respect to the other vectors.
|
||||
if initialStateMode == 0
|
||||
w = rand(vectorSize, 1) - .5;
|
||||
elseif initialStateMode == 1
|
||||
w=whiteningMatrix*guess(:,round);
|
||||
end
|
||||
w = w - B * B' * w;
|
||||
w = w / norm(w);
|
||||
|
||||
wOld = zeros(size(w));
|
||||
wOld2 = zeros(size(w));
|
||||
|
||||
% This is the actual fixed-point iteration loop.
|
||||
% for i = 1 : maxNumIterations + 1
|
||||
i = 1;
|
||||
gabba = 1;
|
||||
while i <= maxNumIterations + gabba
|
||||
%drawnow;
|
||||
|
||||
% Aapo removed::
|
||||
% %%DR also removed from next line " & g_FastICA_interrupt
|
||||
% if (interruptible)
|
||||
% if b_verbose
|
||||
% fprintf('\n\nCalculation interrupted by the user\n');
|
||||
% end
|
||||
% return;
|
||||
% end
|
||||
|
||||
% Project the vector into the space orthogonal to the space
|
||||
% spanned by the earlier found basis vectors. Note that we can do
|
||||
% the projection with matrix B, since the zero entries do not
|
||||
% contribute to the projection.
|
||||
w = w - B * B' * w;
|
||||
w = w / norm(w);
|
||||
|
||||
if notFine
|
||||
if i == maxNumIterations + 1
|
||||
if b_verbose
|
||||
fprintf('\nComponent number %d did not converge in %d iterations.\n', round, maxNumIterations);
|
||||
end
|
||||
round = round - 1;
|
||||
numFailures = numFailures + 1;
|
||||
if numFailures > failureLimit
|
||||
if b_verbose
|
||||
fprintf('Too many failures to converge (%d). Giving up.\n', numFailures);
|
||||
end
|
||||
if round == 0
|
||||
A=[];
|
||||
W=[];
|
||||
end
|
||||
return;
|
||||
end
|
||||
% numFailures > failurelimit
|
||||
break;
|
||||
end
|
||||
% i == maxNumIterations + 1
|
||||
else
|
||||
% if notFine
|
||||
if i >= endFinetuning
|
||||
wOld = w; % So the algorithm will stop on the next test...
|
||||
end
|
||||
end
|
||||
% if notFine
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Show the progress...
|
||||
if b_verbose, fprintf('.'); fflush(stdout); end;
|
||||
|
||||
|
||||
% Test for termination condition. Note that the algorithm has
|
||||
% converged if the direction of w and wOld is the same, this
|
||||
% is why we test the two cases.
|
||||
if norm(w - wOld) < epsilon | norm(w + wOld) < epsilon
|
||||
if finetuningEnabled & notFine
|
||||
if b_verbose, fprintf('Initial convergence, fine-tuning: '); end;
|
||||
notFine = 0;
|
||||
gabba = maxFinetune;
|
||||
wOld = zeros(size(w));
|
||||
wOld2 = zeros(size(w));
|
||||
usedNlinearity = gFine;
|
||||
myy = myyK * myyOrig;
|
||||
|
||||
endFinetuning = maxFinetune + i;
|
||||
|
||||
else
|
||||
numFailures = 0;
|
||||
% Save the vector
|
||||
B(:, round) = w;
|
||||
|
||||
% Calculate the de-whitened vector.
|
||||
A(:,round) = dewhiteningMatrix * w;
|
||||
% Calculate ICA filter.
|
||||
W(round,:) = w' * whiteningMatrix;
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Show the progress...
|
||||
if b_verbose, fprintf('computed ( %d steps ) \n', i); end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Also plot the current state...
|
||||
switch usedDisplay
|
||||
case 1
|
||||
if rem(round, displayInterval) == 0,
|
||||
% There was and may still be other displaymodes...
|
||||
% 1D signals
|
||||
temp = X'*B;
|
||||
icaplot(plottype,temp(:,1:numOfIC)');
|
||||
%drawnow;
|
||||
end
|
||||
case 2
|
||||
if rem(round, displayInterval) == 0,
|
||||
% ... and now there are :-)
|
||||
% 1D basis
|
||||
icaplot(plottype,A');
|
||||
%drawnow;
|
||||
end
|
||||
case 3
|
||||
if rem(round, displayInterval) == 0,
|
||||
% ... and now there are :-)
|
||||
% 1D filters
|
||||
icaplot(plottype,W);
|
||||
%drawnow;
|
||||
end
|
||||
end
|
||||
% switch usedDisplay
|
||||
break; % IC ready - next...
|
||||
end
|
||||
%if finetuningEnabled & notFine
|
||||
elseif stabilizationEnabled
|
||||
if (~stroke) & (norm(w - wOld2) < epsilon | norm(w + wOld2) < ...
|
||||
epsilon)
|
||||
stroke = myy;
|
||||
if b_verbose, fprintf('Stroke!'); end;
|
||||
myy = .5*myy;
|
||||
if mod(usedNlinearity,2) == 0
|
||||
usedNlinearity = usedNlinearity + 1;
|
||||
end
|
||||
elseif stroke
|
||||
myy = stroke;
|
||||
stroke = 0;
|
||||
if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
|
||||
usedNlinearity = usedNlinearity - 1;
|
||||
end
|
||||
elseif (notFine) & (~long) & (i > maxNumIterations / 2)
|
||||
if b_verbose, fprintf('Taking long (reducing step size) '); end;
|
||||
long = 1;
|
||||
myy = .5*myy;
|
||||
if mod(usedNlinearity,2) == 0
|
||||
usedNlinearity = usedNlinearity + 1;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
wOld2 = wOld;
|
||||
wOld = w;
|
||||
|
||||
switch usedNlinearity
|
||||
% pow3
|
||||
case 10
|
||||
w = (X * ((X' * w) .^ 3)) / numSamples - 3 * w;
|
||||
case 11
|
||||
EXGpow3 = (X * ((X' * w) .^ 3)) / numSamples;
|
||||
Beta = w' * EXGpow3;
|
||||
w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
|
||||
case 12
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
w = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2) - 3 * w;
|
||||
case 13
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
EXGpow3 = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2);
|
||||
Beta = w' * EXGpow3;
|
||||
w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
|
||||
% tanh
|
||||
case 20
|
||||
hypTan = tanh(a1 * X' * w);
|
||||
w = (X * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / numSamples;
|
||||
case 21
|
||||
hypTan = tanh(a1 * X' * w);
|
||||
Beta = w' * X * hypTan;
|
||||
w = w - myy * ((X * hypTan - Beta * w) / ...
|
||||
(a1 * sum((1-hypTan .^2)') - Beta));
|
||||
case 22
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
hypTan = tanh(a1 * Xsub' * w);
|
||||
w = (Xsub * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / size(Xsub, 2);
|
||||
case 23
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
hypTan = tanh(a1 * Xsub' * w);
|
||||
Beta = w' * Xsub * hypTan;
|
||||
w = w - myy * ((Xsub * hypTan - Beta * w) / ...
|
||||
(a1 * sum((1-hypTan .^2)') - Beta));
|
||||
% gauss
|
||||
case 30
|
||||
% This has been split for performance reasons.
|
||||
u = X' * w;
|
||||
u2=u.^2;
|
||||
ex=exp(-a2 * u2/2);
|
||||
gauss = u.*ex;
|
||||
dGauss = (1 - a2 * u2) .*ex;
|
||||
w = (X * gauss - sum(dGauss)' * w) / numSamples;
|
||||
case 31
|
||||
u = X' * w;
|
||||
u2=u.^2;
|
||||
ex=exp(-a2 * u2/2);
|
||||
gauss = u.*ex;
|
||||
dGauss = (1 - a2 * u2) .*ex;
|
||||
Beta = w' * X * gauss;
|
||||
w = w - myy * ((X * gauss - Beta * w) / ...
|
||||
(sum(dGauss)' - Beta));
|
||||
case 32
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
u = Xsub' * w;
|
||||
u2=u.^2;
|
||||
ex=exp(-a2 * u2/2);
|
||||
gauss = u.*ex;
|
||||
dGauss = (1 - a2 * u2) .*ex;
|
||||
w = (Xsub * gauss - sum(dGauss)' * w) / size(Xsub, 2);
|
||||
case 33
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
u = Xsub' * w;
|
||||
u2=u.^2;
|
||||
ex=exp(-a2 * u2/2);
|
||||
gauss = u.*ex;
|
||||
dGauss = (1 - a2 * u2) .*ex;
|
||||
Beta = w' * Xsub * gauss;
|
||||
w = w - myy * ((Xsub * gauss - Beta * w) / ...
|
||||
(sum(dGauss)' - Beta));
|
||||
% skew
|
||||
case 40
|
||||
w = (X * ((X' * w) .^ 2)) / numSamples;
|
||||
case 41
|
||||
EXGskew = (X * ((X' * w) .^ 2)) / numSamples;
|
||||
Beta = w' * EXGskew;
|
||||
w = w - myy * (EXGskew - Beta*w)/(-Beta);
|
||||
case 42
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
w = (Xub * ((Xub' * w) .^ 2)) / size(Xsub, 2);
|
||||
case 43
|
||||
Xsub=X(:,getSamples(numSamples, sampleSize));
|
||||
EXGskew = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
|
||||
Beta = w' * EXGskew;
|
||||
w = w - myy * (EXGskew - Beta*w)/(-Beta);
|
||||
|
||||
otherwise
|
||||
error('Code for desired nonlinearity not found!');
|
||||
end
|
||||
|
||||
% Normalize the new w.
|
||||
w = w / norm(w);
|
||||
i = i + 1;
|
||||
end
|
||||
round = round + 1;
|
||||
end
|
||||
if b_verbose, fprintf('Done.\n'); end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Also plot the ones that may not have been plotted.
|
||||
if (usedDisplay > 0) & (rem(round-1, displayInterval) ~= 0)
|
||||
switch usedDisplay
|
||||
case 1
|
||||
% There was and may still be other displaymodes...
|
||||
% 1D signals
|
||||
temp = X'*B;
|
||||
icaplot(plottype,temp(:,1:numOfIC)');
|
||||
drawnow;
|
||||
case 2
|
||||
% ... and now there are :-)
|
||||
% 1D basis
|
||||
icaplot(plottype,A');
|
||||
%drawnow;
|
||||
case 3
|
||||
% ... and now there are :-)
|
||||
% 1D filters
|
||||
icaplot(plottype,W);
|
||||
%drawnow;
|
||||
otherwise
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% In the end let's check the data for some security
|
||||
if imag(A) ~= 0
|
||||
if b_verbose, fprintf('Warning: removing the imaginary part from the result.\n'); end
|
||||
A = real(A);
|
||||
W = real(W);
|
||||
end
|
||||
|
||||
endfunction
|
||||
|
||||
Loading…
Reference in New Issue