Conjugate transpose: Difference between revisions
(→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
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.
- MathWorld: conjugate transpose, Hermitian matrix, normal matrix, 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
<lang ruby>require 'matrix'
- Start with some matrix.
i = Complex::I matrix = Matrix[[i, 0, 0],
[0, i, 0], [0, 0, i]]
- Find the conjugate transpose.
- 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
<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>