On Sun, Jun 14, 2009 at 10:20 AM, Tom K.<t...@kraussfamily.org> wrote: > > > > jseabold wrote: >> >> On Mon, Jun 8, 2009 at 3:33 PM, Robert Kern<robert.k...@gmail.com> wrote: >>> On Mon, Jun 8, 2009 at 14:10, Alan G Isaac<ais...@american.edu> wrote: >>>>>> Going back to Alan Isaac's example: >>>>>> 1) beta = (X.T*X).I * X.T * Y >>>>>> 2) beta = np.dot(np.dot(la.inv(np.dot(X.T,X)),X.T),Y) >>>> >>>> >>>> Robert Kern wrote: >>>>> 4) beta = la.lstsq(X, Y)[0] >>>>> >>>>> I really hate that example. >> > > I propose the following alternative for discussion: > > U, s, Vh = np.linalg.svd(X, full_matrices=False) > 1) beta = Vh.T*np.asmatrix(np.diag(1/s))*U.T*Y > 2) beta = np.dot(Vh.T, np.dot(np.diag(1/s), np.dot(U.T, Y))) > > 1) is with X and Y starting out as matrices, 2) is with . >. Sigh. >
The problem is that I left the diagonal returned by svd as an array rather than a matrix for backward compatibility. Diag returns the diagonal when presented with a 2d matrix (array). I went back and forth on that, but not doing so would have required everyone to change their code to use diagflat instead of diag. I do note that diag doesn't preserve matrices and that looks like a bug to me. This is also an argument against over-overloaded functions such as diag. Such functions are one of the blemishes of MatLab, IMHO. OTOH, there is something of a shortage of short, meaningful names Chuck _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
