%% Demonstrace výpočtu singulárního rozkladu

A=[1 2 3;4 5 6;7 8 9;10 11 12]

%
% Najdeme diagonalizaci matice C=A^TA: V_1^TCV_1=D, na diagonále D jsou druhé mocniny singulárních šísel.

C=A'*A;
[V1,D,ni]=QRA_sym_diag(C)
[V2,D2]=eig(A'*A)

%
% Na matici AV_1 aplikujeme QR rozklad s výběrem vedoucích sloupců, prvky na diagonále R volíme nezáporné.

[Q,R,P]=qr(A*V1);
s=find(diag(R)<0);
if length(s)>0 % v R na diagonale jsou zaporna cisla
    R(s,:)=-R(s,:);
    Q(:,s)=-Q(:,s);
end % otoceni znamenek radku R a sloupcu Q

%
% Matice R je diagonální, v praxi počítáme se zaokrouhlovacími chybami.
% V_1P=V, U=Q, R=\Sigma, U^TAV=\Sigma je singulární rozklad matice A.

V=V1*P;
U=Q;
Sigma=diag(diag(R));

% Matici \Sigma musíme doplnit nulami na stejný rozměr jako má R.

[m,n]=size(R);
if m~=n
    Sigma(m,n)=0;
end

% Výsledek

U,Sigma,V

%
% Test výsledku

U*Sigma*V'-A

% Porovnání se standardním výsledkem Matlabu

[U1,S1,V1]=svd(A)