<br><br><div class="gmail_quote">On Tue, Jan 19, 2010 at 8:02 PM,  <span dir="ltr"><<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
On Tue, Jan 19, 2010 at 9:47 PM, Charles R Harris<br>
<div><div></div><div class="h5"><<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br>
><br>
><br>
> On Tue, Jan 19, 2010 at 6:08 PM, <<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>> wrote:<br>
>><br>
>> On Tue, Jan 19, 2010 at 6:48 PM, Charles R Harris<br>
>> <<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br>
>> ><br>
>> ><br>
>> > On Tue, Jan 19, 2010 at 2:34 PM, <<a href="mailto:josef.pktd@gmail.com">josef.pktd@gmail.com</a>> wrote:<br>
>> >><br>
>> >> On Tue, Jan 19, 2010 at 4:29 PM, Charles R Harris<br>
>> >> <<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br>
>> >> ><br>
>> >> ><br>
>> >> > On Tue, Jan 19, 2010 at 1:47 PM, Gael Varoquaux<br>
>> >> > <<a href="mailto:gael.varoquaux@normalesup.org">gael.varoquaux@normalesup.org</a>> wrote:<br>
>> >> >><br>
>> >> >> On Tue, Jan 19, 2010 at 02:22:30PM -0600, Robert Kern wrote:<br>
>> >> >> > > y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0])<br>
>> >> >><br>
>> >> >> > > with np = numpy and linalg = scipy.linalg where scipy calls<br>
>> >> >> > > ATLAS.<br>
>> >> >><br>
>> >> >> > For clarification, are you trying to find the components of the y<br>
>> >> >> > vectors that are perpendicular to the space spanned by the 10<br>
>> >> >> > orthonormal vectors in confounds?<br>
>> >> >><br>
>> >> >> Yes. Actually, what I am doing is calculating partial correlation<br>
>> >> >> between<br>
>> >> >> x and y conditionally to confounds, with the following code:<br>
>> >> >><br>
>> >> >> def cond_partial_cor(y, x, confounds=[]):<br>
>> >> >>    """ Returns the partial correlation of y and x, conditionning on<br>
>> >> >>        confounds.<br>
>> >> >>    """<br>
>> >> >>    # First orthogonalise y and x relative to confounds<br>
>> >> >>    if len(confounds):<br>
>> >> >>        y = y - np.dot(confounds.T, linalg.lstsq(confounds.T, y)[0])<br>
>> >> >>        x = x - np.dot(confounds.T, linalg.lstsq(confounds.T, x)[0])<br>
>> >> >>    return np.dot(x, y)/sqrt(np.dot(y, y)*np.dot(x, x))<br>
>> >> >><br>
>> >> >> I am not sure that what I am doing is optimal.<br>
>> >> >><br>
>> >> >> > > Most of the time is spent in linalg.lstsq. The length of the<br>
>> >> >> > > vectors<br>
>> >> >> > > is<br>
>> >> >> > > 810, and there are about 10 confounds.<br>
>> >> >><br>
>> >> >> > Exactly what are the shapes? y.shape = (810, N); confounds.shape =<br>
>> >> >> > (810,<br>
>> >> >> > 10)?<br>
>> >> >><br>
>> >> >> Sorry, I should have been more precise:<br>
>> >> >><br>
>> >> >> y.shape = (810, )<br>
>> >> >> confounds.shape = (10, 810)<br>
>> >> >><br>
>> >> ><br>
>> >> > Column stack the bunch so that the last column is y, then do a qr<br>
>> >> > decomposition. The last column of q is the (normalized) orthogonal<br>
>> >> > vector<br>
>> >> > and its amplitude is the last (bottom right) component of r.<br>
>> >><br>
>> >> do you have to do qr twice, once with x and once with y in the last<br>
>> >> column or can this be combined?<br>
>> >><br>
>> >> I was trying to do something similar for partial autocorrelation for<br>
>> >> timeseries but didn't manage or try anything better than repeated<br>
>> >> leastsq or a variant.<br>
>> >><br>
>> ><br>
>> > Depends on what you want to do. The QR decomposition is essentially<br>
>> > Gram-Schmidt on the columns. So if you just want an orthonormal basis<br>
>> > for<br>
>> > the subspace spanned by a bunch of columns, the columns of Q are they.<br>
>> > To<br>
>> > get the part of y orthogonal to that subspace you can do y - Q*Q.T*y,<br>
>> > which<br>
>> > is probably what you want if the x's are fixed and the y's vary. If<br>
>> > there is<br>
>> > just one y, then putting it as the last column lets the QR algorithm do<br>
>> > that<br>
>> > last bit of projection.<br>
>><br>
>> Gram-Schmidt (looking at it for the first time) looks a lot like<br>
>> sequential least squares projection. So, I'm trying to figure out if I<br>
>> can use the partial results up to a specific column as partial least<br>
>> squares and then work my way to the end by including/looking at more<br>
>> columns.<br>
>><br>
><br>
> I don't the QR factorization would work for normal PLS. IIRC, one of the<br>
> algorithms does a svd of the cross correlation matrix. The difference is<br>
> that in some sense the svd picks out the best linear combination of columns,<br>
> while the qr factorization without column pivoting just takes them in order.<br>
> The QR factorization used to be the method of choice for least squares<br>
> because it is straight forward to compute, no iterating needed as in svd,<br>
> but these days that advantage is pretty much gone. It is still a common<br>
> first step in the svd, however. The matrix is factored to Q*R, then the svd<br>
> of R is computed.<br>
<br>
</div></div>I (finally) figured out svd and eigenvalue decomposition for this purpose.<br>
<br>
But from your description of QR, I thought specifically of the case<br>
where we have a "natural" ordering of the regressors, similar to the<br>
polynomial case of you and Anne. In the timeseries case it would be by<br>
increasing lags<br>
<br>
yt on y_{t-1}<br>
yt on y_{t-1}, y_{t-2}<br>
...<br>
...<br>
yt on y_{t-k} for k= 1,...,K<br>
<br>
or yt on xt and the lags of xt<br>
<br>
This is really sequential LS with a predefined sequence, not PLS or<br>
PCA/PCR or similar orthogonalization by "importance".<br>
The usual procedure for deciding on the appropriate number of lags<br>
usually loops over OLS with increasing number of regressors.<br>
>From the discussion, I thought there might be a way to "cheat" in this<br>
using QR and Gram-Schmidt<br>
<br></blockquote><div><br>Ah, then I think your idea would work. The norms of the residuals at each step would be along the diagonal of the R matrix. They won't necessarily decrease monotonically, however.<br><br>Chuck  <br>
</div><br></div>