Polar decomposition software

Unit Wiskunde;
{
This software has been designed and is CopyLefted
by Han de Bruijn:
                          (===)
                         @-O^O-@
                          #/_\#
                           ###

Polar decomposition matrix manipulation routines
================================================
}
INTERFACE

type
  matrix = array of array of double;

procedure Inverse(var B : matrix); { B := B^(-1) }
procedure Kwadraat(B : matrix; var P: matrix); { P := B^T*B }
procedure Produkt(A,B : matrix; var AB: matrix); { AB := A*B }
procedure Kopie(B : matrix; var P : matrix); { P := B }
procedure Drukaf(Q : matrix); { Print Q }
procedure Eenheid(NDM : integer; var E : matrix); { E := unity }
procedure Newton(P : matrix; var Q : matrix); { Q := sqrt(P) }
procedure Algemeen(B : matrix); { Testing in general }

IMPLEMENTATION

procedure Inverse(var B : matrix);
{
  Matrix inversion
  pivots on diagonal
}
var
  pe : double;
  NDM,i,j,k : integer;
begin
  NDM := Length(B);
  for k := 0 to NDM-1 do
  begin
    pe := B[k,k];
    B[k,k] := 1;
    for i := 0 to NDM-1 do
    begin
      B[i,k] := B[i,k]/pe;
      if (i = k) then Continue;
      for j := 0 to NDM-1 do
      begin
        if (j = k) then Continue;
        B[i,j] := B[i,j] - B[i,k]*B[k,j];
      end;
    end;
    for j := 0 to NDM-1 do
      B[k,j] := - B[k,j]/pe;
    B[k,k] := 1/pe;
  end;
end;

procedure Kwadraat(B : matrix; var P: matrix);
{
  P = B^T.B
}
var
  NDM,i,j,k : integer;
  h : double;
begin
  NDM := Length(B);
  SetLength(P,NDM,NDM);
{ Transpose * Original }
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      h := 0;
      for k := 0 TO NDM-1 do
        h := h + B[k,i]*B[k,j];
      P[i,j] := h;
    end;
  end;
end;

procedure Produkt(A,B : matrix; var AB: matrix);
{
  AB = A.B
}
var
  NDM,i,j,k : integer;
  h : double;
begin
  NDM := Length(A);
  SetLength(AB,NDM,NDM);
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      h := 0;
      for k := 0 TO NDM-1 do
        h := h + A[i,k]*B[k,j];
      AB[i,j] := h;
    end;
  end;
end;

procedure Kopie(B : matrix; var P : matrix);
{
  P = B
}
var
  NDM,i,j : integer;
begin
  NDM := Length(B);
  SetLength(P,NDM,NDM);
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      P[i,j] := B[i,j]
    end;
  end;
end;

procedure Drukaf(Q : matrix);
{
  Print
}
var
  NDM,i,j : integer;
begin
  NDM := Length(Q);
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      Write(Q[i,j]:12:6);
    end;
    Writeln;
  end;
  Writeln;
end;

procedure Eenheid(NDM : integer; var E : matrix);
{
  Unit matrix
}
var
  i,j : integer;
begin
  SetLength(E,NDM,NDM);
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      E[i,j] := 0;
    end;
    E[i,i] := 1;
  end;
end;

function norm(A : matrix) : double;
var
  i,j,NDM : integer;
  som : double;
begin
  NDM := Length(A);
  som := 0;
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      som := som + sqr(A[i,j]);
    end;
  end;
  norm := sqrt(som);
end;

procedure Newton(P : matrix; var Q : matrix);
{
  Newton-Raphson method computing
  the square root Q of a matrix P
}
const
  veel : integer = 25;
  eps : double = 1.E-10;
var
  keer,NDM,i,j,maal : integer;
  B,PB,R,QQ : matrix;
  fout : double;
begin
  NDM := Length(P);
  Eenheid(NDM,Q);
  Drukaf(Q);
  SetLength(R,NDM,NDM);
  maal := 0;
  for keer := 0 to veel do
  begin
    Kopie(Q,B);
    Inverse(B);
    Produkt(P,B,PB);
    Produkt(Q,Q,QQ);
    for i := 0 to NDM-1 do
    begin
      for j := 0 to NDM-1 do
      begin
        R[i,j] := QQ[i,j]-P[i,j];
        Q[i,j] := (Q[i,j] + PB[i,j])/2;
      end;
    end;
    fout := norm(R);
    Writeln('    Error =',fout);
    Writeln;
    Drukaf(Q);
    maal := keer;
    if fout < eps then Break;
  end;
  Writeln('    Iterations = ',maal);
  Writeln; 
end;

procedure Algemeen(B : matrix);
{
  Testing in general
}
var
  Q,R,BB,U,E : matrix;
begin
  Kwadraat(B,BB);
  Writeln('Newton Raphson -> Q');
  Writeln;
  Newton(BB,Q);
  Produkt(Q,Q,R);
  Writeln('Test Q^2 = B^T.B');
  Writeln;
  Drukaf(R);
  Drukaf(BB);
  Writeln('Orthogonal matrix U + test');
  Writeln;
  Kopie(Q,R);
  Inverse(R);
  Produkt(B,R,U);
  Drukaf(U);
  Kwadraat(U,E);
  Drukaf(E);
  Writeln('Polar Decomposition + test');
  Writeln;
  Drukaf(B);
  Drukaf(U);
  Drukaf(Q);
  Produkt(U,Q,R);
  Drukaf(R);
end;

END.

Main program

program Polair;
{
https://math.stackexchange.com/questions/1456320/polar-decomposition-of-real-matrices
}
Uses Wiskunde;

procedure Test1;
{
  Test random
  3x3 matrix
}
const
  NDM : integer = 3;
var
  B : matrix;
  i,j : integer;
begin
  SetLength(B,NDM,NDM);
  Random; Random;
  for i := 0 to NDM-1 do
  begin
    for j := 0 to NDM-1 do
    begin
      B[i,j] := Random-0.5;
    end;
  end;
  Writeln('Random matrix B');
  Writeln;
  Drukaf(B);
  Algemeen(B);
end;

procedure Test2;
{
https://math.stackexchange.com/questions/252940/polar-decomposition
}
const
  NDM : integer = 2;
var
  B : matrix;
begin
  SetLength(B,NDM,NDM);
  B[0,0] := 1;
  B[0,1] := -2;
  B[1,0] := 3;
  B[1,1] := -1;
  Writeln('Sample matrix');
  Writeln;
  Drukaf(B);
  Algemeen(B);
end;

procedure Test3;
{
https://math.stackexchange.com/questions/706262/find-the-polar-decomposition
}
const
  NDM : integer = 2;
var
  B,Q,U : matrix;
begin
  SetLength(B,NDM,NDM);
  SetLength(Q,NDM,NDM);
  SetLength(U,NDM,NDM);
  B[0,0] := 1;
  B[0,1] := 0;
  B[1,0] := 1;
  B[1,1] := 1;
  Writeln('Sample matrix');
  Writeln;
  Drukaf(B);
  Algemeen(B);
{
https://socratic.org/questions/what-is-1-over-square-root-of-5
}
  U[0,0] := 2/sqrt(5);
  U[0,1] := -1/sqrt(5);
  U[1,0] := 1/sqrt(5);
  U[1,1] := 2/sqrt(5);
  Drukaf(U);
  Q[0,0] := 3/sqrt(5);
  Q[0,1] := 1/sqrt(5);
  Q[1,0] := 1/sqrt(5);
  Q[1,1] := 2/sqrt(5);
  Drukaf(Q);
end;

procedure Test4;
{
https://math.stackexchange.com/questions/1534751/calculating-polar-decomposition
}
const
  NDM : integer = 2;
var
  B,U,Q : matrix;
begin
  SetLength(B,NDM,NDM);
  SetLength(U,NDM,NDM);
  SetLength(Q,NDM,NDM);
  B[0,0] := 2;
  B[0,1] := -3;
  B[1,0] := 1;
  B[1,1] := 6;
  Writeln('Sample matrix');
  Writeln;
  Drukaf(B);
  Algemeen(B);
  U[0,0] := 2/sqrt(5);
  U[0,1] := -1/sqrt(5);
  U[1,0] := 1/sqrt(5);
  U[1,1] := 2/sqrt(5);
  Drukaf(U);
  Q[0,0] := sqrt(5);
  Q[0,1] := 0;
  Q[1,0] := 0;
  Q[1,1] := 3*sqrt(5);
  Drukaf(Q);
end;

begin
  Test1;
end.

Output

Random matrix B

    0.361048   -0.297419   -0.227079
    0.171654   -0.181309   -0.338205
   -0.127762   -0.074326   -0.417988

Newton Raphson -> Q

    1.000000    0.000000    0.000000
    0.000000    1.000000    0.000000
    0.000000    0.000000    1.000000

    Error = 1.40545314281996E+0000

    0.588072   -0.064505   -0.043319
   -0.064505    0.563428    0.079962
   -0.043319    0.079962    0.670330

    Error = 3.13954995174322E-0001

    0.430353   -0.124777   -0.066436
   -0.124777    0.369094    0.142866
   -0.066436    0.142866    0.574098

    Error = 6.67068255114933E-0002

    0.384513   -0.168507   -0.062670
   -0.168507    0.280613    0.174530
   -0.062670    0.174530    0.554677

    Error = 1.53751088229495E-0002

    0.369956   -0.192516   -0.055290
   -0.192516    0.238488    0.188300
   -0.055290    0.188300    0.549922

    Error = 3.64905916733788E-0003

    0.363640   -0.203485   -0.051739
   -0.203485    0.219437    0.194468
   -0.051739    0.194468    0.547924

    Error = 7.48757484221924E-0004

    0.361442   -0.207302   -0.050503
   -0.207302    0.212807    0.196615
   -0.050503    0.196615    0.547229

    Error = 9.06824086032184E-0005

    0.361091   -0.207912   -0.050306
   -0.207912    0.211748    0.196958
   -0.050306    0.196958    0.547118

    Error = 2.31633678002109E-0006

    0.361082   -0.207928   -0.050301
   -0.207928    0.211719    0.196967
   -0.050301    0.196967    0.547115

    Error = 1.67844089318140E-0009

    0.361082   -0.207928   -0.050301
   -0.207928    0.211719    0.196967
   -0.050301    0.196967    0.547115

    Error = 4.98075748207982E-0012

    0.361082   -0.207928   -0.050301
   -0.207928    0.211719    0.196967
   -0.050301    0.196967    0.547115

    Iterations = 9

Test Q^2 = B^T.B

    0.176144   -0.129009   -0.086638
   -0.129009    0.126855    0.159924
   -0.086638    0.159924    0.340661

    0.176144   -0.129009   -0.086638
   -0.129009    0.126855    0.159924
   -0.086638    0.159924    0.340661

Orthogonal matrix U + test

    0.491571   -0.868969   -0.057015
    0.599413    0.385125   -0.701701
   -0.631714   -0.310760   -0.710187

    1.000000    0.000000   -0.000000
    0.000000    1.000000   -0.000000
   -0.000000   -0.000000    1.000000

Polar Decomposition + test

    0.361048   -0.297419   -0.227079
    0.171654   -0.181309   -0.338205
   -0.127762   -0.074326   -0.417988

    0.491571   -0.868969   -0.057015
    0.599413    0.385125   -0.701701
   -0.631714   -0.310760   -0.710187

    0.361082   -0.207928   -0.050301
   -0.207928    0.211719    0.196967
   -0.050301    0.196967    0.547115

    0.361048   -0.297419   -0.227079
    0.171654   -0.181309   -0.338205
   -0.127762   -0.074326   -0.417988