plotting the eigenvector of covariance matrix using matplotlib and np.linalg

Question

I am trying to draw eigenvector and of covariance matrix received from a bunch of points (polyhedron in 3D). Here is what i do.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import linalg as la
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d

class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)
##################################################################################################
#here i start with drawing the actual polyhedron and the vector 
##################################################################################################
#generate num random points in 3d
num = 5
#coord = 10*np.random.rand(3,num)#num points in 3D #first axis is x, second = y, third = z
#xcod = np.array([1,2,3,2.7,2.4,1])
xcod = np.array([1,1,1,1,1,1])
ycod = np.array([1,1,4.5,5.,6,1])
zcod = np.array([1,-2,0,2,3,1])
#coord = np.concatenate(coord,coord[0])
#####plotting in 3d
fig = plt.figure()
ax = fig.add_subplot(111,projection = '3d')
#plotting all the points
ax.plot(xcod,ycod,zcod,'x-')
#adding labels for vertice
for i in range(num):
    ax.text(xcod[i],ycod[i],zcod[i],'%d(%.2f,%.2f,%.2f)'%(i,xcod[i],ycod[i],zcod[i]))
#supposed centroid
centroid = np.array([np.mean(xcod),np.mean(ycod),np.mean(zcod)])
ax.scatter(centroid[0],centroid[1],centroid[2],marker = 'o',color='r')
#labelling the axes
ax.set_xlabel("x axis")
ax.set_ylabel("y axis")
ax.set_zlabel("z axis")
#getting a stack of all vertices, while removing last repeat vertex
cod = np.vstack((np.delete(xcod,-1),np.delete(ycod,-1),np.delete(zcod,-1)))
#caculating covariance matrix
#ddof = 0 is using simple averages or normalising with N ; ddof = 1 means normalising with N-1
covmat = np.cov(cod,ddof=0)
#computing eigen values and eigen vectors
eigval,eigvec = la.eig(covmat)
#multiplying eigen value and eigen vec
#for counter in range(len(eigval)):
#    eigvec[counter]= eigval[counter]*eigvec[counter]
#####################################################################################
#plotting Eigen vectors
#####################################################################################
for vec in eigvec:#fetching one vector from list of eigvecs
    #drawing the vec, basically drawing a arrow form centroid to the end point of vec
    drawvec = Arrow3D([centroid[0],vec[0]],[centroid[1],vec[1]],[centroid[2],vec[2]],mutation_scale=20,lw=3,arrowstyle="-|>",color='r')
    #adding the arrow to the plot
    ax.add_artist(drawvec)
#plot show
plt.show()  

The plot I get by doing this is less than satisfactory. The view of the eigenvectors from two different angles.

I was expecting something like this. Vectors popping form the centroid to give direction of largest variance. But it seems it is not working, perhaps the eigenvectors are not correctly calculated by np.linalg?
Can you suggest me what am i missing?

Also, I am trying to draw the ellipsoid after I have the eigenvectors. If you can suggest me on that too, it would be great :)

Edit: bit of progress I think the np.linalg is just giving me position vector of eigenvectors from origin, so i just transformed them to be from centroid,

#getting tuples of x,y,z
verts = [zip(xcod,ycod,zcod)]
#plotting polyhedron surface
ax.add_collection3d(Poly3DCollection(verts,alpha=0.5))
#changing eigvec from origin to centroid
for counteri in range(len(eigvec)):
    eigvec[counteri][0]+=centroid[0]
    eigvec[counteri][1]+=centroid[1]
    eigvec[counteri][2]+=centroid[2]

adding above code before going into the plotting part of #plotting Eigen vectors

This now gives me something like below

Answer 1

The eigenvectors in eigvec are column vectors. Therefore, to retrieve the eigenvectors through iteration, you need to transpose eigvec:

for vec in eigvec.T:  

Coupling that with your observation that vec needs to be shifted by centroid:

vec += centroid

yields

For completeness,

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import linalg as LA
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d


class Arrow3D(FancyArrowPatch):

    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        FancyArrowPatch.draw(self, renderer)
##########################################################################
# here i start with drawing the actual polyhedron and the vector
##########################################################################
# generate num random points in 3d
num = 5
# coord = 10*np.random.rand(3,num)#num points in 3D #first axis is x, second = y, third = z
#xcod = np.array([1,2,3,2.7,2.4,1])
xcod = np.array([1, 1, 1, 1, 1, 1])
ycod = np.array([1, 1, 4.5, 5., 6, 1])
zcod = np.array([1, -2, 0, 2, 3, 1])
#coord = np.concatenate(coord,coord[0])
# plotting in 3d
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# plotting all the points
ax.plot(xcod, ycod, zcod, 'x-')
# adding labels for vertice
for i in range(num):
    ax.text(xcod[i], ycod[i], zcod[i], '%d(%.2f,%.2f,%.2f)' %
            (i, xcod[i], ycod[i], zcod[i]))
# supposed centroid
centroid = np.array([np.mean(xcod), np.mean(ycod), np.mean(zcod)])
ax.scatter(centroid[0], centroid[1], centroid[2], marker='o', color='r')
# labelling the axes
ax.set_xlabel("x axis")
ax.set_ylabel("y axis")
ax.set_zlabel("z axis")
# getting a stack of all vertices, while removing last repeat vertex
cod = np.vstack(
    (np.delete(xcod, -1), np.delete(ycod, -1), np.delete(zcod, -1)))
# caculating covariance matrix
# ddof = 0 is using simple averages or normalising with N ; ddof = 1 means
# normalising with N-1
covmat = np.cov(cod, ddof=0)
# computing eigen values and eigen vectors
eigval, eigvec = LA.eig(covmat)
# multiplying eigen value and eigen vec
# for counter in range(len(eigval)):
#    eigvec[counter]= eigval[counter]*eigvec[counter]
##########################################################################
# plotting Eigen vectors
##########################################################################
for vec in eigvec.T:  # fetching one vector from list of eigvecs
    # drawing the vec, basically drawing a arrow form centroid to the end
    # point of vec
    vec += centroid
    drawvec = Arrow3D([centroid[0], vec[0]], [centroid[1], vec[1]], [centroid[2], vec[2]],
                      mutation_scale=20, lw=3, arrowstyle="-|>", color='r')
    # adding the arrow to the plot
    ax.add_artist(drawvec)
# plot show
plt.show()