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:

\[ \|A\|_p = \left( \sum_{i=1}^n | a_i | ^p \right)^{1/p} \]

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\):

\[ C_{ij} = \sum_{\ell=1}^k A_{i\ell} B_{\ell j},\qquad\text{ for }i=1,\ldots,m\text{ and }j=1,\ldots,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));
../../_images/Matrix_Operations_39_0.png

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