is_q_matrix.py

from itertools import product
import numpy as np
from numpy import linalg as la
from matplotlib import pyplot as plt
import matplotlib.patches as patches
import pulp

def is_q_matrix(M):
  
  # 問題のサイズを取得
  n = M.shape[0]

  # 所与の行列と単位行列を連結
  I = np.identity(n)
  A = np.array([I,-M])

  # 相補錐の半空間表示を取得
  H = np.zeros((2**n, n, n))
  for k, pi in enumerate(product([0,1], repeat=n)):
    # 相補錐のベクトル表示を取得
    A_pi = np.concatenate([ A[pi[i],:,i].reshape(-1,1) for i in range(n) ], axis=1)
    # 相補錐の半空間表示を取得
    A_pi_inv = la.inv(A_pi)
    # 行列を格納
    H[k] = A_pi_inv

  # 分割後の凸多面体に対して空集合であるかを判定
  result = True
  for sigma in product(range(n), repeat=2**n):
    # 分割後の凸多面体を取得
    H_sigma = np.concatenate([ H[k,sigma[k]].reshape(1,-1) for k in range(2**n) ], axis=0)
    # 空集合であるかを判定
    result &= is_empty(H_sigma)

  return result

def is_empty(H_sigma):
  # LPライブラリの許容解探索を使って凸多面体が空集合であるかを判定

  # 問題のサイズを取得
  m, n = H_sigma.shape

  # 問題インスタンスを生成
  lp = pulp.LpProblem()

  # 変数を生成
  q = [ pulp.LpVariable('q_'+str(i)) for i in range(n) ]

  # 目的関数を設定
  # 実行可能性のみを知りたいのでゼロベクトルを設定
  c = [0]*n
  lp += pulp.lpDot(c, q)

  # 不等式制約を設定
  # 等号は付かないため十分に小さい値epsを設定
  # epsを使いたくなければファルカスの補題を使って等号付きの制約に変形可能
  # 今回は記事と同様の形式を使用
  eps = 10**(-5)
  for i in range(m):
    lp += pulp.lpDot(H_sigma[i], q) <= -eps

  # 求解
  status = lp.solve(pulp.PULP_CBC_CMD(msg=False))

  return status != 1

def show_complementary_cone(M):

  # 行列のサイズを取得
  n = M.shape[0]

  # 各列ベクトルを単位ベクトルに変換
  for i in range(n):
    M[:,i] /= la.norm(M[:,i])

  # 所与の行列と単位行列を連結
  I = np.identity(n)
  A = np.array([I,-M])

  # 描画インスタンスを生成
  ax = plt.axes()
  ax.set_aspect('equal')

  # 単位円を描画
  c = patches.Circle(xy=(0, 0), radius=1, fill=False)
  ax.add_patch(c)

  # 列ベクトルを描画
  for k in range(2):
    for i in range(n):
      a = A[k,:,i]
      ax.quiver(0, 0, a[0], a[1], angles='xy', scale_units='xy', scale=1)
      ax.text(a[0], a[1], 'a_{'+str(i)+','+str(k)+'}')

  # 相補錐を描画
  e1 = [1,0] # x軸
  for idx, (k1, k2) in enumerate(product([0,1], repeat=n)):
    a1 = A[k1,:,0]
    a2 = A[k2,:,1]

    # x軸とa1の成す角
    cos1 = np.inner(e1,a1)/(la.norm(e1)*la.norm(a1))
    rad1 = np.arccos(cos1)
    deg1 = np.rad2deg(rad1)
    
    # a1とa2の成す角
    cos2 = np.inner(a1,a2)/(la.norm(a1)*la.norm(a2))
    rad2 = np.arccos(cos2)
    deg2 = np.rad2deg(rad2)

    # x軸とa1の符号
    det1 = la.det([e1,a1])
    sign1 = np.sign(det1)

    # a1とa2の符号
    det2 = la.det([a1,a2])
    sign2 = np.sign(det2)

    # a1とa2のx軸からの符号付き角度
    theta1 = sign1*deg1
    theta2 = theta1+sign2*deg2

    # 符号付き角度を昇順に並び替え（180度以下の角度で描画するため）
    theta1, theta2 = sorted([theta1, theta2])

    # 相補錐を描画
    r = 2*(idx+1)/(2**n+1)
    a = patches.Arc((0, 0), height=r, width=r, theta1=theta1, theta2=theta2)
    ax.add_patch(a)

  # 描画
  plt.show()

  return

if __name__ == '__main__':

  n = 2
  M = np.random.rand(n,n)*2-1

  print(is_q_matrix(M))
  
  show_complementary_cone(M)