Chapter 40 Introduction to Linear Algebra
Linear Algebra is probably the most important mathematical field for Scientific Computation. This chapter covers the basics of Linear Algebra in Julia and how it can be used to solve scientific problems.
1
And if you ask someone in the field, they will say most important field PERIOD!
Also, a lot of these functions described here are standard functions (that is, part of the
Base package which is automatically loaded), however many are also part of the LinearAlgebra package, which should be downloaded, but you will need to perform using LinearAlgebra. It is pointed out when you need this package.
Section 40.1 Vectors and Matrices
Vector and matrices are the building blocks of linear algebra. In short a Vector is simply a 1D array and in julia the type
Vector<T> is shorthand for Array{T,1}. We can make a vector just like we did with an array. For example, x = [1,2,3] makes an array of length three and the output:
3-element Vector{Int64}:
1
2
3
shows that it is a 1D integer array of length 3.
Similarly, there is a
Matrix type which is an alias for a 2D array. If
A = [i+j for i=1:3,j=1:3]
then we get the results
3×3 Matrix{Int64}:
2 3 4
3 4 5
4 5 6
Subsection 40.1.1 Addition and Subtraction of Vectors and Matrices
All of the construction methods for arrays of both 1D and 2D as seen in Chapter \ref{ch:arrays} can be used. If we have
3×3 Matrix{Int64}:
5 5 5
3 6 5
5 7 9
and
A-B returns
3×3 Matrix{Int64}:
-1 1 3
3 2 5
3 3 3
5-element Vector{Int64}:
2
5
8
11
14
and
x-y returns
5-element Vector{Int64}:
0
-1
-2
-3
-4
Subsection 40.1.2 scalar multiplication
Scalar multiplication of vectors and matrix work as expected using the
* operator, which is implied if multiplying on the left by a number. Using A and x above, 3x returns
5-element Vector{Int64}:
3
6
9
12
15
and
-4A returns
3×3 Matrix{Int64}:
-8 -12 -16
-12 -16 -20
-16 -20 -24
Subsection 40.1.3 Matrix-vector products
3-element Vector{Int64}:
34
44
54
and the matrix-matrix products
A*B returns
3×3 Matrix{Int64}:
10 18 14
14 24 18
18 30 22
Note that if either scalar or matrix multiplication involves variables, then you will need to explicitly put in a
*, since Julia can have variables of any length, if for example, you entered AB expecting the product, you will get UndefVarError: `AB` not defined in `Main`.
Subsection 40.1.4 Matrix Transpose
Recall that the tranpose of a matrix \(A\) is a matrix where the \(i\)th row and \(j\) column of A becomes the \(i\)th column and \(j\)th row of the transpose. In julia, the
transpose function performs like. For example,
transpose(B)
returns
3×3 transpose(::Matrix{Int64}) with eltype Int64:
3 0 1
2 2 2
1 0 3
Notice that the type of a transpose is not a simple
Matrix. It has transpose in the type. This is mainly for efficiency. For large matrices (say ones that take up a Gb of memory or more), simply a copy of a matrix is more efficient than doing the transpose operation. Julia will then known if it is a transpose then perform the correct operation using the transpose when needed.
Subsection 40.1.5 Trace and Determinant of a Matrix
The trace of a matrix \(A\) is the sum of the elements along the diagonal. In julia, you can use the
tr function in the LinearAlgebra package. The functiontr(A)
returns
12
Also, the determinant of a matrix can be easily found in julia with the
det function in the LinearAlgebra package.
det(A)
returns
0.0
Subsection 40.1.6 Identity Matrix
Recall that the identity matrix, \(I\text{,}\) is the matrix that satisfies \(AI=A\) for any matrix \(A\text{.}\) It is matrix with ones along the diagonal and zeros elsewhere. If needed in julia, you can use
Matrix(I,3,3)
which returns
3×3 Matrix{Bool}:
1 0 0
0 1 0
0 0 1
Note: the
I is part of the LinearAlgebra package and can be used without the size often. For example, I*B returns
3×3 Matrix{Int64}:
3 2 1
0 2 0
1 2 3
Section 40.2 Matrix Inverse
If you have a matrix \(A\text{,}\) then if it exists, the inverse of \(A\text{,}\) denoted \(A^{-1}\) satisfies \(AA^{-1}=I\) and \(A^{-1}A=I\text{.}\)
3×3 Matrix{Float64}:
0.375 -0.25 -0.125
0.0 0.5 -0.0
-0.125 -0.25 0.375
If we try to find the inverse of \(A\text{,}\) above,
inv(A) we get SingularException(3) which occurs because the matrix \(A\) is singular, that is, it doesn’t have a matrix inverse.
Section 40.3 Solving Linear Systems
One of the most important aspects of linear algebra is that of solving a linear system. In general, this is solving the linear system:
\begin{equation*}
A \mathbf{x} = \mathbf{b}
\end{equation*}
A = [0 1 1; 1 2 1; 1 1 -1] b = [1,2,3]
then the solution can be found with
x = A \ b
which returns
3-element Vector{Float64}:
-2.0
3.0
-2.0
which is the same as \(A^{-1} \mathbf{b}\text{,}\) but is more efficient.
Section 40.4 Eigenvalues and Eigenvectors
Recall that a scalar \(\lambda\) and corresponding vector \(\mathbf{v}\) form an eigenvalue-eigenvector pair if
\begin{equation*}
A \mathbf{v} = \lambda \mathbf{v}
\end{equation*}
To find the eigenvalues of a matrix \(A\) in julia,
eigvals(A)
which returns
3-element Vector{Float64}:
-1.6554423815498303
-0.21075588095919165
2.8661982625090223
The eigenvectors can be found with
eigvecs(A)
which returns
3×3 Matrix{Float64}:
-0.462272 0.785797 -0.410888
-0.114102 -0.51223 -0.851235
0.879367 0.346619 -0.326451
where the eigenvectors are the columns of the resulting matrix. The eigenvectors are in the same order as the eigenvalues, which is important in that they form pairs.
To get a particular eigenvector, you can use the standard array notation to pull out a particular column. For example,
eigvecs(A)[:,2] will return the 2nd eigenvector.
Section 40.5 Special Matrices
Special matrices arise throughout Linear Algebra.
Subsection 40.5.1 Diagonal Matrices
For example, the identity matrix, can be thought of as a diagonal matrix and in julia there is a matrix type called
Diagonal. For example,
D=Diagonal([1,2,3,4])
returns
4×4 Diagonal{Int64, Vector{Int64}}:
1 ⋅ ⋅ ⋅
⋅ 2 ⋅ ⋅
⋅ ⋅ 3 ⋅
⋅ ⋅ ⋅ 4
and you’ll notice that it doesn’t just create a matrix where there are zeros on the off-diagonal. It will print out like a regular matrix, but it will only store the diagonal elements. This is precisely for saving on storage. However, you can still do operations with this matrix. For example, finding the eigenvalues of this with
eigvals(D)
and this returns
4-element Vector{Int64}:
1
2
3
4
which is what we expect as the eigenvalues of a diagonal matrix are the diagonal elements.
Subsection 40.5.2 Tridiagonal Matrices
Another matrix that shows up often are tridiagonal matrices. For example, the following arises from discretizes a differential equation:
T = Tridiagonal([1 for i=1:9],[-2 for i=1:10],[1 for i=1:9])
10×10 Tridiagonal{Int64, Vector{Int64}}:
-2 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
1 -2 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1 -2 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1 -2 1 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1 -2 1 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 1 -2 1 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1 -2 1 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 -2 1 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 -2 1
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 -2
Section 40.6 Abstract Matrices
As with number types that we discussed in Section 21.3, matrices in julia are subtypes of the abstract type
AbstractMatrix. There are a few functions that are common to all matrices and understanding this is helpful in working with different types of special matrices.
In short, an
AbstractMatrix is the supertype of all matrices in julia. For example the multiplication of two matrices with *. If
T = Tridiagonal([-1 for i=1:4],[2 for j=1:5],[-1 for j=1:4]) D = Diagonal([1,2,3,4,5])
then
D*T returns
5×5 Tridiagonal{Int64, Vector{Int64}}:
2 -2 ⋅ ⋅ ⋅
-1 4 -3 ⋅ ⋅
⋅ -2 6 -4 ⋅
⋅ ⋅ -3 8 -5
⋅ ⋅ ⋅ -4 10
and notice that the type is a
Tridiagonal. Depending on the types of the two input matrices, the multiplication will result in the appropriate output type.
Subsection 40.6.1 Solving with Matrices/Division
Another common
