Skip to main content

Subsection 6.3.3 Dot product: error results

It is of interest to accumulate the roundoff error encountered during computation as a perturbation of input and/or output parameters:

  • \(\check \kappa = ( x + \deltax )^T y \text{;}\)

    (\(\check \kappa\) is the exact output for a slightly perturbed \(x \))

  • \(\check \kappa = x^T ( y + \deltay ) \text{;}\)

    (\(\check \kappa\) is the exact output for a slightly perturbed \(y \))

  • \(\check \kappa = x^T y + \delta\!\kappa \text{.}\)

    (\(\check \kappa\) equals the exact result plus an error)

The first two are backward error results (error is accumulated onto input parameters, showing that the algorithm is numerically stable since it yields the exact output for a slightly perturbed input) while the last one is a forward error result (error is accumulated onto the answer). We will see that in different situations, a different error result may be needed by analyses of operations that require a dot product.

Let us focus on the second result. Ideally one would show that each of the entries of \(y \) is slightly perturbed relative to that entry:

\begin{equation*} \deltay = \left( \begin{array}{c} \sigma_0 \psi_0 \\ \vdots \\ \sigma_{n-1} \psi_{n-1} \end{array} \right) = \left( \begin{array}{ c c c } \sigma_0 \amp \cdots \amp 0 \\ \vdots \amp \ddots \amp \vdots \\ 0 \amp \cdots \amp \sigma_{n-1} \end{array} \right) \left( \begin{array}{c} \psi_0 \\ \vdots \\ \psi_{n-1} \end{array} \right) = \Sigma y, \end{equation*}

where each \(\sigma_i \) is "small" and \(\Sigma = \diag{ \sigma_0, \ldots, \sigma_{n-1} }\text{.}\) The following special structure of \(\Sigma \text{,}\) inspired by (6.3.3) will be used in the remainder of this note:

\begin{equation*} \Sigma^{(n)} = \left\{ \begin{array}{l l} 0 \times 0 \mbox{ matrix } \amp\mbox{ if } n = 0 \\ \theta_1 \amp\mbox{ if } n = 1 \\ \diag{\theta_{n}, \theta_{n}, \theta_{n-1}, \ldots , \theta_2 } \amp \mbox{ otherwise }. \end{array} \right. \end{equation*}

Recall that \(\theta_j \) is an order of magnitude variable with \(\vert \theta_j \vert \leq \gamma_j \text{.}\)

Homework 6.3.3.1.

Let \(k \geq 0 \) and assume that \(\vert \epsilon_1 \vert, \vert \epsilon_2 \vert \leq \meps \text{,}\) with \(\epsilon_1 = 0 \) if \(k = 0 \text{.}\) Show that

\begin{equation*} \FlaTwoByTwoSingleLine { I + \Sigma^{(k)} }{ 0 } { 0 }{ ( 1 + \epsilon_1 ) }( 1 + \epsilon_2 ) = ( I + \Sigma^{(k+1)} ) . \end{equation*}

Hint: Reason the cases where \(k = 0 \) and \(k = 1\) separately from the case where \(k > 1 \text{.}\)

Solution

Case: \(k = 0 \text{.}\)

Then

\begin{equation*} \begin{array}{l} \FlaTwoByTwoSingleLine { I + \Sigma^{(k)} }{ 0 } { 0 }{ ( 1 + \epsilon_1 ) }( 1 + \epsilon_2 ) \\ ~~~=~~~~ \lt k = 0 \mbox{ means } (I + \Sigma^{(k)}) \mbox{ is } 0 \times 0 \mbox{ and } (1+\epsilon_1) = ( 1 + 0 ) \gt \\ { ( 1 + 0 ) }( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ ( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ ( 1 + \theta_1 ) \\ ~~~=~~~~\\ ( I + \Sigma^{(1)} ) . \end{array} \end{equation*}

Case: \(k = 1 \text{.}\)

Then

\begin{equation*} \begin{array}{l} \FlaTwoByTwoSingleLine { I + \Sigma^{(k)} }{ 0 } { 0 }{ ( 1 + \epsilon_1 ) }( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ \FlaTwoByTwoSingleLine { 1 + \theta_1 }{ 0 } { 0 }{ ( 1 + \epsilon_1 ) }( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ \FlaTwoByTwoSingleLine { (1 + \theta_1)( 1 + \epsilon_2 ) }{ 0 } { 0 }{ ( 1 + \epsilon_1 )( 1 + \epsilon_2 ) } \\ ~~~=~~~~\\ \FlaTwoByTwoSingleLine { (1 + \theta_2) }{ 0 } { 0 }{ ( 1 + \theta_2 )} \\ ~~~=~~~~\\ ( I + \Sigma^{(2)} ) . \end{array} \end{equation*}

Case: \(k > 1 \text{.}\)

Notice that

\begin{equation*} \begin{array}{l} ( I + \Sigma^{(k)} )( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ \left( \begin{array}{c c c c c} 1+\theta_k \amp 0 \amp 0 \amp \cdots \amp 0 \\ 0 \amp 1+\theta_k \amp 0 \amp \cdots \amp 0 \\ 0 \amp 0 \amp 1+\theta_{k-1} \amp \cdots \amp 0 \\ \vdots \amp \vdots \amp \vdots \amp \ddots \amp \vdots \\ 0 \amp 0 \amp 0 \amp \cdots \amp 1+\theta_2 \end{array} \right) ( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ \left( \begin{array}{c c c c c} 1+\theta_{k+1} \amp 0 \amp 0 \amp \cdots \amp 0 \\ 0 \amp 1+\theta_{k+1} \amp 0 \amp \cdots \amp 0 \\ 0 \amp 0 \amp 1+\theta_{k} \amp \cdots \amp 0 \\ \vdots \amp \vdots \amp \vdots \amp \ddots \amp \vdots \\ 0 \amp 0 \amp 0 \amp \cdots \amp 1+\theta_3 \end{array} \right) \end{array} \end{equation*}

Then

\begin{equation*} \begin{array}{l} \FlaTwoByTwoSingleLine { I + \Sigma^{(k)} }{ 0 } { 0 }{ ( 1 + \epsilon_1 ) }( 1 + \epsilon_2 ) \\ ~~~=~~~~\\ \FlaTwoByTwoSingleLine { ( I + \Sigma^{(k)})( 1 + \epsilon_2 ) }{ 0 } { 0 }{ ( 1 + \epsilon_1 )( 1 + \epsilon_2 ) } \\ ~~~=~~~~\\ \FlaTwoByTwoSingleLine { \left( \begin{array}{c c c c c} 1+\theta_{k+1} \amp 0 \amp 0 \amp \cdots \amp 0 \\ 0 \amp 1+\theta_{k+1} \amp 0 \amp \cdots \amp 0 \\ 0 \amp 0 \amp 1+\theta_{k} \amp \cdots \amp 0 \\ \vdots \amp \vdots \amp \vdots \amp \ddots \amp \vdots \\ 0 \amp 0 \amp 0 \amp \cdots \amp 1+\theta_3 \end{array} \right) }{ 0 } { 0 }{ ( 1 + \theta_2 ) } \\ ~~~=~~~~\\ ( I + \Sigma^{(k+1)} ) . \end{array} \end{equation*}

We state a theorem that captures how error is accumulated by the algorithm.

Proof by Mathematical Induction on \(n\text{,}\) the size of vectors \(x\) and \(y\text{.}\)

When going through the below proof, it will help to look at Section 3 in

  • [7] Paolo Bientinesi, Robert A. van de Geijn, The Science of Deriving Stability Analyses, FLAME Working Note #33. Aachen Institute for Computational Engineering Sciences, RWTH Aachen. TR AICES-2008-2. November 2008.

which relates the algorithm for computing a dot product using our "partition-repartition" FLAME notation to the choices \(x_T \text{,}\) \(y_T \text{,}\) etc. It also links the inductive proof to that algorithm.

  • Base case.

    \(m( x ) = m( y ) = 0 \text{.}\) Trivial!

  • Inductive Step.

    Inductive Hypothesis (I.H.): Assume that if \(x_T, y_T \in \R^k \text{,}\) \(k > 0\text{,}\) then

    \begin{equation*} \fl{ x_T^T y_T } = x_T^T (I + \Sigma_T ) y_T , \mbox{ where } \Sigma_T = \Sigma^{(k)}. \end{equation*}

    We will show that when \(x_T, y_T \in \R^{k+1}\text{,}\) the equality \(\fl{ x_T^T y_T } = x_T^T (I + \Sigma_T ) y_T\) holds true again. Assume that \(x_T, y_T \in \R^{k+1}\text{,}\) and partition \(x_T \rightarrow \FlaTwoByOneSingleLine{ x_0 }{ \chi_1 }\) and \(y_T \rightarrow \FlaTwoByOneSingleLine{ y_0 }{ \psi_1 }\text{.}\) Then

    \begin{equation*} \begin{array}{l} \fl{ \FlaTwoByOneSingleLine{ x_0 }{ \chi_1 }^T \FlaTwoByOneSingleLine{ y_0 }{ \psi_1 } } \\ ~~~=~~~~ \lt {\rm definition} \gt \\ \fl{( \fl{ x_0^T y_0 } + \fl{\chi_1 \psi_1} }\\ ~~~=~~~~ \lt {\rm I.H.~with~} x_T = x_0, y_T = y_0, {\rm ~and~} \Sigma_0 = \Sigma^{(k)} \gt \\ \fl{ x_0^T ( I + \Sigma_0) y_0 + \fl{ \chi_1 \psi_1 }} \\ ~~~=~~~~ \lt {\rm SCM,~twice} \gt \\ \left( x_0^T ( I + \Sigma_0 ) y_0 + \chi_1 \psi_1 (1+\epsilon_* )\right) ( 1 + \epsilon_+ ) \\ ~~~=~~~~ \lt {\rm rearrangement} \gt \\ \FlaTwoByOneSingleLine{ x_0 }{ \chi_1 }^T \FlaTwoByTwoSingleLine { ( I +\Sigma_0 ) }{ 0 } { 0 }{ ( 1+ \epsilon_* ) }( 1+ \epsilon_+ ) \FlaTwoByOneSingleLine{ y_0 }{ \psi_1 } \\ ~~~=~~~~ \lt {\rm renaming} \gt \\ x_T^T ( I + \Sigma_T ) y_T \\ \end{array} \end{equation*}

    where \(\vert \epsilon_* \vert, \vert \epsilon_+ \vert \leq \meps \text{,}\) \(\epsilon_+ = 0 \) if \(k = 0 \text{,}\) and

    \begin{equation*} ( I + \Sigma_T ) = \FlaTwoByTwoSingleLine { ( I +\Sigma_0 ) }{ 0 } { 0 }{ ( 1+ \epsilon_*) }( 1+ \epsilon_+) \end{equation*}

    so that \(\Sigma_T = \Sigma^{(k+1)} \text{.}\)

  • By the Principle of Mathematical Induction, the result holds.

A number of useful consequences of Theorem 6.3.3.1 follow. These will be used later as an inventory (library) of error results from which to draw when analyzing operations and algorithms that utilize a dot product.

R-1B

We leave the proof of Corollary 6.3.3.2 R-1B as an exercise.

R-2B

The proof of Corollary 6.3.3.2 R-2B is, of course, just a minor modification of the proof of Corollary 6.3.3.2 R-1B.

R-1F

For Corollary 6.3.3.2 R-1F, let \(\delta\!\kappa = x^T \Sigma^{(n)} y \text{,}\) where \(\Sigma^{(n)} \) is as in Theorem 6.3.3.1. Then

\begin{equation*} \begin{array}{rcl} \vert \delta\!\kappa \vert \amp=\amp \vert x^T \Sigma^{(n)} y \vert \\ \amp\leq \amp \vert \chi_0 \vert \vert \theta_{n} \vert \vert \psi_0 \vert + \vert \chi_1 \vert \vert \theta_{n} \vert \vert \psi_1 \vert + \cdots + \vert \chi_{n-1} \vert \vert \theta_{2} \vert \vert \psi_{n-1} \vert \\ \amp\leq \amp \gamma_n \vert \chi_0 \vert \vert \psi_0 \vert + \gamma_n \vert \chi_1 \vert \vert \psi_1 \vert + \cdots + \gamma_2 \vert \chi_{n-1} \vert \vert \psi_{n-1} \vert \\ \amp\leq \amp \gamma_n \vert x \vert^T \vert y \vert. \end{array} \end{equation*}
Homework 6.3.3.2.

Prove Corollary 6.3.3.2 R1-B.

Solution

From Theorem 6.3.3.1 we know that

\begin{equation*} \check \kappa = x^T ( I + \Sigma^{(n)} ) y = ( x + \begin{array}[t]{c} \underbrace{\Sigma^{(n)} x}\\ \deltax \end{array} )^T y . \end{equation*}

Then

\begin{equation*} \begin{array}{rcl} \vert \deltax \vert \amp=\amp \vert \Sigma^{(n)} x \vert = \left\vert \left( \begin{array}{c} \theta_n \chi_0 \\ \theta_n \chi_1 \\ \theta_{n-1} \chi_2 \\ \vdots \\ \theta_2 \chi_{n-1} \end{array} \right) \right\vert = \left( \begin{array}{c} \vert \theta_n \chi_0 \vert \\ \vert \theta_n \chi_1 \vert \\ \vert \theta_{n-1} \chi_2 \vert \\ \vdots \\ \vert \theta_2 \chi_{n-1} \vert \end{array} \right) = \left( \begin{array}{c} \vert \theta_n \vert \vert \chi_0 \vert \\ \vert \theta_n \vert \vert \chi_1 \vert \\ \vert \theta_{n-1} \vert \vert \chi_2 \vert \\ \vdots \\ \vert \theta_2 \vert \vert \chi_{n-1} \vert \end{array} \right) \\ \amp\leq\amp \left( \begin{array}{c} \vert \theta_n \vert \vert \chi_0 \vert \\ \vert \theta_n \vert \vert \chi_1 \vert \\ \vert \theta_n \vert \vert \chi_2 \vert \\ \vdots \\ \vert \theta_n \vert \vert \chi_{n-1} \vert \end{array} \right) \leq \left( \begin{array}{c} \gamma_n \vert \chi_0 \vert \\ \gamma_n \vert \chi_1 \vert \\ \gamma_n \vert \chi_2 \vert \\ \vdots \\ \gamma_n \vert \chi_{n-1} \vert \end{array} \right) = \gamma_n \vert x \vert. \end{array} \end{equation*}

(Note: strictly speaking, one should probably treat the case \(n = 1 \) separately.)

login