normalize Does Not Need An Estimated Zero

[dmbaker]

This code will work:

def normalize(v, axis=-1, order=2):
    """
    Normalize a vector for arbitrary axes, and giving optimal performance.
    """
    one_dim = np.ndim(v) == 0 or np.ndim(v) == 1
    l2 = np.atleast_1d(np.linalg.norm(v, order, axis))
    l2[l2==0] = 1
    ret = a / np.expand_dims(l2, axis)
    return ret[0] if one_dim else ret

[jokva]

Do you have some example input that would trigger the instability that can be used to build a test case?

[dmbaker]

@jokva who are you asking this question of?

[jokva]

You, dmbaker. Also applies to the other issues, really, with angles and numerical instability.

[dmbaker]

This one is not in not really about the stability. More that you don't need to use a tolerance.

I missed a variable name change. Here is the fixed version.

def normalized(a, axis=-1, order=2):
    """
    Normalize a vector for arbitrary axes, and giving optimal performance.
    """
    one_dim = np.ndim(a) == 0 or np.ndim(a) == 1
    l2 = np.atleast_1d(np.linalg.norm(a, order, axis))
    l2[l2==0] = 1
    ret = a / np.expand_dims(l2, axis)
    return ret[0] if one_dim else ret

[dmbaker]

Here is a link to a very good document on calculating angles with numeric stability. Read page 15 for why you should not use the inverse cosine.

Computing Cross-Products and Rotations in 2- and 3-Dimensional Euclidean Spaces

Also see this one of Kahan's papers starting a page 45

How Futile are Mindless Assessments of Roundoff in Floating-Point Computation ?

From Wikipedia:
Kahan was the primary architect behind the IEEE 754-1985 standard for floating-point computation (and its radix-independent follow-on, IEEE 854). He has been called "The Father of Floating Point", since he was instrumental in creating the original IEEE 754 specification.[1] Kahan continued his contributions to the IEEE 754 revision that led to the current IEEE 754 standard.

[jokva]

I don't see why your normalized() function would perform better as it's doing the same calls, with the exception of you explicitly using the ord = 2 argument. Interestingly, with the ord = 2 argument (which I though was the default for vectors?), the close-to-zero issue seems to without having to filter out close-to-zero values, which is nice.

edit: I guess it interprets the list-of-vectors as a matrix and uses the Frobenius norm instead. Using explicit ord = 2 seems to fix this issue.

[jokva]

#44 The fix

[dmbaker]

@jokva my normalized() is a general function that can handle norms of vectors in dimensions other than 3, norms other than 2 (Euclidian) and the axis along which to calculate the norm. The default handling is the Euclidian norm of 3D vectors.

My initial comment was that you probably don't need the epsilon test.

[dmbaker]

Also, I forgot to give you this link on calculating angles:

On the Angle Between Vectors

[jokva]

My initial comment was that you probably don't need the epsilon test.

Both approaches uses a assign-1-where-zero, but the eps was there because (fuzzed) testing showed that the function was sensitive to sufficiently small values. This goes away when using ord = 2, which is nicer anyway since this is always for 2D vectors and not ND arrays.