[SciPy-user] how to find projection of a point to subspace?
Sun May 31 03:12:53 CDT 2009
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.
Robert Kern wrote:
> On Fri, May 29, 2009 at 13:51, Dmitrey <email@example.com> 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 < 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 * eps)
> i = bad.searchsorted(1)
> if i < m:
> rcond = s[i]
> rcond = -1
> x0 = np.linalg.lstsq(A, b, rcond=rcond)
> null_space = vh[i:]
> y_proj = x0 + (null_space * np.dot(null_space,
> return y_proj
More information about the SciPy-user