[Numpy-discussion] Efficient orthogonalisation with scipy/numpy

[Robert Kern]

On Tue, Jan 19, 2010 at 14:47, Gael Varoquaux
<gael.varoquaux@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.

If confounds is orthonormal, then there is no need to use lstsq().

  y = y - np.dot(np.dot(confounds, y), confounds)

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco