At a reductive level, a matrix is a table of numbers that obeys certain algebraic laws. But matrices are pervasive in scientific computation, mainly because they represent linear operations on vectors. Moreover, vectors go far beyond the three-dimensional representations of physical quantities you learned about in calculus.
2.2.1 Notation ¶ We use capital letters in bold to refer to matrices, and lowercase bold letters for vectors. All named vectors in this book are column vectors . The bold symbol 0 \boldsymbol{0} 0
To refer to a specific element of a matrix, we use the uppercase name of the matrix without boldface, as in A 24 A_{24} A 24 ( 2 , 4 ) (2,4) ( 2 , 4 ) A \mathbf{A} A [1] x 3 x_3 x 3
We will have frequent need to refer to the individual columns of a matrix as vectors. Our convention is to use a lowercase bold version of the matrix name with a subscript to represent the column number. Thus, a 1 , a 2 , … , a n \mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_n a 1 , a 2 , … , a n m × n m\times n m × n A \mathbf{A} A v 1 , … , v p \mathbf{v}_1,\ldots,\mathbf{v}_p v 1 , … , v p V \mathbf{V} V V = [ v j ] \mathbf{V}=\bigl[ \mathbf{v}_j \bigr] V = [ v j ]
The notation A T \mathbf{A}^T A T adjoint A ∗ \mathbf{A}^* A ∗ [2] A \mathbf{A} A A ∗ = A T \mathbf{A}^*=\mathbf{A}^T A ∗ = A T symmetric matrix is a square matrix such that A T = A \mathbf{A}^T=\mathbf{A} A T = A
The identity matrix of size n n n I \mathbf{I} I I n \mathbf{I}_n I n e j \mathbf{e}_j e j
2.2.2 Block matrix expressions ¶ We will often find it useful to break a matrix into separately named pieces. For example, we might write
A = [ A 11 A 12 A 13 A 21 A 22 A 23 ] , B = [ B 1 B 2 B 3 ] . \mathbf{A} =
\begin{bmatrix}
\mathbf{A}_{11} & \mathbf{A}_{12} & \mathbf{A}_{13} \\
\mathbf{A}_{21} & \mathbf{A}_{22} & \mathbf{A}_{23}
\end{bmatrix}, \qquad
\mathbf{B} =
\begin{bmatrix}
\mathbf{B}_1 \\ \mathbf{B}_2 \\ \mathbf{B}_3
\end{bmatrix}. A = [ A 11 A 21 A 12 A 22 A 13 A 23 ] , B = ⎣ ⎡ B 1 B 2 B 3 ⎦ ⎤ . It’s understood that blocks that are on top of one another have the same number of columns, and blocks that are side by side have the same number of rows. Typically, if the blocks all have compatible dimensions, then they can be multiplied as though the blocks were scalars. For instance, continuing with the definitions above, we say that A \mathbf{A} A 2 × 3 2\times 3 2 × 3 B \mathbf{B} B 3 × 1 3\times 1 3 × 1
A B = [ A 11 B 1 + A 12 B 2 + A 13 B 3 A 21 B 1 + A 22 B 2 + A 23 B 3 ] , \mathbf{A} \mathbf{B} =
\begin{bmatrix}
\mathbf{A}_{11}\mathbf{B}_1 + \mathbf{A}_{12}\mathbf{B}_2 + \mathbf{A}_{13}\mathbf{B}_3 \\
\mathbf{A}_{21}\mathbf{B}_1 + \mathbf{A}_{22}\mathbf{B}_2 + \mathbf{A}_{23}\mathbf{B}_3
\end{bmatrix}, AB = [ A 11 B 1 + A 12 B 2 + A 13 B 3 A 21 B 1 + A 22 B 2 + A 23 B 3 ] , provided that the individual block products are well-defined. For transposes we have, for example,
A T = [ A 11 T A 21 T A 12 T A 22 T A 13 T A 23 T ] . \mathbf{A}^T =
\begin{bmatrix}
\mathbf{A}_{11}^T & \mathbf{A}_{21}^T \\[2mm]
\mathbf{A}_{12}^T & \mathbf{A}_{22}^T \\[2mm]
\mathbf{A}_{13}^T & \mathbf{A}_{23}^T
\end{bmatrix}. A T = ⎣ ⎡ A 11 T A 12 T A 13 T A 21 T A 22 T A 23 T ⎦ ⎤ . 2.2.3 Vector and matrix basics ¶ Vectors and matrices are integral to scientific computing. All modern languages provide ways to work with them beyond manipulation of individual elements.
In Julia, vectors and matrices are one-dimensional and two-dimensional arrays, respectively. Square brackets are used to enclose elements of a matrix or vector. Use spaces for horizontal concatenation, and semicolons or new lines to indicate vertical concatenation.
The size function returns the number of rows and columns in a matrix. Use length to get the number of elements in a vector or matrix.
A = [ 1 2 3 4 5; 50 40 30 20 10
π sqrt(2) exp(1) (1+sqrt(5))/2 log(3) ]
A vector is not quite the same thing as a matrix: it has only one dimension, not two. Separate its elements by commas or semicolons:
x = [ 3, 3, 0, 1, 0 ]
size(x)
For some purposes, however, an n n n
y = [ 3 3 0 1 0 ]
size(y)
This 1 × 5 1\times 5 1 × 5
Concatenated elements within brackets may be matrices or vectors for a block representation, as long as all the block sizes are compatible.
The zeros and ones functions construct matrices with entries all zero or one, respectively.
B = [ zeros(3, 2) ones(3, 1) ]
A single quote ' after a matrix returns its adjoint. For real matrices, this is the transpose; for complex-valued matrices, the elements are also conjugated.
If x is simply a vector, then its transpose has a row shape.
There are many convenient shorthand ways of building vectors and matrices other than entering all of their entries directly or in a loop. To get a range with evenly spaced entries between two endpoints, you have two options. One is to use a colon :.
z = 0:3:12 # start:step:stop
(Ranges are not strictly considered vectors, but they behave identically in most circumstances.) Instead of specifying the step size, you can give the number of points in the range if you use range.
Accessing an element is done by giving one (for a vector) or two (for a matrix) index values within square brackets.
The end keyword refers to the last element in a dimension. It saves you from having to compute and store the size of the matrix first.
The indices can be vectors or ranges, in which case a block of the matrix is accessed.
A[1:2, end-2:end] # first two rows, last three columns
If a dimension has only the index : (a colon), then it refers to all the entries in that dimension of the matrix.
A[:, 1:2:end] # all of the odd columns
The matrix and vector senses of addition, subtraction, scalar multiplication, multiplication, and power are all handled by the usual symbols.
Use diagm to construct a matrix by its diagonals. A more general syntax puts elements on super- or subdiagonals.
B = diagm( [-1, 0, -5] ) # create a diagonal matrix
@show size(A), size(B);
BA = B * A # matrix product
A * B causes an error here, because the dimensions aren’t compatible.
Errors are formally called exceptions in Julia.
A square matrix raised to an integer power is the same as repeated matrix multiplication.
Sometimes one instead wants to treat a matrix or vector as a mere array and simply apply a single operation to each element of it. For multiplication, division, and power, the corresponding operators start with a dot.
Because both matrices are 3 × 5 3\times 5 3 × 5 A * C would be an error here, but elementwise operations are fine.
The two operands of a dot operator have to have the same size—unless one is a scalar, in which case it is expanded or broadcast to be the same size as the other operand.
Most of the mathematical functions, such as cos, sin, log, exp, and sqrt, expect scalars as operands. However, you can broadcast any function, including ones that you have defined, across a vector or array by using a special dot syntax.
A dot added to the end of a function name means to apply the function elementwise to an array.
show(cos.(π*x)) # broadcast to a function
Rather than dotting multiple individual functions, you can use @. before an expression to broadcast everything within it.
show(@. cos(π*(x+1)^3)) # broadcast an entire expression
In MATLAB, every numerical value is treated like a matrix. A matrix with one row or one column is interpreted as a vector, and a 1 × 1 1\times 1 1 × 1
Square brackets are used to enclose elements of a matrix or vector. Use spaces for horizontal concatenation, and semicolons or new lines to indicate vertical concatenation.
The size function returns the number of rows and columns in a matrix. Use length to get the number of elements in a vector or matrix.
A = [
1 2 3 4 5;
50 40 30 20 10
pi sqrt(2) exp(1) (1+sqrt(5))/2 log(3)
]
x = [ 3, 3, 0, 1, 0 ]; % row vector
size(x)
Concatenated elements within brackets may be matrices or vectors for a block representation, as long as all the block sizes are compatible.
The zeros and ones functions construct matrices with entries all zero or one, respectively.
B = [ zeros(3, 2) ones(3, 1) ]
A single quote ' after a matrix returns its adjoint. For real matrices, this is the transpose; for complex-valued matrices, the elements are also conjugated.
There are many convenient shorthand ways of building vectors and matrices other than entering all of their entries directly or in a loop. To get a range with evenly spaced entries between two endpoints, you have two options. One is to use a colon :.
z = 0:3:12 % start:step:stop
Instead of specifying the step size, you can give the number of points in the range if you use linspace.
s = linspace(-1, 1, 5) % row result
Accessing an element is done by giving one (for a vector) or two (for a matrix) index values within parentheses.
The end keyword refers to the last element in a dimension. It saves you from having to compute and store the size of the matrix first.
The indices can be vectors or ranges, in which case a block of the matrix is accessed.
A(1:2, end-2:end) % first two rows, last three columns
If a dimension has only the index : (a colon), then it refers to all the entries in that dimension of the matrix.
A(:, 1:2:end) % all of the odd columns
The matrix and vector senses of addition, subtraction, scalar multiplication, multiplication, and power are all handled by the usual symbols.
Use diag to construct a matrix by its diagonals. A more general syntax puts elements on super- or subdiagonals.
B = diag([-1, 0, -5]) % create a diagonal matrix
BA = B * A % matrix product
A * B causes an error here, because the dimensions aren’t compatible.
Errors are formally called exceptions in Julia.
A square matrix raised to an integer power is the same as repeated matrix multiplication.
Sometimes one instead wants to treat a matrix or vector as a mere array and simply apply a single operation to each element of it. For multiplication, division, and power, the corresponding operators start with a dot.
Because both matrices are 3 × 5 3\times 5 3 × 5 A * C would be an error here, but elementwise operations are fine.
The two operands of a dot operator have to have the same size—unless one is a scalar, in which case it is expanded or broadcast to be the same size as the other operand.
Most of the mathematical functions, such as cos, sin, log, exp, and sqrt, can operate elementwise on vectors and matrices.
A vector is created using square brackets and commas to enclose and separate its entries.
x = array([3, 3, 0, 1, 0 ])
print(x.shape)
To construct a matrix, you nest the brackets to create a “vector of vectors”. The inner vectors are the rows.
A = array([
[1, 2, 3, 4, 5],
[50, 40, 30, 20, 10],
[pi, sqrt(2), exp(1), (1+sqrt(5))/2, log(3)]
])
print(A)
print(A.shape)
In this text, we treat all vectors as equivalent to matrices with a single column. That isn’t true in NumPy, because even an n × 1 n \times 1 n × 1
array([[3], [1], [2]]).shape
You can concatenate arrays with compatible dimensions using hstack and vstack.
Transposing a matrix is done by appending .T to it.
For matrices with complex values, we usually want instead the adjoint or hermitian, which is .conj().T.
There are many convenient shorthand ways of building vectors and matrices other than entering all of their entries directly or in a loop. To get a vector with evenly spaced entries between two endpoints, you have two options.
print(arange(1, 7, 2)) # from 1 to 7 (not inclusive), step by 2
print(linspace(-1, 1, 5)) # from -1 to 1 (inclusive), with 5 total values
The practical difference between these is whether you want to specify the step size in arange or the number of points in linspace.
Accessing an element is done by giving one (for a vector) or two index values in square brackets. In Python, indexing always starts with zero, not 1.
A = array([
[1, 2, 3, 4, 5],
[50, 40, 30, 20, 10],
linspace(-5, 5, 5)
])
x = array([3, 2, 0, 1, -1 ])
print("row 2, col 3 of A:", A[1, 2])
print("first element of x:", x[0])
The indices can be ranges, in which case a slice or block of the matrix is accessed. You build these using a colon in the form start:stop. However, the last value of this range is stop-1, not stop.
print(A[1:3, 0:2]) # rows 2 and 3, cols 1 and 2
If start or stop is omitted, the range extends to the first or last index.
print(x[1:]) # elements 2 through the end
print(A[:2, 0]) # first two rows in column 1
Notice in the last case above that even when the slice is in the shape of a column vector, the result is just a vector with one dimension and neither row nor column shape.
There are more variations on the colon ranges. A negative value means to count from the end rather than the beginning. And a colon by itself means to include everything from the relevant dimension.
print(A[:-1, :]) # all rows up to the last, all columns
Finally, start:stop:step means to step size or stride other than one. You can mix this with the other variations.
print(x[::2]) # all the odd indexes
print(A[:, ::-1]) # reverse the columns
The matrix and vector senses of addition, subtraction, and scalar multiplication and division are all handled by the usual symbols. Two matrices of the same size (what NumPy calls shape) are operated on elementwise.
print(A - 2 * ones([3, 5])) # subtract two from each element
If one operand has a smaller number of dimensions than the other, Python tries to broadcast it in the “missing” dimension(s), and the operation proceeds if the resulting shapes are identical.
print(A - 2) # subtract two from each element
u = array([1, 2, 3, 4, 5])
print(A - u) # repeat this row for every row of A
v = array([1, 2, 3])
print(A - v) # broadcasting this would be 3x3, so it's an error
print(A - v.reshape([3, 1])) # broadcasts to each column of A
Matrix–matrix and matrix–vector products are computed using @ or matmul.
B = diag([-1, 0, -5]) # create a diagonal 3x3
print(B @ A) # matrix product
A B AB A B
print(A @ B) # incompatible sizes
The multiplication operator * is reserved for elementwise multiplication. Both operands have to be the same size, after any potential broadcasts.
print(B * A) # not the same size, so it's an error
print((A / 2) * A) # elementwise
To raise to a power elementwise, use a double star. This will broadcast as well.
Most of the mathematical functions, such as cos, sin, log, exp and sqrt, expecting scalars as operands will be broadcast to arrays.
2.2.4 Row and column operations ¶ A critical identity in matrix multiplication is
A e j = a j . \mathbf{A} \mathbf{e}_j = \mathbf{a}_j. A e j = a j . Furthermore, the expression
A [ e 1 e 3 e 5 ] \mathbf{A}
\begin{bmatrix}
\mathbf{e}_1 & \mathbf{e}_3 & \mathbf{e}_5
\end{bmatrix} A [ e 1 e 3 e 5 ] reproduces three columns. An equivalent expression in Julia would be A[:,1:2:5].
We can extend the same idea to rows by using the general identity ( R S ) T = S T R T (\mathbf{R}\mathbf{S})^T=\mathbf{S}^T\mathbf{R}^T ( RS ) T = S T R T B = A T \mathbf{B}=\mathbf{A}^T B = A T [ b j ] \bigl[ \mathbf{b}_j \bigr] [ b j ]
( b j ) T = ( B e j ) T = e j T B T = e j T A . (\mathbf{b}_j)^T = (\mathbf{B} \mathbf{e}_j)^T = \mathbf{e}_j^T \mathbf{B}^T = \mathbf{e}_j^T \mathbf{A}. ( b j ) T = ( B e j ) T = e j T B T = e j T A . But e j T \mathbf{e}_j^T e j T j j j row of I \mathbf{I} I b j T \mathbf{b}_j^T b j T j j j B \mathbf{B} B j j j row of A \mathbf{A} A B = A T \mathbf{B}=\mathbf{A}^T B = A T multiplication on the left by row j j j j j j . Extracting the single element ( i , j ) (i,j) ( i , j ) e i T A e j \mathbf{e}_i^T \mathbf{A} \mathbf{e}_j e i T A e j
Being able to extract specific rows and columns of a matrix via algebra makes it straightforward to do row- and column-oriented operations, such as linear combinations.
Say that A \mathbf{A} A A \mathbf{A} A
A ( e 1 + 2 e 3 ) . \mathbf{A}(\mathbf{e}_1+2\mathbf{e}_3). A ( e 1 + 2 e 3 ) . Suppose we want to do this operation “in place,” meaning replacing the first column of A \mathbf{A} A A \mathbf{A} A A \mathbf{A} A
A [ e 1 + 2 e 3 e 2 e 3 e 4 e 5 ] . \mathbf{A}
\begin{bmatrix}
\mathbf{e}_1+2\mathbf{e}_3 & \mathbf{e}_2 & \mathbf{e}_3 & \mathbf{e}_4 & \mathbf{e}_5
\end{bmatrix}. A [ e 1 + 2 e 3 e 2 e 3 e 4 e 5 ] . The Julia equivalent is
The += operator means to increment the item on the left-hand side. There are similar interpretations for -= and *=.
The MATLAB equivalent is
A(:, 1) = A(:, 1) + 2 * A(:, 3)The NumPy equivalent is
The += operator means to increment the item on the left-hand side. There are similar interpretations for -= and *=
2.2.5 Exercises ¶ ✍ Suppose
C = [ I A − I B ] . \mathbf{C} =
\begin{bmatrix}
\mathbf{I} & \mathbf{A} \\ -\mathbf{I} & \mathbf{B}
\end{bmatrix}. C = [ I − I A B ] . Using block notation, find C 2 \mathbf{C}^2 C 2 C 3 \mathbf{C}^3 C 3
⌨ Let
A = [ 2 1 1 0 0 − 1 4 1 2 2 0 − 2 1 3 − 1 5 ] , B = [ 3 − 1 0 2 7 1 0 2 ] , \mathbf{A} =
\begin{bmatrix}
2 & 1 & 1 & 0 \\ 0 & -1 & 4 & 1 \\ 2 & 2 & 0 & -2 \\ 1 & 3 & -1
& 5
\end{bmatrix},
\quad
\mathbf{B} =
\begin{bmatrix}
3 & -1 & 0 & 2 \\ 7 & 1 & 0 & 2
\end{bmatrix}, A = ⎣ ⎡ 2 0 2 1 1 − 1 2 3 1 4 0 − 1 0 1 − 2 5 ⎦ ⎤ , B = [ 3 7 − 1 1 0 0 2 2 ] , u = [ 2 − 1 3 1 ] , v = [ π e ] . \mathbf{u} =
\begin{bmatrix}
2 \\ -1 \\ 3 \\ 1
\end{bmatrix},
\quad
\mathbf{v} =
\begin{bmatrix}
\pi \\ e
\end{bmatrix}. u = ⎣ ⎡ 2 − 1 3 1 ⎦ ⎤ , v = [ π e ] . (Do not round off the values in v \mathbf{v} v
(a) A B , \mathbf{A}\mathbf{B},\quad AB , (b) B A , \mathbf{B} \mathbf{A},\quad BA , (c) v T B , \mathbf{v}^T \mathbf{B},\quad v T B , (d) B u , \mathbf{B} \mathbf{u},\quad Bu , (e) [ u A u A 2 u A 3 u ] \bigl[ \, \mathbf{u}\:\: \mathbf{A}\mathbf{u} \:\: \mathbf{A}^2 \mathbf{u} \:\: \mathbf{A}^3 \mathbf{u} \bigr] [ u Au A 2 u A 3 u ]
⌨ Let
u = [ 1 3 5 7 9 11 ] , v = [ − 60 − 50 − 40 − 30 − 20 − 10 ] . \mathbf{u} =
\begin{bmatrix}
1\\3\\5\\7\\9\\11
\end{bmatrix}, \qquad
\mathbf{v} =
\begin{bmatrix}
-60 \\ -50 \\ -40 \\ -30 \\ -20 \\ -10
\end{bmatrix}. u = ⎣ ⎡ 1 3 5 7 9 11 ⎦ ⎤ , v = ⎣ ⎡ − 60 − 50 − 40 − 30 − 20 − 10 ⎦ ⎤ . Find the inner products u T v \mathbf{u}^T\mathbf{v} u T v v T u \mathbf{v}^T\mathbf{u} v T u u v T \mathbf{u}\mathbf{v}^T u v T v u T \mathbf{v}\mathbf{u}^T v u T
⌨ In Julia, give a demonstration of the identity ( A B ) T = B T A T (\mathbf{A}\mathbf{B})^T=\mathbf{B}^T\mathbf{A}^T ( AB ) T = B T A T 3 × 4 3\times 4 3 × 4 A \mathbf{A} A 4 × 2 4\times 2 4 × 2 B \mathbf{B} B
✍ Prove that if A \mathbf{A} A B \mathbf{B} B ( A B ) − 1 = B − 1 A − 1 (\mathbf{A}\mathbf{B})^{-1}=\mathbf{B}^{-1}\mathbf{A}^{-1} ( AB ) − 1 = B − 1 A − 1 A B \mathbf{A}\mathbf{B} AB
✍ Suppose B \mathbf{B} B 4 × 3 4\times 3 4 × 3 A \mathbf{A} A B \mathbf{B} B A \mathbf{A} A B \mathbf{B} B
(a) A ∈ R 4 × 1 \mathbf{A}\in\mathbb{R}^{4 \times 1} A ∈ R 4 × 1 B \mathbf{B} B B \mathbf{B} B
(b) The rows of A ∈ R 4 × 3 \mathbf{A}\in\mathbb{R}^{4 \times 3} A ∈ R 4 × 3 B \mathbf{B} B
(c) The first column of A ∈ R 4 × 3 \mathbf{A}\in\mathbb{R}^{4 \times 3} A ∈ R 4 × 3 B \mathbf{B} B A \mathbf{A} A B \mathbf{B} B A \mathbf{A} A B \mathbf{B} B
(d) A A A B \mathbf{B} B
(a) ✍ Prove that for real vectors v \mathbf{v} v w \mathbf{w} w v T w \mathbf{v}^T\mathbf{w} v T w w T v \mathbf{w}^T\mathbf{v} w T v
(b) ✍ Prove true, or give a counterexample for, the equivalent statement about outer products, v w T \mathbf{v}\mathbf{w}^T v w T w v T \mathbf{w}\mathbf{v}^T w v T