What does $\mbox{cholesky}(P)^{-1} [\cos(\theta), \sin(\theta)]$ correspond to?

25 Views Asked by At

I'm trying to understand a function, and its name isn't really helping me (I've found resources similar to this and this but this is either R or not exactly "ELLIPLOT" and I'm not sure I fully understand the concept of an elliplot).

I know what a Cholesky decomposition is but I can't understand the sort of link it has with ellipses. I guess I'm looking for resources or explanation as to what pan of mathematics this is related to. I've encountered this function in a Linear Programming course but we haven't done much more than writing it and I'm a bit confused as to why we did.

def ellipplot(P,gamma=None,xc=None):
    """Matlab ELLIPPLOT Internal function for plotting ellipsoid"""
    if xc is None:
        xc = np.zeros(len(P))
    if gamma is None:
        gamma = 1 
    P=P/gamma
    P=np.linalg.cholesky(P)
    theta = np.linspace(-3.14, 3.14, 1000)
    z = np.array([np.cos(theta), np.sin(theta)])
    x = np.linalg.inv(P)@z
    for n in range(x.shape[1]):
        x[:,n]=x[:,n]+xc
        xe = x[0,:]
        ye = x[1,:]
    return xe, ye

xe, ye = ellipplot(
    np.array(
        [[2, 1], 
         [1, 2]]
    )
)
plt.plot(xe, ye)
plt.show()

This is the matlab version of the function, I believe.

1

There are 1 best solutions below

5
On BEST ANSWER

The ellipse is defined by $$ p_{11}x_1^2+p_{22}x_2^2+p_{12}x_1x_2+p_{21}x_1x_2=1\,. $$ In matrix/vector notation this is $$\tag{1} x^\top Px=1 $$ where $$ P=\left(\begin{matrix}p_{11}&p_{12}\\p_{21}&p_{22}\end{matrix}\right)\,,~~x=\left(\begin{matrix}x_1\\x_2\end{matrix}\right)\,. $$ To plot the points that satisfy (1) the code takes the Cholesky decomposition (square root) of $P=LL^\top$ so that (1) can be written as $$\tag{2} y^\top y=1\,~~\text{ where }y=L^\top x\,. $$ Clearly, $$ y=\left(\begin{matrix}\cos\theta\\\sin\theta\end{matrix}\right)\,\quad\quad\theta\in[-\pi,\pi] $$ satisfy (2). Therefore, to find the vector $x$ we onyly have to invert $y=L^\top x$ and that's what the code does.

Strictly speaking the code inverts $L$ (not $L^\top$) but I checked that it plots an ellipse.