[Numpy-discussion] Cross-covariance function

[Sturla Molden]

Den 20.01.2012 13:39, skrev Pierre Haessig:
> I don't see how does your function relates to numpy.cov [1]. Is it an
> "extended case" function or is there a difference in the underlying math ?
>

If X is rank n x p, then np.cov(X, rowvar=False) is equal
to S after

     cX = X - X.mean(axis=0)[np.newaxis,:]
     S = np.dot(cX.T, cX)/(n-1.)

If we also have Y rank n x p, then the upper p x p
quadrant of

     np.cov(X, y=Y, rowvar=False)

is equal to S after

     XY = np.hstack(X,Y)
     cXY = XY - XY.mean(axis=0)[np.newaxis,:]
S = np.dot(cXY.T, cXY)/(n-1.)

Thus we can see thatthe total cocariance is composed
of four parts:

     S[:p,:p]  = np.dot(cX.T, cX)/(n-1.)   # upper left
     S[:p,p:]  = np.dot(cXY.T, cYY)/(n-1.) # upper right
     S[p:,:p]  = np.dot(cY.T, cX)/(n-1.)   # lower left
     S[p:,:p]  = np.dot(cYX.T, cYX)/(n-1.) # lower right

Often we just want the upper-right p x p quadrant. Thus
we can save 75 % of the cpu time by not computing the rest.

Sturla