[Numpy-discussion] Efficient orthogonalisation with scipy/numpy
Pauli Virtanen
pav at iki.fi
Tue Jan 19 16:17:12 EST 2010
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
More information about the NumPy-Discussion
mailing list