Skip to main content

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.)

Consider 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*}
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.

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.

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).

login