
%porownanie 4 metod rozwiazania LZNK dla macierzy 3x2
A=[1, 1; 1, 0; 1, -1]
b=A(:,1)
[Q,R]=qr(A)
norm(Q*R-A,1)
Q'*Q
s=[1;0];
g=Q'*b;
x=R(1:2,:)\g(1:2)
y=A\b
B=A'*A;
z=B\(A'*b)
%zad 3 rozwiazujemy uklad [aI,A;A',0][r;x]=[b;0] a>0
a=1;B=[a*eye(3),A;A',zeros(2,2)];
w=B\[b;zeros(2,1)];
r=w(1:3)
w=w(4:5)

b=randn(3,1);
g=Q'*b;
x=R(1:2,:)\g(1:2)
er=abs(g(3))
y=A\b
norm(A*y-b)
B=A'*A;
z=B\(A'*b)
%zad 3 rozwiazujemy uklad [aI,A;A',0][r;x]=[b;0] a>0
a=1;B=[a*eye(3),A;A',zeros(2,2)];
w=B\[b;zeros(2,1)];
r=w(1:3)
w=w(4:5)

%testy \ (QR) i ukl r normalnych 
%dla podmacierzy macierzy Hilberta

S=6;
N=3;
er=zeros(S,7);
for k=1:S,
  H=hilb(2*N);
  H=H(:,1:N);
  b=H(:,1);
  s=zeros(N,1);s(1)=1;
  x=H\b; %rozwiazujemy poprzez rozklad QR
  y=(H'*H)\(H'*b);%zad 3 rozwiazujemy uklad [aI,A;A',0][r;x]=[b;0] a>0
  a=1; %mzona wziac rozne a np a=1000 lu 1e-4...
  B=[a*eye(2*N),H;H',zeros(N,N)];
  w=B\[b;zeros(N,1)];
  r=w(1:2*N);
  w=w(2*N+1:end);
  er(k,:)=[N,norm(x-s),norm(H*x-b),norm(y-s),norm(H*y-b),norm(w-s),a*norm(r)];
  N+=2;
endfor
er


%nieregularne LZNK z A i b regularyzujemy tzn
%szukamy xe=arg min (||Ax-b||_2^2 +||eps x||_2^2)
%zregularyzowane LZNK to LZNk z macierza 
%[A;eps I] i [b;0] ma uklad r normalnych (A'A+eps*I)xe=A'b
%rozwiazemy to zreg LZNk na 2 sposoby (\ czyli QR i ukl r normalnych)
%i porownamy z rozwiazaniem LZNK z A i b otrzymanego backslashem
%dla A=[B,2*B] B losowa

N=20;
M=12;
eps=1e-2;
er=zeros(M,5);
for k=1:M,
  B=rand(N,N);
  A=[B,2*B];
  b=rand(N,1);
  x=A\b;
  Ae=[A;eps*eye(2*N)];be=[b;zeros(2*N,1)];
  xe=Ae\be;
  ye=(C=A'*A+(eps*eps)*eye(2*N))\(A'*b);
  er(k,:)=[eps,norm(x-xe),cond(Ae),norm(ye-x),cond(C)];
  eps*=0.1;
endfor
er


function Y=HouseTX(X,v,knv=v'*v)
  %mnozymy macierz X Mx NRHS przez macierz Householdera z lewej strony
  %majac tylko wektor Househ. tej macierzy 
 Y=X-v*((2/knv)*(v'*X)); 
endfunction

v=randn(5,1);
x=randn(5,1);
Hx=HouseTX(x,v);
norm(x)-norm(Hx)
norm(x-HouseTX(Hx,v))
X=randn(5,4);
HX=HouseTX(X,v);
norm(X)-norm(HX)
norm(X,"fro")-norm(HX,"fro")
%sprawdzmy tw z wykladu u<>w 
%||u||_2=||w||_2 to m Householdera z wektrem v=u-w 
%przeprowadza u na v a v na u
v=v/norm(v);
x=x/norm(x);
w=x-v;
x-HouseTX(v,w)
v-HouseTX(x,w)

A=randn(5,3);
w=A(:,1);
w(1)=w(1)-norm(w);
HouseTX(A(:,1),w)
HA=HouseTX(A,w)
norm(A-HouseTX(HA,w),1)


%porownanie czasu mnozenie Hx dla H m Householdera 
%zadanej przez wektor Hous. w
%3 sposoby: tworzymy H i przez nia mnozymy (nie wliczamy czasu utworzenia)
%2/ ze wzoru Hx=x-(2/(knw))*(w*w')*x
%3/ ze wzoru Hx=x-(2/(knw))*w*(w'*x) (tylko inaczej postawione nawiasy)

N=10;
M=4;
time=zeros(M,4);
for k=1:M,
  w=rand(N,1);
  x=rand(N,50);
  knw=w'*w;
  H=eye(N)-(2/knw)*w*w';
  time(k,1)=N;
  tic;H*x;time(k,2)=toc;
  tic;x-(2/knw)*(w*w')*x;time(k,3)=toc;
  tic;x-w*((2/knw)*(w'*x));time(k,4)=toc;
  N*=10;
endfor
time