Hi,

I'm wondering what is the best way to generate uniform random rotation (orthonormal) matrix with a dimension higher than 4 using matlab.

Thanks.

"Martin " <joonyhur@gmail.com> wrote in message <idr87f$ho7$1@fred.mathworks.com>...
> Hi,
> 
> I'm wondering what is the best way to generate uniform random rotation (orthonormal) matrix with a dimension higher than 4 using matlab.
> 

[Q,R] = qr(rand(5));

Consider the properties of Q.

HTH,
John

"Martin " <joonyhur@gmail.com> wrote in message <idr87f$ho7$1@fred.mathworks.com>...
> Hi,
> 
> I'm wondering what is the best way to generate uniform random rotation (orthonormal) matrix with a dimension higher than 4 using matlab.
> 
> Thanks.

What about
[Q, ~] = qr(randn(ndim))

Bruno

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idr8mi$q71$1@fred.mathworks.com>...
> 
> [Q,R] = qr(rand(5));
> 
> Consider the properties of Q.

I'm almost sure using RAND will not produce "uniform" distribution of Q. RANDN will.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrfee$3iu$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idr8mi$q71$1@fred.mathworks.com>...
> > 
> > [Q,R] = qr(rand(5));
> > 
> > Consider the properties of Q.
> 
> I'm almost sure using RAND will not produce "uniform" distribution of Q. RANDN will.
> 
> Bruno

True. The first vector will be a problem.
However, with this slight change, it will do
quite well...

[q,r]  = qr(rand(5)-.5);

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrg69$hrj$1@fred.mathworks.com>...

> 
> [q,r]  = qr(rand(5)-.5);

I still think it's not fix the problem entirely John. I believe the density of q(:,1)=(1,1,1,1,1)'/sqrt(5) is larger than q(:,1)=(1,0,0,0,0)'.

Bruno

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrg69$hrj$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrfee$3iu$1@fred.mathworks.com>...
> > "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idr8mi$q71$1@fred.mathworks.com>...
> > > 
> > > [Q,R] = qr(rand(5));
> > > 
> > > Consider the properties of Q.
> > 
> > I'm almost sure using RAND will not produce "uniform" distribution of Q. RANDN will.
> > 
> > Bruno
> 
> True. The first vector will be a problem.
> However, with this slight change, it will do
> quite well...
> 
> [q,r]  = qr(rand(5)-.5);
> 
> John

On the other hand, is either of these methods
truly random in their result? Try this...

n = 100000;
ang = zeros(n,2);
for i= 1:n
  [q,~] = qr(rand(2) - 0.5);
  ang(i,1) = atan2(q(1,1),q(2,1));
  [q,~] = qr(randn(2) - 0.5);
  ang(i,2) = atan2(q(1,1),q(2,1));
end
hist(ang,100)

John

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrghh$o4t$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrg69$hrj$1@fred.mathworks.com>...
> 
> > 
> > [q,r]  = qr(rand(5)-.5);
> 
> I still think it's not fix the problem entirely John. I believe the density of q(:,1)=(1,1,1,1,1)'/sqrt(5) is larger than q(:,1)=(1,0,0,0,0)'.
> 
> Bruno

True. However, the randn solution also has flaws. See
my other response.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrh6r$7ap$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrghh$o4t$1@fred.mathworks.com>...
> > "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrg69$hrj$1@fred.mathworks.com>...
> > 
> > > 
> > > [q,r]  = qr(rand(5)-.5);
> > 
> > I still think it's not fix the problem entirely John. I believe the density of q(:,1)=(1,1,1,1,1)'/sqrt(5) is larger than q(:,1)=(1,0,0,0,0)'.
> > 
> > Bruno
> 
> True. However, the randn solution also has flaws. See
> my other response.
> 
> John
- - - - - - - -
  Bruno didn't subtract .5 in qr(rand(5)-.5).  That would spoil the uniformity.

Roger Stafford

"John D'Errico" <woodchips@rochester.rr.com> wrote in message 
> On the other hand, is either of these methods
> truly random in their result? Try this...
> > ...
>   [q,~] = qr(randn(2) - 0.5);

The above has flaw, but my original method without shifting

> [q,~] = qr(randn(2));

has not.

Bruno

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <idrhsn$jeg$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrh6r$7ap$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrghh$o4t$1@fred.mathworks.com>...
> > > "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idrg69$hrj$1@fred.mathworks.com>...
> > > 
> > > > 
> > > > [q,r]  = qr(rand(5)-.5);
> > > 
> > > I still think it's not fix the problem entirely John. I believe the density of q(:,1)=(1,1,1,1,1)'/sqrt(5) is larger than q(:,1)=(1,0,0,0,0)'.
> > > 
> > > Bruno
> > 
> > True. However, the randn solution also has flaws. See
> > my other response.
> > 
> > John
> - - - - - - - -
>   Bruno didn't subtract .5 in qr(rand(5)-.5).  That would spoil the uniformity.
> 
> Roger Stafford

Even if you do, neither solution seems to be
truly uniform. Look at my example in 2-d from
my last post.

Perhaps better is this solution...

n = 1000000;
ang = zeros(n,1);
for i= 1:n
  A = randn(2);
  [V,D] = eig(A+A');
  ang(i,1) = atan(V(1,1)/V(2,1));
end
hist(ang,1000)

John

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idri4m$nuc$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message 
> > On the other hand, is either of these methods
> > truly random in their result? Try this...
> > > ...
> >   [q,~] = qr(randn(2) - 0.5);
> 
> The above has flaw, but my original method without shifting
> 
> > [q,~] = qr(randn(2));
> 

Yes it does have a flaw. Have you tried my example
for randn?

LOOK AT THE RESULT. Does it look unbiased?

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idriop$5q2$1@fred.mathworks.com>...

> 
> Yes it does have a flaw. Have you tried my example
> for randn?
> 
> LOOK AT THE RESULT. Does it look unbiased?
> 

Yes I look at the result, unbiased AFAICS, which is what I expected with the theory.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrj41$buh$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <idriop$5q2$1@fred.mathworks.com>...
> 
> > 
> > Yes it does have a flaw. Have you tried my example
> > for randn?
> > 
> > LOOK AT THE RESULT. Does it look unbiased?
> > 
> 
> Yes I look at the result, unbiased AFAICS, which is what I expected with the theory.
> 
> Bruno

Here is what I get: http://yfrog.com/baqrrandnp

Seems pretty uniform to me.

Bruno

Just wonder John, does the QR of your Matlab always output diag(R)>=0 ?

If not, a change should perhaps is nevessary:

[Q, R] = qr(randn(ndim))
Q = Q*diag(sign(diag(R)))

Bruno

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <idriop$5q2$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idri4m$nuc$1@fred.mathworks.com>...
> > The above has flaw, but my original method without shifting
> > 
> > > [q,~] = qr(randn(2));
> > 
> 
> Yes it does have a flaw. Have you tried my example
> for randn?
> 
> LOOK AT THE RESULT. Does it look unbiased?
> 
> John
- - - - - - -
  John, I think there is a communication difficulty here.  In the example you showed with randn you subtracted .5 from it.  That is NOT what Bruno suggested.  He wrote: "[Q, ~] = qr(randn(ndim))" in his first posting.

  In your example where you have ang(i,2) = atan2(q(1,1),q(2,1)) this angle will be the same as atan2(r2,r1) where [r1;r2] are the first column in the qr argument.  With randn alone without .5 subtracted that angle should be uniformly distributed in [0,2*pi].

Roger Stafford

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <idrk4e$t80$1@fred.mathworks.com>...
> Just wonder John, does the QR of your Matlab always output diag(R)>=0 ?
> 
> If not, a change should perhaps is nevessary:
> 
> [Q, R] = qr(randn(ndim))
> Q = Q*diag(sign(diag(R)))
> 
> Bruno

Hi Bruno,

How did you come up with this formula for uniform random rotation matrices? Is there a proof or some literature available?

Thanks!

"G" wrote in message <jl9ktr$f4c$1@newscl01ah.mathworks.com>...
> 
> How did you come up with this formula for uniform random rotation matrices? Is there a proof or some literature available?

A proof, probably. The fact is straight from the property of rotational invariant of canonical multivariate Gaussian pdf.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jlab6e$s2m$1@newscl01ah.mathworks.com>...
> "G" wrote in message <jl9ktr$f4c$1@newscl01ah.mathworks.com>...
> > 
> > How did you come up with this formula for uniform random rotation matrices? Is there a proof or some literature available?
> 
> A proof, probably. The fact is straight from the property of rotational invariant of canonical multivariate Gaussian pdf.
> 
> Bruno

Hi,

The code provided in this thread:

[Q, R] = qr(randn(n))
Q = Q*diag(sign(diag(R)))

Is a correct way to generate random matrices from O(n), the orthogonal group. For reference, see "How to generate random matrices from the classical compact groups", Francesco Mezzadri, Notices of ACM, Volume 54, Number 5, 2007.

The problem is that det(Q) may be -1, and therefore sometimes Q is not a rotation matrix as desired. I'm thinking of two possibilities for det(Q) < 0:

1) exchange any two columns of Q, or
2) compute the svd of Q = USV^T, and then compute Q2 := U diag([ones(1, n - 1), -1]) V^T. The Q2 is then the rotation matrix closest to Q in the Frobenius norm.

I am unable to say what effect these have on the probability measure. Anyone?

Kalle

Pick any column of Q then reverse the sign when det(Q)=-1. You'll be fine.

Bruno