Lesson 11 - Matrix and Vector Operations

Learning goals. By the end of this lecture, students should be able to:

• interpret vectors as ordered collections of unknowns, data, or coefficients,

• compute and interpret dot products as weighted sums,

• compute matrix-vector products using both row and column viewpoints,

• explain the loop structure and approximate cost of a dense matrix-vector product,

• use matrix-matrix products to represent many related matrix-vector products at once.

Key terms. vector, dot product, matrix-vector product, row view, column view.

Why this chapter follows balance laws

In Lesson 10, we assembled a model in the form

\[\mathbf{A}\mathbf{x}=\mathbf{b}\]

and immediately solved it with a black-box routine. The solution was useful, but the notation hid several operations. The expression \(\mathbf{A}\mathbf{x}\) is not a single multiplication. It is a structured collection of many multiplications and additions.

This chapter opens that box. We will treat vectors and matrices both as mathematical objects and as arrays of numbers stored in a computer.

Vectors as ordered data

A vector \(\mathbf{x}\in\mathbb{R}^n\) is an ordered list of \(n\) real numbers, written as

\[\begin{split}\mathbf{x}=\begin{bmatrix}x_1\\x_2\\ \vdots \\x_n\end{bmatrix}.\end{split}\]

Here \(\mathbb{R}^n\) means the set of all length-\(n\) real-valued vectors. An entry \(x_i\) usually has meaning: it might represent the temperature, neutron flux, or mechanical stress in the \(i\)th region of a reactor model. It could also be the \(i\)th coefficient in a fitted heat-transfer model. The order matters because entry \(x_i\) is tied to the \(i\)th feature, cell, basis function, or equation of the model.

A vector can also be scaled and added:

\[\begin{split}\alpha \mathbf{u}+\beta \mathbf{v} = \begin{bmatrix} \alpha u_1+\beta v_1\\ \alpha u_2+\beta v_2\\ \vdots\\ \alpha u_n+\beta v_n \end{bmatrix}.\end{split}\]

Here \(\alpha,\beta\in\mathbb{R}\) are scalars. This operation appears whenever we combine sources, update iterates, blend states, or form linear combinations of basis vectors. The phrase linear combination means exactly this: multiply vectors by scalars, then add the results.

Worked example: combining source patterns.

Problem. Suppose a four-cell model has basis source patterns in cells 1 and 3:

\[\begin{split}\mathbf{s}_1=\begin{bmatrix}1\\0\\0\\0\end{bmatrix} \qquad \text{and} \qquad \mathbf{s}_3=\begin{bmatrix}0\\0\\1\\0\end{bmatrix}.\end{split}\]

Build the source vector for amplitudes \(2~\mathrm{n/s}\) in cell 1 and \(0.5~\mathrm{n/s}\) in cell 3.

Solution. The patterns are dimensionless. The amplitudes carry the physical source units, so the integrated source vector is the linear combination

\[\begin{split}\mathbf{b}=2\mathbf{s}_1+\frac{1}{2}\mathbf{s}_3 = \begin{bmatrix}2\\0\\1/2\\0\end{bmatrix}.\end{split}\]

The entries of \(\mathbf{b}\) have units \(\mathrm{n/s}\). If the same algebra is used with normalized amplitudes, then the entries are dimensionless fractions instead. The linear-combination rule is the same; the units ride on the amplitudes.

Dot products as weighted sums

The dot product of two vectors \(\mathbf{u},\mathbf{v}\in\mathbb{R}^n\) is

\[\mathbf{u}^{\mathsf{T}}\mathbf{v}=\sum_{j=1}^n u_j v_j.\]

The dot product is a weighted sum: multiply corresponding entries and add the results. If \(\mathbf{u}\) contains weights and \(\mathbf{v}\) contains state values, then \(\mathbf{u}^{\mathsf{T}}\mathbf{v}\) is the weighted total.

Worked example: reactor power as a dot product.

Problem. Suppose a nuclear reactor model has three regions. Let

\[\begin{split}\boldsymbol{\phi}=\begin{bmatrix}1.2\\0.9\\0.4\end{bmatrix} \times 10^{12}\ \mathrm{n/(cm^2\,s)}\end{split}\]

be region-averaged scalar neutron fluxes. Also suppose the region volumes and macroscopic fission cross sections combine to give

\[\begin{split}\mathbf{w}=\begin{bmatrix}0.50\\0.75\\1.00\end{bmatrix} \ \mathrm{cm^2}, \qquad w_i=V_i\Sigma_{f,i}.\end{split}\]

Estimate the total fission rate and corresponding power using \(3.2\times 10^{-11}\ \mathrm{J/fission}\).

Solution. Here \(V_i\) has units \(\mathrm{cm^3}\), \(\Sigma_{f,i}\) has units \(\mathrm{cm^{-1}}\), and therefore \(w_i=V_i\Sigma_{f,i}\) has units \(\mathrm{cm^2}\). This is a reaction-rate model: the macroscopic fission cross section converts scalar-flux track-length rate into fission reaction rate. The total fission rate is

\[\mathbf{w}^{\mathsf{T}}\boldsymbol{\phi} = \left[0.50(1.2)+0.75(0.9)+1.00(0.4)\right]\times 10^{12} = 1.675\times 10^{12}\ \mathrm{fissions/s}.\]

The unit check is part of the calculation:

\[(\mathrm{cm^2})\,\frac{\mathrm{n}}{\mathrm{cm^2\,s}} = \frac{\mathrm{fissions}}{\mathrm{s}}.\]

The entries of \(\mathbf{w}\) are physical weights. The entries of \(\boldsymbol{\phi}\) are state values. The dot product forms the weighted total over the regions.

Multiplying that fission rate by \(3.2\times 10^{-11}\ \mathrm{J/fission}\) gives an estimated power of

\[\left(1.675\times 10^{12}\ \frac{\mathrm{fissions}}{\mathrm{s}}\right) \left(3.2\times 10^{-11}\ \frac{\mathrm{J}}{\mathrm{fission}}\right) \approx 53.6\ \mathrm{W}.\]

For example,

\[\begin{split}\begin{bmatrix}2&-1&0\end{bmatrix} \begin{bmatrix}1\\2\\3\end{bmatrix} =2(1)-1(2)+0(3)=0.\end{split}\]

For abstract vectors, or for vectors whose entries have been nondimensionalized or share compatible units, the same operation is also a geometric measure of alignment. That interpretation becomes central in Lesson 15.

A quick preview of that geometry is useful now. The Euclidean length of a vector is

\[\lVert \mathbf{v} \rVert_2=\left(\sum_{j=1}^n v_j^2\right)^{1/2}.\]

For two nonzero vectors, the dot product determines the angle \(\theta\) between them:

\[\cos\theta=\frac{\mathbf{u}^{\mathsf{T}}\mathbf{v}}{\lVert \mathbf{u} \rVert_2 \lVert \mathbf{v} \rVert_2}.\]

If \(\cos\theta=1\), the vectors point in the same direction. If \(\cos\theta=0\), the vectors are orthogonal. Lesson 15 develops this geometric interpretation more carefully.

import numpy as np

u = np.array([2.0, -1.0, 0.0])
v = np.array([1.0, 2.0, 3.0])

s_loop = 0.0
for j in range(len(v)):
    s_loop += u[j] * v[j]

s_numpy = u @ v
print(s_loop, s_numpy)

0.0 0.0

Matrix-vector products: the row view

Let

\[\begin{split}\mathbf{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & a_{22} & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix}, \qquad \mathbf{x}=\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}.\end{split}\]

We write \(\mathbf{A}_{i,:}\) for row \(i\) of \(\mathbf{A}\), and \(\mathbf{A}_{:,j}\) for column \(j\). In NumPy slicing, these are written as A[i, :] and A[:, j], except that Python indices start at \(0\) rather than \(1\).

The matrix-vector product \(\mathbf{y}=\mathbf{A}\mathbf{x}\) is defined row by row:

\[\begin{split}\mathbf{A}\mathbf{x} = \begin{bmatrix} \mathbf{A}_{1,:}\\ \mathbf{A}_{2,:}\\ \vdots\\ \mathbf{A}_{m,:} \end{bmatrix} \mathbf{x} = \begin{bmatrix} \mathbf{A}_{1,:}\mathbf{x}\\ \mathbf{A}_{2,:}\mathbf{x}\\ \vdots\\ \mathbf{A}_{m,:}\mathbf{x} \end{bmatrix}.\end{split}\]

Entry \(i\) is

\[y_i = \mathbf{A}_{i,:}\mathbf{x} = \sum_{j=1}^{n} a_{ij}x_j.\]

Thus each entry of \(\mathbf{y}\) is a dot product between one row of \(\mathbf{A}\) and the vector \(\mathbf{x}\). A physical unit check applies row by row: every term \(a_{ij}x_j\) that is added in row \(i\) must have the same units as \(y_i\) and, in a system \(\mathbf{A}\mathbf{x}=\mathbf{b}\), the same units as \(b_i\).

Matrix-vector products: the column view

The same product can be read as a linear combination of columns. Since \(\mathbf{A}_{:,j}\) is column \(j\) of \(\mathbf{A}\),

\[\begin{split}\mathbf{A}\mathbf{x}= \begin{bmatrix}\mathbf{A}_{:,1}&\mathbf{A}_{:,2}&\cdots&\mathbf{A}_{:,n}\end{bmatrix} \begin{bmatrix} x_1\\x_2\\\vdots\\x_n \end{bmatrix} = x_1\mathbf{A}_{:,1}+x_2\mathbf{A}_{:,2}+\cdots+x_n\mathbf{A}_{:,n}.\end{split}\]

The row view emphasizes one output entry at a time. The column view emphasizes how the entries of \(\mathbf{x}\) weight whole column patterns. Both descriptions produce the same vector.

Worked example: a cumulative-sum matrix.

Problem. Let

\[\begin{split}\mathbf{A}=\begin{bmatrix} 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 1 & 0\\ 1 & 1 & 1 & 1 \end{bmatrix}, \qquad \mathbf{x}=\begin{bmatrix}3\\-1\\4\\2\end{bmatrix}.\end{split}\]

Compute \(\mathbf{A}\mathbf{x}\) using both the row view and the column view.

Solution. The row view shows why this matrix forms cumulative sums:

\[\begin{split}\mathbf{A}\mathbf{x} = \begin{bmatrix} 3\\ 3+(-1)\\ 3+(-1)+4\\ 3+(-1)+4+2 \end{bmatrix} = \begin{bmatrix}3\\2\\6\\8\end{bmatrix}.\end{split}\]

The column view gives the same result by adding step patterns:

\[\begin{split}\mathbf{A}\mathbf{x} = 3\begin{bmatrix}1\\1\\1\\1\end{bmatrix} -1\begin{bmatrix}0\\1\\1\\1\end{bmatrix} +4\begin{bmatrix}0\\0\\1\\1\end{bmatrix} +2\begin{bmatrix}0\\0\\0\\1\end{bmatrix} = \begin{bmatrix}3\\2\\6\\8\end{bmatrix}.\end{split}\]

This matrix is not a transport stencil; it is a running-total operator.

The loop pattern behind matrix-vector multiplication

The examples above used small matrices whose entries could be shown on the page. A real engineering matrix may have thousands or millions of entries, so we also care about the computational work. The direct algorithm for a dense matrix-vector product stores all entries of \(\mathbf{A}\) and loops over every row-column pair:

def matvec(A, x):
    m, n = A.shape
    y = np.zeros(m)
    for i in range(m):
        for j in range(n):
            y[i] += A[i, j] * x[j]
    return y

This loop is the computational form of

\[y_i = \sum_{j=1}^n a_{ij}x_j.\]

In practice, we usually write

A = np.array([[2.0, -1.0, 0.0],
              [-1.0, 2.0, -1.0],
              [0.0, -1.0, 2.0]])
x = np.array([1.0, 2.0, 3.0])
y = A @ x

because NumPy calls highly optimized linear algebra routines. The mathematical operation is the same as the nested loop.

Cost of a dense matrix-vector product

For an \(n\times n\) dense matrix, each of the \(n\) output entries requires approximately \(n\) multiplications and \(n-1\) additions. The total cost is therefore roughly

\[n^2+n(n-1)=2n^2-n\]

floating-point operations. For large \(n\), we summarize this by saying the cost is \(O(n^2)\).

Students can see this cost experimentally in a Jupyter notebook. The exact timings depend on the machine, but the comparison is usually dramatic:

import numpy as np

def matvec(A, x):
    m, n = A.shape
    y = np.zeros(m)
    for i in range(m):
        for j in range(n):
            y[i] += A[i, j] * x[j]
    return y

n = 300
rng = np.random.default_rng(1)
A = rng.random((n, n))
x = rng.random(n)

# In Jupyter/IPython:
%timeit A @ x
%timeit matvec(A, x)

Both lines compute the same mathematical product. The second version makes the \(n^2\) loop structure visible, while the first uses optimized compiled routines.

This cost estimate matters because real engineering matrices can be very large. A dense \(1000\times 1000\) matrix has one million entries. A dense \(100000\times 100000\) matrix has ten billion entries, which is already too large to treat casually. Lesson 14 uses this cost argument to motivate sparse storage and iterative methods.

Matrix-matrix products

A matrix-matrix product can be understood as several matrix-vector products at once. If

\[\mathbf{X}=\begin{bmatrix}\mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_p\end{bmatrix},\]

then

\[\mathbf{A}\mathbf{X}=\begin{bmatrix} \mathbf{A}\mathbf{x}_1 & \mathbf{A}\mathbf{x}_2 & \cdots & \mathbf{A}\mathbf{x}_p \end{bmatrix}.\]

This is useful when the same matrix acts on several vectors. In Lesson 12, the same idea appears when the same coefficient matrix is reused for several right-hand sides.

Entry by entry,

\[c_{ij}=\sum_{k} a_{ik} b_{kj}.\]

Thus each entry of \(\mathbf{C}=\mathbf{A}\mathbf{B}\) is a dot product between a row of \(\mathbf{A}\) and a column of \(\mathbf{B}\).

The same physical unit rule applies inside each dot product. The products \(a_{ik}b_{kj}\) must be compatible before they are added. When the columns of \(\mathbf{B}\) are several right-hand sides, they should represent scenarios with the same physical units unless the units are being tracked deliberately.

Array shapes and noncommutativity

Many linear algebra errors are shape errors. If \(\mathbf{A}\) is \(m\times n\), then \(\mathbf{A}\mathbf{x}\) is defined only when \(\mathbf{x}\in\mathbb{R}^n\). The result is a vector in \(\mathbb{R}^m\).

The same shape rule applies to matrix-matrix products. If \(\mathbf{A}\) is \(m\times n\) and \(\mathbf{B}\) is \(n\times p\), then \(\mathbf{A}\mathbf{B}\) is defined and has shape \(m\times p\). The reverse product \(\mathbf{B}\mathbf{A}\) may not even be defined. Even when both products are defined, they are usually not equal.

For example,

\[\begin{split}\mathbf{A}= \begin{bmatrix} 1&2\\ 0&1 \end{bmatrix}, \qquad \mathbf{B}= \begin{bmatrix} 1&0\\ 3&1 \end{bmatrix}.\end{split}\]

Then

\[\begin{split}\mathbf{A}\mathbf{B}= \begin{bmatrix} 7&2\\ 3&1 \end{bmatrix}, \qquad \mathbf{B}\mathbf{A}= \begin{bmatrix} 1&2\\ 3&7 \end{bmatrix}.\end{split}\]

Both products exist, but they are different. Matrix multiplication is not commutative.

Transposes

The transpose swaps rows and columns:

\[(\mathbf{A}^{\mathsf{T}})_{ij}=a_{ji}.\]

If \(\mathbf{A}\) is \(m\times n\), then \(\mathbf{A}^{\mathsf{T}}\) is \(n\times m\). Transposes appear in dot products, least-squares problems, orthogonality, symmetric matrices, and eigenvalue computations.

For example, if

\[\begin{split}\mathbf{A}= \begin{bmatrix} 1&2&3\\ 4&5&6\\ 7&8&9 \end{bmatrix},\end{split}\]

then

\[\begin{split}\mathbf{A}^{\mathsf{T}}= \begin{bmatrix} 1&4&7\\ 2&5&8\\ 3&6&9 \end{bmatrix}.\end{split}\]

The first row of \(\mathbf{A}\) becomes the first column of \(\mathbf{A}^{\mathsf{T}}\), and similarly for the other rows.

Permutation matrices

Some matrices perform simple bookkeeping operations. A permutation matrix reorders entries. For example,

\[\begin{split}\mathbf{P}= \begin{bmatrix} 0&1&0\\ 1&0&0\\ 0&0&1 \end{bmatrix}, \qquad \mathbf{x}= \begin{bmatrix} 4\\5\\6 \end{bmatrix}.\end{split}\]

Then

\[\begin{split}\mathbf{P}\mathbf{x}= \begin{bmatrix} 5\\4\\6 \end{bmatrix}.\end{split}\]

The first two entries have been swapped. This matters in practical linear solvers because pivoting swaps rows to put better pivot entries in useful positions.

The same permutation matrix swaps rows when it multiplies a matrix on the left. If

\[\begin{split}\mathbf{A}= \begin{bmatrix} 1&2&3\\ 4&5&6\\ 7&8&9 \end{bmatrix},\end{split}\]

then

\[\begin{split}\mathbf{P}\mathbf{A}= \begin{bmatrix} 4&5&6\\ 1&2&3\\ 7&8&9 \end{bmatrix}.\end{split}\]

In numerical linear solves, reordering the equations can sometimes make the computed solution more accurate. Pivoting is the systematic version of this row-reordering idea.

Residuals revisited

The residual from Lesson 10 is

\[\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}.\]

Now we can read this expression operationally:

• compute \(\mathbf{A}\mathbf{x}\), which is one dot product per row,

• subtract the result from \(\mathbf{b}\), entry by entry,

• interpret each component of \(\mathbf{r}\) as the defect in one equation.

A residual check is therefore not mysterious. It is a set of defects in the balance equations computed after the solve.

Exercises

1. Let

   \[\begin{split}\mathbf{A}=\begin{bmatrix}1&0&0\\1&1&0\\1&1&1\end{bmatrix}, \qquad \mathbf{x}=\begin{bmatrix}2\\-3\\5\end{bmatrix}.\end{split}\]

   Compute \(\mathbf{A}\mathbf{x}\).

2. For the same \(\mathbf{A}\), identify \(\mathbf{A}_{2,:}\) and \(\mathbf{A}_{:,2}\). Then write the NumPy slices that would extract them.

3. Compute the transpose of

   \[\begin{split}\mathbf{C}=\begin{bmatrix}1&2&3\\4&5&6\end{bmatrix}.\end{split}\]

4. Let

   \[\begin{split}\mathbf{P}=\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\end{bmatrix}, \qquad \mathbf{x}=\begin{bmatrix}10\\20\\30\end{bmatrix}.\end{split}\]

   Compute \(\mathbf{P}\mathbf{x}\).

Answers.

1. \(\mathbf{A}\mathbf{x}=\begin{bmatrix}2\\-1\\4\end{bmatrix}\).

2. \(\mathbf{A}_{2,:}=\begin{bmatrix}1&1&0\end{bmatrix}\), \(\mathbf{A}_{:,2}=\begin{bmatrix}0\\1\\1\end{bmatrix}\). In NumPy, these are A[1, :] and A[:, 1].

3. \(\mathbf{C}^{\mathsf{T}}=\begin{bmatrix}1&4\\2&5\\3&6\end{bmatrix}\).

4. \(\mathbf{P}\mathbf{x}=\begin{bmatrix}30\\20\\10\end{bmatrix}\).

End-of-lecture summary

Takeaway.

• A dot product is a weighted sum and later becomes a geometric measure of alignment.

• A matrix-vector product is many dot products, one for each row of the matrix; each row can be read as one balance equation.

• The column view says \(\mathbf{A}\mathbf{x}\) is a linear combination of the columns of \(\mathbf{A}\), which can be read as influence or contribution patterns.

• Dense matrix-vector products cost \(O(n^2)\), which motivates sparse methods for large models.

• Matrix-matrix multiplication collects many related matrix-vector products into one operation, preparing the many-right-hand-side solves in Lesson 12.

• Lesson 12 uses these operations to explain elimination, triangular solves, and reusable factorizations.