annotate solvers/SMALL_chol.m @ 247:ecce33192fcc tip

Added tag ver_2.1 for changeset cef4500b936f
author luisf <luis.figueira@eecs.qmul.ac.uk>
date Wed, 31 Oct 2012 12:24:44 +0000
parents 8e660fd14774
children
rev   line source
idamnjanovic@1 1 function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin)
ivan@128 2 %% Implementation of OMP with Cholesky factorisation
ivan@128 3 % Sparse coding of a group of signals based on a given
ivan@128 4 % dictionary and specified number of atoms to use.
ivan@128 5 % input arguments: Dict - the dictionary
ivan@128 6 % X - the signals to represent
ivan@128 7 % m - number of atoms in Dictionary
ivan@128 8 % errorGoal - the maximal allowed representation error for
ivan@128 9 % each signal.
ivan@128 10 %
ivan@128 11 % optional: if Dict is function handle then Transpose Dictionary
ivan@128 12 % handle needs to be specified.
ivan@128 13 %
ivan@128 14 % output arguments: A - sparse coefficient matrix.
ivan@128 15 %
ivan@128 16
idamnjanovic@22 17 %
idamnjanovic@22 18 % Centre for Digital Music, Queen Mary, University of London.
idamnjanovic@22 19 % This file copyright 2009 Ivan Damnjanovic.
idamnjanovic@22 20 %
idamnjanovic@22 21 % This program is free software; you can redistribute it and/or
idamnjanovic@22 22 % modify it under the terms of the GNU General Public License as
idamnjanovic@22 23 % published by the Free Software Foundation; either version 2 of the
idamnjanovic@22 24 % License, or (at your option) any later version. See the file
idamnjanovic@22 25 % COPYING included with this distribution for more information.
idamnjanovic@22 26 %
idamnjanovic@1 27 %%
idamnjanovic@1 28 % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers
idamnjanovic@1 29 explicitD=0;
idamnjanovic@1 30 if isa(Dict,'float')
idamnjanovic@1 31 D =@(z) Dict*z;
idamnjanovic@1 32 Dt =@(z) Dict'*z;
idamnjanovic@63 33 explicitD=0;
idamnjanovic@1 34 elseif isobject(Dict)
idamnjanovic@1 35 D =@(z) Dict*z;
idamnjanovic@1 36 Dt =@(z) Dict'*z;
idamnjanovic@1 37 elseif isa(Dict,'function_handle')
idamnjanovic@1 38 try
idamnjanovic@1 39 DictT=varargin{1};
idamnjanovic@1 40 if isa(DictT,'function_handle');
idamnjanovic@74 41 D=Dict;
idamnjanovic@74 42 Dt=DictT;
idamnjanovic@1 43 else
idamnjanovic@1 44 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1 45 end
idamnjanovic@1 46 catch
idamnjanovic@1 47 error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.');
idamnjanovic@1 48 end
idamnjanovic@1 49 else
idamnjanovic@1 50 error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1 51 end
idamnjanovic@1 52 %%
idamnjanovic@1 53 [n,P]=size(X);
idamnjanovic@1 54
idamnjanovic@1 55
idamnjanovic@1 56
idamnjanovic@1 57 global opts opts_tr machPrec
idamnjanovic@1 58 opts.UT = true;
idamnjanovic@1 59 opts_tr.UT = true; opts_tr.TRANSA = true;
idamnjanovic@1 60 machPrec = 1e-5;
idamnjanovic@1 61
idamnjanovic@1 62 A = sparse(m,size(X,2));
idamnjanovic@1 63 for k=1:1:P,
idamnjanovic@1 64
idamnjanovic@1 65 R_I = [];
idamnjanovic@1 66 x=X(:,k);
idamnjanovic@1 67 residual=x;
idamnjanovic@1 68 indx = [];
idamnjanovic@1 69 a = zeros(m,1);
idamnjanovic@1 70 currResNorm = norm(residual);
idamnjanovic@63 71 errorGoal1=errorGoal*currResNorm;
idamnjanovic@1 72 j = 0;
idamnjanovic@1 73 while currResNorm>errorGoal & j < maxNumCoef,
idamnjanovic@1 74 j = j+1;
idamnjanovic@1 75 dir=Dt(residual);
idamnjanovic@1 76
idamnjanovic@1 77 [tmp__, pos]=max(abs(dir));
idamnjanovic@1 78
idamnjanovic@1 79 [R_I, flag] = updateChol(R_I, n, m, D, explicitD, indx, pos, Dt);
idamnjanovic@1 80
idamnjanovic@1 81
idamnjanovic@1 82 indx(j)=pos;
idamnjanovic@1 83 dx=zeros(m,1);
idamnjanovic@1 84
idamnjanovic@1 85 z = linsolve(R_I,dir(indx),opts_tr);
idamnjanovic@1 86
idamnjanovic@1 87 dx(indx) = linsolve(R_I,z,opts);
idamnjanovic@1 88 a(indx) = a(indx) + dx(indx);
idamnjanovic@1 89
idamnjanovic@1 90 residual=x-D(a);
idamnjanovic@1 91 currResNorm = norm(residual);
idamnjanovic@1 92
idamnjanovic@1 93
idamnjanovic@1 94 end;
idamnjanovic@1 95 if (~isempty(indx))
idamnjanovic@1 96 A(indx,k)=a(indx);
idamnjanovic@1 97 end
idamnjanovic@1 98 end;
idamnjanovic@1 99 return;
idamnjanovic@1 100
idamnjanovic@1 101
idamnjanovic@1 102 function [R, flag] = updateChol(R, n, N, A, explicitA, activeSet, newIndex, varargin)
idamnjanovic@1 103
idamnjanovic@1 104 % updateChol: Updates the Cholesky factor R of the matrix
idamnjanovic@1 105 % A(:,activeSet)'*A(:,activeSet) by adding A(:,newIndex)
idamnjanovic@1 106 % If the candidate column is in the span of the existing
idamnjanovic@1 107 % active set, R is not updated, and flag is set to 1.
idamnjanovic@1 108
idamnjanovic@1 109 global opts_tr machPrec
idamnjanovic@1 110 flag = 0;
idamnjanovic@1 111
idamnjanovic@1 112 if (explicitA)
idamnjanovic@1 113 newVec = A(:,newIndex);
idamnjanovic@1 114 else
idamnjanovic@1 115 At=varargin{1};
idamnjanovic@1 116 e = zeros(N,1);
idamnjanovic@1 117 e(newIndex) = 1;
idamnjanovic@1 118 newVec = A(e);%feval(A,1,n,N,e,1:N,N);
idamnjanovic@1 119 end
idamnjanovic@1 120
idamnjanovic@1 121 if isempty(activeSet),
idamnjanovic@1 122 R = sqrt(sum(newVec.^2));
idamnjanovic@1 123 else
idamnjanovic@1 124 if (explicitA)
idamnjanovic@1 125 p = linsolve(R,A(:,activeSet)'*A(:,newIndex),opts_tr);
idamnjanovic@1 126 else
idamnjanovic@1 127 AnewVec = At(newVec);%feval(A,2,n,length(activeSet),newVec,activeSet,N);
idamnjanovic@1 128 p = linsolve(R,AnewVec(activeSet),opts_tr);
idamnjanovic@1 129 end
idamnjanovic@1 130 q = sum(newVec.^2) - sum(p.^2);
idamnjanovic@1 131 if (q <= machPrec) % Collinear vector
idamnjanovic@1 132 flag = 1;
idamnjanovic@1 133 else
idamnjanovic@1 134 R = [R p; zeros(1, size(R,2)) sqrt(q)];
idamnjanovic@1 135 end
idamnjanovic@1 136 end