ti, 2010-01-19 kello 22:12 +0100, Gael Varoquaux kirjoitti: > On Tue, Jan 19, 2010 at 02:58:32PM -0600, Robert Kern wrote: > > > 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) > > Unfortunately, confounds is not orthonormal, and as it is different at > each call, I cannot orthogonalise it as a preprocessing.
You orthonormalize it on each call, then. It's quite likely cheaper to do than the SVD that lstsq does. {{{ import numpy as np def gram_schmid(V): """ Gram-Schmid orthonormalization of a set of `M` vectors, in-place. Parameters ---------- V : array, shape (N, M) """ # XXX: speed can be improved by using routines from scipy.lib.blas # XXX: maybe there's an orthonormalization routine in LAPACK, too, # apart from QR. too lazy to check... n = V.shape[1] for k in xrange(n): V[:,k] /= np.linalg.norm(V[:,k]) for j in xrange(k+1, n): V[:,j] -= np.vdot(V[:,j], V[:,k]) * V[:,k] return V def relative_ortho(x, V, V_is_orthonormal=False): """ Relative orthogonalization of vector `x` versus vectors in `V`. """ if not V_is_orthonormal: gram_schmid(V) for k in xrange(V.shape[1]): x -= np.vdot(x, V[:,k])*V[:,k] return x V = np.array([[1,0,1], [1,1,0]], dtype=float).T x = np.array([1,1,1], dtype=float) relative_ortho(x, V) print x print np.dot(x, V[:,0]) print np.dot(x, V[:,1]) }}} -- Pauli Virtanen _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion