annotate solvers/SMALL_chol.m @ 5:f44689e95ea4

(none)
author idamnjanovic
date Mon, 22 Mar 2010 10:43:01 +0000
parents 7750624e0c73
children 524cc3fff5ac
rev   line source
idamnjanovic@1 1 function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin)
idamnjanovic@1 2 %%
idamnjanovic@1 3 %=============================================
idamnjanovic@1 4 % Sparse coding of a group of signals based on a given
idamnjanovic@1 5 % dictionary and specified number of atoms to use.
idamnjanovic@1 6 % input arguments: Dict - the dictionary
idamnjanovic@1 7 % X - the signals to represent
idamnjanovic@1 8 % m - number of atoms in Dictionary
idamnjanovic@1 9 % errorGoal - the maximal allowed representation error for
idamnjanovic@1 10 % each signal.
idamnjanovic@1 11 %
idamnjanovic@1 12 % optional: if Dict is function handle then Transpose Dictionary
idamnjanovic@1 13 % handle needs to be specified.
idamnjanovic@1 14 %
idamnjanovic@1 15 % output arguments: A - sparse coefficient matrix.
idamnjanovic@1 16 %
idamnjanovic@1 17 % based on KSVD toolbox solver found on Miki Elad webpage (finding inverse
idamnjanovic@1 18 % with pinv() is changed with OMP Cholesky update)
idamnjanovic@1 19 % Ivan Damnjanovic 2009
idamnjanovic@1 20 %=============================================
idamnjanovic@1 21 %%
idamnjanovic@1 22 % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers
idamnjanovic@1 23 explicitD=0;
idamnjanovic@1 24 if isa(Dict,'float')
idamnjanovic@1 25 D =@(z) Dict*z;
idamnjanovic@1 26 Dt =@(z) Dict'*z;
idamnjanovic@1 27 explicitD=1;
idamnjanovic@1 28 elseif isobject(Dict)
idamnjanovic@1 29 D =@(z) Dict*z;
idamnjanovic@1 30 Dt =@(z) Dict'*z;
idamnjanovic@1 31 elseif isa(Dict,'function_handle')
idamnjanovic@1 32 try
idamnjanovic@1 33 DictT=varargin{1};
idamnjanovic@1 34 if isa(DictT,'function_handle');
idamnjanovic@1 35 D=Dict;
idamnjanovic@1 36 Dt=DictT;
idamnjanovic@1 37 else
idamnjanovic@1 38 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1 39 end
idamnjanovic@1 40 catch
idamnjanovic@1 41 error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.');
idamnjanovic@1 42 end
idamnjanovic@1 43 else
idamnjanovic@1 44 error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1 45 end
idamnjanovic@1 46 %%
idamnjanovic@1 47 [n,P]=size(X);
idamnjanovic@1 48
idamnjanovic@1 49
idamnjanovic@1 50
idamnjanovic@1 51 global opts opts_tr machPrec
idamnjanovic@1 52 opts.UT = true;
idamnjanovic@1 53 opts_tr.UT = true; opts_tr.TRANSA = true;
idamnjanovic@1 54 machPrec = 1e-5;
idamnjanovic@1 55
idamnjanovic@1 56 A = sparse(m,size(X,2));
idamnjanovic@1 57 for k=1:1:P,
idamnjanovic@1 58
idamnjanovic@1 59 R_I = [];
idamnjanovic@1 60 x=X(:,k);
idamnjanovic@1 61 residual=x;
idamnjanovic@1 62 indx = [];
idamnjanovic@1 63 a = zeros(m,1);
idamnjanovic@1 64 currResNorm = norm(residual);
idamnjanovic@1 65 errorGoal=errorGoal*currResNorm;
idamnjanovic@1 66 j = 0;
idamnjanovic@1 67 while currResNorm>errorGoal & j < maxNumCoef,
idamnjanovic@1 68 j = j+1;
idamnjanovic@1 69 dir=Dt(residual);
idamnjanovic@1 70
idamnjanovic@1 71 [tmp__, pos]=max(abs(dir));
idamnjanovic@1 72
idamnjanovic@1 73 [R_I, flag] = updateChol(R_I, n, m, D, explicitD, indx, pos, Dt);
idamnjanovic@1 74
idamnjanovic@1 75
idamnjanovic@1 76 indx(j)=pos;
idamnjanovic@1 77 dx=zeros(m,1);
idamnjanovic@1 78
idamnjanovic@1 79 z = linsolve(R_I,dir(indx),opts_tr);
idamnjanovic@1 80
idamnjanovic@1 81 dx(indx) = linsolve(R_I,z,opts);
idamnjanovic@1 82 a(indx) = a(indx) + dx(indx);
idamnjanovic@1 83
idamnjanovic@1 84 residual=x-D(a);
idamnjanovic@1 85 currResNorm = norm(residual);
idamnjanovic@1 86
idamnjanovic@1 87
idamnjanovic@1 88 end;
idamnjanovic@1 89 if (~isempty(indx))
idamnjanovic@1 90 A(indx,k)=a(indx);
idamnjanovic@1 91 end
idamnjanovic@1 92 end;
idamnjanovic@1 93 return;
idamnjanovic@1 94
idamnjanovic@1 95
idamnjanovic@1 96 function [R, flag] = updateChol(R, n, N, A, explicitA, activeSet, newIndex, varargin)
idamnjanovic@1 97
idamnjanovic@1 98 % updateChol: Updates the Cholesky factor R of the matrix
idamnjanovic@1 99 % A(:,activeSet)'*A(:,activeSet) by adding A(:,newIndex)
idamnjanovic@1 100 % If the candidate column is in the span of the existing
idamnjanovic@1 101 % active set, R is not updated, and flag is set to 1.
idamnjanovic@1 102
idamnjanovic@1 103 global opts_tr machPrec
idamnjanovic@1 104 flag = 0;
idamnjanovic@1 105
idamnjanovic@1 106 if (explicitA)
idamnjanovic@1 107 newVec = A(:,newIndex);
idamnjanovic@1 108 else
idamnjanovic@1 109 At=varargin{1};
idamnjanovic@1 110 e = zeros(N,1);
idamnjanovic@1 111 e(newIndex) = 1;
idamnjanovic@1 112 newVec = A(e);%feval(A,1,n,N,e,1:N,N);
idamnjanovic@1 113 end
idamnjanovic@1 114
idamnjanovic@1 115 if isempty(activeSet),
idamnjanovic@1 116 R = sqrt(sum(newVec.^2));
idamnjanovic@1 117 else
idamnjanovic@1 118 if (explicitA)
idamnjanovic@1 119 p = linsolve(R,A(:,activeSet)'*A(:,newIndex),opts_tr);
idamnjanovic@1 120 else
idamnjanovic@1 121 AnewVec = At(newVec);%feval(A,2,n,length(activeSet),newVec,activeSet,N);
idamnjanovic@1 122 p = linsolve(R,AnewVec(activeSet),opts_tr);
idamnjanovic@1 123 end
idamnjanovic@1 124 q = sum(newVec.^2) - sum(p.^2);
idamnjanovic@1 125 if (q <= machPrec) % Collinear vector
idamnjanovic@1 126 flag = 1;
idamnjanovic@1 127 else
idamnjanovic@1 128 R = [R p; zeros(1, size(R,2)) sqrt(q)];
idamnjanovic@1 129 end
idamnjanovic@1 130 end