A lot of the Data Science methods we will see in this tutorial require some understanding of linear algebra, and in this notebook we will focus on how Julia handles matrices, the types that exist, and how to call basic linear algebra tasks.
# some packages we will use
using LinearAlgebra
using SparseArrays
using Images
using MAT
┌ Info: Precompiling MAT [23992714-dd62-5051-b70f-ba57cb901cac] └ @ Base loading.jl:1317
We will get started with creating a random matrix.
A = rand(10,10); # created a random matrix of size 10-by-10
Atranspose = A' # matrix transpose
A = A*Atranspose; # matrix multiplication
@show A[11] == A[1,2];
A[11] == A[1, 2] = true
b = rand(10); #created a random vector of size 10
x = A\b; #x is the solutions to the linear system Ax=b
@show norm(A*x-b)
;
norm(A * x - b) = 3.689263419778939e-12
A few things that are noteworthy:
A
is a Matrix
type, and b
is a Vector
type.Adjoint
.\
is always the recommended way to solve a linear system. You almost never want to call the inv
function@show typeof(A)
@show typeof(b)
@show typeof(rand(1,10))
@show typeof(Atranspose)
;
typeof(A) = Matrix{Float64} typeof(b) = Vector{Float64} typeof(rand(1, 10)) = Matrix{Float64} typeof(Atranspose) = Adjoint{Float64, Matrix{Float64}}
Matrix{Float64} == Array{Float64,2}
true
Vector{Float64} == Array{Float64,1}
true
Atranspose
10×10 adjoint(::Matrix{Float64}) with eltype Float64: 0.249899 0.332817 0.376245 … 0.196724 0.533504 0.545816 0.51139 0.291047 0.0906729 0.0955851 0.145757 0.328184 0.709529 0.27097 0.0946509 0.661237 0.785675 0.723368 0.93995 0.491115 0.303725 0.705366 0.52196 0.999402 0.10152 0.819124 0.997739 0.616617 0.985117 0.588296 0.556655 0.0694413 0.229682 … 0.402893 0.0629227 0.145643 0.701979 0.743977 0.275638 0.258079 0.389031 0.737054 0.903376 0.746583 0.096613 0.411208 0.0589579 0.463674 0.409278 0.890039 0.789442 0.361242 0.456934 0.839767 0.381081 0.110646 0.654256 0.0364732 0.277699 0.366428
adjoint
in julia is a lazy adjoint -- often, we can easily perform Linear Algebra operations such as A*A'
without actually transposing the matrix.
?adjoint
search: adjoint adjoint! Adjoint
A'
adjoint(A)
Lazy adjoint (conjugate transposition). Note that adjoint
is applied recursively to elements.
For number types, adjoint
returns the complex conjugate, and therefore it is equivalent to the identity function for real numbers.
This operation is intended for linear algebra usage - for general data manipulation see [permutedims
](@ref Base.permutedims).
jldoctest
julia> A = [3+2im 9+2im; 8+7im 4+6im]
2×2 Matrix{Complex{Int64}}:
3+2im 9+2im
8+7im 4+6im
julia> adjoint(A)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
3-2im 8-7im
9-2im 4-6im
julia> x = [3, 4im]
2-element Vector{Complex{Int64}}:
3 + 0im
0 + 4im
julia> x'x
25 + 0im
Atranspose.parent
10×10 Matrix{Float64}: 0.249899 0.51139 0.709529 … 0.903376 0.409278 0.381081 0.332817 0.291047 0.27097 0.746583 0.890039 0.110646 0.376245 0.0906729 0.0946509 0.096613 0.789442 0.654256 0.866039 0.00363302 0.0391991 0.283632 0.524616 0.182196 0.54766 0.403914 0.0984212 0.805436 0.219447 0.384202 0.798657 0.602425 0.630361 … 0.0758388 0.791843 0.0837277 0.0147225 0.0898888 0.39432 0.8799 0.369629 0.30327 0.196724 0.0955851 0.661237 0.411208 0.361242 0.0364732 0.533504 0.145757 0.785675 0.0589579 0.456934 0.277699 0.545816 0.328184 0.723368 0.463674 0.839767 0.366428
sizeof(A)
800
That's because it's an array of Float64's, each is of size 8 bytes, and there are 10*10 numbers.
# To actually copy the matrix:
B = copy(Atranspose)
10×10 Matrix{Float64}: 0.249899 0.332817 0.376245 … 0.196724 0.533504 0.545816 0.51139 0.291047 0.0906729 0.0955851 0.145757 0.328184 0.709529 0.27097 0.0946509 0.661237 0.785675 0.723368 0.93995 0.491115 0.303725 0.705366 0.52196 0.999402 0.10152 0.819124 0.997739 0.616617 0.985117 0.588296 0.556655 0.0694413 0.229682 … 0.402893 0.0629227 0.145643 0.701979 0.743977 0.275638 0.258079 0.389031 0.737054 0.903376 0.746583 0.096613 0.411208 0.0589579 0.463674 0.409278 0.890039 0.789442 0.361242 0.456934 0.839767 0.381081 0.110646 0.654256 0.0364732 0.277699 0.366428
sizeof(B)
800
The \
operator allows you to solve a system of linear equations, and often uses a suitable matrix factorization to solve the problem. We will cover factorizations next.
?\
search: \
\(x, y)
Left division operator: multiplication of y
by the inverse of x
on the left. Gives floating-point results for integer arguments.
jldoctest
julia> 3 \ 6
2.0
julia> inv(3) * 6
2.0
julia> A = [4 3; 2 1]; x = [5, 6];
julia> A \ x
2-element Vector{Float64}:
6.5
-7.0
julia> inv(A) * x
2-element Vector{Float64}:
6.5
-7.0
\(A, B)
Matrix division using a polyalgorithm. For input matrices A
and B
, the result X
is such that A*X == B
when A
is square. The solver that is used depends upon the structure of A
. If A
is upper or lower triangular (or diagonal), no factorization of A
is required and the system is solved with either forward or backward substitution. For non-triangular square matrices, an LU factorization is used.
For rectangular A
the result is the minimum-norm least squares solution computed by a pivoted QR factorization of A
and a rank estimate of A
based on the R factor.
When A
is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt
factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices.
jldoctest
julia> A = [1 0; 1 -2]; B = [32; -4];
julia> X = A \ B
2-element Vector{Float64}:
32.0
18.0
julia> A * X == B
true
(\)(F::QRSparse, B::StridedVecOrMat)
Solve the least squares problem $\min\|Ax - b\|^2$ or the linear system of equations $Ax=b$ when F
is the sparse QR factorization of $A$. A basic solution is returned when the problem is underdetermined.
jldoctest
julia> A = sparse([1,2,4], [1,1,1], [1.0,1.0,1.0], 4, 2)
4×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
1.0 ⋅
1.0 ⋅
⋅ ⋅
1.0 ⋅
julia> qr(A)\fill(1.0, 4)
2-element Vector{Float64}:
1.0
0.0
A common tool used in Linear Algebra is matrix factorizations. These factorizations are often used to solve linear systems like Ax=b
, and as we will see later in this tutorial... Ax=b
comes up in a lot of Data Science problems
L*U = P*A
luA = lu(A)
LU{Float64, Matrix{Float64}} L factor: 10×10 Matrix{Float64}: 1.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.714777 1.0 0.0 0.0 0.0 0.0 0.0 0.431294 0.879205 1.0 0.0 0.0 0.0 0.0 0.402595 0.494058 0.42929 1.0 0.0 0.0 0.0 0.55784 0.0226054 0.41306 0.372797 0.0 0.0 0.0 0.63672 0.756747 0.509309 0.628019 … 0.0 0.0 0.0 0.610917 0.308904 0.236169 0.0853205 0.0 0.0 0.0 0.550327 0.675756 0.622177 0.13992 1.0 0.0 0.0 0.908177 0.621365 0.340249 0.274056 0.267463 1.0 0.0 0.717008 0.291136 0.291335 -0.268112 -0.146209 -0.292422 1.0 U factor: 10×10 Matrix{Float64}: 3.65267 2.61085 1.57538 1.47055 … 2.23148 2.01016 3.31727 0.0 1.23501 1.08582 0.610166 0.381498 0.834564 0.76739 0.0 0.0 0.801694 0.344159 0.189336 0.498796 0.272776 0.0 0.0 0.0 0.572375 0.0488353 0.0800869 0.156863 0.0 0.0 0.0 0.0 -0.120383 -0.353458 -0.455433 0.0 0.0 0.0 0.0 … 0.185715 0.362473 0.131452 0.0 0.0 0.0 0.0 0.292697 0.176947 0.00272506 0.0 0.0 0.0 0.0 0.023562 0.185526 0.0335832 0.0 0.0 0.0 0.0 0.0 -0.0124478 0.010775 0.0 0.0 0.0 0.0 0.0 0.0 0.00230443
norm(luA.L*luA.U - luA.P*A)
6.280369834735101e-16
Q*R = A
qrA = qr(A)
LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}} Q factor: 10×10 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}: -0.466731 0.556721 0.116392 … -0.0400424 0.508724 -0.134302 -0.333609 -0.420934 0.685112 0.108506 0.0656284 -0.184771 -0.201298 -0.47984 -0.535905 -0.287502 0.108374 -0.211555 -0.187904 -0.180434 -0.134551 0.0339839 0.242066 0.0759608 -0.260361 0.29205 -0.382156 0.216634 -0.283205 -0.204056 -0.297177 -0.265198 -0.0384313 … -0.414573 0.0930191 0.259962 -0.33465 0.160772 -0.0662353 -0.154548 0.049564 0.628225 -0.285134 0.0871597 -0.0063624 -0.258897 -0.158312 -0.598196 -0.256855 -0.246975 -0.22312 0.767019 0.198405 0.0427175 -0.423874 -0.00321333 0.100231 0.0691344 -0.71471 0.183707 R factor: 10×10 Matrix{Float64}: -7.82607 -7.39288 -5.76651 -4.70661 … -6.45377 -8.68162 0.0 -1.50819 -1.85438 -1.14894 -1.7525 -1.48307 0.0 0.0 -0.718806 -0.473242 -0.457496 -0.145766 0.0 0.0 0.0 -0.631469 -0.157768 -0.266435 0.0 0.0 0.0 0.0 0.677052 0.616209 0.0 0.0 0.0 0.0 … -0.361714 -0.107839 0.0 0.0 0.0 0.0 -0.217956 -0.00306423 0.0 0.0 0.0 0.0 0.148502 0.0280145 0.0 0.0 0.0 0.0 0.00907696 -0.00774295 0.0 0.0 0.0 0.0 0.0 0.0014477
norm(qrA.Q*qrA.R - A)
1.1605946732881715e-14
L*L' = A
isposdef(A)
true
cholA = cholesky(A)
Cholesky{Float64, Matrix{Float64}} U factor: 10×10 UpperTriangular{Float64, Matrix{Float64}}: 1.9112 1.36608 0.824288 0.769438 … 1.16758 1.05178 1.7357 ⋅ 1.11131 0.977068 0.549051 0.343287 0.750974 0.690528 ⋅ ⋅ 0.895374 0.384375 0.21146 0.557081 0.30465 ⋅ ⋅ ⋅ 0.756554 0.0645496 0.105857 0.207339 ⋅ ⋅ ⋅ ⋅ -0.124073 -0.364293 -0.469395 ⋅ ⋅ ⋅ ⋅ … 0.233683 0.456096 0.165405 ⋅ ⋅ ⋅ ⋅ 0.537544 0.281694 -0.00656903 ⋅ ⋅ ⋅ ⋅ 0.061188 0.41714 0.102246 ⋅ ⋅ ⋅ ⋅ ⋅ 0.157784 -0.0366897 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0887722
norm(cholA.L*cholA.U - A)
1.4043333874306805e-15
cholA.L
10×10 LowerTriangular{Float64, Matrix{Float64}}: 1.9112 ⋅ ⋅ … ⋅ ⋅ ⋅ 1.36608 1.11131 ⋅ ⋅ ⋅ ⋅ 0.824288 0.977068 0.895374 ⋅ ⋅ ⋅ 0.769438 0.549051 0.384375 ⋅ ⋅ ⋅ 1.06614 0.0251216 0.369843 ⋅ ⋅ ⋅ 1.2169 0.840979 0.456022 … ⋅ ⋅ ⋅ 1.37034 0.323543 0.260853 ⋅ ⋅ ⋅ 1.16758 0.343287 0.21146 0.061188 ⋅ ⋅ 1.05178 0.750974 0.557081 0.41714 0.157784 ⋅ 1.7357 0.690528 0.30465 0.102246 -0.0366897 0.0887722
cholA.U
10×10 UpperTriangular{Float64, Matrix{Float64}}: 1.9112 1.36608 0.824288 0.769438 … 1.16758 1.05178 1.7357 ⋅ 1.11131 0.977068 0.549051 0.343287 0.750974 0.690528 ⋅ ⋅ 0.895374 0.384375 0.21146 0.557081 0.30465 ⋅ ⋅ ⋅ 0.756554 0.0645496 0.105857 0.207339 ⋅ ⋅ ⋅ ⋅ -0.124073 -0.364293 -0.469395 ⋅ ⋅ ⋅ ⋅ … 0.233683 0.456096 0.165405 ⋅ ⋅ ⋅ ⋅ 0.537544 0.281694 -0.00656903 ⋅ ⋅ ⋅ ⋅ 0.061188 0.41714 0.102246 ⋅ ⋅ ⋅ ⋅ ⋅ 0.157784 -0.0366897 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0887722
factorize(A)
Cholesky{Float64, Matrix{Float64}} U factor: 10×10 UpperTriangular{Float64, Matrix{Float64}}: 1.9112 1.36608 0.824288 0.769438 … 1.16758 1.05178 1.7357 ⋅ 1.11131 0.977068 0.549051 0.343287 0.750974 0.690528 ⋅ ⋅ 0.895374 0.384375 0.21146 0.557081 0.30465 ⋅ ⋅ ⋅ 0.756554 0.0645496 0.105857 0.207339 ⋅ ⋅ ⋅ ⋅ -0.124073 -0.364293 -0.469395 ⋅ ⋅ ⋅ ⋅ … 0.233683 0.456096 0.165405 ⋅ ⋅ ⋅ ⋅ 0.537544 0.281694 -0.00656903 ⋅ ⋅ ⋅ ⋅ 0.061188 0.41714 0.102246 ⋅ ⋅ ⋅ ⋅ ⋅ 0.157784 -0.0366897 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0887722
?factorize
search: factorize Factorization factorial
factorize(A)
Compute a convenient factorization of A
, based upon the type of the input matrix. factorize
checks A
to see if it is symmetric/triangular/etc. if A
is passed as a generic matrix. factorize
checks every element of A
to verify/rule out each property. It will short-circuit as soon as it can rule out symmetry/triangular structure. The return value can be reused for efficient solving of multiple systems. For example: A=factorize(A); x=A\b; y=A\C
.
Properties of A |
type of factorization |
---|---|
Positive-definite | Cholesky (see cholesky ) |
Dense Symmetric/Hermitian | Bunch-Kaufman (see bunchkaufman ) |
Sparse Symmetric/Hermitian | LDLt (see ldlt ) |
Triangular | Triangular |
Diagonal | Diagonal |
Bidiagonal | Bidiagonal |
Tridiagonal | LU (see lu ) |
Symmetric real tridiagonal | LDLt (see ldlt ) |
General square | LU (see lu ) |
General non-square | QR (see qr ) |
If factorize
is called on a Hermitian positive-definite matrix, for instance, then factorize
will return a Cholesky factorization.
jldoctest
julia> A = Array(Bidiagonal(fill(1.0, (5, 5)), :U))
5×5 Matrix{Float64}:
1.0 1.0 0.0 0.0 0.0
0.0 1.0 1.0 0.0 0.0
0.0 0.0 1.0 1.0 0.0
0.0 0.0 0.0 1.0 1.0
0.0 0.0 0.0 0.0 1.0
julia> factorize(A) # factorize will check to see that A is already factorized
5×5 Bidiagonal{Float64, Vector{Float64}}:
1.0 1.0 ⋅ ⋅ ⋅
⋅ 1.0 1.0 ⋅ ⋅
⋅ ⋅ 1.0 1.0 ⋅
⋅ ⋅ ⋅ 1.0 1.0
⋅ ⋅ ⋅ ⋅ 1.0
This returns a 5×5 Bidiagonal{Float64}
, which can now be passed to other linear algebra functions (e.g. eigensolvers) which will use specialized methods for Bidiagonal
types.
?diagm
search: diagm spdiagm diag diagind Diagonal isdiag Bidiagonal Tridiagonal
diagm(kv::Pair{<:Integer,<:AbstractVector}...)
diagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)
Construct a matrix from Pair
s of diagonals and vectors. Vector kv.second
will be placed on the kv.first
diagonal. By default the matrix is square and its size is inferred from kv
, but a non-square size m
×n
(padded with zeros as needed) can be specified by passing m,n
as the first arguments.
diagm
constructs a full matrix; if you want storage-efficient versions with fast arithmetic, see Diagonal
, Bidiagonal
Tridiagonal
and SymTridiagonal
.
jldoctest
julia> diagm(1 => [1,2,3])
4×4 Matrix{Int64}:
0 1 0 0
0 0 2 0
0 0 0 3
0 0 0 0
julia> diagm(1 => [1,2,3], -1 => [4,5])
4×4 Matrix{Int64}:
0 1 0 0
4 0 2 0
0 5 0 3
0 0 0 0
diagm(v::AbstractVector)
diagm(m::Integer, n::Integer, v::AbstractVector)
Construct a matrix with elements of the vector as diagonal elements. By default, the matrix is square and its size is given by length(v)
, but a non-square size m
×n
can be specified by passing m,n
as the first arguments.
jldoctest
julia> diagm([1,2,3])
3×3 Matrix{Int64}:
1 0 0
0 2 0
0 0 3
# convert(Diagonal{Int64,Array{Int64,1}},diagm([1,2,3]))
Diagonal([1,2,3])
3×3 Diagonal{Int64, Vector{Int64}}: 1 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 3
I
is a function
I(3)
3×3 Diagonal{Bool, Vector{Bool}}: 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ 1
Sparse matrices are stored in Compressed Sparse Column (CSC) form
using SparseArrays
S = sprand(5,5,2/5)
5×5 SparseMatrixCSC{Float64, Int64} with 12 stored entries: 0.27572 0.905187 0.102799 ⋅ ⋅ 0.789701 0.140221 0.763729 ⋅ ⋅ ⋅ ⋅ 0.209574 ⋅ 0.856283 0.0953587 ⋅ 0.577905 ⋅ ⋅ ⋅ ⋅ ⋅ 0.251486 0.238713
S.rowval
12-element Vector{Int64}: 1 2 4 1 2 1 2 3 4 5 3 5
Matrix(S)
5×5 Matrix{Float64}: 0.27572 0.905187 0.102799 0.0 0.0 0.789701 0.140221 0.763729 0.0 0.0 0.0 0.0 0.209574 0.0 0.856283 0.0953587 0.0 0.577905 0.0 0.0 0.0 0.0 0.0 0.251486 0.238713
S.colptr
6-element Vector{Int64}: 1 4 6 10 11 13
S.m
5
Let's get to the more "data science-y" side. We will do so by working with images (which can be viewed as matrices), and we will use the SVD
decomposition.
First let's load an image. I chose this image as it has a lot of details.
X1 = load("data/khiam-small.jpg")
@show typeof(X1)
X1[1,1] # this is pixel [1,1]
typeof(X1) = Matrix{RGB{N0f8}}
We can easily convert this image to gray scale.
Xgray = Gray.(X1)
We can easily extract the RGB layers from the image. We will make use of the reshape
function below to reshape a vector to a matrix.
R = map(i->X1[i].r,1:length(X1))
R = Float64.(reshape(R,size(X1)...))
G = map(i->X1[i].g,1:length(X1))
G = Float64.(reshape(G,size(X1)...))
B = map(i->X1[i].b,1:length(X1))
B = Float64.(reshape(B,size(X1)...))
;
Z = zeros(size(R)...) # just a matrix of all zeros of equal size as the image
RGB.(Z,G,Z)
We can easily obtain the Float64
values of the grayscale image.
Xgrayvalues = Float64.(Xgray)
283×283 Matrix{Float64}: 0.101961 0.0627451 0.0784314 0.0941176 … 0.509804 0.552941 0.666667 0.0666667 0.0980392 0.0745098 0.054902 0.505882 0.584314 0.501961 0.0784314 0.0862745 0.0784314 0.0862745 0.6 0.701961 0.611765 0.0862745 0.0666667 0.0745098 0.0941176 0.658824 0.705882 0.145098 0.0784314 0.101961 0.0901961 0.0745098 0.713725 0.682353 0.231373 0.0745098 0.0745098 0.0784314 0.0862745 … 0.729412 0.701961 0.168627 0.12549 0.0980392 0.0862745 0.0862745 0.627451 0.466667 0.192157 0.439216 0.447059 0.305882 0.137255 0.231373 0.184314 0.137255 0.458824 0.454902 0.45098 0.454902 0.196078 0.101961 0.117647 0.45098 0.466667 0.458824 0.45098 0.584314 0.121569 0.137255 0.458824 0.458824 0.458824 0.454902 … 0.521569 0.513725 0.12549 0.466667 0.45098 0.458824 0.47451 0.576471 0.741176 0.117647 0.45098 0.45098 0.462745 0.458824 0.560784 0.67451 0.117647 ⋮ ⋱ ⋮ 0.494118 0.47451 0.47451 0.462745 0.427451 0.435294 0.443137 0.47451 0.482353 0.470588 0.470588 0.439216 0.431373 0.431373 0.494118 0.501961 0.470588 0.45098 0.447059 0.447059 0.45098 0.470588 0.494118 0.490196 0.482353 0.431373 0.419608 0.419608 0.458824 0.482353 0.47451 0.466667 … 0.454902 0.435294 0.423529 0.443137 0.458824 0.45098 0.45098 0.309804 0.32549 0.34902 0.47451 0.478431 0.462745 0.462745 0.341176 0.345098 0.360784 0.482353 0.478431 0.458824 0.458824 0.419608 0.372549 0.321569 0.552941 0.552941 0.541176 0.533333 0.447059 0.411765 0.372549 0.552941 0.545098 0.576471 0.552941 … 0.435294 0.423529 0.407843 0.564706 0.552941 0.54902 0.505882 0.439216 0.431373 0.419608 0.568627 0.552941 0.517647 0.462745 0.439216 0.431373 0.427451
Next, we will downsample this image using the SVD. First, let's obtain the SVD decomposition.
SVD_V = svd(Xgrayvalues)
SVD{Float64, Float64, Matrix{Float64}} U factor: 283×283 Matrix{Float64}: -0.0635349 -0.131929 -0.116512 … 0.00868688 -0.0118621 -0.0626293 -0.150522 -0.112161 -0.00150511 0.0575172 -0.0616805 -0.165404 -0.105726 -0.0434962 0.0138229 -0.0618268 -0.1749 -0.0962379 0.00219005 0.0054939 -0.0626651 -0.178715 -0.0859759 -0.0395294 -0.0529901 -0.0624273 -0.18066 -0.0777654 … 0.00598019 0.0735455 -0.0617703 -0.1838 -0.0647771 -0.0149741 -0.015221 -0.0612071 -0.189903 -0.0440992 0.0615694 -0.0389833 -0.0600312 -0.193717 -0.026464 -0.0129492 -0.0310368 -0.0595516 -0.193191 -0.0082394 -0.0041925 0.0186706 -0.0592992 -0.189681 0.0127781 … 0.00350977 -0.00980804 -0.0588581 -0.184106 0.0357618 0.0221717 -0.0327037 -0.0588589 -0.173952 0.0602371 -0.0466892 0.0965736 ⋮ ⋱ -0.0569682 -0.0202558 0.0659664 -0.137251 0.0495428 -0.0568661 -0.0124937 0.0650277 -0.0297994 -0.0666082 -0.0570419 -0.00968205 0.0710139 -0.0263281 -0.107249 -0.0569222 -0.0224364 0.080825 0.0962447 -0.0373149 -0.0558972 -0.0287199 0.0772194 … 0.0333739 -0.0503132 -0.0556624 -0.0315443 0.0755169 0.00224254 0.131499 -0.0564903 -0.0249832 0.0791765 -0.0486639 -0.101717 -0.0570535 -0.023077 0.0950146 0.0120903 0.166347 -0.0569607 -0.0341167 0.103528 -0.0729744 -0.00130042 -0.0564621 -0.0405538 0.0920885 … -0.0511204 0.0274548 -0.056417 -0.0347061 0.0899392 0.106521 -0.0937512 -0.0567368 -0.0285247 0.0840819 -0.00502509 0.141953 singular values: 283-element Vector{Float64}: 118.02470275379117 16.329470293845436 14.475717736464638 12.775444083672541 11.69266465798729 9.21226889101219 8.392542050051626 7.46454576274542 6.548820440012604 5.435428419444214 5.289305141102475 4.827268476923795 4.5778154866519385 ⋮ 0.03090885712575563 0.02799539349985343 0.02404889024142825 0.021245492946620915 0.019886446902621976 0.015105271901297816 0.01372268127428849 0.011652834020733208 0.0088554722463145 0.006128400462062909 0.003035083065666921 0.002022289284342762 Vt factor: 283×283 Matrix{Float64}: -0.0423701 -0.0427979 -0.0411351 … -0.068855 -0.0668341 -0.0171889 -0.00656414 -0.0199103 0.0209562 0.0628906 0.0466408 0.0331267 0.0305505 -0.0440781 -0.0465621 -0.0198712 -0.0381761 -0.0419185 -0.0260347 -0.0132403 0.0444547 0.0248013 0.0172284 0.0521957 0.0696417 -0.0524257 -0.0480213 -0.0647193 … 0.0136645 -0.035287 0.044339 0.0570753 0.0552864 0.01628 0.0301614 -0.0278711 -0.020096 -0.0207029 0.0563081 0.0330192 0.0786039 0.0862486 0.132617 0.0382286 0.0364599 0.118329 0.134319 0.110004 -0.0535676 -0.0457819 -0.0599732 -0.0444854 -0.0355576 … -0.0484189 -0.0464691 0.116675 0.100495 0.0815002 0.0208476 -0.0066859 -0.0150851 -0.0391767 -0.0298541 -0.0318803 0.0565046 ⋮ ⋱ 0.00952726 -0.0150849 -0.0535822 0.066615 -0.0853499 0.00522386 0.0304889 0.00156495 -0.0610843 -0.0156719 0.0634549 0.00279093 0.0878891 -0.0532867 0.0457019 0.0105359 -0.104826 0.0599095 0.0228524 -0.0766208 -0.0634161 0.0557495 0.00316178 … -0.00494265 -0.00759012 -0.0390847 0.0359929 -0.0412283 -0.0133298 -0.0410064 -0.0299082 0.0574482 0.0300026 0.0126923 0.00129676 -0.0320997 0.0656215 -0.016768 0.0195575 0.00424165 -0.00949598 -0.00387311 0.124957 0.0344757 -0.0255612 -0.0480881 0.131576 -0.00266735 … 0.0501369 -0.0121825 -0.0371619 0.0111629 0.0412245 0.104353 -0.0941486 -0.00565901 -0.0243119 0.0488415 -0.00606318 0.045091
norm(SVD_V.U*diagm(SVD_V.S)*SVD_V.V' - Xgrayvalues)
8.985416478474782e-13
# use the top 4 singular vectors/values to form a new image
u1 = SVD_V.U[:,1]
v1 = SVD_V.V[:,1]
img1 = SVD_V.S[1]*u1*v1'
i = 2
u1 = SVD_V.U[:,i]
v1 = SVD_V.V[:,i]
img1 += SVD_V.S[i]*u1*v1'
i = 3
u1 = SVD_V.U[:,i]
v1 = SVD_V.V[:,i]
img1 += SVD_V.S[i]*u1*v1'
i = 4
u1 = SVD_V.U[:,i]
v1 = SVD_V.V[:,i]
img1 += SVD_V.S[i]*u1*v1'
283×283 Matrix{Float64}: 0.29657 0.318549 0.343035 0.342863 … 0.552066 0.572353 0.457859 0.301302 0.320177 0.348938 0.349233 0.537413 0.557303 0.429426 0.304972 0.31986 0.351791 0.353049 0.521999 0.540254 0.40225 0.313416 0.323556 0.356918 0.359548 0.515346 0.530358 0.386351 0.32363 0.32932 0.362593 0.366602 0.515531 0.526726 0.380811 0.326611 0.328579 0.361655 0.366891 … 0.508294 0.516375 0.370125 0.330997 0.328022 0.361055 0.367526 0.494981 0.499081 0.351644 0.341715 0.331637 0.364936 0.372786 0.477945 0.476419 0.325571 0.347785 0.33259 0.3661 0.374478 0.458421 0.45296 0.299809 0.355413 0.334757 0.367165 0.376502 0.44522 0.434822 0.282863 0.364931 0.338526 0.368979 0.379199 … 0.433394 0.417385 0.268694 0.3731 0.339849 0.367643 0.379321 0.419181 0.396356 0.253065 0.382235 0.341772 0.365584 0.379015 0.408898 0.378267 0.243974 ⋮ ⋱ ⋮ 0.354242 0.358328 0.352706 0.341149 0.446787 0.439007 0.396856 0.34905 0.35294 0.345329 0.334487 0.447845 0.438984 0.403409 0.351978 0.354083 0.345375 0.335005 0.446849 0.435978 0.402846 0.361234 0.35888 0.352543 0.343088 0.437318 0.423923 0.381958 0.361572 0.364431 0.360956 0.347908 … 0.431406 0.42378 0.373937 0.362478 0.367414 0.365122 0.350681 0.430516 0.425185 0.371956 0.371064 0.37998 0.376679 0.359154 0.4395 0.436824 0.385301 0.379477 0.381453 0.376175 0.361214 0.43536 0.425971 0.377987 0.385815 0.382321 0.378765 0.365761 0.426375 0.413328 0.358817 0.376774 0.373795 0.372047 0.360174 … 0.425234 0.413546 0.355563 0.3738 0.372571 0.369703 0.35723 0.427778 0.417003 0.362891 0.372627 0.376319 0.372809 0.358087 0.436367 0.429276 0.377653
Gray.(img1)
As you can see, it's still far away from the original image. Let's try using 100 singular vectors/values.
i = 1:100
u1 = SVD_V.U[:,i]
v1 = SVD_V.V[:,i]
img1 = u1*spdiagm(0=>SVD_V.S[i])*v1'
Gray.(img1)
This looks almost identical to the original image, even though it's not identical to the original image (and we can see that from the norm difference).
norm(Xgrayvalues-img1)
5.528187095958997
Our next problem will still be related to images, but this time we will solve a simple form of the face recognition problem. Let's get the data first.
M = matread("data/face_recog_qr.mat")
Dict{String, Any} with 1 entry: "V2" => [0.08103 0.0729089 … 0.0529805 0.0594823; 0.089725 0.082329 … 0.05180…
Each vector in M["V2"]
is a fase image. Let's reshape the first one and take a look.
q = reshape(M["V2"][:,1],192,168)
Gray.(q)
Now we will go back to the vectorized version of this image, and try to select the images that are most similar to it from the "dictionary" matrix. Let's use b = q[:]
to be the query image. Note that the notation [:]
vectorizes a matrix column wise.
b = q[:]
32256-element Vector{Float64}: 0.0810300261581645 0.08972497847335578 0.0873803822279144 0.0841839709221837 0.08923537113064388 0.08721125794763387 0.09061926007535659 0.09874041385097843 0.08932030421406685 0.08584587638761629 0.09319940629980634 0.08927783767235538 0.0994655711726987 ⋮ 0.10334467325843819 0.09275226549880594 0.0846311117231841 0.07561567634556146 0.06385097046320852 0.06110169995847822 0.05023127567812609 0.03478055114733056 0.027342088151717514 0.028961527075438608 0.024761941927267744 0.027064071630997635
We will remove the first image from the dictionary. The goal is to find the solution of the linear system Ax=b
where A
is the dictionary of all images. In face recognition problem we really want to minimize the norm differece norm(Ax-b)
but the \
actually solves a least squares problem when the matrix at hand is not invertible.
A = M["V2"][:,2:end]
x = A\b #Ax=b
Gray.(reshape(A*x,192,168))
norm(A*x-b)
9.32163450371629e-14
This was an easy problem. Let's try to make the picture harder to recover. We will add some random error.
qv = q+rand(size(q,1),size(q,2))*0.5
qv = qv./maximum(qv)
Gray.(qv)
b = qv[:];
x = A\b
norm(A*x-b)
22.71317697959601
The error is so much bigger this time.
Gray.(reshape(A*x,192,168))
After finishing this notebook, you should be able to:
We can solve a simple form of the face recognition problem even when a face image has been distorted with wrong pixels. Example, one of our inputs was this image:
And we were able to detect this face to be closest to the input image: