On 03/02/2021 23:23, Henry Baker wrote:
Input: an nxn(*) matrix M of floating point numbers (positive, negative or zero) which is *supposed* to have rank exactly one, ... Unfortunately, due to floating point approximations and and other issues, the values of M have small errors,
and asks: How can you find the vectors whose outer product M approximately is? The obvious approach is to remark that if M = uv'+h (h small) then Mb = (v.b)u + hb which is approximately a multiple of u, and M'a = (u.a)v + h'a which is approximately a multiple of v. So starting with any vector b not perpendicular to v, we can iterate a := (Mb)/b^2 b := (M'a)/a^2. (Why the b^2 and a^2 scaling factors? See the end of the next paragraph.) Before asking how fast this process converges (on the face of it is seems obviously first-order) we should ask what we expect it to converge _to_. Clearly not, in general, the true u,v, but we might hope for convergence to the u,v for which uv' is as close as possible to M. If we take this in the least-squares sense then the resulting derivative=0 equations say that u=Mv/v^2 and v=M'u/u^2; that is, that hv=0 and h'u=0. We still don't expect convergence to exactly u,v because we could equally converge to ku,v/k for some scalar k. But suppose that, given our current a,b, we pick the scaling of u,v to make a=u+d with d as small as possible. In particular, this implies u.d=0. Then our new approximation to v is M'(u+d)/(u+d)^2 = (vu'+h')(u+d)/(u+d)^2 = (u^2 v + h'd)/(u^2+d^2) where the various cross-terms vanish because u.d=0 and h'u=0. So to first order in d, this is v + h'd / u^2 = v + e, say. A little calculation shows that |e|/|v| is at most R/|u||v| times |d|/|u|, where R is the 2-norm of h; of course |u||v| is the 2-norm of uv'. So our errors go down by at worst a factor of about |h|/|uv'| at each iteration. If the errors in M are as small as the words "floating point approximations and other issues" suggest, it shouldn't take more than a few iterations for this to converge to machine accuracy. (Again, convergence will be not to the true u,v but to whichever u,v actually yield the best approximation to the perturbed M.) -- g