how to solve not enough argument on the bold one?
1 view (last 30 days)
Show older comments
function [B,vset,ops,ct] = bbdireg(Y,X,n,nc)
% BBDIREG:
% Improved downwards branch and bound (BBDi) to select subset for least
% squares regression problems, Y = X*p, where X is independent, Y is
% dependent and p is parameter.
% [B,vset,ops,ct] = bbdireg(y,x,n,nc) finds the best n-variable subset of
% X to minimize J = 0.5*(Y-A*p)'*(Y-A*p)
% Inputs:
% y: each row is a dependent smaple, m columns
% x: each row is an independent sample, r columns
% n: subset size to be selected
% nc: the number of top solutions to be selected
% Outputs:
% B: The minimum J of nc top solutions
% vset: subset indices of nc top solutions
% ops: total number of operations spent.
% ct: total computation time spent
%
% By Yi Cao (c) at Cranfield University 17th February 2013
%
% Example
%{
r = 41;
m = 2;
n = 20;
N = 1000;
Y = randn(N,m);
X = randn(N,r);
nc = 10;
[B,vset,ops,ct] = bbdireg(Y,X,n,nc);
%}
% Reference: Kariwala, V., Ye, L. and Cao, Y, "Branch and Bound Method for
% Regression-Based Controlled Variable Selection", Computers & Chemical
% Engineering, to appear 2013
%
% determin size of the problem
[N,r]=size(X);
if size(Y,1)~=N
error('Y and X should have the same number of rows.')
end
% convert to standard problem J = 0.5*x'*inv(C)x;
C = X'*X; % should be rxr
invCo=inv(C);
cmin=min(eig(C));
Z = X'*Y;
bound=0;
fx=false(1,r);
rem=true(1,r);
crem=rem;
B=zeros(nc,1);
vset=zeros(nc,n);
ib=1;
ops=zeros(1,5);
invC=invCo;
p2=zeros(1,r);
y2=p2;
nf=0;
n2=n;
m2=r;
f=false;
downf=false;
t0=cputime;
bbm2sub(fx,rem);
[B,idx]=sort(B);
B = 0.5*(Y(:)'*Y(:) - B);
ct=cputime-t0;
vset=sort(vset(idx,:),2);
% Main recursive function
function bbm2sub(fx0,rem0)
ops(4)=ops(4)+1;
fx=fx0;
rem=rem0;
nf=sum(fx);
m2=sum(rem);
n2=n-nf;
while ~f && m2>=n2 && n2>=0
if ~f && n2>0 && m2>n2 && ~downf %&& (bound>0 || m2-n2==1)
downprune; % downward pruning
else
downf=true; % no need to check pruning conditions
end
if f || m2<n2 || n2<0
return % no feasible solution
elseif (m2==n2 && n2>1) || (~n2 && (m2-n2>1 || ~downf))
break % no need to branch further
end
% downwards branching
fx1=fx;
rem1=rem;
b0=bound;
p0=p2;
invC0=invC;
p2(~rem)=Inf;
[~,id]=min(p2);
fx1(id)=true;
rem1(id)=false;
downf=false;
if sum(rem1)+sum(fx)>n
updateC(fx|rem1,id); % update C inverse required
end
bbm2sub(fx,rem1); % recursive call
downf=any(p0);
invC=invC0;
fx=fx1;
rem=rem1;
if b0<bound % bound improved check pruning
p0(~rem)=Inf;
L=p0<=bound;
if any(L)
rem=rem&~L;
downf=true;
fx=fx|L;
end
end
f=false;
p2=p0;
nf=sum(fx);
n2=n-nf;
m2=sum(rem);
end
if ~f % terminal node update
if m2==n2
update(fx|rem);
elseif ~n2
update(fx);
end
end
end
% Update C inverse
function updateC(s,id)
u = invC(s,id);
d = invC(id,id);
invC(s,s) = invC(s,s) - u*(d\u');
end
% Downwards pruning function
function downprune
while ~downf
f=false;
downf=true;
ops(2)=ops(2)+1; % operation count
ops(3)=ops(3)+m2;
s=fx|rem; % fix + candidate set
rem0=rem;
D=invC(s,s);
y=Z(s,:);
dy=D*y;
b0=dy(:)'*y(:);
if b0 < bound % current node pruning
f = true;
return;
end
dy2 = sum(dy.*dy,2);
d=diag(D);
alpha = (dy2./d)';
p2(s)=b0-alpha; % subnode pruning
if m2-n2==1 % special case, discard one
ops(5)=ops(5)+1;
p2(~rem)=0;
[b,id]=max(p2);
while b>bound % find the optimal solution
s0=s;
s0(id)=false;
B(ib)=b; %avoid sorting
vset(ib,:)=find(s0);
[bound,ib]=min(B);
p2(id)=0;
[b,id]=max(p2);
end
f=true; % terminated branch serach
return;
end
if bound>0 && m2-n2>1
p2(~rem)=Inf;
L=p2<bound; % subnode pruning
if any(L)
rem = rem & ~L;
fx = fx | L;
m2 = sum(rem);
nf = sum(fx);
n2 = n - nf;
if n2<=0
return
end
end
end
rems=rem(s);
if any(crem~=rem)
cmin = 1/max(eig(D(rems,rems))); % update cmin
crem=rem;
end
[sy,idx]=sort(dy2(rems));
ind = find(rem);
idx = ind(idx);
pred=sum(sy(1:(m2-n2)))*cmin; % predictive downwards pruning
if b0-pred<bound
f=true;
return
end
if bound>0 && m2-n2>1 % predictive subnode pruning downwards
y2=p2;
y2(idx(m2-n2+1:m2))=b0-pred-(sy(m2-n2+1:m2)-sy(m2-n2))*cmin;
L=y2<bound;
L(idx(1:m2-n2))=false;
L(~rem)=false;
if any(L)
rem = rem & ~L;
fx = fx | L;
m2 = sum(rem);
nf = sum(fx);
n2 = n - nf;
if n2<=0
return
end
end
end
if bound>0 && n2>0 % predictive subnode pruning upwards
y2(idx(1:m2-n2))=(sy(1:m2-n2)-sy(m2-n2+1))*cmin;
L0=b0-pred+y2<=bound;
L0(idx(m2-n2+1:m2))=false;
L0(~rem0)=false;
if any(L0(L))
f=true;
return
end
L0(~rem)=false;
if any(L0)
downf=false;
rem(L0)=false;
m2=sum(rem);
if m2<=n2
return
end
updateC(rem|fx,L0);
end
end
if downf
return
end
end
end
% Update function
function update(s)
z=Z(s,:)'/chol(C(s,s));
ops(1)=ops(1)+1;
lambda=z(:)'*z(:);
if lambda>=bound
B(ib)=lambda; %avoid sorting
vset(ib,:)=find(s);
[bound,ib]=min(B);
end
end
end
0 Comments
Answers (1)
Cris LaPierre
on 18 Jan 2021
Edited: Cris LaPierre
on 18 Jan 2021
You can delete the bolded line.
[N,r]=size(X);
It is not necessary. N and r are defined about 10-15 row above, and those are the values used to create X in the first place.
X = randn(N,r);
X, N and r are not modified between these lines, so the bolded line of code is not needed.
4 Comments
See Also
Categories
Find more on Regression in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!