Question: Efficient parametrization of unitary matrices

Hello Maple community, although I'm using Maple (now version 11) for quite a while now, I still have problems in designing and writing efficient procedures. Within the context of quantum information theory, I want to write a procedure that parametrizes unitary matrices of arbitrary dimension NxN (in order to perform an optimization over this space). For the beginning, I just want to create random unitary matrices of a given dimension. For this, I implemented the formulas for the "Generalized Euler angle parametrization" (math-ph/0205016). This involves first the creation of N^2-1 matrices that serve as generators for SU(N) and, in a second step, several matrix exponentials and matrix products. These two steps I put in two procedures ("hermitian_basis" and "parametrize_unitary_Euler") (see the linked worksheet). Download 711_random unitary.mws
View file details In principle, this seems to work but, unfortunately, it is really slow (so that one cannot use it inside some optimization loop later!). Interestingly, I saw a very similar implementation for Matlab which also uses this parametrization for unitary matrices (available from http://www.qlib.info). On my laptop (Intel CoreDuo 2.3GHz, 2GB RAM, running Win XP), the Matlab code seems to return a random unitary almost instantly where my Maple code takes several seconds, i.e. we are talking about roughly one order of magnitude difference in speed! So the question is: What did I do wrong? I already used caching of the generator matrices (and store them as sparse). I used evalf wherever it seemed appropriate for me. Since everything involving matrices cannot be compiled so easily, I guess this is also no solution... So, how can I remove stupid bottlenecks? I'm grateful for any hint. Thanks in advance!
Please Wait...