最小二乘法及带约束条件的最小二乘法(空间约束和L2约束)

clc;clear;
n = 50; N = 1000;
x = linspace(-3,3,n)'; X = linspace(-3,3,N)';
pix = pi*x;
y = sin(pix)./(pix) + 0.1*x + 0.05*randn(n,1);
%%
%3.2
% p(:,1) = ones(n,1);P(:,1)= ones(N,1);
% for j= 1:15
%     p(:,2*j) = sin(j/2*x);p(:,2*j+1) = cos(j/2*x);
%     P(:,2*j) = sin(j/2*X);P(:,2*j+1) = cos(j/2*X);
% end
% t = p\y;
% F = P*t;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -0.5 1.2]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%3.8
% hh = 2*0.3^2; t0 = randn(n,1); e = 0.1;
% for o = 1:n*1000
%     i = ceil(rand*n);
%     ki = exp(-(x-x(i)).^2/hh);
%     t = t0 - e*ki*(ki'*t0-y(i));
%     if norm(t - t0) < 0.000001,break,end
%     t0 = t;
% end
% K = exp(-(repmat(X.^2,1,n) + repmat(x.^2',N,1) - 2*X*x')/hh);
% F = K*t;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -0.5 1.2]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%4.3
% y = sin(pix)./(pix) + 0.1*x + 0.05*randn(n,1);
% p(:,1) = ones(n,1);P(:,1)= ones(N,1);
% for j= 1:15
%     p(:,2*j) = sin(j/2*x);p(:,2*j+1) = cos(j/2*x);
%     P(:,2*j) = sin(j/2*X);P(:,2*j+1) = cos(j/2*X);
% end
% t1 = p\y;
% F1 = p*t1;
% t2 = (p*diag([ones(1,11) zeros(1,20)]))\y;
% F2 = P*t2;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -0.8 1.2]);
% plot(x,F1,'g-');plot(X,F2,'r--');plot(x,y,'bo');
% legend('LS','Subspace-constrained LS');
%%
%4.7
% y = sin(pix)./(pix) + 0.1*x + 0.2*randn(n,1);%
% 
% x2 = x.^2;X2=X.^2;hh = 2*0.3^2;l=0.1;
% 
% k = exp(-(repmat(x2,1,n)+repmat(x2',n,1)-2*x*x')/hh);
% K = exp(-(repmat(X2,1,n)+repmat(x2',N,1)-2*X*x')/hh);
% t1 = k\y;
% F1=K*t1;
% t2=(k^2+l*eye(n))\(k*y);
% F2 = K*t2;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -1 1.5]);
% plot(X,F1,'g-');plot(X,F2,'r--');plot(x,y,'bo');
% legend('LS','LS-constrained LS');
%%
%4.14
% y = sin(pix)./(pix) + 0.1*x + 0.2*randn(n,1);%
% 
% x2 = x.^2;xx=repmat(x2,1,n)+repmat(x2',n,1)-2*x*x';
% hhs=2*[0.03 0.3 3].^2;ls=[0.0001 0.1 100];
% m=5;u=floor(m*[0:n-1]/n)+1;u=u(randperm(n));
% for hk=1:length(hhs)
%     hh = hhs(hk);k=exp(-xx/hh);
%     for i=1:m
%         ki = k(u~=i,:);kc=k(u==i,:);yi=y(u~=i);yc=y(u==i);
%         for lk = 1:length(ls)
%             l=ls(lk);t=(ki'*ki*+l*eye(n))\(ki'*yi);fc=kc*t;
%             g(hk,lk,i) = mean((fc-yc).^2);
%         end,end,end
% [gl,ggl] = min(mean(g,3),[],2);[gh1l,gghl] = min(gl);
% L=ls(ggl(gghl));HH = hhs(gghl);
% 
% K = exp(-(repmat(X.^2,1,n)+repmat(x2',N,1)-2*X*x')/HH);
% k = exp(-xx/HH);t = (k^2+L*eye(n))\(k*y);F=K*t;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -0.7 1.7]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%5.8 奇异矩阵行列式 = 0；
% hh = 2*0.3^2; l = 0.1;t0 = randn(n,1); x2 = x.^2;
% k = exp(-(repmat(x2,1,n) + repmat(x2',n,1) - 2*x*x')/hh);
% k2 = k^2;ky = k*y;
% for o = 1:1000
%     t = (k2+1*pinv(diag(abs(t0))))\ky;%diag：取对角元素[1,n];pinv:奇异矩阵或者方阵求伪逆；
%inv:非奇异，则存在逆矩阵；
%     if norm(t - t0) < 0.001,break,end
%     t0 = t;
% end
% K = exp(-(repmat(X.^2,1,n) + repmat(x2',N,1) - 2*X*x')/hh);
% F = K*t;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -1 1.5]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%6.9
% n = 10;N = 1000;x = linspace(-3,3,n)';X = linspace(-4,4,N)';
% y = x + 0.2*randn(n,1);y(n) = -4;
% 
% p(:,1) = ones(n,1);p(:,2) = x;t0 = p\y;e = 1;
% for o = 1:1000
%     r = abs(p*t0-y);w = ones(n,1);w(r>e) = e./r(r>e);
%     t = (p' * (repmat(w,1,2).*p))\(p'*(w.*y));
%     if norm(t - t0) < 0.001,break,end
%     t0 = t;
% end
% P(:,1) = ones(N,1);P(:,2) = X;F = P*t;
% 
% figure(1);clf;hold on;axis([-4 4 -4.5 3.5]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%6.14
% n = 50;N = 1000;x = linspace(-3,3,n)';X = linspace(-3,3,N)';
% pix = pi*x;y = sin(pix)./(pix) + 0.1*x + 0.2*randn(n,1);
% y(n/2) = -5;
% 
% hh = 2*0.3^2; l = 0.1;e = 0.1;t0 = randn(n,1); x2 = x.^2;
% k = exp(-(repmat(x2,1,n) + repmat(x2',n,1) - 2*x*x')/hh);
% for o = 1:1000
%     r = abs(k*t0-y);w = ones(n,1);w(r>e) = e./r(r>e);
%     Z = k*(repmat(w,1,n).*k)+l*pinv(diag(abs(t0)));
%     t = (Z + 0.000001*eye(n))\(k*(w.*y));
%     if norm(t - t0) < 0.001,break,end %norm 二范数
%     t0 = t;
% end
% K = exp(-(repmat(X.^2,1,n) + repmat(x2',N,1) - 2*X*x')/hh);
% F = K*t;
% 
% figure(1);clf;hold on;axis([-2.8 2.8 -1 1.5]);
% plot(X,F,'g-');plot(x,y,'bo');
%%
%7.4
% clear;
% n = 200;a = linspace(0,4*pi,n/2);
% u = [a.*cos(a) (a+pi).*cos(a)]'+1*rand(n,1);%看做两维特征；
% v = [a.*sin(a) (a+pi).*sin(a)]'+1*rand(n,1);
% x = [u v];y = [ones(1,n/2) -ones(1,n/2)]';
% 
% x2 = sum(x.^2,2);hh = 2*1^2;l = 0.01;
% k = exp(-(repmat(x2,1,n)+ repmat(x2',n,1) -2*x*x')/hh);
% t = (k^2 +l*eye(n))\(k*y);
% 
% m = 100;X = linspace(-15,15,m)';X2 = X.^2;
% U = exp(-(repmat(u.^2,1,m)+ repmat(X2',n,1) - 2*u*X')/hh);
% V = exp(-(repmat(v.^2,1,m)+ repmat(X2',n,1) - 2*v*X')/hh);
% figure(1);clf;hold on;axis([-15 15 -15 15]);
% contourf(X,X,sign(V'*(U.*repmat(t,1,m))));
% plot(x(y == 1,1),x(y == 1,2),'bo');
% plot(x(y == -1,1),x(y == -1,2),'rx');
% colormap([1 0.7 1; 0.7 1 1]);
%%
%8.20
% clear;
n = 40;x = [randn(1,n/2)-15 randn(1,n/2)-5;randn(1,n)]';
y = [ones(n/2,1);-ones(n/2,1)];
x(1:2,1) = x(1:2,1)+60;x(:,3) = 1;

l = 0.01;e = 0.01;t0 = zeros(3,1);
for o = 1:1000
        m = (x*t0).*y; v = m + min(1,max(0,1-m));
        a = abs(v-m);w = ones(size(y));w(a>e) = e./a(a>e);
        t = (x'*(repmat(w,1,3).*x)+l*eye(3))\(x'*(w.*v.*y));
        if norm(t - t0) < 0.001,break,end
        t0 = t;
end
figure(1);clf;hold on;z = [-20 50];axis([z -2 2]);
plot(x(y==1,1),x(y==1,2),'bo');
plot(x(y == -1,1),x(y==-1,2),'rx');
plot(z,-(t(3)+z*t(1))/t(2),'k-');