[Numpy-discussion] Efficient orthogonalisation with scipy/numpy

[Charles R Harris]

On Tue, Jan 19, 2010 at 1:47 PM, Gael Varoquaux <
gael.varoquaux at normalesup.org> wrote:

> 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)
>
>
Column stack the bunch so that the last column is y, then do a qr
decomposition. The last column of q is the (normalized) orthogonal vector
and its amplitude is the last (bottom right) component of r.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20100119/7b87e2d9/attachment.html>