

%lab 3 animal breeding
disp("animal breeding");

load rozd5data.txt

%direct solver
n=length(A);m=rows(X);

tic; G=[X'*X,X'*Z;Z'*X,Z'*Z+k*inv(A)];
r=[X'*y;Z'*y];
s1=G\r;b1=s1(1:2);a1=s1(3:end)
t1=toc

%2nd: solve LLSP with [X,Z,0,S] with S=sqrt(k)A^(-0.5)
tic;
[Q,D]=eig(A);
S=sqrt(k)*Q*diag(1./sqrt(diag(D)))*Q';
B2=[X,Z;zeros(n,2),S];
s2=B2\[y;zeros(n,1)];
b2=s2(1:2);a2=s2(3:end)
t2=toc
er2=norm(a2-a1)

%3rd way: solve LLSP with [X,Z,0,S] with S=sqrt(k)inv(L) with L'*L=A so S'*S=k*A^(-1)
tic;
L=chol(A)';
S=sqrt(k)*inv(L);
B3=[X,Z;zeros(n,2),S];
s3=B3\[y;zeros(n,1)];
b3=s3(1:2);a3=s3(3:end)
t3=toc
er3=norm(a3-a1)




%4th way solve LLSP with [X,Z*L;0,sqrt(k)*Id] and [y;0]
%getting [b;g] with g=L^(-1)a i.e. a=L*g L'L=A
tic;
L=chol(A)';
S=sqrt(k)*eye(size(A));
B4=[X,Z*L;zeros(n,2),S];
s4=B4\[y;zeros(n,1)];
b4=s4(1:2);a4=L*s4(3:end)
t4=toc
er4=norm(a4-a1)
t=[t1,t2,t3,t4]
er=[er2,er3,er4]