---
title: "Calculus7 Jacobian Matrix and Coordinate Transformation ヤコビ行列と座標変換"
author: "ryamada"
date: "2017年2月15日"
output:
html_document:
toc: true
toc_depth: 6
number_section: true
---
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```
$y = f(x)$という１変数(xの数が1)のスカラー関数(yが1変数)を$x$で微分すると、それは、$y=f(x)$という曲線の傾きを表す。
ヤコビ行列とヘッセ行列はその多変数版。
The derivative of a one-variable scaler function $y = f(x)$ indicates its slope.
Jacobian and Hessian Matrices are multi-variable versions of the slope.
遺伝統計学分野にとってのヤコビ行列は、傾きを表す行列としての意味と、変数変換・座標変換に関する情報を表す行列としての意味との2つの面で重要である。
There are two important ways to use Jacobian matrix in statistical genetics; one is the indicator of slope and the other is for variable transformation/coordinate transformation.
まず、傾きを表す行列としてのヤコビ行列(とヘッセ行列)を扱い、ついで、変数変換・座標変換でのヤコビ行列を扱う。
Its use for slope will be handled first and its use for coordinate transformation will come later.
正規分布のパラメタを最尤推定することにする。
平均$m$、SD $s$の正規分布からのサンプル $X=\{x_1,x_2,...,x_n\}$がある。
Assume the task of maximum likelihood estimation of parameters of normal distribution, mean $m$ and sd $s$ for samples $X=\{x_1,x_2,...,x_n\}$.
尤度関数は The likelihood function is;
$$
L(m,s|X) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}s}e^{-\frac{(x_i-m)^2}{2 s^2}} .
$$
対数尤度関数は The logarithm of it is;
$$
LL(m,s|X) = n \times (-\frac{1}{2}\log{2\pi} - \log{s}) - \sum_{i=1}^n \frac{(x_i-m)^2}{2 s^2}\\
=-\frac{1}{2s^2}(n\times m^2 -2\sum_{i=1}^n x_i \times m + \sum_{i=1}^n x_i^2) -n \log{s} -\frac{n}{2}\log{2\pi}
$$
この尤度関数と対数尤度関数は、$m,s$という2変数スカラー値関数(スカラー値が尤度)。
These functions are scaler functions with two variables $m$ and $s$.
$X= \{0,1,2\}$という観察をしたとする。
Assume $X= \{0,1,2\}$.
$$
LL(m,s|X) = -\frac{1}{2s^2} (3m^2-6m+5) -3\log{s} -\frac{3}{2}\log{2\pi}
$$
$(m,s)$の最尤推定解 $(\hat{m},\hat{s})$では、$LL$の偏微分が0($\frac{\partial LL(\hat{m},\hat{s})}{\partial m} =0,\frac{\partial LL(\hat{m},\hat{s})}{\partial s} =0$)であるから、以下の連立方程式をとけばよい。
We will get the maximum likelihood estimates $(\hat{m},\hat{s})$, by solving the system of two equations that are partial derivatives of $LL(m,s|X)$.
1つの(対数)尤度関数を各変数で偏微分して連立方程式を作ると、変数の数の方程式が得られる。
この例では、2つの変数、2つの方程式。
When we have one scaler (logarithm) likelihood function and partially differentiate it with $n$ parameters, we will have a system of $n$ equations.
In this case, two parameters and the equations.
$$
\frac{\partial LL(m,s)}{\partial m} = -\frac{1}{2s^2}(6m-6) = 0\\
\frac{\partial LL(m,s)}{\partial s} = \frac{1}{s^3}(3m^2-6m+5) - \frac{3}{s} = 0
$$
第１式より$\hat{m}=1$、それを用いて第２式より$\hat{s}^2 = \frac{2}{3}$を得る。
This can be solved easily as $\hat{m}=1$ and $\hat{s}^2 = \frac{2}{3}$.
この例のように、
$n$個の変数 $x_1,...,x_n$があって、同数$n$個の関数 $f_1(x_1,...,x_n),f_2(x_1,...,x_n),...,f_n(x_1,...,x_n)$ があるとき、
各関数を各変数で偏微分し、その偏微分関数を$n \times n $行列に並べたもの
$$
\partial F = \begin{pmatrix} \frac{\partial}{\partial x_1}f_1, ...,\frac{\partial}{\partial x_1}f_n\\
...,...,...\\
\frac{\partial}{\partial x_n}f_1, ...,\frac{\partial}{\partial x_n}f_n \end{pmatrix} = \begin{pmatrix}\frac{\partial}{\partial x_i}f_j \end{pmatrix}
$$
を$n$個の変数と$n$個の関数のヤコビ行列と言う。
When we have $n$ functions of $n$ variables, we can make $n\times n$ matrix whose elements are partial derivatives of each equation by each variable, as above. This $n\times n$ matrix whose elements are functions is called Jacobian matrix for the system of equations.
上の例では、1つの$n=2$変数スカラー関数を各変数で偏微分することによって$n$個の関数を得たが、それを$n$個の変数でされに偏微分することができるが、そのようにしてできた、$n\times n$行列を、元のスカラー関数にとってのヘッセ行列と言う。
We made the system of two equations by partially differentiate one scaler function of two variables and then made its Jacobian matrix. This Jacobian matrix is called Hessian matrix of the initial scaler function.
2つの方程式があった。
Recall the equations.
$$
f_m(m,s) = -\frac{1}{2s^2}(6m-6) = 0\\
f_s(m,s) = \frac{1}{s^3}(3m^2-6m+5) - \frac{3}{s} = 0
$$
$f_m(m,s) = u,f_s(m,s)=v$と書き直すと、
$m,s$が作る2次元平面の上に$(u,v)$という2次元空間上の点が対応づいて広がっているので、この様子を表示するには$(m,s,u,v)$という4次元空間が必要である。
Rewriting as $f_m(m,s) = u,f_s(m,s)=v$, each point $(m,s)$ in a two dimensional space has a value set $(u,v)$. This corresponds to a point $(m,s,u,v)$ which is a point in 4 dimensional space. It is not possible to visualize the 4 dimensional data as they are but we can draw something in 3D space as below.
少し工夫をして描いてみる。
u,vをそれぞれ別々に描けば
```{r,rgl=TRUE}
M <- seq(from=-1,to=2,length=80)
S <- seq(from=0.5,to=1.5,length=80)
ms <- expand.grid(M,S)
m <- ms[,1]
s <- ms[,2]
u <- -1/(2*s^2)*(6*m-6)
v <- 1/(s^3)*(3*m^2-6*m+5)-3/s
plot3d(m,s,u,main="u over m-s plane")
```
```{r rgl=TRUE}
plot3d(m,s,v,main = "v over m-s plane")
```
s=1に固定して、mと$f_m=u,f_s=v$とを描くと
Fix s at 1 and draw (u,v) along m axis,
```{r,rgl=TRUE}
s1 <- 1
us1 <- -1/(2*s1^2)*(6*M-6)
vs1 <- 1/s1^3*(3*M^2-6*M+5)-3/s1
plot3d(M,us1,vs1,type="l",main="(u,v) along m when s=1")
```
m=0.8に固定して、sと$f_m=u,f_s=v$とを描くと
Fix m at 0.8 and draw (u,v) for s axis,
```{r,rgl=TRUE}
m0 <- 0.8
um0 <- -1/(2*S^2)*(6*m0-6)
vm0 <- 1/S^3*(3*m0^2-6*m0+5)-3/S
plot3d(S,um0,vm0,type="l",main="(u,v) along s when m=0.8")
```
点$(m,s)$において、$(u=f_m(m,s),v=f_s(m,s))$ がどのように傾いているのかを調べる。
Let's evaluate tangents at $(m,s,u=f_m(m,s),v=f_s(m,s))$.
$s$を固定し、$m$方向での$u=f_m$の変化と$v=f_s$の変化を評価すると、それぞれ
Fix s and evaluate change of $u=f_m$ and $v=f_s$ along $m$;
$$
\frac{\partial}{\partial m}f_m(m,s) = -\frac{3}{s^2}\\
\frac{\partial}{\partial m}f_s(m,s) = \frac{6(m-1)}{s^3}
$$
$m$を固定し、$s$方向での$u=f_m$の変化と$v=f_s$の変化を評価すると、それぞれ
Fix m and evaluate change of $u=f_m$ and $v=f_s$ along $s$;
$$
\frac{\partial}{\partial s}f_m(m,s) = \frac{1}{s^3}(6m-6)\\
\frac{\partial}{\partial s}f_s(m,s) = -\frac{3}{s^4}(3m^2-6m+5) + \frac{3}{s^2}
$$
であり、結局、
These are summarized as;
$$
\begin{pmatrix} \frac{\partial}{\partial m}f_m(m,s),\frac{\partial}{\partial s}f_m(m,s)\\
\frac{\partial}{\partial m}f_s(m,s), \frac{\partial}{\partial s}f_s(m,s) \end{pmatrix} \\
=\begin{pmatrix} \frac{\partial ^2 }{\partial m\partial m}f(m,s),\frac{\partial^2}{\partial s\partial m}f(m,s)\\
\frac{\partial^2}{\partial m\partial s}f(m,s), \frac{\partial^2}{\partial s\partial s}f(m,s) \end{pmatrix}\\
= \begin{pmatrix} -\frac{3}{s^2}, \frac{1}{s^3}(6m-6) \\ \frac{6(m-1)}{s^3}, -\frac{3}{s^4}(3m^2-6m+5)+\frac{3}{s^2}\end{pmatrix}
$$
となる。
$f_m,f_s$に関するヤコビ行列であり、対数尤度関数に関するヘッセ行列である。
This is the Jacobian matrix for $f_m,f_s$ and the Hessian matrix for the log-likelihood function.
点$(mx,sx)$における接面を描いてみる。
Let's draw the tangents plane(s) at $(mx,sx)$.
```{r,rgl=TRUE}
mx <- m0
sx <- s1
fm.mxsx <- -1/(2*sx^2)*(6*mx-6)
fs.mxsx <- 1/(sx^3) *(3*mx^2-6*mx+5)-3/sx
dfmdm <- -3/(sx^2)
dfmds <- 1/(sx^3)*(6*mx-6)
dfsdm <- 6*(mx-1)/(sx^3)
dfsds <- -3/(sx^4)*(3*mx^2-6*mx+5)+3/(sx^2)
```
```{r rgl=TRUE}
plot3d(ms[,1],ms[,2],u)
spheres3d(mx,sx,fm.mxsx,radius = 1,col=3)
u.tan <- fm.mxsx + dfmdm * (m-mx) + dfmds * (s-sx)
spheres3d(m,s,u.tan,radius = 0.5,col=2)
```
```{r,rgl=TRUE}
plot3d(m,s,v)
spheres3d(mx,sx,fs.mxsx,radius = 10,col=3)
v.tan <- fs.mxsx + dfsdm * (m-mx) + dfsds * (s-sx)
spheres3d(m,s,v.tan,radius = 5,col=2)
```
s=1に固定して、mと$f_m=u,f_s=v$とを描くと
Fix s at 1 and draw (u,v) along m;
```{r,rgl=TRUE}
plot3d(M,us1,vs1,type="l")
spheres3d(mx,fm.mxsx,fs.mxsx,radius = 1,col=3)
us1.tan <- fm.mxsx + dfmdm * (M-mx) + dfmds * (s1-sx)
vs1.tan <- fs.mxsx + dfsdm * (M-mx) + dfsds * (s1-sx)
spheres3d(M,us1.tan,vs1.tan,radius = 0.1,col=2)
```
m=0に固定して、sと$f_m=u,f_s=v$とを描くと
Fix m at 0 and draw (u,v) along s;
```{r,rgl=TRUE}
plot3d(S,um0,vm0,type="l")
spheres3d(sx,fm.mxsx,fs.mxsx,radius = 0.5,col=3)
um0.tan <- fm.mxsx + dfmdm * (m0-mx) + dfmds * (S-sx)
vm0.tan <- fs.mxsx + dfsdm * (m0-mx) + dfsds * (S-sx)
spheres3d(S,um0.tan,vm0.tan,radius = 0.1,col=2)
```
上記の正規分布変数推定の例では、連立方程式は、簡単に解けてしまったが、ニュートン法を用いて、数値計算的に解を求める練習をすることにする。
The problem to estimate parameters of normal distribution was easy and we already know the answer. The followings are its numeric solution based on the Newton method as practice.
1変数関数 $f(x)=0$のニュートン法では、適当なxを初期値と定め、そこの接線を取り、接線が0と交わる値を次のxの値として取り、それを繰り返した。
In case of the Newton method for one parameter scaler function, we set an initial value and drew the tangent line then renewed the value with the coordinate where the tangent line crosses the line $0$ and repeated the procedures.
今回は(m,s)の2変数があり、2つの式があり、$(f_m(m,s),f_s(m,s))=(0,0)$となるように$(m,s)$の値の組を更新したい。
In our case we have two parameters to be estimated and we have two functions to handle $(f_m(m,s),f_s(m,s))$. We want to repeat procedures to renew $(m,s)$.
正規分布のパラメタ推定は、式がややこしいので、簡易な例をとることにする。
Here we will take a simpler function rather than the complicated function from normal distribution example.
最小化したい2変数スカラー関数を考える。
The following is the function.
$$
f(x,y) = (x+y)^4 + (x-y)^4
$$
```{r}
fr <- function(w) {
x <- w[1]
y <- w[2]
(x+y)^4 + (x-y)^4
}
t <- seq(from=-2,to=2,length=50)
xy <- as.matrix(expand.grid(t,t))
x <- xy[,1]
y <- xy[,2]
z <- apply(xy,1,fr)
plot3d(x,y,z)
```
答えはわかっていて、$x=y=0$。
The answer is known as $x=y=0$.
それを偏微分を用いて計算機的に求めてみる。
Let's solve it numerically with partial derivatives.
$$
\frac{\partial f}{\partial x} = f_x = 4(x+y)^3 + 4(x-y)^3\\
\frac{\partial f}{\partial y} = f_y = 4(x+y)^3 - 4(x-y)^3
$$
最小値をとる$(\hat{x},\hat{y})$では、$f_x(\hat{x},\hat{y})= 0, f_y(\hat{x},\hat{y})=0$。
ニュートン法を実施するために
$f_x,f_y$をさらに偏微分する。
Partially differentiate $f_x,f_y$ further with x and y.
$$
\frac{\partial f_x}{\partial x} = \frac{\partial ^2 f}{\partial x\partial x} = f_{xx} = 12(x+y)^2 + 12 (x-y)^2\\
\frac{\partial f_x}{\partial y} = \frac{\partial ^2 f}{\partial x\partial y} = f_{xy} = 12(x+y)^2 - 12 (x-y)^2\\
\frac{\partial f_y}{\partial x} = \frac{\partial ^2 f}{\partial y\partial x} = f_{yx} = 12(x+y)^2 - 12 (x-y)^2\\
\frac{\partial f_y}{\partial y} = \frac{\partial ^2 f}{\partial y\partial y} = f_{yy} = 12(x+y)^2 + 12 (x-y)^2
$$
$(x_0,y_0,f_x(x_0,y_0),f_y(x_0,y_0))$における接面と$f_x=0,f_y=0$との交点座標を$(x_1,y_1,0,0)$とする。
The tangent plane at $(x_0,y_0,f_x(x_0,y_0),f_y(x_0,y_0))$ crosses with the plane $f_x=0,f_y=0$; the crossing point is now called as $(x_1,y_1,0,0)$ and $(x_1,y_1)$ should replace $(x_0,y_0)$.
接面が$f_1=f_2=0$をよぎる$(m,s)$の座標を取り出すことは、
$(x_1,y_1)$ can be obtained by soving the equations below;
$$
f_x(x_1,y_1) = f_x(x_0,y_0) + \frac{\partial f_x}{\partial x} (x_1-x_0) + \frac{\partial f_x}{\partial y} (y_1-y_0) = 0\\
f_y(x_1,y_1) = f_y(x_0,y_0) + \frac{\partial f_y}{\partial x} (x_1-x_0) + \frac{\partial f_y}{\partial y} (y_1-y_0) = 0\\
$$
を解くことであるから、
行列
$$
\partial F = \begin{pmatrix} \frac{\partial}{\partial x}f_x(x,y),\frac{\partial}{\partial y}f_x(x,y)\\
\frac{\partial}{\partial y}f_y(x,y), \frac{\partial}{\partial y}f_y(x,y) \end{pmatrix} = \begin{pmatrix} f_{xx},f_{xy}\\ f_{yx},f_{yy}\end{pmatrix} = 12 \begin{pmatrix} (x+y)^2+(x-y)^2,(x+y)^2-(x-y)^2\\(x+y)^2-(x-y)^2,(x+y)^2+(x-y)^2 \end{pmatrix}
$$
を用いると、
Using this matrix, the Hessian matrix, the system is expressed as below;
$$
\begin{pmatrix} f_x(x_0,y_0) \\ f_y(x_0,y_0)\end{pmatrix} + \partial F \begin{pmatrix} x_1- x_0 \\ y_1- y_0 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}
$$
と表せる。
この行列 $\partial F$は式$f(x,y) = (x+y)^4 + (x-y)^4$のヘッセ行列である。
$$
\begin{pmatrix} x\\y \end{pmatrix} = w
$$
とすれば
Using this notation, the equation is expressed as below.
$$
F(w_0) + \partial F (w_1-w_0) = 0
$$
となり、結局
$$
\partial F w_1 = \partial F w_0 - F(w_0)\\
w_1 = w_0 - (\partial F)^{-1} F(w_0)
$$
によって求まる。
たとえば、
The following is an example.
$$
(x_0,y_0) = (3,4)
$$
のとき
$$
\partial F(3,4) = 12 \begin{pmatrix} (x+y)^2+(x-y)^2,(x+y)^2-(x-y)^2\\(x+y)^2-(x-y)^2,(x+y)^2+(x-y)^2 \end{pmatrix} = 12 \begin{pmatrix} 50, 48\\48,50\end{pmatrix}
$$
であるから、
$$
(\partial F(3,4))^{-1} = (12 \begin{pmatrix} 50, 48\\48,50\end{pmatrix})^{-1} .
$$
実際に計算すれば
The inverse matrix can be calculated;
```{r}
fx.F <- function(x,y){
tmp1 <- (x+y)^2
tmp2 <- (x-y)^2
12 * matrix(c(tmp1+tmp2,tmp1-tmp2,tmp1-tmp2,tmp1+tmp2),2,2)
}
x0 <- 1
y0 <- 2
F.inv <- solve(fx.F(x0,y0))
F.inv
```
となり、
```{r}
w0 <- matrix(c(x0,y0),nrow=2)
fxfy <- function(w){
c(4*(w[1]+w[2])^3+4*(w[1]-w[2])^3,4*(w[1]+w[2])^3-4*(w[1]-w[2])^3)
}
w1 = w0 - F.inv %*% matrix(fxfy(w0),nrow=2)
w1
```
この$w_1=(x_1,y_1)$を使って$(x_1,y_1,f_x(x_1,y_1),f_y(x_1,y_1))$の接面と$f_x=0,f_y=0$との交点座標を$(x_2,y_2,0,0)$とすれば、
The $(x,y)$ can be renewed further.
```{r}
w2 <- w1 - solve(fx.F(w1[1],w1[2])) %*% matrix(fxfy(w1),nrow=2)
w2
```
これを繰り返すことにする。
Repeat the procedures as below.
```{r}
x0 <- 3
y0 <- 4
w <- matrix(c(x0,y0),nrow=2)
xy <- matrix(c(x0,y0),nrow=1)
n.iter <- 10
for(i in 1:n.iter){
new.w <- w - solve(fx.F(w[1],w[2])) %*% matrix(fxfy(w),nrow=2)
w <- new.w
xy <- rbind(xy,c(w))
}
plot(xy,type="b",pch=20)
matplot(xy,type="l")
```
正規分布の変数$m,s$にニュートン法を適用してみる。以下の説明、コードを理解せよ。
The followings are Newton's method applied to the problem of normal distribution. Read and understand the steps and the codes.
正規分布の変数推定対数尤度関数から、2つの方程式$f_m(m,s)=0$, $f_s(m,s) = 0$が得られる。
Two equations were obtained from Log-likelihood function.
このヤコビ行列を作る。
Make the Jacobian matrix.
ヤコビ行列を利用して、
接面が$u=f_1(m,s)=0$, $v=f_2(m,s)=0$と交わる$(m,s)$を採用する。
Renew the initial (m,s) using the Jacobian matrix.
```{r rgl=TRUE}
mx <- 0.6
sx <- 1
plot3d(m,s,u,col=gray(0.3),size=0.01)
spheres3d(mx,sx,fm.mxsx,radius = 1,col=3)
u.tan <- fm.mxsx + dfmdm * (m-mx) + dfmds * (s-sx)
spheres3d(m,s,u.tan,radius = 0.5,col=2)
spheres3d(m,s,rep(0,length(m)),radius = 0.5,col=4)
```
接面が$u=0$平面と交わって直線が現れる。
The tangent plane corsses with the plane $u=0$, showing a straight line.
```{r,rgl=TRUE}
plot3d(m,s,v,col=gray(0.3),size=0.01)
spheres3d(mx,sx,fs.mxsx,radius = 10,col=3)
v.tan <- fs.mxsx + dfsdm * (m-mx) + dfsds * (s-sx)
spheres3d(m,s,v.tan,radius = 1,col=2)
spheres3d(m,s,rep(0,length(m)),radius = 1,col=4)
```
接面が$v=0$平面と交わって直線が現れる。
The tangent plane corsses with the plane $v=0$, showing a straight line.
この２つの直線の交点が、ニュートン法が示す、次の$(m1,s1)$の位置。
$(m1,s1)$ is the crossing point of these two lines.
```{r,rgl=TRUE}
plot3d(mx,sx,fm.mxsx,radius = 1,col=3,zlim = c(-10,10))
spheres3d(mx,sx,fs.mxsx,radius = 1,col=5,zlim = c(-10,10))
spheres3d(m,s,u.tan,radius = 0.2,col=3,zlim = c(-10,10))
spheres3d(m,s,v.tan,radius = 0.2,col=5,zlim = c(-10,10))
spheres3d(m,s,rep(0,length(m)),radius = 0.2,col=4,zlim = c(-10,10))
```
```{r}
matrix(c(dfmdm,dfmds,dfsdm,dfsds),byrow=TRUE,nrow=2)
```
```{r}
new.ms <- matrix(c(mx,sx),nrow=2) - solve(matrix(c(dfmdm,dfmds,dfsdm,dfsds),byrow=TRUE,nrow=2)) %*% matrix(c(fm.mxsx,fs.mxsx),nrow=2)
new.ms
```
```{r,rgl=TRUE}
plot3d(mx,sx,fm.mxsx,radius = 0.5,col=3)
spheres3d(mx,sx,fs.mxsx,radius = 0.5,col=5)
spheres3d(m,s,u.tan,radius = 0.2,col=3)
spheres3d(m,s,v.tan,radius = 0.2,col=5)
spheres3d(m,s,rep(0,length(m)),radius = 0.2,col=1)
spheres3d(new.ms[1],new.ms[2,],0,radius = 0.5,col=6)
```
$\partial F$を計算する関数を使って
```{r}
partialF <- function(m,s){
matrix(c(-3/(s^2),1/(s^3)*(6*m-6),6*(m-1)/(s^3),-3/(s^4)*(3*m^2-6*m+5) + 3/(s^2)),byrow=TRUE,2,2)
}
m0 <- mx
s0 <- sx
n.iter <- 10
ms <- matrix(c(m0,s0),nrow=1)
solve(partialF(ms[1,1],ms[1,2]))
for(i in 1:n.iter){
pF <- partialF(ms[i,1],ms[i,2])
tmpf1.mxsx <- -1/(2*ms[i,2]^2)*(6*ms[i,1]-6)
tmpf2.mxsx <- 1/(ms[i,2]^3) *(3*ms[i,1]^2-6*ms[i,1]+5)-3/ms[i,2]
new.ms <- matrix(ms[i,],nrow=2) - solve(pF) %*% matrix(c(tmpf1.mxsx,tmpf2.mxsx),nrow=2)
ms <- rbind(ms,c(new.ms))
}
plot(ms,type="b")
points(1,sqrt(2/3),pch=20,col=2)
matplot(ms,type="l")
abline(h=c(1,sqrt(2/3)),col=3)
```
$$
u = \frac{\partial}{\partial x} f(x,y) = f_{x}(x,y)\\
v = \frac{\partial}{\partial x} f(x,y) = f_{y}(x,y)
$$
この関係を
These equations can be considered as coordinate transformation.
$$
u = f(x,y)\\
v = g(x,y)
$$
と捉えなおす。
逆に The opposite transformation;
$$
x = h(u,v)\\
y = k(u,v)
$$
も成り立つとする。
先ほどまで$(x,y)$と言う座標で呼ばれてた点が、$(u,v)$という座標で呼ばれるようになった。言い換えれば、そのような座標を引きなおした、とも読み取れる。
The point who had $(x,y)$ coordinate, has now different coordinate $(u,v)$. On the same space, new coordinate-system lattice appeared, in other words.
$$
u = 2x+\frac{1}{2}y\\
v = x+y
$$
逆に(u,v)という座標でよばれていた点を(x,y)という座標で呼ぶように変えるのが、以下の処理である。
The opposite relation is;
$$
x = \frac{2u-v}{3}\\
y = \frac{-2u + 4v}{3}
$$
上記の関係だと、新たな座標系による格子は次のようになる。
The lattice with new coordinate system is drawn as below.
左のパネルでは$x$の値が等しい点を結んだ線と、$y$の値が等しい点を結んだ線とが現れている。
右のパネルでは$u$の値が等しい点を結んだ線と、$v$の値が等しい点を結んだ線が現れている。
In the left panel, lines indicate sets of points with same x values or same y values.
In the right panel, lines indicate sets of points with same u values or same v values.
元の座標で点を打ち、元の座標で格子を描くと、マス目状の格子が描かれ、元の座標で点を打ち、新たな座標で格子を描くと、ゆがんだ格子が描かれる。
When points are placed in the original (x,y) coordinates and lattice of the original coordinate, then the lattice appears as perpendicular lattice.
When lattice is drawn with (u,v) system, then the lattice lines are curved.
```{r}
n <- 10^4
XY <- matrix(runif(n*2)*5,ncol=2)
x <- XY[,1]
y <- XY[,2]
col.ori <- abs(floor(x-min(x)+1)*floor(y-min(y)+1))
par(mfcol=c(1,2))
plot(x,y,col=col.ori)
u <- 2*x+1/2*y
v <- x+y
col <- abs(floor(u-min(u)+1) * floor(v-min(v)+1))
plot(x,y,col=col,pch=20)
par(mfcol=c(1,1))
```
逆に、新たな座標で点を打ち、それに元の座標で格子を描くと、ゆがみ方が逆になったパターンが現れる。
The opposite;
```{r}
n <- 10^4
UV <- matrix(runif(n*2)*5,ncol=2)
u <- UV[,1]
v <- UV[,2]
col.ori <- abs(floor(u-min(u)+1)*floor(v-min(v)+1))
par(mfcol=c(1,2))
plot(u,v,col=col.ori)
x <- (2*u-v)/3
y <- (-2*u+4*v)/3
col <- abs(floor(x-min(x)+1) * floor(y-min(y)+1))
plot(u,v,col=col,pch=20)
par(mfcol=c(1,1))
```
$0 \le x \le 1, 0 \le y \le 1$という正方領域に
$f(x,y) = x^2+y$という関数があって
Assume a function $f(x,y) = x^2+y$ in the area $0 \le x \le 1, 0 \le y \le 1$, that is a square and we want to know the integral below.
$$
\int_{x=0}^{x=1} \int_{y=0}^{y=1} f(x,y) dxdy
$$
という積分をしたいとする。
正方領域を十分に小さな正方形に分割して、足し合わせることで近似してみる。
小さな正方形の面積で重みづけをした和になる。
Divide the square into small squares and sum the value of the points multiplied by their area, then it is good approximation of the integral.
```{r}
n <- 10
tmp <- seq(from=0,to=1,length=n+1)
tmp <- tmp[-(n+1)]
xy <- as.matrix(expand.grid(tmp,tmp))
x <- xy[,1]
y <- xy[,2]
s <- (tmp[2]-tmp[1])*(tmp[2]-tmp[1])
fxy <- x^2+y
sum(fxy * s)
```
分割を細かくすると、真の値に近づくようだ。
With the division finer, the approximation gets closer to the true value.
```{r}
ns <- 2^(1:10)
ret <- rep(0,length(ns))
for(i in 1:length(ns)){
n <- ns[i]
tmp <- seq(from=0,to=1,length=n+1)
tmp <- tmp[-(n+1)]
xy <- as.matrix(expand.grid(tmp,tmp))
x <- xy[,1]
y <- xy[,2]
deltax <- deltay <- tmp[2]-tmp[1]
s <- deltax * deltay
fxy <- x^2+y
ret[i] <- sum(fxy * s)
}
plot(ns,ret)
```
ここで、小さい面積を$\delta x \times \delta y$として計算した。
We used $\delta x \times \delta y$ as the area of small square.
これは、横軸方向の小さいベクトル$(\delta x,0)$と縦軸方向の小さいベクトル$(0,\delta y)$との作る小さな正方形の面積が$\delta x \times \delta y$と計算できることを利用している。
This squre is consisted of two mutually perpendicular small vectors, $(\delta x,0)$ and $(0,\delta y)$.
では、同様に、$(u,v)$を使うとどうなるだろうか。
Now we want to do the similar for (u,v).
$0 \le u \le 1, 0 \le v \le 1$の領域について、(u,v)の細かい格子を作って
$f(x,y) =f((2*u-v)/3, (-2*u+4*v)/3) = ((2*u-v)/3)^2 + (-2*u+4*v)/3$ の値を
足し合わせればよい。
Make fine lattice and sum up the values $f(x,y) =f((2*u-v)/3, (-2*u+4*v)/3) = ((2*u-v)/3)^2 + (-2*u+4*v)/3$ multiplied by their area.
このときの小さい格子の作る小さい領域は
$\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ と言う小さいベクトルと
$\delta v = \delta x - \delta y$と言う小さいベクトルが作る小領域である。
The small area is approxiameted as paralleogram between $\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ and $\delta v = \delta x - \delta y$.
この面積は$|\delta| u \times |\delta v|$というように小ベクトルの長さの積にはならない。
この小面積は平行四辺形なので、平行四辺形の面積の求め方を確認することにする。
The area of parallelogram is not $|\delta| u \times |\delta v|$.
We should know how to calculate the area of parallelogram.
ベクトル$A=(a,b)$とベクトル$X=(x,y)$とが作る平行四辺形の面積は、２つのベクトルの成す角を$\theta$とすれば
Assume two vecors $A=(a,b)$ and $X=(x,y)$, the angle between which is $\theta$, that make a parallelogram.
Its area is;
$$
S = |A| |X| \sin{\theta} \\
= |A||X|\sqrt{1-\cos^2{\theta}} \\
= \sqrt{|A|^2|X|^2-(A \cdot X)^2}
$$
ただし、$A\cdot X$は内積。
$(A \cdot X)$ is inner product.
これを成分表示すれば Express the area with elements of vectors;
$$
S = \sqrt{(a^2+b^2)\times (x^2+y^2) - (ax+by)^2}\\
= \sqrt{a^2y^2+b^2x^2 - 2abxy}
= |ay-bx|
$$
この$|ay-bx|$は、行列
$$
\begin{pmatrix}a,x\\b,y\end{pmatrix}
$$
の行列式の絶対値。
$|ay-bx|$ is the absolute value of determinant of the matrix above.
変数変換によって求めたい微小面積は
$\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ と言う小さいベクトルと
$\delta v = \delta x + \delta y$とが作る平行四辺形であるから、
The small area for (u,v) coordinate is parallelogram by $\delta u = 2 \times \delta x + \frac{1}{2} \delta y$ and
$\delta v = \delta x + \delta y$.
それは、
$$
\begin{pmatrix} 2 \delta x, \delta x \\ \frac{1}{2} \delta y, \delta y\end{pmatrix}
$$
という行列の行列式の絶対値。
It is the absolute value of determinant of the matrix above.
それは
$$
2\times \delta x \times \delta y - \delta x \times \frac{1}{2} \delta y = (2\times 1 - 1\times \frac{1}{2}) \times \delta x \delta y = \frac{3}{2} \delta x \delta y.
$$
ここで$2\times 1 - 1 \times \frac{1}{2}$というのは
$$
\begin{pmatrix} 2,1\\\frac{1}{2},1 \end{pmatrix} = \begin{pmatrix} \frac{\partial x}{\partial u} , \frac{\partial x}{\partial v}\\ \frac{\partial y}{\partial x},\frac{\partial y}{\partial v} \end{pmatrix}
$$
の行列式。
これは、ヤコビ行列の行列式であることを意味している。
The above formulas have shown, the matrix for the area calculation is Jacobian matrix of coordinate transformation equations.
結局、
$$
dS = |det(Jacobi)| dudv = dxdy
$$
であるから、
$0 \le u \le 1, 0 \le v \le 1$の領域について、
$f(x,y) =f((2u-v)/3, (-2u+4v)/3) = ((2u-v)/3)^2 + (-2u+4v)/3$を積分するには
以下を計算することとなる。
The integral of the problem is expressed as;
$$
\int_{u=0}^{u=1} \int_{v=0}^{v=1} ((2u-v)/3)^2 + (-2u+4v)/3 \times |det \begin{pmatrix} \frac{\partial x}{\partial u} , \frac{\partial x}{\partial v}\\ \frac{\partial y}{\partial x},\frac{\partial y}{\partial v} \end{pmatrix})| dudv \\
= \int_{u=0}^{u=1} \int_{v=0}^{v=1} ((2u-v)/3)^2 + (-2u+4v)/3 \times |det\begin{pmatrix} 2,1\\\frac{1}{2},-1 \end{pmatrix} | dudv
$$
ここで、
$\begin{pmatrix} 2,1\\\frac{1}{2},-1 \end{pmatrix}$は、u,vの関数ではないので、$dS = C dudv = dxdy$のように、微小面積は座標によらず定数倍($C=\frac{3}{2}$倍)であることがわかるが、それは、図を見ると、平行四辺形の格子がどこも同じ大きさであることに対応している。
In this case, the determinant of Jacobian matrix is constant regardless of (u,v), that were shown in the figure where all parallelograms of lattice were identical.
微小面積の変換が定数倍ではない場合、ヤコビ行列の行列式がu,vの関数である例を練習問題として取り上げる。
The following exercise is the case of non-constant Jacobian determinant.
平面の極座標変換
$$
x = r \cos{\theta}\\
y = r \sin{\theta}
$$
のヤコビ行列を求めよ。
This indicates polar coordinate transformation in 2D.
Answer its Jacobian matrix.
その行列式が $r$であることを示せ。
Show its determinant is $r$.
半径1の球の体積の半分を求めることを考える。
We want to know the volume of unit sphere.
$x=r\cos{\theta},y=r\sin{\theta}$における、単位球面のz座標は$\pm \sqrt{1-r^2}$であるから、
z-coordinate for $x=r\cos{\theta},y=r\sin{\theta}$ is $\pm \sqrt{1-r^2}$.
半径rを0から1まで、角$\theta$を0から$2\pi$まで$\sqrt{1-r^2}$について積分すれば、単位球の半分の値となる。
Integral $\sqrt{1-r^2}$ from 0 to 1 for r and from 0 to $2\pi$ for $\theta$ is the half of the volume.
微小面積$dS= |det(Jacobi)| drd\theta = r drd\theta$であるから、半球の体積は
The infinitesimal area is $dS= |det(Jacobi)| drd\theta = r drd\theta$, then;
$$
\int_{\theta=0}^{\theta=2\pi} \int_{r=0}^{r=1} \sqrt{1-r^2} (r dr d\theta)
$$
これを用いて、単位球の体積が$\frac{4}{3}\pi$であることを示せ。
Show that the volume of unit sphere is $\frac{4}{3}\pi$, based on the information above.
