Background

In this note we will be working with matrices and vectors. Simply put, matrices are two dimensional arrays and vectors are one dimensional arrays (or the "usual" notion of arrays). We will be using notation that is consistent with array notation. So e.g. a matrix $\mathbf{A}$ with $m$ rows and $n$ columns (also denoted as an $m\times n$ matrix) will in code be defined as int [][] A = new int[m][n] (assuming the matrix stores integers). Also a vector $\mathbf{x}$ of size $n$ in code will be declared as int [] x = new int[n] (again assuming the vector contains integers). To be consistent with the array notations, we will denote the entry in $\mathbf{A}$ corresponding to the $i$th row and $j$th column as $A[i][j]$ (or A[i][j]). Similarly, the $i$th entry in the vector $\mathbf{x}$ will be denoted as $x[i]$ (or x[i]). We will follow the array convention assume that the indices $i$ and $j$ start at $0$.

If you want a refresher on matrices, you might want to start with this Khan academy video (though if you are comfortable with the array analogy above you should not really need much more for this note):

The Problem

If you have not seen this operation before, this Khan academy video has some examples (or read the next section):

The Inner Product

It turns out that to visualize the matrix-vector multiplication, it is helpful to define a simpler operation between two vectors:

Below is an animation that shows the inner product of two random vectors of size $n=6$ (click on "Next" to continue the computation of the inner product: the animation starts with two new random vectors if "Next" is pressed once the inner product has already been computed):

Armed with the definition of the inner product, here is an equivalent definition of the matrix-vector multiplication problem:

Below is one example of a matrix vector multiplication: \[ \begin{pmatrix} 1 & 2 & -3\\ 2 & 9 &0\\ 6 & -1 & -2 \end{pmatrix} \times \begin{pmatrix} 2\\ 3\\ -1 \end{pmatrix} = \begin{pmatrix} 1\times 2 + 2\times 3 + (-3)\times(-1)\\ 2\times 2+ 9\times 3+0\times (-1)\\ 6\times 2+(-1)\times 3+ (-2)\times(-1) \end{pmatrix} = \begin{pmatrix} 11\\ 31\\ 11 \end{pmatrix}. \]

From the above the following should be easy to show:

In general the above cannot be improved:

Structured Matrices

It turns out that in many practical applications of matrix-vector multiplication, the matrix $\mathbf{A}$ has some structure that potentially allows one to perform matrix vector multiplication (asymptotically) faster than $\Omega(n^2)$ time. We illustrate this point with a specific family of structured matrices:

Next, we present an algorithm to multiply the outer product of two arbitrary vectors $\mathbf{s}$ and $\mathbf{t}$ with a vector $\mathbf{x}$ in time $O(n)$. (Note that this is a significant improvement over the runtime from the trivial algorithm in Exercise 1.) In other words, we want to solve the following problem

                                	//Input: s[i], t[i], x[i] for 0<= i<n

					//Compute <s,t> first
					ip=0;
					for(j=0;j < n; j++)
						ip+=t[j]*x[j];

					//Now compute the output vector
					y = new int[n];
					for(i=0; i< n;i++)
						y[i] = s[i]*ip;

					return y;
				

The runtime analysis is left as an exercise:

Discrete Fourier Transform

In this section, we consider perhaps the world's most widely used matrix-vector multiplication: the Discrete Fourier Transform (or DFT). For applications and more math (than below), look at its Wikipedia page . Here is the definition of the Discrete Fourier transform:

Before going a bit more into the math of DFT, here is a video showing the DFT³ of an image of a table cloth:

And playing with DFTs can be addictive:

The complex number thingy

As you must have noticed the DFT uses complex numbers. Here's a video on why complex numbers are really useful (and why "imaginary number" is one of Math's PR mistakes):

And, sorry could not resist this video on "imaginary Erdos number" (skip this if you do not want to get side-tracked):

So it turns out that the complex numbers used in the discrete Fourier matrix are roots of unity . However, for whatever is coming up next, you only need to know the following two simple facts:

Some Comments on Structure of the DFT

We now present some observation on the discrete Fourier matrix that can be used to speed up the computation of the DFT from the trivial $O(n^2)$ bound. For simplicity we will make the following assumption.

Let $0\le i,j\lt n$ be row and column indices and consider the entry $F[i][j]= e^{-2\pi\iota\cdot i\cdot j/n}$. We will express this entry as a product of new numbers. In particular, let $(i_0,\dots, i_{m-1})$ and $(j_0,\dots,j_{m-1})$ be the binary representation of $i$ and $j$ (note that $m$ are sufficient). In other words, \[i= i_0 + i_1\cdot 2+\cdots +i_{m-1}\cdot 2^{m-1},\] and \[j= j_0 + j_1\cdot 2+\cdots +j_{m-1}\cdot 2^{m-1}.\]

Using the above we get that for $0\le i,j\lt n$ \[\frac{i\cdot j}{n} =2^{-m}\cdot \left(\sum_{k=0}^{m-1} i_k\cdot 2^k\right)\cdot \left(\sum_{k=0}^{m-1} j_k\cdot 2^k\right) = \sum_{k,\ell: 0\le k+\ell\le 2m-2} i_k\cdot j_{\ell}\cdot 2^{k+\ell-m}.\] Using the above identity along with Fact 2, we get that \[ F[i][j] = e^{-2\pi\iota\cdot i\cdot j/n} =\prod_{k,\ell: 0\le k+\ell\le 2m-2} e^{-2\pi\iota\cdot i_k\cdot j_{\ell}\cdot 2^{k+\ell-m}}.\] Using Fact 1, we note that for every $k+\ell\ge m$, we have that $e^{-2\pi\iota\cdot i_k\cdot j_{\ell}\cdot 2^{k+\ell-m}}=1$. This implies that \[ F[i][j] = \prod_{k,\ell: 0\le k+\ell\le m-1} e^{-2\pi\iota\cdot i_k\cdot j_{\ell}\cdot 2^{k+\ell-m}}.\] In other words, we have reduced each entry from being a product of $2m-1$ numbers to being a product of only $m$ numbers. We record this fact

It turns out that the above savings of "only" $m-1$ products (i.e. we went down from a product of $2m-1$ terms to a product of $m$ terms) is enough to drive down the time to compute the DFT from $O(n^2)$ to $O(n\log{n})$:

The Distributive Law

In this section, we will present a simple observation that turns out to be very useful in designing efficient algorithms (especially those involving addition and multiplication of numbers). In particular, this idea (suitably implemented) can also be used to solve Exercise 4 above.

We begin with the statement of the distributive law:

Note that the above implies that \[\sum_{i=0}^{n-1} a\cdot b_i = a\cdot\left(\sum_{i=0}^{n-1} b_i\right).\] The above implies that instead of the $n$ multiplications on the LHS, we can get away with just one multiplication on the RHS. Of course in the case above we still need $\Theta(n)$ arithmetic operation in both cases so we do not get any asymptotic improvement.

However, if we combine the distributive law with some judicious saving of intermediate results then we can result in asymptotic improvement in runtimes. We actually have already seen such an example. Recall the example of multiplying the outer product of two vectors $\mathbf{s}$ and $\mathbf{t}$ with an arbitrary vector $\mathbf{x}$, we saw that $y[i]=s[i]\cdot \left\langle \mathbf{t},\mathbf{x}\right\rangle$. The expression was achieved from the original expression of $\sum_{j=0}^{n-1} s[i]\cdot t[j]\cdot x[j]$ by using the distributive law. Then we noticed that if we stored the "common" term $\left\langle \mathbf{t},\mathbf{x}\right\rangle$, then we could avoid recomputing this term for different $i$. We summarize this strategy as follows:

Computing Sums of Products

It turns out that the distributive law strategy is especially good at efficiently computing sums of products. Instead of defining more formally what I mean by sums of products, let me illustrate my point with a specific problem:⁵

Consider the following trivial algorithm to solve this problem. Compute a $4$-D array defined as (where $0\le i,j,k,\ell\lt n$): \[T[i][j][k][\ell] = A[i][j]\cdot B[j][k]\cdot C[k][\ell]\cdot x[\ell].\] Note that $T$ can be computed in $O(n^4)$ time⁶. Then in further $O(n^3)$ time one can compute $y[i]$ for each $0\le i\lt n$ using the definition of $\mathbf{y}$ above. Thus, we can compute $\mathbf{y}$ in $O(n^4)$ overall. Next, we show that one can compute $\mathbf{y}$ must faster in $O(n^2)$ time using the distributive law strategy.

                                //Input: A[i][j], B[j][k], C[k][l], x[l] for  for 0<= i,j,k,l<n

					z = new int[n];
					u = new int[n];
					y = new int[n];

					//Compute the vector z first
					for( k=0; k< n; k++)
					{
						z[k]=0;
						for(l=0; l< n; l++)
						z[k]+= C[k][l]*x[l];
					}

					//Compute the u vector next
					for( j=0; j< n; j++)
					{
						u[j]=0;
						for(k=0; k< n; k++)
						u[j]+= B[j][k]*z[k];
					}

					//Finally compute the output y
					for(i=0; i< n;i++)
					{
						y[i]=0;
						for(j=0; j< n;j++)
							y[i]+= A[i][j]*u[j];
					}

					return y;
				

We leave the runtime analysis of the algorithm above as an exercise:

What just happened here?

We now return to the unresolved question of why the sums of product considered above is a natural question in itself. It is not too hard to see that the sums of products can also be stated as \[\mathbf{y}= \mathbf{A}\cdot \mathbf{B}\cdot \mathbf{C}\cdot \mathbf{x}.\] Further, the algorithm discussed above basically follows the obvious way to compute $\mathbf{y}$. That is compute the RHS as \[\mathbf{A}\cdot \underbrace{\left(\mathbf{B}\cdot \underbrace{\left(\mathbf{C}\cdot \mathbf{x}\right)}_{\mathbf{z}}\right)}_{\mathbf{u}}.\]