Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Uniform random rotation matrix

Subject: Uniform random rotation matrix

From: Martin

Date: 9 Dec, 2010 18:47:11

Message: 1 of 20

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.

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 18:55:14

Message: 2 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 18:57:05

Message: 3 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 20:50:22

Message: 4 of 20

"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

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 21:03:05

Message: 5 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 21:09:05

Message: 6 of 20

"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

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 21:13:20

Message: 7 of 20

"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

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 21:20:27

Message: 8 of 20

"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

Subject: Uniform random rotation matrix

From: Roger Stafford

Date: 9 Dec, 2010 21:32:07

Message: 9 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 21:36:23

Message: 10 of 20

"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

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 21:45:24

Message: 11 of 20

"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

Subject: Uniform random rotation matrix

From: John D'Errico

Date: 9 Dec, 2010 21:47:05

Message: 12 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 21:53:05

Message: 13 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 21:59:07

Message: 14 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 9 Dec, 2010 22:10:22

Message: 15 of 20

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

Subject: Uniform random rotation matrix

From: Roger Stafford

Date: 9 Dec, 2010 22:12:05

Message: 16 of 20

"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

Subject: Uniform random rotation matrix

From: G

Date: 1 Apr, 2012 13:24:11

Message: 17 of 20

"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!

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 1 Apr, 2012 19:44:15

Message: 18 of 20

"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

Subject: Uniform random rotation matrix

From: Kaba

Date: 25 Oct, 2012 19:29:08

Message: 19 of 20

"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

Subject: Uniform random rotation matrix

From: Bruno Luong

Date: 25 Oct, 2012 22:51:08

Message: 20 of 20

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

Bruno

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us