Conjugate transpose

From Rosetta Code
Revision as of 22:28, 27 March 2012 by rosettacode>Kernigh (→‎{{header|Factor}}: Link to Factor bug #484.)
Conjugate transpose is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Suppose that a matrix contains complex numbers. Then the conjugate transpose of is a matrix containing the complex conjugates of the matrix transposition of .

This means that row , column of the conjugate transpose equals the complex conjugate of row , column of the original matrix.

In the next list, must also be a square matrix.

  • A Hermitian matrix equals its own conjugate transpose: .
  • A normal matrix is commutative in multiplication with its conjugate transpose: .
  • A unitary matrix has its inverse equal to its conjugate transpose: . This is true iff and iff , where is the identity matrix.

Given some matrix of complex numbers, find its conjugate transpose. Also determine if it is a Hermitian matrix, normal matrix, or a unitary matrix.

Ada

<lang Ada>with Ada.Text_IO; use Ada.Text_IO; with Ada.Complex_Text_IO; use Ada.Complex_Text_IO; with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types; with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays; procedure ConTrans is

  subtype CM is Complex_Matrix;
  S2O2 : constant Float := 0.7071067811865;
  procedure Print (mat : CM) is begin
     for row in mat'Range(1) loop for col in mat'Range(2) loop
        Put(mat(row,col), Exp=>0, Aft=>4);
     end loop; New_Line; end loop;
  end Print;
  function almostzero(mat : CM; tol : Float) return Boolean is begin
     for row in mat'Range(1) loop for col in mat'Range(2) loop
        if abs(mat(row,col)) > tol then return False; end if;
     end loop; end loop;
     return True;
  end almostzero;
  procedure Examine (mat : CM) is
     CT : CM := Conjugate (Transpose(mat));
     isherm, isnorm, isunit : Boolean;
  begin
     isherm := almostzero(mat-CT, 1.0e-6);
     isnorm := almostzero(mat*CT-CT*mat, 1.0e-6);
     isunit := almostzero(CT-Inverse(mat), 1.0e-6);
     Print(mat);
     Put_Line("Conjugate transpose:"); Print(CT);
     Put_Line("Hermitian?: " & isherm'Img);
     Put_Line("Normal?: " & isnorm'Img);
     Put_Line("Unitary?: " & isunit'Img);
  end Examine;
  hmat : CM := ((3.0+0.0*i, 2.0+1.0*i), (2.0-1.0*i, 1.0+0.0*i));
  nmat : CM := ((1.0+0.0*i, 1.0+0.0*i, 0.0+0.0*i),
                (0.0+0.0*i, 1.0+0.0*i, 1.0+0.0*i),
                (1.0+0.0*i, 0.0+0.0*i, 1.0+0.0*i));
  umat : CM := ((S2O2+0.0*i, S2O2+0.0*i, 0.0+0.0*i),
                (0.0+S2O2*i, 0.0-S2O2*i, 0.0+0.0*i),
                (0.0+0.0*i, 0.0+0.0*i, 0.0+1.0*i));

begin

  Put_Line("hmat:"); Examine(hmat); New_Line;
  Put_Line("nmat:"); Examine(nmat); New_Line;
  Put_Line("umat:"); Examine(umat);

end ConTrans;</lang>

Output:
hmat:
( 3.0000, 0.0000)( 2.0000, 1.0000)
( 2.0000,-1.0000)( 1.0000, 0.0000)
Conjugate transpose:
( 3.0000,-0.0000)( 2.0000, 1.0000)
( 2.0000,-1.0000)( 1.0000,-0.0000)
Hermitian?: TRUE
Normal?: TRUE
Unitary?: FALSE

nmat:
( 1.0000, 0.0000)( 1.0000, 0.0000)( 0.0000, 0.0000)
( 0.0000, 0.0000)( 1.0000, 0.0000)( 1.0000, 0.0000)
( 1.0000, 0.0000)( 0.0000, 0.0000)( 1.0000, 0.0000)
Conjugate transpose:
( 1.0000,-0.0000)( 0.0000,-0.0000)( 1.0000,-0.0000)
( 1.0000,-0.0000)( 1.0000,-0.0000)( 0.0000,-0.0000)
( 0.0000,-0.0000)( 1.0000,-0.0000)( 1.0000,-0.0000)
Hermitian?: FALSE
Normal?: TRUE
Unitary?: FALSE

umat:
( 0.7071, 0.0000)( 0.7071, 0.0000)( 0.0000, 0.0000)
( 0.0000, 0.7071)( 0.0000,-0.7071)( 0.0000, 0.0000)
( 0.0000, 0.0000)( 0.0000, 0.0000)( 0.0000, 1.0000)
Conjugate transpose:
( 0.7071,-0.0000)( 0.0000,-0.7071)( 0.0000,-0.0000)
( 0.7071,-0.0000)( 0.0000, 0.7071)( 0.0000,-0.0000)
( 0.0000,-0.0000)( 0.0000,-0.0000)( 0.0000,-1.0000)
Hermitian?: FALSE
Normal?: TRUE
Unitary?: TRUE

Factor

This example is incorrect. Please fix the code and remove this message.

Details: By using m. for matrix multiplication, this code encounters bug #484 in Factor: see Talk:Conjugate transpose#Bug in Factor.

<lang factor>USING: kernel math.functions math.matrices sequences ; IN: rosetta.hermitian

conj-t ( matrix -- conjugate-transpose )
   flip [ [ conjugate ] map ] map ;
hermitian-matrix? ( matrix -- ? )
   dup conj-t = ;
normal-matrix? ( matrix -- ? )
   dup conj-t [ m. ] [ swap m. ] 2bi = ;
unitary-matrix? ( matrix -- ? )
   [ dup conj-t m. ] [ length identity-matrix ] bi = ;</lang>

Go

<lang go>package main

import (

   "fmt"
   "math"
   "math/cmplx"

)

// a type to represent matrices type matrix struct {

   ele  []complex128
   cols int

}

// conjugate transpose, implemented here as a method on the matrix type. func (m *matrix) conjTranspose() *matrix {

   r := &matrix{make([]complex128, len(m.ele)), len(m.ele) / m.cols}
   rx := 0
   for _, e := range m.ele {
       r.ele[rx] = cmplx.Conj(e)
       rx += r.cols
       if rx >= len(r.ele) {
           rx -= len(r.ele) - 1
       }
   }
   return r

}

// program to demonstrate capabilites on example matricies func main() {

   show("h", matrixFromRows([][]complex128{
       {3, 2 + 1i},
       {2 - 1i, 1}}))
   show("n", matrixFromRows([][]complex128{
       {1, 1, 0},
       {0, 1, 1},
       {1, 0, 1}}))
   show("u", matrixFromRows([][]complex128{
       {math.Sqrt2 / 2, math.Sqrt2 / 2, 0},
       {math.Sqrt2 / -2i, math.Sqrt2 / 2i, 0},
       {0, 0, 1i}}))

}

func show(name string, m *matrix) {

   m.print(name)
   ct := m.conjTranspose()
   ct.print(name + "_ct")
   fmt.Println("Hermitian:", m.equal(ct, 1e-14))
   mct := m.mult(ct)
   ctm := ct.mult(m)
   fmt.Println("Normal:", mct.equal(ctm, 1e-14))
   i := eye(m.cols)
   fmt.Println("Unitary:", mct.equal(i, 1e-14) && ctm.equal(i, 1e-14))

}

// two constructors func matrixFromRows(rows [][]complex128) *matrix {

   m := &matrix{make([]complex128, len(rows)*len(rows[0])), len(rows[0])}
   for rx, row := range rows {
       copy(m.ele[rx*m.cols:(rx+1)*m.cols], row)
   }
   return m

}

func eye(n int) *matrix {

   r := &matrix{make([]complex128, n*n), n}
   n++
   for x := 0; x < len(r.ele); x += n {
       r.ele[x] = 1
   }
   return r

}

// print method outputs matrix to stdout func (m *matrix) print(heading string) {

   fmt.Print("\n", heading, "\n")
   for e := 0; e < len(m.ele); e += m.cols {
       fmt.Printf("%6.3f ", m.ele[e:e+m.cols])
       fmt.Println()
   }

}

// equal method uses ε to allow for floating point error. func (a *matrix) equal(b *matrix, ε float64) bool {

   for x, aEle := range a.ele {
       if math.Abs(real(aEle)-real(b.ele[x])) > math.Abs(real(aEle))*ε ||
           math.Abs(imag(aEle)-imag(b.ele[x])) > math.Abs(imag(aEle))*ε {
           return false
       }
   }
   return true

}

// mult method taken from matrix multiply task func (m1 *matrix) mult(m2 *matrix) (m3 *matrix) {

   m3 = &matrix{make([]complex128, (len(m1.ele)/m1.cols)*m2.cols), m2.cols}
   for m1c0, m3x := 0, 0; m1c0 < len(m1.ele); m1c0 += m1.cols {
       for m2r0 := 0; m2r0 < m2.cols; m2r0++ {
           for m1x, m2x := m1c0, m2r0; m2x < len(m2.ele); m2x += m2.cols {
               m3.ele[m3x] += m1.ele[m1x] * m2.ele[m2x]
               m1x++
           }
           m3x++
       }
   }
   return m3

}</lang> Output:

h
[( 3.000+0.000i) (+2.000+1.000i)] 
[( 2.000-1.000i) (+1.000+0.000i)] 

h_ct
[( 3.000-0.000i) (+2.000+1.000i)] 
[( 2.000-1.000i) (+1.000-0.000i)] 
Hermitian: true
Normal: true
Unitary: false

n
[( 1.000+0.000i) (+1.000+0.000i) (+0.000+0.000i)] 
[( 0.000+0.000i) (+1.000+0.000i) (+1.000+0.000i)] 
[( 1.000+0.000i) (+0.000+0.000i) (+1.000+0.000i)] 

n_ct
[( 1.000-0.000i) (+0.000-0.000i) (+1.000-0.000i)] 
[( 1.000-0.000i) (+1.000-0.000i) (+0.000-0.000i)] 
[( 0.000-0.000i) (+1.000-0.000i) (+1.000-0.000i)] 
Hermitian: false
Normal: true
Unitary: false

u
[( 0.707+0.000i) (+0.707+0.000i) (+0.000+0.000i)] 
[( 0.000+0.707i) (+0.000-0.707i) (+0.000+0.000i)] 
[( 0.000+0.000i) (+0.000+0.000i) (+0.000+1.000i)] 

u_ct
[( 0.707-0.000i) (+0.000-0.707i) (+0.000-0.000i)] 
[( 0.707-0.000i) (+0.000+0.707i) (+0.000-0.000i)] 
[( 0.000-0.000i) (+0.000-0.000i) (+0.000-1.000i)] 
Hermitian: false
Normal: true
Unitary: true

J

Solution: <lang j> ct =: +@|: NB. Conjugate transpose (ct A is A_ct)</lang> Examples: <lang j> X =: +/ . * NB. Matrix Multiply (x)

  HERMITIAN =:  3 2j1 ,: 2j_1 1  
  (-: ct) HERMITIAN               NB.  A_ct = A

1

  NORMAL    =:  1 1 0 , 0 1 1 ,: 1 0 1
  ((X~ -: X) ct) NORMAL           NB. A_ct x A = A x A_ct

1

  UNITARY   =:  (-:%:2) * 1 1 0 , 0j_1 0j1 0 ,: 0 0 0j1 * %:2
  (ct -: %.)  UNITARY             NB.  A_ct = A^-1

1</lang>

Reference (example matrices for other langs to use):<lang j> HERMITIAN;NORMAL;UNITARY +--------+-----+--------------------------+ | 3 2j1|1 1 0| 0.707107 0.707107 0| |2j_1 1|0 1 1|0j_0.707107 0j0.707107 0| | |1 0 1| 0 0 0j1| +--------+-----+--------------------------+

  NB. In J, PjQ is P + Q*i and the 0.7071... is sqrt(2)</lang>

Mathematica

<lang Mathematica>NormalMatrixQ[a_List?MatrixQ] := Module[{b = Conjugate@Transpose@a},a.b === b.a] UnitaryQ[m_List?MatrixQ] := (Conjugate@Transpose@m.m == IdentityMatrix@Length@m)

m = {{1, 2I, 3}, {3+4I, 5, I}}; m //MatrixForm -> (1 2I 3 3+4I 5 I)

ConjugateTranspose[m] //MatrixForm -> (1 3-4I -2I 5 3 -I)

{HermitianMatrixQ@#, NormalMatrixQ@#, UnitaryQ@#}&@m -> {False, False, False}</lang>


Ruby

This example is incorrect. Please fix the code and remove this message.

Details: Matrix#hermitian? in MRI uses a different definition of Hermitian matrix: it only checks for .

Works with: Ruby version 1.9.3

<lang ruby>require 'matrix'

  1. Start with some matrix.

i = Complex::I matrix = Matrix[[i, 0, 0],

               [0, i, 0],
               [0, 0, i]]
  1. Find the conjugate transpose.
  2. Matrix#conjugate appeared in Ruby 1.9.2.

conjt = matrix.conj.t # aliases for matrix.conjugate.tranpose print 'conjugate tranpose: '; puts conjt

if matrix.square?

 # These predicates appeared in Ruby 1.9.3.
 print 'Hermitian? '; puts matrix.hermitian?
 print '   normal? '; puts matrix.normal?
 print '  unitary? '; puts matrix.unitary?

else

 # Matrix is not square. These predicates would
 # raise ExceptionForMatrix::ErrDimensionMismatch.
 print 'Hermitian? false'
 print '   normal? false'
 print '  unitary? false'

end</lang>

Tcl

Tcl's matrixes (in Tcllib) do not assume that the contents are numeric at all. As such, they do not provide mathematical operations over them and this considerably increases the complexity of the code below. Note the use of lambda terms to simplify access to the complex number package.

Library: Tcllib (Package: math::complexnumbers)
Library: Tcllib (Package: struct::matrix)

<lang tcl>package require struct::matrix package require math::complexnumbers

proc complexMatrix.equal {m1 m2 {epsilon 1e-14}} {

   if {[$m1 rows] != [$m2 rows] || [$m1 columns] != [$m2 columns]} {

return 0

   }
   # Compute the magnitude of the difference between two complex numbers
   set ceq [list apply {{epsilon a b} {

expr {[mod [- $a $b]] < $epsilon}

   } ::math::complexnumbers} $epsilon]
   for {set i 0} {$i<[$m1 columns]} {incr i} {

for {set j 0} {$j<[$m1 rows]} {incr j} { if {![{*}$ceq [$m1 get cell $i $j] [$m2 get cell $i $j]]} { return 0 } }

   }
   return 1

}

proc complexMatrix.multiply {a b} {

   if {[$a columns] != [$b rows]} {
       error "incompatible sizes"
   }
   # Simplest to use a lambda in the complex NS
   set cpm {{sum a b} {

+ $sum [* $a $b]

   } ::math::complexnumbers}
   set c0 [math::complexnumbers::complex 0.0 0.0];   # Complex zero
   set c [struct::matrix]
   $c add columns [$b columns]
   $c add rows [$a rows]
   for {set i 0} {$i < [$a rows]} {incr i} {
       for {set j 0} {$j < [$b columns]} {incr j} {
           set sum $c0

foreach rv [$a get row $i] cv [$b get column $j] { set sum [apply $cpm $sum $rv $cv]

           }

$c set cell $j $i $sum

       }
   }
   return $c

}

proc complexMatrix.conjugateTranspose {matrix} {

   set mat [struct::matrix]
   $mat = $matrix
   $mat transpose
   for {set c 0} {$c < [$mat columns]} {incr c} {

for {set r 0} {$r < [$mat rows]} {incr r} { set val [$mat get cell $c $r] $mat set cell $c $r [math::complexnumbers::conj $val] }

   }
   return $mat

}</lang> Using these tools to test for the properties described in the task: <lang tcl>proc isHermitian {matrix {epsilon 1e-14}} {

   if {[$matrix rows] != [$matrix columns]} {

# Must be square! return 0

   }
   set cc [complexMatrix.conjugateTranspose $matrix]
   set result [complexMatrix.equal $matrix $cc $epsilon]
   $cc destroy
   return $result

}

proc isNormal {matrix {epsilon 1e-14}} {

   if {[$matrix rows] != [$matrix columns]} {

# Must be square! return 0

   }
   set mh [complexMatrix.conjugateTranspose $matrix]
   set mhm [complexMatrix.multiply $mh $matrix]
   set mmh [complexMatrix.multiply $matrix $mh]
   $mh destroy
   set result [complexMatrix.equal $mhm $mmh $epsilon]
   $mhm destroy
   $mmh destroy
   return $result

}

proc isUnitary {matrix {epsilon 1e-14}} {

   if {[$matrix rows] != [$matrix columns]} {

# Must be square! return 0

   }
   set mh [complexMatrix.conjugateTranspose $matrix]
   set mhm [complexMatrix.multiply $mh $matrix]
   set mmh [complexMatrix.multiply $matrix $mh]
   $mh destroy
   set result [complexMatrix.equal $mhm $mmh $epsilon]
   $mhm destroy
   if {$result} {

set id [struct::matrix] $id = $matrix; # Just for its dimensions for {set c 0} {$c < [$id columns]} {incr c} { for {set r 0} {$r < [$id rows]} {incr r} { $id set cell $c $r \ [math::complexnumbers::complex [expr {$c==$r}] 0] } } set result [complexMatrix.equal $mmh $id $epsilon] $id destroy

   }
   $mmh destroy
   return $result

}</lang>