While wandering on the internet, I stumbled upon Paul Hsieh's blog-post, where he demonstrates a way to approximate the norm of a vector without any call to the sqrt function. Let's see if I can reproduce the steps to derive this.

---

Table of contents

---

Setting-up the scene.

Calculating the norm of a vector $(x,y)$, or a complex number $x+iy$ means calculating $\sqrt{x^2+y^2}$. Without loss of generality, we can set $\sqrt{x^2+y^2}=1$. If we draw this, we get the following.

Finding a lower bound to the norm.

Now, the issue with the norm is that the $\sqrt{}$ operation is expensive to compute. That's why we would like another way to approximate the norm. A first idea is to look at other norms available, indeed, what we have called "norm" so far is actually the 2-norm, also named euclidean norm. Let's have a look at two other norms : the infinity norm and the Manhattan norm.

Infinity norm is :

$$\lVert(x,y)\rVert_\infty = \max(x,y)$$

Manhattan norm is :

$$\lVert(x,y)\rVert_1 = |x|+|y|$$

Now we see the Manhattan norm is indeed a lower bound for the 2-norm, even if it's rough. The Infinity norm, however, is too high. But that is not an issue, we could simply scale it up so that it is always higher than the 2-norm. The scaling factor is chosen, such as the yellow curve tangent to the circle. For that, we need it to be equal to $\cos\frac{\pi}{4}=\frac{1}{\sqrt{2}}$.

We have a lower bound! By choosing the closest to the circle between the yellow and green curves, we get an octagon that is very close to the circle. We can define the upper bound of the circle with a function $f$ such as:

$$f(x,y) = \max\left(\max(x,y), \frac{1}{\sqrt{2}}(|x|+|y|)\right)$$

Note that this is different from Paul's article. You do need to take the maximum value of the two norms to select the points that are closest to the center. Generally speaking, for two norms, if one's value is higher than the other, then the former will be drawn closer to the origin when plotting the $\text{norm}(x,y)=1$ curve.

To trace this function, note that Manhattan and infinity norms isolines cross when $|y|=1$ and $|x| = \sqrt{2}-1$ or $|x|=1$ and $|y| = \sqrt{2}-1$.

Finding an upper bound to the norm.

The first idea you can get from the lower bound we found is to scale it up so that the octagon corners touch the circle.

To do so, we need to find the 2-norm of one of the corners and divide $f$ by it.

Let's take the one at $x=1$, $y=\sqrt{2}-1$. We have:

$$\begin{align} \sqrt{x^2+y^2} &=& \sqrt{1 + \left(\sqrt{2}-1\right)^2}\\ &=& \sqrt{1 + 2 - 2\sqrt{2} + 1}\\ &=& \sqrt{4 - 2\sqrt{2}} \end{align}$$

Thus, the upper-bound for the 2-norm with the octagon method is $\sqrt{4 - 2\sqrt{2}}f(x,y)$:

$$f(x,y) \leq \sqrt{x^2+y^2} \leq \sqrt{4 - 2\sqrt{2}}f(x,y)$$

Choosing the best approximation for the norm.

Now, we could stick to Paul Hsieh's choice of taking the middle between the lower and the upper bounds, and it will probably be fine. But come on, let's see if it is the best choice. 😉

Formally, the problem is to find a number $a\in[0,1]$ such as $g$ defined as follows is the closest possible to the norm-2.

$$\begin{align} g(x,y,a) &=& (1-a)f(x,y)+\frac{a}{\sqrt{4 - 2\sqrt{2}}}f(x,y)\\ &=& \left((1-a) + a\sqrt{4 - 2\sqrt{2}}\right)f(x,y) \end{align}$$

Let's plot this function for various values of $a$. To make things easier, I will "unroll" the circle, and plot the norms against $\theta$, the angle between our vector and the $x$ axis.

As expected, we can continuously vary our approximation between the upper and lower bounds. Notice that these functions are periodic and even. We can thus focus on the first half period to minimize the error. The first half period is when the vector is at the first octagon vertices, starting from the $x$ axis and circling anti-clockwise.

To minimize the error with our approximation, we want to minimize the square error. That is:

$$\begin{align} e(a) &=& \int_0^{\arctan\left(\sqrt{2}-1\right)}(g(x,y,a)-1)^2\text{d}\theta \end{align}$$

Thankfully, the expression of $f(x,y)$ and thus of $g(x,y,a)$ should simplify a lot on the given interval. You can see on schematic above that on this interval we have, $f(x,y)=max(|x|,|y|)=|x|=x=\cos\theta$. We can thus rewrite $e(a)$ as follows.

$$\begin{align} e(a) &=& \int_0^{\arctan\left(\sqrt{2}-1\right)}(g(x,y,a)-1)^2\text{d}\theta\\ &=& \int_0^{\arctan\left(\sqrt{2}-1\right)}\left(\left(1-a + a\sqrt{4-2\sqrt{2}}\right)\cos\theta-1\right)^2\text{d}\theta\\ &=& \int_0^{\arctan\left(\sqrt{2}-1\right)}\left(h(a)\cos\theta-1\right)^2\text{d}\theta \end{align}$$

Where $h(a)=\left(1-a + a\sqrt{4-2\sqrt{2}}\right)$ and $\arctan\left(\sqrt{2}-1\right)=\frac{\pi}{8}$.

As we can see from these plots, there is a minimal error, and though 0.5 is a reasonable choice for $a$, we can do slightly better around 0.3.

We can explicitly calculate $e(a)$. Let $h(a)=(1+a(A-1))$. We have

$$\begin{align} e(a) &=& \int_0^{\pi/8}(h(a)\cos\theta-1)^2\text{d}\theta\\ &=&h^2(a)\int_0^{\pi/8}\cos^2\theta\text{d}\theta-2h(a)\int_0^{\pi/8}\cos\theta\text{d}\theta + \frac{\pi}{8}\\ &=& h^2(a)B-2h(a)\sin\frac{\pi}{8} + \frac{\pi}{8} \end{align}$$

Where $B=\frac{\pi}{16}+\frac{1}{4\sqrt2}$. Thus, we look for the position of the minimum, that is where $e'(a)=0$.

$$\begin{align} 0 &=& 2B(A-1)(1+a(A-1))-\sin\frac{\pi}{8}\\ 0 &=& 2B(A-1)(1+a(A-1)) - \frac{A}{2\sqrt2}\\ a &=& \left(\frac{A}{2B\sqrt2}-1\right)\times\frac{1}{A-1}\\ a &\approx& 0.311 \end{align}$$

Not that far from 0.3!

The maximum deviation from the result is then $\max_\theta{|h(a)\cos\theta-1|}$. Looking for that maximum is like looking for the maximum of $\left(h(a)\cos\theta-1\right)^2$. Long story short, the maxima can only occur on the boundaries of the allowed domain for $\theta$, that is $\theta=0$ or $\theta=\pi/8$, meaning

$$\max_\theta{|h(a)\cos\theta-1|} = \max\left(h(a)-1, \left|h(a)\frac{\sqrt{2-\sqrt{2}}}{2}-1\right|\right)$$

With our choice for $a$, we get $h(a)\approx 1.026$, so the maximum deviation is 0.052. That is, we have at most a 5.3% deviation from the norm-2!

Conclusion

That was a fun Sunday project! Originally this was intended to be included in a longer blog-post that is yet to be finished, but I figured it was interesting enough to have its own post. The take-home message being, you can approximate the Euclidean norm of a vector with:

$$\begin{align} \text{norm}(x,y) &=& \frac{\sqrt{2-\sqrt{2}}}{\frac{\pi}{8}+\frac{1}{2\sqrt{2}}}\max\left(\max(|x|,|y|), \frac{1}{\sqrt{2}}(|x|+|y|)\right)\\ &\approx& 1.026\max\left(\max(|x|,|y|), \frac{1}{\sqrt{2}}(|x|+|y|)\right) \end{align}$$

You'll get at most a 5.3% error. This is a bit different from what's proposed on Paul Hsieh's blog-post. Unless I made a mistake, there might be a typo on his blog!

If you are interested in playing with the code used to generate the figures in this article, have a look at the companion notebook!

As always, if you have any question, or want to add something to this post, you can leave me comment or ping me on Twitter or Mastodon.