타원 상의 한 점과 타원 밖의 한 점 사이의 거리 구하기

이번 포스팅은 다소 간단한 내용을 담고 있는 포스팅인데, 제가 실제로 하고 있는 없무와 약간 관련성이 있기 때문에 하는 포스팅입니다. 

 

문제는 타원과 타워 밖의 한 점이 주어져 있을 때, 타원과 점 사이의 거리를 구하는 것 입니다. 문제를 최대한 간단하게 하기 위해서 2차원 공간을 생각하겠습니다. 2차원 공간에서 타원의 방정식은 

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

로 주어집니다. 타원 위의 일반적인 점을 $\vec{x} = (x, y)$라 놓겠습니다. 타원 밖의 한 점의 좌표를 $\vec{p} = (p, q)$ 라고 하겠습니다. 이 점은 타원의 밖에 있습니다. 일반화를 잃지 않고 타원과 타원 밖의 점의 배치를 위 그림과 같이 놓을 수 있습니다. 타원의 중심이 원점이 아니라면, 평행 이동을 통해 타원의 중심을 원점으로 놓을 수 있습니다. 타원의 축 방향이 X, Y 축과 나란하지 않다면 좌표 회전을 통해서 타원의 축 방향이 X, Y 축과 나란하게 놓을 수 있습니다 .또한 점 $(p, q)$가 제 1사분면에 놓여 있지 않은 경우,  X축 훅은 Y축을 기준으로 거울 반사 이동 시켜서 제 1사분면에 놓을 수 있습니다. 평행이동, 회전 변환, 거울 반사 이동은 점과 점 사이의 거리를 보존하는 변환이기 때문에, 이와 같은 변환을 하더라도 문제의 정답은 전혀 달라지지 않습니다. 

 

우리 문제는 $d(\vec{x}) = |\vec{x} - \vec{p}|$를 최소화하는 타원 위의 점 $\vec{x}$를 찾는 것 입니다. 좀 더 명확하게 수학적으로 표현하면,

$$\vec{x}_{\text{opt}} = \arg \min |\vec{x} - \vec{p}|, \text{ subject to } \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

입니다. 매우 단순한 형태의 구속 조건이 있는 하의 최적화 문제 입니다. 

 

1. 타원 상의 점을 매개변수화를 이용한 풀이

 

(가장?) 간단하게 이 문제에 접근 할 수 있는 방법은 타원 상의 점을 매개변수화 하여 표현하고, $|\vec{x} - \vec{p}|$ 를 매개변수로 표현하고, 이를 매개변수로 미분하여 미분이 0이 되는 값을 찾는 것 입니다. 타원상의 점은 $\theta = [0, 2 \pi]$라 했을 때,

$$x = a \cos \theta$$

$$y = b \sin \theta$$

으로 표현할 수 있습니다. 따라서 $d(\vec{x}) = d(\theta) = \sqrt{(a \cos \theta - p)^2 + (b \sin \theta -b)^2}$ 가 됩니다. $d(\theta)$는 표현식에 근호가 있어 복잡하니 이를 제곱한 값을 사용하도록 하겠습니다. 즉, $f(\theta) = d(\theta)^2 = (a \cos \theta - p)^2 + (b \sin \theta -b)^2$ 라고 놓고, 

$$\frac{df}{d\theta} = 0$$

이 되게 하는 $\theta$를 찾으면 됩니다.

$$\frac{df}{d\theta} = -2a\sin \theta (a \cos \theta - x) + 2 b \cos \theta (b \sin \theta -y) $$

인데, $\cot \theta$를 $t$라고 놓으면, $\sin \theta = \sqrt{1- t^2}$이 되고(물론 $\pm$ 부호가 앞에 붙어야 하지만, 양수만을 채택 하도록 하겠습니다. 그 이유는 나중에 설명합니다), 위 식을 $t$에 대해서 정리하면,

$$f(t) = -2a \sqrt{1-t^2} (at -x) + 2bt (b\sqrt{1-t^2}-y)$$

가 됩니다. $f(t) = 0$을 구해야 하는 것이 목적이니, 이는 

$$\sqrt{1-t^2}(t(-a^2 + b^2) + at) = bty$$

와 같고, 루트를 없애기 위해서 양변을 제곱하여 정리하면,

$$(1-t^2)\{ t(-a^2 + b^2) +ax\}^2 - b^2y^2t^2 = 0$$

이 됩니다. 이 식을 $t$에 대해서 정리하면,

$$-(a^2 - b^2)^2t^4 + 2ax(a^2 - b^2)t^3 + \{a^4 - a^2 (2b^2 + x^2) + b^4 - b^2y^2\}t^2 + 2ax(b^2 - a^2)t + a^2x^2= 0$$

와 같이 $t$에 대한 4차 방정식이 됩니다. 4차 방정식은 근의 공식이 있어서 이 식을 해석적으로도 풀 수 있습니다. 

 

$t = \cos \theta$ 라고 두었고, $\vec{p}$는 제 1사분면에 있기 때문에 $\cos \theta$는 양수이며, $\sin \theta$ 역시 양수가 됩니다. 따라서 위 4차 방정식의 해 중에서, 양수이면서 0보다 크거나 같고 1보다 작거나 같은 $t$에 해당하는 값이 우리가 찾는 $\cos \theta$ 값이 됩니다. 이 조건을 만족 시키는 4차 방정식의 해를 $t_0$ 이라고 놓으면, 

$$x = a t_0, y = b \sqrt{1-t_0^2}$$

가 우리가 찾는 최종 해가 됩니다. 4차 방정식의 해를 해석적으로 쓰면 너무 길어지기 때문에, 그냥 $t_0$로 놓고 더 이상 풀어쓰지 않도록 하겠습니다. 

 

이 방법은 직관적이면서 간단한 방법이지만, 공간의 차원이 높아진다면 시도하기 어려운 방법입니다. 예를들어, 같은 문제를 3차원 공간에서 푼다면, 타원을 표현하기 위한 매개변수는 2개가 필요합니다. 예를들어 구면 좌표계와 같은 변환을 사용하면, $$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$$

를 기술하는 매개 변수화는

$$x = a \sin \theta \cos \phi$$

$$y = a \sin \theta \sin \phi$$

$$z = a \cos \theta$$

이 있습니다. 이를 이용하면 $|\vec{x} - \vec{p}|$는 $\theta, \phi$의 함수가 됩니다. $\theta, \phi$에 대한 연립 방정식을 풀어야 하는데, 이 과정이 어려울 수 도 있습니다. 그나마 공간 차원이 3차원이면 그런대로 식이라도 전개할 수 있는데, 일반적으로 $n$ 차원 공간에서의 문제를 풀려고 하려면 한도 끝도 없이 식이 복잡해 지고, 결국에는 문제를 풀 수 없어 집니다.

 

2. 라그랑지 승수를 이용한 풀이

 

구속 조건이 있는 최적화 문제라면 라그랑지 승수Lagrange multiiplier를 이용하여 최적화 문제를 해결 할 수 있습니다. 위 문제에 대한 라그랑지안은

$$\mathcal{L} = |\vec{x} - \vec{p}|^2 - \lambda (\vec{x} \Sigma \vec{x} - 1)$$

입니다. 여기서 $\lambda$는 라그랑지 승수 입니다. 

$$\Sigma =  \begin{pmatrix}\frac{1}{a^2}&0\\0&\frac{1}{b^2}\end{pmatrix} $$

입니다. 따라서 $\vec{x} \Sigma \vec{x} - 1$ 를 풀어서 쓰면, 타원의 방정식 $\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$이 됩니다. 2차 형식의 경우 위와 같이 $\Sigma$ 행렬을 사용하면 간단하게 표현할 수 있습니다. 위 라그랑지안에 대한 최적화 해는

$$\frac{\partial \mathcal{L}}{\partial \vec{x}} = 2(\vec{x} - \vec{p}) - 2 \lambda \Sigma \vec{x} = 0$$

$$\frac{\partial \mathcal{L}}{\partial \lambda} = \vec{x} \Sigma \vec{x} -1 = 0$$

을 동시에 만족하는 $\vec{x}, \lambda$ 입니다. 첫 번째 식을 $\vec{x}$에 대해서 정리하면, 

$$\vec{x} = (1- \lambda \Sigma)^{-1}\vec{p}$$

가 되고, 이를 두 번째 식에 대입하여 정리하면

$$f(\lambda) = (1- \lambda \Sigma)^{-1}\vec{p} \Sigma (1- \lambda \Sigma)^{-1}\vec{p} - 1 = 0$$

이 됩니다. 위 식에서 $\Sigma, \vec{p}$는 이미 주어진 값이니, 위 식은 $\lambda$에 대한 방정식이 됩니다. 이 $\lambda$에 대한 해를 $\lambda_0$이라고 하면, 우리가 찾고자 하는 해는 $\vec{x} = (1 - \lambda \Sigma)^{-1} \vec{p}$ 가 됩니다. 위 식에서 나오는 행렬은 모두 $2 \times 2$ 행렬이기 때문에 역행렬을 해석적으로 구할 수 있습니다. 그런데 전체 식은 $\lambda$에 대한 조금은 복잡한 식이 되니, 이를 해석적으로 계산을 할 수는 있으나 조금은 어렵고 귀찮습니다. 

 

이런 경우에는 위 $f(\lambda)$의 해를 수치적인 방법으로 구할 수 있습니다. 아래는 $\lambda$를 계산하고, 최종적으로는 $\vec{x}$를 계산하는 파이썬 코드 입니다.

import numpy as np
from scipy import optimize

def function(lam):
    A = np.identity(2) - lam* sigma
    Ainv = np.linalg.inv(A)
    y = np.matmul(Ainv, x)
    f = y.transpose() @ sigma @ y-1
    return f

def return_p(lam):
    A = np.identity(2) - lam * sigma
    Ainv = np.linalg.inv(A)
    y = np.matmul(Ainv, x)
    return y


def lagrangian():
    a, b = 2, 1
    x = np.array([5, 5])

    sigma = np.zeros([2, 2])
    sigma[0, 0] = 1 / a ** 2.0
    sigma[1, 1] = 1 / b ** 2.0

    sol = optimize.brentq(function, -1e8, 0)
    p1 = return_p(sol)
    return p1

if __name__ == '__main__':
	lagrangian()

위 코드에서는 scipy의 라이브러리인 optimize.brentq 알고리듬을 활용합니다. 

 

이 방식은 다로 복잡해 보이기는 하지만, 일반적인 공간 차원에서의 문제를 풀 수 있다는 장점이 있습니다. 위에서 문제를 설명하기 위해서 $\Sigma$를 $2 \times 2$ 행렬로 두었지만, 일반적인 $n$차원 공간의 문제는 $\Sigma$가 $n \times n$ 행렬로 바뀌는 것 밖에는 차이가 없습니다. 코딩을 해 둔다면 숫자만 몇 개 더 넣으면 해결되는 문제 입니다.