# Generating unitary matrices numerically - “close” to the identity element

EDIT: broke this into two parts - for these were two different questions.

For numerically obtaining the stabilities of a matricial equation, i need to generate an ensemble of matrices that are members of SU(n), and strategies for uniform random generation of those are outlined here. But I also need these unitary matrices to be clustered close to the SU(n) identity operation, at some fixed "distance". Thus my question distills to -

1. What might be a good metric for SU(n) matrices? EDIT: one option has been provided by @Jim_Belk downstairs (:D) - certainly feel free to add more!

2. How does one generate SU(n) matrices close to the Identity?

EDIT: are there any alternatives other than computing the generators first and then exponentiating them?

If you find the questions framed improperly or dubious, do tell!

-

If you want to generate random elements of $SU(n)$ close to the identity, it seems to me that the most logical approach would be to generate random elements of the Lie algebra $\mathfrak{su}(n)$, and then exponentiate them.

The Lie algebra $\mathfrak{su}(n)$ consists of all $n\times n$ matrices $A$ with trace $0$ satisfying $A^\dagger = -A$. In particular, the diagonal entries of $A$ must be pure imaginary and add up to $0$, and the entries below the diagonal are entirely determined by the entries above the diagonal. Given a matrix $A\in\mathfrak{su}(n)$, the exponential $\exp(tA)$ is a unitary matrix for all $t\in\mathbb{R}$.

There is a standard inner product on $\mathfrak{su}(n)$, namely the negative of the Killing form, which is invariant under the adjoint action. This is given by the formula $$\langle A, B\rangle = -2n\,\mathrm{tr}(AB) = -2n\sum_{i=1}^n\sum_{j=1}^n a_{ij} b_{ji}$$ We will drop the factor of $2n$. Since $b_{ji} = -b_{ij}^*$, this can also be written $$\langle A, B\rangle = \sum_{i=1}^n\sum_{j=1}^n a_{ij} b_{ij}^* \;=\; \sum_{i=1}^n \mathrm{Im}(a_{ii})\mathrm{Im}(b_{ii}) + \sum_{i<j} 2 \bigl[\mathrm{Re}(a_{ij})\mathrm{Re}(b_{ij})+\mathrm{Im}(a_{ij})\mathrm{Im}(b_{ij})\bigr]$$ In particular, $$\langle A,A\rangle \;=\; \sum_{i=1}^n |a_{ii}|^2 + \sum_{i<j} 2|a_{ij}|^2$$ Presumably the right thing to do is to choose random elements $A$ of $\mathfrak{su}(n)$ having unit norm with respect to this inner product, and then compute $\exp(tA)$ for some fixed small value of $t$. To choose random elements with unit norm, you should start by finding an orthonormal basis for $\mathfrak{su}(n)$ (which should be fairy easy given the above formulas, except for the added complication that the sum of the diagonal entries is $0$, which can be solved with Gram-Schmidt), and then choose coefficients for these basis vectors in the same way that you would choose the coordinates of a random point on the unit sphere in $\mathbb{R}^{n^2-1}$.

-
Apparently I dont have enough reps to vote this useful - but it is! – Debanjan Basu Apr 10 '12 at 8:57
However, could you elucidate a bit on that simplification made from $$\sum_{i.j=1}^n a_{ij}b_{ij}= \sum_{i=1}^n \Im(a_{ii})^2 + \sum_{i\le j}2\left[\Re(a_{ij})\Re(b_{ij})+\Im(a_{ij})\Im(b_{ij})\right]$$? For starters i am pretty sure that the first term in the RHS would contain one each of a and b. – Debanjan Basu Apr 10 '12 at 9:14
But however I have already indicated a paper that generates SU(n) matrices for Uniform Random Sampling. That approach employs generating random matrices or elements of a Ginibre Ensemble from a Gaussian Variate(0,1) and Gram-Schmidts it subsequently. Why that works is explained in this essay‌​. – Debanjan Basu Apr 10 '12 at 9:14
@Debanjan Quite right -- I have fixed the formula. The main trick is to combine the $a_{ij}b_{ij}^*$ term with the $a_{ji}b_{ji}^*$ term of the sum. These two are complex conjugates so the result is $2\mathrm{Re}(a_{ij}b_{ij}^*)$, which is equal to the indicated expansion. I wrote it the way I did to emphasize the similarity to dot product if you think of the real and imaginary parts of the entries above the diagonal as separate coordinates. – Jim Belk Apr 10 '12 at 16:03

A cheap and cheerful way to get a (special) unitary matrix close to the identity is to use the product of two Householder reflections, where the unit vectors defining them are close. I reckon the Frobenius norm of the difference between the product and the identity is proportional to the absolute value of the sine of the angle between the vectors.

-
okay! once i get it read, i will edit this answer - others may find it a more ready reference! – Debanjan Basu Apr 12 '12 at 13:51
so i know what householder reflections are - but i dont get how one generates unitary matrices from it! – Debanjan Basu Apr 17 '12 at 6:03