Skip to content

Instantly share code, notes, and snippets.

@wsuzume

wsuzume/ishitah.ipynb

Last active June 15, 2018 16:04
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save wsuzume/0003d764d53b1cabdd0616f21446376d to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save wsuzume/0003d764d53b1cabdd0616f21446376d to your computer and use it in GitHub Desktop.
Download ZIP
Solution for ishitah's problem
Display the source blob
Display the rendered blob
Raw
ishitah.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# (M^T)PM = I を満たすPを求めよ（ただしMは3x2行列）"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"np.random.seed(0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 適当なMを生成"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.5488135 , 0.71518937],\n",
" [ 0.60276338, 0.54488318],\n",
" [ 0.4236548 , 0.64589411]])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = np.random.random((3,2))\n",
"M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 適当な3x3行列を特異値分解して直交行列Vを得る（解の不定性）"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[-0.40343826, -0.76554323, -0.50117974],\n",
" [ 0.34473823, 0.38020044, -0.85825589],\n",
" [-0.84758074, 0.51902909, -0.11052461]])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U, s, V = np.linalg.svd(np.random.random((3,3)))\n",
"V"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.00000000e+00, 6.83388169e-17, 4.78915481e-17],\n",
" [ 6.83388169e-17, 1.00000000e+00, 1.24204090e-16],\n",
" [ 4.78915481e-17, 1.24204090e-16, 1.00000000e+00]])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V.dot(V.T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mに適当なベクトルを1本追加（解の不定性）"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.5488135 , 0.71518937, 0.0871293 ],\n",
" [ 0.60276338, 0.54488318, 0.0202184 ],\n",
" [ 0.4236548 , 0.64589411, 0.83261985]])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M_a = np.concatenate([M, np.random.random(3)[:, np.newaxis]], axis=1)\n",
"M_a"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 適当なベクトルを追加したMをVの空間に写像"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"array([[-0.89518099, 0.05476367, -0.19913632],\n",
" [-1.02937543, -0.10062448, -0.39475771],\n",
" [-0.46792155, -0.67687704, -0.15538017]])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = V.dot(M_a).T\n",
"A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 以下で表される特性行列Bを計算"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"B = np.array([[A[0,0] ** 2, A[0,1] ** 2, A[0,2] ** 2],\n",
" [A[1,0] ** 2, A[1,1] ** 2, A[1,2] ** 2],\n",
" [ (A[0,0] * A[1,0]), (A[0,1] * A[1,1]), (A[0,2] * A[1,2]) ]])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bの逆行列をベクトル (1,1,0)^T にかけるとPの固有値が求まる"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1.32600029, 95.83012756, -8.82579404])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lam = np.linalg.inv(B).dot(np.array([1,1,0]))\n",
"lam"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### (V^T)ΛV でPが計算できる"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 5.26431201, 16.85258856, -28.91228336],\n",
" [ 16.85258856, 12.25199231, -30.25520993],\n",
" [-28.91228336, -30.25520993, 70.81402949]])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = V.T.dot(np.diag(lam)).dot(V)\n",
"P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 本当に解けているかどうか確認"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.00000000e+00, -1.03039624e-15],\n",
" [ -1.40402649e-15, 1.00000000e+00]])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.T.dot(P).dot(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### (Q^T)Q = P なる Q が欲しければ以下のように求めればよいがけっこう固有値が負になりエラーを吐く"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/yoshinobu/.pyenv/versions/anaconda3-5.0.1/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt\n",
" \"\"\"Entry point for launching an IPython kernel.\n"
]
}
],
"source": [
"Q = V.dot(np.diag(1 / np.sqrt(lam)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 複素数の範囲でよければ上の式も計算できることが多い"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@wsuzume
Copy link
Author

wsuzume commented Jun 15, 2018

解法 https://github.com/wsuzume/notes/blob/master/ishitah.pdf

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment