notebooks/LA_integer_matrix.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 整数行列\n",
    "\n",
    "各成分が整数の行列を整数行列と言います。\n",
    "\n",
    "$A$を$m\\times n$の整数行列とします。\n",
    "整数行列の基本変形は次のように整数行列を保ちます：\n",
    "\n",
    "* ある行（列）を−１倍する\n",
    "* ２行（列）を入れ替える\n",
    "* ある行（列）の整数倍を他の行（列）に加える\n",
    "\n",
    "$Q$を$m\\times n$の整数行列とします。\n",
    "このとき、\n",
    "$\\det Q=\\pm 1$である$Q$をユニモジュラ行列と言います。\n",
    "\n",
    "また、\n",
    "* $\\det Q =\\pm 1$と$Q^{-1}$が整数行列であることは同値です。（[証明](https://shakayami-math.hatenablog.com/entry/2019/12/27/212816)）\n",
    "\n",
    "$A$に列基本変形を繰り返し適用するのは、ユニモジュラ行列$Q$を使って$AQ$とすることに対応します。\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 行列式因子\n",
    "\n",
    "$A$の$k$次章行列式の最大公約数$d_k(A)$は行列式因子と呼ばれます。\n",
    "$d_k(A)$はユニモジュラ行列による基本変形によって不変であることに注意しましょう。\n",
    "実際、\n",
    "\n",
    "* ある行を$-1$倍したり、行を入れ替えても、新たにできる小行列式は元の$-1$倍\n",
    "* ある行の整数倍を他の行に加えた場合、新たにできる少行列式は元の少行列式と既存の少行列式の整数倍の和"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hermite標準形\n",
    "\n",
    "rank $m$の整数行列$A: m \\times n$に対して、次のような形に変形できるユニモジュラ行列$Q$が存在します。\n",
    "\n",
    "![hermite](figs/hermite_noraml.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "### 元の行列 ###\n",
      "[[1 4 4 5]\n",
      " [1 6 8 2]\n",
      " [6 8 8 9]]\n",
      "### １行目と２行目を入れ替え ###\n",
      "[[1. 4. 4. 5.]\n",
      " [6. 8. 8. 9.]\n",
      " [1. 6. 8. 2.]]\n",
      "### １列目と２列目を入れ替え ###\n",
      "[[1. 4. 4. 5.]\n",
      " [1. 8. 6. 2.]\n",
      " [6. 8. 8. 9.]]\n",
      "### ２行目を２倍 ###\n",
      "[[ 1.  4.  4.  5.]\n",
      " [ 1.  6.  8.  2.]\n",
      " [12. 16. 16. 18.]]\n",
      "### ２行目を３倍して１行目に追加 ###\n",
      "[[ 1.  4.  4.  5.]\n",
      " [19. 30. 32. 29.]\n",
      " [ 6.  8.  8.  9.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 基本行列を作ります\n",
    "\n",
    "def make_transposition_matrix(size, i, j):\n",
    "    # 行列のiとj番目の行（列）を入れ替える行列を作ります\n",
    "    P = np.eye(size)\n",
    "    P[i, i] = 0\n",
    "    P[j, j] = 0\n",
    "    P[i, j] = 1\n",
    "    P[j, i] = 1\n",
    "    return P\n",
    "\n",
    "M, N = 3, 4\n",
    "A = np.random.randint(1, 10, (M, N))\n",
    "print(\"### 元の行列 ###\")\n",
    "print(A)\n",
    "\n",
    "print(\"### １行目と２行目を入れ替え ###\")\n",
    "P = make_transposition_matrix(M, 1, 2)\n",
    "print(P @ A)\n",
    "\n",
    "print(\"### １列目と２列目を入れ替え ###\")\n",
    "P = make_transposition_matrix(N, 1, 2)\n",
    "print(A @ P)\n",
    "\n",
    "\n",
    "def make_scale_matrix(size, i, s):\n",
    "    # 行列のi行（列）をs倍する行列を作ります　\n",
    "    T = np.eye(size)\n",
    "    T[i, i] = s\n",
    "    return T\n",
    "\n",
    "print(\"### ２行目を２倍 ###\")\n",
    "T =  make_scale_matrix(M, 2, 2)\n",
    "print(T @ A)\n",
    "\n",
    "\n",
    "def make_add_matrix(size, i, j, s):\n",
    "    # 行列のj行（列）をs倍してi行（列）に足す行列を作ります　\n",
    "    E = np.eye(size)\n",
    "    E[i, j] = s\n",
    "    return E\n",
    "\n",
    "print(\"### ２行目を３倍して１行目に追加 ###\")\n",
    "E =  make_add_matrix(M, 1, 2, 3)\n",
    "print(E @ A)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 9  4  0  1  9  0  1]\n",
      " [ 8  9  0  8  6  4  3]\n",
      " [ 0  4  6  8  1  8  4]\n",
      " [ 9  4  0  1  9  0  1]\n",
      " [27 12  0  3 27  0  3]]\n",
      "[[1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0.]\n",
      " [1. 0. 0. 0. 0. 0. 0.]\n",
      " [3. 0. 0. 0. 0. 0. 0.]]\n",
      "[[1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0.]\n",
      " [1. 0. 0. 0. 0. 0. 0.]\n",
      " [3. 0. 0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Hermite標準形\n",
    "\n",
    "# M x N でrank=Rの行列\n",
    "R, M, N = 3, 5, 7\n",
    "np.random.seed(10)\n",
    "A = np.random.randint(0, 10, (R, N))\n",
    "A_left = np.repeat(A[0][None, :], M - R, axis=0)\n",
    "A_left = np.random.randint(0, 10, M - R).reshape(-1, 1) * A_left\n",
    "A = np.vstack([A, A_left])\n",
    "A_origin = A.copy()\n",
    "\n",
    "\n",
    "Q = np.eye(N)  # 右からかける正則行列\n",
    "\n",
    "assert A.shape == (M, N)\n",
    "print(A)\n",
    "\n",
    "all_zero_count = 0\n",
    "while True:\n",
    "    if all_zero_count == M:\n",
    "        break\n",
    "\n",
    "    all_zero_count = 0\n",
    "    for m in range(M):\n",
    "        if len(np.nonzero(A[m, m+1:])[0]) == 0:\n",
    "            all_zero_count += 1\n",
    "            continue\n",
    "        # ある行について、n以上の列の非零の絶対値最小成分を見つけます\n",
    "        abs_min_idx = np.argmin(np.where(A[m, m:]==0, np.infty, np.abs(A[m, m:])))\n",
    "\n",
    "        # 絶対値最小の要素が対角成分にくるように入れ替えます\n",
    "        idx = abs_min_idx + m\n",
    "        P_col = make_transposition_matrix(N, m, idx)\n",
    "        A = A @ P_col\n",
    "        Q = Q @ P_col\n",
    "\n",
    "        # (m, m)番目の項を使って他の行の絶対値を小さくします\n",
    "        for n in range(N):\n",
    "            if n == m:\n",
    "                continue\n",
    "            scale = A[m, n] // A[m, m]\n",
    "            E = make_add_matrix(N, m, n, -scale)\n",
    "            A = A @ E\n",
    "            Q = Q @ E\n",
    "\n",
    "\n",
    "print(np.where(np.abs(A) < 1e-6, 0, A))\n",
    "\n",
    "# 基本変形をまとめたやつで変形してみます\n",
    "A = A_origin @ Q\n",
    "\n",
    "print(np.where(np.abs(A) < 1e-6, 0, A))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hermite標準形について、次の４つの条件は同値です\n",
    "\n",
    "* $A$のHermite標準形が$[I O]$\n",
    "* 行列式因子が$d_m(A)=1$\n",
    "* 任意の整数ベクトル$b$について、$Ax=b$が整数解を持つ\n",
    "* $yA$が整数ベクトルとなる任意の$y$は整数ベクトル"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smith標準形\n",
    "\n",
    "Hermite標準形は左下に非零要素が残っちゃいますが、更に変形すれば、より行列式因子が計算しやすいSmith標準形に変形できます。\n",
    "\n",
    "![smith](figs/smith_normal.png)\n",
    "\n",
    "ここで、$\\alpha_2$は$\\alpha_1$で割り切れる、$\\alpha_3$は$\\alpha_2$で割り切れる、... のような形になっています。\n",
    "これを利用すると、行列式因子はこの最初の左上の成分だけを調べるだけで解決できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}