Conjugate transpose

From Rosetta Code
Task
Conjugate transpose
You are encouraged to solve this task according to the task description, using any language you may know.

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.



Task

Given some matrix of complex numbers, find its conjugate transpose.

Also determine if the matrix is a:

  • Hermitian matrix,
  • normal matrix, or
  • unitary matrix.


See also



Ada[edit]

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;
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

C[edit]

/* Uses C99 specified complex.h, complex datatype has to be defined and operation provided if used on non-C99 compilers */
 
#include<stdlib.h>
#include<stdio.h>
#include<complex.h>
 
typedef struct
{
int rows, cols;
complex **z;
} matrix;
 
matrix
transpose (matrix a)
{
int i, j;
matrix b;
 
b.rows = a.cols;
b.cols = a.rows;
 
b.z = malloc (b.rows * sizeof (complex *));
 
for (i = 0; i < b.rows; i++)
{
b.z[i] = malloc (b.cols * sizeof (complex));
for (j = 0; j < b.cols; j++)
{
b.z[i][j] = conj (a.z[j][i]);
}
}
 
return b;
}
 
int
isHermitian (matrix a)
{
int i, j;
matrix b = transpose (a);
 
if (b.rows == a.rows && b.cols == a.cols)
{
for (i = 0; i < b.rows; i++)
{
for (j = 0; j < b.cols; j++)
{
if (b.z[i][j] != a.z[i][j])
return 0;
}
}
}
 
else
return 0;
 
return 1;
}
 
matrix
multiply (matrix a, matrix b)
{
matrix c;
int i, j;
 
if (a.cols == b.rows)
{
c.rows = a.rows;
c.cols = b.cols;
 
c.z = malloc (c.rows * (sizeof (complex *)));
 
for (i = 0; i < c.rows; i++)
{
c.z[i] = malloc (c.cols * sizeof (complex));
c.z[i][j] = 0 + 0 * I;
for (j = 0; j < b.cols; j++)
{
c.z[i][j] += a.z[i][j] * b.z[j][i];
}
}
 
}
 
return c;
}
 
int
isNormal (matrix a)
{
int i, j;
matrix a_ah, ah_a;
 
if (a.rows != a.cols)
return 0;
 
a_ah = multiply (a, transpose (a));
ah_a = multiply (transpose (a), a);
 
for (i = 0; i < a.rows; i++)
{
for (j = 0; j < a.cols; j++)
{
if (a_ah.z[i][j] != ah_a.z[i][j])
return 0;
}
}
 
return 1;
}
 
int
isUnitary (matrix a)
{
matrix b;
int i, j;
if (isNormal (a) == 1)
{
b = multiply (a, transpose(a));
 
for (i = 0; i < b.rows; i++)
{
for (j = 0; j < b.cols; j++)
{
if ((i == j && b.z[i][j] != 1) || (i != j && b.z[i][j] != 0))
return 0;
}
}
return 1;
}
return 0;
}
 
 
int
main ()
{
complex z = 3 + 4 * I;
matrix a, aT;
int i, j;
printf ("Enter rows and columns :");
scanf ("%d%d", &a.rows, &a.cols);
 
a.z = malloc (a.rows * sizeof (complex *));
printf ("Randomly Generated Complex Matrix A is : ");
for (i = 0; i < a.rows; i++)
{
printf ("\n");
a.z[i] = malloc (a.cols * sizeof (complex));
for (j = 0; j < a.cols; j++)
{
a.z[i][j] = rand () % 10 + rand () % 10 * I;
printf ("\t%f + %fi", creal (a.z[i][j]), cimag (a.z[i][j]));
}
}
 
aT = transpose (a);
 
printf ("\n\nTranspose of Complex Matrix A is : ");
for (i = 0; i < aT.rows; i++)
{
printf ("\n");
aT.z[i] = malloc (aT.cols * sizeof (complex));
for (j = 0; j < aT.cols; j++)
{
aT.z[i][j] = rand () % 10 + rand () % 10 * I;
printf ("\t%f + %fi", creal (aT.z[i][j]), cimag (aT.z[i][j]));
}
}
 
printf ("\n\nComplex Matrix A %s hermitian",
isHermitian (a) == 1 ? "is" : "is not");
printf ("\n\nComplex Matrix A %s unitary",
isUnitary (a) == 1 ? "is" : "is not");
printf ("\n\nComplex Matrix A %s normal",
isNormal (a) == 1 ? "is" : "is not");
 
 
 
return 0;
}
Output:
Enter rows and columns :3 3
Randomly Generated Complex Matrix A is :
        3.000000 + 6.000000i    7.000000 + 5.000000i    3.000000 + 5.000000i
        6.000000 + 2.000000i    9.000000 + 1.000000i    2.000000 + 7.000000i
        0.000000 + 9.000000i    3.000000 + 6.000000i    0.000000 + 6.000000i

Transpose of Complex Matrix A is :
        2.000000 + 6.000000i    1.000000 + 8.000000i    7.000000 + 9.000000i
        2.000000 + 0.000000i    2.000000 + 3.000000i    7.000000 + 5.000000i
        9.000000 + 2.000000i    2.000000 + 8.000000i    9.000000 + 7.000000i

Complex Matrix A is not hermitian

Complex Matrix A is not unitary

Complex Matrix A is not normal

Common Lisp[edit]

 
(defun matrix-multiply (m1 m2)
(mapcar
(lambda (row)
(apply #'mapcar
(lambda (&rest column)
(apply #'+ (mapcar #'* row column))) m2)) m1))
 
(defun identity-p (m &optional (tolerance 1e-6))
"Is m an identity matrix?"
(loop for row in m
for r = 1 then (1+ r) do
(loop for col in row
for c = 1 then (1+ c) do
(if (eql r c)
(unless (< (abs (- col 1)) tolerance) (return-from identity-p nil))
(unless (< (abs col) tolerance) (return-from identity-p nil)) )))
T )
 
(defun conjugate-transpose (m)
(apply #'mapcar #'list (mapcar #'(lambda (r) (mapcar #'conjugate r)) m)) )
 
(defun hermitian-p (m)
(equalp m (conjugate-transpose m)))
 
(defun normal-p (m)
(let ((m* (conjugate-transpose m)))
(equalp (matrix-multiply m m*) (matrix-multiply m* m)) ))
 
(defun unitary-p (m)
(identity-p (matrix-multiply m (conjugate-transpose m))) )
 
Output:
(hermitian-p
  '((3        #C(2 1))
    (#C(2 -1) 1) ))
=> T

(normal-p
  '((#C(0 1) 0)
    (0       #C(3 -5)) ))
==> T

(unitary-p
  '((0.70710677        0.70710677       0)
    (#C(0 -0.70710677) #C(0 0.70710677) 0)
    (0                 0                1) ))
==> T

D[edit]

Translation of: Python
A well typed and mostly imperative version:
import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric;
 
T[][] conjugateTranspose(T)(in T[][] m) pure nothrow @safe {
auto r = new typeof(return)(m[0].length, m.length);
foreach (immutable nr, const row; m)
foreach (immutable nc, immutable c; row)
r[nc][nr] = c.conj;
return r;
}
 
bool isRectangular(T)(in T[][] M) pure nothrow @safe @nogc {
return M.all!(row => row.length == M[0].length);
}
 
T[][] matMul(T)(in T[][] A, in T[][] B) pure nothrow /*@safe*/
in {
assert(A.isRectangular && B.isRectangular &&
!A.empty && !B.empty && A[0].length == B.length);
} body {
auto result = new T[][](A.length, B[0].length);
auto aux = new T[B.length];
 
foreach (immutable j; 0 .. B[0].length) {
foreach (immutable k, const row; B)
aux[k] = row[j];
foreach (immutable i, const ai; A)
result[i][j] = dotProduct(ai, aux);
}
 
return result;
}
 
/// Check any number of complex matrices for equality within
/// some bits of mantissa.
bool areEqual(T)(in Complex!T[][][] matrices, in size_t nBits=20)
pure nothrow /*@safe*/ {
static bool allSame(U)(in U[] v) pure nothrow @nogc {
return v[1 .. $].all!(c => c == v[0]);
}
 
bool allNearSame(in Complex!T[] v) pure nothrow @nogc {
auto v0 = v[0].Complex!T; // To avoid another cast.
return v[1 .. $].all!(c => feqrel(v0.re, c.re) >= nBits &&
feqrel(v0.im, c.im) >= nBits);
}
 
immutable x = matrices.map!(m => m.length).array;
if (!allSame(x))
return false;
immutable y = matrices.map!(m => m[0].length).array;
if (!allSame(y))
return false;
foreach (immutable s; 0 .. x[0])
foreach (immutable t; 0 .. y[0])
if (!allNearSame(matrices.map!(m => m[s][t]).array))
return false;
return true;
}
 
bool isHermitian(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow /*@safe*/ {
return [m, ct].areEqual;
}
 
bool isNormal(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow /*@safe*/ {
return [matMul(m, ct), matMul(ct, m)].areEqual;
}
 
auto complexIdentitymatrix(in size_t side) pure nothrow /*@safe*/ {
return side.iota.map!(r => side.iota.map!(c => complex(r == c)).array).array;
}
 
bool isUnitary(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow /*@safe*/ {
immutable mct = matMul(m, ct);
immutable ident = mct.length.complexIdentitymatrix;
return [mct, matMul(ct, m), ident].areEqual;
}
 
void main() /*@safe*/ {
alias C = complex;
immutable x = 2 ^^ 0.5 / 2;
 
immutable data = [[[C(3.0, 0.0), C(2.0, 1.0)],
[C(2.0, -1.0), C(1.0, 0.0)]],
 
[[C(1.0, 0.0), C(1.0, 0.0), C(0.0, 0.0)],
[C(0.0, 0.0), C(1.0, 0.0), C(1.0, 0.0)],
[C(1.0, 0.0), C(0.0, 0.0), C(1.0, 0.0)]],
 
[[C(x, 0.0), C(x, 0.0), C(0.0, 0.0)],
[C(0.0, -x), C(0.0, x), C(0.0, 0.0)],
[C(0.0, 0.0), C(0.0, 0.0), C(0.0, 1.0)]]];
 
foreach (immutable mat; data) {
enum mFormat = "[%([%(%1.3f, %)],\n %)]]";
writefln("Matrix:\n" ~ mFormat, mat);
immutable ct = conjugateTranspose(mat);
"Its conjugate transpose:".writeln;
writefln(mFormat, ct);
writefln("Hermitian? %s.", isHermitian(mat, ct));
writefln("Normal?  %s.", isNormal(mat, ct));
writefln("Unitary?  %s.\n", isUnitary(mat, ct));
}
}
Output:
Matrix:
[[3.000+0.000i, 2.000+1.000i],
 [2.000-1.000i, 1.000+0.000i]]
Its conjugate transpose:
[[3.000-0.000i, 2.000+1.000i],
 [2.000-1.000i, 1.000-0.000i]]
Hermitian? true.
Normal?    true.
Unitary?   false.

Matrix:
[[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]]
Its conjugate transpose:
[[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.

Matrix:
[[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]]
Its conjugate transpose:
[[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.

Alternative Version[edit]

A more functional version that contains some typing problems (same output).

import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric, std.exception, std.traits;
 
// alias CM(T) = Complex!T[][]; // Not yet useful.
 
auto conjugateTranspose(T)(in Complex!T[][] m) pure nothrow /*@safe*/
if (!hasIndirections!T) {
return iota(m[0].length).map!(i => m.transversal(i).map!conj.array).array;
}
 
T[][] matMul(T)(immutable T[][] A, immutable T[][] B) pure nothrow /*@safe*/ {
immutable Bt = B[0].length.iota.map!(i => B.transversal(i).array).array;
return A.map!(a => Bt.map!(b => a.dotProduct(b)).array).array;
}
 
/// Check any number of complex matrices for equality within
/// some bits of mantissa.
bool areEqual(T)(in Complex!T[][][] matrices, in size_t nBits=20)
pure nothrow /*@safe*/ {
static bool allSame(U)(in U[] v) pure nothrow @nogc @safe {
return v[1 .. $].all!(c => c == v[0]);
}
 
bool allNearSame(in Complex!T[] v) pure nothrow @nogc @safe {
auto v0 = v[0].Complex!T; // To avoid another cast.
return v[1 .. $].all!(c => feqrel(v0.re, c.re) >= nBits &&
feqrel(v0.im, c.im) >= nBits);
}
 
immutable x = matrices.map!(m => m.length).array;
if (!allSame(x))
return false;
immutable y = matrices.map!(m => m[0].length).array;
if (!allSame(y))
return false;
foreach (immutable s; 0 .. x[0])
foreach (immutable t; 0 .. y[0])
if (!allNearSame(matrices.map!(m => m[s][t]).array))
return false;
return true;
}
 
bool isHermitian(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow /*@safe*/ {
return [m, ct].areEqual;
}
 
bool isNormal(T)(immutable Complex!T[][] m, immutable Complex!T[][] ct)
pure nothrow /*@safe*/ {
return [matMul(m, ct), matMul(ct, m)].areEqual;
}
 
auto complexIdentitymatrix(in size_t side) pure nothrow /*@safe*/ {
return side.iota.map!(r => side.iota.map!(c => complex(r == c)).array).array;
}
 
bool isUnitary(T)(immutable Complex!T[][] m, immutable Complex!T[][] ct)
pure nothrow /*@safe*/ {
immutable mct = matMul(m, ct);
immutable ident = mct.length.complexIdentitymatrix;
return [mct, matMul(ct, m), ident].areEqual;
}
 
void main() {
alias C = complex;
immutable x = 2 ^^ 0.5 / 2;
 
foreach (/*immutable*/ const matrix;
[[[C(3.0, 0.0), C(2.0, 1.0)],
[C(2.0, -1.0), C(1.0, 0.0)]],
 
[[C(1.0, 0.0), C(1.0, 0.0), C(0.0, 0.0)],
[C(0.0, 0.0), C(1.0, 0.0), C(1.0, 0.0)],
[C(1.0, 0.0), C(0.0, 0.0), C(1.0, 0.0)]],
 
[[C(x, 0.0), C(x, 0.0), C(0.0, 0.0)],
[C(0.0, -x), C(0.0, x), C(0.0, 0.0)],
[C(0.0, 0.0), C(0.0, 0.0), C(0.0, 1.0)]]]) {
immutable mat = matrix.assumeUnique; //*
 
enum mFormat = "[%([%(%1.3f, %)],\n %)]]";
writefln("Matrix:\n" ~ mFormat, mat);
immutable ct = conjugateTranspose(mat);
"Its conjugate transpose:".writeln;
writefln(mFormat, ct);
writefln("Hermitian? %s.", isHermitian(mat, ct));
writefln("Normal?  %s.", isNormal(mat, ct));
writefln("Unitary?  %s.\n", isUnitary(mat, ct));
}
}

Factor[edit]

Before the fix to Factor bug #484, m. gave the wrong answer and this code failed. Factor 0.94 is too old to work.

Works with: Factor version development (future 0.95)
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 = ;

Usage:

USE: rosetta.hermitian
IN: scratchpad { { C{ 1 2 } 0 }
                 { 0 C{ 3 4 } } }
               [ hermitian-matrix? . ]
               [ normal-matrix? . ]
               [ unitary-matrix? . ] tri
f
t
f

Fortran[edit]

The examples and algorithms are taken from the j solution, except for UnitaryQ. The j solution uses the matrix inverse verb. Compilation on linux, assuming the program is file f.f08 :
gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none f.f08 -o f
 
program conjugate_transpose
 
complex, dimension(3, 3) :: a
integer :: i
a = reshape((/ (i, i=1,9) /), shape(a))
call characterize(a)
a(:2,:2) = reshape((/cmplx(3,0),cmplx(2,-1),cmplx(2,1),cmplx(1,0)/),(/2,2/))
call characterize(a(:2,:2))
call characterize(cmplx(reshape((/1,0,1,1,1,0,0,1,1/),(/3,3/)),0))
a(3,:) = (/cmplx(0,0), cmplx(0,0), cmplx(0,1)/)*sqrt(2.0)
a(2,:) = (/cmplx(0,-1),cmplx(0,1),cmplx(0,0)/)
a(1,:) = (/1,1,0/)
a = a * sqrt(2.0)/2.0
call characterize(a)
 
contains
 
subroutine characterize(a)
complex, dimension(:,:), intent(in) :: a
integer :: i, j
do i=1, size(a,1)
print *,(a(i, j), j=1,size(a,1))
end do
print *,'Is Hermitian? ',HermitianQ(a)
print *,'Is normal? ',NormalQ(a)
print *,'Unitary? ',UnitaryQ(a)
print '(/)'
end subroutine characterize
 
function ct(a) result(b) ! return the conjugate transpose of a matrix
complex, dimension(:,:), intent(in) :: a
complex, dimension(size(a,1),size(a,1)) :: b
b = conjg(transpose(a))
end function ct
 
function identity(n) result(b) ! return identity matrix
integer, intent(in) :: n
real, dimension(n,n) :: b
integer :: i
b = 0
do i=1, n
b(i,i) = 1
end do
end function identity
 
logical function HermitianQ(a)
complex, dimension(:,:), intent(in) :: a
HermitianQ = all(a .eq. ct(a))
end function HermitianQ
 
logical function NormalQ(a)
complex, dimension(:,:), intent(in) :: a
NormalQ = all(matmul(ct(a),a) .eq. matmul(a,ct(a)))
end function NormalQ
 
logical function UnitaryQ(a)
! if A inverse equals A star
! then multiplying each side by A should result in the identity matrix
! Thus show that A times A star is sufficiently close to I .
complex, dimension(:,:), intent(in) :: a
UnitaryQ = all(abs(matmul(a,ct(a)) - identity(size(a,1))) .lt. 1e-6)
end function UnitaryQ
 
end program conjugate_transpose
 
-*- mode: compilation; default-directory: "/tmp/" -*-
Compilation started at Fri Jun  7 16:31:38

a=./f && make $a && time $a
gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none f.f08 -o f
 (  1.00000000    ,  0.00000000    ) (  4.00000000    ,  0.00000000    ) (  7.00000000    ,  0.00000000    )
 (  2.00000000    ,  0.00000000    ) (  5.00000000    ,  0.00000000    ) (  8.00000000    ,  0.00000000    )
 (  3.00000000    ,  0.00000000    ) (  6.00000000    ,  0.00000000    ) (  9.00000000    ,  0.00000000    )
 Is Hermitian?   F
 Is normal?   F
 Unitary?   F


 (  3.00000000    ,  0.00000000    ) (  2.00000000    ,  1.00000000    )
 (  2.00000000    , -1.00000000    ) (  1.00000000    ,  0.00000000    )
 Is Hermitian?   T
 Is normal?   T
 Unitary?   F


 (  1.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    )
 (  1.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    )
 Is Hermitian?   F
 Is normal?   T
 Unitary?   F


 ( 0.707106769    ,  0.00000000    ) ( 0.707106769    ,  0.00000000    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,-0.707106769    ) (  0.00000000    , 0.707106769    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    ) (  0.00000000    , 0.999999940    )
 Is Hermitian?   F
 Is normal?   T
 Unitary?   T



real	0m0.002s
user	0m0.000s
sys	0m0.000s

Compilation finished at Fri Jun  7 16:31:38

Go[edit]

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
}

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

Haskell[edit]

Slow implementation using lists.

import Data.List (transpose)
import Data.Complex
 
type Matrix a = [[a]]
 
main :: IO ()
main =
mapM_ (\a -> do
putStrLn "\nMatrix:"
mapM_ print a
putStrLn "Conjugate Transpose:"
mapM_ print (conjTranspose a)
putStrLn $ "Hermitian? " ++ show (isHermitianMatrix a)
putStrLn $ "Normal? " ++ show (isNormalMatrix a)
putStrLn $ "Unitary? " ++ show (isUnitaryMatrix a))
([[[3, 2:+1],
[2:+(-1), 1 ]],
 
[[1, 1, 0],
[0, 1, 1],
[1, 0, 1]],
 
[[sqrt 2/2:+0, sqrt 2/2:+0, 0 ],
[0:+sqrt 2/2, 0:+ (-sqrt 2/2), 0 ],
[0, 0, 0:+1]]] :: [Matrix (Complex Double)])
 
isHermitianMatrix, isNormalMatrix, isUnitaryMatrix :: RealFloat a => Matrix (Complex a) -> Bool
isHermitianMatrix a = a `approxEqualMatrix` conjTranspose a
isNormalMatrix a = (a `mmul` conjTranspose a) `approxEqualMatrix` (conjTranspose a `mmul` a)
isUnitaryMatrix a = (a `mmul` conjTranspose a) `approxEqualMatrix` ident (length a)
 
approxEqualMatrix :: (Fractional a, Ord a) => Matrix (Complex a) -> Matrix (Complex a) -> Bool
approxEqualMatrix a b = length a == length b && length (head a) == length (head b) &&
and (zipWith approxEqualComplex (concat a) (concat b))
where approxEqualComplex (rx :+ ix) (ry :+ iy) = abs (rx - ry) < eps && abs (ix - iy) < eps
eps = 1e-14
 
mmul :: Num a => Matrix a -> Matrix a -> Matrix a
mmul a b = [[sum (zipWith (*) row column) | column <- transpose b] | row <- a]
 
ident :: Num a => Int -> Matrix a
ident size = [[fromIntegral $ div a b * div b a | a <- [1..size]] | b <- [1..size]]
 
conjTranspose :: Num a => Matrix (Complex a) -> Matrix (Complex a)
conjTranspose = map (map conjugate) . transpose

Output:

Matrix:
[3.0 :+ 0.0,2.0 :+ 1.0]
[2.0 :+ (-1.0),1.0 :+ 0.0]
Conjugate Transpose:
[3.0 :+ (-0.0),2.0 :+ 1.0]
[2.0 :+ (-1.0),1.0 :+ (-0.0)]
Hermitian? True
Normal? True
Unitary? False

Matrix:
[1.0 :+ 0.0,1.0 :+ 0.0,0.0 :+ 0.0]
[0.0 :+ 0.0,1.0 :+ 0.0,1.0 :+ 0.0]
[1.0 :+ 0.0,0.0 :+ 0.0,1.0 :+ 0.0]
Conjugate Transpose:
[1.0 :+ (-0.0),0.0 :+ (-0.0),1.0 :+ (-0.0)]
[1.0 :+ (-0.0),1.0 :+ (-0.0),0.0 :+ (-0.0)]
[0.0 :+ (-0.0),1.0 :+ (-0.0),1.0 :+ (-0.0)]
Hermitian? False
Normal? True
Unitary? False

Matrix:
[0.7071067811865476 :+ 0.0,0.7071067811865476 :+ 0.0,0.0 :+ 0.0]
[0.0 :+ 0.7071067811865476,0.0 :+ (-0.7071067811865476),0.0 :+ 0.0]
[0.0 :+ 0.0,0.0 :+ 0.0,0.0 :+ 1.0]
Conjugate Transpose:
[0.7071067811865476 :+ (-0.0),0.0 :+ (-0.7071067811865476),0.0 :+ (-0.0)]
[0.7071067811865476 :+ (-0.0),0.0 :+ 0.7071067811865476,0.0 :+ (-0.0)]
[0.0 :+ (-0.0),0.0 :+ (-0.0),0.0 :+ (-1.0)]
Hermitian? False
Normal? True
Unitary? True

J[edit]

Solution:
   ct =: [email protected]|:                      NB.  Conjugate transpose (ct A is A_ct)
Examples:
   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
Reference (example matrices for other langs to use):
   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)
 
hermitian=: -: ct
normal =: (X~ -: X) ct
unitary=: ct -: %.
 
(hermitian,normal,unitary)&.>HERMITIAN;NORMAL;UNITARY
+-----+-----+-----+
|1 1 0|0 1 0|0 1 1|
+-----+-----+-----+

jq[edit]

Works with: jq version 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[edit]

(1) transpose/0:

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

# 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 ;

(2) Operations on real/complex numbers

# 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;

(3) Array operations

# 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]) ));

(4) Matrix operations

# 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;

Conjugate transposition[edit]

# (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)  ;

Examples[edit]

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
Output:
$ 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


Julia[edit]

Julia has a built-in matrix type, and the conjugate-transpose of a complex matrix A is simply:

A'

(similar to Matlab). You can check whether A is Hermitian via the built-in function

ishermitian(A)

Ignoring the possibility of roundoff errors for floating-point matrices (like most of the examples in the other languages), you can check whether a matrix is normal or unitary by the following functions

isnormal(A) = size(A,1) == size(A,2) && A'*A == A*A'
isunitary(A) = size(A,1) == size(A,2) && A'*A == eye(A)

Kotlin[edit]

As Kotlin doesn't have built in classes for complex numbers or matrices, some basic functionality needs to be coded in order to tackle this task:

// version 1.1.3
 
typealias C = Complex
typealias Vector = Array<C>
typealias Matrix = Array<Vector>
 
class Complex(val real: Double, val imag: Double) {
 
operator fun plus(other: Complex) =
Complex(this.real + other.real, this.imag + other.imag)
 
operator fun times(other: Complex) =
Complex(this.real * other.real - this.imag * other.imag,
this.real * other.imag + this.imag * other.real)
 
fun conj() = Complex(this.real, -this.imag)
 
/* tolerable equality allowing for rounding of Doubles */
infix fun teq(other: Complex) =
Math.abs(this.real - other.real) <= 1e-14 &&
Math.abs(this.imag - other.imag) <= 1e-14
 
override fun toString() = "${"%.3f".format(real)} " + when {
imag > 0.0 -> "+ ${"%.3f".format(imag)}i"
imag == 0.0 -> "+ 0.000i"
else -> "- ${"%.3f".format(-imag)}i"
}
}
 
fun Matrix.conjTranspose(): Matrix {
val rows = this.size
val cols = this[0].size
return Matrix(cols) { i -> Vector(rows) { j -> this[j][i].conj() } }
}
 
operator fun Matrix.times(other: Matrix): Matrix {
val rows1 = this.size
val cols1 = this[0].size
val rows2 = other.size
val cols2 = other[0].size
require(cols1 == rows2)
val result = Matrix(rows1) { Vector(cols2) { C(0.0, 0.0) } }
for (i in 0 until rows1) {
for (j in 0 until cols2) {
for (k in 0 until rows2) {
result[i][j] += this[i][k] * other[k][j]
}
}
}
return result
}
 
/* tolerable matrix equality using the same concept as for complex numbers */
infix fun Matrix.teq(other: Matrix): Boolean {
if (this.size != other.size || this[0].size != other[0].size) return false
for (i in 0 until this.size) {
for (j in 0 until this[0].size) if (!(this[i][j] teq other[i][j])) return false
}
return true
}
 
fun Matrix.isHermitian() = this teq this.conjTranspose()
 
fun Matrix.isNormal(): Boolean {
val ct = this.conjTranspose()
return (this * ct) teq (ct * this)
}
 
fun Matrix.isUnitary(): Boolean {
val ct = this.conjTranspose()
val prod = this * ct
val ident = identityMatrix(prod.size)
val prod2 = ct * this
return (prod teq ident) && (prod2 teq ident)
}
 
fun Matrix.print() {
val rows = this.size
val cols = this[0].size
for (i in 0 until rows) {
for (j in 0 until cols) {
print(this[i][j])
print(if(j < cols - 1) ", " else "\n")
}
}
println()
}
 
fun identityMatrix(n: Int): Matrix {
require(n >= 1)
val ident = Matrix(n) { Vector(n) { C(0.0, 0.0) } }
for (i in 0 until n) ident[i][i] = C(1.0, 0.0)
return ident
}
 
fun main(args: Array<String>) {
val x = Math.sqrt(2.0) / 2.0
val matrices = arrayOf(
arrayOf(
arrayOf(C(3.0, 0.0), C(2.0, 1.0)),
arrayOf(C(2.0, -1.0), C(1.0, 0.0))
),
arrayOf(
arrayOf(C(1.0, 0.0), C(1.0, 0.0), C(0.0, 0.0)),
arrayOf(C(0.0, 0.0), C(1.0, 0.0), C(1.0, 0.0)),
arrayOf(C(1.0, 0.0), C(0.0, 0.0), C(1.0, 0.0))
),
arrayOf(
arrayOf(C(x, 0.0), C(x, 0.0), C(0.0, 0.0)),
arrayOf(C(0.0, -x), C(0.0, x), C(0.0, 0.0)),
arrayOf(C(0.0, 0.0), C(0.0, 0.0), C(0.0, 1.0))
)
)
 
for (m in matrices) {
println("Matrix:")
m.print()
val mct = m.conjTranspose()
println("Conjugate transpose:")
mct.print()
println("Hermitian? ${mct.isHermitian()}")
println("Normal? ${mct.isNormal()}")
println("Unitary? ${mct.isUnitary()}\n")
}
}
Output:
Matrix:
3.000 + 0.000i,  2.000 + 1.000i
2.000 - 1.000i,  1.000 + 0.000i

Conjugate transpose:
3.000 + 0.000i,  2.000 + 1.000i
2.000 - 1.000i,  1.000 + 0.000i

Hermitian? true
Normal?    true
Unitary?   false

Matrix:
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

Conjugate transpose:
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

Matrix:
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

Conjugate transpose:
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

Maple[edit]

The commands HermitianTranspose and IsUnitary are provided by the LinearAlgebra package.

M:=<<3|2+I>,<2-I|1>>:
 
with(LinearAlgebra):
IsNormal:=A->EqualEntries(A^%H.A,A.A^%H):
 
M^%H;
HermitianTranspose(M);
type(M,'Matrix'(hermitian));
IsNormal(M);
IsUnitary(M);

Output:

                               [  3    2 + I]
                               [            ]
                               [2 - I    1  ]

                               [  3    2 + I]
                               [            ]
                               [2 - I    1  ]

                                    true

                                    true

                                    false

Mathematica / Wolfram Language[edit]

NormalMatrixQ[a_List?MatrixQ] := Module[{b = [email protected]@a},a.b === b.a]
UnitaryQ[m_List?MatrixQ] := ([email protected]@m.m == [email protected]@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)
 
{[email protected]#, [email protected]#, [email protected]#}&@m
-> {False, False, False}

PARI/GP[edit]

conjtranspose(M)=conj(M~)
isHermitian(M)=M==conj(M~)
isnormal(M)=my(H=conj(M~));H*M==M*H
isunitary(M)=M*conj(M~)==1

Perl[edit]

In general, using two or more modules which overload operators can be problematic. For this task, using both Math::Complex and Math::MatrixReal gives us the behavior we want for everything except matrix I/O, i.e. parsing and stringification.

use strict;
use English;
use Math::Complex;
use Math::MatrixReal;
 
my @examples = (example1(), example2(), example3());
foreach my $m (@examples) {
print "Starting matrix:\n", cmat_as_string($m), "\n";
my $m_ct = conjugate_transpose($m);
print "Its conjugate transpose:\n", cmat_as_string($m_ct), "\n";
print "Is Hermitian? ", (cmats_are_equal($m, $m_ct) ? 'TRUE' : 'FALSE'), "\n";
my $product = $m_ct * $m;
print "Is normal? ", (cmats_are_equal($product, $m * $m_ct) ? 'TRUE' : 'FALSE'), "\n";
my $I = identity(($m->dim())[0]);
print "Is unitary? ", (cmats_are_equal($product, $I) ? 'TRUE' : 'FALSE'), "\n";
print "\n";
}
exit 0;
 
sub cmats_are_equal {
my ($m1, $m2) = @ARG;
my $max_norm = 1.0e-7;
return abs($m1 - $m2) < $max_norm; # Math::MatrixReal overloads abs().
}
 
# Note that Math::Complex and Math::MatrixReal both overload '~', for
# complex conjugates and matrix transpositions respectively.
sub conjugate_transpose {
my $m_T = ~ shift;
my $result = $m_T->each(sub {~ $ARG[0]});
return $result;
}
 
sub cmat_as_string {
my $m = shift;
my $n_rows = ($m->dim())[0];
my @row_strings = map { q{[} . join(q{, }, $m-