１０次元立方体を眺める

pythonで１０次元立方体を２次元に射影して回転アニメーションを生成します．

結果

10D cube projected into 2D pic.twitter.com/T4sUQgqm7b

— S.Tabayashi (@tabayashi0117) 2020年4月5日

コード

（矢印を押すと展開します）

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.collections as mc
# from mpl_toolkits.mplot3d import Axes3D 
import itertools
import time

def generate_rotation_2d(dim=3, axis1=0, axis2=1):
    """回転行列を生成する
    
    Args:
        dim(int): 次元
        axis1(int): 回転面の軸のインデックス
        axis2(int): 回転面の軸のインデックス
        
    Returns:
        function: 回転行列
    """
    def rotation_2d(theta):
        """回転行列
        
        Args:
            theta(float): 回転角
            
        Returns:
            np.2darray: 回転行列
        """
        rotation_matrix = np.zeros((dim, dim))
        for i in range(dim):
            rotation_matrix[i,i] = 1
        rotation_matrix[axis1, axis1] = np.cos(theta)
        rotation_matrix[axis1, axis2] = -np.sin(theta)
        rotation_matrix[axis2, axis1] = np.sin(theta)
        rotation_matrix[axis2, axis2] = np.cos(theta)
        return rotation_matrix
    return rotation_2d


def prod_rotation(theta, rotation_2d_list, dim):
    """回転行列の積をとる"""
    # 単位行列を生成
    rotation_matrix = np.identity(dim)
    # すべての回転行列の積をとる
    for rotation in rotation_2d_list:
        rotation_matrix = rotation_matrix @ rotation(theta)
    return rotation_matrix

def plot_cube(dim=3):
    """立方体を図示する"""
    # 回転行列の全体
    rotation_2d_list = [generate_rotation_2d(dim, axis1, axis2) for axis1, axis2 in itertools.combinations(range(dim),2)]
    
    # 頂点全体を生成
    vertex_list = [np.array(list(coordinate)).reshape(-1,1) for coordinate in itertools.product(*[[-1/2,1/2] for _ in range(dim)])]

    # エッジのインデックス全体を生成
    edge_index_list = []
    for comb in itertools.combinations(range(len(vertex_list)),2):
        if np.sum(np.abs(np.array(vertex_list[comb[0]]) - np.array(vertex_list[comb[1]]))) == 1:
            edge_index_list.append(comb)

    # 速度（設定用）
    speed = 4/dim
    # 回転角
    theta = (np.pi /180) * speed 
    
    # 回転行列を生成
    rotation_matrix = prod_rotation(theta, rotation_2d_list, dim)
    
    # 回転角で何回回転したかカウント
    cnt = -1

    def update(frame):
        plt.cla()

        nonlocal vertex_list, edge_index_list, rotation_matrix, dim, theta, cnt

        if cnt > 0:
            vertex_list = [rotation_matrix @ v for v in vertex_list]

        x = np.concatenate([v[0] for v in vertex_list])
        y = np.concatenate([v[1] for v in vertex_list])
#         z = np.concatenate([v[2] for v in vertex_list_rotated])
        
        lim = 1. + 0.1 * dim
        
        ax.set_xlim([-lim, lim])
        ax.set_xlabel("X")
        ax.set_ylim([-lim, lim])
        ax.set_ylabel("Y")
        ax.set_axis_off()

        ax.scatter(x,y, s=40/dim, c="red", alpha=0.6)
        
        # 線のリスト．　[(x0, y0), (x1, y1)]が1つの線
        lines = [[vertex_list[i1][:2].flatten(), vertex_list[i2][:2].flatten()] for i1, i2 in edge_index_list]
        colors = ["blue" for i in range(len(edge_index_list))]
        
        plt.text( 
            0.96, 0.96, f"$\Theta$: {theta*cnt*(180/np.pi)%360.0:.2f}[deg]",
            ha='right', va='top',
            transform=ax.transAxes
        )
        
        # 複数の線を追加
        lc = mc.LineCollection(lines, colors=colors, linewidths=1.0/dim)
        ax.add_collection(lc)
        
        cnt += 1

    fig, ax = plt.subplots(figsize=(5,5))
    ani = animation.FuncAnimation(fig, update, interval = 20, frames = int((np.pi*4)/theta)+1)
    ani.save(f'cube/cube_{dim}d_to_2d_720.mp4', writer="ffmpeg", dpi=145) 

if __name__ == "__main__":
    for dim in range(2,11):
        start = time.time()
        plot_cube(dim=dim)
        elapsed_time = time.time() - start
        print(f"dim: {dim}, elapsed_time:{elapsed_time}[sec]")

数学的なうんちく

２次元の回転行列は

\begin{pmatrix}
\cos\theta & -\sin\theta \\
\sin\theta & \cos\theta
\end{pmatrix}

です．

3次元の回転は

\begin{pmatrix}
\cos\theta & -\sin\theta & 0 \\
\sin\theta & \cos\theta & 0 \\
0 & 0 & 1
\end{pmatrix} , \qquad 

\begin{pmatrix}
\cos\theta & 0 & -\sin\theta \\
0 & 1 & 0 \\
\sin\theta & 0 & \cos\theta 
\end{pmatrix} , \qquad 
\begin{pmatrix}
1 & 0 & 0 \\
0 & \cos\theta & -\sin\theta \\
0 & \sin\theta & \cos\theta
\end{pmatrix} 

の組み合わせで表現できます．

なので，自然な拡張として，N次元の回転を以下の行列の組み合わせで表現することができます：

\begin{align}
& \qquad \qquad \small{\text{i-th}} \qquad \qquad \small{\text{j-th}} \\
& \begin{pmatrix}
1 & \cdots & 0 & \cdots & 0 & \cdots &  0 \\
\vdots & \ddots & \vdots & & \vdots && \vdots \\
 0 & \cdots & \cos\theta & \cdots & -\sin\theta & \cdots & 0 \\
\vdots & & \vdots & & \vdots && \vdots \\
 0 & \cdots & \sin\theta & \cdots & \cos\theta & \cdots & 0 \\
\vdots & & \vdots & & \vdots &\ddots & \vdots \\
0 & \cdots & 0 & \cdots & 0 & \cdots &  1
\end{pmatrix} \begin{matrix} \\ \\  \small{\text{i-th}} \\ \\ \small{\text{j-th}}  \\ \\ \\ \end{matrix} ,  \qquad 1\le i < j \le N .
\end{align}

これを実装したのが上のコードです．

課題

１４次元までは上のコードで30分程度の計算時間で済んだのですが，それ以上次元を上げようとすると処理が終わらなくなりました．

処理が重くなりすぎたのは，回転行列による積の演算ではなく，頂点と辺のプロットでした．

N次元立方体の頂点の数は$2^N$，辺の数は$2^{N-1}N$です．

頂点の数は，10次元で1024, 20次元で1048576, 30次元で1073741824.

辺の数は，10次元で5120, 20次元で10485760, 30次元で16106127360.

恐ろしい増えっぷりですね．

その他の次元（観賞用）

8D pic.twitter.com/oEeI76nMu6

— S.Tabayashi (@tabayashi0117) April 5, 2020

6D pic.twitter.com/IMwKWSfwnz

— S.Tabayashi (@tabayashi0117) April 5, 2020

4D pic.twitter.com/k5HHkLK9Ak

— S.Tabayashi (@tabayashi0117) April 5, 2020

12D pic.twitter.com/0u2YvMkxpI

— S.Tabayashi (@tabayashi0117) April 5, 2020

14D pic.twitter.com/aNgFnMZfXY

— S.Tabayashi (@tabayashi0117) April 5, 2020

無限次元へ

本記事で行ったプロットから察するに，無限次元立方体も２次元に射影すると，ほとんどいたる角度から見たとき**「だいたい円盤」**っぽく見えるんだろうとは思うのですが，実際のところはどうなんでしょうね？

高次元は不思議なことが多いのですが，それを遥かに超えた無限次元はまったく想像がつきません．