function x =LARS(y, D, N)
[nd, md] = size(D);
x = zeros(md, 1);
d = ones(md, 1);
for i = 1:md
d(i) = norm(D(:, i));
if d(i) > 0
D(:, i) = D(:, i)/d(i);
else
d(i) = 1;
end
end
activeSet = zeros(N, 1);
s = zeros(md, 1);
flags = zeros(md, 1);
c = D'*y;
absc = abs(c);
[C, j] = max(absc);
nj = 1;
activeSet(nj) = j;
s(j) = sign(c(j));
flags(j) = 1;
for i = 1:N
if nj == 0
break
end
Xa = ones(nd, 1)*s(activeSet(1:nj))'.*D(:, activeSet(1:nj));
Ga = Xa'*Xa;
e = ones(nj, 1);
Aa = (e'/Ga*e)^-0.5;
w = Aa*Ga\e;
a = D'*Xa*w;
rs = zeros(2*md, 1);
for g = 1:md
if (flags(g) == 1)
continue;
end
rs(2*g-1) = (C-c(g))/(Aa-a(g));
rs(2*g) = (C+c(g))/(Aa+a(g));
end
r = min(rs(rs>0));
if isempty(r)
break
end
j = find(rs==r,1);
j = floor((j+1)/2);
wa = zeros(md, 1);
wa(activeSet(1:nj)) = s(activeSet(1:nj)).*w;
rj = -x./wa;
rh = min(rj(rj>0));
if (~isempty(rh) && rh<=r)
C = C-rh*Aa;
jj = find(rj==rh, 1);
x = x + rh*wa;
flags(jj) = 0;
jj = find(activeSet==jj);
activeSet(jj:nj-1) = activeSet(jj+1:nj);
nj = nj-1;
else
if (r*Aa < 1e-8)
break
end
C = C-r*Aa;
x = x + r*wa;
nj = nj+1;
activeSet(nj) = j;
flags(j) = 1;
end
c = D'*(y-D*x);
s(j) = sign(c(j));
d1 = norm(y - D*x);
if (C<1e-6 || d1 < 1e-6)
break
end
end
x = x./d;
uj5u.com熱心網友回復：
再某些資料下，d1先減少后變大，x的值會突然變大，所以我加了if (r*Aa < 1e-8)
break
end
但是效果不是很好。
uj5u.com熱心網友回復：
help！頂轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/172542.html
標籤：其他開發語言
上一篇：回圈一條條處理資料，不報錯，程式突然像死機，錯誤不可復現
下一篇：python 兩表合并