PfHP Matrix-vector multiplication via dot products

Unit 1.3.3 Matrix-vector multiplication via dot products

Homework 1.3.3.1.

Compute $\left(\begin{array}{rrrr} 2 \amp 2 \amp -1 \amp 2\\ 2 \amp 1 \amp 0 \amp -2\\ -2 \amp -2 \amp 2 \amp 2 \end{array}\right) \left(\begin{array}{r} 2\\ -1\\ 0\\ -1 \end{array}\right) =$

Solution

\begin{equation*} \begin{array}{rcl} \left(\begin{array}{rrrr} 2 \amp 2 \amp -1 \amp 2\\ \hline 2 \amp 1 \amp 0 \amp -2\\ \hline -2 \amp -2 \amp 2 \amp 2 \end{array}\right) \left(\begin{array}{r} 2\\ -1\\ 0\\ -1 \end{array}\right) \amp=\amp \left( \begin{array}{r} \left(\begin{array}{rrrr} 2 \amp 2 \amp -1 \amp 2 \end{array}\right) \left(\begin{array}{r} 2\\ -1\\ 0\\ -1 \end{array}\right) \\ \hline \left(\begin{array}{rrrr} 2 \amp 1 \amp 0 \amp -2 \end{array}\right) \left(\begin{array}{c} 2\\ -1\\ 0\\ -1 \end{array}\right) \\ \hline \left(\begin{array}{rrrr} -2 \amp -2 \amp 2 \amp 2 \end{array}\right) \left(\begin{array}{c} 2\\ -1\\ 0\\ -1 \end{array}\right) \end{array} \right) \\ \amp=\amp \left( \begin{array}{r} 0 \\ 5 \\ -4 \end{array} \right) \end{array} \end{equation*}

(Here we really care more about exposing the dot products of rows with the vector than the final answer.)

permalink

permalinkConsider the matrix-vector multiplication

$\begin{equation*} y := A x + y . \end{equation*}$

The way one is usually taught to compute this operation is that each element of $y \text{,}$ $\psi_i \text{,}$ is updated with the dot product of the corresponding row of $A \text{,}$ $\widetilde a_i^T \text{,}$ with vector $x \text{.}$ With our notation, we can describe this as

$\begin{equation*} \begin{array}{rcl} \left( \begin{array}{c} \psi_0 \\ \hline \psi_1 \\ \hline \vdots \\ \hline \psi_{m-1} \end{array} \right) \amp:=\amp \left( \begin{array}{c} \widetilde a_0^T \\ \hline \widetilde a_1^T \\ \hline \vdots \\ \hline \widetilde a_{m-1}^T \end{array} \right) x + \left( \begin{array}{c} \psi_0 \\ \hline \psi_1 \\ \hline \vdots \\ \hline \psi_{m-1} \end{array} \right) \\ \amp = \amp \left( \begin{array}{c} \widetilde a_0^T x\\ \hline \widetilde a_1^T x\\ \hline \vdots \\ \hline \widetilde a_{m-1}^T x \end{array} \right) + \left( \begin{array}{c} \psi_0 \\ \hline \psi_1 \\ \hline \vdots \\ \hline \psi_{m-1} \end{array} \right) = \left( \begin{array}{c} \widetilde a_0^T x + \psi_0 \\ \hline \widetilde a_1^T x + \psi_1 \\ \hline \vdots \\ \hline \widetilde a_{m-1}^T x + \psi_{m-1} \end{array} \right). \end{array} \end{equation*}$

permalink

Remark 1.3.3.

The way we remember how to do the partitioned matrix-vector multiplication

$\begin{equation*} \left( \begin{array}{c} \widetilde a_0^T \\ \hline \widetilde a_1^T \\ \hline \vdots \\ \hline \widetilde a_{m-1}^T \end{array} \right) x = \left( \begin{array}{c} \widetilde a_0^T x \\ \hline \widetilde a_1^T x \\ \hline \vdots \\ \hline \widetilde a_{m-1}^T x \end{array} \right) \end{equation*}$

is to replace the $m$ with a specific integer, $3 \text{,}$ and each symbol with a number:

$\begin{equation*} \left( \begin{array}{r} 1 \\ \hline -2 \\ \hline 3 \end{array} \right) (-1) . \end{equation*}$

If you view the first as a $3 \times 1$ matrix and the second as a $1 \times 1$ matrix, the result is

$\begin{equation} \left( \begin{array}{r} 1 \\ \hline -2 \\ \hline 3 \end{array} \right) (-1) = \left( \begin{array}{r} (1) \times (-1)\\ \hline (-2) \times (-1) \\ \hline (3) \times (-1) \end{array} \right) .\label{rowpartition-mv1}\tag{1.3.1} \end{equation}$

You should then notice that multiplication with the partitioned matrix and vector works the same way:

$\begin{equation} \left( \begin{array}{c} \widetilde a_0^T \\ \hline \widetilde a_1^T \\ \hline \vdots \\ \hline \widetilde a_{m-1}^T \end{array} \right) x = \left( \begin{array}{c} (\widetilde a_0^T) \times ( x )\\ \hline (\widetilde a_1^T) \times ( x ) \\ \hline \vdots \\ \hline ( \widetilde a_{m-1}^T ) \times ( x ) \end{array} \right) \label{rowpartition-mv2}\tag{1.3.2} \end{equation}$

with the change that in (1.3.1) the multiplication with numbers commutes while in (1.3.1) the multiplication does not (in general) commute:

$\begin{equation*} (-2) \times (-1) = (-1) \times (-2) \end{equation*}$

but, in general,

$\begin{equation*} \widetilde a_{i}^T x \neq x \widetilde a_i^T. \end{equation*}$

Indeed, $\widetilde a_i^T x$ is a dot product (inner product) while we will later be reminded that $x \widetilde a_i^T$ is an outerproduct of vectors $\widetilde a_i$ and $x \text{.}$

If we then expose the individual elements of $A$ and $y$ we get

$\begin{equation*} \begin{array}{l} \left( \begin{array}{c} \psi_0 \\ \psi_1 \\ \vdots \\ \psi_{m-1} \end{array}\right) \\ ~~~:= \left( \begin{array}{c @{~} c @{~} c @{~} c @{~} c @{~} c @{~} c @{~} c @{~} c @{~} c @{~} c} \alpha_{0,0} \amp \chi_0 \amp + \amp \alpha_{0,1} \amp \chi_1 \amp + \amp \cdots \amp \alpha_{0,n-1} \amp \chi_{n-1} \amp+ \amp \psi_0 \\ \alpha_{1,0} \amp \chi_0 \amp + \amp \alpha_{1,1} \amp \chi_1 \amp + \amp \cdots \amp \alpha_{1,n-1} \amp \chi_{n-1} \amp + \amp \psi_1 \\ \amp\amp\amp\amp\amp\vdots\amp\amp\amp \\ \alpha_{m-1,0} \amp \chi_0 \amp + \amp \alpha_{m-1,1} \amp \chi_1 \amp + \amp \cdots \amp \alpha_{m-1,n-1} \amp \chi_{n-1} \amp + \amp \psi_{m-1} \end{array} \right) \\ \end{array} \end{equation*}$

This discussion explains the IJ loop for computing $y := A x + y \text{:}$

$\begin{equation*} \begin{array}{l} {\bf for~} i := 0, \ldots , m-1 \\[0.05in] ~~~ \left. \begin{array}{l} {\bf for~} j := 0, \ldots , n-1 \\ ~~~ ~~~ \psi_i := \alpha_{i,j} \chi_{j} + \psi_i \\ {\bf end} \end{array} \right\} ~~~ \psi_i := \widetilde a_i^T x + \psi_i \\[0.1in] {\bf end} \end{array} \end{equation*}$

where we notice that again the $i$ index ranges over the rows of the matrix and $j$ index ranges over the columns of the matrix.

permalink

Homework 1.3.3.2.

In directory Assignments/Week1/C complete the implementation of matrix-vector multiplication in terms of dot operations

#define alpha( i,j ) A[ (j)*ldA + i ]   // map alpha( i,j ) to array A 
#define chi( i )  x[ (i)*incx ]            // map chi( i )  to array x
#define psi( i )  y[ (i)*incy ]            // map psi( i )  to array y

void Dots( int, const double *, int, const double *, int, double * );

void MyGemv( int m, int n, double *A, int ldA,
           double *x, int incx, double *y, int incy )
{
  for ( int i=0; i<m; i++ )
    Dots(   ,      ,     ,     ,    ,      );
}

in file Assignments/Week1/C/Gemv_I_Dots.c. You can test it by executing

make I_Dots

Upon completion, this should print

It appears all is well

Admittedly, this is a very incomplete test of the routine. However, it doesn't play much of a role in the course, so we choose to be sloppy.

Solution

Assignments/Week1/Answers/Gemv_I_Dots.c.

permalink

Remark 1.3.4.

The file name Gemv_I_Dots.c can be decoded as follows:

The Gemv stands for "GEneral Matrix-Vector multiplication."
The I indicates a loop indexed by $i$ (a loop over rows of the matrix).
The Dots indicates that the operation performed in the body of the loop is a dot product (hence the Dot) with the result added to a scalar (hence the s).