Conjugate transpose: Difference between revisions

=={{header|jq}}==

{{works with|jq|1.4}}

In the following, we use the array [x,y] to represent the complex number x + iy,

but the following functions also accept a number wherever a complex number is acceptable.

====Infrastructure====

'''(1) transpose/0''':

If your jq does not have "transpose" then the following may be used:

<lang jq># transpose/0 expects its input to be a rectangular matrix

# (an array of equal-length arrays):

def transpose:

if (.[0] | length) == 0 then []

else [map(.[0])] + (map(.[1:]) | transpose)

end ;</lang>

'''(2) Operations on real/complex numbers'''

<lang jq># x must be real or complex, and ditto for y;

# always return complex

def plus(x; y):

if (x|type) == "number" then

if (y|type) == "number" then [ x+y, 0 ]

else [ x + y[0], y[1]]

end

elif (y|type) == "number" then plus(y;x)

else [ x[0] + y[0], x[1] + y[1] ]

end;

# x must be real or complex, and ditto for y;

# always return complex

def multiply(x; y):

if (x|type) == "number" then

if (y|type) == "number" then [ x*y, 0 ]

else [x * y[0], x * y[1]]

end

elif (y|type) == "number" then multiply(y;x)

else [ x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]]

end;

# conjugate of a real or complex number

def conjugate:

if type == "number" then [.,0]

else [.[0], -(.[1]) ]

end;</lang>

'''(3) Array operations'''

<lang jq># a and b are arrays of real/complex numbers

def dot_product(a; b):

a as $a | b as $b

| reduce range(0;$a|length) as $i

(0; . as $s | plus($s; multiply($a[$i]; $b[$i]) ));</lang>

'''(4) Matrix operations'''

<lang jq># convert a matrix of mixed real/complex entries to all complex entries

def to_complex:

def toc: if type == "number" then [.,0] else . end;

map( map(toc) );

# simple matrix pretty-printer

def pp(wide):

def pad: tostring | (wide - length) * " " + .;

def row: reduce .[] as $x (""; . + ($x|pad));

reduce .[] as $row (""; . + "\n\($row|row)");

# Matrix multiplication

# A and B should both be real/complex matrices,

# A being m by n, and B being n by p.

def matrix_multiply(A; B):

A as $A | B as $B

| ($B[0]|length) as $p

| ($B|transpose) as $BT

| reduce range(0; $A|length) as $i

([]; reduce range(0; $p) as $j

(.; .[$i][$j] = dot_product( $A[$i]; $BT[$j] ) )) ;

# Complex identity matrix of dimension n

def complex_identity(n):

def indicator(i;n): [range(0;n)] | map( [0,0]) | .[i] = [1,0];

reduce range(0; n) as $i ([]; . + [indicator( $i; n )] );

# Approximate equality of two matrices

# Are two real/complex matrices essentially equal

# in the sense that the sum of the squared element-wise differences

# is less than or equal to epsilon?

# The two matrices must be conformal.

def approximately_equal(M; N; epsilon):

def norm: multiply(. ; conjugate ) | .[0];

def sqdiff( x; y): plus(x; multiply(y; -1)) | norm;

reduce range(0;M|length) as $i

(0; reduce range(0; M[0]|length) as $j

(.; 0 + sqdiff( M[$i][$j]; N[$i][$j] ) ) ) <= epsilon;</lang>

====Conjugate transposition====

<lang jq># (entries may be real and/or complex)

def conjugate_transpose:

map( map(conjugate) ) | transpose;

# A Hermitian matrix equals its own conjugate transpose

def is_hermitian:

to_complex == conjugate_transpose;

# A matrix is normal if it commutes multiplicatively

# with its conjugate transpose

def is_normal:

. as $M

| conjugate_transpose as $H

| matrix_multiply($H; $M) == matrix_multiply($H; $M);

# A unitary matrix (U) has its inverse equal to its conjugate transpose (T)

# i.e. U^-1 == T; NASC is I == UT == TU

def is_unitary:

. as $M

| conjugate_transpose as $H

| complex_identity(length) as $I

| approximately_equal( $I; matrix_multiply($H;$M); 1e-10)

and approximately_equal( $I ; matrix_multiply($M;$H); 1e-10) ; </lang>

====Examples====

<lang jq>def hermitian_example:

[ [ 3, [2,1]],

[[2,-1], 1 ] ];

def normal_example:

[ [1, 1, 0],

[0, 1, 1],

[1, 0, 1] ];

def unitary_example:

0.707107

| [ [ [., 0], [., 0], 0 ],

[ [0, -.], [0, .], 0 ],

[ 0, 0, [0,1] ] ];

def demo:

hermitian_example

| ("Hermitian example:", pp(8)),

"",

("Its conjugate transpose is:", (to_complex | conjugate_transpose | pp(8))),

"",

"Hermitian example: \(hermitian_example | is_hermitian )",

"",

"Normal example: \(normal_example | is_normal )",

"",

"Unitary example: \(unitary_example | is_unitary)"

;

demo</lang>

{{out}}

<lang sh>$ jq -r -c -n -f Conjugate_transpose.jq

Hermitian example:

3 [2,1]

[2,-1] 1

Conjugate transpose:

[3,-0] [2,1]

[2,-1] [1,-0]

Hermitian example: true

Normal example: true

Unitary example: true</lang>