function [B,sset,ops]=b3msv(Q,nc)
%B3MSV Bidirectional Branch and Bound with Minimum Singular Value Criterion for Subset Selection
%
% Consider the following subset selection problem:
% Given a tall (m x n, m>n) matrix, A, to find n rows of A
% such that the resulted n x n square submatrix has the
% largest MSV among all possible n x n submatrices.
% This problem has many applications, where one wishes to
% square down a non-square matrix to get a well-posted problem
% (non-singular, hence maximizing the MSV).
% This problem has been studied in Linear Algebra for many decades.
% Many intuitive solutions have been proposed either analytically or
% numerically. But none of them, except two approaches, the
% exhaustive search and the branch and bound (BB) can guarantee the global
% optimality.
% Exhaustive search can only be used for very small m and n.
% Traditional BB is unidirectional, either upwards, where the subset
% is gragually expending until reaching the desired size, or
% downwards, where subset is shrinking one by one until the target
% size. The performance of both are very limited.
% In B3, the search is carried out in both directions hence is much
% more efficient than unidirectional BB. Moreover, a novel
% determinant based pruning algorithm is implemented to replace
% time-consuming singular value decomposition so that the overall
% efficincy is about several orders of magnitude faster than
% unidirectional approaches.
%Usage:
% [B,subset,ops]=b3msv(A,r)
%Input:
% A, m x n matrix, m>n (row selection) or n>m (column selection)
% r, number of best subsets required, default is 1
%Output:
% B, a r-length vector contains the best MSVs
% subset, a r x n index matrix for the r best n-element subsets
% ops, various counters for total number of different evaluations
% for algorithm developers.
%
% See also b3wc, pb3wc
% Reference
% Y. Cao and V. Kariwala, Bidirectional Branch and Bound for Controlled
% Variable Selection Part I: Principles and Minimum Singular Value
% Criterion, Computers and Chemical Engineering, 32(2008), 2306-2319
%
% By Yi Cao at Cranfield University, 07 November 2007,
% revised on 09/01/2009 for large size problems
%Examples:
%Ex1 check correctness
%{
A = randn(10,5); % Consider a 5 out of 10 problem:
V = nchoosek(1:10,5); % All 252 combinations
B = zeros(252,1);
for k=1:252 % exhaustive search
B(k) = min(svd(A(V(k,:),:)));
end
[B1,k]=sort(B,'descend'); % the solution
[B2,subset]=b3msv(A,10);
norm(B1(1:10)-B2) % should be zero or very small
isequal(subset,V(k(1:10),:)) % should be 1
%}
%Ex2 check efficiency
%{
A = randn(40,20); % 20 out of 40 problem
tic, b3msv(A); toc % 7.5 s with 2.0GHz Core Duo T2500 processor
tic,
for k=1:100, % test time required for MSV
B=min(svd(A(1:20,:)));
end
t1=toc % 0.2 seconds
year=nchoosek(40,20)*t1/100/3600/24/365 % over 2 years!
%}
%Ex3 check special cases
%{
A = randn(400,5); % 5 out of 400
tic, b3msv(A); toc % 0.9 seconds
A = randn(400,395); % 395 out of 400
tic, b3msv(A); toc % 11.4 seconds
%}
%Ex4 a very large case
%{
n = 100000;
m = 20;
A = 1./randn(n,m);
tic,
[B,sset,op] = b3msv(A); % 15 + 142 nodes evaluated
toc % about 3 seconds %
%}
%Error checking
if nargin > 2
error('Too many input arguments');
end
% if call without input to display help
if ~nargin
eval(['help ',eval('mfilename')])
return
end
if nargin<2, % number of best subsets
nc=1;
end
[r,n]=size(Q);
if r==n % trivial case
B=min(svd(Q));
sset=1:n;
ops=1;
return
end
if r<n % column selection
Q=Q';
[r,n]=size(Q);
end
bi=2;
Q2=Q.*Q; % preparation
q2=sum(Q2,2);
h2=q2';
g2=h2;
B=zeros(nc,1);
sset=zeros(nc,n);
ib=1;
bound=0;
ops=[0 0 0 0 0]; % ops count for research
fx=false(1,r); % fixed set
rem=true(1,r); % candidate set
nf=0; % size of fixed set
n2=n; % n2 = n - nf
m2=r; % size of candidate set
upf=false; % upwards pruning flag
downf=false; % downward pruning flag
f=false; % fail flag
bf=false; % bound upodate flag
bbmsvsub(fx,rem); % main recursive solver
[B,idx]=sort(sqrt(B),'descend'); % post processing
sset=sort(sset(idx,:),2);
function bbmsvsub(fx0,rem0)
ops(4)=ops(4)+1;
fx=fx0;
rem=rem0;
nf=sum(fx);
m2=sum(rem);
n2=n-nf;
% bf=true;
while ~f && m2>=n2 && n2>=0 % Pruning loops
while ~f && m2 >= n2 && n2>=0 && (~downf || ~upf || ~bf)
if ~upf || ~bf
upprune; % upwards pruning
end
if ~f && m2 >= n2 && (~downf || ~bf)
downprune; % downward pruning
end
bf=true;
end %pruning loop
if f || m2<=n2 || n2<=0
break % invalid or trivial problems
end
if n2==1 % select one problems
% not very efficient for some problems
% if bound
% h2(~rem)=Inf;
% [b,idk]=min(h2);
% else
g2(~rem)=0;
[b,idk]=max(g2);
% end
s=fx;
s(idk)=true;
update(s);
downf=false;
rem(idk)=false;
m2=m2-1;
continue
end
if m2-n2==1, % discard one problems
h2(~rem)=0;
[b,idk]=max(h2);
rem(idk)=false;
s=fx|rem;
update(s);
fx(idk)=true;
upf=false;
nf=nf+1;
n2=n2-1;
m2=m2-1;
continue
end
fx1=fx;
rem1=rem;
b0=bound;
h0=h2;
g0=g2;
if (n2<0.5*m2 && bi) || bi==1 %upward branching
if bound
h2(~rem)=Inf;
[b,id]=min(h2);
else
g2(~rem)=0;
[b,id]=max(g2);
end
rem1(id)=false;
fx(id)=true;
upf=false;
bbmsvsub(fx,rem1);
upf=true;
downf=false;
rem=rem1;
fx=fx1;
else %downward branching
if bound && n<100
g2(~rem)=Inf;
[b,id]=min(g2);
else
h2(~rem)=0;
[b,id]=max(h2);
end
rem1(id)=false;
fx1(id)=true;
downf=false;
bbmsvsub(fx,rem1);
if g0(id)<=bound
return
end
downf=true;
upf=false;
fx=fx1;
rem=rem1;
end
if b0<bound
bf=false;
L=min(g0,q2')<=bound;
if any(L)
rem=rem&~L;
downf=false;
end
end
f=false;
h2=h0;
g2=g0;
nf=sum(fx);
n2=n-nf;
m2=sum(rem);
end
if ~f % update bound for trivial problems
if m2==n2
update(fx|rem);
elseif ~n2
update(fx);
end
end
end
function upprune
f=false;
if ~nf
g2(rem)=q2(rem);
elseif nf==1
if q2(fx)<bound
f=true;
return;
end
x=Q(rem,:)*Q(fx,:)';
g2(rem)=(x.*x)/(bound-q2(fx))+q2(rem);
else
Z=Q(fx,:)*Q(fx,:)';
[R,f]=chol(Z-bound*eye(nf));
if f,
return;
end
if m2>=n
X=(R'\Q(fx,:))*Q(rem,:)';
else
X=R'\(Q(fx,:)*Q(rem,:)');
end
g2(rem)=q2(rem)-sum(X.*X,1)';
end
ops(2)=ops(2)+1;
ops(3)=ops(3)+m2;
g2(~rem)=Inf;
upf=true;
L=g2<=bound;
if any(L)
downf=false;
rem=rem&~L;
m2=sum(rem);
end
end
function downprune
s=fx|rem;
h2(:)=Inf;
[R,f]=chol(Q(s,:)'*Q(s,:)-bound*eye(n));
if f
return
end
G=R'\Q(rem,:)';
h2(rem)=1-sum(G.*G,1);
ops(2)=ops(2)+1;
ops(3)=ops(3)+m2;
downf=true;
L=h2<0;
if any(L)
upf=false;
fx=fx | L;
rem=rem & ~L;
nf=sum(fx);
m2=sum(rem);
n2=n-nf;
end
end
function bf0=update(s) % update bound
lambda=eig(Q(s,:)*Q(s,:)');
ops(1)=ops(1)+1;
bf0=sum(lambda<bound);
if ~bf0
B(ib)=min(lambda); %avoid sorting
sset(ib,:)=find(s);
[bound,ib]=min(B);
bf=false;
end
end
end
