Matrix Operations
Contents
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.0695614007143797
-0.7354166947850529
0.8053390643232473
A = randn(3,3) # Random matrix
3×3 Matrix{Float64}:
-0.846352 -1.08618 -0.586363
-0.948634 -1.5764 1.11086
0.945739 1.06962 0.391572
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.481854082620411
0.15475338642151515
0.960525077075872
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
1.4427084998285056
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.27700504162095096
1.0086915305143034
0.5532262258777343
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.38307 -0.214799 -0.951445
1.03021 1.41365 -0.236946
-1.02927 -1.47234 -0.612455
1.43243 2.74791 -3.03221
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.475837 -0.824565 1.28374
-2.44897 -2.21858 1.99392
2.3392 3.43864 -0.633649
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}:
-0.0541881-5.93359im -0.0107824-8.71349im -0.070167+2.09157im
0.0414486-12.9819im 0.00744095-19.064im 0.052079+4.57609im
0.0189485+6.57725im 0.000410366+9.65871im 0.0179041-2.31846im
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.750819 0.306796 -2.06616 2.47226
-0.212385 1.41742 -1.19989 2.78556
-0.584636 -0.411776 2.1116 -2.30109
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
1.4427084998285056
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.695701 0.107967 -1.91251
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)
-0.5327813132526754
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) = -0.2002573534834639
tr(A) = 2.6936818694474303
inv(A) = [-7.036833017953787 -0.616015709095034 -2.159320083999361; -4.374825769962134 0.25677908032874996 -1.4123386690116264; -10.292797047476093 -0.5968291949694127 -2.082452509471869]
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}:
-1.6574103438863186 + 0.0im
0.9594431310210694 - 0.4460157310613137im
0.9594431310210694 + 0.4460157310613137im
X = eigvecs(A) # Eigenvectors
3×3 Matrix{ComplexF64}:
0.613502+0.0im 0.741082-0.0im 0.741082+0.0im
0.379175+0.0im -0.500486-0.339601im -0.500486+0.339601im
0.692706+0.0im -0.192428+0.218984im -0.192428-0.218984im
norm(A*X - X*Diagonal(lambda)) # Confirm A*X = X*LAMBDA
2.1689344116778317e-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}}:
-5//3 -9//1 -3//4
1//3 -2//9 9//1
9//5 -3//7 -5//3
3×3 Matrix{Rational{Int64}}:
-47940//1792159 166455//1792159 920430//1792159
-190008//1792159 -46809//1792159 -167265//1792159
-2916//1792159 191808//1792159 -38220//1792159