Repository: pytorch
Issue: Add numerically-stable function for angle between vectors
Labels: module: numerical-stability,triaged,enhancement,module: distance functions

## 🚀 Feature
<!-- A clear and concise description of the feature proposal -->

Add a numerically-stable function that returns the angle between two vectors.

## Motivation

<!-- Please outline the motivation for the proposal. Is your feature request related to a problem? e.g., I'm always frustrated when [...]. If this is related to another GitHub issue, please link here too -->

The most well-known formula for computing the angle between two vectors is
```python3
def angle_v0(a, b, dim=-1):
    c = (a * b).sum(dim) / (a.norm(dim=dim) * b.norm(dim=dim))
    return c.arccos()
```

Unfortunately, this formula is *extremely* numerically unstable when `a` and `b` are close to parallel or antiparallel. In fact, as `c` approaches ±1, the gradient of `c.arccos()` explodes to -∞. This is only an apparent pathology, however, since `∠(a, b)` itself is actually well-defined and well-behaved everywhere besides the origin.

## Alternatives

<!-- A clear and concise description of any alternative solutions or features you've considered, if any. -->

A clunky workaround is to clip the argument of `arccos`, such as
```python3
def angle_v1(a, b, dim=-1, eps=1e-4):
    z = (a * b).sum(dim) / (a.norm(dim=dim) * b.norm(dim=dim))
    return z.clip(-1 + eps, 1 - eps).arccos()
```
However:
1. If `eps` is large, this compromises the correctness of the result.
2. If `eps` is small, this is still numerically unstable.

## Pitch

<!-- A clear and concise description of what you want to happen. -->

See [William Kahan](https://en.wikipedia.org/wiki/William_Kahan)'s paper [*How Futile are Mindless Assessments of Roundoff in Floating-Point Computation? §12: Mangled Angles*](https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf#page=46):

> We shall consider now a simpler example chosen to support the thesis that numerical software is extremely difficult to debug.  Evidence for this thesis is the longevity of inaccurate software in use by vast numbers of numerically active and, in most instances, sophisticated users of MATLAB. Once again we consider the angle between two directions x and y specified now in Cartesian coordinates
> ...
> Here is a better formula less well known than it deserves: ∠(x, y) := 2 arctan( ||x·||y|| – ||x||·y||/||x·||y|| + ||x||·y|| ). When executed in arithmetic rounded to n significant bits its absolute error never much exceeds 1/2ⁿ⁻¹ unless the dimension is gargantuan. And redirected roundoff barely affects this formula.
> ...
> What do mangled angles teach us? The goal of error-analysis is not to find errors but to fix them.  They have to be found first. The embarrassing longevity, over three decades, of inaccurate and/or ugly programs to compute a function so widely used as ∠(X, Y) says something bleak about the difficulty of floating-point error-analysis for experts and nonexperts: Without adequate aids like redirected roundings, diagnosis and cure are becoming practically impossible. 

Here is this formula implemented in PyTorch:

```python3
def angle_v2(a, b, dim=-1):
    a_norm = a.norm(dim=dim, keepdim=True)
    b_norm = b.norm(dim=dim, keepdim=True)
    return 2 * torch.atan2(
        (a * b_norm - a_norm * b).norm(dim=dim),
        (a * b_norm + a_norm * b).norm(dim=dim)
    )
```

The proposal is to add this to the PyTorch standard library so that users can benefit from access to a numerically-stable version of this fundamental geometric function, without having to deal with the issues that arise from the naive formula.

## Additional context

<!-- Add any other context or screenshots about the feature request here. -->

See also the Math StackExchange question [*Numerically stable method for angle between 3D vectors*](https://math.stackexchange.com/questions/1143354/numerically-stable-method-for-angle-between-3d-vectors).