[Numpy-discussion] Efficient orthogonalisation with scipy/numpy

[Pauli Virtanen]

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