Floating point Euclidean norm optimization

Question

Past few days I've been struggling with floating point number exercises, trying to learn how they work and what are their limitations. The exercise I've encountered today is as follows:

The Euclidean norm of an n-dimensional vector $x$ with components $x_{1}...x_{n}$ is defined by: $$||x||_{2}=\left(\sum_{i=1}^{n} x_{i}^{2}\right)^{1/2}$$ How would you avoid overflow and harmful underflow in this computation?

I've come up with this solution, but I'm not sure if it really makes sense:

There are some problems with underflow and overflow that can occur when computing Euclidean norm.Firstly, the component $x_{i}$ of the vector can be in the machine range, however its square root might be not. It is possible, that the square of $x_{i}$ can be either too small (underflow - $x_{i}$ considered to be 0) or too large (overflow) or the sum of the squared components can be too large (overflow again).The solution could be implementing an algorithm, that firstly splits the components of the vector 'x' into 3 parts:

1. small magnitude numbers, for which underflow might occur: $x_{i}<y$
2. medium magnitude numbers, where $y\leq x_{i}<Y$
3. large magnitude numbers, for which $x_{i}>Y$

When $y$ and $Y$ are thresholds set by machine precision limits. Then, small magnitude numbers will be scaled by a scalar $S$ and then their sum of squares 'unscaled' by the same $S$. Similarly, large magnitude numbers will be scaled by an appropriate $s$, and then descaled again.

I've taken the algorithm from this book, but I'm not sure if I understand it correctly and if I explained the solution properly. It would be also helpful to reproduce this algorithm in R language, however, I don't understand the explanation to the point of replicating it.

Answer 1

Roughly right.

You can first calculate the sum of squares. If it overflows, then you find a large power of two M, divide all numbers by M, calculate the norm, and multiply by M. If overflow happens at $2^{1023}$, then you have numbers greater than $2^{500}$. Dividing by M=$2^{700}$ means you get neither overflow nor underflow.

You may have underflow in some of the squares. Check how much error an underflow can produce. Then you check how large a sum wouldn’t be influenced by these errors. If the sum of squares is larger than say $2^{-900}$ then there may be products with underflow but they don’t affect the sum. In that case you are fine.

Otherwise, multiply all numbers by M, calculate the norm, and divide by M.