[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.

    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:
    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