

1==1;

%LZNK - met num 2016/17

%zad 1 
%testy metod rozwiazywania LZNK
A=[2 1 2;1 -1 3]' %jakkakolwiek macierz
s=[1;1];f=A*s;
%backslash - czarna skrzynak octave'a
x=A\f
norm(x-s)
%qr
[q,r]=qr(A);[norm(q'*q-eye(3),1),norm(A-q*r,1)] 
g=q'*f;y=r(1:2,:)\g(1:2)
norm(y-s)
%uklad r. nomrlalnych
z=(A'*A)\(A'*f)
norm(z-s)

%powtrzamy testy metod dla n-1 kolimn macierzy hilberta n x n
ee=[];
for n=5:20,
 H=hilb(n);H=H(:,2:n);
 s=ones(n-1,1);f=H*s;
%backslash - czarna skrzynak octave'a
x=H\f; e(1,1:2)=[norm(x-s), norm(H*x-f)];
%qr
[q,r]=qr(H);
[norm(q'*q-eye(n),1),norm(H-q*r,1)]
g=q'*f;y=r(1:n-1,:)\g(1:n-1);
e(2,1:2)=[norm(y-s), norm(H*y-f)];
%uklad r. nomrlalnych
z=(H'*H)\(H'*f);
e(3,1:2)=[norm(z-s), norm(H*z-f)];
e
ee=[ee;e];
end

%regresja liniowa dane pkty (xk,yk)
%szukamy a i b prostej takiej ze minimalizuje
%\sum_k|ax(k)+b-y(k)|^2
%czyli regresja liniowa - LZNK z A=[x,ones(size(x))] i f=y
%x,y - wektory pionowe 
%testujemy dla pktow bliskich danej prostej np y=2*x-7+eps
%z eps wektorem z zaburzenaimi jednostajnymi na poziomie 1e-2
disp("regresja liniowa");
x=linspace(0,2,10)';
y=2*x-7+2*(rand(size(x))-0.5);
A=[x,ones(size(x))];
z=A\y
norm(z-[2;-7])
s=linspace(-.5,2.5);
plot(x,y,"+",s,z(1)*s+z(2));pause(1);
%zamieniamy x z y tzn szukamy
%c,d t ze minimalizuje
%\sum_k |cy(k)+d-x(k)|^2
A=[y,ones(size(y))];
w=A\x
ymax=max(z(1)*s+z(2));ymin=min(z(1)*s+z(2));
sy=linspace(ymin,ymax);
plot(x,y,"+",s,z(1)*s+z(2),w(1)*sy+w(2),sy)

%testy Householdera
disp("testy macierzy Householdera");
n=10;w=randn(n,1);
H=eye(n)-(2/(w'*w))*w*w';
norm(H'-H,1)
norm(H*H-eye(n),1)
x=randn(n,1);
norm(H*x)-norm(x)

%LZNk z macierza rzedu 1 z wykorzystaniem p Householdera
disp("LZNK z macierza rzedu jeden");
m=10;n=4;
u=randn(m,1);v=randn(n,1);
f=randn(m,1);
A=u*v';
rank(A) %czy rzedu jeden
w1=u;w1(1)=w1(1)-norm(u)
w2=v;w2(1)=w2(1)-norm(v)
%sprawdzamy
H1u=u-((2/(w1'*w1))*(w1'*u))*w1
H2v=v-((2/(w2'*w2))*(w2'*v))*w2
H1=eye(m)-(2/(w1'*w1))*w1*w1';
H2=eye(n)-(2/(w2'*w2))*w2*w2';
H1AH2=H1*A*H2%to tylko do testow....
norm(u)*norm(v)-H1AH2(1,1)
g=f-((2/(w1'*w1))*(w1'*f))*w1;
y=zeros(n,1);
y(1)=g(1)/(norm(u)*norm(v))
x1=H2*y
%a naprawde nie potrzebujemy ani H2 ani w2
x=v*(g(1)/(norm(u)*(v'*v)))
norm(x-x1)
x2=A\f
norm(x-x2)





