CODEWITHSUNDEEP
pythonalgorithmnumpy
Jared Goguen
Jared Goguen

Reputation: 9010

Denormalize unit vector

Suppose I have the unit vector, u.

from numpy import mat
u = mat([[0.9**0.5], [-(0.1)**0.5]])
# [ 0.9486833  -0.31622777]

The unit vector is an eigenvector of a matrix with integer entries. I also know that the eigenvalues are integers. So, the unit vector will contain irrational decimals that, when squared, are decimal approximations of rational numbers.

Is there any good way to find the smallest value k such that all entries of ku are integers? In general, k will be the square root of an integer.

A naive approach would be to square each entry in the vector and find the smallest ki that produces an integer. Then, k will be the square root the of LCM of all ki. I'm hoping there is a better approach than this.

Note that this is not a duplicate of this question.

Upvotes: 9

Views: 1484

Answers (4)

Jared Goguen
Jared Goguen

Reputation: 9010

I've improved the methodology provided by Christian in order to accept a wider range of fractions. The trick was to "pre-normalize" the unit vector by dividing it by the smallest non-zero entry. This works reliably for all fractions of a specified maximum denominator.

from fractions import Fraction, gcd

def recover_integer_vector(u, denom=50):
    '''
    For a given vector u, return the smallest vector with all integer entries.
    '''

    # make smallest non-zero entry 1
    u /= min(abs(x) for x in u if x)

    # get the denominators of the fractions
    denoms = (Fraction(x).limit_denominator(denom).denominator for x in u)

    # multiply the scaled u by LCM(denominators)
    lcm = lambda a, b: a * b / gcd(a, b)
    return u * reduce(lcm, denoms)

Testing:

The following code was used to ensure that all fractions in the given range work correctly.

import numpy as np

from numpy import array
from itertools import combinations_with_replacement


for t in combinations_with_replacement(range(1, 50), 3):
    if reduce(gcd, t) != 1: continue

    v = array(map(float, t))
    u = v / np.linalg.norm(v)

    w = recover_integer_vector(u)

    assert np.allclose(w, v) or np.allclose(w, -v)

As can be seen by the testing time, this code is not very quick. It took my computer about 6 seconds to test. I'm not sure whether the timing can be improved.

Upvotes: 4

B. M.
B. M.

Reputation: 18638

Let u=[cos(α),sin(α)] the unit vector of v=[p,q], where p and q are little integers, prime together . Then tan(α) is a rational.

So there is a straightforward way to find back v from u : fractions.Fraction(u[1]/u[0]).limit_denominator()

Example :

v=array([-1,3])
u=v/np.linalg.norm(v)
# array([-0.31622777,  0.9486833 ])
f=Fraction(u[1]/u[0]).limit_denominator()
# Fraction(-1, 3)

if you want k, norm([f.numerator,f.denominator]) give it.

p=0 is a trivial particular case.

Upvotes: 2

Christian
Christian

Reputation: 729

Ok, here's my best stab at it. This algorithm seems to work, though I can't guarantee it will be faster than what you originally proposed:

import numpy as np
from fractions import Fraction, gcd

# need some finesse and knowledge of problem domain to fine tune this parameter
# otherwise you'll end up with extremely large k due to precision problems
MAX_DENOMINATOR = 1000

# least common multiple
def lcm(a, b):
    return a * b / gcd(a, b)

# returns the denominator of the square of a number
def get_square_denominator(x):
    return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator

# returns the smallest multiplier, k, that can be used to scale a vector to
# have integer coefficients. Assumes vector has the property that it can be 
# expressed as [x**0.5, y**0.5, ...] where x, y, ... are rational with 
# denominators <= MAX_DENOMINATOR
def get_k(v):
    denominators = [get_square_denominator(i.item()) for i in v]
    return reduce(lcm, denominators) ** 0.5

Here it is in use:

>>> u = matrix([[0.9486833], [-0.31622777]])
>>> get_k(u)
3.1622776601683795
>>> get_k(u) * u
matrix([[ 3.00000001],
        [-1.00000001]])

Upvotes: 2

hobbs
hobbs

Reputation: 239920

You could square the terms, find the sequence of continued fraction convergents of each, replace the numerator of each convergent with the nearest square number, and look for the smallest denominator that appears in all of the lists of convergents (as part of convergents with tolerable accuracy). Your definition of "tolerable accuracy" allows the algorithm to deal with a bit of rounding error without producing an insanely large result. In the case of your example, 9/10 and 1/10 would be among the first few convergents, with relative errors less than 10^-9, and sqrt(10) would be your scaling factor (multiplying out to [3, -1] which you can read directly off of the square roots of the numerators of the convergents, if you make sure to recover the signs). It's unclear what to do if this algorithm doesn't find a common denominator, though.

Upvotes: 2

Related Questions