Skip to main content

Subsection 1.5.2 Practical computation of the vector 2-norm

Consider the computation \(\gamma = \sqrt{ \alpha^2 + \beta^2} \) where \(\alpha, \beta \in \R\text{.}\) When computing this with floating point numbers, a fundamental problem is that the intermediate values \(\alpha^2 \) and \(\beta^2 \) may overflow (become larger than the largest number that can be stored) or underflow (become smaller than the smallest positive number that can be stored), even if the resulting \(\gamma \) itself does not overflow or underflow.

The solution is to first determine the largest value, \(\mu = \max( \vert \alpha \vert, \vert \beta \vert ) \) and then compute

\begin{equation*} \gamma = \mu \sqrt{ \left( \frac{\alpha}{\mu} \right)^2 + \left( \frac{\beta}{\mu} \right)^2} \end{equation*}

instead. A careful analysis shows that if \(\gamma \) does not overflow, neither do any of the intermediate values encountered during it computation. While one of the terms \(\left( \frac{\alpha}{\mu} \right)^2 \) or \(\left( \frac{\beta}{\mu} \right)^2 \) may underflow, the other one equals one and hence the overall result does not underflow. A complete discussion of all the intracacies go beyond this note.

This insight generalizes to the computation of \(\| x \|_2 \) where \(x \in \Cn \text{.}\) Rather than computing it as

\begin{equation*} \| x \|_2 = \sqrt{ \vert \chi_0 \vert^2 + \vert \chi_1 \vert^2 + \cdots + \vert \chi_{n-1} \vert^2 } \end{equation*}

and risk overflow or underflow, instead the following computation is used:

\begin{equation*} \begin{array}{lcl} \mu \amp =\amp \| x \|_\infty \\ \| x \|_2 \amp =\amp \mu \sqrt{ \left( \frac{\vert \chi_0 \vert}{\mu} \right)^2 + \left( \frac{\vert \chi_1 \vert}{\mu} \right)^2 + \cdots + \left( \frac{\vert \chi_{n-1} \vert}{\mu} \right)^2 }. \end{array} \end{equation*}

What we describe avoids overflow and underflow when not required. The problem is that is not efficient because most of the cost in performing such vector operations is in bringing data from main memory into registers so that computation can occur.

The described method requires the vector to be read from mainn memory three times:

  • Once to determine \(\mu\text{.}\)

  • Once to compute the sum of the squares of the elements of \(x / \mu \text{.}\)

In a recent talk at an annual workshop we organize, the BLIS Retreat 2023, how to implement the computation of the norm in a way that requires the data to only be read once was discussed by Eleni Vlachopoulou from AMD. We believe you may enjoy her talk!

login