9.1. Matrix Operations#
9.1.1. Creating Vectors and Matrices#
We have already seen how to create 1D and 2D arrays in Julia. These can be naturally used to represent vectors and matrices:
x = randn(3) # Random vector
3-element Vector{Float64}:
1.170089928604851
-1.7581890050001607
0.17850340495995906
A = randn(3,3) # Random matrix
3×3 Matrix{Float64}:
0.100159 0.924001 0.437762
-0.0532033 -1.5979 1.63883
-0.961968 -0.957267 -0.321329
9.1.2. Basic Matrix Operations#
Julia also defines common matrix operations for arrays, some in the base library and some in the LinearAlgebra package.
using LinearAlgebra
For example addition, subtraction, and scalar multiplication work as expected:
y = randn(3)
x - 3y
3-element Vector{Float64}:
-2.715017441848958
-0.617033622804521
0.7797533050482875
Note that this happens to be the same operation as the corresponding elementwise operation x .- 3y, because of the definition of vector addition. However, if you try to multiply two vectors
z = x * y # Error - cannot multiply two vectors
you get an error. Use can use the dot function or the \(\cdot\) syntax (type \cdot and tab) to compute dot products:
x_dot_y_1 = dot(x,y) # Dot product
x_dot_y_2 = x ⋅ y # Same thing
2.1483222324306683
For the special case of vectors of length 3, Julia also defines the cross product:
x_cross_y = cross(x,y)
3-element Vector{Float64}:
0.42027036163434983
0.4656737823096879
1.8318328807655448
The 2-norm of a vector can be calculated using the norm function:
a = norm(x) # 2-norm of x
b = sqrt(sum(abs.(x).^2)) # Should be the same
a - b # Confirm
0.0
More generally, the norm function takes a second argument p in which case it computes the \(p\)-norm:
This includes the so-called max-norm, by setting p to Inf.
9.1.3. Matrix Multiplication#
Matrix multiplication is performed using the * operator. Recall the definition of the product of two matrices \(C=AB\) where \(A\) is \(m\)-by-\(k\), \(B\) is \(k\)-by-\(n\), and \(C\) is \(m\)-by-\(n\):
If the “middle dimensions” (\(k\) in the example above) do not match, Julia gives an error.
B = randn(4,3)
C = B*A # OK, since B is 4-by-3 and A is 3-by-3
4×3 Matrix{Float64}:
-0.509578 -0.793306 -0.26053
0.0224918 -0.194257 -0.155398
-0.210142 -2.39151 2.20166
0.115612 -0.829318 0.202901
Note that unlike addition and subtraction, matrix multiplication is completely different from elementwise multiplication:
AA = A*A # Square of matrix
A2 = A.*A # Square of each entry in matrix
A2 - AA # These are not the same
3×3 Matrix{Float64}:
0.470273 2.65674 -1.22583
1.49965 1.61796 5.85436
0.661693 -0.0319904 1.98992
Similarly, the power ^ operator will compute matrix powers, not elementwise powers:
A_to_the_power_of_2 = A^2 # Same as A*A
A_to_the_power_of_3_5 = A^3.5 # A^3.5, much harder to compute!
3×3 Matrix{ComplexF64}:
1.22298-0.976917im 1.31013+0.98218im -1.51813-1.36301im
1.43161+1.23172im 0.620983-1.23835im 2.29619+1.71851im
0.155062+0.275974im -0.49154-0.277461im 2.74272+0.385044im
Note that Julia automatically returns a matrix of complex numbers when needed (even if the input matrix was real).
9.1.4. Matrix Transpose#
Matrices can be transposed using the transpose function, or conjugate transposed using the adjoint function or the convenient ' syntax:
BT = transpose(B) # B is 4-by-3, so BT is 3-by-4
BT2 = B' # Same thing (since B is real so conjugate does not matter)
A * B' # Well-defined, since A is 3-by-3 and B' is 3-by-4
3×4 Matrix{Float64}:
0.205614 -0.076061 1.33351 0.08526
0.796908 -0.0347191 -1.9857 -0.713121
0.117239 0.334943 -1.31801 0.449975
Since the dot product between two vectors can be written as \(z = x^*y\) (conjugate transpose of \(x\)), matrix multiplication can be used to provide an alternative syntax:
x_dot_y_3 = x' * y
2.1483222324306683
In all these examples, the vectors x and y have been 1D arrays. It is also possible to use 2D arrays (or matrices) to represent vectors, which allows for both column- and row-vectors:
a = randn(3,1) # Column vector
b = randn(1,3) # Row vector
1×3 Matrix{Float64}:
0.477496 1.19731 -0.189257
Note that a' * b is now invalid because of the dimensions. One option is to do a' * b', but it is generally safer to use the dot function to compute the dot product:
a_dot_b = dot(a,b)
1.7676934685493575
9.1.5. Other Matrix Operations#
The LinearAlgebra package defines many other common matrix operations. For example, determinant, trace, and inverse matrix (but note that inverse matrices are often inefficient and ill-conditioned, use \ instead for solving linear systems):
A = randn(3,3)
println("det(A) = ", det(A))
println("tr(A) = ", tr(A))
println("inv(A) = ", inv(A))
det(A) = -5.7051110913542
tr(A) = 1.3910049113959708
inv(A) = [0.012099035501579595 -0.13029261793896074 -0.44297849809270884; 0.5860275009391351 1.0481059549952612 0.26242732525290596; -0.35898961734377355 0.026100367589614572 -0.1574972022909212]
9.1.6. Eigenvalues and Eigenvectors#
Eigenvalues and eigenvectors can be computed using the functions eigvals and eigvects. Note that the output of these functions is in general complex, even though the matrices are real.
A = randn(3,3)
lambda = eigvals(A) # Eigenvalues
3-element Vector{ComplexF64}:
0.6095685531281074 - 0.402127648755167im
0.6095685531281074 + 0.402127648755167im
1.5021885497949794 + 0.0im
X = eigvecs(A) # Eigenvectors
3×3 Matrix{ComplexF64}:
-0.459935+0.379712im -0.459935-0.379712im -0.980182+0.0im
0.174134-0.373901im 0.174134+0.373901im 0.19777+0.0im
0.688588-0.0im 0.688588+0.0im -0.011381+0.0im
norm(A*X - X*Diagonal(lambda)) # Confirm A*X = X*LAMBDA
2.3050943929982774e-15
The Diagonal function above is an example of a special matrix in Julia discussed in the next section.
9.1.7. Example: Eigenvalues of Random Matrices#
Random matrix theory studies the distribution of the eigenvalues of certain classes of randomly generated matrices, which have many important applications. For example, matrices with entries from the normal distribution can be shown to have eigenvalues concentrated inside a circle in the complex plane:
using PyPlot
n = 1000
A = randn(n,n)
l = eigvals(A)
plot(real(l), imag(l), ".")
axis("equal");
# Draw circle
radius = sqrt(n)
phi = 2π*(0:100)/100
plot(radius*cos.(phi), radius*sin.(phi));
9.1.8. Other data types#
We already saw how matrix operations can be performed on complex numbers. Many operations are also defined for rational numbers, which allows for exact linear algebra computations:
num = rand(-10:10, 3, 3)
den = rand(1:10, 3, 3)
rat = num .// den
display(rat)
inv(rat)
3×3 Matrix{Rational{Int64}}:
4//7 0//1 -3//1
1//1 6//5 7//1
-7//10 -9//2 6//5
3×3 Matrix{Rational{Int64}}:
21//19 525//1159 140//1159
-35//171 -55//1159 -2450//10431
-7//57 100//1159 80//3477