General expression of a unitary matrix

Below please find functions that attempt to test for whether a matrix 
is unitary and special unitary (SU) and generate an SU matrix from a 
3-vector and convert a 2x2 SU matrix to a 3-vector.  These are not 
extensively debugged, so they may not be correct.  However, they passed 
a few simple tests, including the following:

	  SU2.4 <- SU2(1:3)%*%SU2(2:4)
	  all.equal(SU2.4, SU2(Re(SU2.vec(SU2.4))))

	  An SU2 matrix has lots of structure.  I computed xi from the ratio of 
the diagonals and zeta from the ratio of the off-diagonals.  Then I used 
those to compute eta = atan(x[1,2]/exp(zeta*1i), x[1,1]/exp(xi*1i)).

	  If you find a counter example, please post it to the list;  maybe 
someone will fix it or explain why Wolfram was wrong.  If you clean up 
the functions I wrote, e.g, adding checks and returning only the real 
part of the 3-vector, etc., please post that to the list.

	  spencer graves
######################################
is.unitary <- function(x,
        eps=prod(dim(x))*.Machine$double.eps){
   x2 <- (x %*% t(Conj(x)))
   (abs(mean(x2-diag(dim(x)[1])))
       < eps)
}

is.unitary(diag(3)+1e-15)
is.unitary(diag(3)+1e-16)

is.SU <- function(x,
        eps=prod(dim(x))*.Machine$double.eps){
   if(is.unitary(x, eps)){
     eig.x <- eigen(x)
     det.x <- prod(eig.x$values)
     return(abs(det.x-1)<eps)
   }
   else FALSE
}

is.SU(diag(3)+1e-15)
is.SU(diag(3)+1e-16)

SU2 <- function(x){
# x = c(xi, eta, zeta)
   eix <- exp(1i*x[1])
   eiz <- exp(1i*x[3])
   su2 <- array(NA, dim=c(2,2))
   diag(su2) <- (c(eix, Conj(eix))*
             cos(x[2]))
   seta <- sin(x[2])
   su2[1,2] <- eiz*seta
   su2[2,1] <- (-Conj(eiz)*seta)
   su2
}

SU2(1:3)
is.SU(SU2(1:3))

is.SU(SU2(1:3)%*%SU2(2:4))

SU2.vec <- function(x,
        eps=prod(dim(x))*.Machine$double.eps){
   xi <- (log(x[1,1]/x[2,2])/(2i))
   zeta <- (log(-x[1,2]/x[2,1])/(2i))
#
   eixi <- exp(xi*1i)
   eizeta <- exp(zeta*1i)
   eta1 <- atan(x[1,2]/eizeta,
                x[1,1]/eixi )
#  eta2 <- atan(-x[2,1]/eizeta,
#               x[2,2]/eixi)
   vec <- c(xi, eta1, zeta)
   vec
}

x <- SU2(1:3)

SU2.4 <- SU2(1:3)%*%SU2(2:4)

SU2.vec(SU2.4)

SU2.4
all.equal(SU2.4,
SU2(SU2.vec(SU2.4)))

all.equal(SU2.4,
SU2(Re(SU2.vec(SU2.4))))

SU2.vec(SU2(1:3))

######################################
Yeah. Let U=U_1*U_2' where U_1 and U_2 are unitary in that form. My
objection is to write U in that form too. However, I can not find a way
to do it.

On Sun, 14 Aug 2005 09:05:19 -0700
#########################
	  Could you provide an example that can NOT be expressed in that form?

	  spencer graves

General expression of a unitary matrix

Thread (5 messages)