Euclidean Norm of a Vector: Best Routine?

[DavidB]

I am looking for a robust routine for calculating the Euclidean norm of a real vector. I came across SNRM2 and DNRM2, but they seem to be kind of old.

Has anything newer and better been developed since these routines were written? I am wondering, is there a "best" routine out there for calculating the Euclidean norm of a vector?

[Julien Langou]

The four xNRM2 routines that you can find in the reference BLAS are
reference implementations. Therefore they are not meant to be
state-of-the-art implementations but they are meant to be robust and
simple (understandable, readable) pieces of code that execute the
functionnality. Maybe you will find better implementations in some
optimized BLAS libraries (ATLAS, Goto BLAS for the open source ones).

That said, the current implementation of xNRM2 in the reference BLAS
(one of the many contributions of Sven Hammarling) is not that simple.

First, note that a code like: sqrt(x^Tx) is not acceptable. Assume that the
largest number, we can store in our computer is 1.8e308. Assume that the
vector X is made of two entries: 3e200 and 4e200. The Euclidean norm of
X is 5e200 which is a perfectly legitimate number on this system (i.e. you
can store it since it is smaller than 1.8e308). However, the computation
of sqrt(x^Tx) leads in first instance to x^TX=25e400 which would
overflow on our computer. Consequently, our result would be +Inf and not
5e200.

Fine, the solution is easy. In case of overflow, we need to scale our data
accordingly. In our case, simply pre-divide X by 1e200, find 5 as the
norm, and post-multiply by 1e200 to get 5e200. That's what the netlib
code is doing.

Since xNRM2 is a Level-1 BLAS subroutine, the efficiency of this routine is
all about memory access, for each memory reference you will perform
about three operations. Not a lot. Therefore the reference BLAS
implementation of xNRM2 uses a trick to perform as few as possible
passes on the original data. In this matter, the proposed algorithm is
optimal: it provides the solution in only one pass on the original data
(obviously the minimum bound). The trick is to compute the scaling factor
on the fly while computing the norm.

The original code and some comments are in the LAPACK routine xLASSQ,
SSQ standing for 'sum of squares'.

Some variants have recently been proposed by Fred Gustavson. The idea
is to use what Fred calls the 'fast path'. The problem of the current xNRM2
(and xLASSQ) routine is that there is a branching (IF statement) inside the
loop (DO statement) and therefore the code is almost impossible to
optimize correctly (by hand or by a compiler). This is the price for
robustness. Since the overflow cases are not supposed to happen all that
much in practice, Fred Gustavson proposed to use by default sqrt(x^Tx)
code, this code does not have any branching, can be easily vectorized,
performs one path on the data, etc., so it is fast. This code however
might overflow so if this code overflows (result is +Inf) then go back to
the current algorithm.

Why not, looks like a good idea.

There are other cases where the conservativeness of LAPACK has a cost
on the performance of the algorithm, while this conversativeness is only
justified in a few particular cases. Some fast path algorithms could used
there too. (Fred Gustavson told me at some points about DLATRS.) Those
algorithms would be faster in most of the cases, and would fail in some
cases. So they might be considered overall better. On going conversation.

[DavidB]

Hi, Julien.

Thank-you for the informative response.

I suppose, given a choice between speed and robustness, I would choose the routine that is robust, the one that yields accurate results for the widest range of inputs.

I notice that there are two versions of xNRM2 routines posted, early ones and others featuring newer code supplied by Sven Hammarling. The routines supplied by Sven Hammarling do not use the machine constants DBL_EPSILON, DBL_MIN, etc. but they are certainly easier to understand than the earlier code. I tried to understand the earlier code but there are too many GOTO statements; I gave up in frustration.

In any case, I will keep using DNRM2, but I will check out some of the other routines you mentioned in your post.

Thanks,


David