Flux in Julia/Learning Julia (Intro_to_Julia)

12. Factorizations and other fun

딥스탯 2018. 10. 5. 21:30
12_Factorizations_and_other_fun

Factorizations and other fun

Based on work by Andreas Noack

Reference

https://github.com/JuliaComputing/JuliaBoxTutorials/tree/master/introductory-tutorials/intro-to-julia (github : JuliaComputing/JuliaBoxTutorials/introductory-tutorials/intro-to-julia/)

Topics:

  1. Factorizations
  2. Special matrix structures
  3. Generic linear algebra
  4. Exercises

Series

Before we get started, let's set up a linear system and use LinearAlgebra to bring in the factorizations and special matrix structures.

In [1]:
using LinearAlgebra
A = rand(3, 3)
x = fill(1, (3,))
b = A * x
Out[1]:
3-element Array{Float64,1}:
 1.4095257746683503
 1.8634312497010534
 1.3670170455218968

Factorizations

LU factorizations

In Julia we can perform an LU factorization

PA = LU

where P is a permutation matrix, L is lower triangular unit diagonal and U is upper triangular, using lufact.

Julia allows computing the LU factorization and defines a composite factorization type for storing it.

In [2]:
Alu = lu(A)
Out[2]:
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
 1.0       0.0       0.0
 0.615572  1.0       0.0
 0.509803  0.759546  1.0
U factor:
3×3 Array{Float64,2}:
 0.932644   0.799276  0.131511
 0.0       -0.208911  0.428852
 0.0        0.0       0.292487
In [3]:
typeof(Alu)
Out[3]:
LU{Float64,Array{Float64,2}}

The different parts of the factorization can be extracted by accessing their special properties

In [4]:
Alu.P
Out[4]:
3×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0
 1.0  0.0  0.0
In [5]:
Alu.L
Out[5]:
3×3 Array{Float64,2}:
 1.0       0.0       0.0
 0.615572  1.0       0.0
 0.509803  0.759546  1.0
In [6]:
Alu.U
Out[6]:
3×3 Array{Float64,2}:
 0.932644   0.799276  0.131511
 0.0       -0.208911  0.428852
 0.0        0.0       0.292487

Julia can dispatch methods on factorization objects.

For example, we can solve the linear system using either the original matrix or the factorization object.

In [7]:
A\b
Out[7]:
3-element Array{Float64,1}:
 1.0000000000000002
 0.9999999999999998
 1.0
In [8]:
Alu\b
Out[8]:
3-element Array{Float64,1}:
 1.0000000000000002
 0.9999999999999998
 1.0

Similarly, we can calculate the determinant of A using either A or the factorization object

In [9]:
det(A) ≈ det(Alu)
Out[9]:
true

QR factorizations

In Julia we can perform a QR factorization

A=QR

where Q is unitary/orthogonal and R is upper triangular, using qrfact.

In [10]:
Aqr = qr(A)
Out[10]:
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.398232   0.456045  -0.795886 
 -0.781148  -0.62344    0.0336249
 -0.480853   0.635095   0.604512 
R factor:
3×3 Array{Float64,2}:
 -1.19394  -0.859561  -0.620766
  0.0      -0.205043   0.554298
  0.0       0.0       -0.232786

Similarly to the LU factorization, the matrices Q and R can be extracted from the QR factorization object via

In [11]:
Aqr.Q
Out[11]:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.398232   0.456045  -0.795886 
 -0.781148  -0.62344    0.0336249
 -0.480853   0.635095   0.604512
In [12]:
Aqr.R
Out[12]:
3×3 Array{Float64,2}:
 -1.19394  -0.859561  -0.620766
  0.0      -0.205043   0.554298
  0.0       0.0       -0.232786

Eigendecompositions

The results from eigendecompositions, singular value decompositions, Hessenberg factorizations, and Schur decompositions are all stored in Factorization types.

The eigendecomposition can be computed

In [13]:
Asym = A + A'
AsymEig = eigen(Asym)
Out[13]:
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
eigenvalues:
3-element Array{Float64,1}:
 -0.4684087192598678
  0.9149919921470593
  3.122512789652273 
eigenvectors:
3×3 Array{Float64,2}:
 -0.764548   0.181307  -0.618541
  0.325394  -0.719802  -0.613192
  0.556403   0.670084  -0.491327

The values and the vectors can be extracted from the Eigen type by special indexing

In [14]:
AsymEig.values
Out[14]:
3-element Array{Float64,1}:
 -0.4684087192598678
  0.9149919921470593
  3.122512789652273
In [15]:
AsymEig.vectors
Out[15]:
3×3 Array{Float64,2}:
 -0.764548   0.181307  -0.618541
  0.325394  -0.719802  -0.613192
  0.556403   0.670084  -0.491327

Once again, when the factorization is stored in a type, we can dispatch on it and write specialized methods that exploit the properties of the factorization, e.g. that $$A^{-1}=(V\Lambda V^{-1})^{-1}=V\Lambda^{-1}V^{-1}.$$

In [16]:
inv(AsymEig)*Asym
Out[16]:
3×3 Array{Float64,2}:
  1.0           2.77556e-16  -5.9952e-15 
 -1.33227e-15   1.0           4.32987e-15
 -1.65146e-15  -8.04912e-16   1.0

Special matrix structures

Matrix structure is very important in linear algebra. To see how important it is, let's work with a larger linear system

In [17]:
n = 1000
A = randn(n,n);

Julia can often infer special matrix structure

In [18]:
Asym = A + A'
issymmetric(Asym)
Out[18]:
true

but sometimes floating point error might get in the way.

In [19]:
Asym_noisy = copy(Asym)
Asym_noisy[1,2] += 5eps()
Out[19]:
3.03091928625111
In [20]:
issymmetric(Asym_noisy)
Out[20]:
false

Luckily we can declare structure explicitly with, for example, Diagonal, Triangular, Symmetric, Hermitian, Tridiagonal and SymTridiagonal.

In [21]:
Asym_explicit = Symmetric(Asym_noisy);

Let's compare how long it takes Julia to compute the eigenvalues of Asym, Asym_noisy, and Asym_explicit

In [22]:
@time eigvals(Asym);

  0.166887 seconds (112.22 k allocations: 13.368 MiB)
In [23]:
@time eigvals(Asym_noisy);

  1.073496 seconds (18 allocations: 7.921 MiB)
In [24]:
@time eigvals(Asym_explicit);

  0.119242 seconds (7.62 k allocations: 8.378 MiB)

In this example, using Symmetric() on Asym_noisy made our calculations about 5x more efficient :)

A big problem

Using the Tridiagonal and SymTridiagonal types to store tridiagonal matrices makes it possible to work with potentially very large tridiagonal problems. The following problem would not be possible to solve on a laptop if the matrix had to be stored as a (dense) Matrix type.

In [25]:
n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

  0.754640 seconds (484.50 k allocations: 206.799 MiB, 17.43% gc time)
Out[25]:
6.370204015035759

Generic linear algebra

The usual way of adding support for numerical linear algebra is by wrapping BLAS and LAPACK subroutines. For matrices with elements of Float32, Float64, Complex{Float32} or Complex{Float64} this is also what Julia does.

However, Julia also supports generic linear algebra, allowing you to, for example, work with matrices and vectors of rational numbers.

Rational numbers

Julia has rational numbers built in. To construct a rational number, use double forward slashes:

In [26]:
1//2
Out[26]:
1//2

Example: Rational linear system of equations

The following example shows how linear system of equations with rational elements can be solved without promoting to floating point element types. Overflow can easily become a problem when working with rational numbers so we use BigInts.

In [27]:
Arational = Matrix{Rational{BigInt}}(rand(1:10, 3, 3))/10
Out[27]:
3×3 Array{Rational{BigInt},2}:
 2//5   9//10  2//5 
 2//5   1//2   9//10
 3//10  1//5   1//10
In [28]:
x = fill(1, 3)
b = Arational*x
Out[28]:
3-element Array{Rational{BigInt},1}:
 17//10
  9//5 
  3//5
In [29]:
Arational\b
Out[29]:
3-element Array{Rational{BigInt},1}:
 1//1
 1//1
 1//1
In [30]:
lu(Arational)
Out[30]:
LU{Rational{BigInt},Array{Rational{BigInt},2}}
L factor:
3×3 Array{Rational{BigInt},2}:
 1//1   0//1   0//1
 1//1   1//1   0//1
 3//4  19//16  1//1
U factor:
3×3 Array{Rational{BigInt},2}:
 2//5   9//10     2//5  
 0//1  -2//5      1//2  
 0//1   0//1   -127//160

Exercises

11.1

What are the eigenvalues of matrix A?

A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]

and assign it a variable A_eigv

In [31]:
using LinearAlgebra
In [32]:
A =[140   97   74  168  131;
    97  106   89  131   36;
    74   89  152  144   71;
    168  131  144   54  142;
    131   36   71  142   36]
Out[32]:
5×5 Array{Int64,2}:
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
In [33]:
A_eigv = eigen(A).values
Out[33]:
5-element Array{Float64,1}:
 -128.49322764802145 
  -55.887784553057   
   42.752167279318854
   87.16111477514494 
  542.4677301466137
In [34]:
@assert A_eigv ==  [-128.49322764802145, -55.887784553056875, 42.7521672793189, 87.16111477514521, 542.4677301466143]

AssertionError: A_eigv == [-128.49322764802145, -55.887784553056875, 42.7521672793189, 87.16111477514521, 542.4677301466143]

Stacktrace:
 [1] top-level scope at In[34]:1

11.2

Create a Diagonal matrix from the eigenvalues of A.

In [35]:
A_diag = Diagonal(A_eigv)
Out[35]:
5×5 Diagonal{Float64,Array{Float64,1}}:
 -128.493     ⋅        ⋅        ⋅         ⋅   
     ⋅     -55.8878    ⋅        ⋅         ⋅   
     ⋅        ⋅      42.7522    ⋅         ⋅   
     ⋅        ⋅        ⋅      87.1611     ⋅   
     ⋅        ⋅        ⋅        ⋅      542.468
In [36]:
@assert A_diag ==  [-128.493    0.0      0.0      0.0       0.0;
    0.0    -55.8878   0.0      0.0       0.0;
    0.0      0.0     42.7522   0.0       0.0;
    0.0      0.0      0.0     87.1611    0.0;
    0.0 0.0      0.0      0.0     542.468]

AssertionError: A_diag == [-128.493 0.0 0.0 0.0 0.0; 0.0 -55.8878 0.0 0.0 0.0; 0.0 0.0 42.7522 0.0 0.0; 0.0 0.0 0.0 87.1611 0.0; 0.0 0.0 0.0 0.0 542.468]

Stacktrace:
 [1] top-level scope at In[36]:1

11.3

Create a LowerTriangular matrix from A and store it in A_lowertri

In [37]:
A_lowertri = LowerTriangular(A)
Out[37]:
5×5 LowerTriangular{Int64,Array{Int64,2}}:
 140    ⋅    ⋅    ⋅   ⋅
  97  106    ⋅    ⋅   ⋅
  74   89  152    ⋅   ⋅
 168  131  144   54   ⋅
 131   36   71  142  36
In [38]:
@assert A_lowertri ==  [140    0    0    0   0;
  97  106    0    0   0;
  74   89  152    0   0;
 168  131  144   54   0;
 131   36   71  142  36]

Please let us know how we're doing!

https://tinyurl.com/introJuliaFeedback


저작자표시 비영리 변경금지