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. 연습문제

함께보기¶

• http://deepstat.tistory.com/45 (01. Getting started)(in English)
• http://deepstat.tistory.com/46 (01. Getting started(한글))
• http://deepstat.tistory.com/47 (02. Strings)(in English)
• http://deepstat.tistory.com/48 (02. Strings(한글))
• http://deepstat.tistory.com/49 (03. Data structures)(in English)
• http://deepstat.tistory.com/50 (03. Data structures(한글))
• http://deepstat.tistory.com/51 (04. Loops)(in English)
• http://deepstat.tistory.com/52 (04. Loops(한글))
• http://deepstat.tistory.com/53 (05. Conditionals)(in English)
• http://deepstat.tistory.com/54 (05. Conditionals(한글))
• http://deepstat.tistory.com/55 (06. Functions)(in English)
• http://deepstat.tistory.com/56 (06. Functions(한글))
• http://deepstat.tistory.com/57 (07. Packages)(in English)
• http://deepstat.tistory.com/58 (07. Packages(한글))
• http://deepstat.tistory.com/59 (08. Plotting)(in English)
• http://deepstat.tistory.com/60 (08. Plotting(한글))
• http://deepstat.tistory.com/61 (09. Julia is fast)(in English)
• http://deepstat.tistory.com/62 (09. Julia is fast(한글))
• http://deepstat.tistory.com/63 (10. Multiple dispatch)(in English)
• http://deepstat.tistory.com/64 (10. Multiple dispatch(한글))
• http://deepstat.tistory.com/65 (11. Basic linear algebra in Julia)(in English)
• http://deepstat.tistory.com/66 (11. Basic linear algebra in Julia(한글))
• http://deepstat.tistory.com/67 (12. Factorizations and other fun)(in English)

시작하기 전에, 선형계를 설정하고, 행렬 분해나 특별한 행렬 구조를 사용하기 위해서 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