Flux in Julia/Learning Julia (Intro_to_Julia)

12. Factorizations and other fun (한글)

딥스탯 2018. 10. 5. 21:32
12_Factorizations_and_other_fun(한글)

Factorizations and other fun

Andreas Noack님의 작성에 기반했습니다.

출처

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

Topics:

  1. 행렬 분해
  2. 특별한 행렬 구조
  3. 일반적인 선형대수
  4. 연습문제

함께보기

시작하기 전에, 선형계를 설정하고, 행렬 분해나 특별한 행렬 구조를 사용하기 위해서 LinearAlgebra 를 불러온다.

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 분해 (LU factorizations)

Julia에서 LU 분해(LU factorizations)를 하려면,

PA = LU

여기서 P는 순열행렬(permutation matrix), L 는 하삼각 단위 대각 행렬 (lower triangular unit diagonal matrix), U는 상삼각행렬(upper triangular matrix)이다.

Julia는 LU 분해를 계산하고, 이를 저장할 때 type을 composite factorization type으로 정의한다.

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}}

특수 속성에 접근해서 분해결과를 각각 추출해낼 수 있다.

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는 분해 object에 있는 method를 빠르게 처리할 수 있다.

예를 들어 원래의 행렬 또는 분해 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

비슷하게, A의 행렬식(determinant)도 계산할 수 있다.

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

QR 분해(QR factorizations)

Julia에서 QR 분해 또한 할 수 있다.

A=QR

여기서 Q는 단위(unitary)/직각(orthogonal) 행렬이고, R은 상삼각(upper triangular) 행렬이다.

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

LU 분해와 마찬가지로, 행렬 Q와 R은 QR 분해 object로부터 추출해낼 수 있다.

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)

고유분해(eigendecompositions), 특이값 분해(singular value decompositions), Hessenberg 분해(Hessenberg factorizations), Schur 분해(Schur decompositions) 등 4가지 분해의 결과는 모두 Factorization type으로 저장된다.

고유분해(eigendecompositions)는 다음과 같이 계산된다.

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

고유값(eigen values)과 고유벡터(eigen vectors)는 Eigen 타입에 특별한 인덱싱(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

또 언급하자면, 분해가 type에 저장되면, 분해의 특성을 이용하는 특수한 method를 작성할 수 있다. 예를 들면, $$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)

행렬구조는 선형대수에서 매우 중요하다. 아래는 이것이 얼마나 중요한 건지 볼 건데, 그 전에 좀 큼 선형계를 설정하자.

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

Julia는 종종 특별한 행렬 구조(special matrix structure)를 추론한다.

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

그러나 가끔 부동소수점 오류가 발생할 수도 있다.

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

다행히도 우리는 이런 구조를 명시적으로 선언할 수도 있는데, Diagonal, Triangular, Symmetric, Hermitian, Tridiagonal, SymTridiagonal 같은 것들이 그 예다.

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

Julia가 Asym, Asym_noisy, Asym_explicit 세 행렬의 eigenvalue를 계산하는데 얼마나 차이가 나는지 비교하자.

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)

이 예에서 Asym_noisy에 Symmetric()을 씌운 것이 5x나 더 효율적이다.

큰 문제

삼중행렬(tridiagonal matrix)를 저장할 때, Tridiagonal나 SymTridiagonal type을 이용하면 잠재적으로 엄청나게 큰 삼중행렬(tridiagonal matrix)도 계산할 수 있다. 다음의 예제는 Matrix type으로 저장된 object라면 일반 PC로 계산할 수 없을 것이다.

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)

보통, 수치 선형대수를 위한 support를 추가하려고 BLAS와 LAPACK subroutine을 묶는다. Float32, Float64, Complex{Float32}, Complex{Float64} type의 구성요소들로 이루어진 행렬은 Julia 또한 지원한다.

그러나, Julia는 일반적인 선형대수(generic linear algebra)또한 지원한다. 예를 들어서, 유리수 벡터와 행렬을 연산할 수 있다.

유리수(Rational numbers)

Julia는 유리수가 기본내장돼있다. 유리수를 입력하기 위해서, 정방향 slash를 2개 쓴다.

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

예제 : 유리선형계 방정식(Rational linear system of equations)

다음의 예제는 유리선형계 방정식(rational linear system of equations)을 부동 소수점 type으로 바꾸지 않고 해결할 수 있는 방법을 보여준다. 유리수로 작업 할 때 오버플로(overflow)가 문제가 쉽게 일어날 수 있으므로 BigInt를 사용한다.

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

연습문제

11.1

다음 행렬 A의 고유값(eigen values)이 뭔지 계산해보자.

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
]

그리고 계산한 값을 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

A의 고유값(eigenvalues)를 이용해서 대각행렬(Diagonal matrix)를 만들자.

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

A의 하삼각(LowerTriangular) 행렬을 만들고, 이를 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]

Julia tutorial 개발 팀이 잘 하고 있는지 평가해 주세요!(Please let us know how we're doing!)

https://tinyurl.com/introJuliaFeedback


저작자표시 비영리 변경금지 (새창열림)