view solvers/SMALL_chol.m @ 128:8e660fd14774 ivand_dev

Feature 186

function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin) 
%% Implementation of OMP with Cholesky factorisation
% 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.
%

%
%   Centre for Digital Music, Queen Mary, University of London.
%   This file copyright 2009 Ivan Damnjanovic.
%
%   This program 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 2 of the
%   License, or (at your option) any later version.  See the file
%   COPYING included with this distribution for more information.
%   
%%
% 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=0;
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);
    errorGoal1=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