NAG Library Routine Document
h05abf
(best_subset_given_size)
1
Purpose
Given a set of $m$ features and a scoring mechanism for any subset of those features, h05abf selects the best $n$ subsets of size $p$ using a direct communication branch and bound algorithm.
2
Specification
Fortran Interface
Subroutine h05abf ( 
mincr,
m,
ip,
nbest,
la,
bscore,
bz,
f,
mincnt,
gamma,
acc,
iuser,
ruser,
ifail) 
Integer, Intent (In)  :: 
mincr,
m,
ip,
nbest,
mincnt  Integer, Intent (Inout)  :: 
iuser(*),
ifail  Integer, Intent (Out)  :: 
la,
bz(mip,nbest)  Real (Kind=nag_wp), Intent (In)  :: 
gamma,
acc(2)  Real (Kind=nag_wp), Intent (Inout)  :: 
ruser(*)  Real (Kind=nag_wp), Intent (Out)  :: 
bscore(nbest)  External  :: 
f 

C Header Interface
#include nagmk26.h
void 
h05abf_ (
const Integer *mincr,
const Integer *m,
const Integer *ip,
const Integer *nbest,
Integer *la,
double bscore[],
Integer bz[],
void (NAG_CALL *f)(
const Integer *m,
const Integer *drop,
const Integer *lz,
const Integer z[],
const Integer *la,
const Integer a[],
double score[],
Integer iuser[],
double ruser[],
Integer *info),
const Integer *mincnt,
const double *gamma,
const double acc[],
Integer iuser[],
double ruser[],
Integer *ifail) 

3
Description
Given $\Omega =\left\{{x}_{\mathrm{i}}:i\in \mathbb{Z},1\le i\le m\right\}$, a set of $m$ unique features and a scoring mechanism $f\left(S\right)$ defined for all $S\subseteq \Omega $ then h05abf is designed to find ${S}_{o1}\subseteq \Omega ,\left{S}_{o1}\right=p$, an optimal subset of size $p$. Here $\left{S}_{o1}\right$ denotes the cardinality of ${S}_{o1}$, the number of elements in the set.
The definition of the optimal subset depends on the properties of the scoring mechanism, if
then the optimal subset is defined as one of the solutions to
else if
then the optimal subset is defined as one of the solutions to
If neither of these properties hold then h05abf cannot be used.
As well as returning the optimal subset,
${S}_{o1}$,
h05abf can return the best
$n$ solutions of size
$p$. If
${S}_{o\mathit{i}}$ denotes the
$\mathit{i}$th best subset, for
$\mathit{i}=1,2,\dots ,n1$, then the
$\left(i+1\right)$th best subset is defined as the solution to either
or
depending on the properties of
$f$.
The solutions are found using a branch and bound method, where each node of the tree is a subset of
$\Omega $. Assuming that
(1) holds then a particular node, defined by subset
${S}_{i}$, can be trimmed from the tree if
$f\left({S}_{i}\right)<\hat{f}\left({S}_{on}\right)$ where
$\hat{f}\left({S}_{on}\right)$ is the
$n$th highest score we have observed so far for a subset of size
$p$, i.e., our current best guess of the score for the
$n$th best subset. In addition, because of
(1) we can also drop all nodes defined by any subset
${S}_{j}$ where
${S}_{j}\subseteq {S}_{i}$, thus avoiding the need to enumerate the whole tree. Similar short cuts can be taken if
(2) holds. A full description of this branch and bound algorithm can be found in
Ridout (1988).
Rather than calculate the score at a given node of the tree
h05abf utilizes the fast branch and bound algorithm of
Somol et al. (2004), and attempts to estimate the score where possible. For each feature,
${x}_{i}$, two values are stored, a count
${c}_{i}$ and
${\hat{\mu}}_{i}$, an estimate of the contribution of that feature. An initial value of zero is used for both
${c}_{i}$ and
${\hat{\mu}}_{i}$. At any stage of the algorithm where both
$f\left(S\right)$ and
$f\left(S\left\{{x}_{i}\right\}\right)$ have been calculated (as opposed to estimated), the estimated contribution of the feature
${x}_{i}$ is updated to
and
${c}_{i}$ is incremented by
$1$, therefore at each stage
${\hat{\mu}}_{i}$ is the mean contribution of
${x}_{i}$ observed so far and
${c}_{i}$ is the number of observations used to calculate that mean.
As long as ${c}_{i}\ge k$, for the usersupplied constant $k$, then rather than calculating $f\left(S\left\{{x}_{\mathrm{i}}\right\}\right)$ this routine estimates it using $\hat{f}\left(S\left\{{x}_{\mathrm{i}}\right\}\right)=f\left(S\right)\gamma {\hat{\mu}}_{i}$ or $\hat{f}\left(S\right)\gamma {\hat{\mu}}_{i}$ if $f\left(S\right)$ has been estimated, where $\gamma $ is a usersupplied scaling factor. An estimated score is never used to trim a node or returned as the optimal score.
Setting
$k=0$ in this routine will cause the algorithm to always calculate the scores, returning to the branch and bound algorithm of
Ridout (1988). In most cases it is preferable to use the fast branch and bound algorithm, by setting
$k>0$, unless the score function is iterative in nature, i.e.,
$f\left(S\right)$ must have been calculated before
$f\left(S\left\{{x}_{i}\right\}\right)$ can be calculated.
h05abf is a direct communication version of
h05aaf.
4
References
Narendra P M and Fukunaga K (1977) A branch and bound algorithm for feature subset selection IEEE Transactions on Computers 9 917–922
Ridout M S (1988) Algorithm AS 233: An improved branch and bound algorithm for feature subset selection Journal of the Royal Statistics Society, Series C (Applied Statistics) (Volume 37) 1 139–147
Somol P, Pudil P and Kittler J (2004) Fast branch and bound algorithms for optimal feature selection IEEE Transactions on Pattern Analysis and Machine Intelligence (Volume 26) 7 900–912
5
Arguments
 1: $\mathbf{mincr}$ – IntegerInput

On entry: 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: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of features in the full feature set.
Constraint:
${\mathbf{m}}\ge 2$.
 3: $\mathbf{ip}$ – IntegerInput

On entry: $p$, the number of features in the subset of interest.
Constraint:
$1\le {\mathbf{ip}}\le {\mathbf{m}}$.
 4: $\mathbf{nbest}$ – IntegerInput

On entry:
$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: $\mathbf{la}$ – IntegerOutput

On exit: the number of best subsets returned.
 6: $\mathbf{bscore}\left({\mathbf{nbest}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: holds the score for the
la best subsets returned in
bz.
 7: $\mathbf{bz}\left({\mathbf{m}}{\mathbf{ip}},{\mathbf{nbest}}\right)$ – Integer arrayOutput

On exit: 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)$.
 8: $\mathbf{f}$ – Subroutine, supplied by the user.External Procedure

f must evaluate the scoring function
$f$.
The specification of
f is:
Fortran Interface
Integer, Intent (In)  :: 
m,
drop,
lz,
z(lz),
la,
a(la)  Integer, Intent (Inout)  :: 
iuser(*),
info  Real (Kind=nag_wp), Intent (Inout)  :: 
ruser(*)  Real (Kind=nag_wp), Intent (Out)  :: 
score(max(la,1)) 

C Header Interface
#include nagmk26.h
void 
f (
const Integer *m,
const Integer *drop,
const Integer *lz,
const Integer z[],
const Integer *la,
const Integer a[],
double score[],
Integer iuser[],
double ruser[],
Integer *info) 

 1: $\mathbf{m}$ – IntegerInput

On entry: $m=\left\Omega \right$, the number of features in the full feature set.
 2: $\mathbf{drop}$ – IntegerInput

On entry: 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: $\mathbf{lz}$ – IntegerInput

On entry: the number of features stored in
z.
 4: $\mathbf{z}\left({\mathbf{lz}}\right)$ – Integer arrayInput

On entry:
${\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: $\mathbf{la}$ – IntegerInput

On entry: 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: $\mathbf{a}\left({\mathbf{la}}\right)$ – Integer arrayInput

On entry:
${\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: $\mathbf{score}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{la}},1\right)\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: 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
h05abf. 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 $.
 8: $\mathbf{iuser}\left(*\right)$ – Integer arrayUser Workspace
 9: $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) arrayUser Workspace

f is called with the arguments
iuser and
ruser as supplied to
h05abf. You should use the arrays
iuser and
ruser to supply information to
f.
 10: $\mathbf{info}$ – IntegerInput/Output

On entry: ${\mathbf{info}}=0$.
On exit: set
info to a nonzero value if you wish
h05abf to terminate with
${\mathbf{ifail}}={\mathbf{82}}$.
f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
h05abf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: f should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
h05abf. If your code inadvertently
does return any NaNs or infinities,
h05abf is likely to produce unexpected results.
 9: $\mathbf{mincnt}$ – IntegerInput

On entry:
$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$.
 10: $\mathbf{gamma}$ – Real (Kind=nag_wp)Input

On entry: $\gamma $, the scaling factor used when estimating scores. If ${\mathbf{gamma}}<0$ then $\gamma =1$ is used.
 11: $\mathbf{acc}\left(2\right)$ – Real (Kind=nag_wp) arrayInput

On entry: a measure of the accuracy of the scoring function,
$f$.
Letting
${a}_{i}={\epsilon}_{1}\leftf\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.
 12: $\mathbf{iuser}\left(*\right)$ – Integer arrayUser Workspace
 13: $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) arrayUser Workspace

iuser and
ruser are not used by
h05abf, but are passed directly to
f and may be used to pass information to this routine.
 14: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=11$

On entry, ${\mathbf{mincr}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{mincr}}=0$ or $1$.
 ${\mathbf{ifail}}=21$

On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}\ge 2$.
 ${\mathbf{ifail}}=31$

On entry, ${\mathbf{ip}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $1\le {\mathbf{ip}}\le {\mathbf{m}}$.
 ${\mathbf{ifail}}=41$

On entry, ${\mathbf{nbest}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{nbest}}\ge 1$.
 ${\mathbf{ifail}}=42$

On entry, ${\mathbf{nbest}}=\u2329\mathit{\text{value}}\u232a$.
But only $\u2329\mathit{\text{value}}\u232a$ best subsets could be calculated.
 ${\mathbf{ifail}}=81$

On exit from
f,
${\mathbf{score}}\left(\u2329\mathit{\text{value}}\u232a\right)=\u2329\mathit{\text{value}}\u232a$, which is inconsistent with the score for the parent node. Score for the parent node is
$\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=82$

A nonzero value for
info has been returned:
${\mathbf{info}}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The subsets returned by h05abf are guaranteed to be optimal up to the accuracy of the calculated scores.
8
Parallelism and Performance
h05abf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
The maximum number of unique subsets of size
$p$ from a set of
$m$ features is
$N=\frac{m!}{\left(mp\right)!p!}$. The efficiency of the branch and bound algorithm implemented in
h05abf comes from evaluating subsets at internal nodes of the tree, that is subsets with more than
$p$ features, and where possible trimming branches of the tree based on the scores at these internal nodes as described in
Narendra and Fukunaga (1977). Because of this it is possible, in some circumstances, for more than
$N$ subsets to be evaluated. This will tend to happen when most of the features have a similar effect on the subset score.
If multiple optimal subsets exist with the same score, and
nbest is too small to return them all, then the choice of which of these optimal subsets is returned is arbitrary.
10
Example
This example finds the three linear regression models, with five variables, that have the smallest residual sums of squares when fitted to a supplied dataset. The data used in this example was simulated.
10.1
Program Text
Program Text (h05abfe.f90)
10.2
Program Data
Program Data (h05abfe.d)
10.3
Program Results
Program Results (h05abfe.r)