

A=[1, 1; 1, 0; 1, -1]
[Q,R]=qr(A)
Q*R-A
norm(Q*R-A,1)

norm(A*x-b)
s=randn(2,1);b=A*s;
x=A\b
norm(x-s)
b=randn(3,1);
x=A\b
y=R(1:2,:)\(Q'*b)(1:2)
norm(x-y)
z=(A'*A)\(A'*b)
norm(x-z)

%zad 2 porownac 3 metody rozwiazania LZNK

n=3;
MM=6;
er=zeros(MM,4);
for k=1:MM,
  m=2*n;
  H=hilb(m);
  A=H(:,1:n);
  b=A(:,1);
  x=A\b;
  y=(A'*A)\(A'*b);
  B=[eye(m,m),A;A',zeros(n,n)];
  z=B\[b;zeros(n,1)];
  [x,y,zz=z((m+1):end)]
  s=zeros(n,1);s(1)=1;
  er(k,:)=[m,norm(x-s),norm(y-s),norm(zz-s)];
  n++;
endfor
er

%regularyzacja zamist LZNK z A i b 
%rozwiazujemy zregularyzowane LZNK 
%xe=arg min (||Ax-b||_2^2+eps||x||_2^2)
m=5;n=m-2;;
B=randn(m,n);
%B=single(B);
A=[B,2*B];
b=rand(m,1);
%b=single(b);
MM=7;
eps=0.1;
s=A\b;
er=zeros(MM,3);
for k=1:MM,
 Ae=[A;eps*eye(columns(A))];
 bb=[b;zeros(columns(A),1)];
 xe=Ae\bb;
 ye=(A'*A+eps*eps*eye(columns(A)))\(A'*b);
 [s,xe,ye]
 er(k,:)=[eps,norm(xe-s),norm(ye-s)];
 eps*=0.1;
end 
er


%znalezc krzywa zadana rowaniem postaci ax^2+by^2-1=0
%najlepiej pasujaca do zadanych punktow z plaszczyzny
%tzn (a,b)=arg min \sum_k|alx(k)^2+be y(k)^2-1|^2

function ab=krzywapas(x,y)
%znalezc krzywa zadana rowaniem postaci ax^2+by^2-1=0
%najlepiej pasujaca do zadanych punktow z plaszczyzny
%tzn (a,b)=arg min \sum_k|alx(k)^2+be y(k)^2-1|^2

  A=[x.*x,y.*y];
  f=ones(length(x),1);
  ab=A\f;
endfunction

%punkty na elpisie 0.5x^2+y^2=1
x=rand(10,1);
y=sqrt(1-0.5*x.*x);
ab=krzywapas(x,y)

y1=y+0.01*randn(size(y));
ab=krzywapas(x,y1)
ms=linspace(-1,1,200);
plot(ms,sqrt((1-ab(1)*ms.*ms))/ab(2),x,y1,"+")


%funkcja mnozaca macierz przez H m Householdera

function Y=HtX(X,v,kd=v'*v)
  Y=X- ((2/kd)*v)*(v'*X);
endfunction

v=randn(4,1);
x=randn(4,1);
y=HtX(x,v);
w=HtX(y,v);
norm(w-x)
norm(y)-norm(x)
A=randn(4,4);
HA=HtX(A,v);
norm(HA)-norm(A)
y=HtX(x,v);

%testujemy tw z wykladu
u=x; 
l=norm(v)/norm(u);
v-HtX(l*u,v-l*u)
u=zeros(size(v));
u(1)=norm(v);
HtX(v,v-u)

%mnozymy przez H m Householdera wektor lub macierz A na 3 sposoby

N=10;M=4;
t=zeros(M,3);
for k=1:M,
  v=rand(N,1);X=randn(N,4);
 
  tic; H=eye(N)-(2/(v'*v)*v)*v';H*X;t(k,1)=toc;
  tic;H*X;t(k,2)=toc;
  tic;X- (((2/kd)*v)*v')*X;t(k,3)=toc;
  tic;X- ((2/kd)*v)*(v'*X);t(k,4)=toc;
  N*=10;
endfor
