Re: [R] General expression of a unitary matrix

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 15 Aug 2005 - 04:51:39 EST

          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

J. Liu wrote:

> Thank you, Spencer. I read through the websites you suggested. What I
> need is how to parameterize a 2\times 2 unitary matrix. Generally,
> since for a complex 2\times 2 matrix, there are 8 free variables, and
> for it to be unitary, there are four constraints (unit norm and
> orthogonality), hence I think there are four free variables left for a
> 2\times 2unitary matrix. The form I found can not decribe all the
> unitary matrix, that is why I suspect that it is not the most general
> one. The form in the second web you suggested is an interesting one,
> however, since only 3 variables invovled, it may not be the most
> general expression. 
> 
> Jing  
> 
> 
> On Sat, 13 Aug 2005 09:06:23 -0700
>  Spencer Graves <spencer.graves@pdf.com> wrote:
> 
>>	  Google led me to 
>>"http://mathworld.wolfram.com/SpecialUnitaryMatrix.html", where I 
>>learned that a "special unitary matrix" U has det(U) = 1 in addition
>>to 
>>the "unitary matrix" requirement that
>>
>>	  U %*% t(Conj(U)) == diag(dim(U)[1]).
>>
>>	  Thus, if U is a k x k unitary matrix with det(U) = exp(th*1i), 
>>exp(-th*1i/k)*U is a special unitary matrix.  Moreover, the special 
>>unitary matrices are a group under multiplication.
>>
>>	  Another Google query led me to  	 
>>"http://mathworld.wolfram.com/SpecialUnitaryGroup.html", which gives
>>a 
>>general expression for a special unitary matrix, which seems to
>>require 
>>three real numbers, not four;  with a fourth, you could get a general
>>
>>unitary matrix.
>>
>>	  spencer graves
>>
>>J. Liu wrote:
>>
>>
>>>Hi, all,
>>>
>>>Does anybody got the most general expression of a unitary matrix?
>>>I found one in the book, four entries of the matrix are:
>>> 
>>>(cos\theta) exp(j\alpha);     -(sin\theta)exp(j(\alpha-\Omega));
>>>(sin\theta)exp(j(\beta+\Omega));   (cos\theta) exp(j\beta);
>>> 
>>>where "j" is for complex. 
>>>However, since for any two unitary matrices, their product should
>>
>>also
>>
>>>be a unitary matrix. When I try to use the above expression to
>>>calculate the product, I can not derive the product into the same
>>
>>form.
>>
>>>Therefore, I suspect that this may not be the most general
>>
>>expression. 
>>
>>>Could you help me out of this? Thanks...
>>>
>>>______________________________________________
>>>R-help@stat.math.ethz.ch mailing list
>>>https://stat.ethz.ch/mailman/listinfo/r-help
>>>PLEASE do read the posting guide!
>>
>>http://www.R-project.org/posting-guide.html
>>
>>-- 
>>Spencer Graves, PhD
>>Senior Development Engineer
>>PDF Solutions, Inc.
>>333 West San Carlos Street Suite 700
>>San Jose, CA 95110, USA
>>
>>spencer.graves@pdf.com
>>www.pdf.com <http://www.pdf.com>
>>Tel:  408-938-4420
>>Fax: 408-280-7915
>>
>>______________________________________________
>>R-help@stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide!
>>http://www.R-project.org/posting-guide.html
> 
> 

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Mon Aug 15 05:08:49 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 15:19:33 EST