function Blatt11

% Matrix A
A = [1 2 3 ; 4 5 6; 7 8 9; 10 0 0];
% Test
Aplus = Moore_Penrose(A);
% Moore-Penrose Axiome überprüfen - via Frobeniusnorm
disp(norm(A*Aplus*A -A,'fro'))
disp(norm(Aplus*A*Aplus -Aplus,'fro'))
disp(norm((A*Aplus)' -A*Aplus,'fro'))
disp(norm((Aplus*A)' -Aplus*A,'fro'))

end

function Aplus = Moore_Penrose(A)
    % Funktion, welche die Moore-Penrose-Inverse berechnet. 
    % Dazu ist nur erlaubt, die Funktion eig zu verwenden.
    % svd ist nicht erlaubt. 
    % Daher berechnen wir zuerst die Singulärwertzerlegung von A*A (es genügt die economy sized SVD zu berechnen, 
    % vgl. was nötig ist für die Moore-Penrose Inverse (MPI) und help svd).
    % Für die MPI reicht es, dass \Sigma^+ eine pxp Matrix ist, V eine nxp
    % und U eine mxp Matrix, denn die anderen Spalten von V und U sind
    % wegen der Nullblöcken in der vollen \Sigma^+-Matrix unnötig.
    AstarA = A'*A;
    [V,EW] = eig(AstarA); % die Eigenvektoren sind gerade die ONB V der SVD
    % filtere die nichtnull Eigenwerte
    [EW,Ind] = sort(diag(EW),1,'descend'); % Diagonale als Vektor schreiben und absteigend sortieren
    % aus dem Vektor die Elemente > 1e-12 herausfiltern, d.h. EW > 0.
    EW = EW(EW > 1e-12); 
    % nun die Singulärwerte
    SV = sqrt(EW);
    p = length(SV);
    % V sortieren - in Ind ist die Sortierung festgehalten
    Vsort = V(:,Ind);
    Veco = Vsort(:,1:p);
    % Konstruiere damit die Matrix U (economy version - die Spalten zu den 0 Eigenwerten werden nicht benötigt, da diese in der späteren Matrixmultiplikation herausfallen)
    for i = 1:p
        Ueco(:,i) = 1/SV(i)*A*Veco(:,i); % entsprechend der Sortierung U aufstellen
    end
    % Konstruiere Sigma^+
    Sigmaplus = zeros(p);
    Sigmaplus(1:p,1:p) = diag(1 ./ SV);
    % Und damit A^+
    Aplus = Vsort*Sigmaplus*Ueco';
end

