Odwracanie macierzy

foster

Algorytm do odwracania macierzy przedstawiony poniżej pochodzi z książki Donalda E. Knuth'a "Sztuka programowania" tom 1. Odwraca on macierz bardzo efektywnie nie korzystając przy tym ze znanego z matematyki wzoru wykorzystującego wyznacznik macierzy i dopełnienia algebraiczne.

Algorytm zastepujący macierz jej macierzą odwrotną przy założeniu, że elementami macierzy są A[i,j] dla 1<=i,j<=n :

1) Dla k=1,2,...,n wykonuj co następuje: przejrzyj wszytkie kolumny wiersza k nie wykorzystane do tej pory jako kolumny osiowe w poszukiwaniu elementu o największej wartości bezwzglednej. Przypisz do C[k] numer kolumny, w której znaleziono taki element, i wykonaj operacje osiowania, przyjmując ten element za element osiowy. Jeśli wszystkie takie elementy są zerami, to macierz jest osobliwa i nie można jej odwrócić.)
2) Dokonaj permutacji wierszy i kolumn w ten sposob, żeby to, co było wierszem k, stało sie wierszem C[k], a to co było kolumną C[k], stało się kolumną k.

 Operacja osiowania przebiega w następujący sposób:
                   Przed osiowaniem         Po osiowaniu

                  kolumna     inna       kolumna     inna
                  osiowa     kolumna     osiowa     kolumna
                /   .           .   \ /    .          .     \
                |   .           .   | |    .          .     |
                |                   | |                     |
  wiersz osiowy |.. a ......... b ..| |.. 1/a ...... b/a ...|
                |                   | |                     |
                |   .           .   | |    .          .     |
                |   .           .   | |    .          .     |
                |   .           .   | |    .          .     |
                |                   | |                     |
                |.. c ......... d ..| |..-c/a .... d-bc/a ..|
  inny wiersz   |                   | |                     |
                |   .           .   | |    .          .     |
                \   .           .   / \    .          .     /

Poniżej przedstawiam kod źródłowy programu realizujacego powyższy algorytm w języku Pascal.

{$N+}
Program Macierz_Odwrotna;
Const
  Max  = 50;      { Maksymalny stopien macierzy }
  Zero = 1.0E-7;  { Umowne zero }
Type
  TMacierz= Array[1..Max,1..Max] of Real;
 
Function OdwracanieMacierzy(Var Macierz: TMacierz;Stopien:Byte): Boolean;
Var
  w,k,w1,k1,i,                { zmienne wykorzystywane przy iteracji }
  p,q              : Integer;  { wspolrzedne elementu osiowego }
  c               : Array[1..Max] of Byte;
  BylaJuzKolOsiowa: Array[1..Max]of boolean;
  Max             : Real;
  MacPom          : TMacierz;
begin
FillChar(BylaJuzKolOsiowa,SizeOf(BylaJuzKolOsiowa),False);
for w:=1 to Stopien do
  begin
    c[w]:=0;
    Max:=0.0;
    for k:=1 to Stopien do
      begin
        if not(BylaJuzKolOsiowa[k]) then
        if Max < Abs(Macierz[w,k]) then
          begin
            Max:=Abs(Macierz[w,k]);
            p:=w;
            q:=k;
          end;
      end;
    if Max<Zero then
      begin
        OdwracanieMacierzy:= False;
        Exit;
      end;
    for w1:=1 to Stopien do
      for k1:=1 to Stopien do
        if (w1<>p)and(k1<>q)then
          Macierz[w1,k1]:= Macierz[w1,k1]-(Macierz[p,k1]*Macierz[w1,q]/Macierz[p,q]);
    for w1:=1 to Stopien do
      if w1<>p then Macierz[w1,q]:= -Macierz[w1,q]/Macierz[p,q];
    for k1:=1 to Stopien do
      if k1<>q then Macierz[p,k1]:= Macierz[p,k1]/Macierz[p,q];
    Macierz[p,q]:= 1/Macierz[p,q];
 
    c[w]:= q;
    BylaJuzKolOsiowa[q]:= True;
  end;
 
(* Permutowanie macierzy *)
MacPom:=Macierz;
for k:=1 to Stopien do
  for i:=1 to Stopien do
    Macierz[c[k],i]:= MacPom[k,i];
 
MacPom:=Macierz;
for k:=1 to Stopien do
  for i:=1 to Stopien do
    Macierz[i,k]:= MacPom[i,c[k]];
 
OdwracanieMacierzy:= True;
end;
 
Procedure Wczytaj(var Macierz:TMacierz;var Stopien:Byte);
Var
  f   : Text;
  i,j : Integer;
Begin
  Assign(f,'.\in.txt');
  Reset(f);
  ReadLn(f,Stopien);
  if Stopien>max then
    begin
      WriteLn('Macierz za duza!!!');
      Halt;
    end;
  for i:=1 to Stopien do
    begin
      for j:=1 to Stopien do
        Read(f,Macierz[i,j]);
      ReadLn(f);
    end;
  Close(f);
end;
 
Procedure Zapisz(Var Macierz:TMacierz; Stopien:Integer; OK:Boolean);
Var
  f   : Text;
  i,j : Integer;
Begin
  Assign(f,'.\out.txt');
  Rewrite(f);
  if OK then
    for i:=1 to Stopien do
      begin
        for j:=1 to Stopien do
          Write(f,Macierz[i,j],' ');
        WriteLn(f);
      end
    else
    WriteLn(f,'Macierz osobliwa!!!');
  Close(f);
end;
 
Var
  Macierz : TMacierz;
  Stopien : Byte;
Begin
  Wczytaj(Macierz,Stopien);
  if OdwracanieMacierzy(Macierz,Stopien) then
    Zapisz(Macierz,Stopien,True)
    else
    Zapisz(Macierz,Stopien,False);
end.

Program można ulepszyć pisząc dobra funkcję permutowania wierszy i kolumn, która nie będzie wykorzystywała macierzy pomocniczej. Stopień macierzy, którą można będzie odwrócić zwiększy się wtedy dwukrotnie (czyli do 100).


Metoda eliminacji Gaussa-Jordana

Dane wejsciowe 
Dwie macierze -jedna odwracana,druga jednostkowa
Na obu macierzach wykonujemy te same operacje
Gdy pierwsza macierz sprowadzimy do macierzy jednostkowej 
druga przyjmie postac macierzy odwrotnej

1. Zamiana wierszy nie wplywa na macierz odwrotna
2. Dodanie wierszy nie wplywa na macierz odwrotna
3. Pomnozenie wierszy przez skalar rozny od zera nie wplywa na macierz odwrotna

 
uses
  Crt;
 
{ Rozmiary macierzy. }
 
const
  rozn = 20;
  rozb = 40;
 
{ Macierz }
 
type
  macierz = array [1..rozn,1..rozb] of real;
 
{ dane:                                             }
{   A - macierz do odwr˘cenia (pierwszych n kolumn) }
{   n - rozmiar macierzy                            }
{ wynik:                                            }
{   A - kolumny od n + 1 do 2n                      }
 
procedure odwmac (var a: macierz; n: integer);
var
  i,j,k,m,wm: Integer;
  s,t,w,p,max: Real;
  b:macierz;
begin
for i:=1 to n do
for j:=1 to n do
 b[i,j]:=a[i,j];
        w:= 1;
      for i:= 1 to n-1 do
         begin
            max:= abs (b[i,i]);
            wm:= i;
            for j:= i+1 to n do         { wybor maxymalnego elementu }
               if abs (b[j,i]) > max then
                  begin
                     max:= abs (b[j,i]);
                     wm:= j
                  end;
            if wm > i then              { zamiana wierszy }
               begin
                  for k:= i to n do
                     begin
                        p:= b[i,k];
                        b[i,k]:= b[wm,k];
                        b[wm,k]:= p
                     end;
                  w:= -w                { zamiana wierszy zmienia znak }
               end;
            if b[i,i] = 0 then          { macierz osobliwa }
               w:= 0;
            if w <> 0 then              { eliminacja Gaussa }
               for j:= i+1 to n do
                  begin
                     p:= b[j,i] / b[i,i];
                     for k:= n downto i+1 do
                        b[j,k]:= b[j,k] - p * b[i,k]
                  end;
            w:= w * b[i,i]              { iloczyn elementow na przekatnej }
         end;
      w:= w * b[n,n];
      if w=0 then
       begin
       writeln('Macierz osobliwa');
       exit;
       end;
  for i := 1 to n do
    for j := n + 1 to 2 * n do
      if i = j - n then
        A [i,j] := 1
      else
        A [i,j] := 0;
    for i := 1 to n do
      begin
             max:= abs (a[i,i]);
            wm:= i;
            for j:= i+1 to n do         { wybor maxymalnego elementu }
               if abs (a[j,i]) > max then
                  begin
                     max:= abs (a[j,i]);
                     wm:= j
                  end;
            if wm > i then              { zamiana wierszy }
               begin
                  for k:= i to 2*n do
                     begin
                        p:= a[i,k];
                        a[i,k]:= a[wm,k];
                        a[wm,k]:= p
                     end;
               end;
        s := A [i,i];
        for j := 1 to 2 * n do
          A [i,j] := A [i,j] / s;
        for j := i + 1 to n do
          begin
            t := A [j,i];
            for k := i to 2 * n do
              A [j,k] := - A [i,k] * t + A [j,k];
          end;
      end;
    for k := n + 1 to 2 * n do
      begin
        for i := n downto 1 do
          begin
            s := A [i,k];
            for j := i + 1 to n do
              s := s - A [i,j] * A [j,k];
            A [i,k] := s;
          end;
      end;
end;
 
var
  A: macierz;
  s,t: Real;
  j,i,n: Integer;
 
begin
  ClrScr;
  WriteLn ('Odwracanie macierzy');
  WriteLn;
  Write ('Podaj rozmiar macierzy n=');
  ReadLn (n);
  for i := 1 to n  do
    begin
      WriteLn ('Wprowad« ',i,'-ty wiersz macierzy');
      for j := 1 to n do
        ReadLn  (A [i,j]);
    end;
  odwmac (A,n);
  for i := 1 to n do
    begin
      for j := n + 1 to 2 * n do
        Write (A [i,j]:4:8,' ');
      WriteLn;
    end;
  ReadKey;
end.
 

6 komentarzy

Uszanowanie dla autora. Algorytm oczywiście działa. Przeliczyłem nie macierz a potem sprawdziłem dwoma sposobami. Mam prośbę, moglibyście wyjaśnić dlaczego to działa?
Próbowałem z książki Knuth-a ale nie zrozumiałem jaki jest związek z "przeosiowaniem macierzy" a jej macierzą odwrotną.Pozdrawiam serdecznie

Nawet fajny algorytm. Tylko ciekawy jestem czy jest on zawsze stabilny. Metoda Gaussa-Jordana moze nie jest najszybsza ale zawsze doprowadzi nas do poprawnego wyniku. Oczywiscie gdy macierz jest nie osobliwa :)

Ładnie, nie czytałem tekstu tylko kod więc jeśli to Twój algorytm to bardzo ładny. Chłopaki - bez przesady Guass był matematykiem a foster jest nie ujmując mu ambitnym programistą, więc o co chodzi ;)

Gauss na pewno jest prostszy, a co za tym idzie, w praktyce moze okazac sie szybszy, zwlaszcza dla mniejszych macierzy (stale pomijane przy szacowaniu zlozonosci, czesto nie sa takie male jakby sie wydawac moglo). A jeszcze lepiej Gaussa-Jordana...
Ale swoja droga, to ciekawa metoda :)

A czy przypadkiem algorytm Gaussa nie jest szybszy?

  1. Tag < delphi >
  2. Formatowanie kodu (wszystkie słowa kluczowe - function, procedure, var, begin - małą literą i wcięcia)
  3. Ten "rysunek" na początku się ciut rozpadł