Download All

Code covered by the BSD License  

Highlights from
Statistical Learning Toolbox

from Statistical Learning Toolbox by Dahua Lin
Functions for statistical learning, pattern recognition and computer vision, covering many topics.

Description of slnbreconweights
Home > sltoolbox > manifold > slnbreconweights.m

slnbreconweights

PURPOSE ^

SLNBRECONWEIGHTS Solve the optimal reconstruction weights on given neighbors

SYNOPSIS ^

function WG = slnbreconweights(X0, X, G, varargin)

DESCRIPTION ^

SLNBRECONWEIGHTS Solve the optimal reconstruction weights on given neighbors

 $ Syntax $
   - WG = slnbreconweights(X0, X, G, ...)

 $ Arguments $
   - X0:       The reference samples to reconstruct the query samples
   - X:        The query samples
   - G:        The graph giving the neighborhood relations 
   - WG:       The weighted graph giving the solved weights

 $ Description $
   - WG = slnbreconweights(X0, X, G, ...) solves the optimal weights to
     reconstruct the samples in X from those in X0. If X is empty, then
     it would use X0 as X. The graph G indicates the neighborhood 
     relation, having n0 sources and n targets. The WG is a graph with
     of the same size as G, and the reconstuction weights are placed
     in the positions corresponding to those in G. 
     You can specify the following properties to control the solving:
     \*
     \t    Table. The Properties of Reconstruction Weights Solving
     \h       name        &     description
            'constraint'  & The constraint on the solution, it can be
                            one of the following string to indicate a
                            single constraint or a cell array of multiple
                            strings to indicate compound constaints.
                            - 'nonneg':  non-negative
                            - 's1':      the weights sum to 1
                            (default = 's1')
            'delta'       & The value of regularization. In practice, 
                            regularization is essential to guarantee the
                            stability of the solution. In implementation,
                            the diagonal elements of the gram matrix will
                            be added with a value:
                               (delta^2) * trace(G) / K
                            here G is X^T * X, K is the neighbor number.
                            (default = 0.1)
            'solver'      & The solver offered by user (function handle).  
                            If the user specify a non-empty solver, then 
                            it will use the user's solver to solve 
                            weights. The solver is like the form:
                               w = f(X, y)
                            Here X is d x K neighbor sample matrix, y is
                            a d x 1 vector representing the target sample.
                            It should output a K x 1 vector giving the
                            reconstruction weights. 
                            By default, solver = [], indicating to use
                            internal solver based on constraint and delta.   
            'thres'       & The thres, if the ratio of a weight value
                            to the average weight for that reconstruction
                            is lower than thres, the weight is set 
                            to strictly zeros. This would significantly
                            reduces the near-zero weights, and thus
                            reduces the complexity of the graph.
                            (default = 1e-8)
     \*

 $ Remarks $
   - When the user specify a non-empty solver, the internal solver will
     not be used, thus constraint and delta will not take effect.

   - G would be in all acceptable graph form. WG will always be a 
     numeric matrix. If G is a sparse adjmat, then WG would be sparse,
     otherwise WG is full.

   - With the WG solved, to reconstruct by weighed combination of
     neighbors, you can simply write it as: Xr = X0 * WG, then Xr
     is a d x n matrix with the j-th column reconstructed from the 
     referenced samples in Xr using the j-th column's weights in WG.

 $ History $
   - Created by Dahua Lin, on Sep 11st, 2006

CROSS-REFERENCE INFORMATION ^

This function calls:
  • sladjmat SLADJMAT Constructs the adjacency matrix representation of a graph
  • slgraphinfo SLGRAPHINFO Extracts basic information of a given graph representation
  • raise_lackinput RAISE_LACKINPUT Raises an error indicating lack of input argument
  • slparseprops SLPARSEPROPS Parses input parameters
This function is called by:
  • sllle SLLLE Performs Locally Linear Embedding

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function WG = slnbreconweights(X0, X, G, varargin)
0002 %SLNBRECONWEIGHTS Solve the optimal reconstruction weights on given neighbors
0003 %
0004 % $ Syntax $
0005 %   - WG = slnbreconweights(X0, X, G, ...)
0006 %
0007 % $ Arguments $
0008 %   - X0:       The reference samples to reconstruct the query samples
0009 %   - X:        The query samples
0010 %   - G:        The graph giving the neighborhood relations
0011 %   - WG:       The weighted graph giving the solved weights
0012 %
0013 % $ Description $
0014 %   - WG = slnbreconweights(X0, X, G, ...) solves the optimal weights to
0015 %     reconstruct the samples in X from those in X0. If X is empty, then
0016 %     it would use X0 as X. The graph G indicates the neighborhood
0017 %     relation, having n0 sources and n targets. The WG is a graph with
0018 %     of the same size as G, and the reconstuction weights are placed
0019 %     in the positions corresponding to those in G.
0020 %     You can specify the following properties to control the solving:
0021 %     \*
0022 %     \t    Table. The Properties of Reconstruction Weights Solving
0023 %     \h       name        &     description
0024 %            'constraint'  & The constraint on the solution, it can be
0025 %                            one of the following string to indicate a
0026 %                            single constraint or a cell array of multiple
0027 %                            strings to indicate compound constaints.
0028 %                            - 'nonneg':  non-negative
0029 %                            - 's1':      the weights sum to 1
0030 %                            (default = 's1')
0031 %            'delta'       & The value of regularization. In practice,
0032 %                            regularization is essential to guarantee the
0033 %                            stability of the solution. In implementation,
0034 %                            the diagonal elements of the gram matrix will
0035 %                            be added with a value:
0036 %                               (delta^2) * trace(G) / K
0037 %                            here G is X^T * X, K is the neighbor number.
0038 %                            (default = 0.1)
0039 %            'solver'      & The solver offered by user (function handle).
0040 %                            If the user specify a non-empty solver, then
0041 %                            it will use the user's solver to solve
0042 %                            weights. The solver is like the form:
0043 %                               w = f(X, y)
0044 %                            Here X is d x K neighbor sample matrix, y is
0045 %                            a d x 1 vector representing the target sample.
0046 %                            It should output a K x 1 vector giving the
0047 %                            reconstruction weights.
0048 %                            By default, solver = [], indicating to use
0049 %                            internal solver based on constraint and delta.
0050 %            'thres'       & The thres, if the ratio of a weight value
0051 %                            to the average weight for that reconstruction
0052 %                            is lower than thres, the weight is set
0053 %                            to strictly zeros. This would significantly
0054 %                            reduces the near-zero weights, and thus
0055 %                            reduces the complexity of the graph.
0056 %                            (default = 1e-8)
0057 %     \*
0058 %
0059 % $ Remarks $
0060 %   - When the user specify a non-empty solver, the internal solver will
0061 %     not be used, thus constraint and delta will not take effect.
0062 %
0063 %   - G would be in all acceptable graph form. WG will always be a
0064 %     numeric matrix. If G is a sparse adjmat, then WG would be sparse,
0065 %     otherwise WG is full.
0066 %
0067 %   - With the WG solved, to reconstruct by weighed combination of
0068 %     neighbors, you can simply write it as: Xr = X0 * WG, then Xr
0069 %     is a d x n matrix with the j-th column reconstructed from the
0070 %     referenced samples in Xr using the j-th column's weights in WG.
0071 %
0072 % $ History $
0073 %   - Created by Dahua Lin, on Sep 11st, 2006
0074 %
0075 
0076 %% parse and verify input arguments
0077 
0078 if nargin < 3
0079     raise_lackinput('slnbreconweights', 3);
0080 end
0081 
0082 if ~isnumeric(X0) || ndims(X0) ~= 2
0083     error('sltoolbox:invalidarg', ...
0084         'X should be a 2D numeric matrix');
0085 end
0086 [d, n0] = size(X0);
0087 
0088 
0089 if isempty(X)
0090     X = X0;    
0091 else
0092     if ~isnumeric(X) || ndims(X) ~= 2
0093         error('sltoolbox:invalidarg', ...
0094             'X0 should be a 2D numeric matrix');
0095     end        
0096     if size(X, 1) ~= d
0097         error('sltoolbox:sizmismatch', ...
0098             'The sample dimension in X is not the same as that in X0');
0099     end
0100 end
0101 n = size(X, 2);
0102 
0103 gi = slgraphinfo(G);
0104 if gi.n ~= n0 || gi.nt ~= n
0105     error('sltoolbox:sizmismatch', ...
0106         'The size of the graph is not consisitent with the sample set');
0107 end
0108 
0109 opts.constraint = 's1';
0110 opts.delta = 0.1;
0111 opts.solver = [];
0112 opts.thres = 1e-8;
0113 opts = slparseprops(opts, varargin{:});
0114 
0115 thres = opts.thres;
0116 
0117 %% Prepare parameters
0118 
0119 % prepare graph
0120 if ~strcmp(gi.form, 'adjmat')
0121     G = sladjmat(G, 'sparse', true, 'valtype', 'logical');
0122 end
0123     
0124 % prepare solver
0125 if isempty(opts.solver)
0126        
0127     % parse constraints
0128     cs = opts.constraint;
0129     if ~isempty(cs)
0130         if ~iscell(cs)
0131             cs = {cs};
0132         end
0133         constraint = parse_constraints(cs);
0134     else
0135         constraint = parse_constraints({});
0136     end
0137     
0138     % decide solver
0139     if ~constraint.nonneg
0140         delta2 = opts.delta^2;
0141         if ~constraint.s1
0142             wsolver = @(X, y) internal_wsolver_unc(X, y, delta2);
0143         else
0144             wsolver = @(X, y) internal_wsolver_s1(X, y, delta2);
0145         end
0146     else
0147         optimopts = optimset('Display', 'off', 'LargeScale', 'off');
0148         if ~constraint.s1
0149             wsolver = @(X, y) internal_wsolver_nonneg(X, y, opts.delta, optimopts);
0150         else
0151             wsolver = @(X, y) internal_wsolver_nonneg_s1(X, y, opts.delta, optimopts);
0152         end
0153     end                           
0154     
0155 else
0156     if ~isa(opts.solver, 'function_handle')
0157         error('The weight solver should be a function handle');
0158     end
0159     wsolver = opts.solver;
0160 end
0161 
0162 
0163 %% main skeleton
0164 
0165 % init WG
0166 if issparse(G)
0167     WG = spalloc(n0, n, nnz(G));
0168 else
0169     WG = zeros(n0, n);
0170 end
0171 
0172 % solve weights
0173 for i = 1 : n
0174     nbinds = find(G(:,i));
0175     if ~isempty(nbinds)        
0176         Xnb = X0(:, nbinds);
0177         y = X(:,i);
0178         w = wsolver(Xnb, y);
0179         if thres > 0
0180             absw = abs(w);
0181             curthres = thres * sum(absw) / length(w);
0182             w(absw < curthres) = 0;
0183         end        
0184         WG(nbinds, i) = w;
0185     end    
0186 end
0187 
0188 
0189 %% constraint parsing function
0190 
0191 function c = parse_constraints(cs)
0192 
0193 ncs = length(cs);
0194 
0195 c = struct('nonneg', false, 's1', false);
0196 
0197 for i = 1 : ncs
0198     cname = cs{i};
0199     if ~ischar(cname)
0200         error('sltoolbox:invalidarg', ...
0201             'The constraint should be given in char string');
0202     end
0203     switch cname
0204         case 'nonneg'
0205             c.nonneg = true;
0206         case 's1'
0207             c.s1 = true;
0208         otherwise
0209             error('sltoolbox:invalidarg', ...
0210                 'Invalid constraint name for weight solving: %s', cname);
0211     end
0212 end
0213 
0214 
0215 
0216 %% The internal weight solvers
0217 
0218 % unconstrained solver
0219 function w = internal_wsolver_unc(X, y, delta2)
0220 
0221 [G, Xty] = compute_G_Xty(X, y, delta2);
0222 w = G \ Xty; 
0223 
0224 % solver with s1 constraint
0225 function w = internal_wsolver_s1(X, y, delta2)
0226 
0227 [G, Xty, K] = compute_G_Xty(X, y, delta2);
0228 wu = G \ [Xty, ones(K, 1)];
0229 
0230 w = wu(:,1);
0231 u = wu(:,2);
0232 lambda = (1 - sum(w)) / sum(u);
0233 w = w + lambda * u;
0234 
0235 % solver with nonnegative constraint
0236 function w = internal_wsolver_nonneg(X, y, delta, optimopts)
0237 
0238 [Xa, ya, K] = augformulate(X, y, delta);
0239 if K <= 20
0240     w = lsqnonneg(Xa, ya, [], optimopts);
0241 else
0242     lb = zeros(K, 1);
0243     w = lsqlin(Xa, ya, [], [], [], [], lb, [], [], optimopts);
0244 end
0245 
0246 % solver with nonnegative and s1 constraint
0247 function w = internal_wsolver_nonneg_s1(X, y, delta, optimopts)
0248 
0249 [Xa, ya, K] = augformulate(X, y, delta);
0250 
0251 Aeq = ones(1, K);
0252 beq = 1;
0253 lb = zeros(K, 1);
0254 w = lsqlin(Xa, ya, [], [], Aeq, beq, lb, [], [], optimopts);
0255 
0256 
0257 % solver preparation function
0258 
0259 function [G, Xty, K] = compute_G_Xty(X, y, delta2)
0260 
0261 % compute Xt, G, and Xty
0262 K = size(X, 2);
0263 Xt = X';
0264 G = Xt * X;
0265 Xty = Xt * y;
0266 
0267 % regularize
0268 if delta2 > 0
0269     diaginds = (1:K)*(K+1) - K;
0270     rv = delta2 * sum(G(diaginds)) / K;
0271     G(diaginds) = G(diaginds) + rv;
0272 end
0273 
0274 function [Xa, ya, K] = augformulate(X, y, delta)
0275 
0276 K = size(X, 2);
0277 if delta ~= 0
0278     Xa = [X; delta * eye(K)];
0279     ya = [y; zeros(K, 1)];
0280 else
0281     Xa = X;
0282     ya = y;
0283 end
0284
Generated on Wed 20-Sep-2006 12:43:11 by m2html © 2003

Contact us at files@mathworks.com