Mercurial > hg > smallbox
view solvers/SMALL_chol.m @ 59:23f9dd7b9d78
(none)
author | idamnjanovic |
---|---|
date | Mon, 14 Mar 2011 17:25:38 +0000 (2011-03-14) |
parents | 524cc3fff5ac |
children | e08af264ef93 |
line wrap: on
line source
function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin) %% % % 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. % % 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. % %% % 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