怎么随机生成正交矩阵$M\in O(n)$

hbghlyj

怎么随机生成正交矩阵$M\in O(n)$


web.archive.org/web/20220121021352/https://st … ubgroup-rand-var.pdf
We suggest a simple algorithm for Monte Carlo generation of uniformly distributed variables on a compact group. Examples include random permutations, Rubik's cube positions, orthogonal, unitary, and symplectic matrices, and elements of $G L_n$ over a finite field. The algorithm reduces to the "standard" fast algorithm when there is one, but many new examples are included.

hbghlyj

概述了一种均匀随机 3D 旋转矩阵的方法，步骤如下：

• 绕 z 轴随机旋转（z 单位向量未移动）
• 将z 单位向量旋转到单位球面上的随机位置

其中第一步是具有随机角度的简单 2D 旋转矩阵，参数为$x_1 \sim U(0, 1)$
$$R=\left[\begin{array}{ccc}
\cos \left(2 \pi x_1\right) & \sin \left(2 \pi x_1\right) & 0 \\
-\sin \left(2 \pi x_1\right) & \cos \left(2 \pi x_1\right) & 0 \\
0 & 0 & 1
\end{array}\right]$$
第二步比较棘手：它使用某些 Householder 矩阵，如下所示：
$$H=I-2 v v^{\mathrm{\top}}$$
$$v=\left[\begin{array}{c}
\cos \left(2 \pi x_2\right) \sqrt{x_3} \\
\sin \left(2 \pi x_2\right) \sqrt{x_3} \\
\sqrt{1-x_3}
\end{array}\right]$$
其中 $x_2$ 和 $x_3$ 也取自 $U(0, 1)$。最终完成的随机 3D 旋转矩阵为（包括通过原点的反射，将 Householder 矩阵变换从反射转换为旋转）：$$M=-H R$$

hbghlyj

Python

1. import numpy as np
3. def uniform_random_rotation(x):
4.     """Apply a random rotation in 3D, with a distribution uniform over the
5.     sphere.
7.     Arguments:
8.         x: vector or set of vectors with dimension (n, 3), where n is the
9.             number of vectors
11.     Returns:
12.         Array of shape (n, 3) containing the randomly rotated vectors of x,
13.         about the mean coordinate of x.
15.     Algorithm taken from "Fast Random Rotation Matrices" (James Avro, 1992):
16.     https://doi.org/10.1016/B978-0-08-050755-2.50034-8
17.     """
19.     def generate_random_z_axis_rotation():
20.         """Generate random rotation matrix about the z axis."""
21.         R = np.eye(3)
22.         x1 = np.random.rand()
23.         R[0, 0] = R[1, 1] = np.cos(2 * np.pi * x1)
24.         R[0, 1] = -np.sin(2 * np.pi * x1)
25.         R[1, 0] = np.sin(2 * np.pi * x1)
26.         return R
28.     # There are two random variables in [0, 1) here (naming is same as paper)
29.     x2 = 2 * np.pi * np.random.rand()
30.     x3 = np.random.rand()
32.     # Rotation of all points around x axis using matrix
33.     R = generate_random_z_axis_rotation()
34.     v = np.array([
35.         np.cos(x2) * np.sqrt(x3),
36.         np.sin(x2) * np.sqrt(x3),
37.         np.sqrt(1 - x3)
38.     ])
39.     H = np.eye(3) - (2 * np.outer(v, v))
40.     M = -(H @ R)
41.     x = x.reshape((-1, 3))
42.     mean_coord = np.mean(x, axis=0)
43.     return ((x - mean_coord) @ M) + mean_coord @ M

Copy the Code

Mathematica

1. uniformRandomRotation[x_] := Module[
2.   {generateRandomZAxisRotation, x2, x3, R, v, H, M, meanCoord},
3.   
4.   generateRandomZAxisRotation[] := Module[{m, r},
5.     r = RandomReal[];
6.     m = IdentityMatrix[3];
7.     m[[1,1]] = m[[2,2]] = Cos[2 Pi r];
8.     m[[1,2]] = -Sin[2 Pi r];
9.     m[[2,1]] = Sin[2 Pi r];
10.     m
11.   ];
12.   
13.   x2 = 2 Pi RandomReal[];
14.   x3 = RandomReal[];
15.   R = generateRandomZAxisRotation[];
16.   v = {Cos[x2] Sqrt[x3], Sin[x2] Sqrt[x3], Sqrt[1 - x3]};
17.   H = IdentityMatrix[3] - 2 Outer[Times, v, v];
18.   M = - (H . R);
19.   meanCoord = Mean[x];
20.   
21.   ((x - meanCoord) M) + meanCoord M
22. ]

Copy the Code

hbghlyj

Sampling a Uniformly Random Rotation
这个是什么原理呢

1. GenUniformRandRotation[dp_,dd_,xn_]:=Module[{x,\[Theta],\[Phi],z,r,Vx,Vy,Vz,st,ct,Sx,Sy},
2. x=xn;(*RandomReal[1,3];*)
3. \[Theta]=(x[[1]]-1/2)dp 2\[Pi];(*Rotation about the pole (Z). This is improved from Arvo's implementation by subtracting 1/2 so that the rotation is unbiased. *)
4. \[Phi]=x[[2]]2\[Pi];(*For direction of pole deflection. *)
5. z=x[[3]]dd 2;(*For magnitude of pole deflection. *)(*Compute a vector V used for distributing points over the sphere via the reflection I-V Transpose(V).This formulation of V will guarantee that if x[1] and x[2] are uniformly distributed, the reflected points will be uniform on the sphere.Note that V has length Sqrt[2] to eliminate the 2 in the Householder matrix. *)
6. r=Sqrt[z];
7. Vx=Sin[\[Phi]]r;
8. Vy=Cos[\[Phi]]r;
9. Vz=Sqrt[2-z];
10. (*Compute the row vector S=Transpose(V)*R, where R is a simple rotation by theta about the z-axis. No need to compute Sz since it's just Vz. *)
11. st=Sin[\[Theta]];
12. ct=Cos[\[Theta]];
13. Sx=Vx ct-Vy st;
14. Sy=Vx st+Vy ct;
16. (*Construct the rotation matrix (V Transpose(V)-I) R *R_z[\[Pi]], which is equivalent to V S-R.
17. Note that Arvo's code is missing the multiplication by R_z[\[Pi]], which is unnoticeable for random rotations, but causes problems for random perturbations.*)
18. M=(-Vx Sx+ct        -Vx Sy+st        Vx Vz
19. -Vy Sx-st        -Vy Sy+ct        Vy Vz
20. -Vz Sx        -Vz Sy        1-z
22. )]
23. sc = 0.02;
24. maxNumSamples = 1000;
25. xns = RandomReal[1,{3,maxNumSamples}]; (*precompute random numbers so that the Manipulate does not 'bounce'*)
26. (*Test code:
27. MatrixForm[GenUniformRandRotation[0,0,RandomReal[1,{3,1}]]]
28. Rz[\[Theta]_] := (\[NoBreak]Cos[\[Theta]]        Sin[\[Theta]]        0
29. -Sin[\[Theta]]        Cos[\[Theta]]        0
30. 0        0        1
32. \[NoBreak])
33. A = (\[NoBreak]a        b        c
34. d        e        f
35. g        h        j
37. \[NoBreak]);
38. A.Rz[\[Pi]]
39. Rz[\[Pi]].A*)

Copy the Code