comparison solvers/SMALL_chol.m @ 1:7750624e0c73 version0.5

(none)

+function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin) 
+%%
+%=============================================
+% Sparse coding of a group of signals based on a given 
+% dictionary and specified number of atoms to use. 
+% input arguments: Dict - the dictionary
+%                   X - the signals to represent
+%                   m - number of atoms in Dictionary
+%                   errorGoal - the maximal allowed representation error for
+%                   each signal.
+%
+% optional:         if Dict is function handle then Transpose Dictionary
+%                   handle needs to be specified.
+%
+% output arguments: A - sparse coefficient matrix.
+%
+% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
+% with pinv() is changed with OMP Cholesky update)
+% Ivan Damnjanovic 2009
+%=============================================
+%%
+% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers 
+explicitD=0;
+if     isa(Dict,'float')      
+            D =@(z) Dict*z;  
+            Dt =@(z) Dict'*z;
+            explicitD=1;
+elseif isobject(Dict) 
+            D =@(z) Dict*z;  
+            Dt =@(z) Dict'*z;
+elseif isa(Dict,'function_handle') 
+    try
+        DictT=varargin{1};
+        if isa(DictT,'function_handle'); 
+                D=Dict; 
+                Dt=DictT;
+        else
+                error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
+        end
+    catch
+        error('If Dictionary is a function handle, Transpose Dictionary  needs to be specified. Exiting.');
+    end
+else
+    error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
+end
+%%
+[n,P]=size(X);
+
+
+
+global opts opts_tr machPrec
+opts.UT = true; 
+opts_tr.UT = true; opts_tr.TRANSA = true;
+machPrec = 1e-5;
+
+A = sparse(m,size(X,2));
+for k=1:1:P,
+    
+    R_I = [];
+    x=X(:,k);
+    residual=x;
+	indx = [];
+	a = zeros(m,1);
+	currResNorm = norm(residual);
+    errorGoal=errorGoal*currResNorm;
+	j = 0;
+    while currResNorm>errorGoal & j < maxNumCoef,
+		j = j+1;
+        dir=Dt(residual);
+        
+        [tmp__, pos]=max(abs(dir));
+        
+        [R_I, flag] = updateChol(R_I, n, m, D, explicitD, indx, pos, Dt);
+        
+       
+            indx(j)=pos;
+            dx=zeros(m,1);
+        
+            z = linsolve(R_I,dir(indx),opts_tr);
+        
+            dx(indx) = linsolve(R_I,z,opts);
+            a(indx) = a(indx) + dx(indx);
+                
+            residual=x-D(a);
+            currResNorm = norm(residual);
+        
+    
+   end;
+   if (~isempty(indx))
+       A(indx,k)=a(indx);
+   end
+end;
+return;
+
+
+function [R, flag] = updateChol(R, n, N, A, explicitA, activeSet, newIndex, varargin)
+
+% updateChol: Updates the Cholesky factor R of the matrix 
+% A(:,activeSet)'*A(:,activeSet) by adding A(:,newIndex)
+% If the candidate column is in the span of the existing 
+% active set, R is not updated, and flag is set to 1.
+
+global opts_tr machPrec
+flag = 0;
+
+if (explicitA)
+    newVec = A(:,newIndex);
+else
+    At=varargin{1};
+    e = zeros(N,1);
+    e(newIndex) = 1;
+    newVec = A(e);%feval(A,1,n,N,e,1:N,N); 
+end
+
+if isempty(activeSet),
+    R = sqrt(sum(newVec.^2));
+else
+    if (explicitA)
+        p = linsolve(R,A(:,activeSet)'*A(:,newIndex),opts_tr);
+    else
+        AnewVec = At(newVec);%feval(A,2,n,length(activeSet),newVec,activeSet,N);
+        p = linsolve(R,AnewVec(activeSet),opts_tr);
+    end
+    q = sum(newVec.^2) - sum(p.^2);
+    if (q <= machPrec) % Collinear vector
+        flag = 1;
+    else
+        R = [R p; zeros(1, size(R,2)) sqrt(q)];
+    end
+end