invdet.m

Feature Selection using matlab

function [InvM, detM, RankM]= InvDet(M);

            n = size(M,1);      maxpivot = 0;          npivot = 0;    r= 0;
           eps = 2.2 * 10e-10;
           
           if sum(sum(isnan(M))) > 0
                    RankM = 0;
                    InvM = zeros(n,n);
                    detM = 0;
                    return
           end
              
           
        InvM = zeros(n,n);   
        B = M;                           
		A = zeros(n, 2*n);
        A(1:n,1:n) = B;
        for i=1:n, A(i,i+n) = 1; end
        
		for k=1:n
			maxpivot=abs(A(k,k));
            if maxpivot == 0, 
                maxpivot = eps;
            end
			npivot=k;
			for i=k:n
                   if maxpivot == 0, 
                        maxpivot = eps;
                   end
                         
                   if abs(A(i,k)) == 0, 
                        A(i,k) = eps;
                   end
                   
				if maxpivot < abs(A(i,k))
         		    maxpivot=abs(A(i,k));
					npivot=i;
                end
            end
            
			
				if npivot~=k
				   for j=k:2*n
					    temp = A(npivot,j);
					    A(npivot,j)=A(k,j);
					    A(k,j) = temp;	
                    end
              	   r=r+1;
				end
								
				for i=1:n
				     if i > k 
                         
                         if A(k,k) == 0
                              disp('Singular CovMatrix');
                              [ InvM, detM, RankM ] = UseOnlyNonZeroEigValues(M,n);
                                 return
                         end
                         
                         
					   D  = - A(i,k)/ A(k,k);
						  for j=k:2*n
							A(i,j) = A(i,j) + D*A(k,j) ;
                          end
                     end
                end

     end
           
		   detM = 1;
           for k=1:n
			  D = A(k,k);
		       detM = detM * A(k,k);
	    		for j=k:2*n, 	A(k,j)=A(k,j) / D ; end
           end
              detM = abs(detM); % Rem: Why ?
             
           for k=n:-1:2
	    	   for i=1: k-1
				 D = A(i,k);
					for j=(i+1): 2*n,
                        A(i,j) = A(i,j) - D*A(k,j);
                    end   
               end
           end

          InvM(1:n, 1:n) = A(1:n, (n+1):(2*n) );
          RankM = n;

          if abs(detM) < 10^(-30) || isnan(InvM(1,1))
                      [ InvM, detM, RankM ] = UseOnlyNonZeroEigValues(M,n);
          end
          
          
%-------------------------------------------------------- 
function  [ InvM, detM, RankM ] = UseOnlyNonZeroEigValues(M,n)


                [EigVec, EigVal] = eig(M); 
                EigVal = diag(EigVal);
                IndxNonZeroEigVals = find(abs(EigVal) > 1e-4);
  
                if isempty( IndxNonZeroEigVals)
                        RankM = 0;
                        InvM = zeros(n,n);
                        detM = 0;
                        return
                end
                
                LIndxNonZeroEigVals  = length( IndxNonZeroEigVals );
                InvEigCovs = zeros(n,n);
  
                % Inverse = Vectors' * Values^(-1) * Vectors
                
                
                for IndexEigVals  =  IndxNonZeroEigVals'
                        
                    InvEigCovs = InvEigCovs + ...
                        1 / EigVal(IndexEigVals)  * ...
                                          (EigVec(:, IndexEigVals ) * EigVec(:, IndexEigVals )' );
                end
  
                    detM      = prod( EigVal(   IndxNonZeroEigVals ) );
                    RankM = LIndxNonZeroEigVals;
                    InvM     = InvEigCovs;
return