# [Numpy-discussion] Efficient orthogonalisation with scipy/numpy

Pauli Virtanen pav@iki...
Tue Jan 19 15:17:12 CST 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

```