www.pudn.com > SURE_SVT_code.zip > sure_spec_fun.m, change:2012-10-15,size:5122b


function    R   =   sure_spec_fun(f, df, sigma, varargin)
    %
    %   SURE_SVT:   Stein's Unbiased Risk Estimate (SURE) for matrix-valued
    %               spectral functions
    %
    %   USAGE:
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, Y)
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, S, [M N])    
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, S, [M N], IS_REAL)
    %
    %   DESCRIPTION:
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, Y)
    %       takes as inputs: 
    %       (i)     The function handles F and DF, that use as input the 
    %               singular value (column) vector S and returns the vector 
    %               F(S) of values of F at S, and the vector DF(S) of the 
    %               values of the derivative of F at S,
    %       (ii)    The standard deviation of the noise SIGMA > 0,
    %       (iii)   The observed matrix Y, which can be real- or 
    %               complex-valued,
    %       and it returns R, the value of SURE evaluated at Y.
    %
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, S, [M N])
    %       R   =   SURE_SPEC_FUN(F, DF, SIGMA, S, [M N], IS_REAL)
    %       accepts the vector S of observed singular values, and the 
    %       vector [M N] containing the number of rows and columns of the 
    %       observed matrix. The optional flag IS_REAL indicates whether 
    %       the observed matrix is real-valued (if IS_REAL = 1) or not. 
    %       The default is IS_REAL = 1.
    %
    %   REFERENCES:
    %
    %       Candes, E.J., Sing-Long, C.A., and Trzasko, J.D.,
    %       ``Unbiased Risk Estimates for Singular Value Thresholding''
    %
    %   V1.0:   Oct. 2012.
    %
    if( nargin == 4 ),
        [M N]       =   size(varargin{1});
        is_real     =   isreal(varargin{1});
        s           =   svd(varargin{1});
    elseif( nargin == 5 ),
        s           =   varargin{1};
        M           =   varargin{2}(1);
        N           =   varargin{2}(2);
        is_real     =   1;
    elseif( nargin == 6 ),
        s           =   varargin{1};
        M           =   varargin{2}(1);
        N           =   varargin{2}(2);
        is_real     =   varargin{3};
    end
    
    R           =   sure_div_spec_fun(f, df, s, is_real, [M N]);
    
    %   additional terms
    R           =   -M*N*sigma^2 + sum((f(s) - s).^2) + 2*sigma^2*R;    
    if( ~is_real ),
        R       =   R - M*N*sigma^2;
    end
end
function    x   =   sure_div_spec_fun(f, df, s, is_real, sz)
    svThreshold 	=   1E-8;       %   threshold to determine whether two singular
                                    %   values are the same
    M               =   sz(1);
    N               =   sz(2);
    
    %   *** safeguard
    %       check multiplicities of singular values in a robust manner
    z           =   s(2:end);
    s           =   [s(1) 1];
    Is          =   1;
    while( ~isempty(z) ),
        idx         =   find(abs(z - s(Is, 1)) < svThreshold );
        if( isempty(idx) )
            s           =   [s; [z(1) 1]];
            z(1)        =   [];
            Is          =   Is + 1;
        end
        z(idx)      =   [];
        s(Is, 2)    =   s(Is, 2) + numel(idx);
    end
    clear z
    
    %       warns the user about using SURE with a non-simple, not-full
    %       rank matrix
    if( any( s(:, 1) < svThreshold ) ),
        fprintf('   +   [SURE_SPEC_FUN] Warning: argument might be rank-deficient.\n')
    end
    if( any( s(:, 2) > 1 ) ),
        fprintf('   +   [SURE_SPEC_FUN] Warning: argument might have repeated singular values.\n')
    end

    %       find positive singular values
    idx_p   =   ( s(:, 1) > svThreshold );
    
    if( is_real ),
        x   =   div_spec_fun_real(f, df, s, idx_p, M, N);
    else
        x   =   div_spec_fun_complex(f, df, s, idx_p, M, N);
    end
end
function    x   =   div_spec_fun_real(f, df, s, idx_p, M, N)
    x  	=   0;
    if( any(idx_p)  ),
        x   =   x + sum( 0.5*s(idx_p, 2).*(s(idx_p, 2) + 1).*df(s(idx_p, 1)) );
        x   =   x + sum( s(idx_p, 2).*(abs(M-N) + 0.5*(s(idx_p, 2) - 1)).*(f(s(idx_p, 1))./s(idx_p, 1)) );
    end
    if( any(~idx_p) ),
        x   =   x + df(0)*sum( (abs(M - N) + s(~idx_p, 2)).*s(~idx_p, 2) );
    end
    
    D  	=   zeros(size(s, 1));
    for Ik = 1:size(s, 1),
        D(:, Ik)    =   s(Ik, 2)*s(:, 2).*s(:, 1).*f( s(:, 1) )./(s(:, 1).^2 - s(Ik, 1).^2);
    end
    D( isnan(D) | isinf(D) | abs(D) > 1E6 )     =   0;
    
    x 	=   x + 2*sum(D(:));
end
function    x   =   div_spec_fun_complex(f, df, s, idx_p, M, N)
    x  	=   0;
    if( any(idx_p)  ),
        x   =   x + sum( (s(idx_p, 2).^2).*df(s(idx_p, 1)) );
        x   =   x + sum( s(idx_p, 2).*(2*abs(M-N) + s(idx_p, 2)).*(f(s(idx_p, 1))./s(idx_p, 1)) );
    end
    if( any(~idx_p) ),
        x   =   x + 2*df(0)*sum( (abs(M - N) + s(~idx_p, 2)).*s(~idx_p, 2) );
    end
    
    D  	=   zeros(size(s, 1));
    for Ik = 1:size(s, 1),
        D(:, Ik)    =   s(Ik, 2)*s(:, 2).*s(:, 1).*f( s(:, 1) )./(s(:, 1).^2 - s(Ik, 1).^2);
    end
    D( isnan(D) | isinf(D) | abs(D) > 1E6 )     =   0;
    
    x 	=   x + 4*sum(D(:));
end