On 9/22/09, Fred lunnon <fred.lunnon@gmail.com> wrote:
... PSLQ is implemented in both Mathematica and Maple. I compared the Maple version with my home-grown geometrically-inspired "Umbrella" algorithm, only to discover that PSLQ detected the polynomial satisfied by a quadratic surd, using only half the precision required by Umbrella --- Umbrella retires, blown inside-out.
All these algorithms have an inherent limitation: they detect only homogeneous relations (i.e. with RHS = 0) with integer coefficients, between a given bunch of reals; whereas the Lego orrery appears to demand a RHS equal to the gear ratio required to be approximated. I'm still thinking about this --- if only as a last-ditch attempt to salvage a little pride.
Mature reflection suggests that GAMP [as I'm thinking of re-christening the umbrella algorithm, just to keep up with MDCF, PSOS, HJLS, PSLQ, ...] only needed twice as many digits as PSLQ in tests, because I hadn't bothered to evaluate inner products with double precision, as any serious implementation should. Casual inspection might suggest that the inhomogeneous problem (the RHS might just as well be unity) should only involve some fairly trivial modification to the homogeneous problem. It seems to me that geometry gives some insight into this situation, so I'll take the liberty of describing GAMP briefly. The input target vector of real components to be rationally approximated is interpreted as a line in projective n-space joining the origin to the point at infinity with those coordinates. A `ratio' basis of n integer lattice points, regarded as vectors rooted at the origin, is initially set equal to the unit points on the Cartesian axes; as the algorithm progresses, its vectors fold ever more tightly around the line, like an umbrella closing around its shaft, in such a way that the volume of the enclosed parallelepiped remains unity, and their convex hull (as it were) continues to contains the line. Each member of the basis constitutes a rational approximation to the target vector; furthermore, since (Minkowski?) there can be no lattice points within the parallelepiped, the point at which the line emerges from it gives a lower bound on the length of any closer lattice point. At each iteration, the next vector in cyclic order is selected as candidate for improvement; its components with respect to the (non-Cartesian) basis comprising the target point, the origin, and the other n-1 vectors, are computed by inverting the matrix of that basis; then they are truncated to integers, and the corresponding linear combination subtracted from the candidate. Meanwhile the corresponding `relation' basis, the (Hodge) dual and matrix inverse of the ratio basis, yields progressively better approximate null linear combinations of the target reals. Termination is in practice easy to detect: either the candidate fails to improve, when an exact solution, or (nearly) best approximations given the precision, have been attained; or else it improves and lengthens wildly, when a relation will have been detected. Convergence takes place within at most d*log(10)/log(x) iterations, where d equals the decimal precision of the input, and x^n = x^(n-1) + ... + 1, with x ~ 2 - 1/2^n. Now then, an integer relation has a straightforward geometrical interpretation, as a prime flat (hyperplane) with integer coefficients, passing close to or meeting the target point; if the relation happens to be inhomogeneous, the prime simply fails to meet the origin as well. But there is then no apparent geometrical equivalent associated with its corresponding dual: rational approximate ratios are meaningful only for the homogeneous problem! This suggests that no obvious strategy, such as recasting GAMP to operate directly with the relation basis instead, can possibly succeed. I'm not at this stage sufficiently familiar with the other algorithms to comment on them; but I begin strongly to suspect that this difficulty is inherent. Perhaps the intelligent thing to do now would be to stop blathering on, and go and consult Ferguson ... Fred Lunnon