NAG Toolbox: nag_best_subset_given_size (h05ab)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{mincr}$ – int64int32nag_int scalar

Flag indicating whether the scoring function $f$ is increasing or decreasing.

$\mathbf{mincr}=1$

$f\left(S_i\right)\le f\left(S_j\right)$, i.e., the subsets with the largest score will be selected.

$\mathbf{mincr}=0$

$f\left(S_i\right)\ge f\left(S_j\right)$, i.e., the subsets with the smallest score will be selected.

For all $S_j\subseteq \Omega$ and $S_i\subseteq S_j$

Constraint: $\mathbf{mincr}=0$ or $1$.

2:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of features in the full feature set.

Constraint: $\mathbf{m}\ge 2$.

3:     $\mathrm{ip}$ – int64int32nag_int scalar

$p$, the number of features in the subset of interest.

Constraint: $1\le \mathbf{ip}\le \mathbf{m}$.

4:     $\mathrm{nbest}$ – int64int32nag_int scalar

$n$, the maximum number of best subsets required. The actual number of subsets returned is given by la on final exit. If on final exit $\mathbf{la}\ne \mathbf{nbest}$ then $\mathbf{ifail}=\mathbf{42}$ is returned.

Constraint: $\mathbf{nbest}\ge 1$.

5:     $\mathrm{f}$ – function handle or string containing name of m-file

f must evaluate the scoring function $f$.

[score, user, info] = f(m, drop, lz, z, la, a, user, info)

Input Parameters

1:     $\mathrm{m}$ – int64int32nag_int scalar

$m=\left|\Omega \right|$, the number of features in the full feature set.

2:     $\mathrm{drop}$ – int64int32nag_int scalar

Flag indicating whether the intermediate subsets should be constructed by dropping features from the full set ($\mathbf{drop}=1$) or adding features to the empty set ($\mathbf{drop}=0$). See score for additional details.

3:     $\mathrm{lz}$ – int64int32nag_int scalar

The number of features stored in z.

4:     $\mathrm{z}\left(\mathbf{lz}\right)$ – int64int32nag_int array

$\mathbf{z}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{lz}$, contains the list of features which, along with those specified in a, define the subsets whose score is required. See score for additional details.

5:     $\mathrm{la}$ – int64int32nag_int scalar

If $\mathbf{la}>0$, the number of subsets for which a score must be returned.

If $\mathbf{la}=0$, the score for a single subset should be returned. See score for additional details.

6:     $\mathrm{a}\left(\mathbf{la}\right)$ – int64int32nag_int array

$\mathbf{a}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{la}$, contains the list of features which, along with those specified in z, define the subsets whose score is required. See score for additional details.

7:     $\mathrm{user}$ – Any MATLAB object

f is called from nag_best_subset_given_size (h05ab) with the object supplied to nag_best_subset_given_size (h05ab).

8:     $\mathrm{info}$ – int64int32nag_int scalar

$\mathbf{info}=0$.

Output Parameters

1:     $\mathrm{score}\left(\max \left(\mathbf{la},1\right)\right)$ – double array

The value $f\left(S_{\mathit{j}}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{la}$, the score associated with the $j$th subset. $S_j$ is constructed as follows:

$\mathbf{drop}=1$

$S_j$ is constructed by dropping the features specified in the first lz elements of z and the single feature given in $\mathbf{a}\left(j\right)$ from the full set of features, $\Omega$. The subset will therefore contain $\mathbf{m}-\mathbf{lz}-1$ features.

$\mathbf{drop}=0$

$S_j$ is constructed by adding the features specified in the first lz elements of z and the single feature specified in $\mathbf{a}\left(j\right)$ to the empty set, $\varnothing$. The subset will therefore contain $\mathbf{lz}+1$ features.

In both cases the individual features are referenced by the integers $1$ to m with $1$ indicating the first feature, $2$ the second, etc., for some arbitrary ordering of the features, chosen by you prior to calling nag_best_subset_given_size (h05ab). For example, $1$ might refer to the first variable in a particular set of data, $2$ the second, etc..

If $\mathbf{la}=0$, the score for a single subset should be returned. This subset is constructed by adding or removing only those features specified in the first lz elements of z. If $\mathbf{lz}=0$, this subset will either be $\Omega$ or $\varnothing$.

2:     $\mathrm{user}$ – Any MATLAB object

3:     $\mathrm{info}$ – int64int32nag_int scalar

Set info to a nonzero value if you wish nag_best_subset_given_size (h05ab) to terminate with $\mathbf{ifail}=\mathbf{82}$.

6:     $\mathrm{mincnt}$ – int64int32nag_int scalar

$k$, the minimum number of times the effect of each feature, $x_i$, must have been observed before $f\left(S-\left\{x_i\right\}\right)$ is estimated from $f\left(S\right)$ as opposed to being calculated directly.

If $k=0$ then $f\left(S-\left\{x_i\right\}\right)$ is never estimated. If $\mathbf{mincnt}<0$ then $k$ is set to $1$.

7:     $\mathrm{gamma}$ – double scalar

$\gamma$, the scaling factor used when estimating scores. If $\mathbf{gamma}<0$ then $\gamma =1$ is used.

8:     $\mathrm{acc}\left(2\right)$ – double array

A measure of the accuracy of the scoring function, $f$.

Letting $a_i={\epsilon }_1\left|f\left(S_i\right)\right|+{\epsilon }_2$, then when confirming whether the scoring function is strictly increasing or decreasing (as described in mincr), or when assessing whether a node defined by subset $S_i$ can be trimmed, then any values in the range $f\left(S_i\right)\pm a_i$ are treated as being numerically equivalent.

If $0\le \mathbf{acc}\left(1\right)\le 1$ then ${\epsilon }_1=\mathbf{acc}\left(1\right)$, otherwise ${\epsilon }_1=0$.

If $\mathbf{acc}\left(2\right)\ge 0$ then ${\epsilon }_2=\mathbf{acc}\left(2\right)$, otherwise ${\epsilon }_2=0$.

In most situations setting both ${\epsilon }_1$ and ${\epsilon }_2$ to zero should be sufficient. Using a nonzero value, when one is not required, can significantly increase the number of subsets that need to be evaluated.

Optional Input Parameters

1:     $\mathrm{user}$ – Any MATLAB object

user is not used by nag_best_subset_given_size (h05ab), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Output Parameters

1:     $\mathrm{la}$ – int64int32nag_int scalar

The number of best subsets returned.

2:     $\mathrm{bscore}\left(\mathbf{nbest}\right)$ – double array

Holds the score for the la best subsets returned in bz.

3:     $\mathrm{bz}\left(\mathbf{m}-\mathbf{ip},\mathbf{nbest}\right)$ – int64int32nag_int array

The $j$th best subset is constructed by dropping the features specified in $\mathbf{bz}\left(\mathit{i},\mathit{j}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{m}-\mathbf{ip}$ and $\mathit{j}=1,2,\dots ,\mathbf{la}$, from the set of all features, $\Omega$. The score for the $j$th best subset is given in $\mathbf{bscore}\left(j\right)$.

4:     $\mathrm{user}$ – Any MATLAB object

5:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

   $\mathbf{ifail}=11$

Constraint: $\mathbf{mincr}=0$ or $1$.

   $\mathbf{ifail}=21$

Constraint: $\mathbf{m}\ge 2$.

   $\mathbf{ifail}=31$

Constraint: $1\le \mathbf{ip}\le \mathbf{m}$.

   $\mathbf{ifail}=41$

Constraint: $\mathbf{nbest}\ge 1$.

W  $\mathbf{ifail}=42$

On entry, $\mathbf{nbest}=\_$.
But only $\_$ best subsets could be calculated.

   $\mathbf{ifail}=81$

On exit from f, $\mathbf{score}\left(\_\right)=\_$, which is inconsistent with the score for the parent node.

   $\mathbf{ifail}=82$

A nonzero value for info has been returned:

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: h05ab_example

function h05ab_example


fprintf('h05ab example results\n\n');

m      = int64(14);
ip     = int64(5);
nbest  = int64(3);
mincr  = int64(0);
mincnt = int64(-1);
gamma  = -1;
acc    = [-1, -1];
mip    = m - ip;
z      = zeros(mip, 1, 'int64');
a      = zeros(max(nbest, m), 1, 'int64');
bz     = zeros(mip, nbest, 'int64');

% Data required by the scoring function
n = int64(40);
[x,y] = gen_data(n,m);

% Put x and y in user to be passed to f, along with a placeholder for the
% number of subsets evaluated
user = {x; y; 0};

% Call the forward communication best subset routine
[la, bscore, bz, user, ifail] = ...
    h05ab(...
          mincr, m, ip, nbest, @f, mincnt, gamma, acc, 'user', user);

fprintf('\n    Score        Feature Subset\n');
fprintf('    -----        --------------\n');
% Display the best subsets and corresponding scores. h05aa returns
% a list of features excluded from the best subsets, so this is inverted
% to give the set of features included in each subset
ibz = 1:m;
for i = 1:la
  mask = ones(1, m, 'int64');
  for j = 1:mip
    mask(bz(j, i)) = 0;
  end
  fprintf('%12.5e %5d %5d %5d %5d %5d\n', bscore(i), ibz(logical(mask)));
end
fprintf('\n%d subsets evaluated in total\n', user{3});



function [score, user, info] = f(m, drop, lz, z, la, a, user, info)
  % Set up the initial feature set.
  % If drop = 0, this is the Null set (i.e. no features).
  % If drop = 1 then this is the full set (i.e. all features)
  if drop == 0
    isx = zeros(m, 1, 'int64');
  else
    isx = ones(m, 1, 'int64');
  end

  % Add (if drop = 0) or remove (if drop = 1) all the features specified in z
  inv_drop = not(drop);
  for i = 1:lz
    isx(z(i)) = inv_drop;
  end

  score = zeros(max(la, 1), 1);
  for i=1:max(la, 1)
    if (la > 0)
      if (i > 1)
        % Reset the feature altered at the last iteration
        isx(a(i-1)) = drop;
      end

      %  Add or drop the i'th feature in a
      isx(a(i)) = inv_drop;
    end

    ip = int64(sum(isx));

    % Fit the regression model
    rss = g02da('z', user{1}, isx, ip, user{2});

    % Return the score (the residual sums of squares)
    score(i) = rss;
  end

  % Keep track of the number of subsets evaluated
  user{3} = user{3} + max(1, la);

function [x, y] = gen_data(n,m)
  x = zeros(n,m);
  genid = int64(3);
  subid = int64(1);
  seed(1) = int64(23124124);
  [state, ifail] = g05kf(genid, subid, seed);
  for i = 1:m
    [state, x(1:n,i), ifail] = g05sk( ...
                                      n, 0, sqrt(3), state);
  end
  [state, b, ifail] = g05sk(m, 1.5, 3, state);
  [state, y, ifail] = g05sk(n, 0, 1, state);
  y = x*b + y;

h05ab example results


    Score        Feature Subset
    -----        --------------
 1.21567e+03     3     4     6     7    14
 1.32495e+03     3     5     6     7    14
 1.37244e+03     3     5     6    12    14

337 subsets evaluated in total