Computing with matrices
import Pkg; Pkg.activate("/Users/driscoll/Documents/GitHub/fnc")
using FNCFunctions
using Plots
default(
titlefont=(11,"Helvetica"),
guidefont=(11,"Helvetica"),
linewidth = 2,
markersize = 3,
msa = 0,
size=(500,320),
label="",
html_output_format = "svg"
)
using PrettyTables, LaTeXStrings, Printf
using LinearAlgebra Activating project at `~/Documents/GitHub/fnc`
2.2. Computing with matrices ¶ 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 may refer to a vector of all zeros or to a zero matrix, depending on context; we use 0 as the scalar zero only.
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 to mean the ( 2 , 4 ) (2,4) ( 2 , 4 ) element of A \mathbf{A} A .[1] To refer to an element of a vector, we use just one subscript, as in x 3 x_3 x 3 . If you see a boldface character with one or more subscripts, then you know that it is a matrix or vector that belongs to a sequence or indexed collection.
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 are the columns of the m × n m\times n m × n matrix A \mathbf{A} A . Conversely, whenever we define a sequence of vectors v 1 , … , v p \mathbf{v}_1,\ldots,\mathbf{v}_p v 1 , … , v p , we can implicitly consider them to be columns of a matrix V \mathbf{V} V . Sometimes we might write V = [ v j ] \mathbf{V}=\bigl[ \mathbf{v}_j \bigr] V = [ v j ] to emphasize the connection.
The notation A T \mathbf{A}^T A T is used for the transpose of a matrix, in which the rows and columns switch places. In the case of complex matrices, the complex conjugate[2] becomes involved with this operation.
The adjoint or hermitian of a matrix, denoted, A ∗ \mathbf{A}^* A ∗ , is the elementwise complex conjugate of the transpose of A \mathbf{A} A .
If A \mathbf{A} A is real, then A ∗ = A T \mathbf{A}^*=\mathbf{A}^T A ∗ = A T .
The identity matrix of size n n n is denoted I \mathbf{I} I , or sometimes I n \mathbf{I}_n I n , is an n × n n\times n n × n matrix with ones on the main diagonal and zeros elsewhere. It is the multiplicative identity for matrix multiplication. For columns of the identity, we break with our usual naming convention and denote them by 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 is block-2 × 3 2\times 3 2 × 3 and B \mathbf{B} B is block-3 × 1 3\times 1 3 × 1 , so we can write
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) ]3×5 Matrix{Float64}:
1.0 2.0 3.0 4.0 5.0
50.0 40.0 30.0 20.0 10.0
3.14159 1.41421 2.71828 1.61803 1.09861
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 -vector in Julia is treated like having a column shape. Note the difference if we use spaces instead of commas inside the brackets:
y = [ 3 3 0 1 0 ]
size(y)This 1 × 5 1\times 5 1 × 5 matrix is not equivalent to a vector.
Concatenated elements within brackets may be matrices or vectors for a block representation, as long as all the block sizes are compatible.
5×2 Matrix{Int64}:
3 3
3 3
0 0
1 1
0 0
10-element Vector{Int64}:
3
3
0
1
0
3
3
0
1
0
The zeros and ones functions construct matrices with entries all zero or one, respectively.
B = [ zeros(3, 2) ones(3, 1) ]3×3 Matrix{Float64}:
0.0 0.0 1.0
0.0 0.0 1.0
0.0 0.0 1.0
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.
5×3 adjoint(::Matrix{Float64}) with eltype Float64:
1.0 50.0 3.14159
2.0 40.0 1.41421
3.0 30.0 2.71828
4.0 20.0 1.61803
5.0 10.0 1.09861
If x is simply a vector, then its transpose has a row shape.
1×5 adjoint(::Vector{Int64}) with eltype Int64:
3 3 0 1 0
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.
┌ Warning: Timed out waiting for `Base.active_repl_backend.ast_transforms` to become available. Autoloads will not work.
└ @ BasicAutoloads ~/.julia/packages/BasicAutoloads/08hIo/src/BasicAutoloads.jl:117
[ Info: If you have a slow startup file, consider moving `register_autoloads` to the end of it.
20.0
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 columns2×3 Matrix{Float64}:
3.0 4.0 5.0
30.0 20.0 10.0
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 columns3×3 Matrix{Float64}:
1.0 3.0 5.0
50.0 30.0 10.0
3.14159 2.71828 1.09861
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 matrix3×3 Matrix{Int64}:
-1 0 0
0 0 0
0 0 -5
@show size(A), size(B);
BA = B * A # matrix product(size(A), size(B)) = ((3, 5), (3, 3)) 3×5 Matrix{Float64}:
-1.0 -2.0 -3.0 -4.0 -5.0
0.0 0.0 0.0 0.0 0.0
-15.708 -7.07107 -13.5914 -8.09017 -5.49306
A * B causes an error here, because the dimensions aren’t compatible.
Errors are formally called exceptions in Julia.
DimensionMismatch: incompatible dimensions for matrix multiplication: tried to multiply a matrix of size (3, 5) with a matrix of size (3, 3). The second dimension of the first matrix: 5, does not match the first dimension of the second matrix: 3.
Stacktrace:
[1] matmul_size_check_error ( sizeA :: Tuple {Int64, Int64} , sizeB :: Tuple {Int64, Int64} )
@ LinearAlgebra ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:431
[2] matmul_size_check
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:418 [inlined]
[3] matmul_size_check
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:441 [inlined]
[4] _generic_matmatmul! ( C :: Matrix {Float64} , A :: Matrix {Float64} , B :: Matrix {Int64} , alpha :: Bool, beta :: Bool )
@ LinearAlgebra ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:1020
[5] generic_matmatmul!
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:1011 [inlined]
[6] generic_matmatmul_wrapper!
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:343 [inlined]
[7] _mul!
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:328 [inlined]
[8] mul!
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:297 [inlined]
[9] mul!
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:265 [inlined]
[10] mul
@ ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:118 [inlined]
[11] * ( A :: Matrix {Float64} , B :: Matrix {Int64} )
@ LinearAlgebra ~/.julia/juliaup/julia-1.12.3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.12/LinearAlgebra/src/ matmul.jl:114
[12] top-level scope
@ In[20]:1
[13] eval ( m :: Module, e :: Any )
@ Core ./ boot.jl:489 A square matrix raised to an integer power is the same as repeated matrix multiplication.
3×3 Matrix{Int64}:
-1 0 0
0 0 0
0 0 -125
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.
3×5 Matrix{Float64}:
-1.0 -4.0 -9.0 -16.0 -25.0
-2500.0 -1600.0 -900.0 -400.0 -100.0
-9.8696 -2.0 -7.38906 -2.61803 -1.20695
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.
5-element Vector{Int64}:
9
9
0
1
0
5-element Vector{Int64}:
8
8
1
2
1
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[-1.0, -1.0, 1.0, -1.0, 1.0] Rather than dotting multiple individual functions, you can use @. before an expression to broadcast everything within it.
show(@. cospi( (x + 1)^3) ) # broadcast an entire expression[1.0, 1.0, -1.0, 1.0, -1.0] 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 . Let B = A T \mathbf{B}=\mathbf{A}^T B = A T have columns [ b j ] \bigl[ \mathbf{b}_j \bigr] [ b j ] , and note
( 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 is the j j j th row of I \mathbf{I} I , and b j T \mathbf{b}_j^T b j T is the transpose of the j j j th column of B \mathbf{B} B , which is the j j j th row of A \mathbf{A} A by B = A T \mathbf{B}=\mathbf{A}^T B = A T . Thus, multiplication on the left by row j j j of the identity extracts the j j j th row . Extracting the single element ( i , j ) (i,j) ( i , j ) from the matrix is, therefore, 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 has five columns. Adding twice the third column of A \mathbf{A} A to its first column is done by
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 with this value and leaving the other four columns of A \mathbf{A} A alone. We can replace A \mathbf{A} A with
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 *=.
2.2.5 Exercises ¶ ⌨ 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 —find them using native commands.) For each expression below, use the computer to find the result, or explain why the result does not exist.
(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 ] .
✍ Prove that if A \mathbf{A} A and B \mathbf{B} B are invertible and have the same size, then ( A B ) − 1 = B − 1 A − 1 (\mathbf{A}\mathbf{B})^{-1}=\mathbf{B}^{-1}\mathbf{A}^{-1} ( AB ) − 1 = B − 1 A − 1 . (In producing the inverse, it follows that A B \mathbf{A}\mathbf{B} AB is invertible as well.)
✍ Suppose B \mathbf{B} B is an arbitrary 4 × 3 4\times 3 4 × 3 matrix. In each part below a matrix A \mathbf{A} A is described in terms of B \mathbf{B} B . Express A \mathbf{A} A as a product of B \mathbf{B} B with one or more other matrices.
(a) A ∈ R 4 × 1 \mathbf{A}\in\real^{4 \times 1} A ∈ R 4 × 1 is the result of adding the first column of B \mathbf{B} B to -2 times the last column of B \mathbf{B} B .
(b) The rows of A ∈ R 4 × 3 \mathbf{A}\in\real^{4 \times 3} A ∈ R 4 × 3 are the rows of B \mathbf{B} B in order 4,3,2,1.
(c) The first column of A ∈ R 4 × 3 \mathbf{A}\in\real^{4 \times 3} A ∈ R 4 × 3 is 1 times the first column of B \mathbf{B} B , the second column of A \mathbf{A} A is 2 times the second column of B \mathbf{B} B ,
and the third column of A \mathbf{A} A is 3 times the third column of B \mathbf{B} B .
(d) A A A is the scalar sum of all elements of B \mathbf{B} B .
(a) ✍ Prove that for real vectors v \mathbf{v} v and w \mathbf{w} w of the same length, the inner products v T w \mathbf{v}^T\mathbf{w} v T w and w T v \mathbf{w}^T\mathbf{v} w T v are equal.
(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 and w v T \mathbf{w}\mathbf{v}^T w v T .