Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_best_subset_given_size (h05ab)

Purpose

Given a set of m$m$ features and a scoring mechanism for any subset of those features, nag_best_subset_given_size (h05ab) selects the best n$n$ subsets of size p$p$ using a direct communication branch and bound algorithm.

Syntax

[la, bscore, bz, user, ifail] = h05ab(mincr, m, ip, nbest, f, mincnt, gamma, acc, 'user', user)
[la, bscore, bz, user, ifail] = nag_best_subset_given_size(mincr, m, ip, nbest, f, mincnt, gamma, acc, 'user', user)

Description

Given Ω = {xi : i,1im}$\Omega =\left\{{x}_{i}:i\in ℤ,1\le i\le m\right\}$, a set of m$m$ unique features and a scoring mechanism f(S) $f\left(S\right)$ defined for all S Ω $S\subseteq \Omega$ then nag_best_subset_given_size (h05ab) is designed to find So1 Ω , |So1| = p ${S}_{o1}\subseteq \Omega ,|{S}_{o1}|=p$, an optimal subset of size p$p$. Here |So1|$|{S}_{o1}|$ denotes the cardinality of So1${S}_{o1}$, the number of elements in the set.
The definition of the optimal subset depends on the properties of the scoring mechanism, if
 f(Si) ≤ f(Sj) , for all ​ Sj ⊆ Ω ​ and ​ Si ⊆ Sj
$f(Si) ≤ f(Sj) , for all ​ Sj ⊆ Ω ​ and ​ Si ⊆ Sj$
(1)
then the optimal subset is defined as one of the solutions to
 maximize f(S)  subject to  |S| = p S ⊆ Ω
$maximize S⊆Ω f(S) subject to |S| = p$
else if
 f(Si) ≥ f(Sj) , for all ​ Sj ⊆ Ω ​ and ​ Si ⊆ Sj
$f( Si ) ≥ f(Sj) , for all ​ Sj ⊆ Ω ​ and ​ Si ⊆ Sj$
(2)
then the optimal subset is defined as one of the solutions to
 minimize f(S)  subject to  |S| = p. S ⊆ Ω
$minimize S ⊆ Ω f(S) subject to | S | = p .$
If neither of these properties hold then nag_best_subset_given_size (h05ab) cannot be used.
As well as returning the optimal subset, So1${S}_{o1}$, nag_best_subset_given_size (h05ab) can return the best n$n$ solutions of size p$p$. If Soi${S}_{o\mathit{i}}$ denotes the i$\mathit{i}$th best subset, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,n-1$, then the (i + 1)$\left(i+1\right)$th best subset is defined as the solution to either
 maximize f(S)  subject to  |S| = p S ⊆ Ω − {Soj : j ∈ ℤ,1 ≤ j ≤ i}
$maximize S ⊆ Ω - { Soj : j∈ℤ , 1≤j≤i } f(S) subject to | S | = p$
or
 minimize f(S)  subject to  |S| = p S ⊆ Ω − {Soj : j ∈ ℤ,1 ≤ j ≤ i}
$minimize S ⊆ Ω - { Soj : j∈ℤ,1≤j≤i } f(S) subject to | S | = p$
depending on the properties of f$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 Si${S}_{i}$, can be trimmed from the tree if f(Si) < (Son) $f\left({S}_{i}\right)<\stackrel{^}{f}\left({S}_{on}\right)$ where (Son) $\stackrel{^}{f}\left({S}_{on}\right)$ is the n$n$th highest score we have observed so far for a subset of size p$p$, i.e., our current best guess of the score for the n$n$th best subset. In addition, because of (1) we can also drop all nodes defined by any subset Sj${S}_{j}$ where SjSi${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 nag_best_subset_given_size (h05ab) utilizes the fast branch and bound algorithm of Somol et al. (2004), and attempts to estimate the score where possible. For each feature, xi${x}_{i}$, two values are stored, a count ci${c}_{i}$ and μ̂i${\stackrel{^}{\mu }}_{i}$, an estimate of the contribution of that feature. An initial value of zero is used for both ci${c}_{i}$ and μ̂i${\stackrel{^}{\mu }}_{i}$. At any stage of the algorithm where both f(S) $f\left(S\right)$ and f(S{xi}) $f\left(S-\left\{{x}_{i}\right\}\right)$ have been calculated (as opposed to estimated), the estimated contribution of the feature xi${x}_{i}$ is updated to
 ( ciμ̂i + [f(S) − f(S − {xj})] )/( ci + 1 ) $ciμ^i + [ f( S ) - f( S - { xj } ) ] ci+1$
and ci${c}_{i}$ is incremented by 1$1$, therefore at each stage μ̂i${\stackrel{^}{\mu }}_{i}$ is the mean contribution of xi${x}_{i}$ observed so far and ci${c}_{i}$ is the number of observations used to calculate that mean.
As long as cik${c}_{i}\ge k$, for the user-supplied constant k$k$, then rather than calculating f(S{xi}) $f\left(S-\left\{{x}_{i}\right\}\right)$ this function estimates it using (S{xi}) = f(S) γ μ̂i $\stackrel{^}{f}\left(S-\left\{{x}_{i}\right\}\right)=f\left(S\right)-\gamma {\stackrel{^}{\mu }}_{i}$ or (S) γ μ̂i $\stackrel{^}{f}\left(S\right)-\gamma {\stackrel{^}{\mu }}_{i}$ if f(S) $f\left(S\right)$ has been estimated, where γ$\gamma$ is a user-supplied scaling factor. An estimated score is never used to trim a node or returned as the optimal score.
Setting k = 0$k=0$ in this function 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$k>0$, unless the score function is iterative in nature, i.e., f(S) $f\left(S\right)$ must have been calculated before f(S{xi}) $f\left(S-\left\{{x}_{i}\right\}\right)$ can be calculated.
nag_best_subset_given_size (h05ab) is a direct communication version of nag_best_subset_given_size_revcomm (h05aa).

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

Parameters

Compulsory Input Parameters

1:     mincr – int64int32nag_int scalar
Flag indicating whether the scoring function f$f$ is increasing or decreasing.
mincr = 1${\mathbf{mincr}}=1$
f(Si) f(Sj) $f\left({S}_{i}\right)\le f\left({S}_{j}\right)$, i.e., the subsets with the largest score will be selected.
mincr = 0${\mathbf{mincr}}=0$
f(Si) f(Sj) $f\left({S}_{i}\right)\ge f\left({S}_{j}\right)$, i.e., the subsets with the smallest score will be selected.
For all SjΩ${S}_{j}\subseteq \Omega$ and SiSj${S}_{i}\subseteq {S}_{j}$
Constraint: mincr = 0${\mathbf{mincr}}=0$ or 1$1$.
2:     m – int64int32nag_int scalar
m$m$, the number of features in the full feature set.
Constraint: m2${\mathbf{m}}\ge 2$.
3:     ip – int64int32nag_int scalar
p$p$, the number of features in the subset of interest.
Constraint: 1ipm$1\le {\mathbf{ip}}\le {\mathbf{m}}$.
4:     nbest – int64int32nag_int scalar
n$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: nbest1${\mathbf{nbest}}\ge 1$.
5:     f – function handle or string containing name of m-file
f must evaluate the scoring function f$f$.
[score, user, info] = f(m, drop, lz, z, la, a, user, info)

Input Parameters

1:     m – int64int32nag_int scalar
m = |Ω|$m=|\Omega |$, the number of features in the full feature set.
2:     drop – int64int32nag_int scalar
Flag indicating whether the intermediate subsets should be constructed by dropping features from the full set (drop = 1${\mathbf{drop}}=1$) or adding features to the empty set (drop = 0${\mathbf{drop}}=0$). See score for additional details.
3:     lz – int64int32nag_int scalar
The number of features stored in z.
4:     z(lz) – int64int32nag_int array
z(i)${\mathbf{z}}\left(\mathit{i}\right)$, for i = 1,2,,lz$\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:     la – int64int32nag_int scalar
If la > 0${\mathbf{la}}>0$, the number of subsets for which a score must be returned.
If la = 0${\mathbf{la}}=0$, the score for a single subset should be returned. See score for additional details.
6:     a(la) – int64int32nag_int array
a(j)${\mathbf{a}}\left(\mathit{j}\right)$, for j = 1,2,,la$\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:     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:     info – int64int32nag_int scalar
info = 0${\mathbf{info}}=0$.

Output Parameters

1:     score(max (la,1)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{la}},1\right)$) – double array
The value f(Sj) $f\left({S}_{\mathit{j}}\right)$, for j = 1,2,,la$\mathit{j}=1,2,\dots ,{\mathbf{la}}$, the score associated with the j$j$th subset. Sj${S}_{j}$ is constructed as follows:
drop = 1${\mathbf{drop}}=1$
Sj${S}_{j}$ is constructed by dropping the features specified in the first lz elements of z and the single feature given in a(j)${\mathbf{a}}\left(j\right)$ from the full set of features, Ω$\Omega$. The subset will therefore contain mlz1${\mathbf{m}}-{\mathbf{lz}}-1$ features.
drop = 0${\mathbf{drop}}=0$
Sj${S}_{j}$ is constructed by adding the features specified in the first lz elements of z and the single feature specified in a(j)${\mathbf{a}}\left(j\right)$ to the empty set, $\varnothing$. The subset will therefore contain lz + 1${\mathbf{lz}}+1$ features.
In both cases the individual features are referenced by the integers 1$1$ to m with 1$1$ indicating the first feature, 2$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$1$ might refer to the first variable in a particular set of data, 2$2$ the second, etc..
If la = 0${\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 lz = 0${\mathbf{lz}}=0$, this subset will either be Ω$\Omega$ or $\varnothing$.
2:     user – Any MATLAB object
3:     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:     mincnt – int64int32nag_int scalar
k$k$, the minimum number of times the effect of each feature, xi${x}_{i}$, must have been observed before f(S{xi}) $f\left(S-\left\{{x}_{i}\right\}\right)$ is estimated from f(S) $f\left(S\right)$ as opposed to being calculated directly.
If k = 0$k=0$ then f(S{xi}) $f\left(S-\left\{{x}_{i}\right\}\right)$ is never estimated. If mincnt < 0${\mathbf{mincnt}}<0$ then k$k$ is set to 1$1$.
7:     gamma – double scalar
γ$\gamma$, the scaling factor used when estimating scores. If gamma < 0${\mathbf{gamma}}<0$ then γ$\gamma$ is set to 1$1$.
8:     acc(2$2$) – double array
A measure of the accuracy of the scoring function, f$f$.
Letting ai = ε1 |f(Si)| + ε2 ${a}_{i}={\epsilon }_{1}|f\left({S}_{i}\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 Si${S}_{i}$ can be trimmed, then any values in the range f(Si) ± ai $f\left({S}_{i}\right)±{a}_{i}$ are treated as being numerically equivalent.
If 0acc(1)1$0\le {\mathbf{acc}}\left(1\right)\le 1$ then ε1 = acc(1)${\epsilon }_{1}={\mathbf{acc}}\left(1\right)$, otherwise ε1 = 0${\epsilon }_{1}=0$.
If acc(2)0${\mathbf{acc}}\left(2\right)\ge 0$ then ε2 = acc(2)${\epsilon }_{2}={\mathbf{acc}}\left(2\right)$, otherwise ε2 = 0${\epsilon }_{2}=0$.
In most situations setting both ε1${\epsilon }_{1}$ and ε2${\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:     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.

iuser ruser

Output Parameters

1:     la – int64int32nag_int scalar
The number of best subsets returned.
2:     bscore(nbest) – double array
Holds the score for the la best subsets returned in bz.
3:     bz(mip${\mathbf{m}}-{\mathbf{ip}}$,nbest) – int64int32nag_int array
The j$j$th best subset is constructed by dropping the features specified in bz(i,j)${\mathbf{bz}}\left(\mathit{i},\mathit{j}\right)$, for i = 1,2,,mip$\mathit{i}=1,2,\dots ,{\mathbf{m}}-{\mathbf{ip}}$ and j = 1,2,,la$\mathit{j}=1,2,\dots ,{\mathbf{la}}$, from the set of all features, Ω$\Omega$. The score for the j$j$th best subset is given in bscore(j)${\mathbf{bscore}}\left(j\right)$.
4:     user – Any MATLAB object
5:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

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

ifail = 11${\mathbf{ifail}}=11$
Constraint: mincr = 0${\mathbf{mincr}}=0$ or 1$1$.
ifail = 21${\mathbf{ifail}}=21$
Constraint: m2${\mathbf{m}}\ge 2$.
ifail = 31${\mathbf{ifail}}=31$
Constraint: 1ipm$1\le {\mathbf{ip}}\le {\mathbf{m}}$. Constraint: 1ipm$1\le {\mathbf{ip}}\le {\mathbf{m}}$.
ifail = 41${\mathbf{ifail}}=41$
Constraint: nbest1${\mathbf{nbest}}\ge 1$.
W ifail = 42${\mathbf{ifail}}=42$
On entry.
But only _$_$ best subsets could be calculated.
ifail = 81${\mathbf{ifail}}=81$
On exit from f, score(_) = _${\mathbf{score}}\left(_\right)=_$, which is inconsistent with the score for the parent node. On exit from f, score(_) = _${\mathbf{score}}\left(_\right)=_$, which is inconsistent with the score for the parent node.
ifail = 82${\mathbf{ifail}}=82$
A nonzero value for info has been returned:
ifail = 999${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

Accuracy

The subsets returned by nag_best_subset_given_size (h05ab) are guaranteed to be optimal up to the accuracy of the calculated scores.

The maximum number of unique subsets of size p$p$ from a set of m$m$ features is N = (m ! )/( (mp) ! p ! ) $N=\frac{m!}{\left(m-p\right)!p!}$. The efficiency of the branch and bound algorithm implemented in nag_best_subset_given_size (h05ab) comes from evaluating subsets at internal nodes of the tree, that is subsets with more than p$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$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.

Example

```function nag_best_subset_given_size_example
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 = 40;
x = [
-1.59, 0.19, 0.40, 0.43,-0.40, 0.79, 0.06, 0.33, 1.60, 0.58,-1.12, 1.23, 1.07,-0.07;
-0.25, 0.61,-0.36, 1.16, 0.61,-2.05,-0.02,-0.04, 0.80,-0.73,-0.63,-0.75,-0.73, 1.43;
-2.28, 0.46,-0.65, 0.33, 0.16,-0.21,-1.61,-0.54, 0.48, 0.37,-0.95,-2.14, 0.48, 2.02;
-0.52, 1.05, 0.64, 0.02,-1.12, 0.23, 0.06,-1.26, 1.40,-0.98, 2.47, 0.49,-0.02,-0.05;
-0.84, 1.86, 0.10, 0.73,-1.41, 0.98, 0.20,-0.89, 1.84, 2.56, 0.60,-0.12, 0.71, 0.23;
1.12,-0.51,-0.58, 0.09,-1.14, 2.11,-0.11,-0.34,-1.04,-0.43,-0.01,-0.38, 1.80, 0.05;
0.06, 0.85,-2.09, 0.22,-1.35,-0.36, 1.20, 0.41, 0.80,-0.28, 0.18, 0.27, 0.92, 0.63;
-0.48,-1.02, 0.08,-0.06, 0.13,-1.18, 2.30, 0.03, 0.45, 0.62,-1.97, 0.97, 0.93,-0.18;
0.08,-0.31, 0.43,-0.38, 0.01, 1.30, 0.66, 0.65,-0.59, 0.76, 0.04, 0.17,-0.76,-0.90;
0.66, 1.14, 0.40, 2.37, 1.10, 0.17,-0.38, 1.15,-1.00,-0.13,-0.69,-0.62,-0.18, 0.00;
-1.08,-0.21,-1.13,-0.79,-0.76,-1.58, 0.38,-0.03, 1.26,-0.51,-0.75, 0.86, 0.29, 0.68;
-0.74,-1.59,-0.58,-1.09, 1.18,-1.70,-1.02, 0.36, 1.05, 1.30,-0.98,-1.36,-1.28,-1.32;
0.40,-1.58,-1.30,-0.10,-1.34, 0.65,-0.56, 0.39,-0.73,-0.32, 2.19,-0.49, 0.69, 0.18;
0.75,-3.09,-0.61,-1.89, 0.15, 0.77,-0.49,-0.63, 1.20,-0.04, 1.02, 0.31, 0.81,-0.45;
-0.65, 1.57,-1.50,-1.45, 0.21, 0.06, 0.24, 2.24,-1.34, 0.30, 1.39,-0.38,-0.71, 0.48;
1.36, 1.40,-1.40,-0.90, 0.36,-0.21,-0.97, 0.36,-0.26, 0.08, 0.06,-1.49, 0.43,-1.61;
-1.01, 1.50,-0.61,-0.25,-1.01,-0.43, 1.90,-1.33,-0.96,-0.02, 0.51,-1.38,-0.78, 1.82;
1.34, 1.02, 3.50, 0.10, 0.50, 0.04, 0.61,-0.57,-2.69,-0.64,-0.34,-0.21,-1.97,-0.19;
0.29, 0.67,-0.38,-0.63,-0.24, 1.21,-0.09, 0.90,-2.20, 1.72, 0.29, 0.66, 0.19,-0.57;
0.67,-0.56,-0.41, 1.22,-0.30, 0.77, 0.82, 0.36, 1.18, 1.87,-1.48, 0.52, 1.35, 0.13;
-0.40,-1.10,-0.83, 0.71, 1.99,-0.24, 1.30,-0.34,-0.70, 0.28, 0.16, 0.27, 0.37,-1.79;
-0.78, 0.60,-0.45,-0.26,-0.23, 0.89, 0.87, 1.01, 1.20, 0.28, 0.79, 2.76, 0.35, 1.31;
-1.29, 0.62,-0.59, 1.52, 0.62, 0.21, 1.31, 1.09,-0.36,-0.34,-0.03,-0.59,-1.70,-0.03;
0.40,-1.45,-0.98, 2.10,-1.09,-0.53,-0.38,-1.36, 0.13, 0.70,-1.51, 0.08,-0.62,-0.64;
0.43,-0.86, 0.70,-1.07,-0.76, 0.72,-0.14,-1.58, 0.00, 0.58,-0.21, 1.30, 2.02, 1.52;
-0.48, 0.01, 1.30, 0.58,-0.54, 1.09, 0.91, 2.90, 1.32,-1.20,-0.59,-0.51, 0.20,-1.74;
-1.32,-1.41,-0.58,-1.29, 1.61,-0.35,-0.72,-1.92,-1.09, 0.56,-0.87,-0.71, 1.25, 0.10;
1.43, 0.69, 1.34,-0.32, 2.84,-1.43,-0.47,-0.01, 0.83,-0.72,-0.78, 0.50,-1.22, 0.54;
0.82, 0.46, 0.15,-0.57, 0.93, 1.33,-0.23,-1.07, 0.76, 0.25,-1.96, 0.39, 0.24,-0.26;
-0.91, 0.23,-0.19, 1.58,-0.27, 0.33,-0.60,-1.39,-0.30,-0.81,-0.95, 0.88,-0.09,-0.35;
0.65,-1.14, 1.18,-1.06,-0.68,-0.22, 0.21, 0.94, 1.08, 0.81,-0.33, 0.42,-0.90, 0.49;
-0.36,-0.50,-0.02,-0.04, 0.77, 0.62,-1.35,-0.64, 1.20, 1.22, 0.18,-1.39,-0.81,-0.99;
-1.82, 1.06, 0.28, 0.14, 0.62,-0.80,-1.08,-2.15, 1.37, 1.57,-1.48,-0.79, 0.28,-0.20;
1.54, 0.50, 0.13,-0.68, 0.26,-1.13, 0.62,-0.43, 0.39, 1.14, 0.15, 1.03, 0.46, 0.40;
-1.61,-0.61, 0.93,-0.37, 0.44,-1.45, 0.58,-1.77, 0.72,-2.05,-0.03,-1.24,-1.40,-0.06;
-0.48, 0.67, 0.04, 0.27,-0.84,-0.06,-3.67, 0.09, 1.66,-0.30, 1.67, 1.08, 0.00, 0.43;
-1.65,-1.16,-1.17, 1.12, 0.11,-0.15, 0.48,-1.72, 1.08,-0.94, 0.49,-0.56, 0.95, 1.09;
-0.85,-0.02, 1.18,-1.16, 0.49, 1.56,-0.60, 0.32, 0.72,-1.20, 2.52, 1.78, 0.16,-0.01;
-0.60,-0.73,-1.23, 1.50, 0.40,-0.20,-0.65, 0.68, 1.09, 0.40,-1.50,-2.10, 0.21,-0.18;
-0.66,-0.01,-0.01, 0.85,-2.04, 1.17,-0.56, 1.72,-0.18, 1.14,-0.96,-0.92,-0.28, 1.58];
y = [-2.44;
-2.97;
7.42;
3.00;
8.83;
0.03;
2.57;
8.31;
4.55;
-23.10;
3.38;
-0.13;
5.47;
13.97;
20.94;
-12.87;
28.26;
6.89;
5.37;
-1.50;
-23.00;
14.09;
-11.05;
-32.04;
23.36;
-5.58;
2.48;
-5.30;
-7.77;
-34.25;
26.78;
-11.85;
-8.62;
12.35;
-1.54;
-16.59;
-8.69;
7.82;
-18.56;
17.21];

% 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] = ...
nag_best_subset_given_size(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. nag_best_subset_given_size_revcomm 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
for j = 1:mip
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 = nag_correg_linregm_fit('z', user{1}, isx, ip, user{2});

% Return the score (the residual sums of squares)
end

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

Score        Feature Subset
-----        --------------
1.04753e+03     4     7     8    10    14
1.05987e+03     4     5     7     8    14
1.07016e+03     4     5     7    10    14

45 subsets evaluated in total

```
```function h05ab_example
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 = 40;
x = [
-1.59, 0.19, 0.40, 0.43,-0.40, 0.79, 0.06, 0.33, 1.60, 0.58,-1.12, 1.23, 1.07,-0.07;
-0.25, 0.61,-0.36, 1.16, 0.61,-2.05,-0.02,-0.04, 0.80,-0.73,-0.63,-0.75,-0.73, 1.43;
-2.28, 0.46,-0.65, 0.33, 0.16,-0.21,-1.61,-0.54, 0.48, 0.37,-0.95,-2.14, 0.48, 2.02;
-0.52, 1.05, 0.64, 0.02,-1.12, 0.23, 0.06,-1.26, 1.40,-0.98, 2.47, 0.49,-0.02,-0.05;
-0.84, 1.86, 0.10, 0.73,-1.41, 0.98, 0.20,-0.89, 1.84, 2.56, 0.60,-0.12, 0.71, 0.23;
1.12,-0.51,-0.58, 0.09,-1.14, 2.11,-0.11,-0.34,-1.04,-0.43,-0.01,-0.38, 1.80, 0.05;
0.06, 0.85,-2.09, 0.22,-1.35,-0.36, 1.20, 0.41, 0.80,-0.28, 0.18, 0.27, 0.92, 0.63;
-0.48,-1.02, 0.08,-0.06, 0.13,-1.18, 2.30, 0.03, 0.45, 0.62,-1.97, 0.97, 0.93,-0.18;
0.08,-0.31, 0.43,-0.38, 0.01, 1.30, 0.66, 0.65,-0.59, 0.76, 0.04, 0.17,-0.76,-0.90;
0.66, 1.14, 0.40, 2.37, 1.10, 0.17,-0.38, 1.15,-1.00,-0.13,-0.69,-0.62,-0.18, 0.00;
-1.08,-0.21,-1.13,-0.79,-0.76,-1.58, 0.38,-0.03, 1.26,-0.51,-0.75, 0.86, 0.29, 0.68;
-0.74,-1.59,-0.58,-1.09, 1.18,-1.70,-1.02, 0.36, 1.05, 1.30,-0.98,-1.36,-1.28,-1.32;
0.40,-1.58,-1.30,-0.10,-1.34, 0.65,-0.56, 0.39,-0.73,-0.32, 2.19,-0.49, 0.69, 0.18;
0.75,-3.09,-0.61,-1.89, 0.15, 0.77,-0.49,-0.63, 1.20,-0.04, 1.02, 0.31, 0.81,-0.45;
-0.65, 1.57,-1.50,-1.45, 0.21, 0.06, 0.24, 2.24,-1.34, 0.30, 1.39,-0.38,-0.71, 0.48;
1.36, 1.40,-1.40,-0.90, 0.36,-0.21,-0.97, 0.36,-0.26, 0.08, 0.06,-1.49, 0.43,-1.61;
-1.01, 1.50,-0.61,-0.25,-1.01,-0.43, 1.90,-1.33,-0.96,-0.02, 0.51,-1.38,-0.78, 1.82;
1.34, 1.02, 3.50, 0.10, 0.50, 0.04, 0.61,-0.57,-2.69,-0.64,-0.34,-0.21,-1.97,-0.19;
0.29, 0.67,-0.38,-0.63,-0.24, 1.21,-0.09, 0.90,-2.20, 1.72, 0.29, 0.66, 0.19,-0.57;
0.67,-0.56,-0.41, 1.22,-0.30, 0.77, 0.82, 0.36, 1.18, 1.87,-1.48, 0.52, 1.35, 0.13;
-0.40,-1.10,-0.83, 0.71, 1.99,-0.24, 1.30,-0.34,-0.70, 0.28, 0.16, 0.27, 0.37,-1.79;
-0.78, 0.60,-0.45,-0.26,-0.23, 0.89, 0.87, 1.01, 1.20, 0.28, 0.79, 2.76, 0.35, 1.31;
-1.29, 0.62,-0.59, 1.52, 0.62, 0.21, 1.31, 1.09,-0.36,-0.34,-0.03,-0.59,-1.70,-0.03;
0.40,-1.45,-0.98, 2.10,-1.09,-0.53,-0.38,-1.36, 0.13, 0.70,-1.51, 0.08,-0.62,-0.64;
0.43,-0.86, 0.70,-1.07,-0.76, 0.72,-0.14,-1.58, 0.00, 0.58,-0.21, 1.30, 2.02, 1.52;
-0.48, 0.01, 1.30, 0.58,-0.54, 1.09, 0.91, 2.90, 1.32,-1.20,-0.59,-0.51, 0.20,-1.74;
-1.32,-1.41,-0.58,-1.29, 1.61,-0.35,-0.72,-1.92,-1.09, 0.56,-0.87,-0.71, 1.25, 0.10;
1.43, 0.69, 1.34,-0.32, 2.84,-1.43,-0.47,-0.01, 0.83,-0.72,-0.78, 0.50,-1.22, 0.54;
0.82, 0.46, 0.15,-0.57, 0.93, 1.33,-0.23,-1.07, 0.76, 0.25,-1.96, 0.39, 0.24,-0.26;
-0.91, 0.23,-0.19, 1.58,-0.27, 0.33,-0.60,-1.39,-0.30,-0.81,-0.95, 0.88,-0.09,-0.35;
0.65,-1.14, 1.18,-1.06,-0.68,-0.22, 0.21, 0.94, 1.08, 0.81,-0.33, 0.42,-0.90, 0.49;
-0.36,-0.50,-0.02,-0.04, 0.77, 0.62,-1.35,-0.64, 1.20, 1.22, 0.18,-1.39,-0.81,-0.99;
-1.82, 1.06, 0.28, 0.14, 0.62,-0.80,-1.08,-2.15, 1.37, 1.57,-1.48,-0.79, 0.28,-0.20;
1.54, 0.50, 0.13,-0.68, 0.26,-1.13, 0.62,-0.43, 0.39, 1.14, 0.15, 1.03, 0.46, 0.40;
-1.61,-0.61, 0.93,-0.37, 0.44,-1.45, 0.58,-1.77, 0.72,-2.05,-0.03,-1.24,-1.40,-0.06;
-0.48, 0.67, 0.04, 0.27,-0.84,-0.06,-3.67, 0.09, 1.66,-0.30, 1.67, 1.08, 0.00, 0.43;
-1.65,-1.16,-1.17, 1.12, 0.11,-0.15, 0.48,-1.72, 1.08,-0.94, 0.49,-0.56, 0.95, 1.09;
-0.85,-0.02, 1.18,-1.16, 0.49, 1.56,-0.60, 0.32, 0.72,-1.20, 2.52, 1.78, 0.16,-0.01;
-0.60,-0.73,-1.23, 1.50, 0.40,-0.20,-0.65, 0.68, 1.09, 0.40,-1.50,-2.10, 0.21,-0.18;
-0.66,-0.01,-0.01, 0.85,-2.04, 1.17,-0.56, 1.72,-0.18, 1.14,-0.96,-0.92,-0.28, 1.58];
y = [-2.44;
-2.97;
7.42;
3.00;
8.83;
0.03;
2.57;
8.31;
4.55;
-23.10;
3.38;
-0.13;
5.47;
13.97;
20.94;
-12.87;
28.26;
6.89;
5.37;
-1.50;
-23.00;
14.09;
-11.05;
-32.04;
23.36;
-5.58;
2.48;
-5.30;
-7.77;
-34.25;
26.78;
-11.85;
-8.62;
12.35;
-1.54;
-16.59;
-8.69;
7.82;
-18.56;
17.21];

% 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
for j = 1:mip
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)
end

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

Score        Feature Subset
-----        --------------
1.04753e+03     4     7     8    10    14
1.05987e+03     4     5     7     8    14
1.07016e+03     4     5     7    10    14

45 subsets evaluated in total

```