Alex Ghitza

Linear algebra in Julia: Matrix operations

A. Ghitza

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 `Float64` before doing the work; as a pure mathematician, it hurts my soul that the conversion is not done to `Rational` by 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.

@ 2023 Alexandru Ghitza · Created with Zola and PureCSS