Linear algebra in Julia: Matrix operations
2021-01-31
Packages
We mentioned the
LinearAlgebra package in the previous
post. Since it is included with the default installation of Julia, we just
needed to activate it with
julia> using LinearAlgebra
There are thousands of other packages adding all sorts of extra functionality.
Julia has a builtin package manager
Pkg that can take care of all the dirty
business of downloading and building packages with their dependencies.
As an example, let's get a package that computes the reduced row echelon form
of a matrix. First we must enter the
Pkg REPL by pressing
] in Julia, whereupon the prompt changes
to something like
(@v1.5) pkg>
Time to get the package:
(@v1.5) pkg> add RowEchelon
(chatter omitted)
We can return to Julia by pressing
Backspace at the
pkg> prompt, then take the new package
for a spin.
Reduced row echelon form
julia> using RowEchelon
[ Info: Precompiling RowEchelon (code omitted)]
julia> A = [[-1, 1, 3] [2, 0, 1] [1, 1, 0] [-3, -1, 3] [1, 3, -1]];
julia> rref(A)
3x5 Array{Float64,2}:
1.0 0.0 0.0 1.0 5.5112e-17
0.0 1.0 0.0 0.0 -1.0
0.0 0.0 1.0 -2.0 3.0
Some discussion is in order. The row reduction algorithm is numerically
unstable (small errors accumulate), which is exemplified by the appearance of
5.5112e-17 in the upper right corner, as a
synonym for 0.0. We can choose to stay in
the realm of exact (as opposed to approximate) computation by converting to
Rational:
julia> B = Array{Rational}(A)
3×5 Array{Rational,2}:
-1//1 2//1 1//1 -3//1 1//1
1//1 0//1 1//1 -1//1 3//1
3//1 1//1 0//1 3//1 -1//1
julia> rref(B)
3×5 Array{Rational,2}:
1//1 0//1 0//1 1//1 0//1
0//1 1//1 0//1 0//1 -1//1
0//1 0//1 1//1 -2//1 3//1
Yes, rational numbers look rather ugly in Julia.
Suppose we wish to solve the system of equations TODO. We can start by computing the rref of the augmented matrix.
julia> C = [[4, -1, 12, -7] [4, -5, 0, 0] [-4, 9, 12, -7] [-2, -2, 13, 11] [5, -5, 9, 14] [-5, 14, -1, -14]]
4×6 Array{Int64,2}:
4 4 -4 -2 5 -5
-1 -5 9 -2 -5 14
12 0 12 13 9 -1
-7 0 -7 11 14 -14
julia> rref(C)
4×6 Array{Float64,2}:
1.0 0.0 1.0 0.0 0.0 0.824724
0.0 1.0 -2.0 0.0 0.0 -2.74216
0.0 0.0 0.0 1.0 0.0 -0.945917
0.0 0.0 0.0 0.0 1.0 0.155582
We can read the (parametrized) solution off the rref.
Another way of solving the system is by treating the homogeneous and non-homogeneous parts separately:
julia> D = C[1:4, 1:5]
4×5 Array{Int64,2}:
4 4 -4 -2 5
-1 -5 9 -2 -5
12 0 12 13 9
-7 0 -7 11 14
julia> z = nullspace(D)
5×1 Array{Float64,2}:
0.40824829046386274
-0.8164965809277263
-0.4082482904638629
-6.095501145100942e-16
4.691634128813368e-16
This vector z forms a basis for the
solution space of the homogeneous system.
It remains to get one solution of the non-homogeneous system:
julia> v = C[1:4, 6]
4-element Array{Int64,1}:
-5
14
-1
-14
julia> x = D\v
5-element Array{Float64,1}:
-0.2267838206332375
-0.639145584002424
1.0515073473716108
-0.9459172852598086
0.15558248750189305
Every solution of the homogeneous system is then a sum of this vector
x and a scalar multiple of the vector
z.
Matrix inverse and determinant
Of course, you need a square matrix for such things:
julia> size(D)
(4, 5)
julia> det(D)
ERROR: DimensionMismatch("matrix is not square: dimensions are (4, 5)")
Stacktrace: (omitted)
Let's try again with a 4x4 matrix
julia> E = D[1:4, 1:4]
4×4 Array{Int64,2}:
4 4 -4 -2
-1 -5 9 -2
12 0 12 13
-7 0 -7 11
julia> det(E)
-3.765876499528531e-13
Also known as zero, but worth noticing a couple of things:
-
Julia likes to do linear algebra in floating point, and if you give it an
integer matrix it will convert it to
Float64before doing the work; as a pure mathematician, it hurts my soul that the conversion is not done toRationalby default, but presumably the performance gain in working with floats is significant enough to justify this choice, given Julia's target audience (not pure mathematicians); - error propagation is a thing, and the sooner you internalize this the less painful your journey into numerical computation will be.
Moving on to inverses:
julia> inv(E)
4×4 Array{Float64,2}:
2.9608e15 2.36864e15 -1.32771e14 1.1259e15
-5.9216e15 -4.73728e15 2.65542e14 -2.2518e15
-2.9608e15 -2.36864e15 1.32771e14 -1.1259e15
-0.141509 -0.113208 0.0377358 0.0
julia> E * inv(E)
4×4 Array{Float64,2}:
0.283019 0.226415 0.424528 0.0
0.283019 0.226415 -0.325472 0.0
-1.83962 -1.4717 1.99057 0.0
-1.5566 -1.24528 -0.459906 0.0
In other words, if you're not careful you can get complete garbage.
The "truth" is of course
julia> F = Array{Rational}(E)
4×4 Array{Rational,2}:
4//1 4//1 -4//1 -2//1
-1//1 -5//1 9//1 -2//1
12//1 0//1 12//1 13//1
-7//1 0//1 -7//1 11//1
julia> det(F)
0//1
julia> inv(F)
ERROR: SingularException(3)
Stacktrace: (omitted)
We'll have another go at this with a different matrix (and another method):
julia> A = [[0, 0, 1, 2] [0, -1, 0, 1] [1, 3, 0, 0] [0, 0, 1, -3]]
4×4 Array{Int64,2}:
0 0 1 0
0 -1 3 0
1 0 0 1
2 1 0 -3
julia> det(A)
-5.0
Looks promising.
julia> inv(A)
4×4 Array{Float64,2}:
-0.6 0.2 0.6 0.2
3.0 -1.0 0.0 0.0
1.0 0.0 0.0 0.0
0.6 -0.2 0.4 -0.2
Now the other approach:
julia> B = Matrix(I, 4, 4)
4×4 Array{Bool,2}:
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
julia> M = [A B]
4×8 Array{Int64,2}:
0 0 1 0 1 0 0 0
0 -1 3 0 0 1 0 0
1 0 0 1 0 0 1 0
2 1 0 -3 0 0 0 1
julia> rref(M)
4×8 Array{Float64,2}:
1.0 0.0 0.0 0.0 -0.6 0.2 0.6 0.2
0.0 1.0 0.0 0.0 3.0 -1.0 -2.22045e-16 1.11022e-16
0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 0.6 -0.2 0.4 -0.2
This consists of two 4x4 matrices stacked horizontally; the fact that the first one is the identity matrix tells us that the second matrix is the inverse of A.
Just for fun, we can try to do the same with the matrix E that bombed spectacularly above:
julia> M = [E B]
4×8 Array{Int64,2}:
4 4 -4 -2 1 0 0 0
-1 -5 9 -2 0 1 0 0
12 0 12 13 0 0 1 0
-7 0 -7 11 0 0 0 1
julia> rref(M)
4×8 Array{Float64,2}:
1.0 0.0 1.0 0.0 0.0 0.0 0.0493274 -0.058296
0.0 1.0 -2.0 0.0 0.0 -0.2 -0.0224215 -0.00986547
0.0 0.0 0.0 1.0 0.0 0.0 0.0313901 0.0538117
0.0 0.0 0.0 0.0 1.0 0.8 -0.044843 0.380269
The first half of the rref is pretty convincingly not the 4x4 identity matrix, so E is not invertible.