Zoltan Szabo avatar Zoltan Szabo committed 5da18f8

Jacobi optimization based ICA solution with general entropy/mutual information estimators: added; see 'estimate_ICA.m', 'estimate_ICA_Jacobi1.m'

Comments (0)

Files changed (7)

+v0.34 (Mar 22, 2013):
+-Jacobi optimization based ICA solution with general entropy/mutual information estimators: added; see 'estimate_ICA.m', 'estimate_ICA_Jacobi1.m'.
+
 v0.33 (Mar 6, 2013):
 -Two one-dimensional Shannon entropy estimators based on the maximum entropy method: added; see 'HShannon_MaxEnt1_initialization.m', 'HShannon_MaxEnt1_estimation.m', 'HShannon_MaxEnt2_initialization.m', 'HShannon_MaxEnt2_estimation.m'.
 
 
 **Download** the latest release: 
 
-- code: [zip](https://bitbucket.org/szzoli/ite/downloads/ITE-0.33_code.zip), [tar.bz2](https://bitbucket.org/szzoli/ite/downloads/ITE-0.33_code.tar.bz2), 
-- [documentation (pdf)](https://bitbucket.org/szzoli/ite/downloads/ITE-0.33_documentation.pdf).
+- code: [zip](https://bitbucket.org/szzoli/ite/downloads/ITE-0.34_code.zip), [tar.bz2](https://bitbucket.org/szzoli/ite/downloads/ITE-0.34_code.tar.bz2), 
+- [documentation (pdf)](https://bitbucket.org/szzoli/ite/downloads/ITE-0.34_documentation.pdf).
 
 

code/IPA/demos/demo_ISA.m

                 ICA.opt_type = 'fastICA'; 
             %II:
                 %ICA.opt_type = 'EASI';
+            %III:
+                %ICA.opt_type = 'Jacobi1';
+                %ICA.cost_type = 'sumH' or 'I', examples:
+                    %i:
+                        %ICA.cost_type = 'I';
+                        %ICA.cost_name = 'KCCA'; %you can change it to any mutual information estimator
+                    %ii:
+                        %ICA.cost_type = 'sumH';
+                        %ICA.cost_name = 'Renyi_kNN_1tok'; %you can change it to any entropy estimator
         %ISA solver (clustering of the ICA elements):
             ISA.cost_type = 'sumH'; %'I','sumH', 'sum-I','Irecursive', 'Ipairwise', 'Ipairwise1d'
             ISA.cost_name = 'Renyi_kNN_1tok'; %example: ISA.cost_type = 'sumH', ISA.cost_name = 'Renyi_kNN_1tok' means that we use an entropy sum ISA formulation ('sumH'), where the entropies are Renyi entropies estimated via kNN methods ('Renyi_kNN_1tok').

code/IPA/demos/estimate_ICA.m

         W_hat = W_ICA * W_whiten;
     case 'EASI'
         [e_hat,W_hat] = estimate_ICA_EASI(x,dim_reduction);
+    case 'Jacobi1'
+        %whitening: see the note at 'fastICA' above.
+            [x,W_whiten] = whiten(x,dim_reduction);
+        [e_hat,W_ICA] = estimate_ICA_Jacobi1(x,ICA);
+        W_hat = W_ICA * W_whiten;
     otherwise
         error('ICA method=?');
 end

code/IPA/optimization/cost_Jacobi_coordinate_pair.m

+function [cost] = cost_Jacobi_coordinate_pair(Y,co,cost_type)
+%Estimates the ICA cost (cost) of signal Y (2xnumber of samples) for a given cost object (co) and cost type ('I' or 'sumH').
+%
+%INPUT:
+%   Y: Y(:,t) is the t^th sample.
+%
+%Copyright (C) 2012 Zoltan Szabo ("http://nipg.inf.elte.hu/szzoli", "szzoli (at) cs (dot) elte (dot) hu")
+%
+%This file is part of the ITE (Information Theoretical Estimators) Matlab/Octave toolbox.
+%
+%ITE is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
+%the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+%
+%This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+%
+%You should have received a copy of the GNU General Public License along with ITE. If not, see <http://www.gnu.org/licenses/>.
+
+%verification:
+    if size(Y,1)~=2;
+        error('The size of Y is not correct, it should be: 2 x number of samples.');
+    end
+
+switch cost_type
+    case 'sumH'
+        H1 = H_estimation(Y(1,:),co);
+        H2 = H_estimation(Y(2,:),co);
+        cost = H1 + H2;
+    case 'I'
+        cost = I_estimation(Y,[1;1],co);
+    otherwise
+        error('cost type=?');
+end

code/IPA/optimization/rotation_matrix.m

+function [R] = rotation_matrix(angle)
+%Returns the 2D rotation matrix (R) corresponding to the given rotation angle: R = [cos(angle), -sin(angle); sin(angle), cos(angle)].
+%
+%Copyright (C) 2012 Zoltan Szabo ("http://nipg.inf.elte.hu/szzoli", "szzoli (at) cs (dot) elte (dot) hu")
+%
+%This file is part of the ITE (Information Theoretical Estimators) Matlab/Octave toolbox.
+%
+%ITE is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
+%the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+%
+%This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+%
+%You should have received a copy of the GNU General Public License along with ITE. If not, see <http://www.gnu.org/licenses/>.
+
+C = cos(angle); 
+S = sin(angle);
+R = [C, -S; S, C];

code/shared/estimate_ICA_Jacobi1.m

+function [e_hat,W_hat] = estimate_ICA_Jacobi1(x,ICA)
+%Performs ICA on the whitened signal 'x' using Jacobi (also known as Givens) rotation based optimization.
+%
+%INPUT:
+%   x: x(:,t) is the t^th sample.
+%   ICA: solver for independent component analysis. 
+%      Mandatory fields:
+%         ICA.cost_type: cost type ('sumH' or 'I'). 
+%         ICA.cost_name: cost name (an entropy/mutual information cost name).
+%      Examples: 
+%         ICA.cost_type = 'sumH'; ICA.cost_name = 'Renyi_kNN_1tok';
+%         ICA.cost_type = 'I'; ICA.cost_name = 'KCCA';
+%OUTPUT:
+%   e_hat: estimated ICA elements; e_hat(:,t) is the t^th sample.
+%   W_hat: estimated ICA demixing matrix.
+%
+%REFERENCE:
+%  Sergey Kirshner and Barnabas Poczos. ICA and ISA Using Schweizer-Wolff Measure of Dependence. International Conference on Machine Learning (ICML), pages 464-471, 2008. (The optimization technique of this paper has been adapted to general entropy/mutual information estimators.)
+%
+%Copyright (C) 2012 Zoltan Szabo ("http://nipg.inf.elte.hu/szzoli", "szzoli (at) cs (dot) elte (dot) hu")
+%
+%This file is part of the ITE (Information Theoretical Estimators) Matlab/Octave toolbox.
+%
+%ITE is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
+%the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+%
+%This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+%
+%You should have received a copy of the GNU General Public License along with ITE. If not, see <http://www.gnu.org/licenses/>.
+
+%initialization:
+    %initialization of the entropy/mutual information estimator:
+        %verification:
+            if ~strcmp(ICA.cost_type,'sumH') && ~strcmp(ICA.cost_type,'I')
+                error('ICA.cost type is not correct, it should be: sumH or I.');
+            end
+        co = co_initialization(ICA.cost_type,ICA.cost_name,0);%'0': multiplicative constant is not important
+    d = size(x,1); %dimension
+    %parameters you may would like to play with (you can also pass them via the ICA input structure):
+        num_of_sweeps = d;  %for fixed level the number of sweeps    
+        num_of_angles = 90; %for fixed sweep the (maximal) number of angles; the maximum is attained for the last level                
+        num_of_levels = 3; %number of levels
+    e_hat = x; %estimated source
+    W_hat = eye(d); %estimated demixing matrix
+        
+for level = 1 : num_of_levels
+    %angles for the actual level:
+        nangles = ceil(num_of_angles / 2^(num_of_levels-level));
+        angles = pi/2 * [0:nangles-1]/nangles;
+    for sweep = 1 : num_of_sweeps
+        %disp:
+            disp(strcat(['Jacobi1 optimization, level: ',num2str(level),'(/',num2str(num_of_levels),')', ' -> sweep: ',num2str(sweep),'(/',num2str(num_of_sweeps),').']));            
+        %sweep through over all (i1,i2) pairs [i1<i2]:
+            for i1 = 1 : d-1
+            for i2 = i1+1 : d
+                cost_opt = Inf;
+                %optimal angle (angle_opt):
+                    for angle = angles
+                        R = rotation_matrix(angle);
+                        %cost of 'R * e_hat([i1,i2],:))':
+                            cost = cost_Jacobi_coordinate_pair(R * e_hat([i1,i2],:),co,ICA.cost_type);
+                        if cost < cost_opt %the rotation seems to be useful
+                            cost_opt = cost;
+                            angle_opt = angle;
+                        end
+                    end
+                %update the estimated demixing matrix (W_hat), source (e_hat):
+                    R = rotation_matrix(angle_opt);%R optimal
+                    W_hat([i1,i2],:) = R * W_hat([i1,i2],:);
+                    e_hat([i1,i2],:) = R * e_hat([i1,i2],:); 
+            end
+            end
+    end
+end
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.