[Numpy-discussion] Efficient orthogonalisation with scipy/numpy

[Gael Varoquaux]

On Tue, Jan 19, 2010 at 02:22:30PM -0600, Robert Kern wrote:
> > y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0])

> > with np = numpy and linalg = scipy.linalg where scipy calls ATLAS.

> For clarification, are you trying to find the components of the y
> vectors that are perpendicular to the space spanned by the 10
> orthonormal vectors in confounds?

Yes. Actually, what I am doing is calculating partial correlation between
x and y conditionally to confounds, with the following code:

def cond_partial_cor(y, x, confounds=[]):
    """ Returns the partial correlation of y and x, conditionning on
        confounds.
    """
    # First orthogonalise y and x relative to confounds
    if len(confounds):
        y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0])
        x = x - np.dot(confounds.T, linalg.lstsq(confounds.T, x)[0])
    return np.dot(x, y)/sqrt(np.dot(y, y)*np.dot(x, x))

I am not sure that what I am doing is optimal.

> > Most of the time is spent in linalg.lstsq. The length of the vectors is
> > 810, and there are about 10 confounds.

> Exactly what are the shapes? y.shape = (810, N); confounds.shape = (810, 10)?

Sorry, I should have been more precise:

y.shape = (810, )
confounds.shape = (10, 810)

Thanks,

Gaël