[SciPy-user] how to find projection of a point to subspace?

[Robert Kern]

On Fri, May 29, 2009 at 13:51, Dmitrey <dmitrey15 at ukr.net> wrote:
> hi all,
>
> Suppose I have point x = [x0, ..., xn-1] and linear subspace defined by
> Ax=b, A is m x n matrix, b is vector of length m.

For clarity, I will use "y" to refer to the vector you want to project
and "x" to refer to any vector in the subspace defined by "dot(A, x) =
b". And I assume that m < n.

> What is the best way to find projection of the point x to the
> subspace? (Suitable for ill-conditioned cases)

The subspace is parallel to the null space of A. You just need to
translate it by a vector, which can be any solution to dot(A,x)=b.
np.linalg.lstsq(A,b) gives you the minimum-norm x that satisfies this
equation, let's call it x0. To find the null space, use the SVD;
namely the vectors will be Vh[m:] where Vh is the third output from
np.linalg.svd(A). Then the projection of y onto your subspace is

  x0 + (Vh[i:] * dot(Vh[i:], y-x0)[:,np.newaxis]).sum(axis=0)

What this does is make x0 the origin temporarily; compute the
decomposition of the projection onto the null space orthonormal frame;
assembles the decomposition back into the original coordinates; then
restores the origin.

For a well-conditioned A, i==m. For ill-conditioned A, find the first
i such that s[i]/s[0] < epsilon, where s is the vector of singular
values (i.e. the second output from np.linalg.svd(A)). Also, use
np.linalg.lstsq(A, b, rcond=s[i]) in order to find x0 in the
ill-conditioned case.


import numpy as np


def project_subspace(y, A, b, eps=np.finfo(float).eps):
    """ Project a vector onto the subspace defined by "dot(A,x) = b".
    """
    m, n = A.shape
    u, s, vh = np.linalg.svd(A)
    # Find the first singular value to drop below the cutoff.
    bad = (s < s[0] * eps)
    i = bad.searchsorted(1)
    if i < m:
        rcond = s[i]
    else:
        rcond = -1
    x0 = np.linalg.lstsq(A, b, rcond=rcond)[0]
    null_space = vh[i:]
    y_proj = x0 + (null_space * np.dot(null_space,
y-x0)[:,np.newaxis]).sum(axis=0)
    return y_proj

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco