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

[Dmitrey]

Thank you,
Stéfan van der Walt  and Andrew York have sent me a letter with that 
article mentioned, where iterative algorithm is proposed and declared as 
more efficient than direct approach.
D.

Robert Kern wrote:
> 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
>
>