Conjugate transpose: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎Tcl: Added implementation)
Line 241: Line 241:
print ' unitary? false'
print ' unitary? false'
end</lang>
end</lang>

=={{header|Tcl}}==
{{tcllib|math::complexnumbers}}
{{tcllib|struct::matrix}}
<lang tcl>package require struct::matrix
package require math::complexnumbers

proc 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
}

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 isHermitian {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
return 0
}
set cc [conjugateTranspose $matrix]
set result [complexMatrix.equal $matrix $cc $epsilon]
$cc destroy
return $result
}

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 isNormal {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
return 0
}
set mh [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 [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>
<!-- Wot, no demonstration? -->

Revision as of 11:16, 14 March 2012

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.

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

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

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

proc 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

}

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 isHermitian {matrix {epsilon 1e-14}} {

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

# Must be square! return 0

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

}

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 isNormal {matrix {epsilon 1e-14}} {

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

# Must be square! return 0

   }
   set mh [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 [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>