

#Householder i LZNK

#zad 1 test czy H Houseloldera ortogonalna
w=randn(10,1);
a=2/(w'*w);
H=eye(10)- a*w*w';

norm(H-H',1)
norm(H*H-eye(size(H)),1)

b=0;
for k=1:100,
 x=randn(size(w));
 if(abs(norm(x)-norm(H*x))>1e-14)
   b++;
   printf("Differecne of norms>1e-14\n");
 endif
endfor
b

#zad 2 test wlasnosci
n=10;
u=randn(10,1);
v=randn(10,1);
v=(norm(u)/norm(v))*v;
w=u-v;
kw=2/(w'*w);

y=u-(kw*(w'*u))*w;
norm(y-v)


z=v-(kw*(w'*v))*w;
norm(z-u)


#zad 3 metoda householdera krok po kroku dla 
#regresji liniowej tzn dane (xk,yk) k=1,..,m
#szukamy a,b takich \sum_k |a*xk+b-yk|^2 minimalna

m=10;
x=linspace(0,1,m)';
y=2*x+1 +0.1*randn(m,1);

A=[x,ones(m,1)];

function [rozw,er,R,w1,w2]=lznkhous2k(A,y)
#funkcja rozwiazujaca LZNK z macierza 
#2 kolumnowa A ui wektorem y metoda householdera
#INput:
#A - macierz 2 kolumnowa (kolumny niezlzne liniowo)
#y - wektor prawej strony w LZNK
#OUTPUT
#rozw wektor dl 2 - rozwiazanie LZNK z A i y
#er = blad LZNK tzn ||A*rozw-y||_2
#R - macierz 2x2 - 2 pierwsze wiersze macierzy QR1
#    z rozkladu QR  zQ ortogonalna mxm i R1 m x 2 
#    otrzymanego m Householdera tzn R=R1(1:2,:)
#w1,w2 - wektory Householdera dla odp macierzy Householdera 
#       tzn w1 wektor Hous. dla H1 a w2 dla H2 z 
#       H2H1A=R1 czyli Q=H1*H2       
w1=A(:,1); n1=norm(w1,2); 
if(w1(1) > 0)
 printf("H1A(1,1) niedodatnie: %f \n",w1(1));
 w1(1)=w1(1)-n1;
else
 printf("H1A(1,1) dodatnie: %f \n",w1(1));
 w1(1)=w1(1)+n1;
endif 
k1=2/(w1'*w1);
#to nieoptymalne mnozymy przez H1 pierwsza kolumne
#a wiemy ze wyjdzie n1 (lub -n1) razy wersor e1
H1A=A-w1*(k1*(w1'*A));

w2=H1A(:,2);w2(1)=0;
n2=norm(w2,2);
if(w2(2)<=0)
 printf("H1A(2,2) niedodatnie: %g \n",w2(2));
 w2(2)=w2(2)-n2;
else
 printf("H1A(2,2) dodatnie: %g \n",w2(2));
 w2(2)=w2(2)+n2;
endif
k2=2/(w2'*w2);
H2H1A=H1A-w2*(k2*(w2'*H1A));
R=H2H1A(1:2,:);
H1y=y-w1*(k1*(w1'*y));
g=H2H1y=H1y-w2*(k2*(w2'*H1y));

rozw=R\g(1:2);
er = norm(g(3:length(g)),2);

endfunction

#wywolujemy 
[rozw,er]=lznkhous2k(A,y)
#porownujemy z \ 
r1=A\y
norm(A*r1-y)





