精华内容
下载资源
问答
  • 旋转矩阵 随机生成
    千次阅读
    2021-01-13 18:46:49

    在下面的代码中,我对一般的平方线性系统Ax=b实现了带有部分旋转的高斯消去。我测试了我的代码,它产生了正确的输出。不过,现在我正在尝试做以下事情,但我不太确定如何编码它,寻找一些帮助与此!

    我想通过求解Ax=b来测试我的实现,其中A是随机的100x100矩阵,b是随机的100x10向量。

    在我的代码中,我把矩阵

    A = np.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])

    b = np.array([[3.],[15.],[14.]])

    得到以下正确的输出:

    [3. 1. 2.]

    [3. 1. 2.]

    但是现在我如何改变它来生成随机矩阵呢?

    下面是我的代码:

    import numpy as np

    def GEPP(A, b, doPricing = True):

    '''

    Gaussian elimination with partial pivoting.

    input: A is an n x n numpy matrix

    b is an n x 1 numpy array

    output: x is the solution of Ax=b

    with the entries permuted in

    accordance with the pivoting

    done by the algorithm

    post-condition: A and b have been modified.

    '''

    n = len(A)

    if b.size != n:

    raise ValueError("Invalid argument: incompatible sizes between"+

    "A & b.", b.size, n)

    # k represents the current pivot row. Since GE traverses the matrix in the

    # upper right triangle, we also use k for indicating the k-th diagonal

    # column index.

    # Elimination

    for k in range(n-1):

    if doPricing:

    # Pivot

    maxindex = abs(A[k:,k]).argmax() + k

    if A[maxindex, k] == 0:

    raise ValueError("Matrix is singular.")

    # Swap

    if maxindex != k:

    A[[k,maxindex]] = A[[maxindex, k]]

    b[[k,maxindex]] = b[[maxindex, k]]

    else:

    if A[k, k] == 0:

    raise ValueError("Pivot element is zero. Try setting doPricing to True.")

    #Eliminate

    for row in range(k+1, n):

    multiplier = A[row,k]/A[k,k]

    A[row, k:] = A[row, k:] - multiplier*A[k, k:]

    b[row] = b[row] - multiplier*b[k]

    # Back Substitution

    x = np.zeros(n)

    for k in range(n-1, -1, -1):

    x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]

    return x

    if __name__ == "__main__":

    A = np.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])

    b = np.array([[3.],[15.],[14.]])

    print (GEPP(np.copy(A), np.copy(b), doPricing = False))

    print (GEPP(A,b))

    更多相关内容
  • (本代码已作废,但无法删除该资源,非常抱歉,请查看新的文件)注:非常抱歉,由于考虑欠周全,本代码适用范围确实较窄,为弥补这一缺陷,认真考虑后重新编写代码,详见旋转矩阵到四元数源代码(新),如有不当之处...
  • 计算向量X=(x1,x2,⋯ ,xn)X=(x_1,x_2,\cdots,x_n)X=(x1​,x2​,⋯,xn​)绕V=(v1,v2,⋯ ,vn)V=(v_1,v_2,\cdots,v_n)V=(v1​,v2​,⋯,vn​)顺时针旋转θ\thetaθ后得到的向量X′X'X′ 任意的旋转都可以看作绕着一个...

    本来想做绕轴旋转,绕轴旋转是在垂直于轴向量的空间里旋转,但是n维空间里与某个向量垂直的空间为n-1维,而旋转只在二维空间里有定义。所以这里就改成了,向任意方向旋转。

    计算单位向量 X = ( x 1 , x 2 , ⋯   , x n ) X=(x_1,x_2,\cdots,x_n) X=(x1,x2,,xn)旋转到单位向量 V = ( v 1 , v 2 , ⋯   , v n ) V=(v_1,v_2,\cdots,v_n) V=(v1,v2,,vn)的旋转矩阵 R R R

    旋转坐标系

    任意的旋转都可以看作绕着一个轴,在某个平面上的旋转。
    不失一般性,假定向量 V = ( v 1 , v 2 , ⋯   , v n ) V=(v_1,v_2,\cdots,v_n) V=(v1,v2,,vn) v 1 × v 2 v_1\times v_2 v1×v2张成的平面上旋转,用矩阵乘法可表示为: V ′ = R 1 V V'=R_1V V=R1V其中旋转矩阵 R 1 R_1 R1定义为
    R 1 = ( cos ⁡ α − sin ⁡ α sin ⁡ α cos ⁡ α 1 ⋱ 1 ) R_1=\left( \begin{array}{ccc} \cos\alpha& -\sin\alpha& & & \\ \sin\alpha& \cos\alpha& & & \\ & &1 & & \\ & & & \ddots & \\ & & & & 1 \\ \end{array} \right) R1=cosαsinαsinαcosα11
    其中 cos ⁡ α = v 2 v 1 2 + v 2 2 \cos\alpha=\frac{v_2}{\sqrt{v_1^2+v_2^2}} cosα=v12+v22 v2 sin ⁡ α = v 1 v 1 2 + v 2 2 \sin\alpha=\frac{v_1}{\sqrt{v_1^2+v_2^2}} sinα=v12+v22 v1
    由此我们将向量 V V V旋转到与 ( v 1 , 0 , ⋯   , 0 ) (v_1,0,\cdots,0) (v1,0,,0)正交。得 V ′ = ( 0 , v 1 2 + v 2 2 , v 3 , v 4 , ⋯   , v n ) V'=(0,\sqrt{v_1^2+v_2^2}, v_3,v_4,\cdots,v_n) V=(0,v12+v22 ,v3,v4,,vn)
    以此类推,将 V V V旋转到 ( 0 , ⋯   , 0 , 1 ) (0,\cdots,0,1) (0,,0,1)需要 n − 1 n-1 n1次旋转。将这些旋转变换组合为一个旋转矩阵 R 0 = R n − 1 R n − 2 ⋯ R 1 R_0=R_{n-1}R_{n-2}\cdots R_1 R0=Rn1Rn2R1(注意,乘法不满足交换律,顺序不能颠倒)。顺便我们可以计算出这个旋转的逆矩阵 R 0 − 1 = R 1 ⊤ ⋯ R n − 2 ⊤ R n − 1 ⊤ R_0^{-1}=R_1^\top\cdots R_{n-2}^\top R_{n-1}^\top R01=R1Rn2Rn1

    旋转向量

    现在我们要做的是将 R 0 X R_0X R0X旋转到坐标轴 ( 0 , ⋯   , 0 , 1 ) (0,\cdots,0,1) (0,,0,1),并求旋转矩阵。做法与上面相同。最后可以求得,在旋转后的坐标系下,旋转矩阵为 R X R_X RX: R 0 V = R X R 0 X R_0V=R_XR_0X R0V=RXR0X

    虽然 R R R R X R_X RX表示的旋转角度相同,但是这两个旋转所在的坐标系是不同的,还需要将 R X R_X RX的坐标轴旋转回去:
    R = R 0 − 1 R X R 0 R=R_0^{-1}R_XR_0 R=R01RXR0

    代码实现

    import numpy as np
    
    def Rot_map(V):
    	assert len(V.shape) == 1
    	assert np.linalg.norm(V) - 1 < 1e-8
    	n_dim = V.shape[0]
        Rot = np.eye(n_dim)
        Rot_inv = np.eye(n_dim)
        for rotate in range(n_dim-1):
            rot_mat = np.eye(n_dim)
            rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
            cos_theta = V[rotate+1]/rot_norm
            sin_theta = V[rotate]/rot_norm
            rot_mat[rotate,rotate] = cos_theta
            rot_mat[rotate,rotate+1] = - sin_theta
            rot_mat[rotate+1,rotate] = sin_theta
            rot_mat[rotate+1,rotate+1] = cos_theta
    
            V = np.dot(rot_mat, V)
    
            Rot = np.dot(rot_mat, Rot)
            Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
        return Rot, Rot_inv
        
    n_dim = 512
    
    X = np.random.rand(n_dim)
    X = X / np.linalg.norm(X)
    V = np.random.rand(n_dim)
    V = V / np.linalg.norm(V)
    
    R_0, R_0_inv = Rot_map(V)
    R_X, _ = Rot_map(np.dot(R_0, X))
    
    R = np.dot(np.dot(R_0_inv, R_X), R_0)
    
    assert np.linalg.norm(np.dot(R,X) - V) < 1e-8
    

    扩展

    这个方法可以用来在正交于某个向量的空间中随机采样。例如,在 n n n维空间中,采样与 V = ( v 1 , v 2 , ⋯   , v n ) V=(v_1,v_2,\cdots,v_n) V=(v1,v2,,vn)正交的向量,则代码如下:

    import numpy as np
    
    def Rot_map(V):
    	assert len(V.shape) == 1
    	assert np.linalg.norm(V) - 1 < 1e-8
        z_dim = V.shape[0]
        Rot = np.eye(z_dim)
        Rot_inv = np.eye(z_dim)
        for rotate in range(z_dim-1):
            rot_mat = np.eye(z_dim)
            rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
            cos_theta = V[rotate+1]/rot_norm
            sin_theta = V[rotate]/rot_norm
            rot_mat[rotate,rotate] = cos_theta
            rot_mat[rotate,rotate+1] = - sin_theta
            rot_mat[rotate+1,rotate] = sin_theta
            rot_mat[rotate+1,rotate+1] = cos_theta
    
            V = np.dot(rot_mat, V)
    
            Rot = np.dot(rot_mat, Rot)
            Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
        return Rot, Rot_inv
        
    z_dim = 512
    n_samples = 10000
    
    V = np.random.rand(z_dim)
    V = V / np.linalg.norm(V)
    
    R_0, R_0_inv = Rot_map(V)
    
    samples = np.random.rand(z_dim-1, n_samples)
    samples = np.concatenate((samples, np.zeros((1, n_samples))))
    samples = np.dot(R_0_inv, samples)
    
    assert np.mean(np.abs(np.dot(V, samples))) < 1e-8
    

    进一步,我们采样与一组向量 V s Vs Vs正交的向量:

    import numpy as np
    
    
    def gs(X, row_vecs=True, norm = True):
    	'''
    	https://gist.github.com/iizukak/1287876/edad3c337844fac34f7e56ec09f9cb27d4907cc7
    	'''
        if not row_vecs:
            X = X.T
        Y = X[0:1,:].copy()
        for i in range(1, X.shape[0]):
            proj = np.diag((X[i,:].dot(Y.T)/np.linalg.norm(Y,axis=1)**2).flat).dot(Y)
            Y = np.vstack((Y, X[i,:] - proj.sum(0)))
        if norm:
            Y = np.diag(1/np.linalg.norm(Y,axis=1)).dot(Y)
        if row_vecs:
            return Y
        else:
            return Y.T
    
    
    def Rot_map(V):
    	assert len(V.shape) == 1
    	assert np.linalg.norm(V) - 1 < 1e-8
        z_dim = V.shape[0]
        Rot = np.eye(z_dim)
        Rot_inv = np.eye(z_dim)
        for rotate in range(z_dim-1):
            rot_mat = np.eye(z_dim)
            rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
            cos_theta = V[rotate+1]/rot_norm
            sin_theta = V[rotate]/rot_norm
            rot_mat[rotate,rotate] = cos_theta
            rot_mat[rotate,rotate+1] = - sin_theta
            rot_mat[rotate+1,rotate] = sin_theta
            rot_mat[rotate+1,rotate+1] = cos_theta
    
            V = np.dot(rot_mat, V)
    
            Rot = np.dot(rot_mat, Rot)
            Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
        return Rot, Rot_inv
        
    z_dim = 10 # number of dimensions
    n_vectors= 2 # number of Vs
    assert n_vectors < z_dim
    n_samples = 10 # number of samples to orthogonal Vs
    
    Vs = np.random.rand(n_vectors, z_dim)
    # Gram-Schmidt Orthogonization
    orth_Vs = gs(Vs)
    orth_Vs = orth_Vs.T
    
    Rots = np.eye(z_dim)
    Rot_invs = np.eye(z_dim)
    for idx in range(n_vectors):
        V = orth_Vs[:, idx]
        R_0 = np.eye(z_dim)
        R_0_inv = np.eye(z_dim)
        R_0[0:z_dim-idx,0:z_dim-idx], R_0_inv[0:z_dim-idx,0:z_dim-idx] = Rot_map(V[range(z_dim-idx)])
        
        orth_Vs = R_0.dot(orth_Vs)
        Rots = np.dot(R_0, Rots)
        Rot_invs = np.dot(Rot_invs, R_0_inv)
    
    samples = np.random.rand(z_dim-n_vectors, n_samples)
    samples = np.concatenate((samples, np.zeros((n_vectors, n_samples))))
    samples = np.dot(Rot_invs, samples)
    
    assert np.mean(np.abs(np.dot(Vs, samples))) < 1e-8
    

    用矩阵旋转太慢了,上述采样操作可以用矩阵分解解决(在任意向量的正交空间中采样):

    import numpy as np
    
    z_dim = 10 # number of dimensions
    n_vectors= 2 # number of Vs
    assert n_vectors < z_dim
    n_samples = 20 # number of samples to orthogonal Vs
    
    Vs = np.random.rand(n_vectors, z_dim)
    Vs = np.concatenate((Vs,np.zeros((z_dim-n_vectors, z_dim))))
    
    U, s, V = np.linalg.svd(Vs)
    
    directions = np.eye(z_dim)
    directions[range(n_vectors),range(n_vectors)] = 0
    directions = np.dot(np.dot(U, directions), V)
    
    samples = np.random.rand(z_dim, n_samples)
    samples = np.dot(directions.T, samples)
    
    assert np.mean(np.abs(np.dot(Vs, samples))) < 1e-8
    
    展开全文
  • 其实您错了,世界上并没有绝对理想的随机数,就算是用电脑也只能产生接近随机的数:在彩票开奖产生号码时,诸如摇奖机的物理特性、每个球的重量和光滑度的差异、空气的流动性等等都会使开奖结果产生一定的偏态,在中...
  • 常用的旋转参数化方式有轴角、旋转矩阵、欧拉角、四元数等,它们之间的转换推导可以查看这里。 不同旋转表示方式之间的转换在网上可以找到很多相关的代码,同时也有一些库帮助我们实现了它们之间的转换,比如python...

      在做三维相关工作的时候,经常会遇到需要在不同旋转表示方式之间进行转换的情况。常用的旋转参数化方式有轴角、旋转矩阵、欧拉角、四元数等,它们之间的转换推导可以查看这里
      不同旋转表示方式之间的转换在网上可以找到很多相关的代码,同时也有一些库帮助我们实现了它们之间的转换,比如python里的scipy包,其旋转相关的转换代码在scipy.spatial.transform里面,C++里则可以使用Eigen库来实现。而在pytorch里,则可以使用pytorch3d包。在pytorch3dtransforms模块里,实现了丰富的(甚至还有2019年CVPR上提出的6d表示方式的转换)旋转表示方式之间的转换,而且都是batch形式的,对于深度学习非常友好方便。
      pytorch3d.transforms模块实现的转换包括:

    axis_angle_to_matrix,            # 轴角转旋转矩阵
    axis_angle_to_quaternion,        # 轴角转四元数
    euler_angles_to_matrix,          # 欧拉角转旋转矩阵
    matrix_to_euler_angles,          # 旋转矩阵转欧拉角
    matrix_to_quaternion,            # 旋转矩阵转四元数
    matrix_to_rotation_6d,           # 旋转矩阵转6d表示
    quaternion_apply,                # 使用四元数对三维点进行旋转
    quaternion_invert,               # 求四元数旋转的逆
    quaternion_multiply,             # 组合两次四元数旋转为一个四元数
    quaternion_raw_multiply,         # 同上一样
    quaternion_to_axis_angle,        # 四元数转轴角
    quaternion_to_matrix,            # 四元数转旋转矩阵
    random_quaternions,              # 生成随机四元数
    random_rotation,                 # 生成随机旋转矩阵
    random_rotations,                # 同上一样
    rotation_6d_to_matrix,           # 6d表示转旋转矩阵
    standardize_quaternion,          # 将一个单位四元数转换为标准形式(实部≥0)
    

      下面贴上pytorch3d.transforms模块的内容,如果只是需要其中某个转换的代码的话,可以直接复制过去当成一个函数使用。

    # Copyright (c) Facebook, Inc. and its affiliates.
    # All rights reserved.
    #
    # This source code is licensed under the BSD-style license found in the
    # LICENSE file in the root directory of this source tree.
    
    import functools
    from typing import Optional
    
    import torch
    import torch.nn.functional as F
    
    from ..common.types import Device
    
    
    """
    The transformation matrices returned from the functions in this file assume
    the points on which the transformation will be applied are column vectors.
    i.e. the R matrix is structured as
    
        R = [
                [Rxx, Rxy, Rxz],
                [Ryx, Ryy, Ryz],
                [Rzx, Rzy, Rzz],
            ]  # (3, 3)
    
    This matrix can be applied to column vectors by post multiplication
    by the points e.g.
    
        points = [[0], [1], [2]]  # (3 x 1) xyz coordinates of a point
        transformed_points = R * points
    
    To apply the same matrix to points which are row vectors, the R matrix
    can be transposed and pre multiplied by the points:
    
    e.g.
        points = [[0, 1, 2]]  # (1 x 3) xyz coordinates of a point
        transformed_points = points * R.transpose(1, 0)
    """
    
    
    def quaternion_to_matrix(quaternions):
        """
        Convert rotations given as quaternions to rotation matrices.
    
        Args:
            quaternions: quaternions with real part first,
                as tensor of shape (..., 4).
    
        Returns:
            Rotation matrices as tensor of shape (..., 3, 3).
        """
        r, i, j, k = torch.unbind(quaternions, -1)
        two_s = 2.0 / (quaternions * quaternions).sum(-1)
    
        o = torch.stack(
            (
                1 - two_s * (j * j + k * k),
                two_s * (i * j - k * r),
                two_s * (i * k + j * r),
                two_s * (i * j + k * r),
                1 - two_s * (i * i + k * k),
                two_s * (j * k - i * r),
                two_s * (i * k - j * r),
                two_s * (j * k + i * r),
                1 - two_s * (i * i + j * j),
            ),
            -1,
        )
        return o.reshape(quaternions.shape[:-1] + (3, 3))
    
    
    def _copysign(a, b):
        """
        Return a tensor where each element has the absolute value taken from the,
        corresponding element of a, with sign taken from the corresponding
        element of b. This is like the standard copysign floating-point operation,
        but is not careful about negative 0 and NaN.
    
        Args:
            a: source tensor.
            b: tensor whose signs will be used, of the same shape as a.
    
        Returns:
            Tensor of the same shape as a with the signs of b.
        """
        signs_differ = (a < 0) != (b < 0)
        return torch.where(signs_differ, -a, a)
    
    
    def _sqrt_positive_part(x: torch.Tensor) -> torch.Tensor:
        """
        Returns torch.sqrt(torch.max(0, x))
        but with a zero subgradient where x is 0.
        """
        ret = torch.zeros_like(x)
        positive_mask = x > 0
        ret[positive_mask] = torch.sqrt(x[positive_mask])
        return ret
    
    
    def matrix_to_quaternion(matrix: torch.Tensor) -> torch.Tensor:
        """
        Convert rotations given as rotation matrices to quaternions.
    
        Args:
            matrix: Rotation matrices as tensor of shape (..., 3, 3).
    
        Returns:
            quaternions with real part first, as tensor of shape (..., 4).
        """
        if matrix.size(-1) != 3 or matrix.size(-2) != 3:
            raise ValueError(f"Invalid rotation matrix  shape f{matrix.shape}.")
    
        batch_dim = matrix.shape[:-2]
        m00, m01, m02, m10, m11, m12, m20, m21, m22 = torch.unbind(
            matrix.reshape(*batch_dim, 9), dim=-1
        )
    
        q_abs = _sqrt_positive_part(
            torch.stack(
                [
                    1.0 + m00 + m11 + m22,
                    1.0 + m00 - m11 - m22,
                    1.0 - m00 + m11 - m22,
                    1.0 - m00 - m11 + m22,
                ],
                dim=-1,
            )
        )
    
        # we produce the desired quaternion multiplied by each of r, i, j, k
        quat_by_rijk = torch.stack(
            [
                torch.stack([q_abs[..., 0] ** 2, m21 - m12, m02 - m20, m10 - m01], dim=-1),
                torch.stack([m21 - m12, q_abs[..., 1] ** 2, m10 + m01, m02 + m20], dim=-1),
                torch.stack([m02 - m20, m10 + m01, q_abs[..., 2] ** 2, m12 + m21], dim=-1),
                torch.stack([m10 - m01, m20 + m02, m21 + m12, q_abs[..., 3] ** 2], dim=-1),
            ],
            dim=-2,
        )
    
        # We floor here at 0.1 but the exact level is not important; if q_abs is small,
        # the candidate won't be picked.
        # pyre-ignore [16]: `torch.Tensor` has no attribute `new_tensor`.
        quat_candidates = quat_by_rijk / (2.0 * q_abs[..., None].max(q_abs.new_tensor(0.1)))
    
        # if not for numerical problems, quat_candidates[i] should be same (up to a sign),
        # forall i; we pick the best-conditioned one (with the largest denominator)
    
        return quat_candidates[
            F.one_hot(q_abs.argmax(dim=-1), num_classes=4) > 0.5, :  # pyre-ignore[16]
        ].reshape(*batch_dim, 4)
    
    
    def _axis_angle_rotation(axis: str, angle):
        """
        Return the rotation matrices for one of the rotations about an axis
        of which Euler angles describe, for each value of the angle given.
    
        Args:
            axis: Axis label "X" or "Y or "Z".
            angle: any shape tensor of Euler angles in radians
    
        Returns:
            Rotation matrices as tensor of shape (..., 3, 3).
        """
    
        cos = torch.cos(angle)
        sin = torch.sin(angle)
        one = torch.ones_like(angle)
        zero = torch.zeros_like(angle)
    
        if axis == "X":
            R_flat = (one, zero, zero, zero, cos, -sin, zero, sin, cos)
        if axis == "Y":
            R_flat = (cos, zero, sin, zero, one, zero, -sin, zero, cos)
        if axis == "Z":
            R_flat = (cos, -sin, zero, sin, cos, zero, zero, zero, one)
    
        return torch.stack(R_flat, -1).reshape(angle.shape + (3, 3))
    
    
    def euler_angles_to_matrix(euler_angles, convention: str):
        """
        Convert rotations given as Euler angles in radians to rotation matrices.
    
        Args:
            euler_angles: Euler angles in radians as tensor of shape (..., 3).
            convention: Convention string of three uppercase letters from
                {"X", "Y", and "Z"}.
    
        Returns:
            Rotation matrices as tensor of shape (..., 3, 3).
        """
        if euler_angles.dim() == 0 or euler_angles.shape[-1] != 3:
            raise ValueError("Invalid input euler angles.")
        if len(convention) != 3:
            raise ValueError("Convention must have 3 letters.")
        if convention[1] in (convention[0], convention[2]):
            raise ValueError(f"Invalid convention {convention}.")
        for letter in convention:
            if letter not in ("X", "Y", "Z"):
                raise ValueError(f"Invalid letter {letter} in convention string.")
        matrices = map(_axis_angle_rotation, convention, torch.unbind(euler_angles, -1))
        return functools.reduce(torch.matmul, matrices)
    
    
    def _angle_from_tan(
        axis: str, other_axis: str, data, horizontal: bool, tait_bryan: bool
    ):
        """
        Extract the first or third Euler angle from the two members of
        the matrix which are positive constant times its sine and cosine.
    
        Args:
            axis: Axis label "X" or "Y or "Z" for the angle we are finding.
            other_axis: Axis label "X" or "Y or "Z" for the middle axis in the
                convention.
            data: Rotation matrices as tensor of shape (..., 3, 3).
            horizontal: Whether we are looking for the angle for the third axis,
                which means the relevant entries are in the same row of the
                rotation matrix. If not, they are in the same column.
            tait_bryan: Whether the first and third axes in the convention differ.
    
        Returns:
            Euler Angles in radians for each matrix in data as a tensor
            of shape (...).
        """
    
        i1, i2 = {"X": (2, 1), "Y": (0, 2), "Z": (1, 0)}[axis]
        if horizontal:
            i2, i1 = i1, i2
        even = (axis + other_axis) in ["XY", "YZ", "ZX"]
        if horizontal == even:
            return torch.atan2(data[..., i1], data[..., i2])
        if tait_bryan:
            return torch.atan2(-data[..., i2], data[..., i1])
        return torch.atan2(data[..., i2], -data[..., i1])
    
    
    def _index_from_letter(letter: str):
        if letter == "X":
            return 0
        if letter == "Y":
            return 1
        if letter == "Z":
            return 2
    
    
    def matrix_to_euler_angles(matrix, convention: str):
        """
        Convert rotations given as rotation matrices to Euler angles in radians.
    
        Args:
            matrix: Rotation matrices as tensor of shape (..., 3, 3).
            convention: Convention string of three uppercase letters.
    
        Returns:
            Euler angles in radians as tensor of shape (..., 3).
        """
        if len(convention) != 3:
            raise ValueError("Convention must have 3 letters.")
        if convention[1] in (convention[0], convention[2]):
            raise ValueError(f"Invalid convention {convention}.")
        for letter in convention:
            if letter not in ("X", "Y", "Z"):
                raise ValueError(f"Invalid letter {letter} in convention string.")
        if matrix.size(-1) != 3 or matrix.size(-2) != 3:
            raise ValueError(f"Invalid rotation matrix  shape f{matrix.shape}.")
        i0 = _index_from_letter(convention[0])
        i2 = _index_from_letter(convention[2])
        tait_bryan = i0 != i2
        if tait_bryan:
            central_angle = torch.asin(
                matrix[..., i0, i2] * (-1.0 if i0 - i2 in [-1, 2] else 1.0)
            )
        else:
            central_angle = torch.acos(matrix[..., i0, i0])
    
        o = (
            _angle_from_tan(
                convention[0], convention[1], matrix[..., i2], False, tait_bryan
            ),
            central_angle,
            _angle_from_tan(
                convention[2], convention[1], matrix[..., i0, :], True, tait_bryan
            ),
        )
        return torch.stack(o, -1)
    
    
    def random_quaternions(
        n: int, dtype: Optional[torch.dtype] = None, device: Optional[Device] = None
    ):
        """
        Generate random quaternions representing rotations,
        i.e. versors with nonnegative real part.
    
        Args:
            n: Number of quaternions in a batch to return.
            dtype: Type to return.
            device: Desired device of returned tensor. Default:
                uses the current device for the default tensor type.
    
        Returns:
            Quaternions as tensor of shape (N, 4).
        """
        o = torch.randn((n, 4), dtype=dtype, device=device)
        s = (o * o).sum(1)
        o = o / _copysign(torch.sqrt(s), o[:, 0])[:, None]
        return o
    
    
    def random_rotations(
        n: int, dtype: Optional[torch.dtype] = None, device: Optional[Device] = None
    ):
        """
        Generate random rotations as 3x3 rotation matrices.
    
        Args:
            n: Number of rotation matrices in a batch to return.
            dtype: Type to return.
            device: Device of returned tensor. Default: if None,
                uses the current device for the default tensor type.
    
        Returns:
            Rotation matrices as tensor of shape (n, 3, 3).
        """
        quaternions = random_quaternions(n, dtype=dtype, device=device)
        return quaternion_to_matrix(quaternions)
    
    
    def random_rotation(
        dtype: Optional[torch.dtype] = None, device: Optional[Device] = None
    ):
        """
        Generate a single random 3x3 rotation matrix.
    
        Args:
            dtype: Type to return
            device: Device of returned tensor. Default: if None,
                uses the current device for the default tensor type
    
        Returns:
            Rotation matrix as tensor of shape (3, 3).
        """
        return random_rotations(1, dtype, device)[0]
    
    
    def standardize_quaternion(quaternions):
        """
        Convert a unit quaternion to a standard form: one in which the real
        part is non negative.
    
        Args:
            quaternions: Quaternions with real part first,
                as tensor of shape (..., 4).
    
        Returns:
            Standardized quaternions as tensor of shape (..., 4).
        """
        return torch.where(quaternions[..., 0:1] < 0, -quaternions, quaternions)
    
    
    def quaternion_raw_multiply(a, b):
        """
        Multiply two quaternions.
        Usual torch rules for broadcasting apply.
    
        Args:
            a: Quaternions as tensor of shape (..., 4), real part first.
            b: Quaternions as tensor of shape (..., 4), real part first.
    
        Returns:
            The product of a and b, a tensor of quaternions shape (..., 4).
        """
        aw, ax, ay, az = torch.unbind(a, -1)
        bw, bx, by, bz = torch.unbind(b, -1)
        ow = aw * bw - ax * bx - ay * by - az * bz
        ox = aw * bx + ax * bw + ay * bz - az * by
        oy = aw * by - ax * bz + ay * bw + az * bx
        oz = aw * bz + ax * by - ay * bx + az * bw
        return torch.stack((ow, ox, oy, oz), -1)
    
    
    def quaternion_multiply(a, b):
        """
        Multiply two quaternions representing rotations, returning the quaternion
        representing their composition, i.e. the versor with nonnegative real part.
        Usual torch rules for broadcasting apply.
    
        Args:
            a: Quaternions as tensor of shape (..., 4), real part first.
            b: Quaternions as tensor of shape (..., 4), real part first.
    
        Returns:
            The product of a and b, a tensor of quaternions of shape (..., 4).
        """
        ab = quaternion_raw_multiply(a, b)
        return standardize_quaternion(ab)
    
    
    def quaternion_invert(quaternion):
        """
        Given a quaternion representing rotation, get the quaternion representing
        its inverse.
    
        Args:
            quaternion: Quaternions as tensor of shape (..., 4), with real part
                first, which must be versors (unit quaternions).
    
        Returns:
            The inverse, a tensor of quaternions of shape (..., 4).
        """
    
        return quaternion * quaternion.new_tensor([1, -1, -1, -1])
    
    
    def quaternion_apply(quaternion, point):
        """
        Apply the rotation given by a quaternion to a 3D point.
        Usual torch rules for broadcasting apply.
    
        Args:
            quaternion: Tensor of quaternions, real part first, of shape (..., 4).
            point: Tensor of 3D points of shape (..., 3).
    
        Returns:
            Tensor of rotated points of shape (..., 3).
        """
        if point.size(-1) != 3:
            raise ValueError(f"Points are not in 3D, f{point.shape}.")
        real_parts = point.new_zeros(point.shape[:-1] + (1,))
        point_as_quaternion = torch.cat((real_parts, point), -1)
        out = quaternion_raw_multiply(
            quaternion_raw_multiply(quaternion, point_as_quaternion),
            quaternion_invert(quaternion),
        )
        return out[..., 1:]
    
    
    def axis_angle_to_matrix(axis_angle):
        """
        Convert rotations given as axis/angle to rotation matrices.
    
        Args:
            axis_angle: Rotations given as a vector in axis angle form,
                as a tensor of shape (..., 3), where the magnitude is
                the angle turned anticlockwise in radians around the
                vector's direction.
    
        Returns:
            Rotation matrices as tensor of shape (..., 3, 3).
        """
        return quaternion_to_matrix(axis_angle_to_quaternion(axis_angle))
    
    
    def matrix_to_axis_angle(matrix):
        """
        Convert rotations given as rotation matrices to axis/angle.
    
        Args:
            matrix: Rotation matrices as tensor of shape (..., 3, 3).
    
        Returns:
            Rotations given as a vector in axis angle form, as a tensor
                of shape (..., 3), where the magnitude is the angle
                turned anticlockwise in radians around the vector's
                direction.
        """
        return quaternion_to_axis_angle(matrix_to_quaternion(matrix))
    
    
    def axis_angle_to_quaternion(axis_angle):
        """
        Convert rotations given as axis/angle to quaternions.
    
        Args:
            axis_angle: Rotations given as a vector in axis angle form,
                as a tensor of shape (..., 3), where the magnitude is
                the angle turned anticlockwise in radians around the
                vector's direction.
    
        Returns:
            quaternions with real part first, as tensor of shape (..., 4).
        """
        angles = torch.norm(axis_angle, p=2, dim=-1, keepdim=True)
        half_angles = 0.5 * angles
        eps = 1e-6
        small_angles = angles.abs() < eps
        sin_half_angles_over_angles = torch.empty_like(angles)
        sin_half_angles_over_angles[~small_angles] = (
            torch.sin(half_angles[~small_angles]) / angles[~small_angles]
        )
        # for x small, sin(x/2) is about x/2 - (x/2)^3/6
        # so sin(x/2)/x is about 1/2 - (x*x)/48
        sin_half_angles_over_angles[small_angles] = (
            0.5 - (angles[small_angles] * angles[small_angles]) / 48
        )
        quaternions = torch.cat(
            [torch.cos(half_angles), axis_angle * sin_half_angles_over_angles], dim=-1
        )
        return quaternions
    
    
    def quaternion_to_axis_angle(quaternions):
        """
        Convert rotations given as quaternions to axis/angle.
    
        Args:
            quaternions: quaternions with real part first,
                as tensor of shape (..., 4).
    
        Returns:
            Rotations given as a vector in axis angle form, as a tensor
                of shape (..., 3), where the magnitude is the angle
                turned anticlockwise in radians around the vector's
                direction.
        """
        norms = torch.norm(quaternions[..., 1:], p=2, dim=-1, keepdim=True)
        half_angles = torch.atan2(norms, quaternions[..., :1])
        angles = 2 * half_angles
        eps = 1e-6
        small_angles = angles.abs() < eps
        sin_half_angles_over_angles = torch.empty_like(angles)
        sin_half_angles_over_angles[~small_angles] = (
            torch.sin(half_angles[~small_angles]) / angles[~small_angles]
        )
        # for x small, sin(x/2) is about x/2 - (x/2)^3/6
        # so sin(x/2)/x is about 1/2 - (x*x)/48
        sin_half_angles_over_angles[small_angles] = (
            0.5 - (angles[small_angles] * angles[small_angles]) / 48
        )
        return quaternions[..., 1:] / sin_half_angles_over_angles
    
    
    def rotation_6d_to_matrix(d6: torch.Tensor) -> torch.Tensor:
        """
        Converts 6D rotation representation by Zhou et al. [1] to rotation matrix
        using Gram--Schmidt orthogonalization per Section B of [1].
        Args:
            d6: 6D rotation representation, of size (*, 6)
    
        Returns:
            batch of rotation matrices of size (*, 3, 3)
    
        [1] Zhou, Y., Barnes, C., Lu, J., Yang, J., & Li, H.
        On the Continuity of Rotation Representations in Neural Networks.
        IEEE Conference on Computer Vision and Pattern Recognition, 2019.
        Retrieved from http://arxiv.org/abs/1812.07035
        """
    
        a1, a2 = d6[..., :3], d6[..., 3:]
        b1 = F.normalize(a1, dim=-1)
        b2 = a2 - (b1 * a2).sum(-1, keepdim=True) * b1
        b2 = F.normalize(b2, dim=-1)
        b3 = torch.cross(b1, b2, dim=-1)
        return torch.stack((b1, b2, b3), dim=-2)
    
    
    def matrix_to_rotation_6d(matrix: torch.Tensor) -> torch.Tensor:
        """
        Converts rotation matrices to 6D rotation representation by Zhou et al. [1]
        by dropping the last row. Note that 6D representation is not unique.
        Args:
            matrix: batch of rotation matrices of size (*, 3, 3)
    
        Returns:
            6D rotation representation, of size (*, 6)
    
        [1] Zhou, Y., Barnes, C., Lu, J., Yang, J., & Li, H.
        On the Continuity of Rotation Representations in Neural Networks.
        IEEE Conference on Computer Vision and Pattern Recognition, 2019.
        Retrieved from http://arxiv.org/abs/1812.07035
        """
        return matrix[..., :2, :].clone().reshape(*matrix.size()[:-2], 6)
    
    展开全文
  • 旋转矩阵的数学原理

    千次阅读 2019-10-02 00:05:57
    旋转矩阵涉及到的是一种组合设计:覆盖设计。而覆盖设计,填装设计,斯坦纳系,t-设计都是离散数学中组合优化问题。它们解决的是如何组合集合中的元素以达到某种特定的要求。 为了使读者更容易明白这些问题,下面先...

    一、从寇克曼女生问题讲起

    旋转矩阵涉及到的是一种组合设计:覆盖设计。而覆盖设计,填装设计,斯坦纳系,t-设计都是离散数学中组合优化问题。它们解决的是如何组合集合中的元素以达到某种特定的要求。

    为了使读者更容易明白这些问题,下面先从一道相当古老的数学名题讲起。

    (一)寇克曼女生问题

    某教员打算这样安排她班上的十五名女生散步:散步时三名女生为一组,共五组。问能否在一周内每日安排一次散步,使得每两名女生在这周内一道散步恰好一次?

    看起来题目似乎很简单,然而它的彻底解决并不容易。事实上,寇克曼于1847年提出了该问题,过了100多年后,对于一般形式的寇克曼问题的存在性才彻底解决。

    用1-15这15个数字分别代表这15个女生,下面给出一组符合要求的分组方法:

    星期日:(1,2,3),(4,8,12),(5,10,15),(6,11,13),(7,9,14)

    星期一:(1,4,5),(2,8,10),(3,13,14),(6,9,15),(7,11,12)

    星期二:(1,6,7),(2,9,11),(3,12,15),(4,10,14),(5,8,13)

    星期三:(1,8,9),(2,12,14),(3,5,6),(4,11,15),(7,10,13)

    星期四:(1,10,11),(2,13,15),(3,4,7),(5,9,12),(6,8,14)

    星期五:(1,12,13),(2,4,6),(3,9,10),(5,11,14),(7,8,15)

    星期六:(1,14,15),(2,5,7),(3,8,11),(4,9,13),(6,10,12)

    该问题就是最典型的组合设计问题。其本质就是如何将一个集合中的元素组合成一定的子集系以满足一定的要求。表面上看起来,寇克曼女生问题是纯粹的数学游戏,然而它的解却在医药试验设计上有很广泛的运用。

    寇克曼女生问题是t-设计中很特殊的一类——可分解斯坦纳设计。下面我会详细解释这几个名词的含义。

    (二)几种组合设计的含义

    所谓t-设计是“策略组态,Tactical Configuration”的简称。

    不妨用数学语言来定义t-设计

    \(S=\{S_1,S_2,...,S_v\}\)是一个包含有\(v\)个元素的集合;

    \(B_1,B_2,...,B_b\)\(S\)\(b\)个子集,而它们包含的元素个数和都是\(k\)个;

    \(B=\{B_1,B_2,...,B_b\}\)是由这\(b\)个子集组成的集合(子集系),

    对于固定整数t,和S的任意一个t元子集(t≥1),如果包含该子集的\(B\)中子集的个数都是同一个常数\(\lambda t\),则称\(B=\{B_1,B_2,...,B_b\}\)是集合\(S\)上的一个\(t-(v,k,\lambda t)\)设计,简称t-设计

    如果\(t-(v,k,\lambda t)\)设计中,\(t=2, \lambda=1\),则称为斯坦纳系(Steiner)。在该领域,我国已故的数学家陆家羲作出的巨大的贡献,如今每一本讲组合设计的书讲到这个问题,就不能不提到他的大名和以他的名字命名的定理。至今为止,斯坦纳系仍然存在着许多未解决的问题,至今还没有解决\((17,5,4=476)\)\(S(18,6,5=1428)\)的存在或不存在。虽然它的参数显得很小。

    而旋转矩阵涉及的则是另一种更加复杂、参数更多的组合设计——覆盖设计。

    覆盖设计是一种经过精心设计的b个区组组成的子集系,其中每个区组都有\(k\)个元素组成。它可以确保如果选出\(k\)个元素,有\(m\)个在其中,至少有\(\lambda\)个区组中的元素有\(t\)个元素符合。区组中元素的顺序与区组的排列顺序不影响覆盖设计本身。

    \((c:v,k,t,m,\lambda=b)\)

    可以用数学语言来定义比较简单的覆盖设计:

    \(S=\{S_1,S_2,...,S_v\}\)是一个包含有\(v\)个元素的集合;

    \(B_1,B_2,...,B_b\)\(S\)\(b\)个子集,而它们包含的元素个数都是\(k\)个;

    \(B=\{B_1,B_2,...,B_b\}\)是由这\(b\)个子集组成的集合(子集系)。

    对于固定整数\(t\),和\(S\)的任意一个\(t\)元子集(\(t\ge 1\)),如果该子集至少包含在\(B\)\(\lambda\)个区组中,则称\(B=\{B_1,B_2,...,B_b\}\)是集合\(S\)上的一个\(c-(v,k,\lambda t)\)设计,简称覆盖设计。

    填装设计是与覆盖设计相反的设计:

    \(S=\{S_1,S_2,...,S_v\}\)是一个包含有\(v\)个元素的集合;

    \(B_1,B_2,...,B_b\)\(S\)\(b\)个子集,而它们包含的元素个数都是\(k\)个;

    \(B=\{B_1,B_2,...,B_b\}\)是由这\(b\)个子集组成的集合(子集系)。

    对于固定整数\(t\),和S的任意一个\(t\)元子集(\(t \ge 1\)),如果该子集至多包含在\(B\)\(\lambda\)个区组中,则称\(B=\{B_1,B_2,...,B_b\}\)是集合\(S\)上的一个\(p-(v,k,\lambda t)\)设计,简称填装设计。

    t-设计又叫恰好覆盖与恰好填装。t-设计不一定存在,而覆盖设计一定存在。t-设计中,\(\lambda=1\),而覆盖设计一般\(\lambda>1\)。此外,t-设计\(m=t\),所以t-设计只是覆盖设计中比较特殊的一种。

    只要\(b\)足够大,显然覆盖设计一定存在。而有意义的是找到\(b\)的最小值,并找出在上最小值下的覆盖设计,此时的覆盖设计叫做最小覆盖。寻找最小覆盖的问题是组合优化问题的一类,被称为集合覆盖问题(SCP,Set covering problem),与著名的推销员旅行问题或成本最小化、利润最大化问题,都是优化问题的一种。

    但是集合覆盖问题往往经这些问题更加困难。因为其它问题往往已经有比较成熟的、固定的方法。而覆盖设计并没有通用的公式,所以大部分的设计即使用如克雷般超级电脑也很难求出,全盘搜索的算法耗用时间将会是一个天文数字。

    这方面,算法就显得相当重要。Oester Grad教授创造出一种全新的模拟算法,它大大提高了求解覆盖设计的速度,但它不能保证找到的覆盖设计一定是最小覆盖设计。它具有很强的通用性。而之前的其他算法往往只能解决固定某些参数的特定问题,解决的往往只是一类问题。

    对覆盖设计的研究始于19世纪,1835年J·Plue Cker和W.S.B.Wool House(1844)开始研究此类问题。

    到了1969年,人们发现它对军队中布阵与战略设计以及计算机芯片设计都大有用途,因此得到了迅速发展。在统计上,医药设计,农业试验,核研究,质量控制甚至在彩票中都大有用途。

    组合设计问题往往来自于智力游戏,对它们的研究也是纯数学的。但是当研究逐渐深入时,人们逐渐地在生产与其他学科中发现了它的用武之地。这样对它的研究就有了更强大的动力,吸引了更多人的注意,成果也就更加丰富。

    在选7的彩票涉及的旋转矩阵中,所有的(6,六)型和(5,五)型旋转矩阵都是t-设计。而一般的旋转矩阵都是覆盖设计。由于数学上对t-设计研究的比较多,所以有时候我们可以利用t-设计生成一些覆盖设计。

    如以下的设计即为一个\(t-(10,3,3)\)设计,它在有限射影几何中有很广泛的运用。

    B:(2,3,4),(1,5,10),(1,6,9)

    (1,7,8),(2,9,10),(3,8,10)

    (4,8,9),(4,6,7),(3,5,7),(2,5,6)

    即1-10每个数字都出现了3次,而且每两个数字恰好一起出现1次。从它可以生成10注10个号(7,六)型矩阵(它相当对称,平衡但不是最优的),具体生成方法很简单,取每一组的剩余的7个数就可以生成对应的一组。

    (三)组合设计的研究内容

    1.存在性问题

    若给出要求,研究符合要求的组合设计是否存在,以及存在的条件问题。比如,彩票中的覆盖设计问题,它的存在性就不是问题,因为只要注数足够多,总是可以覆盖的(它的上限为复式投注即完全组合,有意义的是它的下限)。而t-设计又叫恰好覆盖,它的存在性就是一个很值得研究的问题,也就是说,参数要符合什么条件,才会存在恰好覆盖一次的设计。

    对存在性的研究更多的是从理论上。然而,对于一般情形的t-设计是否存在的问题,还远远没有解决。

    2.构造问题

    如果已知某种组合设计存在时,如何把它们构造出来?这是与实际应用联系最紧的问题。实际上,最终无论在彩票中,还是新药设计中,人们关心的是构造出的组合设计。经过数学家上百年的努力,现在已经有一些构造方法。如利用有限的射影几何,关联矩阵,数论中的差集等构造出大量的设计。用组合论自身也能解决一些构造问题。然而,对于一般情形的组合设计的构造性问题离解决还相当遥远。比如彩票中覆盖设计问题(即旋转矩阵)当参数变大时,设计的难度是几何级数上升。

    对于一般的最小覆盖问题,仍然没有通用的构造方法。也就是说,目前市场上出现的许许多号码比较多的旋转矩阵,都很难保证是最小覆盖设计,也就是无法保证它是最优的。很多旋转矩阵不断地有人刷新它的下限纪录,也就是越来越接近最小覆盖设计。然而,要证明一个旋转矩阵是否已经是最小覆盖设计,是极其困难的,如果号码很少,还可以通过计算机编程用穷举的方式来解决,而号码稍微多一点,用穷举法超级电脑运算所耗用的时间也将是天文数字。

    3.组合设计之间的关系

    例如:一个组合设计是否与另外一个组合设计本质上一样的(同构)。比如把组合中某两个数字互换,这两个设计应该算同一种设计。每一种设计的同构设计是非常多的。有些同构是很难直接看出的,所以就需要研究同构的设计有什么特点,如何准确快速的判断和产生同构设计。

    组合设计还研究如何由一个组合设计构造出另外一个。

    比如旋转矩阵中存在着这样的问题,比如10个号码01-10,开始我先选定3注:

    01,02,03,04,05,06,07,

    01,03,05,07,08,09,10,

    02,04,06,07,08,09,10

    问如何添上尽可能少的注数,使它成为(7,六)型平衡式矩阵。

    又如一个旋转矩阵与另外一个旋转矩阵是否同构。即使两个旋转矩阵所有参数都相同,也不一定同构。然而,在实际运用中,人们并不关心同构问题。因为只要能用就行了。

    又如10个号码(7,六)型的有8注,比如是01-10这10个号码,问能否在这基础上添上尽可能少的注数,使得它成为11个号码的(7,六)型的旋转矩阵(01-11)。

    4.计数问题

    如果已知某类组合设计存在,自然希望知道这类设计的个数。也就是说互不同构的设计的个数。然而,这个问题是一个极其艰难的问题,现在还很少人去研究它。

    比如很简单的10个号码的(7,六)型矩阵,共有多少种。号码一多,这将是一个很困难的问题。

    5.最优设计

    在诸多的满足要求的组合设计中,找到一个最优的设计,这是它研究的内容。比如覆盖设计很多,如何找出最小覆盖设计就是一具艰难的问题。旋转矩阵中需要用到组合优化的算法与组合构造算法。

    二、旋转矩阵的主要算法

    (一)对旋转矩阵做出突出贡献的主要数学家

    旋转矩阵是一个看似简单实际却异常复杂的问题,尽管有许许多多的人对它非常感兴趣,然而真正在这个领域内做出了开创性贡献的人却不是很多。要想在此领域有所作为,不仅要对组合设计的经典理论和常用方法有深入的了解,还要在此基础上有所创新。有许多国外的彩票专家,比如美国的盖尔·霍华德女士,声称旋转矩阵是由她首先提出来的。实际上,所有的旋转矩阵都是组合数学家们经过多年的精心研究得出的,而不是霍华德这样的彩票专家所能研究出来的。

    在此领域内做出了突出贡献的主要组合数学家有以下几位:

    1.Patric Osergard

    他的主要贡献是用了全新的模拟冷却算法解决了旋转矩阵的构造问题,运用他的模拟冷却程序,可以很迅速的产生许许多多的旋转矩阵。

    2.Alex Sidorenko

    他研究出了许多旋转矩阵和几种产生旋转矩阵的基于图灵理论的一般方法。

    3.Greg Kuperberg

    他注意到线性的[v,t]编码的补集可以给出区组长度不定的覆盖设计,而这可以产生对现有的旋转矩阵的一系列改进。

    4.Dan Gordon

    他不仅是覆盖设计领域内多篇经典论文的合作者,而且总结了所有的旋转矩阵的成果,并且时时关注着该领域的最新进展。他收集的旋转矩阵是迄今为止最全面、最权威的。而这一切全凭他个人的兴趣,没有任何经费的支持。

    以下我将对以上的数学家作一些介绍:

    Dan Gordon是圣地亚哥的通讯研究中心的研究员。个人兴趣:计数理论、组合学、代数分析。

    Greg Kuperberg是美国加州大学的数学系的副教授。主要研究方向是复杂性分析和微积分。他在覆盖设计的主要论文有:

    (1)Asymptotically optimal covering designs(with Daniel Gordon,Oren Patashnik,and Joel Spe-ncer).J.Combin.Theory Ser.A 75(1996),page 270-280.

    (2)Highly saturated packings and reduced coverings(with Gabor Fejes Toth and Wlodzimierz Kuperberg).Monats.Math.125(1998),page 127-145.

    Patric Ostergard是芬兰赫尔辛基理工大学计算科学和工程系的教授。

    他的兴趣集中在数学和计算机科学中系列问题。他的主要研究方向可分为以下几类:

    (1)组合结构的设计

    编码(覆盖编码,纠错编码等等)

    组合设计

    几何填装和覆盖问题

    (2)局部搜索的优化

    模拟冷却算法

    禁忌搜索

    全局优化的随机方法

    (3)加密与解密

    他是1996寇克曼奖的得主,这个奖是由国际组合学协会颁发的,以已故的著名组合学家寇克曼的名字来命名,用来奖励对组合学有突出贡献的数学家。除此之外,他还是组合设计杂志的编辑。

    他在覆盖设计的主要贡献是彩了模拟冷却方法研究出了全新的构造覆盖设计的全新方法。

    他在此领域的主要论文:

    Nurmela,K.J.and Ostergard,P.R.J.“Constructing Covering Designs by Simulated Annealing”,Helsinki University of Technology Digital Systems Laboratory Series B:Technical Reports,No.10;January,1993.

    (二)旋转矩阵的主要算法

    旋转矩阵的定义是很容易明白的,一般的业余数学爱好者理解没有任何障碍。然而,如何快速有效的构造旋转矩阵是一个数学家们一直在研究的问题。当然,这其中最关键的就是算法。而近年来最好的算法无疑是模拟冷却算法,它主要是由Patric Ostergard首创,并且得到了许多后来者的发展。

    下面我简要介绍一下他论文所用的算法的主要思想。

    1.Simulated Annealing模拟冷却算法

    模拟冷却算法是一种随机搜索方法,它的主要特点是不用穷遍集合中每一种可能性就可以找到最优或几乎最优的状态。它是通过模拟一个分子系统的自然冷却过程来做到这一点的。在每一种状态,它随机地选择了一种相邻的状态,如这种相邻的状态有一个更低的成本,系统将会转移到该状态。如果这种相邻的状态有一个更高的成本,系统将可能会转移到该状态,也可能不会转移到该状态。转移的概率依赖现在的状态的温度参数(该值越高,转移的概率越大)和两个状态之间的成本的差异(差异越大,转移的概率越大)。温度将会渐渐低下来,最终会达到均衡。模拟冷却算法常常用来尝试发现离散数学中一些问题的几乎最优的解。

    模拟冷却算法的一般步骤如下:

    1. 给定一个初始状态和初始温度
    2. 外部循环
      1. 内部循环
        1. 随机选择一个相邻状态。若相邻状态的成本更低,转移
        2. 若相邻状态的成本更高,转移的概率为exp{-成本差异/温度}
      2. 降低温度
    3. 返回所遇到的最优状态

    模拟冷却算法的设计者需要选择以下6个参数:

    • 初始温度和初始状态
    • 一种状态的成本函数
    • 一种状态的相邻函数
    • 冷却程序
    • 内部循环方法
    • 外部循环方法

    初始状态和初始温度实际上对算法影响不大,成本函数一般来说也比较容易定义,尤其是对覆盖设计来说,成本可以定义成重复数字的总个数。相邻函数也可以随机挑选一个向量来解决。而有效的冷却程序一般用T'=rT,这里T指原来的温度,T'是新的温度,r是常数,也叫冷却因子。

    Patric Ostergard的关于覆盖设计的经典论文基本上就是如此定义模拟算法的参数的。

    运用该算法,可以很容易算出一般的旋转矩阵。

    除了模拟冷却算法之外,还有另外一些构造旋转矩阵的常用方法。

    2.非连通的集合来结合覆盖设计

    如果对某个\(v=v_1+v_2\)和所有的\(t_1+t_2=t\),都有大小为\(N_1\)的覆盖设计\((v_1,k_1,t_1)\)和大小为\(N_2\)的覆盖设计\((v_2,k_2,t_2)\)存在,那么将有大小为\(N=N_1*N_2\)的覆盖设计存在。然而,可以用这种方法产生的旋转矩阵数量很少,而且构造的过程也很复杂。很少的旋转矩阵是用这种方法产生的。

    3.贪婪算法

    这种算法产生了许多许多的旋转矩阵,这种算法的核心思想是:每个区组都尽可能少重复前面区组的数字,一直重复下去,直到你得到一个覆盖设计。你可以用顺序、逆序或灰色、随机的顺序来重复这个过程。或者可以用你所喜欢的其它顺序。事实上,笔者起初的时候正是用这个方法来产生一些比较简单的矩阵,但是这种算法看起来容易,实际上却十分繁琐,如果不用计算机,即使是很简单的矩阵,也要耗费无数的精力。而且,这种算法只能保证可以产生旋转矩阵,却无法保证产生的旋转矩阵一定是最优的。当参数很大时,用它产生的矩阵离最优的矩阵还差的很远。

    但是,可以用这种方法产生旋转矩阵,然后利用其他的优化算法对它再进一步优化,这样可以产生比较优良的旋转矩阵。

    4.诱致算法

    Greg Kuperberg是这种算法的主要创立者和提倡者。

    先利用一个巨大的参数为\((V,K,t)\)的旋转矩阵,从\(V\)个点中按照某种顺序或完全随机的选出\(v\)个点,然后将用他们原来的长度为K的区组隔断,得到了每个区组个数不定的一个覆盖。最后,将这个覆盖进行如下的修补即可:对每一个长度为1的区组,将该区组替换成一个\((1,k,t)\)的覆盖设计。这是一种比较复杂的算法,然而,却是迄今最好的算法之一。

    运用它可以产生优化程度比较高的矩阵。然而,运用这种算法的一个很大的限制是,必须要有一个参数很大的旋转矩阵和许许多多的参数比它小的矩阵。

    这样的条件比较苛刻,所以它的运用不是十分广泛。

    三、旋转矩阵如何提高中奖概率

    (一)对彩票中一些常用的概率的理解

    彩票之所以能够吸引成千上万的彩民,是因为它给许多人提供了一夜暴富的机会。而且成本很小,两元钱的投入就可能带来数以百万的回报。然而,它的每百万的大奖都是由上百万张没有中奖的彩票提供的。中大奖的机会微乎其微,这正是彩票的本质与魅力根源。

    尽管彩票中大奖的机会都十分小,然而各种不同彩票中大奖的概率还是有大差别的,比如32选7的北京风采中一等奖的概率为:

    而36选7的北京体彩中一等奖的概率仅为:

    简而言之,彩票号码球的个数越少,中大奖的机会就越大,而同样数目的号码球,选5型的比选7型的中大奖的机会更大。当然,天下不会有免费的午餐,中大奖的机会越大,大奖累种的奖金额越小。越难中的大奖,它累积的大奖金额越高。

    一些人买彩票的心态可以分为两种,一种是瞄准百万大奖的,梦想着有朝一日可以改变自己的生活。这种彩民一般并不在意中奖的机会有多渺茫,而比较关心本期累积的大奖奖金有多高。另一种则更希望经常在彩市中有所斩获得一些不大不小的奖。这种人比较适合玩中奖机会稍大的彩票,如上海的天天彩与体彩的“四花选四”。这两种彩票中奖机会都比较高。

    除了考虑奖金最高的一等奖,还要考虑次等级的二等奖、三等奖,它们的奖金虽不如一等奖那么高,也颇为丰厚。那么,如果买单注的话,中这些级别的奖的概率都可以算出来。但是,如果买多注的话,如何计算中各个级别的概率呢?如果这些注数胡乱买的,相互之间没有经过科学组合,那么应该如何估计中奖概率呢?

    旋转矩阵利用科学的组合方式提高了这种中奖概率。


    (二)组合投注的中奖概率分析

    不妨举个例子,以北京体彩(36选7)为例进行分析。比如我买了一注,那么中二等奖(6个正选号,无特别号)的概率可以用如下方法计算:

    现在我想买10个号,如果用复式投注的话,需要买120注,中二等奖的概率可以这样计算:因为只有选的10个号中至少有6个正选号才可能出现二等奖。复式投注实际上由于各注之间重复的号码太多可以效率很低。

    如果用旋转矩阵来投注的话,10个号码,需要购买10注。旋转矩阵的含义是,只要选的10个号中包含了7个正选号,一定含有二等奖。如果10个号中包含了6个正选号,也有可能有二等奖。

    各种方式具体中奖概率比较:

    1.假设用复式投注

    若10个号码中含全部正选号;此条件下中二等奖的概率为100%,此条件发生概率为:

    若10个号中含6个正选号,此条件下中二等奖的概率为100%,此条件发生概率为:

    若10个号中含的正号少于6个,此条件下中二等奖的概率为0;

    根据全概率公式,复式投注中二等奖概率为:

    也就是说,旋转矩阵用的注数为复式投注的1/12,而中二等奖的概率却是复式投注的1/3,从数字的利用效率来说,旋转矩阵组合号码的效率大约是复式投注的4倍。从成本收益对比来看,旋转矩阵具有明显的优势。

    (三)旋转矩阵中奖的上下限分析

    可以分析如果你选中某些正选号的情况下中奖情况(具体见平衡式旋转矩阵那一章)。

    旋转矩阵为什么可以提高中奖的概率呢?

    实际上,笔者一再强调天下没有免费的午餐。旋转矩阵之所以可以比复式投注与一般组合平均每注中二等奖概率要高,是以它牺牲中二等奖时的注数为前提的。举个例子,就可以明白这个问题了。

    比如你准备选01,02,03,04,05,06,07,08,09,10这10个号,买10注有几种组合方式,如:

    (1)最极端的方式:01,02,03,04,05,06,07,同样的号连买10注。

    (2)一般组合方式:如轮次矩阵,10注。表2-1:10个号码组合轮次矩阵

    1. 01 02 03 04 05 06 07
    2. 02 03 04 05 06 07 08
    3. 03 04 05 06 07 08 09
    4. 04 05 06 07 08 09 10
    5. 05 06 07 08 09 10 01
    6. 06 07 08 09 10 01 02
    7. 07 08 09 10 01 02 03
    8. 08 09 10 01 02 03 04
    9. 09 10 01 02 03 04 05
    10. 10 01 02 03 04 05 06

    (3)旋转矩阵方式:10注。表2-2:10个号码(7,六)型旋转矩阵(10注)

    1. 01 05 06 07 08 09 10
    2. 02 03 04 06 07 08 09
    3. 02 03 04 05 07 08 10
    4. 02 03 04 05 06 09 10
    5. 01 03 04 05 06 07 08
    6. 01 02 04 05 06 07 09
    7. 01 02 03 05 06 07 10
    8. 01 02 03 05 08 09 10
    9. 01 02 04 06 08 09 10
    10. 01 03 04 07 08 09 10

    假设开出的中奖号码正选号即为01,02,03,04,05,06,07,那么显然第一种方法中的奖最多,为10个特等奖;系二种方法次之,有1个特等奖和2个二等奖;旋转矩阵只有3个二等奖。

    假设开出的中奖号码为02,03,04,05,06,07,08

    第一种方法有10个二等奖;

    第二种方法有1个特等奖和2个二等奖;

    第三种方法有2个二等奖。

    假设开出的号码为01,02,04,05,07,08,09

    第一种方法:没有二等奖;

    第二种方法:没有二等奖;

    第三种方法:1个二等奖。

    从以上可以看出,旋转矩阵的优势在于它比较均匀,只要选对7个正选号,无论是哪7个,它都能保证二等以上的奖。而一般组合方法,在某种情况下,可能会得一堆二等奖,其它情况可能只会得很小的奖。

    当然,从以上例子可以看出,在某些情况下,旋转矩阵得的奖比一般组合方式要少。这是由于旋转矩阵尽量分散了风险所致。实际上,旋转矩阵正如保险一样,它使你的收益更加确定。当然,在某些情况下,你收益可能会减少,但总体来说,稳定的收益更符合人的天性。

    此外,旋转矩阵中蕴含的修习的原理在整个保险的保费确定和保险经营中都得到了应用。实际上,旋转矩阵中蕴含的深刻的组合设计理论近百年来一直吸引着数学家,至今仍有许多这方面未解的难题。旋转矩阵以其无穷的魅力,激励着一代代数学家去突破、解决。

    转载于:https://www.cnblogs.com/hhh5460/p/5677733.html

    展开全文
  • 四元数 (Quaternion)3.1 四元数的运算3.2 四元数默认在极坐标下3.3 四元数的常用插值方法3.4 贝塞尔曲线和 Squad 插值4 欧拉角、旋转矩阵、四元数的互相转换4.1 欧拉角和旋转矩阵4.2 四元数和旋转矩阵4.3 欧拉角和...
  • 双色球彩票生成

    2018-06-25 13:32:56
    c++双色球彩票生成器,可以生成红球篮球,循环输出。。
  • 计算两个对应点集之间的旋转矩阵R和转移矩阵T

    万次阅读 多人点赞 2018-05-05 23:06:35
    为了寻找这两个点集之间的旋转矩阵RRR和转移矩阵ttt。可以将这个问题建模成如下的公式: B=R∗A+tB=R∗A+tB = R*A+t 为了寻找RRR和ttt,通常需要一下三个步骤: 计算点集合的中心点 将点集合移动到原点,...
  • 可以生成M*N矩阵,用于一些特定需要标记的用途,使用MFC编写,界面操作简单直观,可以拖动显示矩阵
  • %这里随机生成两个四元数 q1 和 q2 作为测试数据 q1 = rand(1,4); q2 = rand(1,4); q1 = q1./norm(q1); q2 = q2./norm(q2); %计算俩个四元数之间的旋转角度差值(单位度数) 2*acos(abs(sum(q1.*q2)))*180/pi %...
  • 使用numpy实现矩阵的翻转(flip)与旋转

    千次阅读 2022-03-26 15:55:17
    如果是一维数组,可以使用a = a[:,:,-1]实现反转! numpy.flip(m, axis=None) Reverse the order of elements in an array along the given axis. The shape of the array ...随机生成一个二维数组 1 2.
  • http://simulations.narod.ru/ 生成正则 simpex,函数 simplex.m 单纯形是多维四面体 也random_rotation.m-生成随机旋转矩阵 运行 zz_plot_simplex_in_2d.m 生成单纯形,随机旋转,投影到 2d,在 2d 中绘制单纯形
  • C#实现N*N的顺时针旋转矩阵

    千次阅读 2013-06-19 20:38:15
    看面试算法题突然想到大学时写的一算法题 using System; using System.Collections.Generic;...using System.Linq;...//输出按n时的矩阵 namespace n时矩阵 ...随机生成边长在1~20中的矩阵
  •   记直接线性变换得到的旋转矩阵为RRR,平移向量为ttt,投影方式为, R←(RRT)−12R R \leftarrow (RR^T)^{-\frac{1}{2}}R R←(RRT)−21​R 然后代入线性方程组,利用最小二乘法更新平移向量ttt。   验证: 随机...
  • 随机生成某些稀疏矩阵

    千次阅读 2016-01-11 11:27:39
    1.单位稀疏矩阵 函数 speye 格式 S = speye(m,n) %生成m×n的单位稀疏矩阵 ...格式 R = sprand(S) %生成与S具有相同稀疏结构的均匀分布随机矩阵 R = sprand(m,n,density) %生成一个m×n的服
  • import numpy as np from numpy.linalg import cholesky import matplotlib...R2变换相当于先沿坐标轴伸缩后旋转翻转,va是伸缩比,酉阵vc是由斜转正的右乘旋转矩阵,其列向量为主轴方向(指上图椭圆最长轴和最短轴)。
  • 但是这些方法产生的随机数质量往往不是很高,而今天介绍的梅森旋转算法可以产生高质量的伪随机数,并且效率高效,弥补了传统伪随机生成器的不足。梅森旋转算法的最长周期取自一个梅森素数,由此命名为梅森旋转算法...
  • matlab 稀疏矩阵生成

    千次阅读 2019-02-27 15:01:48
    1.单位稀疏矩阵 函数 speye 格式 S = speye(m,n) %生成m×n的单位稀疏矩阵 ...格式 R = sprand(S) %生成与S具有相同稀疏结构的均匀分布随机矩阵 R = sprand(m,n,density) %生成一个m×n的服...
  • [转载]matlab 生成矩阵

    2021-04-21 07:58:21
    具体方法如下:将矩阵的元素用方括号括起来,按矩阵行的顺序输入各元素,同一行的各元素之间用空格或逗号分隔,不同行的元素之间用分号分隔。+2.利用M文件建立矩阵++++对于比较大且比较复杂的矩阵,可以为它专门...
  • python三维矩阵

    千次阅读 2020-11-29 13:25:49
    np.ones((256, 256, 75))merged = np.zeros(hr_shape)merged += hr_xmerged +=hr_ymerged += hr_zprint(merged)#部分输出]]上面是我敲的简单代码,我想建立一个三维的图片矩阵,具体操作跟代...
  • 然后利用改进的随机抽样一致算法去除误匹配点对,获得最优单应性矩阵;最后,对重叠区域像素进行动态融合处理。实验结果表明,采用本文算法可以有效地消除视频重叠区域的拼接缝和鬼影,同时可保证视频拼接系统的实时性。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 22,828
精华内容 9,131
关键字:

旋转矩阵 随机生成