[SciPy-user] how to find projection of a point to subspace?
Robert Kern
robert.kern@gmail....
Sat May 30 15:36:22 CDT 2009
On Fri, May 29, 2009 at 13:51, Dmitrey <dmitrey15@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
More information about the SciPy-user
mailing list